Skip to content

Commit

Permalink
[optimization][geometry snapper] Get rid of custom tree, use GEOS
Browse files Browse the repository at this point in the history
  • Loading branch information
@nirvn
nirvn committed Apr 19, 2021
1 parent e255d0f commit 2436f82
Show file tree
Hide file tree
Showing 4 changed files with 96 additions and 141 deletions.
1 change: 1 addition & 0 deletions python/analysis/auto_generated/vector/qgsgeometrysnapper.sip.in
View file
Expand Up @@ -10,6 +10,7 @@




class QgsGeometrySnapper : QObject
{
%Docstring(signature="appended")
Expand Down
202 changes: 84 additions & 118 deletions src/analysis/vector/qgsgeometrysnapper.cpp
View file
Expand Up @@ -14,7 +14,6 @@
* *
***************************************************************************/

#include <QtConcurrentMap>
#include "qgsfeatureiterator.h"
#include "qgsgeometry.h"
#include "qgsvectorlayer.h"
Expand All @@ -25,6 +24,10 @@
#include "qgssurface.h"
#include "qgsmultisurface.h"
#include "qgscurve.h"
#include "qgsgeos.h"

#include <QtConcurrentMap>
#include <geos_c.h>

///@cond PRIVATE

Expand Down Expand Up @@ -184,128 +187,48 @@ class Raytracer

///////////////////////////////////////////////////////////////////////////////

QgsSnapIndex::GridRow::~GridRow()
{
const auto constMCells = mCells;
for ( const QgsSnapIndex::Cell &cell : constMCells )
{
qDeleteAll( cell );
}
}

QgsSnapIndex::Cell &QgsSnapIndex::GridRow::getCreateCell( int col )
{
if ( col < mColStartIdx )
{
mCells.reserve( mColStartIdx - col );
for ( int i = col; i < mColStartIdx; ++i )
{
mCells.prepend( Cell() );
}
mColStartIdx = col;
return mCells.front();
}
else if ( col >= mColStartIdx + mCells.size() )
{
mCells.reserve( col - ( mColStartIdx + mCells.size() ) );
for ( int i = mColStartIdx + mCells.size(); i <= col; ++i )
{
mCells.append( Cell() );
}
return mCells.back();
}
else
{
return mCells[col - mColStartIdx];
}
}

const QgsSnapIndex::Cell *QgsSnapIndex::GridRow::getCell( int col ) const
{
if ( col < mColStartIdx || col >= mColStartIdx + mCells.size() )
{
return nullptr;
}
else
{
return &mCells.at( col - mColStartIdx );
}
}

QList<QgsSnapIndex::SnapItem *> QgsSnapIndex::GridRow::getSnapItems( int colStart, int colEnd ) const
{
colStart = std::max( colStart, mColStartIdx );
colEnd = std::min( colEnd, mColStartIdx + mCells.size() - 1 );

QList<SnapItem *> items;
items.reserve( colEnd - colStart + 1 );
for ( int col = colStart; col <= colEnd; ++col )
{
items.append( mCells.at( col - mColStartIdx ) );
}
return items;
}

///////////////////////////////////////////////////////////////////////////////

QgsSnapIndex::QgsSnapIndex( const QgsPoint &origin, double cellSize )
: mOrigin( origin )
, mCellSize( cellSize )
, mRowsStartIdx( 0 )
{
mSTRTree = GEOSSTRtree_create_r( QgsGeos::getGEOSHandler(), ( size_t )10 );
}

QgsSnapIndex::~QgsSnapIndex()
{
qDeleteAll( mCoordIdxs );
}

qDeleteAll( mSnapItems );

const QgsSnapIndex::Cell *QgsSnapIndex::getCell( int col, int row ) const
{
if ( row < mRowsStartIdx || row >= mRowsStartIdx + mGridRows.size() )
{
return nullptr;
}
else
{
return mGridRows.at( row - mRowsStartIdx ).getCell( col );
}
}

QgsSnapIndex::Cell &QgsSnapIndex::getCreateCell( int col, int row )
{
if ( row < mRowsStartIdx )
{
mGridRows.reserve( mRowsStartIdx - row );
for ( int i = row; i < mRowsStartIdx; ++i )
{
mGridRows.prepend( GridRow() );
}
mRowsStartIdx = row;
return mGridRows.front().getCreateCell( col );
}
else if ( row >= mRowsStartIdx + mGridRows.size() )
GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
for ( GEOSGeometry *geom : mSTRTreeItems )
{
mGridRows.reserve( row - ( mRowsStartIdx + mGridRows.size() ) );
for ( int i = mRowsStartIdx + mGridRows.size(); i <= row; ++i )
{
mGridRows.append( GridRow() );
}
return mGridRows.back().getCreateCell( col );
}
else
{
return mGridRows[row - mRowsStartIdx].getCreateCell( col );
GEOSGeom_destroy_r( geosctxt, geom );
}
GEOSSTRtree_destroy_r( geosctxt, mSTRTree );
}

void QgsSnapIndex::addPoint( const CoordIdx *idx, bool isEndPoint )
{
QgsPoint p = idx->point();
int col = std::floor( ( p.x() - mOrigin.x() ) / mCellSize );
int row = std::floor( ( p.y() - mOrigin.y() ) / mCellSize );
getCreateCell( col, row ).append( new PointSnapItem( idx, isEndPoint ) );

GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
#if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, row, col ) );
#else
GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
GEOSCoordSeq_setX_r( geosctxt, seq, 0, row );
GEOSCoordSeq_setY_r( geosctxt, seq, 0, col );
geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
#endif

PointSnapItem *item = new PointSnapItem( idx, isEndPoint );
GEOSSTRtree_insert_r( geosctxt, mSTRTree, point.get(), item );
#if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR<9
mSTRTreeItems << point.release();
#endif
mSnapItems << item;
}

void QgsSnapIndex::addSegment( const CoordIdx *idxFrom, const CoordIdx *idxTo )
Expand All @@ -319,9 +242,24 @@ void QgsSnapIndex::addSegment( const CoordIdx *idxFrom, const CoordIdx *idxTo )
float y1 = ( pTo.y() - mOrigin.y() ) / mCellSize;

Raytracer rt( x0, y0, x1, y1 );
GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();
for ( ; rt.isValid(); rt.next() )
{
getCreateCell( rt.curCol(), rt.curRow() ).append( new SegmentSnapItem( idxFrom, idxTo ) );
#if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
geos::unique_ptr point( GEOSGeom_createPointFromXY_r( geosctxt, rt.curRow(), rt.curCol() ) );
#else
GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
GEOSCoordSeq_setX_r( geosctxt, seq, 0, rt.curRow() );
GEOSCoordSeq_setY_r( geosctxt, seq, 0, rt.curCol() );
geos::unique_ptr point( GEOSGeom_createPoint_r( geosctxt, seq ) );
#endif

SegmentSnapItem *item = new SegmentSnapItem( idxFrom, idxTo );
GEOSSTRtree_insert_r( geosctxt, mSTRTree, point.get(), item );
#if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR<9
mSTRTreeItems << point.release();
#endif
mSnapItems << item;
}
}

