Skip to content
Permalink
Browse files

Modify kdbush to store feature id alonside point, for optimised stora…

…ge/retrieval
  • Loading branch information
nyalldawson committed Jul 7, 2018
1 parent c41f207 commit 66c17880bef7692611f662862a3a6b7a157fb101
@@ -16,7 +16,7 @@ struct nth {
}
};

template <typename TPoint, typename TIndex = std::size_t>
template <typename TPoint, typename TContainer, typename TIndex = std::size_t>
class KDBush {

public:
@@ -54,12 +54,12 @@ class KDBush {
if (size == 0) return;

points.reserve(size);
ids.reserve(size);
//ids.reserve(size);

TIndex i = 0;
for (auto p = points_begin; p != points_end; p++) {
points.emplace_back(nth<0, TPoint>::get(*p), nth<1, TPoint>::get(*p));
ids.push_back(i++);
//ids.push_back(i++);
}

sortKD(0, size - 1, 0);
@@ -71,17 +71,17 @@ class KDBush {
const TNumber maxX,
const TNumber maxY,
const TVisitor &visitor) {
range(minX, minY, maxX, maxY, visitor, 0, static_cast<TIndex>(ids.size() - 1), 0);
range(minX, minY, maxX, maxY, visitor, 0, static_cast<TIndex>(points.size() - 1), 0);
}

template <typename TVisitor>
void within(const TNumber qx, const TNumber qy, const TNumber r, const TVisitor &visitor) {
within(qx, qy, r, visitor, 0, static_cast<TIndex>(ids.size() - 1), 0);
within(qx, qy, r, visitor, 0, static_cast<TIndex>(points.size() - 1), 0);
}

protected:
std::vector<TIndex> ids;
std::vector<std::pair<TNumber, TNumber>> points;
//std::vector<TIndex> ids;
std::vector<TContainer> points;
std::uint8_t nodeSize;

template <typename TVisitor>
@@ -93,21 +93,23 @@ class KDBush {
const TIndex left,
const TIndex right,
const std::uint8_t axis) {
if ( points.empty() )
return;

if (right - left <= nodeSize) {
for (auto i = left; i <= right; i++) {
const TNumber x = std::get<0>(points[i]);
const TNumber y = std::get<1>(points[i]);
if (x >= minX && x <= maxX && y >= minY && y <= maxY) visitor(ids[i]);
const TNumber x = std::get<0>(points[i].coords);
const TNumber y = std::get<1>(points[i].coords);
if (x >= minX && x <= maxX && y >= minY && y <= maxY) visitor(points[i]);
}
return;
}

const TIndex m = (left + right) >> 1;
const TNumber x = std::get<0>(points[m]);
const TNumber y = std::get<1>(points[m]);
const TNumber x = std::get<0>(points[m].coords);
const TNumber y = std::get<1>(points[m].coords);

if (x >= minX && x <= maxX && y >= minY && y <= maxY) visitor(ids[m]);
if (x >= minX && x <= maxX && y >= minY && y <= maxY) visitor(points[m]);

if (axis == 0 ? minX <= x : minY <= y)
range(minX, minY, maxX, maxY, visitor, left, m - 1, (axis + 1) % 2);
@@ -125,22 +127,25 @@ class KDBush {
const TIndex right,
const std::uint8_t axis) {

if ( points.empty() )
return;

const TNumber r2 = r * r;

if (right - left <= nodeSize) {
for (auto i = left; i <= right; i++) {
const TNumber x = std::get<0>(points[i]);
const TNumber y = std::get<1>(points[i]);
if (sqDist(x, y, qx, qy) <= r2) visitor(ids[i]);
const TNumber x = std::get<0>(points[i].coords);
const TNumber y = std::get<1>(points[i].coords);
if (sqDist(x, y, qx, qy) <= r2) visitor(points[i]);
}
return;
}

const TIndex m = (left + right) >> 1;
const TNumber x = std::get<0>(points[m]);
const TNumber y = std::get<1>(points[m]);
const TNumber x = std::get<0>(points[m].coords);
const TNumber y = std::get<1>(points[m].coords);

if (sqDist(x, y, qx, qy) <= r2) visitor(ids[m]);
if (sqDist(x, y, qx, qy) <= r2) visitor(points[m]);

if (axis == 0 ? qx - r <= x : qy - r <= y)
within(qx, qy, r, visitor, left, m - 1, (axis + 1) % 2);
@@ -175,20 +180,20 @@ class KDBush {
select<I>(k, std::max(left, TIndex(r)), std::min(right, TIndex(r + s)));
}

const TNumber t = std::get<I>(points[k]);
const TNumber t = std::get<I>(points[k].coords);
TIndex i = left;
TIndex j = right;

swapItem(left, k);
if (std::get<I>(points[right]) > t) swapItem(left, right);
if (std::get<I>(points[right].coords) > t) swapItem(left, right);

while (i < j) {
swapItem(i++, j--);
while (std::get<I>(points[i]) < t) i++;
while (std::get<I>(points[j]) > t) j--;
while (std::get<I>(points[i].coords) < t) i++;
while (std::get<I>(points[j].coords) > t) j--;
}

if (std::get<I>(points[left]) == t)
if (std::get<I>(points[left].coords) == t)
swapItem(left, j);
else {
swapItem(++j, right);
@@ -200,7 +205,7 @@ class KDBush {
}

void swapItem(const TIndex i, const TIndex j) {
std::iter_swap(ids.begin() + static_cast<std::int32_t>(i), ids.begin() + static_cast<std::int32_t>(j));
// std::iter_swap(ids.begin() + static_cast<std::int32_t>(i), ids.begin() + static_cast<std::int32_t>(j));
std::iter_swap(points.begin() + static_cast<std::int32_t>(i), points.begin() + static_cast<std::int32_t>(j));
}

@@ -0,0 +1,47 @@
/************************************************************************
* This file has been generated automatically from *
* *
* src/core/qgsspatialindexkdbushdata.h *
* *
* Do not edit manually ! Edit header and run scripts/sipify.pl again *
************************************************************************/




class QgsSpatialIndexKDBushData
{
%Docstring

A container for data stored inside a QgsSpatialIndexKDBush index.

.. versionadded:: 3.4
%End

%TypeHeaderCode
#include "qgsspatialindexkdbushdata.h"
%End
public:

QgsSpatialIndexKDBushData( QgsFeatureId id, double x, double y );
%Docstring
Constructor for QgsSpatialIndexKDBushData, for a feature with the
given ``id`` and ``x``, ``y`` coordinate.
%End


QgsPointXY point() const;
%Docstring
Returns the indexed point.
%End

QgsFeatureId id;
};

/************************************************************************
* This file has been generated automatically from *
* *
* src/core/qgsspatialindexkdbushdata.h *
* *
* Do not edit manually ! Edit header and run scripts/sipify.pl again *
************************************************************************/
@@ -113,6 +113,7 @@
%Include auto_generated/qgssnappingutils.sip
%Include auto_generated/qgsspatialindex.sip
%Include auto_generated/qgsspatialindexkdbush.sip
%Include auto_generated/qgsspatialindexkdbushdata.sip
%Include auto_generated/qgssqlstatement.sip
%Include auto_generated/qgsstatisticalsummary.sip
%Include auto_generated/qgsstringstatisticalsummary.sip
@@ -915,6 +915,7 @@ SET(QGIS_CORE_HDRS
qgsspatialindex.h
qgsspatialindexkdbush.h
qgsspatialindexkdbush_p.h
qgsspatialindexkdbushdata.h
qgsspatialiteutils.h
qgssqlstatement.h
qgssqliteutils.h
@@ -16,10 +16,10 @@
***************************************************************************/

#include "qgsspatialindexkdbush.h"
#include "qgsspatialindexkdbush_p.h"
#include "qgsfeatureiterator.h"
#include "qgsfeedback.h"
#include "qgsfeaturesource.h"
#include "qgsspatialindexkdbush_p.h"

QgsSpatialIndexKDBush::QgsSpatialIndexKDBush( QgsFeatureIterator &fi, QgsFeedback *feedback )
: d( new QgsSpatialIndexKDBushPrivate( fi, feedback ) )
@@ -59,11 +59,11 @@ QgsSpatialIndexKDBush::~QgsSpatialIndexKDBush()
QSet<QgsFeatureId> QgsSpatialIndexKDBush::within( const QgsPointXY &point, double radius ) const
{
QSet<QgsFeatureId> result;
d->index->within( point.x(), point.y(), radius, [&result]( const QgsFeatureId id ) { result.insert( id ); } );
d->index->within( point.x(), point.y(), radius, [&result]( const QgsSpatialIndexKDBushData & p ) { result.insert( p.id ); } );
return result;
}

void QgsSpatialIndexKDBush::within( const QgsPointXY &point, double radius, const std::function<void( QgsFeatureId )> &visitor )
void QgsSpatialIndexKDBush::within( const QgsPointXY &point, double radius, const std::function<void( QgsSpatialIndexKDBushData )> &visitor )
{
d->index->within( point.x(), point.y(), radius, visitor );
}
@@ -84,11 +84,11 @@ QSet<QgsFeatureId> QgsSpatialIndexKDBush::intersect( const QgsRectangle &rectang
d->index->range( rectangle.xMinimum(),
rectangle.yMinimum(),
rectangle.xMaximum(),
rectangle.yMaximum(), [&result]( const QgsFeatureId id ) { result << id; } );
rectangle.yMaximum(), [&result]( const QgsSpatialIndexKDBushData & p ) { result << p.id; } );
return result;
}

void QgsSpatialIndexKDBush::intersect( const QgsRectangle &rectangle, const std::function<void ( QgsFeatureId )> &visitor ) const
void QgsSpatialIndexKDBush::intersect( const QgsRectangle &rectangle, const std::function<void( QgsSpatialIndexKDBushData )> &visitor ) const
{
d->index->range( rectangle.xMinimum(),
rectangle.yMinimum(),
@@ -27,6 +27,7 @@ class QgsRectangle;

#include "qgis_core.h"
#include "qgsfeature.h"
#include "qgsspatialindexkdbushdata.h"
#include <memory>
#include <QList>

@@ -92,7 +93,7 @@ class CORE_EXPORT QgsSpatialIndexKDBush
*
* \note Not available in Python bindings
*/
void intersect( const QgsRectangle &rectangle, const std::function<void( QgsFeatureId )> &visitor ) const SIP_SKIP;
void intersect( const QgsRectangle &rectangle, const std::function<void( QgsSpatialIndexKDBushData )> &visitor ) const SIP_SKIP;

/**
* Returns the set of features which are within the given search \a radius
@@ -106,7 +107,7 @@ class CORE_EXPORT QgsSpatialIndexKDBush
*
* \note Not available in Python bindings
*/
void within( const QgsPointXY &point, double radius, const std::function<void( QgsFeatureId )> &visitor ) SIP_SKIP;
void within( const QgsPointXY &point, double radius, const std::function<void( QgsSpatialIndexKDBushData )> &visitor ) SIP_SKIP;

/**
* Fetches the point from the index with matching \a id and stores it in \a point.
@@ -32,15 +32,16 @@
//

#include "qgsfeature.h"

#include "qgsspatialindexkdbushdata.h"
#include "qgsfeatureiterator.h"
#include "qgsfeedback.h"
#include "qgsfeaturesource.h"
#include <memory>
#include <QList>
#include "kdbush.hpp"

class PointXYKDBush : public kdbush::KDBush< std::pair<double, double>, QgsFeatureId >

class PointXYKDBush : public kdbush::KDBush< std::pair<double, double>, QgsSpatialIndexKDBushData, std::size_t >
{
public:

@@ -52,14 +53,13 @@ class PointXYKDBush : public kdbush::KDBush< std::pair<double, double>, QgsFeatu
explicit PointXYKDBush( const QgsFeatureSource &source, QgsFeedback *feedback )
{
points.reserve( source.featureCount() );
ids.reserve( source.featureCount() );
QgsFeatureIterator it = source.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ) );
fillFromIterator( it, feedback );
}

void fillFromIterator( QgsFeatureIterator &fi, QgsFeedback *feedback = nullptr )
{
QgsFeatureId size = 0;
std::size_t size = 0;

QgsFeature f;
while ( fi.nextFeature( f ) )
@@ -73,29 +73,31 @@ class PointXYKDBush : public kdbush::KDBush< std::pair<double, double>, QgsFeatu
if ( QgsWkbTypes::flatType( f.geometry().wkbType() ) == QgsWkbTypes::Point )
{
const QgsPoint *point = qgsgeometry_cast< const QgsPoint * >( f.geometry().constGet() );
points.emplace_back( point->x(), point->y() );
mIdToPoint[ f.id() ] = QgsPointXY( point->x(), point->y() );
points.emplace_back( QgsSpatialIndexKDBushData( f.id(), point->x(), point->y() ) );
}
else
{
// not a point
continue;
}

ids.push_back( f.id() );
size++;
}

if ( size == 0 )
return;

sortKD( 0, size - 1, 0 );
}

bool point( QgsFeatureId id, QgsPointXY &point ) const
{
auto it = mIdToPoint.constFind( id );
if ( it == mIdToPoint.constEnd() )
auto it = std::find_if( points.begin(), points.end(),
[id]( const QgsSpatialIndexKDBushData & d ) { return d.id == id; } );
if ( it == points.end() )
return false;

point = *it;
point = QgsPointXY( it->coords.first, it->coords.second );
return true;
}

@@ -104,10 +106,6 @@ class PointXYKDBush : public kdbush::KDBush< std::pair<double, double>, QgsFeatu
return points.size();
}

private:

QHash< QgsFeatureId, QgsPointXY > mIdToPoint;

};

class QgsSpatialIndexKDBushPrivate

0 comments on commit 66c1788

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