QGIS API Documentation  3.17.0-Master (a035f434f4)
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, int &count, int totalCount )
53 {
54  QgsFeature f;
55  int pntId = 0;
56 
57  while ( fi.nextFeature( f ) )
58  {
59  if ( feedback->isCanceled() )
60  break;
61 
62  QgsGeometry g = f.geometry();
63 
64  for ( auto it = g.vertices_begin(); it != g.vertices_end(); ++it )
65  {
66  QgsPoint pt = *it;
67  QgsRectangle rect( pt.x(), pt.y(), pt.x(), pt.y() );
68 
69  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  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  QgsRectangle rect( x - thresh, y - thresh, x + thresh, y + thresh );
105 
106  const QList<QgsFeatureId> ids = index.intersects( rect );
107  for ( QgsFeatureId pointb : ids )
108  {
109  if ( pointb == point )
110  continue;
111 
112  double dx = pnts[pointb].x - pnts[point].x;
113  double dy = pnts[pointb].y - pnts[point].y;
114  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  double dx2 = pnts[pnts[pointb].anchor].x - pnts[pointb].x;
128  double dy2 = pnts[pnts[pointb].anchor].y - pnts[pointb].y;
129  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  int spoint = fids[0];
145  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  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  double x = linestring->xAt( v );
171  double y = linestring->yAt( v );
172  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  int spoint = fids.first();
179  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  double x1 = linestring->xAt( v );
198  double x2 = linestring->xAt( v + 1 );
199  double y1 = linestring->yAt( v );
200  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  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 ( QgsFeatureId fid : fids )
220  {
221  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  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  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  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  int count = 0;
309  int 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 }
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a feature to the index.
Wrapper for iterator of features from vector data provider or vector layer.
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
A rectangle specified with double values.
Definition: qgsrectangle.h:41
double y
Definition: qgspoint.h:42
void setPoints(const QgsPointSequence &points)
Resets the line string to match the specified list of points.
void setY(double y) SIP_HOLDGIL
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:280
int numPoints() const override SIP_HOLDGIL
Returns the number of points in the curve.
double x
coordinates of the point
A class to represent a 2D point.
Definition: qgspointxy.h:43
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:316
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:62
An interface for objects which accept features via addFeature(s) methods.
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features ...
Definition: qgsfeatureid.h:28
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:123
void setYAt(int index, double y)
Sets the y-coordinate of the specified node in the line string.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
void setXAt(int index, double x)
Sets the x-coordinate of the specified node in the line string.
Base class for feedback objects to be used for cancellation of something running in a worker thread...
Definition: qgsfeedback.h:43
void setX(double x) SIP_HOLDGIL
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:269
Geometry collection.
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
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.
Abstract base class for all geometries.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
record about anchor being along a segment
A spatial index for QgsFeature objects.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:53
An interface for objects which provide features via a getFeatures method.
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags())
Adds a single feature to the sink.
double distance(double x, double y) const SIP_HOLDGIL
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:196
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
int anchor
Index of the anchor point.
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:139
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
double along
Distance of the anchor point along the segment.
QgsGeometry geometry
Definition: qgsfeature.h:67
bool nextFeature(QgsFeature &f)
Polygon geometry type.
Definition: qgspolygon.h:33
int anchor
Anchor information: 0+ - index of anchor to which this point should be snapped -1 - initial value (un...
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry...
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
QgsRectangle boundingBox() const override SIP_HOLDGIL
Returns the minimal bounding box for the geometry.
Definition: qgspoint.cpp:747
virtual long featureCount() const =0
Returns the number of features contained in the source, or -1 if the feature count is unknown...
double x
Definition: qgspoint.h:41
record about vertex coordinates and index of anchor to which it is snapped