improve speed of creation spatial index for meshes

qgis · Jan 24, 2019 · 9d4d1c9 · 9d4d1c9
1 parent a83a1b2
commit 9d4d1c9
Show file tree

Hide file tree

Showing 5 changed files with 464 additions and 20 deletions.
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
@@ -473,13 +473,14 @@ SET(QGIS_CORE_SRCS
 
   mesh/qgsmeshdataprovider.cpp
   mesh/qgsmeshlayer.cpp
+  mesh/qgsmeshlayerinterpolator.cpp
   mesh/qgsmeshlayerrenderer.cpp
   mesh/qgsmeshlayerutils.cpp
   mesh/qgsmeshmemorydataprovider.cpp
-  mesh/qgstriangularmesh.cpp
-  mesh/qgsmeshlayerinterpolator.cpp
-  mesh/qgsmeshvectorrenderer.cpp
   mesh/qgsmeshrenderersettings.cpp
+  mesh/qgsmeshspatialindex.cpp
+  mesh/qgsmeshvectorrenderer.cpp
+  mesh/qgstriangularmesh.cpp
 
   geometry/qgsabstractgeometry.cpp
   geometry/qgsbox3d.cpp
@@ -1104,12 +1105,13 @@ SET(QGIS_CORE_HDRS
   raster/qgssinglebandpseudocolorrenderer.h
   raster/qgshillshaderenderer.h
 
-  mesh/qgstriangularmesh.h
-  mesh/qgsmeshlayerrenderer.h
   mesh/qgsmeshlayerinterpolator.h
+  mesh/qgsmeshlayerrenderer.h
   mesh/qgsmeshlayerutils.h
-  mesh/qgsmeshvectorrenderer.h
   mesh/qgsmeshrenderersettings.h
+  mesh/qgsmeshspatialindex.h
+  mesh/qgsmeshvectorrenderer.h
+  mesh/qgstriangularmesh.h
 
   scalebar/qgsdoubleboxscalebarrenderer.h
   scalebar/qgsnumericscalebarrenderer.h

diff --git a/src/core/mesh/qgsmeshspatialindex.cpp b/src/core/mesh/qgsmeshspatialindex.cpp
@@ -0,0 +1,346 @@
+/***************************************************************************
+    qgsmeshspatialindex.cpp
+    -----------------------
+    begin                : January 2019
+    copyright            : (C) 2019 by Peter Petrik
+    email                : zilolv at gmail dot com
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+///@cond PRIVATE
+
+#include "qgsmeshspatialindex.h"
+
+#include "qgsmeshdataprovider.h"
+#include "qgsrectangle.h"
+#include "qgslogger.h"
+#include "qgsfeedback.h"
+
+#include <spatialindex/SpatialIndex.h>
+#include <QMutex>
+#include <QMutexLocker>
+#include <memory>
+
+using namespace SpatialIndex;
+
+static Region faceToRegion( const QgsMesh &mesh, int id )
+{
+  const QgsMeshFace face = mesh.face( id );
+  const QVector<QgsMeshVertex> &vertices = mesh.vertices;
+  Q_ASSERT( face.size() > 0 );
+  double xMinimum = vertices[face[0]].x();
+  double yMinimum = vertices[face[0]].y();
+  double xMaximum = vertices[face[0]].x();
+  double yMaximum = vertices[face[0]].y();
+
+  for ( int i = 1; i < face.size(); ++i )
+  {
+    xMinimum = std::min( vertices[face[i]].x(), xMinimum );
+    yMinimum = std::min( vertices[face[i]].y(), yMinimum );
+    xMaximum = std::max( vertices[face[i]].x(), xMaximum );
+    yMaximum = std::max( vertices[face[i]].y(), yMaximum );
+  }
+
+  double pt1[2] = { xMinimum, yMinimum };
+  double pt2[2] = { xMaximum, yMaximum };
+  return SpatialIndex::Region( pt1, pt2, 2 );
+}
+
+static Region rectToRegion( const QgsRectangle &rect )
+{
+  double pt1[2] = { rect.xMinimum(), rect.yMinimum() };
+  double pt2[2] = { rect.xMaximum(), rect.yMaximum() };
+  return SpatialIndex::Region( pt1, pt2, 2 );
+}
+
+/**
+ * \ingroup core
+ * \class QgisMeshVisitor
+ * \brief Custom visitor that adds found faces to list.
+ * \note not available in Python bindings
+ */
+class QgisMeshVisitor : public SpatialIndex::IVisitor
+{
+  public:
+    explicit QgisMeshVisitor( QList<int> &list )
+      : mList( list ) {}
+
+    void visitNode( const INode &n ) override
+    { Q_UNUSED( n ); }
+
+    void visitData( const IData &d ) override
+    {
+      mList.append( static_cast<int>( d.getIdentifier() ) );
+    }
+
+    void visitData( std::vector<const IData *> &v ) override
+    { Q_UNUSED( v ); }
+
+  private:
+    QList<int> &mList;
+};
+
+/**
+ * \ingroup core
+ * \class QgsMeshSpatialIndexCopyVisitor
+ * \note not available in Python bindings
+ */
+class QgsMeshSpatialIndexCopyVisitor : public SpatialIndex::IVisitor
+{
+  public:
+    explicit QgsMeshSpatialIndexCopyVisitor( SpatialIndex::ISpatialIndex *newIndex )
+      : mNewIndex( newIndex ) {}
+
+    void visitNode( const INode &n ) override
+    { Q_UNUSED( n ); }
+
+    void visitData( const IData &d ) override
+    {
+      SpatialIndex::IShape *shape = nullptr;
+      d.getShape( &shape );
+      mNewIndex->insertData( 0, nullptr, *shape, d.getIdentifier() );
+      delete shape;
+    }
+
+    void visitData( std::vector<const IData *> &v ) override
+    { Q_UNUSED( v ); }
+
+  private:
+    SpatialIndex::ISpatialIndex *mNewIndex = nullptr;
+};
+
+
+/**
+ * \ingroup core
+ * \class QgsMeshFaceIteratorDataStream
+ * \brief Utility class for bulk loading of R-trees. Not a part of public API.
+ * \note not available in Python bindings
+*/
+class QgsMeshFaceIteratorDataStream : public IDataStream
+{
+  public:
+    //! constructor - needs to load all data to a vector for later access when bulk loading
+    explicit QgsMeshFaceIteratorDataStream( const QgsMesh &triangularMesh, QgsFeedback *feedback = nullptr )
+      : mMesh( triangularMesh )
+      , mFeedback( feedback )
+    {
+      readNextEntry();
+    }
+
+    ~QgsMeshFaceIteratorDataStream() override
+    {
+      delete mNextData;
+    }
+
+    //! returns a pointer to the next entry in the stream or 0 at the end of the stream.
+    IData *getNext() override
+    {
+      if ( mFeedback && mFeedback->isCanceled() )
+        return nullptr;
+
+      RTree::Data *ret = mNextData;
+      mNextData = nullptr;
+      readNextEntry();
+      return ret;
+    }
+
+    //! returns true if there are more items in the stream.
+    bool hasNext() override
+    {
+      return nullptr != mNextData;
+    }
+
+    //! returns the total number of entries available in the stream.
+    uint32_t size() override
+    {
+      return static_cast<uint32_t>( mMesh.faceCount() );
+    }
+
+    //! sets the stream pointer to the first entry, if possible.
+    void rewind() override
+    {
+      mIterator = 0;
+    }
+
+  protected:
+    void readNextEntry()
+    {
+      SpatialIndex::Region r;
+      const int faceCount = mMesh.faceCount();
+      if ( mIterator < faceCount )
+      {
+        r = faceToRegion( mMesh, mIterator );
+        mNextData = new RTree::Data(
+          0,
+          nullptr,
+          r,
+          mIterator );
+        ++mIterator;
+      }
+    }
+
+  private:
+    int mIterator = 0;
+    const QgsMesh &mMesh;
+    RTree::Data *mNextData = nullptr;
+    QgsFeedback *mFeedback = nullptr;
+};
+
+
+/**
+ * \ingroup core
+ *  \class QgsSpatialIndexData
+ * \brief Data of spatial index that may be implicitly shared
+ * \note not available in Python bindings
+*/
+class QgsMeshSpatialIndexData : public QSharedData
+{
+  public:
+    QgsMeshSpatialIndexData()
+    {
+      initTree();
+    }
+
+    /**
+     * Constructor for QgsSpatialIndexData which bulk loads features from the specified feature iterator
+     * \a fi.
+     *
+     * The optional \a feedback object can be used to allow cancelation of bulk feature loading. Ownership
+     * of \a feedback is not transferred, and callers must take care that the lifetime of feedback exceeds
+     * that of the spatial index construction.
+     */
+    explicit QgsMeshSpatialIndexData( const QgsMesh &fi, QgsFeedback *feedback = nullptr )
+    {
+      QgsMeshFaceIteratorDataStream fids( fi, feedback );
+      initTree( &fids );
+    }
+
+    QgsMeshSpatialIndexData( const QgsMeshSpatialIndexData &other )
+      : QSharedData( other )
+    {
+      QMutexLocker locker( &other.mMutex );
+
+      initTree();
+
+      // copy R-tree data one by one (is there a faster way??)
+      double low[]  = { std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest() };
+      double high[] = { std::numeric_limits<double>::max(), std::numeric_limits<double>::max() };
+      SpatialIndex::Region query( low, high, 2 );
+      QgsMeshSpatialIndexCopyVisitor visitor( mRTree.get() );
+      other.mRTree->intersectsWithQuery( query, visitor );
+    }
+
+    ~QgsMeshSpatialIndexData() = default;
+
+    QgsMeshSpatialIndexData &operator=( const QgsMeshSpatialIndexData &rh ) = delete;
+
+    void initTree( IDataStream *inputStream = nullptr )
+    {
+      // for now only memory manager
+      mStorage.reset( StorageManager::createNewMemoryStorageManager() );
+
+      // R-Tree parameters
+      double fillFactor = 0.7;
+      unsigned int indexCapacity = 10;
+      unsigned int leafCapacity = 10;
+      unsigned int dimension = 2;
+      RTree::RTreeVariant variant = RTree::RV_RSTAR;
+
+      // create R-tree
+      SpatialIndex::id_type indexId;
+
+      if ( inputStream && inputStream->hasNext() )
+        mRTree.reset(
+          RTree::createAndBulkLoadNewRTree(
+            RTree::BLM_STR,
+            *inputStream,
+            *mStorage, fillFactor,
+            indexCapacity,
+            leafCapacity,
+            dimension,
+            variant,
+            indexId )
+        );
+      else
+        mRTree.reset(
+          RTree::createNewRTree(
+            *mStorage,
+            fillFactor,
+            indexCapacity,
+            leafCapacity,
+            dimension,
+            variant,
+            indexId )
+        );
+    }
+
+    //! Storage manager
+    std::unique_ptr<SpatialIndex::IStorageManager> mStorage;
+
+    //! R-tree containing spatial index
+    std::unique_ptr<SpatialIndex::ISpatialIndex> mRTree;
+
+    mutable QMutex mMutex;
+
+};
+
+// -------------------------------------------------------------------------
+
+QgsMeshSpatialIndex::QgsMeshSpatialIndex()
+{
+  d = new QgsMeshSpatialIndexData;
+}
+
+QgsMeshSpatialIndex::QgsMeshSpatialIndex( const QgsMesh &triangularMesh, QgsFeedback *feedback )
+{
+  d = new QgsMeshSpatialIndexData( triangularMesh, feedback );
+}
+
+QgsMeshSpatialIndex::QgsMeshSpatialIndex( const QgsMeshSpatialIndex &other ) //NOLINT
+  : d( other.d )
+{
+}
+
+QgsMeshSpatialIndex:: ~QgsMeshSpatialIndex() = default; //NOLINT
+
+QgsMeshSpatialIndex &QgsMeshSpatialIndex::operator=( const QgsMeshSpatialIndex &other )
+{
+  if ( this != &other )
+    d = other.d;
+  return *this;
+}
+
+QList<int> QgsMeshSpatialIndex::intersects( const QgsRectangle &rect ) const
+{
+  QList<int> list;
+  QgisMeshVisitor visitor( list );
+
+  SpatialIndex::Region r = rectToRegion( rect );
+
+  QMutexLocker locker( &d->mMutex );
+  d->mRTree->intersectsWithQuery( r, visitor );
+
+  return list;
+}
+
+QList<int> QgsMeshSpatialIndex::nearestNeighbor( const QgsPointXY &point, int neighbors ) const
+{
+  QList<int> list;
+  QgisMeshVisitor visitor( list );
+
+  double pt[2] = { point.x(), point.y() };
+  Point p( pt, 2 );
+
+  QMutexLocker locker( &d->mMutex );
+  d->mRTree->nearestNeighborQuery( static_cast<uint32_t>( neighbors ), p, visitor );
+
+  return list;
+}
+
+///@endcond