[FEATURE] QgsSpatialIndexKDBush

A very fast static spatial index for 2D points based on a flat KD-tree,
using https://github.com/mourner/kdbush.hpp

Compared to QgsSpatialIndex, this index:
 - supports single point features only (no multipoints)
 - is static (features cannot be added or removed from the index after construction)
 - is much faster!
 - supports true "distance based" searches, i.e. return all points within a radius
from a search point

external/kdbush/LICENSE

@@ -0,0 +1,13 @@
+Copyright (c) 2016, Vladimir Agafonkin
+
+Permission to use, copy, modify, and/or distribute this software for any
+purpose with or without fee is hereby granted, provided that the above
+copyright notice and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

external/kdbush/README.md

@@ -0,0 +1,3 @@
+A C++11 port of [kdbush](https://github.com/mourner/kdbush), a fast static spatial index for 2D points.
+
+[![Build Status](https://travis-ci.org/mourner/kdbush.hpp.svg?branch=master)](https://travis-ci.org/mourner/kdbush.hpp)

external/kdbush/bench.cpp

@@ -0,0 +1,33 @@
+#include "include/kdbush.hpp"
+
+#include <chrono>
+#include <iostream>
+#include <random>
+#include <vector>
+
+int main() {
+    std::mt19937 gen(0);
+    std::uniform_int_distribution<> dis(-10000, 10000);
+
+    using Point = std::pair<int, int>;
+
+    const std::size_t num_points = 1000000;
+
+    std::vector<Point> points;
+    points.reserve(num_points);
+
+    for (std::size_t i = 0; i < num_points; i++) {
+        points.emplace_back(dis(gen), dis(gen));
+    }
+
+    const auto started = std::chrono::high_resolution_clock::now();
+    kdbush::KDBush<Point> index(points);
+    const auto finished = std::chrono::high_resolution_clock::now();
+
+    const auto duration =
+        std::chrono::duration_cast<std::chrono::nanoseconds>(finished - started).count();
+
+    std::cerr << "indexed 1M points in " << (duration / 1e6) << "ms\n";
+
+    return 0;
+}

external/kdbush/include/kdbush.hpp

@@ -0,0 +1,214 @@
+#pragma once
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <tuple>
+#include <vector>
+#include <cassert>
+
+namespace kdbush {
+
+template <std::uint8_t I, typename T>
+struct nth {
+    inline static typename std::tuple_element<I, T>::type get(const T &t) {
+        return std::get<I>(t);
+    }
+};
+
+template <typename TPoint, typename TIndex = std::size_t>
+class KDBush {
+
+public:
+    using TNumber = decltype(nth<0, TPoint>::get(std::declval<TPoint>()));
+    static_assert(
+        std::is_same<TNumber, decltype(nth<1, TPoint>::get(std::declval<TPoint>()))>::value,
+        "point component types must be identical");
+
+    static const std::uint8_t defaultNodeSize = 64;
+
+    KDBush(const std::uint8_t nodeSize_ = defaultNodeSize) : nodeSize(nodeSize_) {
+    }
+
+    KDBush(const std::vector<TPoint> &points_, const std::uint8_t nodeSize_ = defaultNodeSize)
+        : KDBush(std::begin(points_), std::end(points_), nodeSize_) {
+    }
+
+    template <typename TPointIter>
+    KDBush(const TPointIter &points_begin,
+           const TPointIter &points_end,
+           const std::uint8_t nodeSize_ = defaultNodeSize)
+        : nodeSize(nodeSize_) {
+        fill(points_begin, points_end);
+    }
+
+    void fill(const std::vector<TPoint> &points_) {
+        fill(std::begin(points_), std::end(points_));
+    }
+
+    template <typename TPointIter>
+    void fill(const TPointIter &points_begin, const TPointIter &points_end) {
+        assert(points.empty());
+        const TIndex size = static_cast<TIndex>(std::distance(points_begin, points_end));
+
+        if (size == 0) return;
+
+        points.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++);
+        }
+
+        sortKD(0, size - 1, 0);
+    }
+
+    template <typename TVisitor>
+    void range(const TNumber minX,
+               const TNumber minY,
+               const TNumber maxX,
+               const TNumber maxY,
+               const TVisitor &visitor) {
+        range(minX, minY, maxX, maxY, visitor, 0, static_cast<TIndex>(ids.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);
+    }
+
+protected:
+    std::vector<TIndex> ids;
+    std::vector<std::pair<TNumber, TNumber>> points;
+    std::uint8_t nodeSize;
+
+    template <typename TVisitor>
+    void range(const TNumber minX,
+               const TNumber minY,
+               const TNumber maxX,
+               const TNumber maxY,
+               const TVisitor &visitor,
+               const TIndex left,
+               const TIndex right,
+               const std::uint8_t axis) {
+
+        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]);
+            }
+            return;
+        }
+
+        const TIndex m = (left + right) >> 1;
+        const TNumber x = std::get<0>(points[m]);
+        const TNumber y = std::get<1>(points[m]);
+
+        if (x >= minX && x <= maxX && y >= minY && y <= maxY) visitor(ids[m]);
+
+        if (axis == 0 ? minX <= x : minY <= y)
+            range(minX, minY, maxX, maxY, visitor, left, m - 1, (axis + 1) % 2);
+
+        if (axis == 0 ? maxX >= x : maxY >= y)
+            range(minX, minY, maxX, maxY, visitor, m + 1, right, (axis + 1) % 2);
+    }
+
+    template <typename TVisitor>
+    void within(const TNumber qx,
+                const TNumber qy,
+                const TNumber r,
+                const TVisitor &visitor,
+                const TIndex left,
+                const TIndex right,
+                const std::uint8_t axis) {
+
+        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]);
+            }
+            return;
+        }
+
+        const TIndex m = (left + right) >> 1;
+        const TNumber x = std::get<0>(points[m]);
+        const TNumber y = std::get<1>(points[m]);
+
+        if (sqDist(x, y, qx, qy) <= r2) visitor(ids[m]);
+
+        if (axis == 0 ? qx - r <= x : qy - r <= y)
+            within(qx, qy, r, visitor, left, m - 1, (axis + 1) % 2);
+
+        if (axis == 0 ? qx + r >= x : qy + r >= y)
+            within(qx, qy, r, visitor, m + 1, right, (axis + 1) % 2);
+    }
+
+    void sortKD(const TIndex left, const TIndex right, const std::uint8_t axis) {
+        if (right - left <= nodeSize) return;
+        const TIndex m = (left + right) >> 1;
+        if (axis == 0) {
+            select<0>(m, left, right);
+        } else {
+            select<1>(m, left, right);
+        }
+        sortKD(left, m - 1, (axis + 1) % 2);
+        sortKD(m + 1, right, (axis + 1) % 2);
+    }
+
+    template <std::uint8_t I>
+    void select(const TIndex k, TIndex left, TIndex right) {
+
+        while (right > left) {
+            if (right - left > 600) {
+                const double n = right - left + 1;
+                const double m = k - left + 1;
+                const double z = std::log(n);
+                const double s = 0.5 * std::exp(2 * z / 3);
+                const double r =
+                    k - m * s / n + 0.5 * std::sqrt(z * s * (1 - s / n)) * (2 * m < n ? -1 : 1);
+                select<I>(k, std::max(left, TIndex(r)), std::min(right, TIndex(r + s)));
+            }
+
+            const TNumber t = std::get<I>(points[k]);
+            TIndex i = left;
+            TIndex j = right;
+
+            swapItem(left, k);
+            if (std::get<I>(points[right]) > 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--;
+            }
+
+            if (std::get<I>(points[left]) == t)
+                swapItem(left, j);
+            else {
+                swapItem(++j, right);
+            }
+
+            if (j <= k) left = j + 1;
+            if (k <= j) right = j - 1;
+        }
+    }
+
+    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(points.begin() + static_cast<std::int32_t>(i), points.begin() + static_cast<std::int32_t>(j));
+    }
+
+    TNumber sqDist(const TNumber ax, const TNumber ay, const TNumber bx, const TNumber by) {
+        auto dx = ax - bx;
+        auto dy = ay - by;
+        return dx * dx + dy * dy;
+    }
+};
+
+} // namespace kdbush

