QGIS API Documentation 3.37.0-Master (fdefdf9c27f)
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 for ( int point = 0; point < pnts.count(); ++point )
94 {
95 if ( pnts[point].anchor >= 0 )
96 continue;
97
98 pnts[point].anchor = -2; // make it anchor
99
100 // Find points in threshold
101 double x = pnts[point].x, y = pnts[point].y;
102 const QgsRectangle rect( x - thresh, y - thresh, x + thresh, y + thresh );
103
104 const QList<QgsFeatureId> ids = index.intersects( rect );
105 for ( const QgsFeatureId pointb : ids )
106 {
107 if ( pointb == point )
108 continue;
109
110 const double dx = pnts[pointb].x - pnts[point].x;
111 const double dy = pnts[pointb].y - pnts[point].y;
112 const double dist2 = dx * dx + dy * dy;
113 if ( dist2 > thresh2 )
114 continue; // outside threshold
115
116 if ( pnts[pointb].anchor == -1 )
117 {
118 // doesn't have an anchor yet
119 pnts[pointb].anchor = point;
120 }
121 else if ( pnts[pointb].anchor >= 0 )
122 {
123 // check distance to previously assigned anchor
124 const double dx2 = pnts[pnts[pointb].anchor].x - pnts[pointb].x;
125 const double dy2 = pnts[pnts[pointb].anchor].y - pnts[pointb].y;
126 const double dist2_a = dx2 * dx2 + dy2 * dy2;
127 if ( dist2 < dist2_a )
128 pnts[pointb].anchor = point; // replace old anchor
129 }
130 }
131 }
132}
133
134
135static bool snapPoint( QgsPoint *pt, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts )
136{
137 // Find point ( should always find one point )
138 QList<QgsFeatureId> fids = index.intersects( QgsRectangle( pt->x(), pt->y(), pt->x(), pt->y() ) );
139 Q_ASSERT( fids.count() == 1 );
140
141 const int spoint = fids[0];
142 const int anchor = pnts[spoint].anchor;
143
144 if ( anchor >= 0 )
145 {
146 // to be snapped
147 pt->setX( pnts[anchor].x );
148 pt->setY( pnts[anchor].y );
149 return true;
150 }
151
152 return false;
153}
154
155
156static bool snapLineString( QgsLineString *linestring, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
157{
158 QVector<QgsPoint> newPoints;
159 QVector<int> anchors; // indexes of anchors for vertices
160 const double thresh2 = thresh * thresh;
161 double minDistX, minDistY; // coordinates of the closest point on the segment line
162 bool changed = false;
163
164 // snap vertices
165 for ( int v = 0; v < linestring->numPoints(); v++ )
166 {
167 const double x = linestring->xAt( v );
168 const double y = linestring->yAt( v );
169 const QgsRectangle rect( x, y, x, y );
170
171 // Find point ( should always find one point )
172 QList<QgsFeatureId> fids = index.intersects( rect );
173 Q_ASSERT( fids.count() == 1 );
174
175 const int spoint = fids.first();
176 const int anchor = pnts[spoint].anchor;
177 if ( anchor >= 0 )
178 {
179 // to be snapped
180 linestring->setXAt( v, pnts[anchor].x );
181 linestring->setYAt( v, pnts[anchor].y );
182 anchors.append( anchor ); // point on new location
183 changed = true;
184 }
185 else
186 {
187 anchors.append( spoint ); // old point
188 }
189 }
190
191 // Snap all segments to anchors in threshold
192 for ( int v = 0; v < linestring->numPoints() - 1; v++ )
193 {
194 const double x1 = linestring->xAt( v );
195 const double x2 = linestring->xAt( v + 1 );
196 const double y1 = linestring->yAt( v );
197 const double y2 = linestring->yAt( v + 1 );
198
199 newPoints << linestring->pointN( v );
200
201 // Box
202 double xmin = x1, xmax = x2, ymin = y1, ymax = y2;
203 if ( xmin > xmax )
204 std::swap( xmin, xmax );
205 if ( ymin > ymax )
206 std::swap( ymin, ymax );
207
208 const QgsRectangle rect( xmin - thresh, ymin - thresh, xmax + thresh, ymax + thresh );
209
210 // Find points
211 const QList<QgsFeatureId> fids = index.intersects( rect );
212
213 QVector<AnchorAlongSegment> newVerticesAlongSegment;
214
215 // Snap to anchor in threshold different from end points
216 for ( const QgsFeatureId fid : fids )
217 {
218 const int spoint = fid;
219
220 if ( spoint == anchors[v] || spoint == anchors[v + 1] )
221 continue; // end point
222 if ( pnts[spoint].anchor >= 0 )
223 continue; // point is not anchor
224
225 // Check the distance
226 const double dist2 = QgsGeometryUtilsBase::sqrDistToLine( pnts[spoint].x, pnts[spoint].y, x1, y1, x2, y2, minDistX, minDistY, 0 );
227 // skip points that are behind segment's endpoints or extremely close to them
228 double dx1 = minDistX - x1, dx2 = minDistX - x2;
229 double dy1 = minDistY - y1, dy2 = minDistY - y2;
230 const bool isOnSegment = !qgsDoubleNear( dx1 * dx1 + dy1 * dy1, 0 ) && !qgsDoubleNear( dx2 * dx2 + dy2 * dy2, 0 );
231 if ( isOnSegment && dist2 <= thresh2 )
232 {
233 // an anchor is in the threshold
235 item.anchor = spoint;
236 item.along = QgsPointXY( x1, y1 ).distance( minDistX, minDistY );
237 newVerticesAlongSegment << item;
238 }
239 }
240
241 if ( !newVerticesAlongSegment.isEmpty() )
242 {
243 // sort by distance along the segment
244 std::sort( newVerticesAlongSegment.begin(), newVerticesAlongSegment.end(), []( AnchorAlongSegment p1, AnchorAlongSegment p2 )
245 {
246 return p1.along < p2.along;
247 } );
248
249 // insert new vertices
250 for ( int i = 0; i < newVerticesAlongSegment.count(); i++ )
251 {
252 const int anchor = newVerticesAlongSegment[i].anchor;
253 newPoints << QgsPoint( pnts[anchor].x, pnts[anchor].y, 0 );
254 }
255 changed = true;
256 }
257 }
258
259 // append end point
260 newPoints << linestring->pointN( linestring->numPoints() - 1 );
261
262 // replace linestring's points
263 if ( changed )
264 linestring->setPoints( newPoints );
265
266 return changed;
267}
268
269
270static bool snapGeometry( QgsAbstractGeometry *g, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
271{
272 bool changed = false;
273 if ( QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( g ) )
274 {
275 changed |= snapLineString( linestring, index, pnts, thresh );
276 }
277 else if ( QgsPolygon *polygon = qgsgeometry_cast<QgsPolygon *>( g ) )
278 {
279 if ( QgsLineString *exteriorRing = qgsgeometry_cast<QgsLineString *>( polygon->exteriorRing() ) )
280 changed |= snapLineString( exteriorRing, index, pnts, thresh );
281 for ( int i = 0; i < polygon->numInteriorRings(); ++i )
282 {
283 if ( QgsLineString *interiorRing = qgsgeometry_cast<QgsLineString *>( polygon->interiorRing( i ) ) )
284 changed |= snapLineString( interiorRing, index, pnts, thresh );
285 }
286 }
287 else if ( QgsGeometryCollection *collection = qgsgeometry_cast<QgsGeometryCollection *>( g ) )
288 {
289 for ( int i = 0; i < collection->numGeometries(); ++i )
290 changed |= snapGeometry( collection->geometryN( i ), index, pnts, thresh );
291 }
292 else if ( QgsPoint *pt = qgsgeometry_cast<QgsPoint *>( g ) )
293 {
294 changed |= snapPoint( pt, index, pnts );
295 }
296
297 return changed;
298}
299
300
301int QgsGeometrySnapperSingleSource::run( const QgsFeatureSource &source, QgsFeatureSink &sink, double thresh, QgsFeedback *feedback )
302{
303 // the logic here comes from GRASS implementation of Vect_snap_lines_list()
304
305 long long count = 0;
306 const long long totalCount = source.featureCount() * 2;
307
308 // step 1: record all point locations in a spatial index + extra data structure to keep
309 // reference to which other point they have been snapped to (in the next phase).
310
311 QgsSpatialIndex index;
312 QVector<AnchorPoint> pnts;
313 QgsFeatureRequest request;
314 request.setNoAttributes();
315 QgsFeatureIterator fi = source.getFeatures( request );
316 buildSnapIndex( fi, index, pnts, feedback, count, totalCount );
317
318 if ( feedback->isCanceled() )
319 return 0;
320
321 // step 2: go through all registered points and if not yet marked mark it as anchor and
322 // assign this anchor to all not yet marked points in threshold
323
324 assignAnchors( index, pnts, thresh );
325
326 // step 3: alignment of vertices and segments to the anchors
327 // Go through all lines and:
328 // 1) for all vertices: if not anchor snap it to its anchor
329 // 2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course)
330
331 int modified = 0;
332 QgsFeature f;
333 fi = source.getFeatures();
334 while ( fi.nextFeature( f ) )
335 {
336 if ( feedback->isCanceled() )
337 break;
338
339 QgsGeometry geom = f.geometry();
340 if ( snapGeometry( geom.get(), index, pnts, thresh ) )
341 {
342 f.setGeometry( geom );
343 ++modified;
344 }
345
347
348 ++count;
349 feedback->setProgress( 100. * count / totalCount );
350 }
351
352 return modified;
353}
Abstract base class for all geometries.
virtual QgsRectangle boundingBox() const
Returns the minimal bounding box for the geometry.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
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:167
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:44
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:61
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)
Returns the squared distance between a point and a line.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
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
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:60
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:207
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
void setY(double y)
Sets the point's y-coordinate.
Definition: qgspoint.h:343
void setX(double x)
Sets the point's x-coordinate.
Definition: qgspoint.h:332
Q_GADGET double x
Definition: qgspoint.h:52
double y
Definition: qgspoint.h:53
Polygon geometry type.
Definition: qgspolygon.h:33
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:5207
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...