QGIS API Documentation  3.26.3-Buenos Aires (65e4edfdad)
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 }
qgspolygon.h
QgsPointXY::distance
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
QgsFeedback::setProgress
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:76
QgsLineString::pointN
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
Definition: qgslinestring.cpp:996
QgsLineString::setYAt
void setYAt(int index, double y)
Sets the y-coordinate of the specified node in the line string.
Definition: qgslinestring.cpp:1068
QgsLineString::xAt
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
Definition: qgslinestring.cpp:1045
qgsgeometrysnappersinglesource.h
QgsPolygon
Polygon geometry type.
Definition: qgspolygon.h:33
qgslinestring.h
AnchorPoint::x
double x
coordinates of the point
Definition: qgsgeometrysnappersinglesource.cpp:45
QgsPoint
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:48
QgsGeometrySnapperSingleSource::run
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...
Definition: qgsgeometrysnappersinglesource.cpp:304
QgsLineString::yAt
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
Definition: qgslinestring.cpp:1053
QgsFeatureSink::addFeature
virtual bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags())
Adds a single feature to the sink.
Definition: qgsfeaturesink.cpp:20
qgsfeatureiterator.h
QgsGeometryUtils::sqrDistToLine
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.
Definition: qgsgeometryutils.cpp:203
qgsfeaturesink.h
QgsFeature::geometry
QgsGeometry geometry
Definition: qgsfeature.h:71
QgsFeedback::isCanceled
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:67
QgsFeatureSource
An interface for objects which provide features via a getFeatures method.
Definition: qgsfeaturesource.h:37
QgsLineString::setXAt
void setXAt(int index, double x)
Sets the x-coordinate of the specified node in the line string.
Definition: qgslinestring.cpp:1061
AnchorAlongSegment
record about anchor being along a segment
Definition: qgsgeometrysnappersinglesource.cpp:45
QgsLineString
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:44
QgsRectangle
A rectangle specified with double values.
Definition: qgsrectangle.h:41
AnchorAlongSegment::along
double along
Distance of the anchor point along the segment.
Definition: qgsgeometrysnappersinglesource.cpp:48
QgsPoint::y
double y
Definition: qgspoint.h:70
QgsFeatureRequest
This class wraps a request for features to a vector layer (or directly its vector data provider).
Definition: qgsfeaturerequest.h:83
QgsGeometryCollection
Geometry collection.
Definition: qgsgeometrycollection.h:36
AnchorPoint
record about vertex coordinates and index of anchor to which it is snapped
Definition: qgsgeometrysnappersinglesource.cpp:29
QgsFeature::setGeometry
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
Definition: qgsfeature.cpp:170
qgsDoubleNear
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:2265
QgsFeedback
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:44
QgsLineString::setPoints
void setPoints(size_t size, const double *x, const double *y, const double *z=nullptr, const double *m=nullptr)
Resets the line string to match the specified point data.
Definition: qgslinestring.cpp:1092
QgsFeatureSource::getFeatures
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
QgsPoint::setX
void setX(double x) SIP_HOLDGIL
Sets the point's x-coordinate.
Definition: qgspoint.h:297
qgsfeaturesource.h
QgsFeatureSource::featureCount
virtual long long featureCount() const =0
Returns the number of features contained in the source, or -1 if the feature count is unknown.
QgsSpatialIndex
A spatial index for QgsFeature objects.
Definition: qgsspatialindex.h:67
AnchorAlongSegment::anchor
int anchor
Index of the anchor point.
Definition: qgsgeometrysnappersinglesource.cpp:47
QgsFeatureRequest::setNoAttributes
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
Definition: qgsfeaturerequest.cpp:235
QgsAbstractGeometry
Abstract base class for all geometries.
Definition: qgsabstractgeometry.h:79
qgsgeometryutils.h
QgsPointXY
A class to represent a 2D point.
Definition: qgspointxy.h:58
AnchorPoint::anchor
int anchor
Anchor information: 0+ - index of anchor to which this point should be snapped -1 - initial value (un...
Definition: qgsgeometrysnappersinglesource.cpp:53
QgsGeometry::get
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
Definition: qgsgeometry.cpp:141
QgsFeatureIterator::nextFeature
bool nextFeature(QgsFeature &f)
Definition: qgsfeatureiterator.h:399
QgsGeometry
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:124
AnchorPoint::y
double y
Definition: qgsgeometrysnappersinglesource.cpp:45
QgsSpatialIndex::intersects
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
Definition: qgsspatialindex.cpp:492
qgsgeometrycollection.h
qgsspatialindex.h
QgsFeature
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:55
QgsGeometry::vertices_end
QgsAbstractGeometry::vertex_iterator vertices_end() const
Returns STL-style iterator pointing to the imaginary vertex after the last vertex of the geometry.
Definition: qgsgeometry.cpp:1989
QgsPoint::setY
void setY(double y) SIP_HOLDGIL
Sets the point's y-coordinate.
Definition: qgspoint.h:308
QgsGeometry::vertices_begin
QgsAbstractGeometry::vertex_iterator vertices_begin() const
Returns STL-style iterator pointing to the first vertex of the geometry.
Definition: qgsgeometry.cpp:1982
QgsLineString::numPoints
int numPoints() const override SIP_HOLDGIL
Returns the number of points in the curve.
Definition: qgslinestring.cpp:986
QgsFeatureIterator
Wrapper for iterator of features from vector data provider or vector layer.
Definition: qgsfeatureiterator.h:289
qgsfeedback.h
QgsPoint::boundingBox
QgsRectangle boundingBox() const override SIP_HOLDGIL
Returns the minimal bounding box for the geometry.
Definition: qgspoint.cpp:772
QgsSpatialIndex::addFeature
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a feature to the index.
Definition: qgsspatialindex.cpp:407
QgsPoint::x
double x
Definition: qgspoint.h:69
QgsFeatureSink
An interface for objects which accept features via addFeature(s) methods.
Definition: qgsfeaturesink.h:33
QgsFeatureSink::FastInsert
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
Definition: qgsfeaturesink.h:70
QgsFeatureId
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28