external/kdbush/test.cpp

@@ -0,0 +1,55 @@
+#include "include/kdbush.hpp"
+
+#include <cassert>
+#include <iostream>
+#include <iterator>
+#include <vector>
+
+using TPoint = std::pair<int, int>;
+using TIds = std::vector<std::size_t>;
+
+static std::vector<TPoint> points = {
+    { 54, 1 },  { 97, 21 }, { 65, 35 }, { 33, 54 }, { 95, 39 }, { 54, 3 },  { 53, 54 }, { 84, 72 },
+    { 33, 34 }, { 43, 15 }, { 52, 83 }, { 81, 23 }, { 1, 61 },  { 38, 74 }, { 11, 91 }, { 24, 56 },
+    { 90, 31 }, { 25, 57 }, { 46, 61 }, { 29, 69 }, { 49, 60 }, { 4, 98 },  { 71, 15 }, { 60, 25 },
+    { 38, 84 }, { 52, 38 }, { 94, 51 }, { 13, 25 }, { 77, 73 }, { 88, 87 }, { 6, 27 },  { 58, 22 },
+    { 53, 28 }, { 27, 91 }, { 96, 98 }, { 93, 14 }, { 22, 93 }, { 45, 94 }, { 18, 28 }, { 35, 15 },
+    { 19, 81 }, { 20, 81 }, { 67, 53 }, { 43, 3 },  { 47, 66 }, { 48, 34 }, { 46, 12 }, { 32, 38 },
+    { 43, 12 }, { 39, 94 }, { 88, 62 }, { 66, 14 }, { 84, 30 }, { 72, 81 }, { 41, 92 }, { 26, 4 },
+    { 6, 76 },  { 47, 21 }, { 57, 70 }, { 71, 82 }, { 50, 68 }, { 96, 18 }, { 40, 31 }, { 78, 53 },
+    { 71, 90 }, { 32, 14 }, { 55, 6 },  { 32, 88 }, { 62, 32 }, { 21, 67 }, { 73, 81 }, { 44, 64 },
+    { 29, 50 }, { 70, 5 },  { 6, 22 },  { 68, 3 },  { 11, 23 }, { 20, 42 }, { 21, 73 }, { 63, 86 },
+    { 9, 40 },  { 99, 2 },  { 99, 76 }, { 56, 77 }, { 83, 6 },  { 21, 72 }, { 78, 30 }, { 75, 53 },
+    { 41, 11 }, { 95, 20 }, { 30, 38 }, { 96, 82 }, { 65, 48 }, { 33, 18 }, { 87, 28 }, { 10, 10 },
+    { 40, 34 }, { 10, 20 }, { 47, 29 }, { 46, 78 }
+};
+
+static void testRange() {
+    kdbush::KDBush<TPoint> index(points, 10);
+    TIds expectedIds = { 3, 90, 77, 72, 62, 96, 47, 8, 17, 15, 69, 71, 44, 19, 18, 45, 60, 20 };
+    TIds result;
+    index.range(20, 30, 50, 70, [&result](const auto id) { result.push_back(id); });
+
+    assert(std::equal(expectedIds.begin(), expectedIds.end(), result.begin()));
+}
+
+static void testRadius() {
+    kdbush::KDBush<TPoint> index(points, 10);
+    TIds expectedIds = { 3, 96, 71, 44, 18, 45, 60, 6, 25, 92, 42, 20 };
+    TIds result;
+    index.within(50, 50, 20, [&result](const auto id) { result.push_back(id); });
+
+    assert(std::equal(expectedIds.begin(), expectedIds.end(), result.begin()));
+}
+
+static void testEmpty() {
+    auto emptyPoints = std::vector<TPoint>{};
+    kdbush::KDBush<TPoint> index(emptyPoints);
+}
+
+int main() {
+    testRange();
+    testRadius();
+    testEmpty();
+    return 0;
+}

