QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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 
18 #include "qgsmeshdataprovider.h"
19 #include "qgsrectangle.h"
20 #include "qgslogger.h"
21 #include "qgsfeedback.h"
22 
23 #include <spatialindex/SpatialIndex.h>
24 #include <QMutex>
25 #include <QMutexLocker>
26 #include <memory>
27 
28 using namespace SpatialIndex;
29 
31 
32 static Region faceToRegion( const QgsMesh &mesh, int id )
33 {
34  const QgsMeshFace face = mesh.face( id );
35  const QVector<QgsMeshVertex> &vertices = mesh.vertices;
36  Q_ASSERT( face.size() > 0 );
37  double xMinimum = vertices[face[0]].x();
38  double yMinimum = vertices[face[0]].y();
39  double xMaximum = vertices[face[0]].x();
40  double yMaximum = vertices[face[0]].y();
41 
42  for ( int i = 1; i < face.size(); ++i )
43  {
44  xMinimum = std::min( vertices[face[i]].x(), xMinimum );
45  yMinimum = std::min( vertices[face[i]].y(), yMinimum );
46  xMaximum = std::max( vertices[face[i]].x(), xMaximum );
47  yMaximum = std::max( vertices[face[i]].y(), yMaximum );
48  }
49 
50  double pt1[2] = { xMinimum, yMinimum };
51  double pt2[2] = { xMaximum, yMaximum };
52  return SpatialIndex::Region( pt1, pt2, 2 );
53 }
54 
55 static Region rectToRegion( const QgsRectangle &rect )
56 {
57  double pt1[2] = { rect.xMinimum(), rect.yMinimum() };
58  double pt2[2] = { rect.xMaximum(), rect.yMaximum() };
59  return SpatialIndex::Region( pt1, pt2, 2 );
60 }
61 
68 class QgisMeshVisitor : public SpatialIndex::IVisitor
69 {
70  public:
71  explicit QgisMeshVisitor( QList<int> &list )
72  : mList( list ) {}
73 
74  void visitNode( const INode &n ) override
75  { Q_UNUSED( n ) }
76 
77  void visitData( const IData &d ) override
78  {
79  mList.append( static_cast<int>( d.getIdentifier() ) );
80  }
81 
82  void visitData( std::vector<const IData *> &v ) override
83  { Q_UNUSED( v ) }
84 
85  private:
86  QList<int> &mList;
87 };
88 
94 class QgsMeshSpatialIndexCopyVisitor : public SpatialIndex::IVisitor
95 {
96  public:
97  explicit QgsMeshSpatialIndexCopyVisitor( SpatialIndex::ISpatialIndex *newIndex )
98  : mNewIndex( newIndex ) {}
99 
100  void visitNode( const INode &n ) override
101  { Q_UNUSED( n ) }
102 
103  void visitData( const IData &d ) override
104  {
105  SpatialIndex::IShape *shape = nullptr;
106  d.getShape( &shape );
107  mNewIndex->insertData( 0, nullptr, *shape, d.getIdentifier() );
108  delete shape;
109  }
110 
111  void visitData( std::vector<const IData *> &v ) override
112  { Q_UNUSED( v ) }
113 
114  private:
115  SpatialIndex::ISpatialIndex *mNewIndex = nullptr;
116 };
117 
118 
125 class QgsMeshFaceIteratorDataStream : public IDataStream
126 {
127  public:
129  explicit QgsMeshFaceIteratorDataStream( const QgsMesh &triangularMesh, QgsFeedback *feedback = nullptr )
130  : mMesh( triangularMesh )
131  , mFeedback( feedback )
132  {
133  readNextEntry();
134  }
135 
136  ~QgsMeshFaceIteratorDataStream() override
137  {
138  delete mNextData;
139  }
140 
142  IData *getNext() override
143  {
144  if ( mFeedback && mFeedback->isCanceled() )
145  return nullptr;
146 
147  RTree::Data *ret = mNextData;
148  mNextData = nullptr;
149  readNextEntry();
150  return ret;
151  }
152 
154  bool hasNext() override
155  {
156  return nullptr != mNextData;
157  }
158 
160  uint32_t size() override
161  {
162  return static_cast<uint32_t>( mMesh.faceCount() );
163  }
164 
166  void rewind() override
167  {
168  mIterator = 0;
169  }
170 
171  protected:
172  void readNextEntry()
173  {
174  SpatialIndex::Region r;
175  const int faceCount = mMesh.faceCount();
176  if ( mIterator < faceCount )
177  {
178  r = faceToRegion( mMesh, mIterator );
179  mNextData = new RTree::Data(
180  0,
181  nullptr,
182  r,
183  mIterator );
184  ++mIterator;
185  }
186  }
187 
188  private:
189  int mIterator = 0;
190  const QgsMesh &mMesh;
191  RTree::Data *mNextData = nullptr;
192  QgsFeedback *mFeedback = nullptr;
193 };
194 
195 
202 class QgsMeshSpatialIndexData : public QSharedData
203 {
204  public:
205  QgsMeshSpatialIndexData()
206  {
207  initTree();
208  }
209 
218  explicit QgsMeshSpatialIndexData( const QgsMesh &fi, QgsFeedback *feedback = nullptr )
219  {
220  QgsMeshFaceIteratorDataStream fids( fi, feedback );
221  initTree( &fids );
222  }
223 
224  QgsMeshSpatialIndexData( const QgsMeshSpatialIndexData &other )
225  : QSharedData( other )
226  {
227  QMutexLocker locker( &other.mMutex );
228 
229  initTree();
230 
231  // copy R-tree data one by one (is there a faster way??)
232  double low[] = { std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest() };
233  double high[] = { std::numeric_limits<double>::max(), std::numeric_limits<double>::max() };
234  SpatialIndex::Region query( low, high, 2 );
235  QgsMeshSpatialIndexCopyVisitor visitor( mRTree.get() );
236  other.mRTree->intersectsWithQuery( query, visitor );
237  }
238 
239  ~QgsMeshSpatialIndexData() = default;
240 
241  QgsMeshSpatialIndexData &operator=( const QgsMeshSpatialIndexData &rh ) = delete;
242 
243  void initTree( IDataStream *inputStream = nullptr )
244  {
245  // for now only memory manager
246  mStorage.reset( StorageManager::createNewMemoryStorageManager() );
247 
248  // R-Tree parameters
249  double fillFactor = 0.7;
250  unsigned int indexCapacity = 10;
251  unsigned int leafCapacity = 10;
252  unsigned int dimension = 2;
253  RTree::RTreeVariant variant = RTree::RV_RSTAR;
254 
255  // create R-tree
256  SpatialIndex::id_type indexId;
257 
258  if ( inputStream && inputStream->hasNext() )
259  mRTree.reset(
260  RTree::createAndBulkLoadNewRTree(
261  RTree::BLM_STR,
262  *inputStream,
263  *mStorage, fillFactor,
264  indexCapacity,
265  leafCapacity,
266  dimension,
267  variant,
268  indexId )
269  );
270  else
271  mRTree.reset(
272  RTree::createNewRTree(
273  *mStorage,
274  fillFactor,
275  indexCapacity,
276  leafCapacity,
277  dimension,
278  variant,
279  indexId )
280  );
281  }
282 
284  std::unique_ptr<SpatialIndex::IStorageManager> mStorage;
285 
287  std::unique_ptr<SpatialIndex::ISpatialIndex> mRTree;
288 
289  mutable QMutex mMutex;
290 
291 };
292 
294 
296 {
297  d = new QgsMeshSpatialIndexData;
298 }
299 
301 {
302  d = new QgsMeshSpatialIndexData( triangularMesh, feedback );
303 }
304 
306  : d( other.d )
307 {
308 }
309 
310 QgsMeshSpatialIndex:: ~QgsMeshSpatialIndex() = default; //NOLINT
311 
313 {
314  if ( this != &other )
315  d = other.d;
316  return *this;
317 }
318 
319 QList<int> QgsMeshSpatialIndex::intersects( const QgsRectangle &rect ) const
320 {
321  QList<int> list;
322  QgisMeshVisitor visitor( list );
323 
324  SpatialIndex::Region r = rectToRegion( rect );
325 
326  QMutexLocker locker( &d->mMutex );
327  d->mRTree->intersectsWithQuery( r, visitor );
328 
329  return list;
330 }
331 
332 QList<int> QgsMeshSpatialIndex::nearestNeighbor( const QgsPointXY &point, int neighbors ) const
333 {
334  QList<int> list;
335  QgisMeshVisitor visitor( list );
336 
337  double pt[2] = { point.x(), point.y() };
338  Point p( pt, 2 );
339 
340  QMutexLocker locker( &d->mMutex );
341  d->mRTree->nearestNeighborQuery( static_cast<uint32_t>( neighbors ), p, visitor );
342 
343  return list;
344 }
A rectangle specified with double values.
Definition: qgsrectangle.h:41
QgsMeshSpatialIndex()
Constructor for QgsSpatialIndex.
QgsMeshSpatialIndex & operator=(const QgsMeshSpatialIndex &other)
Implement assignment operator.
double y
Definition: qgspointxy.h:48
A class to represent a 2D point.
Definition: qgspointxy.h:43
QList< int > intersects(const QgsRectangle &rectangle) const
Returns a list of face ids with a bounding box which intersects the specified rectangle.
QVector< QgsMeshVertex > vertices
vertices
Base class for feedback objects to be used for cancellation of something running in a worker thread...
Definition: qgsfeedback.h:44
QList< int > nearestNeighbor(const QgsPointXY &point, int neighbors) const
Returns nearest neighbors to a point.
double x
Definition: qgspointxy.h:47
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:177
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:162
A spatial index for QgsMeshFace objects.
Mesh - vertices and faces.
QVector< int > QgsMeshFace
List of vertex indexes.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
~QgsMeshSpatialIndex()
Destructor finalizes work with spatial index.
QgsMeshFace face(int index) const
Returns a face at the index.
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:172