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