QGIS API Documentation  3.21.0-Master (5b68dc587e)
qgsgeometrysnapper.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * qgsgeometrysnapper.cpp *
3  * ------------------- *
4  * copyright : (C) 2014 by Sandro Mani / Sourcepole AG *
5  * email : [email protected] *
6  ***************************************************************************/
7 
8 /***************************************************************************
9  * *
10  * This program is free software; you can redistribute it and/or modify *
11  * it under the terms of the GNU General Public License as published by *
12  * the Free Software Foundation; either version 2 of the License, or *
13  * (at your option) any later version. *
14  * *
15  ***************************************************************************/
16 
17 #include "qgsfeatureiterator.h"
18 #include "qgsgeometry.h"
19 #include "qgsvectorlayer.h"
20 #include "qgsgeometrysnapper.h"
21 #include "qgsvectordataprovider.h"
22 #include "qgsgeometryutils.h"
23 #include "qgsmapsettings.h"
24 #include "qgssurface.h"
25 #include "qgsmultisurface.h"
26 #include "qgscurve.h"
27 #include "qgslinestring.h"
28 
29 #include <QtConcurrentMap>
30 #include <geos_c.h>
31 
33 
34 QgsSnapIndex::PointSnapItem::PointSnapItem( const QgsSnapIndex::CoordIdx *_idx, bool isEndPoint )
35  : SnapItem( isEndPoint ? QgsSnapIndex::SnapEndPoint : QgsSnapIndex::SnapPoint )
36  , idx( _idx )
37 {}
38 
39 QgsPoint QgsSnapIndex::PointSnapItem::getSnapPoint( const QgsPoint &/*p*/ ) const
40 {
41  return idx->point();
42 }
43 
44 QgsSnapIndex::SegmentSnapItem::SegmentSnapItem( const QgsSnapIndex::CoordIdx *_idxFrom, const QgsSnapIndex::CoordIdx *_idxTo )
45  : SnapItem( QgsSnapIndex::SnapSegment )
46  , idxFrom( _idxFrom )
47  , idxTo( _idxTo )
48 {}
49 
50 QgsPoint QgsSnapIndex::SegmentSnapItem::getSnapPoint( const QgsPoint &p ) const
51 {
52  return QgsGeometryUtils::projectPointOnSegment( p, idxFrom->point(), idxTo->point() );
53 }
54 
55 bool QgsSnapIndex::SegmentSnapItem::getIntersection( const QgsPoint &p1, const QgsPoint &p2, QgsPoint &inter ) const
56 {
57  const QgsPoint &q1 = idxFrom->point(), & q2 = idxTo->point();
58  QgsVector v( p2.x() - p1.x(), p2.y() - p1.y() );
59  QgsVector w( q2.x() - q1.x(), q2.y() - q1.y() );
60  const double vl = v.length();
61  const double wl = w.length();
62 
63  if ( qgsDoubleNear( vl, 0, 0.000000000001 ) || qgsDoubleNear( wl, 0, 0.000000000001 ) )
64  {
65  return false;
66  }
67  v = v / vl;
68  w = w / wl;
69 
70  const double d = v.y() * w.x() - v.x() * w.y();
71 
72  if ( d == 0 )
73  return false;
74 
75  const double dx = q1.x() - p1.x();
76  const double dy = q1.y() - p1.y();
77  const double k = ( dy * w.x() - dx * w.y() ) / d;
78 
79  inter = QgsPoint( p1.x() + v.x() * k, p1.y() + v.y() * k );
80 
81  const double lambdav = QgsVector( inter.x() - p1.x(), inter.y() - p1.y() ) * v;
82  if ( lambdav < 0. + 1E-8 || lambdav > vl - 1E-8 )
83  return false;
84 
85  const double lambdaw = QgsVector( inter.x() - q1.x(), inter.y() - q1.y() ) * w;
86  return !( lambdaw < 0. + 1E-8 || lambdaw >= wl - 1E-8 );
87 }
88 
89 bool QgsSnapIndex::SegmentSnapItem::getProjection( const QgsPoint &p, QgsPoint &pProj )
90 {
91  const QgsPoint &s1 = idxFrom->point();
92  const QgsPoint &s2 = idxTo->point();
93  const double nx = s2.y() - s1.y();
94  const double ny = -( s2.x() - s1.x() );
95  const double t = ( p.x() * ny - p.y() * nx - s1.x() * ny + s1.y() * nx ) / ( ( s2.x() - s1.x() ) * ny - ( s2.y() - s1.y() ) * nx );
96  if ( t < 0. || t > 1. )
97  {
98  return false;
99  }
100  pProj = QgsPoint( s1.x() + ( s2.x() - s1.x() ) * t, s1.y() + ( s2.y() - s1.y() ) * t );
101  return true;
102 }
103 
104 bool QgsSnapIndex::SegmentSnapItem::withinDistance( const QgsPoint &p, const double tolerance )
105 {
106  double minDistX, minDistY;
107  const double distance = QgsGeometryUtils::sqrDistToLine( p.x(), p.y(), idxFrom->point().x(), idxFrom->point().y(), idxTo->point().x(), idxTo->point().y(), minDistX, minDistY, 4 * std::numeric_limits<double>::epsilon() );
108  return distance <= tolerance;
109 }
110 
112 
113 QgsSnapIndex::QgsSnapIndex()
114 {
115  mSTRTree = GEOSSTRtree_create_r( QgsGeos::getGEOSHandler(), ( size_t )10 );
116 }
117 
118 QgsSnapIndex::~QgsSnapIndex()
119 {
120  qDeleteAll( mCoordIdxs );
121  qDeleteAll( mSnapItems );
122 
123  GEOSSTRtree_destroy_r( QgsGeos::getGEOSHandler(), mSTRTree );
124 }
125 
126 void QgsSnapIndex::addPoint( const CoordIdx *idx, bool isEndPoint )
127 {
128  const QgsPoint p = idx->point();
129 
130  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
131 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
132  geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, p.x(), p.y() ) );
133 #else
134  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
135  GEOSCoordSeq_setX_r( geosctxt, seq, 0, p.x() );
136  GEOSCoordSeq_setY_r( geosctxt, seq, 0, p.y() );
137  geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
138 #endif
139 
140  PointSnapItem *item = new PointSnapItem( idx, isEndPoint );
141  GEOSSTRtree_insert_r( geosctxt, mSTRTree, point.get(), item );
142 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR<9
143  mSTRTreeItems.emplace_back( std::move( point ) );
144 #endif
145  mSnapItems << item;
146 }
147 
148 void QgsSnapIndex::addSegment( const CoordIdx *idxFrom, const CoordIdx *idxTo )
149 {
150  const QgsPoint pointFrom = idxFrom->point();
151  const QgsPoint pointTo = idxTo->point();
152 
153  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
154 
155  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 2, 2 );
156 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
157  GEOSCoordSeq_setXY_r( geosctxt, coord, 0, pointFrom.x(), pointFrom.y() );
158  GEOSCoordSeq_setXY_r( geosctxt, coord, 1, pointTo.x(), pointTo.y() );
159 #else
160  GEOSCoordSeq_setX_r( geosctxt, coord, 0, pointFrom.x() );
161  GEOSCoordSeq_setY_r( geosctxt, coord, 0, pointFrom.y() );
162  GEOSCoordSeq_setX_r( geosctxt, coord, 1, pointTo.x() );
163  GEOSCoordSeq_setY_r( geosctxt, coord, 1, pointTo.y() );
164 #endif
165  geos::unique_ptr segment( GEOSGeom_createLineString_r( geosctxt, coord ) );
166 
167  SegmentSnapItem *item = new SegmentSnapItem( idxFrom, idxTo );
168  GEOSSTRtree_insert_r( geosctxt, mSTRTree, segment.get(), item );
169 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR<9
170  mSTRTreeItems.push_back( std::move( segment ) );
171 #endif
172  mSnapItems << item;
173 }
174 
175 void QgsSnapIndex::addGeometry( const QgsAbstractGeometry *geom )
176 {
177  for ( int iPart = 0, nParts = geom->partCount(); iPart < nParts; ++iPart )
178  {
179  for ( int iRing = 0, nRings = geom->ringCount( iPart ); iRing < nRings; ++iRing )
180  {
181  int nVerts = geom->vertexCount( iPart, iRing );
182 
183  if ( qgsgeometry_cast< const QgsSurface * >( geom ) )
184  nVerts--;
185  else if ( const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( geom ) )
186  {
187  if ( curve->isClosed() )
188  nVerts--;
189  }
190 
191  for ( int iVert = 0; iVert < nVerts; ++iVert )
192  {
193  CoordIdx *idx = new CoordIdx( geom, QgsVertexId( iPart, iRing, iVert ) );
194  CoordIdx *idx1 = new CoordIdx( geom, QgsVertexId( iPart, iRing, iVert + 1 ) );
195  mCoordIdxs.append( idx );
196  mCoordIdxs.append( idx1 );
197  addPoint( idx, iVert == 0 || iVert == nVerts - 1 );
198  if ( iVert < nVerts - 1 )
199  addSegment( idx, idx1 );
200  }
201  }
202  }
203 }
204 
205 struct _GEOSQueryCallbackData
206 {
207  QList< QgsSnapIndex::SnapItem * > *list;
208 };
209 
210 void _GEOSQueryCallback( void *item, void *userdata )
211 {
212  reinterpret_cast<_GEOSQueryCallbackData *>( userdata )->list->append( static_cast<QgsSnapIndex::SnapItem *>( item ) );
213 }
214 
215 QgsPoint QgsSnapIndex::getClosestSnapToPoint( const QgsPoint &startPoint, const QgsPoint &midPoint )
216 {
217  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
218 
219  // Look for intersections on segment from the target point to the point opposite to the point reference point
220  // p2 = p1 + 2 * (q - p1)
221  const QgsPoint endPoint( 2 * midPoint.x() - startPoint.x(), 2 * midPoint.y() - startPoint.y() );
222 
223  QgsPoint minPoint = startPoint;
224  double minDistance = std::numeric_limits<double>::max();
225 
226  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 2, 2 );
227 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
228  GEOSCoordSeq_setXY_r( geosctxt, coord, 0, startPoint.x(), startPoint.y() );
229  GEOSCoordSeq_setXY_r( geosctxt, coord, 1, endPoint.x(), endPoint.y() );
230 #else
231  GEOSCoordSeq_setX_r( geosctxt, coord, 0, startPoint.x() );
232  GEOSCoordSeq_setY_r( geosctxt, coord, 0, startPoint.y() );
233  GEOSCoordSeq_setX_r( geosctxt, coord, 1, endPoint.x() );
234  GEOSCoordSeq_setY_r( geosctxt, coord, 1, endPoint.y() );
235 #endif
236  geos::unique_ptr searchDiagonal( GEOSGeom_createLineString_r( geosctxt, coord ) );
237 
238  QList<SnapItem *> items;
239  struct _GEOSQueryCallbackData callbackData;
240  callbackData.list = &items;
241  GEOSSTRtree_query_r( geosctxt, mSTRTree, searchDiagonal.get(), _GEOSQueryCallback, &callbackData );
242  for ( const SnapItem *item : std::as_const( items ) )
243  {
244  if ( item->type == SnapSegment )
245  {
246  QgsPoint inter;
247  if ( static_cast<const SegmentSnapItem *>( item )->getIntersection( startPoint, endPoint, inter ) )
248  {
249  const double dist = QgsGeometryUtils::sqrDistance2D( midPoint, inter );
250  if ( dist < minDistance )
251  {
252  minDistance = dist;
253  minPoint = inter;
254  }
255  }
256  }
257  }
258 
259  return minPoint;
260 }
261 
262 QgsSnapIndex::SnapItem *QgsSnapIndex::getSnapItem( const QgsPoint &pos, const double tolerance, QgsSnapIndex::PointSnapItem **pSnapPoint, QgsSnapIndex::SegmentSnapItem **pSnapSegment, bool endPointOnly ) const
263 {
264  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
265 
266  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 2, 2 );
267 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
268  GEOSCoordSeq_setXY_r( geosctxt, coord, 0, pos.x() - tolerance, pos.y() - tolerance );
269  GEOSCoordSeq_setXY_r( geosctxt, coord, 1, pos.x() + tolerance, pos.y() + tolerance );
270 #else
271  GEOSCoordSeq_setX_r( geosctxt, coord, 0, pos.x() - tolerance );
272  GEOSCoordSeq_setY_r( geosctxt, coord, 0, pos.y() - tolerance );
273  GEOSCoordSeq_setX_r( geosctxt, coord, 1, pos.x() + tolerance );
274  GEOSCoordSeq_setY_r( geosctxt, coord, 1, pos.y() + tolerance );
275 #endif
276 
277  geos::unique_ptr searchDiagonal( GEOSGeom_createLineString_r( geosctxt, coord ) );
278 
279  QList<SnapItem *> items;
280  struct _GEOSQueryCallbackData callbackData;
281  callbackData.list = &items;
282  GEOSSTRtree_query_r( geosctxt, mSTRTree, searchDiagonal.get(), _GEOSQueryCallback, &callbackData );
283 
284  double minDistSegment = std::numeric_limits<double>::max();
285  double minDistPoint = std::numeric_limits<double>::max();
286  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
287  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
288 
289  const auto constItems = items;
290  for ( QgsSnapIndex::SnapItem *item : constItems )
291  {
292  if ( ( ! endPointOnly && item->type == SnapPoint ) || item->type == SnapEndPoint )
293  {
294  const double dist = QgsGeometryUtils::sqrDistance2D( item->getSnapPoint( pos ), pos );
295  if ( dist < minDistPoint )
296  {
297  minDistPoint = dist;
298  snapPoint = static_cast<PointSnapItem *>( item );
299  }
300  }
301  else if ( item->type == SnapSegment && !endPointOnly )
302  {
303  if ( !static_cast<SegmentSnapItem *>( item )->withinDistance( pos, tolerance ) )
304  continue;
305 
306  QgsPoint pProj;
307  if ( !static_cast<SegmentSnapItem *>( item )->getProjection( pos, pProj ) )
308  continue;
309 
310  const double dist = QgsGeometryUtils::sqrDistance2D( pProj, pos );
311  if ( dist < minDistSegment )
312  {
313  minDistSegment = dist;
314  snapSegment = static_cast<SegmentSnapItem *>( item );
315  }
316  }
317  }
318  snapPoint = minDistPoint < tolerance * tolerance ? snapPoint : nullptr;
319  snapSegment = minDistSegment < tolerance * tolerance ? snapSegment : nullptr;
320  if ( pSnapPoint ) *pSnapPoint = snapPoint;
321  if ( pSnapSegment ) *pSnapSegment = snapSegment;
322  return minDistPoint < minDistSegment ? static_cast<QgsSnapIndex::SnapItem *>( snapPoint ) : static_cast<QgsSnapIndex::SnapItem *>( snapSegment );
323 }
324 
326 
327 
328 //
329 // QgsGeometrySnapper
330 //
331 
333  : mReferenceSource( referenceSource )
334 {
335  // Build spatial index
336  mIndex = QgsSpatialIndex( *mReferenceSource );
337 }
338 
339 QgsFeatureList QgsGeometrySnapper::snapFeatures( const QgsFeatureList &features, double snapTolerance, SnapMode mode )
340 {
341  QgsFeatureList list = features;
342  QtConcurrent::blockingMap( list, ProcessFeatureWrapper( this, snapTolerance, mode ) );
343  return list;
344 }
345 
346 void QgsGeometrySnapper::processFeature( QgsFeature &feature, double snapTolerance, SnapMode mode )
347 {
348  if ( !feature.geometry().isNull() )
349  feature.setGeometry( snapGeometry( feature.geometry(), snapTolerance, mode ) );
350  emit featureSnapped();
351 }
352 
353 QgsGeometry QgsGeometrySnapper::snapGeometry( const QgsGeometry &geometry, double snapTolerance, SnapMode mode ) const
354 {
355  // Get potential reference features and construct snap index
356  QList<QgsGeometry> refGeometries;
357  mIndexMutex.lock();
358  QgsRectangle searchBounds = geometry.boundingBox();
359  searchBounds.grow( snapTolerance );
360  const QgsFeatureIds refFeatureIds = qgis::listToSet( mIndex.intersects( searchBounds ) );
361  mIndexMutex.unlock();
362 
363  if ( refFeatureIds.isEmpty() )
364  return QgsGeometry( geometry );
365 
366  refGeometries.reserve( refFeatureIds.size() );
367  QgsFeatureIds missingFeatureIds;
368  const QgsFeatureIds cachedIds = qgis::listToSet( mCachedReferenceGeometries.keys() );
369  for ( const QgsFeatureId id : refFeatureIds )
370  {
371  if ( cachedIds.contains( id ) )
372  {
373  refGeometries.append( mCachedReferenceGeometries[id] );
374  }
375  else
376  {
377  missingFeatureIds << id;
378  }
379  }
380 
381  if ( missingFeatureIds.size() > 0 )
382  {
383 
384  mReferenceLayerMutex.lock();
385  const QgsFeatureRequest refFeatureRequest = QgsFeatureRequest().setFilterFids( missingFeatureIds ).setNoAttributes();
386  QgsFeatureIterator refFeatureIt = mReferenceSource->getFeatures( refFeatureRequest );
387  QgsFeature refFeature;
388  while ( refFeatureIt.nextFeature( refFeature ) )
389  {
390  refGeometries.append( refFeature.geometry() );
391  }
392  mReferenceLayerMutex.unlock();
393  }
394 
395  return snapGeometry( geometry, snapTolerance, refGeometries, mode );
396 }
397 
398 QgsGeometry QgsGeometrySnapper::snapGeometry( const QgsGeometry &geometry, double snapTolerance, const QList<QgsGeometry> &referenceGeometries, QgsGeometrySnapper::SnapMode mode )
399 {
401  ( mode == EndPointPreferClosest || mode == EndPointPreferNodes || mode == EndPointToEndPoint ) )
402  return geometry;
403 
404  const QgsPoint center = qgsgeometry_cast< const QgsPoint * >( geometry.constGet() ) ? *static_cast< const QgsPoint * >( geometry.constGet() ) :
405  QgsPoint( geometry.constGet()->boundingBox().center() );
406 
407  QgsSnapIndex refSnapIndex;
408  for ( const QgsGeometry &geom : referenceGeometries )
409  {
410  refSnapIndex.addGeometry( geom.constGet() );
411  }
412 
413  // Snap geometries
414  QgsAbstractGeometry *subjGeom = geometry.constGet()->clone();
415  QList < QList< QList<PointFlag> > > subjPointFlags;
416 
417  // Pass 1: snap vertices of subject geometry to reference vertices
418  for ( int iPart = 0, nParts = subjGeom->partCount(); iPart < nParts; ++iPart )
419  {
420  subjPointFlags.append( QList< QList<PointFlag> >() );
421 
422  for ( int iRing = 0, nRings = subjGeom->ringCount( iPart ); iRing < nRings; ++iRing )
423  {
424  subjPointFlags[iPart].append( QList<PointFlag>() );
425 
426  for ( int iVert = 0, nVerts = polyLineSize( subjGeom, iPart, iRing ); iVert < nVerts; ++iVert )
427  {
428  if ( ( mode == EndPointPreferClosest || mode == EndPointPreferNodes || mode == EndPointToEndPoint ) &&
429  QgsWkbTypes::geometryType( subjGeom->wkbType() ) == QgsWkbTypes::LineGeometry && ( iVert > 0 && iVert < nVerts - 1 ) )
430  {
431  //endpoint mode and not at an endpoint, skip
432  subjPointFlags[iPart][iRing].append( Unsnapped );
433  continue;
434  }
435 
436  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
437  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
438  const QgsVertexId vidx( iPart, iRing, iVert );
439  const QgsPoint p = subjGeom->vertexAt( vidx );
440  if ( !refSnapIndex.getSnapItem( p, snapTolerance, &snapPoint, &snapSegment, mode == EndPointToEndPoint ) )
441  {
442  subjPointFlags[iPart][iRing].append( Unsnapped );
443  }
444  else
445  {
446  switch ( mode )
447  {
448  case PreferNodes:
450  case EndPointPreferNodes:
451  case EndPointToEndPoint:
452  {
453  // Prefer snapping to point
454  if ( snapPoint )
455  {
456  subjGeom->moveVertex( vidx, snapPoint->getSnapPoint( p ) );
457  subjPointFlags[iPart][iRing].append( SnappedToRefNode );
458  }
459  else if ( snapSegment )
460  {
461  subjGeom->moveVertex( vidx, snapSegment->getSnapPoint( p ) );
462  subjPointFlags[iPart][iRing].append( SnappedToRefSegment );
463  }
464  break;
465  }
466 
467  case PreferClosest:
470  {
471  QgsPoint nodeSnap, segmentSnap;
472  double distanceNode = std::numeric_limits<double>::max();
473  double distanceSegment = std::numeric_limits<double>::max();
474  if ( snapPoint )
475  {
476  nodeSnap = snapPoint->getSnapPoint( p );
477  distanceNode = nodeSnap.distanceSquared( p );
478  }
479  if ( snapSegment )
480  {
481  segmentSnap = snapSegment->getSnapPoint( p );
482  distanceSegment = segmentSnap.distanceSquared( p );
483  }
484  if ( snapPoint && distanceNode < distanceSegment )
485  {
486  subjGeom->moveVertex( vidx, nodeSnap );
487  subjPointFlags[iPart][iRing].append( SnappedToRefNode );
488  }
489  else if ( snapSegment )
490  {
491  subjGeom->moveVertex( vidx, segmentSnap );
492  subjPointFlags[iPart][iRing].append( SnappedToRefSegment );
493  }
494  break;
495  }
496  }
497  }
498  }
499  }
500  }
501 
502  // no extra vertices to add for point geometry
503  if ( qgsgeometry_cast< const QgsPoint * >( subjGeom ) )
504  return QgsGeometry( subjGeom );
505 
506  // nor for no extra vertices modes and end point only snapping
508  {
509  QgsGeometry result( subjGeom );
510  result.removeDuplicateNodes();
511  return result;
512  }
513 
514  std::unique_ptr< QgsSnapIndex > subjSnapIndex( new QgsSnapIndex() );
515  subjSnapIndex->addGeometry( subjGeom );
516 
517  std::unique_ptr< QgsAbstractGeometry > origSubjGeom( subjGeom->clone() );
518  std::unique_ptr< QgsSnapIndex > origSubjSnapIndex( new QgsSnapIndex() );
519  origSubjSnapIndex->addGeometry( origSubjGeom.get() );
520 
521  // Pass 2: add missing vertices to subject geometry
522  for ( const QgsGeometry &refGeom : referenceGeometries )
523  {
524  for ( int iPart = 0, nParts = refGeom.constGet()->partCount(); iPart < nParts; ++iPart )
525  {
526  for ( int iRing = 0, nRings = refGeom.constGet()->ringCount( iPart ); iRing < nRings; ++iRing )
527  {
528  for ( int iVert = 0, nVerts = polyLineSize( refGeom.constGet(), iPart, iRing ); iVert < nVerts; ++iVert )
529  {
530  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
531  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
532  const QgsPoint point = refGeom.constGet()->vertexAt( QgsVertexId( iPart, iRing, iVert ) );
533  if ( subjSnapIndex->getSnapItem( point, snapTolerance, &snapPoint, &snapSegment ) )
534  {
535  // Snap to segment, unless a subject point was already snapped to the reference point
536  if ( snapPoint )
537  {
538  const QgsPoint snappedPoint = snapPoint->getSnapPoint( point );
539  if ( QgsGeometryUtils::sqrDistance2D( snappedPoint, point ) < 1E-16 )
540  continue;
541  }
542 
543  if ( snapSegment )
544  {
545  // Look if there is a closer reference segment, if so, ignore this point
546  const QgsPoint pProj = snapSegment->getSnapPoint( point );
547  const QgsPoint closest = refSnapIndex.getClosestSnapToPoint( point, pProj );
548  if ( QgsGeometryUtils::sqrDistance2D( pProj, point ) > QgsGeometryUtils::sqrDistance2D( pProj, closest ) )
549  {
550  continue;
551  }
552 
553  // If we are too far away from the original geometry, do nothing
554  if ( !origSubjSnapIndex->getSnapItem( point, snapTolerance ) )
555  {
556  continue;
557  }
558 
559  const QgsSnapIndex::CoordIdx *idx = snapSegment->idxFrom;
560  subjGeom->insertVertex( QgsVertexId( idx->vidx.part, idx->vidx.ring, idx->vidx.vertex + 1 ), point );
561  subjPointFlags[idx->vidx.part][idx->vidx.ring].insert( idx->vidx.vertex + 1, SnappedToRefNode );
562  subjSnapIndex.reset( new QgsSnapIndex() );
563  subjSnapIndex->addGeometry( subjGeom );
564  }
565  }
566  }
567  }
568  }
569  }
570  subjSnapIndex.reset();
571  origSubjSnapIndex.reset();
572  origSubjGeom.reset();
573 
574  // Pass 3: remove superfluous vertices: all vertices which are snapped to a segment and not preceded or succeeded by an unsnapped vertex
575  for ( int iPart = 0, nParts = subjGeom->partCount(); iPart < nParts; ++iPart )
576  {
577  for ( int iRing = 0, nRings = subjGeom->ringCount( iPart ); iRing < nRings; ++iRing )
578  {
579  const bool ringIsClosed = subjGeom->vertexAt( QgsVertexId( iPart, iRing, 0 ) ) == subjGeom->vertexAt( QgsVertexId( iPart, iRing, subjGeom->vertexCount( iPart, iRing ) - 1 ) );
580  for ( int iVert = 0, nVerts = polyLineSize( subjGeom, iPart, iRing ); iVert < nVerts; ++iVert )
581  {
582  const int iPrev = ( iVert - 1 + nVerts ) % nVerts;
583  const int iNext = ( iVert + 1 ) % nVerts;
584  const QgsPoint pMid = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iVert ) );
585  const QgsPoint pPrev = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iPrev ) );
586  const QgsPoint pNext = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iNext ) );
587 
588  if ( subjPointFlags[iPart][iRing][iVert] == SnappedToRefSegment &&
589  subjPointFlags[iPart][iRing][iPrev] != Unsnapped &&
590  subjPointFlags[iPart][iRing][iNext] != Unsnapped &&
591  QgsGeometryUtils::sqrDistance2D( QgsGeometryUtils::projectPointOnSegment( pMid, pPrev, pNext ), pMid ) < 1E-12 )
592  {
593  if ( ( ringIsClosed && nVerts > 3 ) || ( !ringIsClosed && nVerts > 2 ) )
594  {
595  subjGeom->deleteVertex( QgsVertexId( iPart, iRing, iVert ) );
596  subjPointFlags[iPart][iRing].removeAt( iVert );
597  iVert -= 1;
598  nVerts -= 1;
599  }
600  else
601  {
602  // Don't delete vertices if this would result in a degenerate geometry
603  break;
604  }
605  }
606  }
607  }
608  }
609 
610  QgsGeometry result( subjGeom );
611  result.removeDuplicateNodes();
612  return result;
613 }
614 
615 int QgsGeometrySnapper::polyLineSize( const QgsAbstractGeometry *geom, int iPart, int iRing )
616 {
617  const int nVerts = geom->vertexCount( iPart, iRing );
618 
619  if ( qgsgeometry_cast< const QgsSurface * >( geom ) || qgsgeometry_cast< const QgsMultiSurface * >( geom ) )
620  {
621  const QgsPoint front = geom->vertexAt( QgsVertexId( iPart, iRing, 0 ) );
622  const QgsPoint back = geom->vertexAt( QgsVertexId( iPart, iRing, nVerts - 1 ) );
623  if ( front == back )
624  return nVerts - 1;
625  }
626 
627  return nVerts;
628 }
629 
630 
631 
632 
633 
634 //
635 // QgsInternalGeometrySnapper
636 //
637 
639  : mSnapTolerance( snapTolerance )
640  , mMode( mode )
641 {}
642 
644 {
645  if ( !feature.hasGeometry() )
646  return QgsGeometry();
647 
648  QgsFeature feat = feature;
649  QgsGeometry geometry = feat.geometry();
650  if ( !mFirstFeature )
651  {
652  // snap against processed geometries
653  // Get potential reference features and construct snap index
654  QgsRectangle searchBounds = geometry.boundingBox();
655  searchBounds.grow( mSnapTolerance );
656  const QgsFeatureIds refFeatureIds = qgis::listToSet( mProcessedIndex.intersects( searchBounds ) );
657  if ( !refFeatureIds.isEmpty() )
658  {
659  QList< QgsGeometry > refGeometries;
660  const auto constRefFeatureIds = refFeatureIds;
661  for ( const QgsFeatureId id : constRefFeatureIds )
662  {
663  refGeometries << mProcessedGeometries.value( id );
664  }
665 
666  geometry = QgsGeometrySnapper::snapGeometry( geometry, mSnapTolerance, refGeometries, mMode );
667  }
668  }
669  mProcessedGeometries.insert( feat.id(), geometry );
670  mProcessedIndex.addFeature( feat );
671  mFirstFeature = false;
672  return geometry;
673 }
674 
Abstract base class for all geometries.
virtual int ringCount(int part=0) const =0
Returns the number of rings of which this geometry is built.
virtual bool moveVertex(QgsVertexId position, const QgsPoint &newPos)=0
Moves a vertex within the geometry.
virtual int vertexCount(int part=0, int ring=0) const =0
Returns the number of vertices of which this geometry is built.
virtual QgsRectangle boundingBox() const =0
Returns the minimal bounding box for the geometry.
virtual QgsPoint vertexAt(QgsVertexId id) const =0
Returns the point corresponding to a specified vertex id.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
virtual bool insertVertex(QgsVertexId position, const QgsPoint &vertex)=0
Inserts a vertex into the geometry.
virtual int partCount() const =0
Returns count of parts contained in the geometry.
virtual bool deleteVertex(QgsVertexId position)=0
Deletes a vertex within the geometry.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
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 & setFilterFids(const QgsFeatureIds &fids)
Sets the feature IDs that should be fetched.
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
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.
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
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:205
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
Definition: qgsfeature.cpp:145
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
void featureSnapped()
Emitted each time a feature has been processed when calling snapFeatures()
QgsFeatureList snapFeatures(const QgsFeatureList &features, double snapTolerance, SnapMode mode=PreferNodes)
Snaps a set of features to the reference layer and returns the result.
QgsGeometry snapGeometry(const QgsGeometry &geometry, double snapTolerance, SnapMode mode=PreferNodes) const
Snaps a geometry to the reference layer and returns the result.
SnapMode
Snapping modes.
@ EndPointPreferClosest
Only snap start/end points of lines (point features will also be snapped, polygon features will not b...
@ PreferClosestNoExtraVertices
Snap to closest point, regardless of it is a node or a segment. No new nodes will be inserted.
@ EndPointPreferNodes
Only snap start/end points of lines (point features will also be snapped, polygon features will not b...
@ PreferNodes
Prefer to snap to nodes, even when a segment may be closer than a node. New nodes will be inserted to...
@ PreferClosest
Snap to closest point, regardless of it is a node or a segment. New nodes will be inserted to make ge...
@ EndPointToEndPoint
Only snap the start/end points of lines to other start/end points of lines.
@ PreferNodesNoExtraVertices
Prefer to snap to nodes, even when a segment may be closer than a node. No new nodes will be inserted...
QgsGeometrySnapper(QgsFeatureSource *referenceSource)
Constructor for QgsGeometrySnapper.
static double sqrDistance2D(const QgsPoint &pt1, const QgsPoint &pt2) SIP_HOLDGIL
Returns the squared 2D distance between two points.
static QgsPoint projectPointOnSegment(const QgsPoint &p, const QgsPoint &s1, const QgsPoint &s2) SIP_HOLDGIL
Project the point on a segment.
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:124
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Q_GADGET bool isNull
Definition: qgsgeometry.h:126
bool removeDuplicateNodes(double epsilon=4 *std::numeric_limits< double >::epsilon(), bool useZValues=false)
Removes duplicate nodes from the geometry, wherever removing the nodes does not result in a degenerat...
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:3213
QgsInternalGeometrySnapper(double snapTolerance, QgsGeometrySnapper::SnapMode mode=QgsGeometrySnapper::PreferNodes)
Constructor for QgsInternalGeometrySnapper.
QgsGeometry snapFeature(const QgsFeature &feature)
Snaps a single feature's geometry against all feature geometries already processed by calls to snapFe...
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
Q_GADGET double x
Definition: qgspoint.h:52
QgsPoint vertexAt(QgsVertexId) const override
Returns the point corresponding to a specified vertex id.
Definition: qgspoint.cpp:525
double distanceSquared(double x, double y) const SIP_HOLDGIL
Returns the Cartesian 2D squared distance between this point a specified x, y coordinate.
Definition: qgspoint.h:367
double y
Definition: qgspoint.h:53
A rectangle specified with double values.
Definition: qgsrectangle.h:42
void grow(double delta)
Grows the rectangle in place by the specified amount.
Definition: qgsrectangle.h:296
QgsPointXY center() const SIP_HOLDGIL
Returns the center point of the rectangle.
Definition: qgsrectangle.h:251
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.
A class to represent a vector.
Definition: qgsvector.h:30
static GeometryType geometryType(Type type) SIP_HOLDGIL
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:938
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition: qgsgeos.h:79
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:1234
QList< QgsFeature > QgsFeatureList
Definition: qgsfeature.h:834
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeatureid.h:37
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28
QLineF segment(int index, QRectF rect, double radius)
Utility class for identifying a unique vertex within a geometry.