QGIS API Documentation  3.21.0-Master (402cf1129e)
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 
28 #include <QtConcurrentMap>
29 #include <geos_c.h>
30 
32 
33 QgsSnapIndex::PointSnapItem::PointSnapItem( const QgsSnapIndex::CoordIdx *_idx, bool isEndPoint )
34  : SnapItem( isEndPoint ? QgsSnapIndex::SnapEndPoint : QgsSnapIndex::SnapPoint )
35  , idx( _idx )
36 {}
37 
38 QgsPoint QgsSnapIndex::PointSnapItem::getSnapPoint( const QgsPoint &/*p*/ ) const
39 {
40  return idx->point();
41 }
42 
43 QgsSnapIndex::SegmentSnapItem::SegmentSnapItem( const QgsSnapIndex::CoordIdx *_idxFrom, const QgsSnapIndex::CoordIdx *_idxTo )
44  : SnapItem( QgsSnapIndex::SnapSegment )
45  , idxFrom( _idxFrom )
46  , idxTo( _idxTo )
47 {}
48 
49 QgsPoint QgsSnapIndex::SegmentSnapItem::getSnapPoint( const QgsPoint &p ) const
50 {
51  return QgsGeometryUtils::projectPointOnSegment( p, idxFrom->point(), idxTo->point() );
52 }
53 
54 bool QgsSnapIndex::SegmentSnapItem::getIntersection( const QgsPoint &p1, const QgsPoint &p2, QgsPoint &inter ) const
55 {
56  const QgsPoint &q1 = idxFrom->point(), & q2 = idxTo->point();
57  QgsVector v( p2.x() - p1.x(), p2.y() - p1.y() );
58  QgsVector w( q2.x() - q1.x(), q2.y() - q1.y() );
59  double vl = v.length();
60  double wl = w.length();
61 
62  if ( qgsDoubleNear( vl, 0, 0.000000000001 ) || qgsDoubleNear( wl, 0, 0.000000000001 ) )
63  {
64  return false;
65  }
66  v = v / vl;
67  w = w / wl;
68 
69  double d = v.y() * w.x() - v.x() * w.y();
70 
71  if ( d == 0 )
72  return false;
73 
74  double dx = q1.x() - p1.x();
75  double dy = q1.y() - p1.y();
76  double k = ( dy * w.x() - dx * w.y() ) / d;
77 
78  inter = QgsPoint( p1.x() + v.x() * k, p1.y() + v.y() * k );
79 
80  double lambdav = QgsVector( inter.x() - p1.x(), inter.y() - p1.y() ) * v;
81  if ( lambdav < 0. + 1E-8 || lambdav > vl - 1E-8 )
82  return false;
83 
84  double lambdaw = QgsVector( inter.x() - q1.x(), inter.y() - q1.y() ) * w;
85  return !( lambdaw < 0. + 1E-8 || lambdaw >= wl - 1E-8 );
86 }
87 
88 bool QgsSnapIndex::SegmentSnapItem::getProjection( const QgsPoint &p, QgsPoint &pProj )
89 {
90  const QgsPoint &s1 = idxFrom->point();
91  const QgsPoint &s2 = idxTo->point();
92  double nx = s2.y() - s1.y();
93  double ny = -( s2.x() - s1.x() );
94  double t = ( p.x() * ny - p.y() * nx - s1.x() * ny + s1.y() * nx ) / ( ( s2.x() - s1.x() ) * ny - ( s2.y() - s1.y() ) * nx );
95  if ( t < 0. || t > 1. )
96  {
97  return false;
98  }
99  pProj = QgsPoint( s1.x() + ( s2.x() - s1.x() ) * t, s1.y() + ( s2.y() - s1.y() ) * t );
100  return true;
101 }
102 
104 
105 class Raytracer
106 {
107  // Raytrace on an integer, unit-width 2D grid with floating point coordinates
108  // See http://playtechs.blogspot.ch/2007/03/raytracing-on-grid.html
109  public:
110  Raytracer( float x0, float y0, float x1, float y1 )
111  : m_dx( std::fabs( x1 - x0 ) )
112  , m_dy( std::fabs( y1 - y0 ) )
113  , m_x( std::floor( x0 ) )
114  , m_y( std::floor( y0 ) )
115  , m_n( 1 )
116  {
117  if ( m_dx == 0. )
118  {
119  m_xInc = 0.;
120  m_error = std::numeric_limits<float>::infinity();
121  }
122  else if ( x1 > x0 )
123  {
124  m_xInc = 1;
125  m_n += int( std::floor( x1 ) ) - m_x;
126  m_error = ( std::floor( x0 ) + 1 - x0 ) * m_dy;
127  }
128  else
129  {
130  m_xInc = -1;
131  m_n += m_x - int( std::floor( x1 ) );
132  m_error = ( x0 - std::floor( x0 ) ) * m_dy;
133  }
134  if ( m_dy == 0. )
135  {
136  m_yInc = 0.;
137  m_error = -std::numeric_limits<float>::infinity();
138  }
139  else if ( y1 > y0 )
140  {
141  m_yInc = 1;
142  m_n += int( std::floor( y1 ) ) - m_y;
143  m_error -= ( std::floor( y0 ) + 1 - y0 ) * m_dx;
144  }
145  else
146  {
147  m_yInc = -1;
148  m_n += m_y - int( std::floor( y1 ) );
149  m_error -= ( y0 - std::floor( y0 ) ) * m_dx;
150  }
151  }
152  int curCol() const { return m_x; }
153  int curRow() const { return m_y; }
154  void next()
155  {
156  if ( m_error > 0 )
157  {
158  m_y += m_yInc;
159  m_error -= m_dx;
160  }
161  else if ( m_error < 0 )
162  {
163  m_x += m_xInc;
164  m_error += m_dy;
165  }
166  else
167  {
168  m_x += m_xInc;
169  m_y += m_yInc;
170  m_error += m_dx;
171  m_error -= m_dy;
172  --m_n;
173  }
174  --m_n;
175  }
176 
177  bool isValid() const { return m_n > 0; }
178 
179  private:
180  float m_dx, m_dy;
181  int m_x, m_y;
182  int m_xInc, m_yInc;
183  float m_error;
184  int m_n;
185 };
186 
188 
189 QgsSnapIndex::QgsSnapIndex( const QgsPoint &origin, double cellSize )
190  : mOrigin( origin )
191  , mCellSize( cellSize )
192 {
193  mSTRTree = GEOSSTRtree_create_r( QgsGeos::getGEOSHandler(), ( size_t )10 );
194 }
195 
196 QgsSnapIndex::~QgsSnapIndex()
197 {
198  qDeleteAll( mCoordIdxs );
199  qDeleteAll( mSnapItems );
200 
201  GEOSSTRtree_destroy_r( QgsGeos::getGEOSHandler(), mSTRTree );
202 }
203 
204 void QgsSnapIndex::addPoint( const CoordIdx *idx, bool isEndPoint )
205 {
206  QgsPoint p = idx->point();
207  int col = std::floor( ( p.x() - mOrigin.x() ) / mCellSize );
208  int row = std::floor( ( p.y() - mOrigin.y() ) / mCellSize );
209 
210  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
211 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
212  geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, row, col ) );
213 #else
214  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
215  GEOSCoordSeq_setX_r( geosctxt, seq, 0, row );
216  GEOSCoordSeq_setY_r( geosctxt, seq, 0, col );
217  geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
218 #endif
219 
220  PointSnapItem *item = new PointSnapItem( idx, isEndPoint );
221  GEOSSTRtree_insert_r( geosctxt, mSTRTree, point.get(), item );
222 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR<9
223  mSTRTreeItems.emplace_back( std::move( point ) );
224 #endif
225  mSnapItems << item;
226 }
227 
228 void QgsSnapIndex::addSegment( const CoordIdx *idxFrom, const CoordIdx *idxTo )
229 {
230  QgsPoint pFrom = idxFrom->point();
231  QgsPoint pTo = idxTo->point();
232  // Raytrace along the grid, get touched cells
233  float x0 = ( pFrom.x() - mOrigin.x() ) / mCellSize;
234  float y0 = ( pFrom.y() - mOrigin.y() ) / mCellSize;
235  float x1 = ( pTo.x() - mOrigin.x() ) / mCellSize;
236  float y1 = ( pTo.y() - mOrigin.y() ) / mCellSize;
237 
238  Raytracer rt( x0, y0, x1, y1 );
239  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
240  for ( ; rt.isValid(); rt.next() )
241  {
242 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
243  geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, rt.curRow(), rt.curCol() ) );
244 #else
245  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
246  GEOSCoordSeq_setX_r( geosctxt, seq, 0, rt.curRow() );
247  GEOSCoordSeq_setY_r( geosctxt, seq, 0, rt.curCol() );
248  geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
249 #endif
250 
251  SegmentSnapItem *item = new SegmentSnapItem( idxFrom, idxTo );
252  GEOSSTRtree_insert_r( geosctxt, mSTRTree, point.get(), item );
253 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR<9
254  mSTRTreeItems.push_back( std::move( point ) );
255 #endif
256  mSnapItems << item;
257  }
258 }
259 
260 void QgsSnapIndex::addGeometry( const QgsAbstractGeometry *geom )
261 {
262  for ( int iPart = 0, nParts = geom->partCount(); iPart < nParts; ++iPart )
263  {
264  for ( int iRing = 0, nRings = geom->ringCount( iPart ); iRing < nRings; ++iRing )
265  {
266  int nVerts = geom->vertexCount( iPart, iRing );
267 
268  if ( qgsgeometry_cast< const QgsSurface * >( geom ) )
269  nVerts--;
270  else if ( const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( geom ) )
271  {
272  if ( curve->isClosed() )
273  nVerts--;
274  }
275 
276  for ( int iVert = 0; iVert < nVerts; ++iVert )
277  {
278  CoordIdx *idx = new CoordIdx( geom, QgsVertexId( iPart, iRing, iVert ) );
279  CoordIdx *idx1 = new CoordIdx( geom, QgsVertexId( iPart, iRing, iVert + 1 ) );
280  mCoordIdxs.append( idx );
281  mCoordIdxs.append( idx1 );
282  addPoint( idx, iVert == 0 || iVert == nVerts - 1 );
283  if ( iVert < nVerts - 1 )
284  addSegment( idx, idx1 );
285  }
286  }
287  }
288 }
289 
290 struct _GEOSQueryCallbackData
291 {
292  QList< QgsSnapIndex::SnapItem * > *list;
293 };
294 
295 void _GEOSQueryCallback( void *item, void *userdata )
296 {
297  reinterpret_cast<_GEOSQueryCallbackData *>( userdata )->list->append( static_cast<QgsSnapIndex::SnapItem *>( item ) );
298 }
299 
300 QgsPoint QgsSnapIndex::getClosestSnapToPoint( const QgsPoint &p, const QgsPoint &q )
301 {
302  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
303 
304  // Look for intersections on segment from the target point to the point opposite to the point reference point
305  // p2 = p1 + 2 * (q - p1)
306  QgsPoint p2( 2 * q.x() - p.x(), 2 * q.y() - p.y() );
307 
308  // Raytrace along the grid, get touched cells
309  float x0 = ( p.x() - mOrigin.x() ) / mCellSize;
310  float y0 = ( p.y() - mOrigin.y() ) / mCellSize;
311  float x1 = ( p2.x() - mOrigin.x() ) / mCellSize;
312  float y1 = ( p2.y() - mOrigin.y() ) / mCellSize;
313 
314  Raytracer rt( x0, y0, x1, y1 );
315  double dMin = std::numeric_limits<double>::max();
316  QgsPoint pMin = p;
317  for ( ; rt.isValid(); rt.next() )
318  {
319 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
320  geos::unique_ptr searchPoint( GEOSGeom_createPointFromXY_r( geosctxt, rt.curRow(), rt.curCol() ) );
321 #else
322  GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
323  GEOSCoordSeq_setX_r( geosctxt, seq, 0, rt.curRow() );
324  GEOSCoordSeq_setY_r( geosctxt, seq, 0, rt.curCol() );
325  geos::unique_ptr searchPoint( GEOSGeom_createPoint_r( geosctxt, seq ) );
326 #endif
327  QList<SnapItem *> items;
328  struct _GEOSQueryCallbackData callbackData;
329  callbackData.list = &items;
330  GEOSSTRtree_query_r( geosctxt, mSTRTree, searchPoint.get(), _GEOSQueryCallback, &callbackData );
331  for ( const SnapItem *item : items )
332  {
333  if ( item->type == SnapSegment )
334  {
335  QgsPoint inter;
336  if ( static_cast<const SegmentSnapItem *>( item )->getIntersection( p, p2, inter ) )
337  {
338  double dist = QgsGeometryUtils::sqrDistance2D( q, inter );
339  if ( dist < dMin )
340  {
341  dMin = dist;
342  pMin = inter;
343  }
344  }
345  }
346  }
347  }
348 
349  return pMin;
350 }
351 
352 QgsSnapIndex::SnapItem *QgsSnapIndex::getSnapItem( const QgsPoint &pos, double tol, QgsSnapIndex::PointSnapItem **pSnapPoint, QgsSnapIndex::SegmentSnapItem **pSnapSegment, bool endPointOnly ) const
353 {
354  int colStart = std::floor( ( pos.x() - tol - mOrigin.x() ) / mCellSize );
355  int rowStart = std::floor( ( pos.y() - tol - mOrigin.y() ) / mCellSize );
356  int colEnd = std::floor( ( pos.x() + tol - mOrigin.x() ) / mCellSize );
357  int rowEnd = std::floor( ( pos.y() + tol - mOrigin.y() ) / mCellSize );
358 
359  GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
360 
361  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 2, 2 );
362 #if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
363  GEOSCoordSeq_setXY_r( geosctxt, coord, 0, rowStart, colStart );
364  GEOSCoordSeq_setXY_r( geosctxt, coord, 1, rowEnd, colEnd );
365 #else
366  GEOSCoordSeq_setX_r( geosctxt, coord, 0, rowStart );
367  GEOSCoordSeq_setY_r( geosctxt, coord, 0, colStart );
368  GEOSCoordSeq_setX_r( geosctxt, coord, 1, rowEnd );
369  GEOSCoordSeq_setY_r( geosctxt, coord, 1, colEnd );
370 #endif
371 
372  geos::unique_ptr searchDiagonal( GEOSGeom_createLineString_r( geosctxt, coord ) );
373 
374  QList<SnapItem *> items;
375  struct _GEOSQueryCallbackData callbackData;
376  callbackData.list = &items;
377  GEOSSTRtree_query_r( geosctxt, mSTRTree, searchDiagonal.get(), _GEOSQueryCallback, &callbackData );
378 
379  double minDistSegment = std::numeric_limits<double>::max();
380  double minDistPoint = std::numeric_limits<double>::max();
381  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
382  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
383 
384  const auto constItems = items;
385  for ( QgsSnapIndex::SnapItem *item : constItems )
386  {
387  if ( ( ! endPointOnly && item->type == SnapPoint ) || item->type == SnapEndPoint )
388  {
389  double dist = QgsGeometryUtils::sqrDistance2D( item->getSnapPoint( pos ), pos );
390  if ( dist < minDistPoint )
391  {
392  minDistPoint = dist;
393  snapPoint = static_cast<PointSnapItem *>( item );
394  }
395  }
396  else if ( item->type == SnapSegment && !endPointOnly )
397  {
398  QgsPoint pProj;
399  if ( !static_cast<SegmentSnapItem *>( item )->getProjection( pos, pProj ) )
400  {
401  continue;
402  }
403  double dist = QgsGeometryUtils::sqrDistance2D( pProj, pos );
404  if ( dist < minDistSegment )
405  {
406  minDistSegment = dist;
407  snapSegment = static_cast<SegmentSnapItem *>( item );
408  }
409  }
410  }
411  snapPoint = minDistPoint < tol * tol ? snapPoint : nullptr;
412  snapSegment = minDistSegment < tol * tol ? snapSegment : nullptr;
413  if ( pSnapPoint ) *pSnapPoint = snapPoint;
414  if ( pSnapSegment ) *pSnapSegment = snapSegment;
415  return minDistPoint < minDistSegment ? static_cast<QgsSnapIndex::SnapItem *>( snapPoint ) : static_cast<QgsSnapIndex::SnapItem *>( snapSegment );
416 }
417 
419 
420 
421 //
422 // QgsGeometrySnapper
423 //
424 
426  : mReferenceSource( referenceSource )
427 {
428  // Build spatial index
429  mIndex = QgsSpatialIndex( *mReferenceSource );
430 }
431 
432 QgsFeatureList QgsGeometrySnapper::snapFeatures( const QgsFeatureList &features, double snapTolerance, SnapMode mode )
433 {
434  QgsFeatureList list = features;
435  QtConcurrent::blockingMap( list, ProcessFeatureWrapper( this, snapTolerance, mode ) );
436  return list;
437 }
438 
439 void QgsGeometrySnapper::processFeature( QgsFeature &feature, double snapTolerance, SnapMode mode )
440 {
441  if ( !feature.geometry().isNull() )
442  feature.setGeometry( snapGeometry( feature.geometry(), snapTolerance, mode ) );
443  emit featureSnapped();
444 }
445 
446 QgsGeometry QgsGeometrySnapper::snapGeometry( const QgsGeometry &geometry, double snapTolerance, SnapMode mode ) const
447 {
448  // Get potential reference features and construct snap index
449  QList<QgsGeometry> refGeometries;
450  mIndexMutex.lock();
451  QgsRectangle searchBounds = geometry.boundingBox();
452  searchBounds.grow( snapTolerance );
453  QgsFeatureIds refFeatureIds = qgis::listToSet( mIndex.intersects( searchBounds ) );
454  mIndexMutex.unlock();
455 
456  if ( refFeatureIds.isEmpty() )
457  return QgsGeometry( geometry );
458 
459  refGeometries.reserve( refFeatureIds.size() );
460  QgsFeatureIds missingFeatureIds;
461  const QgsFeatureIds cachedIds = qgis::listToSet( mCachedReferenceGeometries.keys() );
462  for ( QgsFeatureId id : refFeatureIds )
463  {
464  if ( cachedIds.contains( id ) )
465  {
466  refGeometries.append( mCachedReferenceGeometries[id] );
467  }
468  else
469  {
470  missingFeatureIds << id;
471  }
472  }
473 
474  if ( missingFeatureIds.size() > 0 )
475  {
476 
477  mReferenceLayerMutex.lock();
478  QgsFeatureRequest refFeatureRequest = QgsFeatureRequest().setFilterFids( missingFeatureIds ).setNoAttributes();
479  QgsFeatureIterator refFeatureIt = mReferenceSource->getFeatures( refFeatureRequest );
480  QgsFeature refFeature;
481  while ( refFeatureIt.nextFeature( refFeature ) )
482  {
483  refGeometries.append( refFeature.geometry() );
484  }
485  mReferenceLayerMutex.unlock();
486  }
487 
488  return snapGeometry( geometry, snapTolerance, refGeometries, mode );
489 }
490 
491 QgsGeometry QgsGeometrySnapper::snapGeometry( const QgsGeometry &geometry, double snapTolerance, const QList<QgsGeometry> &referenceGeometries, QgsGeometrySnapper::SnapMode mode )
492 {
494  ( mode == EndPointPreferClosest || mode == EndPointPreferNodes || mode == EndPointToEndPoint ) )
495  return geometry;
496 
497  QgsPoint center = qgsgeometry_cast< const QgsPoint * >( geometry.constGet() ) ? *static_cast< const QgsPoint * >( geometry.constGet() ) :
498  QgsPoint( geometry.constGet()->boundingBox().center() );
499 
500  QgsSnapIndex refSnapIndex( center, 10 * snapTolerance );
501  for ( const QgsGeometry &geom : referenceGeometries )
502  {
503  refSnapIndex.addGeometry( geom.constGet() );
504  }
505 
506  // Snap geometries
507  QgsAbstractGeometry *subjGeom = geometry.constGet()->clone();
508  QList < QList< QList<PointFlag> > > subjPointFlags;
509 
510  // Pass 1: snap vertices of subject geometry to reference vertices
511  for ( int iPart = 0, nParts = subjGeom->partCount(); iPart < nParts; ++iPart )
512  {
513  subjPointFlags.append( QList< QList<PointFlag> >() );
514 
515  for ( int iRing = 0, nRings = subjGeom->ringCount( iPart ); iRing < nRings; ++iRing )
516  {
517  subjPointFlags[iPart].append( QList<PointFlag>() );
518 
519  for ( int iVert = 0, nVerts = polyLineSize( subjGeom, iPart, iRing ); iVert < nVerts; ++iVert )
520  {
521  if ( ( mode == EndPointPreferClosest || mode == EndPointPreferNodes || mode == EndPointToEndPoint ) &&
522  QgsWkbTypes::geometryType( subjGeom->wkbType() ) == QgsWkbTypes::LineGeometry && ( iVert > 0 && iVert < nVerts - 1 ) )
523  {
524  //endpoint mode and not at an endpoint, skip
525  subjPointFlags[iPart][iRing].append( Unsnapped );
526  continue;
527  }
528 
529  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
530  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
531  QgsVertexId vidx( iPart, iRing, iVert );
532  QgsPoint p = subjGeom->vertexAt( vidx );
533  if ( !refSnapIndex.getSnapItem( p, snapTolerance, &snapPoint, &snapSegment, mode == EndPointToEndPoint ) )
534  {
535  subjPointFlags[iPart][iRing].append( Unsnapped );
536  }
537  else
538  {
539  switch ( mode )
540  {
541  case PreferNodes:
543  case EndPointPreferNodes:
544  case EndPointToEndPoint:
545  {
546  // Prefer snapping to point
547  if ( snapPoint )
548  {
549  subjGeom->moveVertex( vidx, snapPoint->getSnapPoint( p ) );
550  subjPointFlags[iPart][iRing].append( SnappedToRefNode );
551  }
552  else if ( snapSegment )
553  {
554  subjGeom->moveVertex( vidx, snapSegment->getSnapPoint( p ) );
555  subjPointFlags[iPart][iRing].append( SnappedToRefSegment );
556  }
557  break;
558  }
559 
560  case PreferClosest:
563  {
564  QgsPoint nodeSnap, segmentSnap;
565  double distanceNode = std::numeric_limits<double>::max();
566  double distanceSegment = std::numeric_limits<double>::max();
567  if ( snapPoint )
568  {
569  nodeSnap = snapPoint->getSnapPoint( p );
570  distanceNode = nodeSnap.distanceSquared( p );
571  }
572  if ( snapSegment )
573  {
574  segmentSnap = snapSegment->getSnapPoint( p );
575  distanceSegment = segmentSnap.distanceSquared( p );
576  }
577  if ( snapPoint && distanceNode < distanceSegment )
578  {
579  subjGeom->moveVertex( vidx, nodeSnap );
580  subjPointFlags[iPart][iRing].append( SnappedToRefNode );
581  }
582  else if ( snapSegment )
583  {
584  subjGeom->moveVertex( vidx, segmentSnap );
585  subjPointFlags[iPart][iRing].append( SnappedToRefSegment );
586  }
587  break;
588  }
589  }
590  }
591  }
592  }
593  }
594 
595  //nothing more to do for points
596  if ( qgsgeometry_cast< const QgsPoint * >( subjGeom ) )
597  return QgsGeometry( subjGeom );
598  //or for end point snapping
600  {
601  QgsGeometry result( subjGeom );
602  result.removeDuplicateNodes();
603  return result;
604  }
605 
606  // SnapIndex for subject feature
607  std::unique_ptr< QgsSnapIndex > subjSnapIndex( new QgsSnapIndex( center, 10 * snapTolerance ) );
608  subjSnapIndex->addGeometry( subjGeom );
609 
610  std::unique_ptr< QgsAbstractGeometry > origSubjGeom( subjGeom->clone() );
611  std::unique_ptr< QgsSnapIndex > origSubjSnapIndex( new QgsSnapIndex( center, 10 * snapTolerance ) );
612  origSubjSnapIndex->addGeometry( origSubjGeom.get() );
613 
614  // Pass 2: add missing vertices to subject geometry
615  for ( const QgsGeometry &refGeom : referenceGeometries )
616  {
617  for ( int iPart = 0, nParts = refGeom.constGet()->partCount(); iPart < nParts; ++iPart )
618  {
619  for ( int iRing = 0, nRings = refGeom.constGet()->ringCount( iPart ); iRing < nRings; ++iRing )
620  {
621  for ( int iVert = 0, nVerts = polyLineSize( refGeom.constGet(), iPart, iRing ); iVert < nVerts; ++iVert )
622  {
623 
624  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
625  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
626  QgsPoint point = refGeom.constGet()->vertexAt( QgsVertexId( iPart, iRing, iVert ) );
627  if ( subjSnapIndex->getSnapItem( point, snapTolerance, &snapPoint, &snapSegment ) )
628  {
629  // Snap to segment, unless a subject point was already snapped to the reference point
630  if ( snapPoint && QgsGeometryUtils::sqrDistance2D( snapPoint->getSnapPoint( point ), point ) < 1E-16 )
631  {
632  continue;
633  }
634  else if ( snapSegment )
635  {
636  // Look if there is a closer reference segment, if so, ignore this point
637  QgsPoint pProj = snapSegment->getSnapPoint( point );
638  QgsPoint closest = refSnapIndex.getClosestSnapToPoint( point, pProj );
639  if ( QgsGeometryUtils::sqrDistance2D( pProj, point ) > QgsGeometryUtils::sqrDistance2D( pProj, closest ) )
640  {
641  continue;
642  }
643 
644  // If we are too far away from the original geometry, do nothing
645  if ( !origSubjSnapIndex->getSnapItem( point, snapTolerance ) )
646  {
647  continue;
648  }
649 
650  const QgsSnapIndex::CoordIdx *idx = snapSegment->idxFrom;
651  subjGeom->insertVertex( QgsVertexId( idx->vidx.part, idx->vidx.ring, idx->vidx.vertex + 1 ), point );
652  subjPointFlags[idx->vidx.part][idx->vidx.ring].insert( idx->vidx.vertex + 1, SnappedToRefNode );
653  subjSnapIndex.reset( new QgsSnapIndex( center, 10 * snapTolerance ) );
654  subjSnapIndex->addGeometry( subjGeom );
655  }
656  }
657  }
658  }
659  }
660  }
661  subjSnapIndex.reset();
662  origSubjSnapIndex.reset();
663  origSubjGeom.reset();
664 
665  // Pass 3: remove superfluous vertices: all vertices which are snapped to a segment and not preceded or succeeded by an unsnapped vertex
666  for ( int iPart = 0, nParts = subjGeom->partCount(); iPart < nParts; ++iPart )
667  {
668  for ( int iRing = 0, nRings = subjGeom->ringCount( iPart ); iRing < nRings; ++iRing )
669  {
670  bool ringIsClosed = subjGeom->vertexAt( QgsVertexId( iPart, iRing, 0 ) ) == subjGeom->vertexAt( QgsVertexId( iPart, iRing, subjGeom->vertexCount( iPart, iRing ) - 1 ) );
671  for ( int iVert = 0, nVerts = polyLineSize( subjGeom, iPart, iRing ); iVert < nVerts; ++iVert )
672  {
673  int iPrev = ( iVert - 1 + nVerts ) % nVerts;
674  int iNext = ( iVert + 1 ) % nVerts;
675  QgsPoint pMid = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iVert ) );
676  QgsPoint pPrev = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iPrev ) );
677  QgsPoint pNext = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iNext ) );
678 
679  if ( subjPointFlags[iPart][iRing][iVert] == SnappedToRefSegment &&
680  subjPointFlags[iPart][iRing][iPrev] != Unsnapped &&
681  subjPointFlags[iPart][iRing][iNext] != Unsnapped &&
682  QgsGeometryUtils::sqrDistance2D( QgsGeometryUtils::projectPointOnSegment( pMid, pPrev, pNext ), pMid ) < 1E-12 )
683  {
684  if ( ( ringIsClosed && nVerts > 3 ) || ( !ringIsClosed && nVerts > 2 ) )
685  {
686  subjGeom->deleteVertex( QgsVertexId( iPart, iRing, iVert ) );
687  subjPointFlags[iPart][iRing].removeAt( iVert );
688  iVert -= 1;
689  nVerts -= 1;
690  }
691  else
692  {
693  // Don't delete vertices if this would result in a degenerate geometry
694  break;
695  }
696  }
697  }
698  }
699  }
700 
701  QgsGeometry result( subjGeom );
702  result.removeDuplicateNodes();
703  return result;
704 }
705 
706 int QgsGeometrySnapper::polyLineSize( const QgsAbstractGeometry *geom, int iPart, int iRing )
707 {
708  int nVerts = geom->vertexCount( iPart, iRing );
709 
710  if ( qgsgeometry_cast< const QgsSurface * >( geom ) || qgsgeometry_cast< const QgsMultiSurface * >( geom ) )
711  {
712  QgsPoint front = geom->vertexAt( QgsVertexId( iPart, iRing, 0 ) );
713  QgsPoint back = geom->vertexAt( QgsVertexId( iPart, iRing, nVerts - 1 ) );
714  if ( front == back )
715  return nVerts - 1;
716  }
717 
718  return nVerts;
719 }
720 
721 
722 
723 
724 
725 //
726 // QgsInternalGeometrySnapper
727 //
728 
730  : mSnapTolerance( snapTolerance )
731  , mMode( mode )
732 {}
733 
735 {
736  if ( !feature.hasGeometry() )
737  return QgsGeometry();
738 
739  QgsFeature feat = feature;
740  QgsGeometry geometry = feat.geometry();
741  if ( !mFirstFeature )
742  {
743  // snap against processed geometries
744  // Get potential reference features and construct snap index
745  QgsRectangle searchBounds = geometry.boundingBox();
746  searchBounds.grow( mSnapTolerance );
747  QgsFeatureIds refFeatureIds = qgis::listToSet( mProcessedIndex.intersects( searchBounds ) );
748  if ( !refFeatureIds.isEmpty() )
749  {
750  QList< QgsGeometry > refGeometries;
751  const auto constRefFeatureIds = refFeatureIds;
752  for ( QgsFeatureId id : constRefFeatureIds )
753  {
754  refGeometries << mProcessedGeometries.value( id );
755  }
756 
757  geometry = QgsGeometrySnapper::snapGeometry( geometry, mSnapTolerance, refGeometries, mMode );
758  }
759  }
760  mProcessedGeometries.insert( feat.id(), geometry );
761  mProcessedIndex.addFeature( feat );
762  mFirstFeature = false;
763  return geometry;
764 }
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 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.
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:3166
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:763
QList< QgsFeature > QgsFeatureList
Definition: qgsfeature.h:831
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
Utility class for identifying a unique vertex within a geometry.