QGIS API Documentation 3.27.0-Master (f261cc1f8b)
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"
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
52static 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
90static 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
138static 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
159static 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
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
273static 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
304int 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:170
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:125
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:45
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(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.
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:2396
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...