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