QGIS API Documentation  3.23.0-Master (dd0cd13a00)
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 
26 using namespace SpatialIndex;
27 
29 
30 static 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 
62 static 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 
77 static 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 
90 class 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 
117 class 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 
148 class 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 
236 class 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  {
256  case QgsMesh::ElementType::Edge:
257  {
258  QgsMeshIteratorDataStream fids( fi, fi.edgeCount(), edgeToRegion, feedback );
259  initTree( &fids );
260  }
261  break;
262  case QgsMesh::ElementType::Face:
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  : d( other.d )
358 {
359 }
360 
361 QgsMeshSpatialIndex:: ~QgsMeshSpatialIndex() = default; //NOLINT
362 
364 {
365  if ( this != &other )
366  d = other.d;
367  return *this;
368 }
369 
370 QList<int> QgsMeshSpatialIndex::intersects( const QgsRectangle &rect ) const
371 {
372  QList<int> list;
373  QgisMeshVisitor visitor( list );
374 
375  const SpatialIndex::Region r = rectToRegion( rect );
376 
377  const QMutexLocker locker( &d->mMutex );
378  d->mRTree->intersectsWithQuery( r, visitor );
379 
380  return list;
381 }
382 
383 QList<int> QgsMeshSpatialIndex::nearestNeighbor( const QgsPointXY &point, int neighbors ) const
384 {
385  QList<int> list;
386  QgisMeshVisitor visitor( list );
387 
388  double pt[2] = { point.x(), point.y() };
389  const Point p( pt, 2 );
390 
391  const QMutexLocker locker( &d->mMutex );
392  d->mRTree->nearestNeighborQuery( static_cast<uint32_t>( neighbors ), p, visitor );
393 
394  return list;
395 }
396 
398 {
399  return mElementType;
400 }
401 
402 void QgsMeshSpatialIndex::addFace( int faceIndex, const QgsMesh &mesh )
403 {
404  if ( mesh.face( faceIndex ).isEmpty() )
405  return;
406 
407  bool ok = false;
408  const SpatialIndex::Region r( faceToRegion( mesh, faceIndex, ok ) );
409  if ( !ok )
410  return;
411 
412  const QMutexLocker locker( &d.constData()->mMutex );
413 
414  try
415  {
416  d.constData()->mRTree->insertData( 0, nullptr, r, faceIndex );
417  }
418  catch ( Tools::Exception &e )
419  {
420  Q_UNUSED( e )
421  QgsDebugMsg( QStringLiteral( "Tools::Exception caught: " ).arg( e.what().c_str() ) );
422  }
423  catch ( const std::exception &e )
424  {
425  Q_UNUSED( e )
426  QgsDebugMsg( QStringLiteral( "std::exception caught: " ).arg( e.what() ) );
427  }
428  catch ( ... )
429  {
430  QgsDebugMsg( QStringLiteral( "unknown spatial index exception caught" ) );
431  }
432 }
433 
434 void QgsMeshSpatialIndex::removeFace( int faceIndex, const QgsMesh &mesh )
435 {
436  if ( mesh.face( faceIndex ).isEmpty() )
437  return;
438  const QMutexLocker locker( &d.constData()->mMutex );
439  bool ok = false;
440  d.constData()->mRTree->deleteData( faceToRegion( mesh, faceIndex, ok ), faceIndex );
441 }
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
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:59
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
Q_GADGET double x
Definition: qgspoint.h:52
double y
Definition: qgspoint.h:53
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
#define QgsDebugMsg(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.