Branch data Line data Source code
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 : :
16 : : #include "qgsgeometrysnappersinglesource.h"
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 : :
28 : : //! record about vertex coordinates and index of anchor to which it is snapped
29 : : struct AnchorPoint
30 : : {
31 : : //! coordinates of the point
32 : : double x, y;
33 : :
34 : : /**
35 : : * Anchor information:
36 : : * 0+ - index of anchor to which this point should be snapped
37 : : * -1 - initial value (undefined)
38 : : * -2 - this point is an anchor, i.e. do not snap this point (snap others to this point)
39 : : */
40 : : int anchor;
41 : : };
42 : :
43 : :
44 : : //! record about anchor being along a segment
45 : : struct AnchorAlongSegment
46 : : {
47 : : int anchor; //!< Index of the anchor point
48 : : double along; //!< Distance of the anchor point along the segment
49 : : };
50 : :
51 : :
52 : 0 : static void buildSnapIndex( QgsFeatureIterator &fi, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, QgsFeedback *feedback, int &count, int totalCount )
53 : : {
54 : 0 : QgsFeature f;
55 : 0 : int pntId = 0;
56 : :
57 : 0 : while ( fi.nextFeature( f ) )
58 : : {
59 : 0 : if ( feedback->isCanceled() )
60 : 0 : break;
61 : :
62 : 0 : QgsGeometry g = f.geometry();
63 : :
64 : 0 : for ( auto it = g.vertices_begin(); it != g.vertices_end(); ++it )
65 : : {
66 : 0 : QgsPoint pt = *it;
67 : 0 : QgsRectangle rect( pt.x(), pt.y(), pt.x(), pt.y() );
68 : :
69 : 0 : QList<QgsFeatureId> ids = index.intersects( rect );
70 : 0 : if ( ids.isEmpty() )
71 : : {
72 : : // add to tree and to structure
73 : 0 : index.addFeature( pntId, pt.boundingBox() );
74 : :
75 : : AnchorPoint xp;
76 : 0 : xp.x = pt.x();
77 : 0 : xp.y = pt.y();
78 : 0 : xp.anchor = -1;
79 : 0 : pnts.append( xp );
80 : 0 : pntId++;
81 : 0 : }
82 : 0 : }
83 : :
84 : 0 : ++count;
85 : 0 : feedback->setProgress( 100. * count / totalCount );
86 : 0 : }
87 : 0 : }
88 : :
89 : :
90 : 0 : static void assignAnchors( QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
91 : : {
92 : 0 : double thresh2 = thresh * thresh;
93 : 0 : int nanchors = 0, ntosnap = 0;
94 : 0 : for ( int point = 0; point < pnts.count(); ++point )
95 : : {
96 : 0 : if ( pnts[point].anchor >= 0 )
97 : 0 : continue;
98 : :
99 : 0 : pnts[point].anchor = -2; // make it anchor
100 : 0 : nanchors++;
101 : :
102 : : // Find points in threshold
103 : 0 : double x = pnts[point].x, y = pnts[point].y;
104 : 0 : QgsRectangle rect( x - thresh, y - thresh, x + thresh, y + thresh );
105 : :
106 : 0 : const QList<QgsFeatureId> ids = index.intersects( rect );
107 : 0 : for ( QgsFeatureId pointb : ids )
108 : : {
109 : 0 : if ( pointb == point )
110 : 0 : continue;
111 : :
112 : 0 : double dx = pnts[pointb].x - pnts[point].x;
113 : 0 : double dy = pnts[pointb].y - pnts[point].y;
114 : 0 : double dist2 = dx * dx + dy * dy;
115 : 0 : if ( dist2 > thresh2 )
116 : 0 : continue; // outside threshold
117 : :
118 : 0 : if ( pnts[pointb].anchor == -1 )
119 : : {
120 : : // doesn't have an anchor yet
121 : 0 : pnts[pointb].anchor = point;
122 : 0 : ntosnap++;
123 : 0 : }
124 : 0 : else if ( pnts[pointb].anchor >= 0 )
125 : : {
126 : : // check distance to previously assigned anchor
127 : 0 : double dx2 = pnts[pnts[pointb].anchor].x - pnts[pointb].x;
128 : 0 : double dy2 = pnts[pnts[pointb].anchor].y - pnts[pointb].y;
129 : 0 : double dist2_a = dx2 * dx2 + dy2 * dy2;
130 : 0 : if ( dist2 < dist2_a )
131 : 0 : pnts[pointb].anchor = point; // replace old anchor
132 : 0 : }
133 : : }
134 : 0 : }
135 : 0 : }
136 : :
137 : :
138 : 0 : static bool snapPoint( QgsPoint *pt, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts )
139 : : {
140 : : // Find point ( should always find one point )
141 : 0 : QList<QgsFeatureId> fids = index.intersects( QgsRectangle( pt->x(), pt->y(), pt->x(), pt->y() ) );
142 : : Q_ASSERT( fids.count() == 1 );
143 : :
144 : 0 : int spoint = fids[0];
145 : 0 : int anchor = pnts[spoint].anchor;
146 : :
147 : 0 : if ( anchor >= 0 )
148 : : {
149 : : // to be snapped
150 : 0 : pt->setX( pnts[anchor].x );
151 : 0 : pt->setY( pnts[anchor].y );
152 : 0 : return true;
153 : : }
154 : :
155 : 0 : return false;
156 : 0 : }
157 : :
158 : :
159 : 0 : static bool snapLineString( QgsLineString *linestring, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
160 : : {
161 : 0 : QVector<QgsPoint> newPoints;
162 : 0 : QVector<int> anchors; // indexes of anchors for vertices
163 : 0 : double thresh2 = thresh * thresh;
164 : : double minDistX, minDistY; // coordinates of the closest point on the segment line
165 : 0 : bool changed = false;
166 : :
167 : : // snap vertices
168 : 0 : for ( int v = 0; v < linestring->numPoints(); v++ )
169 : : {
170 : 0 : double x = linestring->xAt( v );
171 : 0 : double y = linestring->yAt( v );
172 : 0 : QgsRectangle rect( x, y, x, y );
173 : :
174 : : // Find point ( should always find one point )
175 : 0 : QList<QgsFeatureId> fids = index.intersects( rect );
176 : : Q_ASSERT( fids.count() == 1 );
177 : :
178 : 0 : int spoint = fids.first();
179 : 0 : int anchor = pnts[spoint].anchor;
180 : 0 : if ( anchor >= 0 )
181 : : {
182 : : // to be snapped
183 : 0 : linestring->setXAt( v, pnts[anchor].x );
184 : 0 : linestring->setYAt( v, pnts[anchor].y );
185 : 0 : anchors.append( anchor ); // point on new location
186 : 0 : changed = true;
187 : 0 : }
188 : : else
189 : : {
190 : 0 : anchors.append( spoint ); // old point
191 : : }
192 : 0 : }
193 : :
194 : : // Snap all segments to anchors in threshold
195 : 0 : for ( int v = 0; v < linestring->numPoints() - 1; v++ )
196 : : {
197 : 0 : double x1 = linestring->xAt( v );
198 : 0 : double x2 = linestring->xAt( v + 1 );
199 : 0 : double y1 = linestring->yAt( v );
200 : 0 : double y2 = linestring->yAt( v + 1 );
201 : :
202 : 0 : newPoints << linestring->pointN( v );
203 : :
204 : : // Box
205 : 0 : double xmin = x1, xmax = x2, ymin = y1, ymax = y2;
206 : 0 : if ( xmin > xmax )
207 : 0 : std::swap( xmin, xmax );
208 : 0 : if ( ymin > ymax )
209 : 0 : std::swap( ymin, ymax );
210 : :
211 : 0 : QgsRectangle rect( xmin - thresh, ymin - thresh, xmax + thresh, ymax + thresh );
212 : :
213 : : // Find points
214 : 0 : const QList<QgsFeatureId> fids = index.intersects( rect );
215 : :
216 : 0 : QVector<AnchorAlongSegment> newVerticesAlongSegment;
217 : :
218 : : // Snap to anchor in threshold different from end points
219 : 0 : for ( QgsFeatureId fid : fids )
220 : : {
221 : 0 : int spoint = fid;
222 : :
223 : 0 : if ( spoint == anchors[v] || spoint == anchors[v + 1] )
224 : 0 : continue; // end point
225 : 0 : if ( pnts[spoint].anchor >= 0 )
226 : 0 : continue; // point is not anchor
227 : :
228 : : // Check the distance
229 : 0 : 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 : 0 : double dx1 = minDistX - x1, dx2 = minDistX - x2;
232 : 0 : double dy1 = minDistY - y1, dy2 = minDistY - y2;
233 : 0 : bool isOnSegment = !qgsDoubleNear( dx1 * dx1 + dy1 * dy1, 0 ) && !qgsDoubleNear( dx2 * dx2 + dy2 * dy2, 0 );
234 : 0 : if ( isOnSegment && dist2 <= thresh2 )
235 : : {
236 : : // an anchor is in the threshold
237 : : AnchorAlongSegment item;
238 : 0 : item.anchor = spoint;
239 : 0 : item.along = QgsPointXY( x1, y1 ).distance( minDistX, minDistY );
240 : 0 : newVerticesAlongSegment << item;
241 : 0 : }
242 : : }
243 : :
244 : 0 : if ( !newVerticesAlongSegment.isEmpty() )
245 : : {
246 : : // sort by distance along the segment
247 : 0 : std::sort( newVerticesAlongSegment.begin(), newVerticesAlongSegment.end(), []( AnchorAlongSegment p1, AnchorAlongSegment p2 )
248 : : {
249 : 0 : return p1.along < p2.along;
250 : : } );
251 : :
252 : : // insert new vertices
253 : 0 : for ( int i = 0; i < newVerticesAlongSegment.count(); i++ )
254 : : {
255 : 0 : int anchor = newVerticesAlongSegment[i].anchor;
256 : 0 : newPoints << QgsPoint( pnts[anchor].x, pnts[anchor].y, 0 );
257 : 0 : }
258 : 0 : changed = true;
259 : 0 : }
260 : 0 : }
261 : :
262 : : // append end point
263 : 0 : newPoints << linestring->pointN( linestring->numPoints() - 1 );
264 : :
265 : : // replace linestring's points
266 : 0 : if ( changed )
267 : 0 : linestring->setPoints( newPoints );
268 : :
269 : 0 : return changed;
270 : 0 : }
271 : :
272 : :
273 : 0 : static bool snapGeometry( QgsAbstractGeometry *g, QgsSpatialIndex &index, QVector<AnchorPoint> &pnts, double thresh )
274 : : {
275 : 0 : bool changed = false;
276 : 0 : if ( QgsLineString *linestring = qgsgeometry_cast<QgsLineString *>( g ) )
277 : : {
278 : 0 : changed |= snapLineString( linestring, index, pnts, thresh );
279 : 0 : }
280 : 0 : else if ( QgsPolygon *polygon = qgsgeometry_cast<QgsPolygon *>( g ) )
281 : : {
282 : 0 : if ( QgsLineString *exteriorRing = qgsgeometry_cast<QgsLineString *>( polygon->exteriorRing() ) )
283 : 0 : changed |= snapLineString( exteriorRing, index, pnts, thresh );
284 : 0 : for ( int i = 0; i < polygon->numInteriorRings(); ++i )
285 : : {
286 : 0 : if ( QgsLineString *interiorRing = qgsgeometry_cast<QgsLineString *>( polygon->interiorRing( i ) ) )
287 : 0 : changed |= snapLineString( interiorRing, index, pnts, thresh );
288 : 0 : }
289 : 0 : }
290 : 0 : else if ( QgsGeometryCollection *collection = qgsgeometry_cast<QgsGeometryCollection *>( g ) )
291 : : {
292 : 0 : for ( int i = 0; i < collection->numGeometries(); ++i )
293 : 0 : changed |= snapGeometry( collection->geometryN( i ), index, pnts, thresh );
294 : 0 : }
295 : 0 : else if ( QgsPoint *pt = qgsgeometry_cast<QgsPoint *>( g ) )
296 : : {
297 : 0 : changed |= snapPoint( pt, index, pnts );
298 : 0 : }
299 : :
300 : 0 : return changed;
301 : : }
302 : :
303 : :
304 : 0 : 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 : 0 : int count = 0;
309 : 0 : 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 : 0 : QgsSpatialIndex index;
315 : 0 : QVector<AnchorPoint> pnts;
316 : 0 : QgsFeatureRequest request;
317 : 0 : request.setNoAttributes();
318 : 0 : QgsFeatureIterator fi = source.getFeatures( request );
319 : 0 : buildSnapIndex( fi, index, pnts, feedback, count, totalCount );
320 : :
321 : 0 : if ( feedback->isCanceled() )
322 : 0 : 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 : 0 : 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 : 0 : int modified = 0;
335 : 0 : QgsFeature f;
336 : 0 : fi = source.getFeatures();
337 : 0 : while ( fi.nextFeature( f ) )
338 : : {
339 : 0 : if ( feedback->isCanceled() )
340 : 0 : break;
341 : :
342 : 0 : QgsGeometry geom = f.geometry();
343 : 0 : if ( snapGeometry( geom.get(), index, pnts, thresh ) )
344 : : {
345 : 0 : f.setGeometry( geom );
346 : 0 : ++modified;
347 : 0 : }
348 : :
349 : 0 : sink.addFeature( f, QgsFeatureSink::FastInsert );
350 : :
351 : 0 : ++count;
352 : 0 : feedback->setProgress( 100. * count / totalCount );
353 : 0 : }
354 : :
355 : 0 : return modified;
356 : 0 : }
|
