Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgspointlocator.cpp
3 : : --------------------------------------
4 : : Date : November 2014
5 : : Copyright : (C) 2014 by Martin Dobias
6 : : Email : wonder dot sk at gmail dot com
7 : : ***************************************************************************
8 : : * *
9 : : * This program is free software; you can redistribute it and/or modify *
10 : : * it under the terms of the GNU General Public License as published by *
11 : : * the Free Software Foundation; either version 2 of the License, or *
12 : : * (at your option) any later version. *
13 : : * *
14 : : ***************************************************************************/
15 : :
16 : : #include "qgspointlocator.h"
17 : :
18 : : #include "qgsfeatureiterator.h"
19 : : #include "qgsgeometry.h"
20 : : #include "qgsvectorlayer.h"
21 : : #include "qgswkbptr.h"
22 : : #include "qgis.h"
23 : : #include "qgslogger.h"
24 : : #include "qgsrenderer.h"
25 : : #include "qgsapplication.h"
26 : : #include "qgsvectorlayerfeatureiterator.h"
27 : : #include "qgsexpressioncontextutils.h"
28 : : #include "qgslinestring.h"
29 : : #include "qgscurvepolygon.h"
30 : : #include "qgspointlocatorinittask.h"
31 : : #include <spatialindex/SpatialIndex.h>
32 : :
33 : : #include <QLinkedListIterator>
34 : : #include <QtConcurrent>
35 : :
36 : : using namespace SpatialIndex;
37 : :
38 : :
39 : :
40 : 0 : static SpatialIndex::Point point2point( const QgsPointXY &point )
41 : : {
42 : 0 : double plow[2] = { point.x(), point.y() };
43 : 0 : return Point( plow, 2 );
44 : : }
45 : :
46 : :
47 : 0 : static SpatialIndex::Region rect2region( const QgsRectangle &rect )
48 : : {
49 : 0 : double pLow[2] = { rect.xMinimum(), rect.yMinimum() };
50 : 0 : double pHigh[2] = { rect.xMaximum(), rect.yMaximum() };
51 : 0 : return SpatialIndex::Region( pLow, pHigh, 2 );
52 : : }
53 : :
54 : :
55 : : // Ahh.... another magic number. Taken from QgsVectorLayer::snapToGeometry() call to closestSegmentWithContext().
56 : : // The default epsilon used for sqrDistToSegment (1e-8) is too high when working with lat/lon coordinates
57 : : // I still do not fully understand why the sqrDistToSegment() code uses epsilon and if the square distance
58 : : // is lower than epsilon it will have a special logic...
59 : : static const double POINT_LOC_EPSILON = 1e-12;
60 : :
61 : : ////////////////////////////////////////////////////////////////////////////
62 : :
63 : :
64 : : /**
65 : : * \ingroup core
66 : : * \brief Helper class for bulk loading of R-trees.
67 : : * \note not available in Python bindings
68 : : */
69 : 0 : class QgsPointLocator_Stream : public IDataStream
70 : : {
71 : : public:
72 : 0 : explicit QgsPointLocator_Stream( const QLinkedList<RTree::Data *> &dataList )
73 : 0 : : mDataList( dataList )
74 : 0 : , mIt( mDataList )
75 : 0 : { }
76 : :
77 : 0 : IData *getNext() override { return mIt.next(); }
78 : 0 : bool hasNext() override { return mIt.hasNext(); }
79 : :
80 : 0 : uint32_t size() override { Q_ASSERT( false && "not available" ); return 0; }
81 : 0 : void rewind() override { Q_ASSERT( false && "not available" ); }
82 : :
83 : : private:
84 : : QLinkedList<RTree::Data *> mDataList;
85 : : QLinkedListIterator<RTree::Data *> mIt;
86 : : };
87 : :
88 : :
89 : : ////////////////////////////////////////////////////////////////////////////
90 : :
91 : :
92 : : /**
93 : : * \ingroup core
94 : : * \brief Helper class used when traversing the index looking for vertices - builds a list of matches.
95 : : * \note not available in Python bindings
96 : : */
97 : 0 : class QgsPointLocator_VisitorNearestVertex : public IVisitor
98 : : {
99 : : public:
100 : 0 : QgsPointLocator_VisitorNearestVertex( QgsPointLocator *pl, QgsPointLocator::Match &m, const QgsPointXY &srcPoint, QgsPointLocator::MatchFilter *filter = nullptr )
101 : 0 : : mLocator( pl )
102 : 0 : , mBest( m )
103 : 0 : , mSrcPoint( srcPoint )
104 : 0 : , mFilter( filter )
105 : 0 : {}
106 : :
107 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
108 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
109 : :
110 : 0 : void visitData( const IData &d ) override
111 : : {
112 : 0 : QgsFeatureId id = d.getIdentifier();
113 : 0 : QgsGeometry *geom = mLocator->mGeoms.value( id );
114 : : int vertexIndex, beforeVertex, afterVertex;
115 : : double sqrDist;
116 : :
117 : 0 : QgsPointXY pt = geom->closestVertex( mSrcPoint, vertexIndex, beforeVertex, afterVertex, sqrDist );
118 : 0 : if ( sqrDist < 0 )
119 : 0 : return; // probably empty geometry
120 : :
121 : 0 : QgsPointLocator::Match m( QgsPointLocator::Vertex, mLocator->mLayer, id, std::sqrt( sqrDist ), pt, vertexIndex );
122 : : // in range queries the filter may reject some matches
123 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
124 : 0 : return;
125 : :
126 : 0 : if ( !mBest.isValid() || m.distance() < mBest.distance() )
127 : 0 : mBest = m;
128 : 0 : }
129 : :
130 : : private:
131 : : QgsPointLocator *mLocator = nullptr;
132 : : QgsPointLocator::Match &mBest;
133 : : QgsPointXY mSrcPoint;
134 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
135 : : };
136 : :
137 : :
138 : :
139 : : /**
140 : : * \ingroup core
141 : : * \brief Helper class used when traversing the index looking for centroid - builds a list of matches.
142 : : * \note not available in Python bindings
143 : : * \since QGIS 3.12
144 : : */
145 : 0 : class QgsPointLocator_VisitorNearestCentroid : public IVisitor
146 : : {
147 : : public:
148 : :
149 : : /**
150 : : * \ingroup core
151 : : * \brief Helper class used when traversing the index looking for centroid - builds a list of matches.
152 : : * \note not available in Python bindings
153 : : * \since QGIS 3.12
154 : : */
155 : 0 : QgsPointLocator_VisitorNearestCentroid( QgsPointLocator *pl, QgsPointLocator::Match &m, const QgsPointXY &srcPoint, QgsPointLocator::MatchFilter *filter = nullptr )
156 : 0 : : mLocator( pl )
157 : 0 : , mBest( m )
158 : 0 : , mSrcPoint( srcPoint )
159 : 0 : , mFilter( filter )
160 : 0 : {}
161 : :
162 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
163 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
164 : :
165 : 0 : void visitData( const IData &d ) override
166 : : {
167 : 0 : QgsFeatureId id = d.getIdentifier();
168 : 0 : QgsGeometry *geom = mLocator->mGeoms.value( id );
169 : :
170 : 0 : QgsPointXY pt = geom->centroid().asPoint();
171 : :
172 : 0 : QgsPointLocator::Match m( QgsPointLocator::Centroid, mLocator->mLayer, id, std::sqrt( mSrcPoint.sqrDist( pt ) ), pt, -1 );
173 : : // in range queries the filter may reject some matches
174 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
175 : 0 : return;
176 : :
177 : 0 : if ( !mBest.isValid() || m.distance() < mBest.distance() )
178 : 0 : mBest = m;
179 : :
180 : 0 : }
181 : :
182 : : private:
183 : : QgsPointLocator *mLocator = nullptr;
184 : : QgsPointLocator::Match &mBest;
185 : : QgsPointXY mSrcPoint;
186 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
187 : : };
188 : :
189 : : ////////////////////////////////////////////////////////////////////////////
190 : :
191 : : /**
192 : : * \ingroup core
193 : : * \brief Helper class used when traversing the index looking for middle segment - builds a list of matches.
194 : : * \note not available in Python bindings
195 : : * \since QGIS 3.12
196 : : */
197 : 0 : class QgsPointLocator_VisitorNearestMiddleOfSegment: public IVisitor
198 : : {
199 : : public:
200 : :
201 : : /**
202 : : * \ingroup core
203 : : * \brief Helper class used when traversing the index looking for middle segment - builds a list of matches.
204 : : * \note not available in Python bindings
205 : : * \since QGIS 3.12
206 : : */
207 : 0 : QgsPointLocator_VisitorNearestMiddleOfSegment( QgsPointLocator *pl, QgsPointLocator::Match &m, const QgsPointXY &srcPoint, QgsPointLocator::MatchFilter *filter = nullptr )
208 : 0 : : mLocator( pl )
209 : 0 : , mBest( m )
210 : 0 : , mSrcPoint( srcPoint )
211 : 0 : , mFilter( filter )
212 : 0 : {}
213 : :
214 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
215 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
216 : :
217 : 0 : void visitData( const IData &d ) override
218 : : {
219 : 0 : QgsFeatureId id = d.getIdentifier();
220 : 0 : QgsGeometry *geom = mLocator->mGeoms.value( id );
221 : 0 : QgsPointXY pt;
222 : : int afterVertex;
223 : 0 : double sqrDist = geom->closestSegmentWithContext( mSrcPoint, pt, afterVertex, nullptr, POINT_LOC_EPSILON );
224 : 0 : if ( sqrDist < 0 )
225 : 0 : return;
226 : :
227 : 0 : QgsPointXY edgePoints[2];
228 : 0 : edgePoints[0] = geom->vertexAt( afterVertex - 1 );
229 : 0 : edgePoints[1] = geom->vertexAt( afterVertex );
230 : 0 : pt = QgsPointXY( ( edgePoints[0].x() + edgePoints[1].x() ) / 2.0, ( edgePoints[0].y() + edgePoints[1].y() ) / 2.0 );
231 : :
232 : 0 : QgsPointLocator::Match m( QgsPointLocator::MiddleOfSegment, mLocator->mLayer, id, std::sqrt( mSrcPoint.sqrDist( pt ) ), pt, afterVertex - 1 );
233 : : // in range queries the filter may reject some matches
234 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
235 : 0 : return;
236 : :
237 : 0 : if ( !mBest.isValid() || m.distance() < mBest.distance() )
238 : 0 : mBest = m;
239 : :
240 : 0 : }
241 : :
242 : : private:
243 : : QgsPointLocator *mLocator = nullptr;
244 : : QgsPointLocator::Match &mBest;
245 : : QgsPointXY mSrcPoint;
246 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
247 : : };
248 : :
249 : : ////////////////////////////////////////////////////////////////////////////
250 : :
251 : : /**
252 : : * \ingroup core
253 : : * \brief Helper class used when traversing the index looking for line endpoints (start or end vertex) - builds a list of matches.
254 : : * \note not available in Python bindings
255 : : * \since QGIS 3.20
256 : : */
257 : 0 : class QgsPointLocator_VisitorNearestLineEndpoint : public IVisitor
258 : : {
259 : : public:
260 : :
261 : : /**
262 : : * \ingroup core
263 : : * \brief Helper class used when traversing the index looking for line endpoints (start or end vertex) - builds a list of matches.
264 : : */
265 : 0 : QgsPointLocator_VisitorNearestLineEndpoint( QgsPointLocator *pl, QgsPointLocator::Match &m, const QgsPointXY &srcPoint, QgsPointLocator::MatchFilter *filter = nullptr )
266 : 0 : : mLocator( pl )
267 : 0 : , mBest( m )
268 : 0 : , mSrcPoint( srcPoint )
269 : 0 : , mFilter( filter )
270 : 0 : {}
271 : :
272 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
273 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
274 : :
275 : 0 : void visitData( const IData &d ) override
276 : : {
277 : 0 : QgsFeatureId id = d.getIdentifier();
278 : 0 : const QgsGeometry *geom = mLocator->mGeoms.value( id );
279 : :
280 : 0 : QgsPointXY bestPoint;
281 : 0 : int bestVertexNumber = -1;
282 : 0 : auto replaceIfBetter = [this, &bestPoint, &bestVertexNumber]( const QgsPoint & candidate, int vertexNumber )
283 : : {
284 : 0 : if ( bestPoint.isEmpty() || candidate.distanceSquared( mSrcPoint.x(), mSrcPoint.y() ) < bestPoint.sqrDist( mSrcPoint ) )
285 : : {
286 : 0 : bestPoint = QgsPointXY( candidate );
287 : 0 : bestVertexNumber = vertexNumber;
288 : 0 : }
289 : 0 : };
290 : :
291 : 0 : switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
292 : : {
293 : : case QgsWkbTypes::PointGeometry:
294 : : case QgsWkbTypes::UnknownGeometry:
295 : : case QgsWkbTypes::NullGeometry:
296 : 0 : return;
297 : :
298 : : case QgsWkbTypes::LineGeometry:
299 : : {
300 : 0 : int partStartVertexNum = 0;
301 : 0 : for ( auto partIt = geom->const_parts_begin(); partIt != geom->const_parts_end(); ++partIt )
302 : : {
303 : 0 : if ( const QgsCurve *curve = qgsgeometry_cast< const QgsCurve * >( *partIt ) )
304 : : {
305 : 0 : replaceIfBetter( curve->startPoint(), partStartVertexNum );
306 : 0 : replaceIfBetter( curve->endPoint(), partStartVertexNum + curve->numPoints() - 1 );
307 : 0 : partStartVertexNum += curve->numPoints();
308 : 0 : }
309 : 0 : }
310 : 0 : break;
311 : : }
312 : :
313 : : case QgsWkbTypes::PolygonGeometry:
314 : : {
315 : 0 : int partStartVertexNum = 0;
316 : 0 : for ( auto partIt = geom->const_parts_begin(); partIt != geom->const_parts_end(); ++partIt )
317 : : {
318 : 0 : if ( const QgsCurvePolygon *polygon = qgsgeometry_cast< const QgsCurvePolygon * >( *partIt ) )
319 : : {
320 : 0 : if ( polygon->exteriorRing() )
321 : : {
322 : 0 : replaceIfBetter( polygon->exteriorRing()->startPoint(), partStartVertexNum );
323 : 0 : partStartVertexNum += polygon->exteriorRing()->numPoints();
324 : 0 : }
325 : 0 : for ( int i = 0; i < polygon->numInteriorRings(); ++i )
326 : : {
327 : 0 : const QgsCurve *ring = polygon->interiorRing( i );
328 : 0 : replaceIfBetter( ring->startPoint(), partStartVertexNum );
329 : 0 : partStartVertexNum += ring->numPoints();
330 : 0 : }
331 : 0 : }
332 : 0 : }
333 : 0 : break;
334 : : }
335 : : }
336 : :
337 : 0 : QgsPointLocator::Match m( QgsPointLocator::LineEndpoint, mLocator->mLayer, id, std::sqrt( mSrcPoint.sqrDist( bestPoint ) ), bestPoint, bestVertexNumber );
338 : : // in range queries the filter may reject some matches
339 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
340 : 0 : return;
341 : :
342 : 0 : if ( !mBest.isValid() || m.distance() < mBest.distance() )
343 : 0 : mBest = m;
344 : 0 : }
345 : :
346 : : private:
347 : : QgsPointLocator *mLocator = nullptr;
348 : : QgsPointLocator::Match &mBest;
349 : : QgsPointXY mSrcPoint;
350 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
351 : : };
352 : :
353 : :
354 : : ////////////////////////////////////////////////////////////////////////////
355 : :
356 : :
357 : : /**
358 : : * \ingroup core
359 : : * \brief Helper class used when traversing the index looking for edges - builds a list of matches.
360 : : * \note not available in Python bindings
361 : : */
362 : 0 : class QgsPointLocator_VisitorNearestEdge : public IVisitor
363 : : {
364 : : public:
365 : 0 : QgsPointLocator_VisitorNearestEdge( QgsPointLocator *pl, QgsPointLocator::Match &m, const QgsPointXY &srcPoint, QgsPointLocator::MatchFilter *filter = nullptr )
366 : 0 : : mLocator( pl )
367 : 0 : , mBest( m )
368 : 0 : , mSrcPoint( srcPoint )
369 : 0 : , mFilter( filter )
370 : 0 : {}
371 : :
372 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
373 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
374 : :
375 : 0 : void visitData( const IData &d ) override
376 : : {
377 : 0 : QgsFeatureId id = d.getIdentifier();
378 : 0 : QgsGeometry *geom = mLocator->mGeoms.value( id );
379 : 0 : QgsPointXY pt;
380 : : int afterVertex;
381 : 0 : double sqrDist = geom->closestSegmentWithContext( mSrcPoint, pt, afterVertex, nullptr, POINT_LOC_EPSILON );
382 : 0 : if ( sqrDist < 0 )
383 : 0 : return;
384 : :
385 : 0 : QgsPointXY edgePoints[2];
386 : 0 : edgePoints[0] = geom->vertexAt( afterVertex - 1 );
387 : 0 : edgePoints[1] = geom->vertexAt( afterVertex );
388 : 0 : QgsPointLocator::Match m( QgsPointLocator::Edge, mLocator->mLayer, id, std::sqrt( sqrDist ), pt, afterVertex - 1, edgePoints );
389 : : // in range queries the filter may reject some matches
390 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
391 : 0 : return;
392 : :
393 : 0 : if ( !mBest.isValid() || m.distance() < mBest.distance() )
394 : 0 : mBest = m;
395 : 0 : }
396 : :
397 : : private:
398 : : QgsPointLocator *mLocator = nullptr;
399 : : QgsPointLocator::Match &mBest;
400 : : QgsPointXY mSrcPoint;
401 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
402 : : };
403 : :
404 : :
405 : : ////////////////////////////////////////////////////////////////////////////
406 : :
407 : : /**
408 : : * \ingroup core
409 : : * \brief Helper class used when traversing the index with areas - builds a list of matches.
410 : : * \note not available in Python bindings
411 : : */
412 : 0 : class QgsPointLocator_VisitorArea : public IVisitor
413 : : {
414 : : public:
415 : : //! constructor
416 : 0 : QgsPointLocator_VisitorArea( QgsPointLocator *pl, const QgsPointXY &origPt, QgsPointLocator::MatchList &list )
417 : 0 : : mLocator( pl )
418 : 0 : , mList( list )
419 : 0 : , mGeomPt( QgsGeometry::fromPointXY( origPt ) )
420 : 0 : {}
421 : :
422 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
423 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
424 : :
425 : 0 : void visitData( const IData &d ) override
426 : : {
427 : 0 : QgsFeatureId id = d.getIdentifier();
428 : 0 : QgsGeometry *g = mLocator->mGeoms.value( id );
429 : 0 : if ( g->intersects( mGeomPt ) )
430 : 0 : mList << QgsPointLocator::Match( QgsPointLocator::Area, mLocator->mLayer, id, 0, mGeomPt.asPoint() );
431 : 0 : }
432 : : private:
433 : : QgsPointLocator *mLocator = nullptr;
434 : : QgsPointLocator::MatchList &mList;
435 : : QgsGeometry mGeomPt;
436 : : };
437 : :
438 : :
439 : : ////////////////////////////////////////////////////////////////////////////
440 : :
441 : : // code adapted from
442 : : // http://en.wikipedia.org/wiki/Cohen%E2%80%93Sutherland_algorithm
443 : : struct _CohenSutherland
444 : : {
445 : 0 : explicit _CohenSutherland( const QgsRectangle &rect ) : mRect( rect ) {}
446 : :
447 : : typedef int OutCode;
448 : :
449 : : static const int INSIDE = 0; // 0000
450 : : static const int LEFT = 1; // 0001
451 : : static const int RIGHT = 2; // 0010
452 : : static const int BOTTOM = 4; // 0100
453 : : static const int TOP = 8; // 1000
454 : :
455 : : QgsRectangle mRect;
456 : :
457 : 0 : OutCode computeOutCode( double x, double y )
458 : : {
459 : 0 : OutCode code = INSIDE; // initialized as being inside of clip window
460 : 0 : if ( x < mRect.xMinimum() ) // to the left of clip window
461 : 0 : code |= LEFT;
462 : 0 : else if ( x > mRect.xMaximum() ) // to the right of clip window
463 : 0 : code |= RIGHT;
464 : 0 : if ( y < mRect.yMinimum() ) // below the clip window
465 : 0 : code |= BOTTOM;
466 : 0 : else if ( y > mRect.yMaximum() ) // above the clip window
467 : 0 : code |= TOP;
468 : 0 : return code;
469 : : }
470 : :
471 : 0 : bool isSegmentInRect( double x0, double y0, double x1, double y1 )
472 : : {
473 : : // compute outcodes for P0, P1, and whatever point lies outside the clip rectangle
474 : 0 : OutCode outcode0 = computeOutCode( x0, y0 );
475 : 0 : OutCode outcode1 = computeOutCode( x1, y1 );
476 : 0 : bool accept = false;
477 : 0 :
478 : 0 : while ( true )
479 : 0 : {
480 : 0 : if ( !( outcode0 | outcode1 ) )
481 : : {
482 : : // Bitwise OR is 0. Trivially accept and get out of loop
483 : 0 : accept = true;
484 : 0 : break;
485 : : }
486 : 0 : else if ( outcode0 & outcode1 )
487 : : {
488 : : // Bitwise AND is not 0. Trivially reject and get out of loop
489 : 0 : break;
490 : : }
491 : : else
492 : : {
493 : : // failed both tests, so calculate the line segment to clip
494 : : // from an outside point to an intersection with clip edge
495 : : double x, y;
496 : :
497 : : // At least one endpoint is outside the clip rectangle; pick it.
498 : 0 : OutCode outcodeOut = outcode0 ? outcode0 : outcode1;
499 : :
500 : : // Now find the intersection point;
501 : : // use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0)
502 : 0 : if ( outcodeOut & TOP )
503 : : {
504 : : // point is above the clip rectangle
505 : 0 : x = x0 + ( x1 - x0 ) * ( mRect.yMaximum() - y0 ) / ( y1 - y0 );
506 : 0 : y = mRect.yMaximum();
507 : 0 : }
508 : 0 : else if ( outcodeOut & BOTTOM )
509 : : {
510 : : // point is below the clip rectangle
511 : 0 : x = x0 + ( x1 - x0 ) * ( mRect.yMinimum() - y0 ) / ( y1 - y0 );
512 : 0 : y = mRect.yMinimum();
513 : 0 : }
514 : 0 : else if ( outcodeOut & RIGHT )
515 : : {
516 : : // point is to the right of clip rectangle
517 : 0 : y = y0 + ( y1 - y0 ) * ( mRect.xMaximum() - x0 ) / ( x1 - x0 );
518 : 0 : x = mRect.xMaximum();
519 : 0 : }
520 : 0 : else if ( outcodeOut & LEFT )
521 : : {
522 : : // point is to the left of clip rectangle
523 : 0 : y = y0 + ( y1 - y0 ) * ( mRect.xMinimum() - x0 ) / ( x1 - x0 );
524 : 0 : x = mRect.xMinimum();
525 : 0 : }
526 : : else
527 : 0 : break;
528 : :
529 : : // Now we move outside point to intersection point to clip
530 : : // and get ready for next pass.
531 : 0 : if ( outcodeOut == outcode0 )
532 : : {
533 : 0 : x0 = x;
534 : 0 : y0 = y;
535 : 0 : outcode0 = computeOutCode( x0, y0 );
536 : 0 : }
537 : : else
538 : : {
539 : 0 : x1 = x;
540 : 0 : y1 = y;
541 : 0 : outcode1 = computeOutCode( x1, y1 );
542 : : }
543 : : }
544 : : }
545 : 0 : return accept;
546 : : }
547 : : };
548 : :
549 : :
550 : 0 : static QgsPointLocator::MatchList _geometrySegmentsInRect( QgsGeometry *geom, const QgsRectangle &rect, QgsVectorLayer *vl, QgsFeatureId fid )
551 : : {
552 : : // this code is stupidly based on QgsGeometry::closestSegmentWithContext
553 : : // we need iterator for segments...
554 : :
555 : 0 : QgsPointLocator::MatchList lst;
556 : :
557 : : // geom is converted to a MultiCurve
558 : 0 : QgsGeometry straightGeom = geom->convertToType( QgsWkbTypes::LineGeometry, true );
559 : : // and convert to straight segemnt / converts curve to linestring
560 : 0 : straightGeom.convertToStraightSegment();
561 : :
562 : : // so, you must have multilinestring
563 : : //
564 : : // Special case: Intersections cannot be done on an empty linestring like
565 : : // QgsGeometry(QgsLineString()) or QgsGeometry::fromWkt("LINESTRING EMPTY")
566 : 0 : if ( straightGeom.isEmpty() || ( ( straightGeom.type() != QgsWkbTypes::LineGeometry ) && ( !straightGeom.isMultipart() ) ) )
567 : 0 : return lst;
568 : :
569 : 0 : _CohenSutherland cs( rect );
570 : :
571 : 0 : int pointIndex = 0;
572 : 0 : for ( auto part = straightGeom.const_parts_begin(); part != straightGeom.const_parts_end(); ++part )
573 : : {
574 : : // Checking for invalid linestrings
575 : : // A linestring should/(must?) have at least two points.
576 : 0 : QgsCurve *curve = qgsgeometry_cast<QgsCurve *>( *part );
577 : : Q_ASSERT( !curve->hasCurvedSegments() );
578 : 0 : if ( curve->numPoints() < 2 )
579 : 0 : continue;
580 : :
581 : 0 : QgsAbstractGeometry::vertex_iterator it = ( *part )->vertices_begin();
582 : 0 : QgsPointXY prevPoint( *it );
583 : 0 : it++;
584 : 0 : while ( it != ( *part )->vertices_end() )
585 : : {
586 : 0 : QgsPointXY thisPoint( *it );
587 : 0 : if ( cs.isSegmentInRect( prevPoint.x(), prevPoint.y(), thisPoint.x(), thisPoint.y() ) )
588 : : {
589 : 0 : QgsPointXY edgePoints[2];
590 : 0 : edgePoints[0] = prevPoint;
591 : 0 : edgePoints[1] = thisPoint;
592 : 0 : lst << QgsPointLocator::Match( QgsPointLocator::Edge, vl, fid, 0, QgsPointXY(), pointIndex - 1, edgePoints );
593 : 0 : }
594 : 0 : prevPoint = QgsPointXY( *it );
595 : 0 : it++;
596 : 0 : pointIndex += 1;
597 : :
598 : : }
599 : 0 : }
600 : 0 : return lst;
601 : 0 : }
602 : :
603 : : /**
604 : : * \ingroup core
605 : : * \brief Helper class used when traversing the index looking for edges - builds a list of matches.
606 : : * \note not available in Python bindings
607 : : */
608 : 0 : class QgsPointLocator_VisitorEdgesInRect : public IVisitor
609 : : {
610 : : public:
611 : 0 : QgsPointLocator_VisitorEdgesInRect( QgsPointLocator *pl, QgsPointLocator::MatchList &lst, const QgsRectangle &srcRect, QgsPointLocator::MatchFilter *filter = nullptr )
612 : 0 : : mLocator( pl )
613 : 0 : , mList( lst )
614 : 0 : , mSrcRect( srcRect )
615 : 0 : , mFilter( filter )
616 : 0 : {}
617 : :
618 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
619 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
620 : :
621 : 0 : void visitData( const IData &d ) override
622 : : {
623 : 0 : QgsFeatureId id = d.getIdentifier();
624 : 0 : QgsGeometry *geom = mLocator->mGeoms.value( id );
625 : :
626 : 0 : const auto segmentsInRect {_geometrySegmentsInRect( geom, mSrcRect, mLocator->mLayer, id )};
627 : 0 : for ( const QgsPointLocator::Match &m : segmentsInRect )
628 : : {
629 : : // in range queries the filter may reject some matches
630 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
631 : 0 : continue;
632 : :
633 : 0 : mList << m;
634 : : }
635 : 0 : }
636 : :
637 : : private:
638 : : QgsPointLocator *mLocator = nullptr;
639 : : QgsPointLocator::MatchList &mList;
640 : : QgsRectangle mSrcRect;
641 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
642 : : };
643 : :
644 : : ////////////////////////////////////////////////////////////////////////////
645 : :
646 : : /**
647 : : * \ingroup core
648 : : * \brief Helper class used when traversing the index looking for vertices - builds a list of matches.
649 : : * \note not available in Python bindings
650 : : * \since QGIS 3.6
651 : : */
652 : 0 : class QgsPointLocator_VisitorVerticesInRect : public IVisitor
653 : : {
654 : : public:
655 : : //! Constructs the visitor
656 : 0 : QgsPointLocator_VisitorVerticesInRect( QgsPointLocator *pl, QgsPointLocator::MatchList &lst, const QgsRectangle &srcRect, QgsPointLocator::MatchFilter *filter = nullptr )
657 : 0 : : mLocator( pl )
658 : 0 : , mList( lst )
659 : 0 : , mSrcRect( srcRect )
660 : 0 : , mFilter( filter )
661 : 0 : {}
662 : :
663 : 0 : void visitNode( const INode &n ) override { Q_UNUSED( n ) }
664 : 0 : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ) }
665 : :
666 : 0 : void visitData( const IData &d ) override
667 : : {
668 : 0 : QgsFeatureId id = d.getIdentifier();
669 : 0 : const QgsGeometry *geom = mLocator->mGeoms.value( id );
670 : :
671 : 0 : for ( QgsAbstractGeometry::vertex_iterator it = geom->vertices_begin(); it != geom->vertices_end(); ++it )
672 : : {
673 : 0 : if ( mSrcRect.contains( *it ) )
674 : : {
675 : 0 : QgsPointLocator::Match m( QgsPointLocator::Vertex, mLocator->mLayer, id, 0, *it, geom->vertexNrFromVertexId( it.vertexId() ) );
676 : :
677 : : // in range queries the filter may reject some matches
678 : 0 : if ( mFilter && !mFilter->acceptMatch( m ) )
679 : 0 : continue;
680 : :
681 : 0 : mList << m;
682 : 0 : }
683 : 0 : }
684 : 0 : }
685 : :
686 : : private:
687 : : QgsPointLocator *mLocator = nullptr;
688 : : QgsPointLocator::MatchList &mList;
689 : : QgsRectangle mSrcRect;
690 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
691 : : };
692 : :
693 : : /**
694 : : * \ingroup core
695 : : * \brief Helper class used when traversing the index looking for centroid - builds a list of matches.
696 : : * \note not available in Python bindings
697 : : * \since QGIS 3.10
698 : : */
699 : : class QgsPointLocator_VisitorCentroidsInRect : public IVisitor
700 : : {
701 : : public:
702 : : //! Constructs the visitor
703 : : QgsPointLocator_VisitorCentroidsInRect( QgsPointLocator *pl, QgsPointLocator::MatchList &lst, const QgsRectangle &srcRect, QgsPointLocator::MatchFilter *filter = nullptr )
704 : : : mLocator( pl )
705 : : , mList( lst )
706 : : , mSrcRect( srcRect )
707 : : , mFilter( filter )
708 : : {}
709 : :
710 : : void visitNode( const INode &n ) override { Q_UNUSED( n ); }
711 : : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ); }
712 : :
713 : : void visitData( const IData &d ) override
714 : : {
715 : : QgsFeatureId id = d.getIdentifier();
716 : : const QgsGeometry *geom = mLocator->mGeoms.value( id );
717 : : const QgsPointXY centroid = geom->centroid().asPoint();
718 : : if ( mSrcRect.contains( centroid ) )
719 : : {
720 : : QgsPointLocator::Match m( QgsPointLocator::Centroid, mLocator->mLayer, id, 0, centroid, -1 );
721 : :
722 : : // in range queries the filter may reject some matches
723 : : if ( !( mFilter && !mFilter->acceptMatch( m ) ) )
724 : : mList << m;
725 : : }
726 : : }
727 : :
728 : : private:
729 : : QgsPointLocator *mLocator = nullptr;
730 : : QgsPointLocator::MatchList &mList;
731 : : QgsRectangle mSrcRect;
732 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
733 : : };
734 : :
735 : : /**
736 : : * \ingroup core
737 : : * \brief Helper class used when traversing the index looking for middle segment - builds a list of matches.
738 : : * \note not available in Python bindings
739 : : * \since QGIS 3.10
740 : : */
741 : : class QgsPointLocator_VisitorMiddlesInRect : public IVisitor
742 : : {
743 : : public:
744 : : //! Constructs the visitor
745 : : QgsPointLocator_VisitorMiddlesInRect( QgsPointLocator *pl, QgsPointLocator::MatchList &lst, const QgsRectangle &srcRect, QgsPointLocator::MatchFilter *filter = nullptr )
746 : : : mLocator( pl )
747 : : , mList( lst )
748 : : , mSrcRect( srcRect )
749 : : , mFilter( filter )
750 : : {}
751 : :
752 : : void visitNode( const INode &n ) override { Q_UNUSED( n ); }
753 : : void visitData( std::vector<const IData *> &v ) override { Q_UNUSED( v ); }
754 : :
755 : : void visitData( const IData &d ) override
756 : : {
757 : : QgsFeatureId id = d.getIdentifier();
758 : : const QgsGeometry *geom = mLocator->mGeoms.value( id );
759 : :
760 : : for ( QgsAbstractGeometry::const_part_iterator itPart = geom->const_parts_begin() ; itPart != geom->const_parts_end() ; ++itPart )
761 : : {
762 : : QgsAbstractGeometry::vertex_iterator it = ( *itPart )->vertices_begin();
763 : : QgsAbstractGeometry::vertex_iterator itPrevious = ( *itPart )->vertices_begin();
764 : : it++;
765 : : for ( ; it != geom->vertices_end(); ++it, ++itPrevious )
766 : : {
767 : : QgsPointXY pt( ( ( *itPrevious ).x() + ( *it ).x() ) / 2.0, ( ( *itPrevious ).y() + ( *it ).y() ) / 2.0 );
768 : : if ( mSrcRect.contains( pt ) )
769 : : {
770 : : QgsPointLocator::Match m( QgsPointLocator::MiddleOfSegment, mLocator->mLayer, id, 0, pt, geom->vertexNrFromVertexId( it.vertexId() ) );
771 : :
772 : : // in range queries the filter may reject some matches
773 : : if ( mFilter && !mFilter->acceptMatch( m ) )
774 : : continue;
775 : :
776 : : mList << m;
777 : : }
778 : : }
779 : : }
780 : : }
781 : :
782 : : private:
783 : : QgsPointLocator *mLocator = nullptr;
784 : : QgsPointLocator::MatchList &mList;
785 : : QgsRectangle mSrcRect;
786 : : QgsPointLocator::MatchFilter *mFilter = nullptr;
787 : : };
788 : :
789 : : ////////////////////////////////////////////////////////////////////////////
790 : : #include <QStack>
791 : :
792 : : /**
793 : : * \ingroup core
794 : : * \brief Helper class to dump the R-index nodes and their content
795 : : * \note not available in Python bindings
796 : : */
797 : : class QgsPointLocator_DumpTree : public SpatialIndex::IQueryStrategy
798 : : {
799 : : private:
800 : : QStack<id_type> ids;
801 : :
802 : : public:
803 : :
804 : : void getNextEntry( const IEntry &entry, id_type &nextEntry, bool &hasNext ) override
805 : : {
806 : : const INode *n = dynamic_cast<const INode *>( &entry );
807 : : if ( !n )
808 : : return;
809 : :
810 : : QgsDebugMsgLevel( QStringLiteral( "NODE: %1" ).arg( n->getIdentifier() ), 4 );
811 : : if ( n->getLevel() > 0 )
812 : : {
813 : : // inner nodes
814 : : for ( uint32_t cChild = 0; cChild < n->getChildrenCount(); cChild++ )
815 : : {
816 : : QgsDebugMsgLevel( QStringLiteral( "- CH: %1" ).arg( n->getChildIdentifier( cChild ) ), 4 );
817 : : ids.push( n->getChildIdentifier( cChild ) );
818 : : }
819 : : }
820 : : else
821 : : {
822 : : // leaves
823 : : for ( uint32_t cChild = 0; cChild < n->getChildrenCount(); cChild++ )
824 : : {
825 : : QgsDebugMsgLevel( QStringLiteral( "- L: %1" ).arg( n->getChildIdentifier( cChild ) ), 4 );
826 : : }
827 : : }
828 : :
829 : : if ( ! ids.empty() )
830 : : {
831 : : nextEntry = ids.back();
832 : : ids.pop();
833 : : hasNext = true;
834 : : }
835 : : else
836 : : hasNext = false;
837 : : }
838 : : };
839 : :
840 : : ////////////////////////////////////////////////////////////////////////////
841 : :
842 : :
843 : 0 : QgsPointLocator::QgsPointLocator( QgsVectorLayer *layer, const QgsCoordinateReferenceSystem &destCRS, const QgsCoordinateTransformContext &transformContext, const QgsRectangle *extent )
844 : 0 : : mLayer( layer )
845 : 0 : {
846 : 0 : if ( destCRS.isValid() )
847 : : {
848 : 0 : mTransform = QgsCoordinateTransform( layer->crs(), destCRS, transformContext );
849 : 0 : }
850 : :
851 : 0 : setExtent( extent );
852 : :
853 : 0 : mStorage.reset( StorageManager::createNewMemoryStorageManager() );
854 : :
855 : 0 : connect( mLayer, &QgsVectorLayer::featureAdded, this, &QgsPointLocator::onFeatureAdded );
856 : 0 : connect( mLayer, &QgsVectorLayer::featureDeleted, this, &QgsPointLocator::onFeatureDeleted );
857 : 0 : connect( mLayer, &QgsVectorLayer::geometryChanged, this, &QgsPointLocator::onGeometryChanged );
858 : 0 : connect( mLayer, &QgsVectorLayer::attributeValueChanged, this, &QgsPointLocator::onAttributeValueChanged );
859 : 0 : connect( mLayer, &QgsVectorLayer::dataChanged, this, &QgsPointLocator::destroyIndex );
860 : 0 : }
861 : :
862 : :
863 : 0 : QgsPointLocator::~QgsPointLocator()
864 : 0 : {
865 : : // don't delete a locator if there is an indexing task running on it
866 : 0 : mIsDestroying = true;
867 : 0 : if ( mIsIndexing )
868 : 0 : waitForIndexingFinished();
869 : :
870 : 0 : destroyIndex();
871 : 0 : }
872 : :
873 : 0 : QgsCoordinateReferenceSystem QgsPointLocator::destinationCrs() const
874 : : {
875 : 0 : return mTransform.isValid() ? mTransform.destinationCrs() : QgsCoordinateReferenceSystem();
876 : : }
877 : :
878 : 0 : void QgsPointLocator::setExtent( const QgsRectangle *extent )
879 : : {
880 : 0 : if ( mIsIndexing )
881 : : // already indexing, return!
882 : 0 : return;
883 : :
884 : 0 : mExtent.reset( extent ? new QgsRectangle( *extent ) : nullptr );
885 : :
886 : 0 : destroyIndex();
887 : 0 : }
888 : :
889 : 0 : void QgsPointLocator::setRenderContext( const QgsRenderContext *context )
890 : : {
891 : 0 : if ( mIsIndexing )
892 : : // already indexing, return!
893 : 0 : return;
894 : :
895 : 0 : disconnect( mLayer, &QgsVectorLayer::styleChanged, this, &QgsPointLocator::destroyIndex );
896 : :
897 : 0 : destroyIndex();
898 : 0 : mContext.reset( nullptr );
899 : :
900 : 0 : if ( context )
901 : : {
902 : 0 : mContext = std::unique_ptr<QgsRenderContext>( new QgsRenderContext( *context ) );
903 : 0 : connect( mLayer, &QgsVectorLayer::styleChanged, this, &QgsPointLocator::destroyIndex );
904 : 0 : }
905 : :
906 : 0 : }
907 : :
908 : 0 : void QgsPointLocator::onInitTaskFinished()
909 : : {
910 : : // Check that we don't call this method twice, when calling waitForFinished
911 : : // for instance (because of taskCompleted signal)
912 : 0 : if ( !mIsIndexing )
913 : 0 : return;
914 : :
915 : 0 : if ( mIsDestroying )
916 : 0 : return;
917 : :
918 : 0 : mIsIndexing = false;
919 : 0 : mRenderer.reset();
920 : 0 : mSource.reset();
921 : :
922 : : // treat added and deleted feature while indexing
923 : 0 : for ( QgsFeatureId fid : mAddedFeatures )
924 : 0 : onFeatureAdded( fid );
925 : 0 : mAddedFeatures.clear();
926 : :
927 : 0 : for ( QgsFeatureId fid : mDeletedFeatures )
928 : 0 : onFeatureDeleted( fid );
929 : 0 : mDeletedFeatures.clear();
930 : :
931 : 0 : emit initFinished( mInitTask->isBuildOK() );
932 : 0 : }
933 : :
934 : 0 : bool QgsPointLocator::init( int maxFeaturesToIndex, bool relaxed )
935 : : {
936 : 0 : const QgsWkbTypes::GeometryType geomType = mLayer->geometryType();
937 : 0 : if ( geomType == QgsWkbTypes::NullGeometry // nothing to index
938 : 0 : || hasIndex()
939 : 0 : || mIsIndexing ) // already indexing, return!
940 : 0 : return true;
941 : :
942 : 0 : if ( !mLayer->dataProvider()
943 : 0 : || !mLayer->dataProvider()->isValid() )
944 : 0 : return false;
945 : :
946 : 0 : mRenderer.reset( mLayer->renderer() ? mLayer->renderer()->clone() : nullptr );
947 : 0 : mSource.reset( new QgsVectorLayerFeatureSource( mLayer ) );
948 : :
949 : 0 : if ( mContext )
950 : : {
951 : 0 : mContext->expressionContext() << QgsExpressionContextUtils::layerScope( mLayer );
952 : 0 : }
953 : :
954 : 0 : mIsIndexing = true;
955 : :
956 : 0 : if ( relaxed )
957 : : {
958 : 0 : mInitTask = new QgsPointLocatorInitTask( this );
959 : 0 : connect( mInitTask, &QgsPointLocatorInitTask::taskTerminated, this, &QgsPointLocator::onInitTaskFinished );
960 : 0 : connect( mInitTask, &QgsPointLocatorInitTask::taskCompleted, this, &QgsPointLocator::onInitTaskFinished );
961 : 0 : QgsApplication::taskManager()->addTask( mInitTask );
962 : 0 : return true;
963 : : }
964 : : else
965 : : {
966 : 0 : const bool ok = rebuildIndex( maxFeaturesToIndex );
967 : 0 : mIsIndexing = false;
968 : 0 : emit initFinished( ok );
969 : 0 : return ok;
970 : : }
971 : 0 : }
972 : :
973 : 0 : void QgsPointLocator::waitForIndexingFinished()
974 : : {
975 : 0 : mInitTask->waitForFinished();
976 : :
977 : 0 : if ( !mIsDestroying )
978 : 0 : onInitTaskFinished();
979 : 0 : }
980 : :
981 : 0 : bool QgsPointLocator::hasIndex() const
982 : : {
983 : 0 : return mIsIndexing || mRTree || mIsEmptyLayer;
984 : : }
985 : :
986 : 0 : bool QgsPointLocator::prepare( bool relaxed )
987 : : {
988 : 0 : if ( mIsIndexing )
989 : : {
990 : 0 : if ( relaxed )
991 : 0 : return false;
992 : : else
993 : 0 : waitForIndexingFinished();
994 : 0 : }
995 : :
996 : 0 : if ( !mRTree )
997 : : {
998 : 0 : init( -1, relaxed );
999 : 0 : if ( ( relaxed && mIsIndexing ) || !mRTree ) // relaxed mode and currently indexing or still invalid?
1000 : 0 : return false;
1001 : 0 : }
1002 : :
1003 : 0 : return true;
1004 : 0 : }
1005 : :
1006 : 0 : bool QgsPointLocator::rebuildIndex( int maxFeaturesToIndex )
1007 : : {
1008 : 0 : QElapsedTimer t;
1009 : 0 : t.start();
1010 : :
1011 : 0 : QgsDebugMsgLevel( QStringLiteral( "RebuildIndex start : %1" ).arg( mSource->id() ), 2 );
1012 : :
1013 : 0 : destroyIndex();
1014 : :
1015 : 0 : QLinkedList<RTree::Data *> dataList;
1016 : 0 : QgsFeature f;
1017 : :
1018 : 0 : QgsFeatureRequest request;
1019 : 0 : request.setNoAttributes();
1020 : :
1021 : 0 : if ( mExtent )
1022 : : {
1023 : 0 : QgsRectangle rect = *mExtent;
1024 : 0 : if ( mTransform.isValid() )
1025 : : {
1026 : : try
1027 : : {
1028 : 0 : rect = mTransform.transformBoundingBox( rect, QgsCoordinateTransform::ReverseTransform );
1029 : 0 : }
1030 : : catch ( const QgsException &e )
1031 : : {
1032 : 0 : Q_UNUSED( e )
1033 : : // See https://github.com/qgis/QGIS/issues/20749
1034 : 0 : QgsDebugMsg( QStringLiteral( "could not transform bounding box to map, skipping the snap filter (%1)" ).arg( e.what() ) );
1035 : 0 : }
1036 : 0 : }
1037 : 0 : request.setFilterRect( rect );
1038 : 0 : }
1039 : :
1040 : 0 : bool filter = false;
1041 : 0 : QgsRenderContext *ctx = nullptr;
1042 : 0 : if ( mContext )
1043 : : {
1044 : 0 : ctx = mContext.get();
1045 : 0 : if ( mRenderer )
1046 : : {
1047 : : // setup scale for scale dependent visibility (rule based)
1048 : 0 : mRenderer->startRender( *ctx, mSource->fields() );
1049 : 0 : filter = mRenderer->capabilities() & QgsFeatureRenderer::Filter;
1050 : 0 : request.setSubsetOfAttributes( mRenderer->usedAttributes( *ctx ), mSource->fields() );
1051 : 0 : }
1052 : 0 : }
1053 : :
1054 : 0 : QgsFeatureIterator fi = mSource->getFeatures( request );
1055 : 0 : int indexedCount = 0;
1056 : :
1057 : 0 : while ( fi.nextFeature( f ) )
1058 : : {
1059 : 0 : if ( !f.hasGeometry() )
1060 : 0 : continue;
1061 : :
1062 : 0 : if ( filter && ctx && mRenderer )
1063 : : {
1064 : 0 : ctx->expressionContext().setFeature( f );
1065 : 0 : if ( !mRenderer->willRenderFeature( f, *ctx ) )
1066 : : {
1067 : 0 : continue;
1068 : : }
1069 : 0 : }
1070 : :
1071 : 0 : if ( mTransform.isValid() )
1072 : : {
1073 : : try
1074 : : {
1075 : 0 : QgsGeometry transformedGeometry = f.geometry();
1076 : 0 : transformedGeometry.transform( mTransform );
1077 : 0 : f.setGeometry( transformedGeometry );
1078 : 0 : }
1079 : : catch ( const QgsException &e )
1080 : : {
1081 : 0 : Q_UNUSED( e )
1082 : : // See https://github.com/qgis/QGIS/issues/20749
1083 : 0 : QgsDebugMsg( QStringLiteral( "could not transform geometry to map, skipping the snap for it (%1)" ).arg( e.what() ) );
1084 : : continue;
1085 : 0 : }
1086 : 0 : }
1087 : :
1088 : 0 : const QgsRectangle bbox = f.geometry().boundingBox();
1089 : 0 : if ( bbox.isFinite() )
1090 : : {
1091 : 0 : SpatialIndex::Region r( rect2region( bbox ) );
1092 : 0 : dataList << new RTree::Data( 0, nullptr, r, f.id() );
1093 : :
1094 : 0 : if ( mGeoms.contains( f.id() ) )
1095 : 0 : delete mGeoms.take( f.id() );
1096 : 0 : mGeoms[f.id()] = new QgsGeometry( f.geometry() );
1097 : 0 : ++indexedCount;
1098 : 0 : }
1099 : :
1100 : 0 : if ( maxFeaturesToIndex != -1 && indexedCount > maxFeaturesToIndex )
1101 : : {
1102 : 0 : qDeleteAll( dataList );
1103 : 0 : destroyIndex();
1104 : 0 : return false;
1105 : : }
1106 : : }
1107 : :
1108 : : // R-Tree parameters
1109 : 0 : double fillFactor = 0.7;
1110 : 0 : unsigned long indexCapacity = 10;
1111 : 0 : unsigned long leafCapacity = 10;
1112 : 0 : unsigned long dimension = 2;
1113 : 0 : RTree::RTreeVariant variant = RTree::RV_RSTAR;
1114 : : SpatialIndex::id_type indexId;
1115 : :
1116 : 0 : if ( dataList.isEmpty() )
1117 : : {
1118 : 0 : mIsEmptyLayer = true;
1119 : 0 : return true; // no features
1120 : : }
1121 : :
1122 : 0 : QgsPointLocator_Stream stream( dataList );
1123 : 0 : mRTree.reset( RTree::createAndBulkLoadNewRTree( RTree::BLM_STR, stream, *mStorage, fillFactor, indexCapacity,
1124 : 0 : leafCapacity, dimension, variant, indexId ) );
1125 : :
1126 : 0 : if ( ctx && mRenderer )
1127 : : {
1128 : 0 : mRenderer->stopRender( *ctx );
1129 : 0 : }
1130 : :
1131 : 0 : QgsDebugMsgLevel( QStringLiteral( "RebuildIndex end : %1 ms (%2)" ).arg( t.elapsed() ).arg( mSource->id() ), 2 );
1132 : :
1133 : 0 : return true;
1134 : 0 : }
1135 : :
1136 : :
1137 : 0 : void QgsPointLocator::destroyIndex()
1138 : : {
1139 : 0 : mRTree.reset();
1140 : :
1141 : 0 : mIsEmptyLayer = false;
1142 : :
1143 : 0 : qDeleteAll( mGeoms );
1144 : :
1145 : 0 : mGeoms.clear();
1146 : 0 : }
1147 : :
1148 : 0 : void QgsPointLocator::onFeatureAdded( QgsFeatureId fid )
1149 : : {
1150 : 0 : if ( mIsIndexing )
1151 : : {
1152 : : // will modify index once current indexing is finished
1153 : 0 : mAddedFeatures << fid;
1154 : 0 : return;
1155 : : }
1156 : :
1157 : 0 : if ( !mRTree )
1158 : : {
1159 : 0 : if ( mIsEmptyLayer )
1160 : : {
1161 : : // layer is not empty any more, let's build the index
1162 : 0 : mIsEmptyLayer = false;
1163 : 0 : init();
1164 : 0 : }
1165 : 0 : return; // nothing to do if we are not initialized yet
1166 : : }
1167 : :
1168 : 0 : QgsFeature f;
1169 : 0 : if ( mLayer->getFeatures( QgsFeatureRequest( fid ) ).nextFeature( f ) )
1170 : : {
1171 : 0 : if ( !f.hasGeometry() )
1172 : 0 : return;
1173 : :
1174 : 0 : if ( mContext )
1175 : : {
1176 : 0 : std::unique_ptr< QgsFeatureRenderer > renderer( mLayer->renderer() ? mLayer->renderer()->clone() : nullptr );
1177 : 0 : QgsRenderContext *ctx = nullptr;
1178 : :
1179 : 0 : mContext->expressionContext() << QgsExpressionContextUtils::layerScope( mLayer );
1180 : 0 : ctx = mContext.get();
1181 : 0 : if ( renderer && ctx )
1182 : : {
1183 : 0 : bool pass = false;
1184 : 0 : renderer->startRender( *ctx, mLayer->fields() );
1185 : :
1186 : 0 : ctx->expressionContext().setFeature( f );
1187 : 0 : if ( !renderer->willRenderFeature( f, *ctx ) )
1188 : : {
1189 : 0 : pass = true;
1190 : 0 : }
1191 : :
1192 : 0 : renderer->stopRender( *ctx );
1193 : 0 : if ( pass )
1194 : 0 : return;
1195 : 0 : }
1196 : 0 : }
1197 : :
1198 : 0 : if ( mTransform.isValid() )
1199 : : {
1200 : : try
1201 : : {
1202 : 0 : QgsGeometry transformedGeom = f.geometry();
1203 : 0 : transformedGeom.transform( mTransform );
1204 : 0 : f.setGeometry( transformedGeom );
1205 : 0 : }
1206 : : catch ( const QgsException &e )
1207 : : {
1208 : 0 : Q_UNUSED( e )
1209 : : // See https://github.com/qgis/QGIS/issues/20749
1210 : 0 : QgsDebugMsg( QStringLiteral( "could not transform geometry to map, skipping the snap for it (%1)" ).arg( e.what() ) );
1211 : : return;
1212 : 0 : }
1213 : 0 : }
1214 : :
1215 : 0 : const QgsRectangle bbox = f.geometry().boundingBox();
1216 : 0 : if ( bbox.isFinite() )
1217 : : {
1218 : 0 : SpatialIndex::Region r( rect2region( bbox ) );
1219 : 0 : mRTree->insertData( 0, nullptr, r, f.id() );
1220 : :
1221 : 0 : if ( mGeoms.contains( f.id() ) )
1222 : 0 : delete mGeoms.take( f.id() );
1223 : 0 : mGeoms[fid] = new QgsGeometry( f.geometry() );
1224 : 0 : }
1225 : 0 : }
1226 : 0 : }
1227 : :
1228 : 0 : void QgsPointLocator::onFeatureDeleted( QgsFeatureId fid )
1229 : : {
1230 : 0 : if ( mIsIndexing )
1231 : : {
1232 : 0 : if ( mAddedFeatures.contains( fid ) )
1233 : : {
1234 : 0 : mAddedFeatures.remove( fid );
1235 : 0 : }
1236 : : else
1237 : : {
1238 : : // will modify index once current indexing is finished
1239 : 0 : mDeletedFeatures << fid;
1240 : : }
1241 : 0 : return;
1242 : : }
1243 : :
1244 : 0 : if ( !mRTree )
1245 : 0 : return; // nothing to do if we are not initialized yet
1246 : :
1247 : 0 : if ( mGeoms.contains( fid ) )
1248 : : {
1249 : 0 : mRTree->deleteData( rect2region( mGeoms[fid]->boundingBox() ), fid );
1250 : 0 : delete mGeoms.take( fid );
1251 : 0 : }
1252 : :
1253 : 0 : }
1254 : :
1255 : 0 : void QgsPointLocator::onGeometryChanged( QgsFeatureId fid, const QgsGeometry &geom )
1256 : : {
1257 : 0 : Q_UNUSED( geom )
1258 : 0 : onFeatureDeleted( fid );
1259 : 0 : onFeatureAdded( fid );
1260 : 0 : }
1261 : :
1262 : 0 : void QgsPointLocator::onAttributeValueChanged( QgsFeatureId fid, int idx, const QVariant &value )
1263 : : {
1264 : : Q_UNUSED( idx )
1265 : 0 : Q_UNUSED( value )
1266 : 0 : if ( mContext )
1267 : : {
1268 : 0 : onFeatureDeleted( fid );
1269 : 0 : onFeatureAdded( fid );
1270 : 0 : }
1271 : 0 : }
1272 : :
1273 : :
1274 : 0 : QgsPointLocator::Match QgsPointLocator::nearestVertex( const QgsPointXY &point, double tolerance, MatchFilter *filter, bool relaxed )
1275 : : {
1276 : 0 : if ( !prepare( relaxed ) )
1277 : 0 : return Match();
1278 : :
1279 : 0 : Match m;
1280 : 0 : QgsPointLocator_VisitorNearestVertex visitor( this, m, point, filter );
1281 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1282 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1283 : 0 : if ( m.isValid() && m.distance() > tolerance )
1284 : 0 : return Match(); // make sure that only match strictly within the tolerance is returned
1285 : 0 : return m;
1286 : 0 : }
1287 : :
1288 : 0 : QgsPointLocator::Match QgsPointLocator::nearestCentroid( const QgsPointXY &point, double tolerance, MatchFilter *filter, bool relaxed )
1289 : : {
1290 : 0 : if ( !prepare( relaxed ) )
1291 : 0 : return Match();
1292 : :
1293 : 0 : Match m;
1294 : 0 : QgsPointLocator_VisitorNearestCentroid visitor( this, m, point, filter );
1295 : :
1296 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1297 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1298 : 0 : if ( m.isValid() && m.distance() > tolerance )
1299 : 0 : return Match(); // make sure that only match strictly within the tolerance is returned
1300 : 0 : return m;
1301 : 0 : }
1302 : :
1303 : 0 : QgsPointLocator::Match QgsPointLocator::nearestMiddleOfSegment( const QgsPointXY &point, double tolerance, MatchFilter *filter, bool relaxed )
1304 : : {
1305 : 0 : if ( !prepare( relaxed ) )
1306 : 0 : return Match();
1307 : :
1308 : 0 : Match m;
1309 : 0 : QgsPointLocator_VisitorNearestMiddleOfSegment visitor( this, m, point, filter );
1310 : :
1311 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1312 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1313 : 0 : if ( m.isValid() && m.distance() > tolerance )
1314 : 0 : return Match(); // make sure that only match strictly within the tolerance is returned
1315 : 0 : return m;
1316 : 0 : }
1317 : :
1318 : 0 : QgsPointLocator::Match QgsPointLocator::nearestLineEndpoints( const QgsPointXY &point, double tolerance, QgsPointLocator::MatchFilter *filter, bool relaxed )
1319 : : {
1320 : 0 : if ( !prepare( relaxed ) )
1321 : 0 : return Match();
1322 : :
1323 : 0 : Match m;
1324 : 0 : QgsPointLocator_VisitorNearestLineEndpoint visitor( this, m, point, filter );
1325 : :
1326 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1327 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1328 : 0 : if ( m.isValid() && m.distance() > tolerance )
1329 : 0 : return Match(); // make sure that only match strictly within the tolerance is returned
1330 : 0 : return m;
1331 : 0 : }
1332 : :
1333 : 0 : QgsPointLocator::Match QgsPointLocator::nearestEdge( const QgsPointXY &point, double tolerance, MatchFilter *filter, bool relaxed )
1334 : : {
1335 : 0 : if ( !prepare( relaxed ) )
1336 : 0 : return Match();
1337 : :
1338 : 0 : QgsWkbTypes::GeometryType geomType = mLayer->geometryType();
1339 : 0 : if ( geomType == QgsWkbTypes::PointGeometry )
1340 : 0 : return Match();
1341 : :
1342 : 0 : Match m;
1343 : 0 : QgsPointLocator_VisitorNearestEdge visitor( this, m, point, filter );
1344 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1345 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1346 : 0 : if ( m.isValid() && m.distance() > tolerance )
1347 : 0 : return Match(); // make sure that only match strictly within the tolerance is returned
1348 : 0 : return m;
1349 : 0 : }
1350 : :
1351 : 0 : QgsPointLocator::Match QgsPointLocator::nearestArea( const QgsPointXY &point, double tolerance, MatchFilter *filter, bool relaxed )
1352 : : {
1353 : 0 : if ( !prepare( relaxed ) )
1354 : 0 : return Match();
1355 : :
1356 : 0 : MatchList mlist = pointInPolygon( point );
1357 : 0 : if ( !mlist.isEmpty() && mlist.at( 0 ).isValid() )
1358 : : {
1359 : 0 : return mlist.at( 0 );
1360 : : }
1361 : :
1362 : 0 : if ( tolerance == 0 )
1363 : : {
1364 : 0 : return Match();
1365 : : }
1366 : :
1367 : : // discard point and line layers to keep only polygons
1368 : 0 : QgsWkbTypes::GeometryType geomType = mLayer->geometryType();
1369 : 0 : if ( geomType == QgsWkbTypes::PointGeometry || geomType == QgsWkbTypes::LineGeometry )
1370 : 0 : return Match();
1371 : :
1372 : : // use edges for adding tolerance
1373 : 0 : Match m = nearestEdge( point, tolerance, filter );
1374 : 0 : if ( m.isValid() )
1375 : 0 : return Match( Area, m.layer(), m.featureId(), m.distance(), m.point() );
1376 : : else
1377 : 0 : return Match();
1378 : 0 : }
1379 : :
1380 : :
1381 : 0 : QgsPointLocator::MatchList QgsPointLocator::edgesInRect( const QgsRectangle &rect, QgsPointLocator::MatchFilter *filter, bool relaxed )
1382 : : {
1383 : 0 : if ( !prepare( relaxed ) )
1384 : 0 : return MatchList();
1385 : :
1386 : 0 : QgsWkbTypes::GeometryType geomType = mLayer->geometryType();
1387 : 0 : if ( geomType == QgsWkbTypes::PointGeometry )
1388 : 0 : return MatchList();
1389 : :
1390 : 0 : MatchList lst;
1391 : 0 : QgsPointLocator_VisitorEdgesInRect visitor( this, lst, rect, filter );
1392 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1393 : :
1394 : 0 : return lst;
1395 : 0 : }
1396 : :
1397 : 0 : QgsPointLocator::MatchList QgsPointLocator::edgesInRect( const QgsPointXY &point, double tolerance, QgsPointLocator::MatchFilter *filter, bool relaxed )
1398 : : {
1399 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1400 : 0 : return edgesInRect( rect, filter, relaxed );
1401 : : }
1402 : :
1403 : 0 : QgsPointLocator::MatchList QgsPointLocator::verticesInRect( const QgsRectangle &rect, QgsPointLocator::MatchFilter *filter, bool relaxed )
1404 : : {
1405 : 0 : if ( !prepare( relaxed ) )
1406 : 0 : return MatchList();
1407 : :
1408 : 0 : MatchList lst;
1409 : 0 : QgsPointLocator_VisitorVerticesInRect visitor( this, lst, rect, filter );
1410 : 0 : mRTree->intersectsWithQuery( rect2region( rect ), visitor );
1411 : :
1412 : 0 : return lst;
1413 : 0 : }
1414 : :
1415 : 0 : QgsPointLocator::MatchList QgsPointLocator::verticesInRect( const QgsPointXY &point, double tolerance, QgsPointLocator::MatchFilter *filter, bool relaxed )
1416 : : {
1417 : 0 : QgsRectangle rect( point.x() - tolerance, point.y() - tolerance, point.x() + tolerance, point.y() + tolerance );
1418 : 0 : return verticesInRect( rect, filter, relaxed );
1419 : : }
1420 : :
1421 : 0 : QgsPointLocator::MatchList QgsPointLocator::pointInPolygon( const QgsPointXY &point, bool relaxed )
1422 : : {
1423 : 0 : if ( !prepare( relaxed ) )
1424 : 0 : return MatchList();
1425 : :
1426 : 0 : QgsWkbTypes::GeometryType geomType = mLayer->geometryType();
1427 : 0 : if ( geomType == QgsWkbTypes::PointGeometry || geomType == QgsWkbTypes::LineGeometry )
1428 : 0 : return MatchList();
1429 : :
1430 : 0 : MatchList lst;
1431 : 0 : QgsPointLocator_VisitorArea visitor( this, point, lst );
1432 : 0 : mRTree->intersectsWithQuery( point2point( point ), visitor );
1433 : 0 : return lst;
1434 : 0 : }
|