Expand Down Expand Up @@ -355,9 +293,20 @@ void QgsSnapIndex::addGeometry( const QgsAbstractGeometry *geom )
}
}

struct _GEOSQueryCallbackData
{
QList< QgsSnapIndex::SnapItem * > *list;
};

void _GEOSQueryCallback( void *item, void *userdata )
{
reinterpret_cast<_GEOSQueryCallbackData *>( userdata )->list->append( static_cast<QgsSnapIndex::SnapItem *>( item ) );
}

QgsPoint QgsSnapIndex::getClosestSnapToPoint( const QgsPoint &p, const QgsPoint &q )
{
GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();

// Look for intersections on segment from the target point to the point opposite to the point reference point
// p2 = p1 + 2 * (q - p1)
QgsPoint p2( 2 * q.x() - p.x(), 2 * q.y() - p.y() );
Expand All @@ -373,12 +322,19 @@ QgsPoint QgsSnapIndex::getClosestSnapToPoint( const QgsPoint &p, const QgsPoint
QgsPoint pMin = p;
for ( ; rt.isValid(); rt.next() )
{
const Cell *cell = getCell( rt.curCol(), rt.curRow() );
if ( !cell )
{
continue;
}
for ( const SnapItem *item : *cell )
#if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
geos::unique_ptr searchPoint( GEOSGeom_createPointFromXY_r( geosctxt, rt.curRow(), rt.curCol() ) );
#else
GEOSCoordSequence *seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
GEOSCoordSeq_setX_r( geosctxt, seq, 0, rt.curRow() );
GEOSCoordSeq_setY_r( geosctxt, seq, 0, rt.curCol() );
geos::unique_ptr searchPoint( GEOSGeom_createPoint_r( geosctxt, seq ) );
#endif
QList<SnapItem *> items;
struct _GEOSQueryCallbackData callbackData;
callbackData.list = &items;
GEOSSTRtree_query_r( geosctxt, mSTRTree, searchPoint.get(), _GEOSQueryCallback, &callbackData );
for ( const SnapItem *item : items )
{
if ( item->type == SnapSegment )
{
Expand Down Expand Up @@ -406,15 +362,25 @@ QgsSnapIndex::SnapItem *QgsSnapIndex::getSnapItem( const QgsPoint &pos, double t
int colEnd = std::floor( ( pos.x() + tol - mOrigin.x() ) / mCellSize );
int rowEnd = std::floor( ( pos.y() + tol - mOrigin.y() ) / mCellSize );

rowStart = std::max( rowStart, mRowsStartIdx );
rowEnd = std::min( rowEnd, mRowsStartIdx + mGridRows.size() - 1 );
GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler();

GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 2, 2 );
#if GEOS_VERSION_MAJOR>3 || GEOS_VERSION_MINOR>=8
GEOSCoordSeq_setXY_r( geosctxt, coord, 0, rowStart, colStart );
GEOSCoordSeq_setXY_r( geosctxt, coord, 1, rowEnd, colEnd );
#else
GEOSCoordSeq_setX_r( geosctxt, coord, 0, rowStart );
GEOSCoordSeq_setY_r( geosctxt, coord, 0, colStart );
GEOSCoordSeq_setX_r( geosctxt, coord, 1, rowEnd );
GEOSCoordSeq_setY_r( geosctxt, coord, 1, colEnd );
#endif

geos::unique_ptr searchDiagonal( GEOSGeom_createLineString_r( geosctxt, coord ) );

QList<SnapItem *> items;
items.reserve( rowEnd - rowStart + 1 );
for ( int row = rowStart; row <= rowEnd; ++row )
{
items.append( mGridRows.at( row - mRowsStartIdx ).getSnapItems( colStart, colEnd ) );
}
struct _GEOSQueryCallbackData callbackData;
callbackData.list = &items;
GEOSSTRtree_query_r( geosctxt, mSTRTree, searchDiagonal.get(), _GEOSQueryCallback, &callbackData );

double minDistSegment = std::numeric_limits<double>::max();
double minDistPoint = std::numeric_limits<double>::max();
Expand Down
29 changes: 8 additions & 21 deletions src/analysis/vector/qgsgeometrysnapper.h
View file
Expand Up @@ -17,15 +17,17 @@
#ifndef QGS_GEOMETRY_SNAPPER_H
#define QGS_GEOMETRY_SNAPPER_H

#include <QMutex>
#include <QFuture>
#include <QStringList>
#include "qgsspatialindex.h"
#include "qgsabstractgeometry.h"
#include "qgspoint.h"
#include "qgsgeometry.h"
#include "qgis_analysis.h"

#include <QMutex>
#include <QFuture>
#include <QStringList>
#include <geos_c.h>

class QgsVectorLayer;

/**
Expand Down Expand Up @@ -226,32 +228,17 @@ class QgsSnapIndex
typedef QList<SnapItem *> Cell;
typedef QPair<QgsPoint, QgsPoint> Segment;

class GridRow
{
public:
GridRow() = default;
~GridRow();
const Cell *getCell( int col ) const;
Cell &getCreateCell( int col );
QList<SnapItem *> getSnapItems( int colStart, int colEnd ) const;

private:
QList<QgsSnapIndex::Cell> mCells;
int mColStartIdx = 0;
};

QgsPoint mOrigin;
double mCellSize;

QList<CoordIdx *> mCoordIdxs;
QList<GridRow> mGridRows;
int mRowsStartIdx;
QList<SnapItem *> mSnapItems;

void addPoint( const CoordIdx *idx, bool isEndPoint );
void addSegment( const CoordIdx *idxFrom, const CoordIdx *idxTo );
const Cell *getCell( int col, int row ) const;
Cell &getCreateCell( int col, int row );

GEOSSTRtree *mSTRTree = nullptr;
QList<GEOSGeometry *> mSTRTreeItems;
};

///@endcond
Expand Down
5 changes: 3 additions & 2 deletions tests/src/analysis/testqgsgeometrysnapper.cpp
View file
Expand Up @@ -20,6 +20,8 @@ Email : nyall dot dawson at gmail dot com
#include <qgsapplication.h>
#include "qgsvectordataprovider.h"
#include "qgsvectorlayer.h"
#include "qgsgeos.h"
#include <geos_c.h>


class TestQgsGeometrySnapper : public QObject
Expand Down Expand Up @@ -51,6 +53,7 @@ class TestQgsGeometrySnapper : public QObject
void insertExtra();
void duplicateNodes();
void snapMultiPolygonToPolygon();

};

void TestQgsGeometrySnapper::initTestCase()
Expand All @@ -68,11 +71,9 @@ void TestQgsGeometrySnapper::cleanupTestCase()
}
void TestQgsGeometrySnapper::init()
{

}
void TestQgsGeometrySnapper::cleanup()
{

}

void TestQgsGeometrySnapper::snapPolygonToPolygon()
Expand Down

0 comments on commit 2436f82

Please sign in to comment.