Branch data Line data Source code
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 : :
26 : : using namespace SpatialIndex;
27 : :
28 : : ///@cond PRIVATE
29 : :
30 : 0 : static Region faceToRegion( const QgsMesh &mesh, int id )
31 : : {
32 : 0 : const QgsMeshFace face = mesh.face( id );
33 : 0 : const QVector<QgsMeshVertex> &vertices = mesh.vertices;
34 : : Q_ASSERT( face.size() > 0 );
35 : 0 : double xMinimum = vertices[face[0]].x();
36 : 0 : double yMinimum = vertices[face[0]].y();
37 : 0 : double xMaximum = vertices[face[0]].x();
38 : 0 : double yMaximum = vertices[face[0]].y();
39 : :
40 : 0 : for ( int i = 1; i < face.size(); ++i )
41 : : {
42 : 0 : xMinimum = std::min( vertices[face[i]].x(), xMinimum );
43 : 0 : yMinimum = std::min( vertices[face[i]].y(), yMinimum );
44 : 0 : xMaximum = std::max( vertices[face[i]].x(), xMaximum );
45 : 0 : yMaximum = std::max( vertices[face[i]].y(), yMaximum );
46 : 0 : }
47 : :
48 : 0 : double pt1[2] = { xMinimum, yMinimum };
49 : 0 : double pt2[2] = { xMaximum, yMaximum };
50 : 0 : return SpatialIndex::Region( pt1, pt2, 2 );
51 : 0 : }
52 : :
53 : 0 : static Region edgeToRegion( const QgsMesh &mesh, int id )
54 : : {
55 : 0 : const QgsMeshEdge edge = mesh.edge( id );
56 : 0 : const QgsMeshVertex firstVertex = mesh.vertices[edge.first];
57 : 0 : const QgsMeshVertex secondVertex = mesh.vertices[edge.second];
58 : 0 : double xMinimum = std::min( firstVertex.x(), secondVertex.x() );
59 : 0 : double yMinimum = std::min( firstVertex.y(), secondVertex.y() );
60 : 0 : double xMaximum = std::max( firstVertex.x(), secondVertex.x() );
61 : 0 : double yMaximum = std::max( firstVertex.y(), secondVertex.y() );
62 : 0 : double pt1[2] = { xMinimum, yMinimum };
63 : 0 : double pt2[2] = { xMaximum, yMaximum };
64 : 0 : return SpatialIndex::Region( pt1, pt2, 2 );
65 : 0 : }
66 : :
67 : 0 : static Region rectToRegion( const QgsRectangle &rect )
68 : : {
69 : 0 : double pt1[2] = { rect.xMinimum(), rect.yMinimum() };
70 : 0 : double pt2[2] = { rect.xMaximum(), rect.yMaximum() };
71 : 0 : return SpatialIndex::Region( pt1, pt2, 2 );
72 : : }
73 : :
74 : : /**
75 : : * \ingroup core
76 : : * \class QgisMeshVisitor
77 : : * \brief Custom visitor that adds found faces or edges to list.
78 : : * \note not available in Python bindings
79 : : */
80 : 0 : class QgisMeshVisitor : public SpatialIndex::IVisitor
81 : : {
82 : : public:
83 : 0 : explicit QgisMeshVisitor( QList<int> &list )
84 : 0 : : mList( list ) {}
85 : :
86 : 0 : void visitNode( const INode &n ) override
87 : 0 : { Q_UNUSED( n ) }
88 : :
89 : 0 : void visitData( const IData &d ) override
90 : : {
91 : 0 : mList.append( static_cast<int>( d.getIdentifier() ) );
92 : 0 : }
93 : :
94 : 0 : void visitData( std::vector<const IData *> &v ) override
95 : 0 : { Q_UNUSED( v ) }
96 : :
97 : : private:
98 : : QList<int> &mList;
99 : : };
100 : :
101 : : /**
102 : : * \ingroup core
103 : : * \class QgsMeshSpatialIndexCopyVisitor
104 : 0 : * \brief A copy visitor for populating a spatial index.
105 : : * \note not available in Python bindings
106 : : */
107 : : class QgsMeshSpatialIndexCopyVisitor : public SpatialIndex::IVisitor
108 : : {
109 : : public:
110 : : explicit QgsMeshSpatialIndexCopyVisitor( SpatialIndex::ISpatialIndex *newIndex )
111 : : : mNewIndex( newIndex ) {}
112 : :
113 : : void visitNode( const INode &n ) override
114 : : { Q_UNUSED( n ) }
115 : :
116 : : void visitData( const IData &d ) override
117 : : {
118 : : SpatialIndex::IShape *shape = nullptr;
119 : : d.getShape( &shape );
120 : : mNewIndex->insertData( 0, nullptr, *shape, d.getIdentifier() );
121 : : delete shape;
122 : : }
123 : :
124 : : void visitData( std::vector<const IData *> &v ) override
125 : : { Q_UNUSED( v ) }
126 : :
127 : : private:
128 : : SpatialIndex::ISpatialIndex *mNewIndex = nullptr;
129 : : };
130 : :
131 : :
132 : : /**
133 : : * \ingroup core
134 : : * \class QgsMeshFaceIteratorDataStream
135 : : * \brief Utility class for bulk loading of R-trees. Not a part of public API.
136 : : * \note not available in Python bindings
137 : : */
138 : : class QgsMeshIteratorDataStream : public IDataStream
139 : : {
140 : : public:
141 : : //! constructor - needs to load all data to a vector for later access when bulk loading
142 : 0 : explicit QgsMeshIteratorDataStream( const QgsMesh &mesh,
143 : : int featuresCount,
144 : : std::function<Region( const QgsMesh &mesh, int id )> featureToRegionFunction,
145 : : QgsFeedback *feedback = nullptr )
146 : 0 : : mMesh( mesh )
147 : 0 : , mFeaturesCount( featuresCount )
148 : 0 : , mFeatureToRegionFunction( featureToRegionFunction )
149 : 0 : , mFeedback( feedback )
150 : 0 : {
151 : 0 : readNextEntry();
152 : 0 : }
153 : :
154 : 0 : ~QgsMeshIteratorDataStream() override
155 : 0 : {
156 : 0 : delete mNextData;
157 : 0 : }
158 : :
159 : : //! returns a pointer to the next entry in the stream or 0 at the end of the stream.
160 : 0 : IData *getNext() override
161 : : {
162 : 0 : if ( mFeedback && mFeedback->isCanceled() )
163 : 0 : return nullptr;
164 : :
165 : 0 : RTree::Data *ret = mNextData;
166 : 0 : mNextData = nullptr;
167 : 0 : readNextEntry();
168 : 0 : return ret;
169 : 0 : }
170 : :
171 : : //! returns true if there are more items in the stream.
172 : 0 : bool hasNext() override
173 : : {
174 : 0 : return nullptr != mNextData;
175 : : }
176 : :
177 : : //! returns the total number of entries available in the stream.
178 : 0 : uint32_t size() override
179 : : {
180 : 0 : return static_cast<uint32_t>( mFeaturesCount );
181 : : }
182 : :
183 : : //! sets the stream pointer to the first entry, if possible.
184 : 0 : void rewind() override
185 : : {
186 : 0 : mIterator = 0;
187 : 0 : }
188 : :
189 : : protected:
190 : 0 : void readNextEntry()
191 : : {
192 : 0 : SpatialIndex::Region r;
193 : 0 : if ( mIterator < mFeaturesCount )
194 : : {
195 : 0 : r = mFeatureToRegionFunction( mMesh, mIterator );
196 : 0 : mNextData = new RTree::Data(
197 : : 0,
198 : : nullptr,
199 : : r,
200 : 0 : mIterator );
201 : 0 : ++mIterator;
202 : 0 : }
203 : 0 : }
204 : :
205 : : private:
206 : 0 : int mIterator = 0;
207 : : const QgsMesh &mMesh;
208 : : int mFeaturesCount = 0;
209 : : std::function<Region( const QgsMesh &mesh, int id )> mFeatureToRegionFunction;
210 : 0 : RTree::Data *mNextData = nullptr;
211 : : QgsFeedback *mFeedback = nullptr;
212 : : };
213 : :
214 : : /**
215 : : * \ingroup core
216 : : * \class QgsSpatialIndexData
217 : : * \brief Data of spatial index that may be implicitly shared
218 : : * \note not available in Python bindings
219 : : */
220 : : class QgsMeshSpatialIndexData : public QSharedData
221 : : {
222 : : public:
223 : 0 : QgsMeshSpatialIndexData()
224 : 0 : {
225 : 0 : initTree();
226 : 0 : }
227 : :
228 : : /**
229 : : * Constructor for QgsSpatialIndexData which bulk loads faces from the specified mesh
230 : : * \a fi.
231 : : *
232 : : * The optional \a feedback object can be used to allow cancellation of bulk face loading. Ownership
233 : : * of \a feedback is not transferred, and callers must take care that the lifetime of feedback exceeds
234 : : * that of the spatial index construction.
235 : : */
236 : 0 : explicit QgsMeshSpatialIndexData( const QgsMesh &fi, QgsFeedback *feedback, QgsMesh::ElementType elementType )
237 : 0 : {
238 : 0 : switch ( elementType )
239 : : {
240 : : case QgsMesh::ElementType::Edge:
241 : : {
242 : 0 : QgsMeshIteratorDataStream fids( fi, fi.edgeCount(), edgeToRegion, feedback );
243 : 0 : initTree( &fids );
244 : 0 : }
245 : 0 : break;
246 : : case QgsMesh::ElementType::Face:
247 : : {
248 : 0 : QgsMeshIteratorDataStream fids( fi, fi.faceCount(), faceToRegion, feedback );
249 : 0 : initTree( &fids );
250 : 0 : }
251 : 0 : break;
252 : : default:
253 : : // vertices are not supported
254 : : Q_ASSERT( false );
255 : 0 : break;
256 : : }
257 : 0 : }
258 : :
259 : : QgsMeshSpatialIndexData( const QgsMeshSpatialIndexData &other )
260 : : : QSharedData( other )
261 : : {
262 : : QMutexLocker locker( &other.mMutex );
263 : :
264 : : initTree();
265 : :
266 : : // copy R-tree data one by one (is there a faster way??)
267 : : double low[] = { std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest() };
268 : : double high[] = { std::numeric_limits<double>::max(), std::numeric_limits<double>::max() };
269 : : SpatialIndex::Region query( low, high, 2 );
270 : : QgsMeshSpatialIndexCopyVisitor visitor( mRTree.get() );
271 : : other.mRTree->intersectsWithQuery( query, visitor );
272 : : }
273 : :
274 : 0 : ~QgsMeshSpatialIndexData() = default;
275 : :
276 : : QgsMeshSpatialIndexData &operator=( const QgsMeshSpatialIndexData &rh ) = delete;
277 : :
278 : 0 : void initTree( IDataStream *inputStream = nullptr )
279 : : {
280 : : // for now only memory manager
281 : 0 : mStorage.reset( StorageManager::createNewMemoryStorageManager() );
282 : :
283 : : // R-Tree parameters
284 : 0 : double fillFactor = 0.7;
285 : 0 : unsigned int indexCapacity = 10;
286 : 0 : unsigned int leafCapacity = 10;
287 : 0 : unsigned int dimension = 2;
288 : 0 : RTree::RTreeVariant variant = RTree::RV_RSTAR;
289 : :
290 : : // create R-tree
291 : : SpatialIndex::id_type indexId;
292 : :
293 : 0 : if ( inputStream && inputStream->hasNext() )
294 : 0 : mRTree.reset(
295 : 0 : RTree::createAndBulkLoadNewRTree(
296 : : RTree::BLM_STR,
297 : 0 : *inputStream,
298 : 0 : *mStorage, fillFactor,
299 : 0 : indexCapacity,
300 : 0 : leafCapacity,
301 : 0 : dimension,
302 : 0 : variant,
303 : : indexId )
304 : : );
305 : : else
306 : 0 : mRTree.reset(
307 : 0 : RTree::createNewRTree(
308 : 0 : *mStorage,
309 : 0 : fillFactor,
310 : 0 : indexCapacity,
311 : 0 : leafCapacity,
312 : 0 : dimension,
313 : 0 : variant,
314 : : indexId )
315 : : );
316 : 0 : }
317 : :
318 : : //! Storage manager
319 : : std::unique_ptr<SpatialIndex::IStorageManager> mStorage;
320 : :
321 : : //! R-tree containing spatial index
322 : : std::unique_ptr<SpatialIndex::ISpatialIndex> mRTree;
323 : :
324 : : mutable QMutex mMutex;
325 : : };
326 : :
327 : : ///@endcond
328 : :
329 : 0 : QgsMeshSpatialIndex::QgsMeshSpatialIndex()
330 : : {
331 : 0 : d = new QgsMeshSpatialIndexData;
332 : 0 : }
333 : :
334 : 0 : QgsMeshSpatialIndex::QgsMeshSpatialIndex( const QgsMesh &mesh, QgsFeedback *feedback, QgsMesh::ElementType elementType )
335 : 0 : : mElementType( elementType )
336 : : {
337 : 0 : d = new QgsMeshSpatialIndexData( mesh, feedback, elementType );
338 : 0 : }
339 : :
340 : 0 : QgsMeshSpatialIndex::QgsMeshSpatialIndex( const QgsMeshSpatialIndex &other ) //NOLINT
341 : 0 : : d( other.d )
342 : : {
343 : 0 : }
344 : :
345 : 0 : QgsMeshSpatialIndex:: ~QgsMeshSpatialIndex() = default; //NOLINT
346 : :
347 : 0 : QgsMeshSpatialIndex &QgsMeshSpatialIndex::operator=( const QgsMeshSpatialIndex &other )
348 : : {
349 : 0 : if ( this != &other )
350 : 0 : d = other.d;
351 : 0 : return *this;
352 : : }
353 : :
354 : 0 : QList<int> QgsMeshSpatialIndex::intersects( const QgsRectangle &rect ) const
355 : : {
356 : 0 : QList<int> list;
357 : 0 : QgisMeshVisitor visitor( list );
358 : :
359 : 0 : SpatialIndex::Region r = rectToRegion( rect );
360 : :
361 : 0 : QMutexLocker locker( &d->mMutex );
362 : 0 : d->mRTree->intersectsWithQuery( r, visitor );
363 : :
364 : 0 : return list;
365 : 0 : }
366 : :
367 : 0 : QList<int> QgsMeshSpatialIndex::nearestNeighbor( const QgsPointXY &point, int neighbors ) const
368 : : {
369 : 0 : QList<int> list;
370 : 0 : QgisMeshVisitor visitor( list );
371 : :
372 : 0 : double pt[2] = { point.x(), point.y() };
373 : 0 : Point p( pt, 2 );
374 : :
375 : 0 : QMutexLocker locker( &d->mMutex );
376 : 0 : d->mRTree->nearestNeighborQuery( static_cast<uint32_t>( neighbors ), p, visitor );
377 : :
378 : 0 : return list;
379 : 0 : }
380 : :
381 : 0 : QgsMesh::ElementType QgsMeshSpatialIndex::elementType() const
382 : : {
383 : 0 : return mElementType;
384 : : }
|
