Skip to content
Permalink
Browse files

improve speed of creation spatial index for meshes

  • Loading branch information
PeterPetrik committed Jan 24, 2019
1 parent a83a1b2 commit 9d4d1c99d4de5a9ac5be4287ad3679e3743ff17e
@@ -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
@@ -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

0 comments on commit 9d4d1c9

Please sign in to comment.
You can’t perform that action at this time.