QGIS API Documentation  3.10.0-A Coruña (6c816b4204)
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 <QtConcurrentMap>
18 #include "qgsfeatureiterator.h"
19 #include "qgsgeometry.h"
20 #include "qgsvectorlayer.h"
21 #include "qgsgeometrysnapper.h"
22 #include "qgsvectordataprovider.h"
23 #include "qgsgeometryutils.h"
24 #include "qgsmapsettings.h"
25 #include "qgssurface.h"
26 #include "qgsmultisurface.h"
27 #include "qgscurve.h"
28 
30 
31 QgsSnapIndex::PointSnapItem::PointSnapItem( const QgsSnapIndex::CoordIdx *_idx, bool isEndPoint )
32  : SnapItem( isEndPoint ? QgsSnapIndex::SnapEndPoint : QgsSnapIndex::SnapPoint )
33  , idx( _idx )
34 {}
35 
36 QgsPoint QgsSnapIndex::PointSnapItem::getSnapPoint( const QgsPoint &/*p*/ ) const
37 {
38  return idx->point();
39 }
40 
41 QgsSnapIndex::SegmentSnapItem::SegmentSnapItem( const QgsSnapIndex::CoordIdx *_idxFrom, const QgsSnapIndex::CoordIdx *_idxTo )
42  : SnapItem( QgsSnapIndex::SnapSegment )
43  , idxFrom( _idxFrom )
44  , idxTo( _idxTo )
45 {}
46 
47 QgsPoint QgsSnapIndex::SegmentSnapItem::getSnapPoint( const QgsPoint &p ) const
48 {
49  return QgsGeometryUtils::projectPointOnSegment( p, idxFrom->point(), idxTo->point() );
50 }
51 
52 bool QgsSnapIndex::SegmentSnapItem::getIntersection( const QgsPoint &p1, const QgsPoint &p2, QgsPoint &inter ) const
53 {
54  const QgsPoint &q1 = idxFrom->point(), & q2 = idxTo->point();
55  QgsVector v( p2.x() - p1.x(), p2.y() - p1.y() );
56  QgsVector w( q2.x() - q1.x(), q2.y() - q1.y() );
57  double vl = v.length();
58  double wl = w.length();
59 
60  if ( qgsDoubleNear( vl, 0, 0.000000000001 ) || qgsDoubleNear( wl, 0, 0.000000000001 ) )
61  {
62  return false;
63  }
64  v = v / vl;
65  w = w / wl;
66 
67  double d = v.y() * w.x() - v.x() * w.y();
68 
69  if ( d == 0 )
70  return false;
71 
72  double dx = q1.x() - p1.x();
73  double dy = q1.y() - p1.y();
74  double k = ( dy * w.x() - dx * w.y() ) / d;
75 
76  inter = QgsPoint( p1.x() + v.x() * k, p1.y() + v.y() * k );
77 
78  double lambdav = QgsVector( inter.x() - p1.x(), inter.y() - p1.y() ) * v;
79  if ( lambdav < 0. + 1E-8 || lambdav > vl - 1E-8 )
80  return false;
81 
82  double lambdaw = QgsVector( inter.x() - q1.x(), inter.y() - q1.y() ) * w;
83  return !( lambdaw < 0. + 1E-8 || lambdaw >= wl - 1E-8 );
84 }
85 
86 bool QgsSnapIndex::SegmentSnapItem::getProjection( const QgsPoint &p, QgsPoint &pProj )
87 {
88  const QgsPoint &s1 = idxFrom->point();
89  const QgsPoint &s2 = idxTo->point();
90  double nx = s2.y() - s1.y();
91  double ny = -( s2.x() - s1.x() );
92  double t = ( p.x() * ny - p.y() * nx - s1.x() * ny + s1.y() * nx ) / ( ( s2.x() - s1.x() ) * ny - ( s2.y() - s1.y() ) * nx );
93  if ( t < 0. || t > 1. )
94  {
95  return false;
96  }
97  pProj = QgsPoint( s1.x() + ( s2.x() - s1.x() ) * t, s1.y() + ( s2.y() - s1.y() ) * t );
98  return true;
99 }
100 
102 
103 class Raytracer
104 {
105  // Raytrace on an integer, unit-width 2D grid with floating point coordinates
106  // See http://playtechs.blogspot.ch/2007/03/raytracing-on-grid.html
107  public:
108  Raytracer( float x0, float y0, float x1, float y1 )
109  : m_dx( std::fabs( x1 - x0 ) )
110  , m_dy( std::fabs( y1 - y0 ) )
111  , m_x( std::floor( x0 ) )
112  , m_y( std::floor( y0 ) )
113  , m_n( 1 )
114  {
115  if ( m_dx == 0. )
116  {
117  m_xInc = 0.;
118  m_error = std::numeric_limits<float>::infinity();
119  }
120  else if ( x1 > x0 )
121  {
122  m_xInc = 1;
123  m_n += int( std::floor( x1 ) ) - m_x;
124  m_error = ( std::floor( x0 ) + 1 - x0 ) * m_dy;
125  }
126  else
127  {
128  m_xInc = -1;
129  m_n += m_x - int( std::floor( x1 ) );
130  m_error = ( x0 - std::floor( x0 ) ) * m_dy;
131  }
132  if ( m_dy == 0. )
133  {
134  m_yInc = 0.;
135  m_error = -std::numeric_limits<float>::infinity();
136  }
137  else if ( y1 > y0 )
138  {
139  m_yInc = 1;
140  m_n += int( std::floor( y1 ) ) - m_y;
141  m_error -= ( std::floor( y0 ) + 1 - y0 ) * m_dx;
142  }
143  else
144  {
145  m_yInc = -1;
146  m_n += m_y - int( std::floor( y1 ) );
147  m_error -= ( y0 - std::floor( y0 ) ) * m_dx;
148  }
149  }
150  int curCol() const { return m_x; }
151  int curRow() const { return m_y; }
152  void next()
153  {
154  if ( m_error > 0 )
155  {
156  m_y += m_yInc;
157  m_error -= m_dx;
158  }
159  else if ( m_error < 0 )
160  {
161  m_x += m_xInc;
162  m_error += m_dy;
163  }
164  else
165  {
166  m_x += m_xInc;
167  m_y += m_yInc;
168  m_error += m_dx;
169  m_error -= m_dy;
170  --m_n;
171  }
172  --m_n;
173  }
174 
175  bool isValid() const { return m_n > 0; }
176 
177  private:
178  float m_dx, m_dy;
179  int m_x, m_y;
180  int m_xInc, m_yInc;
181  float m_error;
182  int m_n;
183 };
184 
186 
187 QgsSnapIndex::GridRow::~GridRow()
188 {
189  const auto constMCells = mCells;
190  for ( const QgsSnapIndex::Cell &cell : constMCells )
191  {
192  qDeleteAll( cell );
193  }
194 }
195 
196 QgsSnapIndex::Cell &QgsSnapIndex::GridRow::getCreateCell( int col )
197 {
198  if ( col < mColStartIdx )
199  {
200  for ( int i = col; i < mColStartIdx; ++i )
201  {
202  mCells.prepend( Cell() );
203  }
204  mColStartIdx = col;
205  return mCells.front();
206  }
207  else if ( col >= mColStartIdx + mCells.size() )
208  {
209  for ( int i = mColStartIdx + mCells.size(); i <= col; ++i )
210  {
211  mCells.append( Cell() );
212  }
213  return mCells.back();
214  }
215  else
216  {
217  return mCells[col - mColStartIdx];
218  }
219 }
220 
221 const QgsSnapIndex::Cell *QgsSnapIndex::GridRow::getCell( int col ) const
222 {
223  if ( col < mColStartIdx || col >= mColStartIdx + mCells.size() )
224  {
225  return nullptr;
226  }
227  else
228  {
229  return &mCells[col - mColStartIdx];
230  }
231 }
232 
233 QList<QgsSnapIndex::SnapItem *> QgsSnapIndex::GridRow::getSnapItems( int colStart, int colEnd ) const
234 {
235  colStart = std::max( colStart, mColStartIdx );
236  colEnd = std::min( colEnd, mColStartIdx + mCells.size() - 1 );
237 
238  QList<SnapItem *> items;
239 
240  for ( int col = colStart; col <= colEnd; ++col )
241  {
242  items.append( mCells[col - mColStartIdx] );
243  }
244  return items;
245 }
246 
248 
249 QgsSnapIndex::QgsSnapIndex( const QgsPoint &origin, double cellSize )
250  : mOrigin( origin )
251  , mCellSize( cellSize )
252  , mRowsStartIdx( 0 )
253 {
254 }
255 
256 QgsSnapIndex::~QgsSnapIndex()
257 {
258  qDeleteAll( mCoordIdxs );
259 }
260 
261 
262 const QgsSnapIndex::Cell *QgsSnapIndex::getCell( int col, int row ) const
263 {
264  if ( row < mRowsStartIdx || row >= mRowsStartIdx + mGridRows.size() )
265  {
266  return nullptr;
267  }
268  else
269  {
270  return mGridRows[row - mRowsStartIdx].getCell( col );
271  }
272 }
273 
274 QgsSnapIndex::Cell &QgsSnapIndex::getCreateCell( int col, int row )
275 {
276  if ( row < mRowsStartIdx )
277  {
278  for ( int i = row; i < mRowsStartIdx; ++i )
279  {
280  mGridRows.prepend( GridRow() );
281  }
282  mRowsStartIdx = row;
283  return mGridRows.front().getCreateCell( col );
284  }
285  else if ( row >= mRowsStartIdx + mGridRows.size() )
286  {
287  for ( int i = mRowsStartIdx + mGridRows.size(); i <= row; ++i )
288  {
289  mGridRows.append( GridRow() );
290  }
291  return mGridRows.back().getCreateCell( col );
292  }
293  else
294  {
295  return mGridRows[row - mRowsStartIdx].getCreateCell( col );
296  }
297 }
298 
299 void QgsSnapIndex::addPoint( const CoordIdx *idx, bool isEndPoint )
300 {
301  QgsPoint p = idx->point();
302  int col = std::floor( ( p.x() - mOrigin.x() ) / mCellSize );
303  int row = std::floor( ( p.y() - mOrigin.y() ) / mCellSize );
304  getCreateCell( col, row ).append( new PointSnapItem( idx, isEndPoint ) );
305 }
306 
307 void QgsSnapIndex::addSegment( const CoordIdx *idxFrom, const CoordIdx *idxTo )
308 {
309  QgsPoint pFrom = idxFrom->point();
310  QgsPoint pTo = idxTo->point();
311  // Raytrace along the grid, get touched cells
312  float x0 = ( pFrom.x() - mOrigin.x() ) / mCellSize;
313  float y0 = ( pFrom.y() - mOrigin.y() ) / mCellSize;
314  float x1 = ( pTo.x() - mOrigin.x() ) / mCellSize;
315  float y1 = ( pTo.y() - mOrigin.y() ) / mCellSize;
316 
317  Raytracer rt( x0, y0, x1, y1 );
318  for ( ; rt.isValid(); rt.next() )
319  {
320  getCreateCell( rt.curCol(), rt.curRow() ).append( new SegmentSnapItem( idxFrom, idxTo ) );
321  }
322 }
323 
324 void QgsSnapIndex::addGeometry( const QgsAbstractGeometry *geom )
325 {
326  for ( int iPart = 0, nParts = geom->partCount(); iPart < nParts; ++iPart )
327  {
328  for ( int iRing = 0, nRings = geom->ringCount( iPart ); iRing < nRings; ++iRing )
329  {
330  int nVerts = geom->vertexCount( iPart, iRing );
331 
332  if ( qgsgeometry_cast< const QgsSurface * >( geom ) )
333  nVerts--;
334  else if ( const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( geom ) )
335  {
336  if ( curve->isClosed() )
337  nVerts--;
338  }
339 
340  for ( int iVert = 0; iVert < nVerts; ++iVert )
341  {
342  CoordIdx *idx = new CoordIdx( geom, QgsVertexId( iPart, iRing, iVert ) );
343  CoordIdx *idx1 = new CoordIdx( geom, QgsVertexId( iPart, iRing, iVert + 1 ) );
344  mCoordIdxs.append( idx );
345  mCoordIdxs.append( idx1 );
346  addPoint( idx, iVert == 0 || iVert == nVerts - 1 );
347  if ( iVert < nVerts - 1 )
348  addSegment( idx, idx1 );
349  }
350  }
351  }
352 }
353 
354 
355 QgsPoint QgsSnapIndex::getClosestSnapToPoint( const QgsPoint &p, const QgsPoint &q )
356 {
357  // Look for intersections on segment from the target point to the point opposite to the point reference point
358  // p2 = p1 + 2 * (q - p1)
359  QgsPoint p2( 2 * q.x() - p.x(), 2 * q.y() - p.y() );
360 
361  // Raytrace along the grid, get touched cells
362  float x0 = ( p.x() - mOrigin.x() ) / mCellSize;
363  float y0 = ( p.y() - mOrigin.y() ) / mCellSize;
364  float x1 = ( p2.x() - mOrigin.x() ) / mCellSize;
365  float y1 = ( p2.y() - mOrigin.y() ) / mCellSize;
366 
367  Raytracer rt( x0, y0, x1, y1 );
368  double dMin = std::numeric_limits<double>::max();
369  QgsPoint pMin = p;
370  for ( ; rt.isValid(); rt.next() )
371  {
372  const Cell *cell = getCell( rt.curCol(), rt.curRow() );
373  if ( !cell )
374  {
375  continue;
376  }
377  for ( const SnapItem *item : *cell )
378  {
379  if ( item->type == SnapSegment )
380  {
381  QgsPoint inter;
382  if ( static_cast<const SegmentSnapItem *>( item )->getIntersection( p, p2, inter ) )
383  {
384  double dist = QgsGeometryUtils::sqrDistance2D( q, inter );
385  if ( dist < dMin )
386  {
387  dMin = dist;
388  pMin = inter;
389  }
390  }
391  }
392  }
393  }
394 
395  return pMin;
396 }
397 
398 QgsSnapIndex::SnapItem *QgsSnapIndex::getSnapItem( const QgsPoint &pos, double tol, QgsSnapIndex::PointSnapItem **pSnapPoint, QgsSnapIndex::SegmentSnapItem **pSnapSegment, bool endPointOnly ) const
399 {
400  int colStart = std::floor( ( pos.x() - tol - mOrigin.x() ) / mCellSize );
401  int rowStart = std::floor( ( pos.y() - tol - mOrigin.y() ) / mCellSize );
402  int colEnd = std::floor( ( pos.x() + tol - mOrigin.x() ) / mCellSize );
403  int rowEnd = std::floor( ( pos.y() + tol - mOrigin.y() ) / mCellSize );
404 
405  rowStart = std::max( rowStart, mRowsStartIdx );
406  rowEnd = std::min( rowEnd, mRowsStartIdx + mGridRows.size() - 1 );
407 
408  QList<SnapItem *> items;
409  for ( int row = rowStart; row <= rowEnd; ++row )
410  {
411  items.append( mGridRows[row - mRowsStartIdx].getSnapItems( colStart, colEnd ) );
412  }
413 
414  double minDistSegment = std::numeric_limits<double>::max();
415  double minDistPoint = std::numeric_limits<double>::max();
416  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
417  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
418 
419  const auto constItems = items;
420  for ( QgsSnapIndex::SnapItem *item : constItems )
421  {
422  if ( ( ! endPointOnly && item->type == SnapPoint ) || item->type == SnapEndPoint )
423  {
424  double dist = QgsGeometryUtils::sqrDistance2D( item->getSnapPoint( pos ), pos );
425  if ( dist < minDistPoint )
426  {
427  minDistPoint = dist;
428  snapPoint = static_cast<PointSnapItem *>( item );
429  }
430  }
431  else if ( item->type == SnapSegment && !endPointOnly )
432  {
433  QgsPoint pProj;
434  if ( !static_cast<SegmentSnapItem *>( item )->getProjection( pos, pProj ) )
435  {
436  continue;
437  }
438  double dist = QgsGeometryUtils::sqrDistance2D( pProj, pos );
439  if ( dist < minDistSegment )
440  {
441  minDistSegment = dist;
442  snapSegment = static_cast<SegmentSnapItem *>( item );
443  }
444  }
445  }
446  snapPoint = minDistPoint < tol * tol ? snapPoint : nullptr;
447  snapSegment = minDistSegment < tol * tol ? snapSegment : nullptr;
448  if ( pSnapPoint ) *pSnapPoint = snapPoint;
449  if ( pSnapSegment ) *pSnapSegment = snapSegment;
450  return minDistPoint < minDistSegment ? static_cast<QgsSnapIndex::SnapItem *>( snapPoint ) : static_cast<QgsSnapIndex::SnapItem *>( snapSegment );
451 }
452 
454 
455 
456 //
457 // QgsGeometrySnapper
458 //
459 
461  : mReferenceSource( referenceSource )
462 {
463  // Build spatial index
464  mIndex = QgsSpatialIndex( *mReferenceSource );
465 }
466 
467 QgsFeatureList QgsGeometrySnapper::snapFeatures( const QgsFeatureList &features, double snapTolerance, SnapMode mode )
468 {
469  QgsFeatureList list = features;
470  QtConcurrent::blockingMap( list, ProcessFeatureWrapper( this, snapTolerance, mode ) );
471  return list;
472 }
473 
474 void QgsGeometrySnapper::processFeature( QgsFeature &feature, double snapTolerance, SnapMode mode )
475 {
476  if ( !feature.geometry().isNull() )
477  feature.setGeometry( snapGeometry( feature.geometry(), snapTolerance, mode ) );
478  emit featureSnapped();
479 }
480 
481 QgsGeometry QgsGeometrySnapper::snapGeometry( const QgsGeometry &geometry, double snapTolerance, SnapMode mode ) const
482 {
483  // Get potential reference features and construct snap index
484  QList<QgsGeometry> refGeometries;
485  mIndexMutex.lock();
486  QgsRectangle searchBounds = geometry.boundingBox();
487  searchBounds.grow( snapTolerance );
488  QgsFeatureIds refFeatureIds = mIndex.intersects( searchBounds ).toSet();
489  mIndexMutex.unlock();
490 
491  QgsFeatureRequest refFeatureRequest = QgsFeatureRequest().setFilterFids( refFeatureIds ).setNoAttributes();
492  mReferenceLayerMutex.lock();
493  QgsFeature refFeature;
494  QgsFeatureIterator refFeatureIt = mReferenceSource->getFeatures( refFeatureRequest );
495 
496  while ( refFeatureIt.nextFeature( refFeature ) )
497  {
498  refGeometries.append( refFeature.geometry() );
499  }
500  mReferenceLayerMutex.unlock();
501 
502  return snapGeometry( geometry, snapTolerance, refGeometries, mode );
503 }
504 
505 QgsGeometry QgsGeometrySnapper::snapGeometry( const QgsGeometry &geometry, double snapTolerance, const QList<QgsGeometry> &referenceGeometries, QgsGeometrySnapper::SnapMode mode )
506 {
508  ( mode == EndPointPreferClosest || mode == EndPointPreferNodes || mode == EndPointToEndPoint ) )
509  return geometry;
510 
511  QgsPoint center = qgsgeometry_cast< const QgsPoint * >( geometry.constGet() ) ? *static_cast< const QgsPoint * >( geometry.constGet() ) :
512  QgsPoint( geometry.constGet()->boundingBox().center() );
513 
514  QgsSnapIndex refSnapIndex( center, 10 * snapTolerance );
515  for ( const QgsGeometry &geom : referenceGeometries )
516  {
517  refSnapIndex.addGeometry( geom.constGet() );
518  }
519 
520  // Snap geometries
521  QgsAbstractGeometry *subjGeom = geometry.constGet()->clone();
522  QList < QList< QList<PointFlag> > > subjPointFlags;
523 
524  // Pass 1: snap vertices of subject geometry to reference vertices
525  for ( int iPart = 0, nParts = subjGeom->partCount(); iPart < nParts; ++iPart )
526  {
527  subjPointFlags.append( QList< QList<PointFlag> >() );
528 
529  for ( int iRing = 0, nRings = subjGeom->ringCount( iPart ); iRing < nRings; ++iRing )
530  {
531  subjPointFlags[iPart].append( QList<PointFlag>() );
532 
533  for ( int iVert = 0, nVerts = polyLineSize( subjGeom, iPart, iRing ); iVert < nVerts; ++iVert )
534  {
535  if ( ( mode == EndPointPreferClosest || mode == EndPointPreferNodes || mode == EndPointToEndPoint ) &&
536  QgsWkbTypes::geometryType( subjGeom->wkbType() ) == QgsWkbTypes::LineGeometry && ( iVert > 0 && iVert < nVerts - 1 ) )
537  {
538  //endpoint mode and not at an endpoint, skip
539  subjPointFlags[iPart][iRing].append( Unsnapped );
540  continue;
541  }
542 
543  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
544  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
545  QgsVertexId vidx( iPart, iRing, iVert );
546  QgsPoint p = subjGeom->vertexAt( vidx );
547  if ( !refSnapIndex.getSnapItem( p, snapTolerance, &snapPoint, &snapSegment, mode == EndPointToEndPoint ) )
548  {
549  subjPointFlags[iPart][iRing].append( Unsnapped );
550  }
551  else
552  {
553  switch ( mode )
554  {
555  case PreferNodes:
557  case EndPointPreferNodes:
558  case EndPointToEndPoint:
559  {
560  // Prefer snapping to point
561  if ( snapPoint )
562  {
563  subjGeom->moveVertex( vidx, snapPoint->getSnapPoint( p ) );
564  subjPointFlags[iPart][iRing].append( SnappedToRefNode );
565  }
566  else if ( snapSegment )
567  {
568  subjGeom->moveVertex( vidx, snapSegment->getSnapPoint( p ) );
569  subjPointFlags[iPart][iRing].append( SnappedToRefSegment );
570  }
571  break;
572  }
573 
574  case PreferClosest:
577  {
578  QgsPoint nodeSnap, segmentSnap;
579  double distanceNode = std::numeric_limits<double>::max();
580  double distanceSegment = std::numeric_limits<double>::max();
581  if ( snapPoint )
582  {
583  nodeSnap = snapPoint->getSnapPoint( p );
584  distanceNode = nodeSnap.distanceSquared( p );
585  }
586  if ( snapSegment )
587  {
588  segmentSnap = snapSegment->getSnapPoint( p );
589  distanceSegment = segmentSnap.distanceSquared( p );
590  }
591  if ( snapPoint && distanceNode < distanceSegment )
592  {
593  subjGeom->moveVertex( vidx, nodeSnap );
594  subjPointFlags[iPart][iRing].append( SnappedToRefNode );
595  }
596  else if ( snapSegment )
597  {
598  subjGeom->moveVertex( vidx, segmentSnap );
599  subjPointFlags[iPart][iRing].append( SnappedToRefSegment );
600  }
601  break;
602  }
603  }
604  }
605  }
606  }
607  }
608 
609  //nothing more to do for points
610  if ( qgsgeometry_cast< const QgsPoint * >( subjGeom ) )
611  return QgsGeometry( subjGeom );
612  //or for end point snapping
614  {
615  QgsGeometry result( subjGeom );
616  result.removeDuplicateNodes();
617  return result;
618  }
619 
620  // SnapIndex for subject feature
621  std::unique_ptr< QgsSnapIndex > subjSnapIndex( new QgsSnapIndex( center, 10 * snapTolerance ) );
622  subjSnapIndex->addGeometry( subjGeom );
623 
624  std::unique_ptr< QgsAbstractGeometry > origSubjGeom( subjGeom->clone() );
625  std::unique_ptr< QgsSnapIndex > origSubjSnapIndex( new QgsSnapIndex( center, 10 * snapTolerance ) );
626  origSubjSnapIndex->addGeometry( origSubjGeom.get() );
627 
628  // Pass 2: add missing vertices to subject geometry
629  for ( const QgsGeometry &refGeom : referenceGeometries )
630  {
631  for ( int iPart = 0, nParts = refGeom.constGet()->partCount(); iPart < nParts; ++iPart )
632  {
633  for ( int iRing = 0, nRings = refGeom.constGet()->ringCount( iPart ); iRing < nRings; ++iRing )
634  {
635  for ( int iVert = 0, nVerts = polyLineSize( refGeom.constGet(), iPart, iRing ); iVert < nVerts; ++iVert )
636  {
637 
638  QgsSnapIndex::PointSnapItem *snapPoint = nullptr;
639  QgsSnapIndex::SegmentSnapItem *snapSegment = nullptr;
640  QgsPoint point = refGeom.constGet()->vertexAt( QgsVertexId( iPart, iRing, iVert ) );
641  if ( subjSnapIndex->getSnapItem( point, snapTolerance, &snapPoint, &snapSegment ) )
642  {
643  // Snap to segment, unless a subject point was already snapped to the reference point
644  if ( snapPoint && QgsGeometryUtils::sqrDistance2D( snapPoint->getSnapPoint( point ), point ) < 1E-16 )
645  {
646  continue;
647  }
648  else if ( snapSegment )
649  {
650  // Look if there is a closer reference segment, if so, ignore this point
651  QgsPoint pProj = snapSegment->getSnapPoint( point );
652  QgsPoint closest = refSnapIndex.getClosestSnapToPoint( point, pProj );
653  if ( QgsGeometryUtils::sqrDistance2D( pProj, point ) > QgsGeometryUtils::sqrDistance2D( pProj, closest ) )
654  {
655  continue;
656  }
657 
658  // If we are too far away from the original geometry, do nothing
659  if ( !origSubjSnapIndex->getSnapItem( point, snapTolerance ) )
660  {
661  continue;
662  }
663 
664  const QgsSnapIndex::CoordIdx *idx = snapSegment->idxFrom;
665  subjGeom->insertVertex( QgsVertexId( idx->vidx.part, idx->vidx.ring, idx->vidx.vertex + 1 ), point );
666  subjPointFlags[idx->vidx.part][idx->vidx.ring].insert( idx->vidx.vertex + 1, SnappedToRefNode );
667  subjSnapIndex.reset( new QgsSnapIndex( center, 10 * snapTolerance ) );
668  subjSnapIndex->addGeometry( subjGeom );
669  }
670  }
671  }
672  }
673  }
674  }
675  subjSnapIndex.reset();
676  origSubjSnapIndex.reset();
677  origSubjGeom.reset();
678 
679  // Pass 3: remove superfluous vertices: all vertices which are snapped to a segment and not preceded or succeeded by an unsnapped vertex
680  for ( int iPart = 0, nParts = subjGeom->partCount(); iPart < nParts; ++iPart )
681  {
682  for ( int iRing = 0, nRings = subjGeom->ringCount( iPart ); iRing < nRings; ++iRing )
683  {
684  bool ringIsClosed = subjGeom->vertexAt( QgsVertexId( iPart, iRing, 0 ) ) == subjGeom->vertexAt( QgsVertexId( iPart, iRing, subjGeom->vertexCount( iPart, iRing ) - 1 ) );
685  for ( int iVert = 0, nVerts = polyLineSize( subjGeom, iPart, iRing ); iVert < nVerts; ++iVert )
686  {
687  int iPrev = ( iVert - 1 + nVerts ) % nVerts;
688  int iNext = ( iVert + 1 ) % nVerts;
689  QgsPoint pMid = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iVert ) );
690  QgsPoint pPrev = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iPrev ) );
691  QgsPoint pNext = subjGeom->vertexAt( QgsVertexId( iPart, iRing, iNext ) );
692 
693  if ( subjPointFlags[iPart][iRing][iVert] == SnappedToRefSegment &&
694  subjPointFlags[iPart][iRing][iPrev] != Unsnapped &&
695  subjPointFlags[iPart][iRing][iNext] != Unsnapped &&
696  QgsGeometryUtils::sqrDistance2D( QgsGeometryUtils::projectPointOnSegment( pMid, pPrev, pNext ), pMid ) < 1E-12 )
697  {
698  if ( ( ringIsClosed && nVerts > 3 ) || ( !ringIsClosed && nVerts > 2 ) )
699  {
700  subjGeom->deleteVertex( QgsVertexId( iPart, iRing, iVert ) );
701  subjPointFlags[iPart][iRing].removeAt( iVert );
702  iVert -= 1;
703  nVerts -= 1;
704  }
705  else
706  {
707  // Don't delete vertices if this would result in a degenerate geometry
708  break;
709  }
710  }
711  }
712  }
713  }
714 
715  QgsGeometry result( subjGeom );
716  result.removeDuplicateNodes();
717  return result;
718 }
719 
720 int QgsGeometrySnapper::polyLineSize( const QgsAbstractGeometry *geom, int iPart, int iRing )
721 {
722  int nVerts = geom->vertexCount( iPart, iRing );
723 
724  if ( qgsgeometry_cast< const QgsSurface * >( geom ) || qgsgeometry_cast< const QgsMultiSurface * >( geom ) )
725  {
726  QgsPoint front = geom->vertexAt( QgsVertexId( iPart, iRing, 0 ) );
727  QgsPoint back = geom->vertexAt( QgsVertexId( iPart, iRing, nVerts - 1 ) );
728  if ( front == back )
729  return nVerts - 1;
730  }
731 
732  return nVerts;
733 }
734 
735 
736 
737 
738 
739 //
740 // QgsInternalGeometrySnapper
741 //
742 
744  : mSnapTolerance( snapTolerance )
745  , mMode( mode )
746 {}
747 
749 {
750  if ( !feature.hasGeometry() )
751  return QgsGeometry();
752 
753  QgsFeature feat = feature;
754  QgsGeometry geometry = feat.geometry();
755  if ( !mFirstFeature )
756  {
757  // snap against processed geometries
758  // Get potential reference features and construct snap index
759  QgsRectangle searchBounds = geometry.boundingBox();
760  searchBounds.grow( mSnapTolerance );
761  QgsFeatureIds refFeatureIds = mProcessedIndex.intersects( searchBounds ).toSet();
762  if ( !refFeatureIds.isEmpty() )
763  {
764  QList< QgsGeometry > refGeometries;
765  const auto constRefFeatureIds = refFeatureIds;
766  for ( QgsFeatureId id : constRefFeatureIds )
767  {
768  refGeometries << mProcessedGeometries.value( id );
769  }
770 
771  geometry = QgsGeometrySnapper::snapGeometry( geometry, mSnapTolerance, refGeometries, mMode );
772  }
773  }
774  mProcessedGeometries.insert( feat.id(), geometry );
775  mProcessedIndex.addFeature( feat );
776  mFirstFeature = false;
777  return geometry;
778 }
QgsFeatureId id
Definition: qgsfeature.h:64
Wrapper for iterator of features from vector data provider or vector layer.
A rectangle specified with double values.
Definition: qgsrectangle.h:41
double y
Definition: qgspoint.h:42
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeatureid.h:34
virtual bool deleteVertex(QgsVertexId position)=0
Deletes a vertex within the geometry.
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...
virtual bool insertVertex(QgsVertexId position, const QgsPoint &vertex)=0
Inserts a vertex into the geometry.
QList< QgsFeature > QgsFeatureList
Definition: qgsfeature.h:571
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:280
Only snap the start/end points of lines to other start/end points of lines.
QgsGeometrySnapper(QgsFeatureSource *referenceSource)
Constructor for QgsGeometrySnapper.
qint64 QgsFeatureId
Definition: qgsfeatureid.h:25
void featureSnapped()
Emitted each time a feature has been processed when calling snapFeatures()
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:122
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
virtual QgsRectangle boundingBox() const =0
Returns the minimal bounding box for the geometry.
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:197
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
QgsPoint vertexAt(QgsVertexId) const override
Returns the point corresponding to a specified vertex id.
Definition: qgspoint.cpp:480
QgsFeatureList snapFeatures(const QgsFeatureList &features, double snapTolerance, SnapMode mode=PreferNodes)
Snaps a set of features to the reference layer and returns the result.
Utility class for identifying a unique vertex within a geometry.
void grow(double delta)
Grows the rectangle in place by the specified amount.
Definition: qgsrectangle.h:275
Only snap start/end points of lines (point features will also be snapped, polygon features will not b...
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
static GeometryType geometryType(Type type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:812
This class wraps a request for features to a vector layer (or directly its vector data provider)...
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
virtual int ringCount(int part=0) const =0
Returns the number of rings of which this geometry is built.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
Abstract base class for all geometries.
QgsWkbTypes::Type wkbType() const
Returns the WKB type of the geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
double length() const
Returns the length of the vector.
Definition: qgsvector.cpp:71
A class to represent a vector.
Definition: qgsvector.h:29
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsGeometry snapGeometry(const QgsGeometry &geometry, double snapTolerance, SnapMode mode=PreferNodes) const
Snaps a geometry to the reference layer and returns the result.
Prefer to snap to nodes, even when a segment may be closer than a node. New nodes will be inserted to...
static double sqrDistance2D(const QgsPoint &pt1, const QgsPoint &pt2)
Returns the squared 2D distance between two points.
A spatial index for QgsFeature objects.
QgsFeatureRequest & setFilterFids(const QgsFeatureIds &fids)
Sets feature IDs that should be fetched.
Snap to closest point, regardless of it is a node or a segment. New nodes will be inserted to make ge...
An interface for objects which provide features via a getFeatures method.
SnapMode
Snapping modes.
virtual bool moveVertex(QgsVertexId position, const QgsPoint &newPos)=0
Moves a vertex within the geometry.
QgsGeometry snapFeature(const QgsFeature &feature)
Snaps a single feature&#39;s geometry against all feature geometries already processed by calls to snapFe...
virtual int vertexCount(int part=0, int ring=0) const =0
Returns the number of vertices of which this geometry is built.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:137
static QgsPoint projectPointOnSegment(const QgsPoint &p, const QgsPoint &s1, const QgsPoint &s2)
Project the point on a segment.
QgsInternalGeometrySnapper(double snapTolerance, QgsGeometrySnapper::SnapMode mode=QgsGeometrySnapper::PreferNodes)
Constructor for QgsInternalGeometrySnapper.
Prefer to snap to nodes, even when a segment may be closer than a node. No new nodes will be inserted...
QgsGeometry geometry
Definition: qgsfeature.h:67
QgsPointXY center() const
Returns the center point of the rectangle.
Definition: qgsrectangle.h:230
Only snap start/end points of lines (point features will also be snapped, polygon features will not b...
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=nullptr) override
Adds a feature to the index.
bool nextFeature(QgsFeature &f)
double distanceSquared(double x, double y) const
Returns the Cartesian 2D squared distance between this point a specified x, y coordinate.
Definition: qgspoint.h:348
virtual QgsPoint vertexAt(QgsVertexId id) const =0
Returns the point corresponding to a specified vertex id.
Snap to closest point, regardless of it is a node or a segment. No new nodes will be inserted...
QList< QgsFeatureId > intersects(const QgsRectangle &rectangle) const
Returns a list of features with a bounding box which intersects the specified rectangle.
virtual int partCount() const =0
Returns count of parts contained in the geometry.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
double x
Definition: qgspoint.h:41