python/core/auto_generated/qgsspatialindex.sip.in

@@ -27,6 +27,8 @@ QgsSpatialIndex objects are implicitly shared and can be inexpensively copied.
    While the underlying libspatialindex is not thread safe on some platforms, the QgsSpatialIndex
    class implements its own locks and accordingly, a single QgsSpatialIndex object can safely
    be used across multiple threads.
+
+.. seealso:: :py:class:`QgsSpatialIndexKDBush`
 %End
 
 %TypeHeaderCode

python/core/core_auto.sip

@@ -112,6 +112,7 @@
 %Include auto_generated/qgssimplifymethod.sip
 %Include auto_generated/qgssnappingutils.sip
 %Include auto_generated/qgsspatialindex.sip
+%Include auto_generated/qgsspatialindexkdbush.sip
 %Include auto_generated/qgssqlstatement.sip
 %Include auto_generated/qgsstatisticalsummary.sip
 %Include auto_generated/qgsstringstatisticalsummary.sip

scripts/astyle.sh

@@ -85,7 +85,7 @@ astyleit() {
 
 for f in "$@"; do
 	case "$f" in
-                src/plugins/grass/qtermwidget/*|external/o2/*|external/astyle/*|python/ext-libs/*|ui_*.py|*.astyle|tests/testdata/*|editors/*)
+                src/plugins/grass/qtermwidget/*|external/o2/*|external/astyle/*|external/kdbush/*|python/ext-libs/*|ui_*.py|*.astyle|tests/testdata/*|editors/*)
 			echo -ne "$f skipped $elcr"
 			continue
 			;;

src/core/CMakeLists.txt

@@ -2,6 +2,8 @@
 # sources
 
 SET(QGIS_CORE_SRCS
+  ${CMAKE_SOURCE_DIR}/external/kdbush/include/kdbush.hpp
+
   ${CMAKE_SOURCE_DIR}/external/nmea/context.c
   ${CMAKE_SOURCE_DIR}/external/nmea/gmath.c
   ${CMAKE_SOURCE_DIR}/external/nmea/info.c
@@ -288,6 +290,7 @@ SET(QGIS_CORE_SRCS
   qgssimplifymethod.cpp
   qgssnappingutils.cpp
   qgsspatialindex.cpp
+  qgsspatialindexkdbush.cpp
   qgssqlexpressioncompiler.cpp
   qgssqliteexpressioncompiler.cpp
   qgssqlstatement.cpp
@@ -910,6 +913,7 @@ SET(QGIS_CORE_HDRS
   qgssimplifymethod.h
   qgssnappingutils.h
   qgsspatialindex.h
+  qgsspatialindexkdbush.h
   qgsspatialiteutils.h
   qgssqlstatement.h
   qgssqliteutils.h
@@ -1184,6 +1188,7 @@ INCLUDE_DIRECTORIES(
   symbology
   metadata
   mesh
+  ${CMAKE_SOURCE_DIR}/external/kdbush/include
   ${CMAKE_SOURCE_DIR}/external/nmea
   ${CMAKE_SOURCE_DIR}/external/poly2tri
 )

src/core/qgsspatialindex.h

@@ -60,6 +60,8 @@ class QgsFeatureSource;
  * \note While the underlying libspatialindex is not thread safe on some platforms, the QgsSpatialIndex
  * class implements its own locks and accordingly, a single QgsSpatialIndex object can safely
  * be used across multiple threads.
+ *
+ * \see QgsSpatialIndexKDBush, which is an optimised non-mutable index for point geometries only.
  */
 class CORE_EXPORT QgsSpatialIndex
 {