QGIS API Documentation  3.21.0-Master (5b68dc587e)
qgsgeometrysnappersinglesource.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgeometrysnappersinglesource.cpp
3  ---------------------
4  Date : May 2018
5  Copyright : (C) 2018 by Martin Dobias
6  Email : wonder dot sk at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
17 
18 #include "qgsfeatureiterator.h"
19 #include "qgsfeaturesink.h"
20 #include "qgsfeaturesource.h"
21 #include "qgsfeedback.h"
22 #include "qgsgeometrycollection.h"
23 #include "qgsgeometryutils.h"
24 #include "qgslinestring.h"
25 #include "qgspolygon.h"
26 #include "qgsspatialindex.h"
27 
30 {
32  double x, y;
33 
40  int anchor;
41 };
42 
43 
46 {
47  int anchor;
48  double along;
49 };
50 
51 
52 static void buildSnapIndex( QgsFeatureIterator &fi, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, QgsFeedback *feedback, long long &count, long long totalCount )
53 {
54  QgsFeature f;
55  int pntId = 0;
56 
57  while ( fi.nextFeature( f ) )
58  {
59  if ( feedback->isCanceled() )
60  break;
61 
62  const QgsGeometry g = f.geometry();
63 
64  for ( auto it = g.vertices_begin(); it != g.vertices_end(); ++it )
65  {
66  const QgsPoint pt = *it;
67  const QgsRectangle rect( pt.x(), pt.y(), pt.x(), pt.y() );
68 
69  const QList<QgsFeatureId> ids = index.intersects( rect );
70  if ( ids.isEmpty() )
71  {
72  // add to tree and to structure
73  index.addFeature( pntId, pt.boundingBox() );
74 
75  AnchorPoint xp;
76  xp.x = pt.x();
77  xp.y = pt.y();
78  xp.anchor = -1;
79  pnts.append( xp );
80  pntId++;
81  }
82  }
83 
84  ++count;
85  feedback->setProgress( 100. * count / totalCount );
86  }
87 }
88 
89 
90 static void assignAnchors( QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
91 {
92  const double thresh2 = thresh * thresh;
93  int nanchors = 0, ntosnap = 0;
94  for ( int point = 0; point < pnts.count(); ++point )
95  {
96  if ( pnts[point].anchor >= 0 )
97  continue;
98 
99  pnts[point].anchor = -2; // make it anchor
100  nanchors++;
101 
102  // Find points in threshold
103  double x = pnts[point].x, y = pnts[point].y;
104  const QgsRectangle rect( x - thresh, y - thresh, x + thresh, y + thresh );
105 
106  const QList<QgsFeatureId> ids = index.intersects( rect );
107  for ( const QgsFeatureId pointb : ids )
108  {
109  if ( pointb == point )
110  continue;
111 
112  const double dx = pnts[pointb].x - pnts[point].x;
113  const double dy = pnts[pointb].y - pnts[point].y;
114  const double dist2 = dx * dx + dy * dy;
115  if ( dist2 > thresh2 )
116  continue; // outside threshold
117 
118  if ( pnts[pointb].anchor == -1 )
119  {
120  // doesn't have an anchor yet
121  pnts[pointb].anchor = point;
122  ntosnap++;
123  }
124  else if ( pnts[pointb].anchor >= 0 )
125  {
126  // check distance to previously assigned anchor
127  const double dx2 = pnts[pnts[pointb].anchor].x - pnts[pointb].x;
128  const double dy2 = pnts[pnts[pointb].anchor].y - pnts[pointb].y;
129  const double dist2_a = dx2 * dx2 + dy2 * dy2;
130  if ( dist2 < dist2_a )
131  pnts[pointb].anchor = point; // replace old anchor
132  }
133  }
134  }
135 }
136 
137 
138 static bool snapPoint( QgsPoint *pt, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts )
139 {
140  // Find point ( should always find one point )
141  QList<QgsFeatureId> fids = index.intersects( QgsRectangle( pt->x(), pt->y(), pt->x(), pt->y() ) );
142  Q_ASSERT( fids.count() == 1 );
143 
144  const int spoint = fids[0];
145  const int anchor = pnts[spoint].anchor;
146 
147  if ( anchor >= 0 )
148  {
149  // to be snapped
150  pt->setX( pnts[anchor].x );
151  pt->setY( pnts[anchor].y );
152  return true;
153  }
154 
155  return false;
156 }
157 
158 
159 static bool snapLineString( QgsLineString *linestring, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
160 {
161  QVector<QgsPoint> newPoints;
162  QVector<int> anchors; // indexes of anchors for vertices
163  const double thresh2 = thresh * thresh;
164  double minDistX, minDistY; // coordinates of the closest point on the segment line
165  bool changed = false;
166 
167  // snap vertices
168  for ( int v = 0; v < linestring->numPoints(); v++ )
169  {
170  const double x = linestring->xAt( v );
171  const double y = linestring->yAt( v );
172  const QgsRectangle rect( x, y, x, y );
173 
174  // Find point ( should always find one point )
175  QList<QgsFeatureId> fids = index.intersects( rect );
176  Q_ASSERT( fids.count() == 1 );
177 
178  const int spoint = fids.first();
179  const int anchor = pnts[spoint].anchor;
180  if ( anchor >= 0 )
181  {
182  // to be snapped
183  linestring->setXAt( v, pnts[anchor].x );
184  linestring->setYAt( v, pnts[anchor].y );
185  anchors.append( anchor ); // point on new location
186  changed = true;
187  }
188  else
189  {
190  anchors.append( spoint ); // old point
191  }
192  }
193 
194  // Snap all segments to anchors in threshold
195  for ( int v = 0; v < linestring->numPoints() - 1; v++ )
196  {
197  const double x1 = linestring->xAt( v );
198  const double x2 = linestring->xAt( v + 1 );
199  const double y1 = linestring->yAt( v );
200  const double y2 = linestring->yAt( v + 1 );
201 
202  newPoints << linestring->pointN( v );
203 
204  // Box
205  double xmin = x1, xmax = x2, ymin = y1, ymax = y2;
206  if ( xmin > xmax )
207  std::swap( xmin, xmax );
208  if ( ymin > ymax )
209  std::swap( ymin, ymax );
210 
211  const QgsRectangle rect( xmin - thresh, ymin - thresh, xmax + thresh, ymax + thresh );
212 
213  // Find points
214  const QList<QgsFeatureId> fids = index.intersects( rect );
215 
216  QVector<AnchorAlongSegment> newVerticesAlongSegment;
217 
218  // Snap to anchor in threshold different from end points
219  for ( const QgsFeatureId fid : fids )
220  {
221  const int spoint = fid;
222 
223  if ( spoint == anchors[v] || spoint == anchors[v + 1] )
224  continue; // end point
225  if ( pnts[spoint].anchor >= 0 )
226  continue; // point is not anchor
227 
228  // Check the distance
229  const double dist2 = QgsGeometryUtils::sqrDistToLine( pnts[spoint].x, pnts[spoint].y, x1, y1, x2, y2, minDistX, minDistY, 0 );
230  // skip points that are behind segment's endpoints or extremely close to them
231  double dx1 = minDistX - x1, dx2 = minDistX - x2;
232  double dy1 = minDistY - y1, dy2 = minDistY - y2;
233  const bool isOnSegment = !qgsDoubleNear( dx1 * dx1 + dy1 * dy1, 0 ) && !qgsDoubleNear( dx2 * dx2 + dy2 * dy2, 0 );
234  if ( isOnSegment && dist2 <= thresh2 )
235  {
236  // an anchor is in the threshold
237  AnchorAlongSegment item;
238  item.anchor = spoint;
239  item.along = QgsPointXY( x1, y1 ).distance( minDistX, minDistY );
240  newVerticesAlongSegment << item;
241  }
242  }
243 
244  if ( !newVerticesAlongSegment.isEmpty() )
245  {
246  // sort by distance along the segment
247  std::sort( newVerticesAlongSegment.begin(), newVerticesAlongSegment.end(), []( AnchorAlongSegment p1, AnchorAlongSegment p2 )
248  {
249  return p1.along < p2.along;
250  } );
251 
252  // insert new vertices
253  for ( int i = 0; i < newVerticesAlongSegment.count(); i++ )
254  {
255  const int anchor = newVerticesAlongSegment[i].anchor;
256  newPoints << QgsPoint( pnts[anchor].x, pnts[anchor].y, 0 );
257  }
258  changed = true;
259  }
260  }
261 
262  // append end point
263  newPoints << linestring->pointN( linestring->numPoints() - 1 );
264 
265  // replace linestring's points
266  if ( changed )
267  linestring->setPoints( newPoints );
268 
269  return changed;
270 }
271 
272 
273 static bool snapGeometry( QgsAbstractGeometry *g, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
274 {
275  bool changed = false;
276  if ( QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( g ) )
277  {
278  changed |= snapLineString( linestring, index, pnts, thresh );
279  }
280  else if ( QgsPolygon *polygon = qgsgeometry_cast<QgsPolygon *>( g ) )
281  {
282  if ( QgsLineString *exteriorRing = qgsgeometry_cast<QgsLineString *>( polygon->exteriorRing() ) )
283  changed |= snapLineString( exteriorRing, index, pnts, thresh );
284  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
285  {
286  if ( QgsLineString *interiorRing = qgsgeometry_cast<QgsLineString *>( polygon->interiorRing( i ) ) )
287  changed |= snapLineString( interiorRing, index, pnts, thresh );
288  }
289  }
290  else if ( QgsGeometryCollection *collection = qgsgeometry_cast<QgsGeometryCollection *>( g ) )
291  {
292  for ( int i = 0; i < collection->numGeometries(); ++i )
293  changed |= snapGeometry( collection->geometryN( i ), index, pnts, thresh );
294  }
295  else if ( QgsPoint *pt = qgsgeometry_cast<QgsPoint *>( g ) )
296  {
297  changed |= snapPoint( pt, index, pnts );
298  }
299 
300  return changed;
301 }
302 
303 
304 int QgsGeometrySnapperSingleSource::run( const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback )
305 {
306  // the logic here comes from GRASS implementation of Vect_snap_lines_list()
307 
308  long long count = 0;
309  const long long totalCount = source.featureCount() * 2;
310 
311  // step 1: record all point locations in a spatial index + extra data structure to keep
312  // reference to which other point they have been snapped to (in the next phase).
313 
314  QgsSpatialIndex index;
315  QVector<AnchorPoint> pnts;
316  QgsFeatureRequest request;
317  request.setNoAttributes();
318  QgsFeatureIterator fi = source.getFeatures( request );
319  buildSnapIndex( fi, index, pnts, feedback, count, totalCount );
320 
321  if ( feedback->isCanceled() )
322  return 0;
323 
324  // step 2: go through all registered points and if not yet marked mark it as anchor and
325  // assign this anchor to all not yet marked points in threshold
326 
327  assignAnchors( index, pnts, thresh );
328 
329  // step 3: alignment of vertices and segments to the anchors
330  // Go through all lines and:
331  // 1) for all vertices: if not anchor snap it to its anchor
332  // 2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course)
333 
334  int modified = 0;
335  QgsFeature f;
336  fi = source.getFeatures();
337  while ( fi.nextFeature( f ) )
338  {
339  if ( feedback->isCanceled() )
340  break;
341 
342  QgsGeometry geom = f.geometry();
343  if ( snapGeometry( geom.get(), index, pnts, thresh ) )
344  {
345  f.setGeometry( geom );
346  ++modified;
347  }
348 
350 
351  ++count;
352  feedback->setProgress( 100. * count / totalCount );
353  }
354 
355  return modified;
356 }
Abstract base class for all geometries.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
This class wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
An interface for objects which accept features via addFeature(s) methods.
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags())
Adds a single feature to the sink.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
An interface for objects which provide features via a getFeatures method.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
virtual long long featureCount() const =0
Returns the number of features contained in the source, or -1 if the feature count is unknown.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QgsGeometry geometry
Definition: qgsfeature.h:67
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
Definition: qgsfeature.cpp:145
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
Geometry collection.
static int run(const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback)
Run the algorithm on given source and output results to the sink, using threshold value in the source...
static double sqrDistToLine(double ptX, double ptY, double x1, double y1, double x2, double y2, double &minDistX, double &minDistY, double epsilon) SIP_HOLDGIL
Returns the squared distance between a point and a line.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:124
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:44
int numPoints() const override SIP_HOLDGIL
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
void setPoints(const QgsPointSequence &points)
Resets the line string to match the specified list of points.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
void setYAt(int index, double y)
Sets the y-coordinate of the specified node in the line string.
void setXAt(int index, double x)
Sets the x-coordinate of the specified node in the line string.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
A class to represent a 2D point.
Definition: qgspointxy.h:59
double distance(double x, double y) const SIP_HOLDGIL
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:211
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
void setX(double x) SIP_HOLDGIL
Sets the point's x-coordinate.
Definition: qgspoint.h:280
Q_GADGET double x
Definition: qgspoint.h:52
void setY(double y) SIP_HOLDGIL
Sets the point's y-coordinate.
Definition: qgspoint.h:291
QgsRectangle boundingBox() const override SIP_HOLDGIL
Returns the minimal bounding box for the geometry.
Definition: qgspoint.cpp:772
double y
Definition: qgspoint.h:53
Polygon geometry type.
Definition: qgspolygon.h:34
A rectangle specified with double values.
Definition: qgsrectangle.h:42
A spatial index for QgsFeature objects.
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a feature to the index.
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:1234
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28
record about anchor being along a segment
int anchor
Index of the anchor point.
double along
Distance of the anchor point along the segment.
record about vertex coordinates and index of anchor to which it is snapped
double x
coordinates of the point
int anchor
Anchor information: 0+ - index of anchor to which this point should be snapped -1 - initial value (un...