QGIS API Documentation 3.36.0-Maidenhead (09951dc0acf)
Loading...
Searching...
No Matches
qgsmeshspatialindex.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsmeshspatialindex.cpp
3 -----------------------
4 begin : January 2019
5 copyright : (C) 2019 by Peter Petrik
6 email : zilolv 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 "qgsmeshspatialindex.h"
17#include "qgsrectangle.h"
18#include "qgslogger.h"
19#include "qgsfeedback.h"
20
21#include <spatialindex/SpatialIndex.h>
22#include <QMutex>
23#include <QMutexLocker>
24#include <memory>
25
26using namespace SpatialIndex;
27
29
30static Region faceToRegion( const QgsMesh &mesh, int id, bool &ok )
31{
32 const QgsMeshFace face = mesh.face( id );
33
34 if ( face.isEmpty() )
35 {
36 ok = false;
37 return Region();
38 }
39
40 const QVector<QgsMeshVertex> &vertices = mesh.vertices;
41
42 double xMinimum = vertices[face[0]].x();
43 double yMinimum = vertices[face[0]].y();
44 double xMaximum = vertices[face[0]].x();
45 double yMaximum = vertices[face[0]].y();
46
47 for ( int i = 1; i < face.size(); ++i )
48 {
49 xMinimum = std::min( vertices[face[i]].x(), xMinimum );
50 yMinimum = std::min( vertices[face[i]].y(), yMinimum );
51 xMaximum = std::max( vertices[face[i]].x(), xMaximum );
52 yMaximum = std::max( vertices[face[i]].y(), yMaximum );
53 }
54
55 double pt1[2] = { xMinimum, yMinimum };
56 double pt2[2] = { xMaximum, yMaximum };
57
58 ok = true;
59 return SpatialIndex::Region( pt1, pt2, 2 );
60}
61
62static Region edgeToRegion( const QgsMesh &mesh, int id, bool &ok )
63{
64 const QgsMeshEdge edge = mesh.edge( id );
65 const QgsMeshVertex firstVertex = mesh.vertices[edge.first];
66 const QgsMeshVertex secondVertex = mesh.vertices[edge.second];
67 const double xMinimum = std::min( firstVertex.x(), secondVertex.x() );
68 const double yMinimum = std::min( firstVertex.y(), secondVertex.y() );
69 const double xMaximum = std::max( firstVertex.x(), secondVertex.x() );
70 const double yMaximum = std::max( firstVertex.y(), secondVertex.y() );
71 double pt1[2] = { xMinimum, yMinimum };
72 double pt2[2] = { xMaximum, yMaximum };
73 ok = true;
74 return SpatialIndex::Region( pt1, pt2, 2 );
75}
76
77static Region rectToRegion( const QgsRectangle &rect )
78{
79 double pt1[2] = { rect.xMinimum(), rect.yMinimum() };
80 double pt2[2] = { rect.xMaximum(), rect.yMaximum() };
81 return SpatialIndex::Region( pt1, pt2, 2 );
82}
83
90class QgisMeshVisitor : public SpatialIndex::IVisitor
91{
92 public:
93 explicit QgisMeshVisitor( QList<int> &list )
94 : mList( list ) {}
95
96 void visitNode( const INode &n ) override
97 { Q_UNUSED( n ) }
98
99 void visitData( const IData &d ) override
100 {
101 mList.append( static_cast<int>( d.getIdentifier() ) );
102 }
103
104 void visitData( std::vector<const IData *> &v ) override
105 { Q_UNUSED( v ) }
106
107 private:
108 QList<int> &mList;
109};
110
117class QgsMeshSpatialIndexCopyVisitor : public SpatialIndex::IVisitor
118{
119 public:
120 explicit QgsMeshSpatialIndexCopyVisitor( SpatialIndex::ISpatialIndex *newIndex )
121 : mNewIndex( newIndex ) {}
122
123 void visitNode( const INode &n ) override
124 { Q_UNUSED( n ) }
125
126 void visitData( const IData &d ) override
127 {
128 SpatialIndex::IShape *shape = nullptr;
129 d.getShape( &shape );
130 mNewIndex->insertData( 0, nullptr, *shape, d.getIdentifier() );
131 delete shape;
132 }
133
134 void visitData( std::vector<const IData *> &v ) override
135 { Q_UNUSED( v ) }
136
137 private:
138 SpatialIndex::ISpatialIndex *mNewIndex = nullptr;
139};
140
141
148class QgsMeshIteratorDataStream : public IDataStream
149{
150 public:
152 explicit QgsMeshIteratorDataStream( const QgsMesh &mesh,
153 int featuresCount,
154 std::function<Region( const QgsMesh &mesh, int id, bool &ok )> featureToRegionFunction,
155 QgsFeedback *feedback = nullptr )
156 : mMesh( mesh )
157 , mFeaturesCount( featuresCount )
158 , mFeatureToRegionFunction( featureToRegionFunction )
159 , mFeedback( feedback )
160 {
161 readNextEntry();
162 }
163
164 ~QgsMeshIteratorDataStream() override
165 {
166 delete mNextData;
167 }
168
170 IData *getNext() override
171 {
172 if ( mFeedback && mFeedback->isCanceled() )
173 return nullptr;
174
175 RTree::Data *ret = mNextData;
176 mNextData = nullptr;
177 readNextEntry();
178 return ret;
179 }
180
182 bool hasNext() override
183 {
184 return nullptr != mNextData;
185 }
186
188 uint32_t size() override
189 {
190 return static_cast<uint32_t>( mFeaturesCount );
191 }
192
194 void rewind() override
195 {
196 mIterator = 0;
197 }
198
199 protected:
200 void readNextEntry()
201 {
202 SpatialIndex::Region r;
203 while ( mIterator < mFeaturesCount )
204 {
205 bool ok = false;
206 r = mFeatureToRegionFunction( mMesh, mIterator, ok );
207 if ( ok )
208 {
209 mNextData = new RTree::Data( 0, nullptr, r, mIterator );
210 ++mIterator;
211 return;
212 }
213 else
214 {
215 ++mIterator;
216 continue;
217 }
218 }
219 }
220
221 private:
222 int mIterator = 0;
223 const QgsMesh &mMesh;
224 int mFeaturesCount = 0;
225 std::function<Region( const QgsMesh &mesh, int id, bool &ok )> mFeatureToRegionFunction;
226 RTree::Data *mNextData = nullptr;
227 QgsFeedback *mFeedback = nullptr;
228};
229
236class QgsMeshSpatialIndexData : public QSharedData
237{
238 public:
239 QgsMeshSpatialIndexData()
240 {
241 initTree();
242 }
243
252 explicit QgsMeshSpatialIndexData( const QgsMesh &fi, QgsFeedback *feedback, QgsMesh::ElementType elementType )
253 {
254 switch ( elementType )
255 {
257 {
258 QgsMeshIteratorDataStream fids( fi, fi.edgeCount(), edgeToRegion, feedback );
259 initTree( &fids );
260 }
261 break;
263 {
264 QgsMeshIteratorDataStream fids( fi, fi.faceCount(), faceToRegion, feedback );
265 initTree( &fids );
266 }
267 break;
268 default:
269 // vertices are not supported
270 Q_ASSERT( false );
271 break;
272 }
273 }
274
275 QgsMeshSpatialIndexData( const QgsMeshSpatialIndexData &other )
276 : QSharedData( other )
277 {
278 const QMutexLocker locker( &other.mMutex );
279
280 initTree();
281
282 // copy R-tree data one by one (is there a faster way??)
283 double low[] = { std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest() };
284 double high[] = { std::numeric_limits<double>::max(), std::numeric_limits<double>::max() };
285 const SpatialIndex::Region query( low, high, 2 );
286 QgsMeshSpatialIndexCopyVisitor visitor( mRTree.get() );
287 other.mRTree->intersectsWithQuery( query, visitor );
288 }
289
290 ~QgsMeshSpatialIndexData() = default;
291
292 QgsMeshSpatialIndexData &operator=( const QgsMeshSpatialIndexData &rh ) = delete;
293
294 void initTree( IDataStream *inputStream = nullptr )
295 {
296 // for now only memory manager
297 mStorage.reset( StorageManager::createNewMemoryStorageManager() );
298
299 // R-Tree parameters
300 const double fillFactor = 0.7;
301 const unsigned int indexCapacity = 10;
302 const unsigned int leafCapacity = 10;
303 const unsigned int dimension = 2;
304 const RTree::RTreeVariant variant = RTree::RV_RSTAR;
305
306 // create R-tree
307 SpatialIndex::id_type indexId;
308
309 if ( inputStream && inputStream->hasNext() )
310 mRTree.reset(
311 RTree::createAndBulkLoadNewRTree(
312 RTree::BLM_STR,
313 *inputStream,
314 *mStorage, fillFactor,
315 indexCapacity,
316 leafCapacity,
317 dimension,
318 variant,
319 indexId )
320 );
321 else
322 mRTree.reset(
323 RTree::createNewRTree(
324 *mStorage,
325 fillFactor,
326 indexCapacity,
327 leafCapacity,
328 dimension,
329 variant,
330 indexId )
331 );
332 }
333
335 std::unique_ptr<SpatialIndex::IStorageManager> mStorage;
336
338 std::unique_ptr<SpatialIndex::ISpatialIndex> mRTree;
339
340 mutable QMutex mMutex;
341};
342
344
346{
347 d = new QgsMeshSpatialIndexData;
348}
349
351 : mElementType( elementType )
352{
353 d = new QgsMeshSpatialIndexData( mesh, feedback, elementType );
354}
355
357 : mElementType( other.mElementType )
358 , d( other.d )
359{
360}
361
363
365{
366 if ( this != &other )
367 {
368 mElementType = other.mElementType;
369 d = other.d;
370 }
371 return *this;
372}
373
374QList<int> QgsMeshSpatialIndex::intersects( const QgsRectangle &rect ) const
375{
376 QList<int> list;
377 QgisMeshVisitor visitor( list );
378
379 const SpatialIndex::Region r = rectToRegion( rect );
380
381 const QMutexLocker locker( &d->mMutex );
382 d->mRTree->intersectsWithQuery( r, visitor );
383
384 return list;
385}
386
387QList<int> QgsMeshSpatialIndex::nearestNeighbor( const QgsPointXY &point, int neighbors ) const
388{
389 QList<int> list;
390 QgisMeshVisitor visitor( list );
391
392 double pt[2] = { point.x(), point.y() };
393 const Point p( pt, 2 );
394
395 const QMutexLocker locker( &d->mMutex );
396 d->mRTree->nearestNeighborQuery( static_cast<uint32_t>( neighbors ), p, visitor );
397
398 return list;
399}
400
402{
403 return mElementType;
404}
405
406void QgsMeshSpatialIndex::addFace( int faceIndex, const QgsMesh &mesh )
407{
408 if ( mesh.face( faceIndex ).isEmpty() )
409 return;
410
411 bool ok = false;
412 const SpatialIndex::Region r( faceToRegion( mesh, faceIndex, ok ) );
413 if ( !ok )
414 return;
415
416 const QMutexLocker locker( &d.constData()->mMutex );
417
418 try
419 {
420 d.constData()->mRTree->insertData( 0, nullptr, r, faceIndex );
421 }
422 catch ( Tools::Exception &e )
423 {
424 Q_UNUSED( e )
425 QgsDebugError( QStringLiteral( "Tools::Exception caught: " ).arg( e.what().c_str() ) );
426 }
427 catch ( const std::exception &e )
428 {
429 Q_UNUSED( e )
430 QgsDebugError( QStringLiteral( "std::exception caught: " ).arg( e.what() ) );
431 }
432 catch ( ... )
433 {
434 QgsDebugError( QStringLiteral( "unknown spatial index exception caught" ) );
435 }
436}
437
438void QgsMeshSpatialIndex::removeFace( int faceIndex, const QgsMesh &mesh )
439{
440 if ( mesh.face( faceIndex ).isEmpty() )
441 return;
442 const QMutexLocker locker( &d.constData()->mMutex );
443 bool ok = false;
444 d.constData()->mRTree->deleteData( faceToRegion( mesh, faceIndex, ok ), faceIndex );
445}
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition qgsfeedback.h:44
A spatial index for QgsMeshFace or QgsMeshEdge objects.
QList< int > intersects(const QgsRectangle &rectangle) const
Returns a list of face ids with a bounding box which intersects the specified rectangle.
QgsMesh::ElementType elementType() const
Returns the type of mesh elements that are indexed.
QgsMeshSpatialIndex()
Constructor for QgsSpatialIndex.
void addFace(int faceIndex, const QgsMesh &mesh)
Adds a face with faceIndex from the mesh in the spatial index.
void removeFace(int faceIndex, const QgsMesh &mesh)
Removes a face with faceIndex from the mesh in the spatial index.
QgsMeshSpatialIndex & operator=(const QgsMeshSpatialIndex &other)
Implement assignment operator.
QList< int > nearestNeighbor(const QgsPointXY &point, int neighbors) const
Returns nearest neighbors to a point.
~QgsMeshSpatialIndex()
Destructor finalizes work with spatial index.
A class to represent a 2D point.
Definition qgspointxy.h:60
double y
Definition qgspointxy.h:64
double x
Definition qgspointxy.h:63
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
double x
Definition qgspoint.h:52
double y
Definition qgspoint.h:53
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
#define QgsDebugError(str)
Definition qgslogger.h:38
QVector< int > QgsMeshFace
List of vertex indexes.
QPair< int, int > QgsMeshEdge
Edge is a straight line seqment between 2 points.
Mesh - vertices, edges and faces.
QVector< QgsMeshVertex > vertices
QgsMeshFace face(int index) const
Returns a face at the index.
int faceCount() const
Returns number of faces.
ElementType
Defines type of mesh elements.
QgsMeshEdge edge(int index) const
Returns an edge at the index.
int edgeCount() const
Returns number of edge.