c++boostgeometryr-tree

C++ Boost RTree overlaps does not work correctly

I am kinda confused and think I get something wrong but I really fail to see it.

I have a boost::geometry::index::rtree that stores 2D boxes of 2D geographic points. Now I try to check if a new box I add to that rtree overlaps with any box already in that rtree. And this check somehow fails for one test and I really don't get why because I do not believe the error is in the rtree/overlaps implementation.

My code is the following (in Visual Studio test environment):

using CoordinateSystem = boost::geometry::cs::geographic<boost::geometry::degree>;
using Point = boost::geometry::model::point<double, 2, CoordinateSystem>;

using RTreeBox = boost::geometry::model::box<Point>;
using RTreeEntry = std::pair<RTreeBox, uint64_t>;

constexpr static auto kRTreeMaxElementsPerNode = 4;
using RTreeAlgorithm = boost::geometry::index::rstar<kRTreeMaxElementsPerNode>;

using RTree = boost::geometry::index::rtree<RTreeEntry, RTreeAlgorithm>;

bool TestAddTreeEntry(RTree& tree, uint64_t index, RTreeBox box)
{
    if (!boost::geometry::is_valid(box)) {
        boost::geometry::correct(box);
    }

    std::vector<RTreeEntry> query_results;
    tree.query(boost::geometry::index::overlaps(box), std::back_inserter(query_results));
    if (query_results.size() > 0)
    {
        return false;
    }

    tree.insert(std::make_pair(box, index));
    return true;
}

TEST_METHOD(test_rtree_mapping) {

    RTree tree;

    Assert::IsTrue(TestAddTreeEntry(tree, 1, RTreeBox({ 1, 1 }, { 3, 3 })));
    Assert::IsTrue(TestAddTreeEntry(tree, 1, RTreeBox({ 4, 1 }, { 9, 5 })));
    Assert::IsTrue(TestAddTreeEntry(tree, 1, RTreeBox({ 1, 4 }, { 2, 9 })));

    Assert::IsFalse(TestAddTreeEntry(tree, 1, RTreeBox({ 1, 2.75 }, { 2, 9 })));
    Assert::IsFalse(TestAddTreeEntry(tree, 1, RTreeBox({ 1, 4 }, { 3.5, 9 })));
}

The first Assert::IsFalse works - but also unexpectedly only overlaps with the first box ({1, 1}, {3, 3}) and not with the third one ({1, 4}, {2, 9}). The second Assert::IsFalse does not work because the entry is successfully added.

Anyone knows a reason behind this? Has this something to do with geographic coordinates that I do not understand yet?

Kind regards

Solution

  • Overlaps doesn't do what you expect. It's too strict. The docs state:

    The function overlaps implements function Overlaps from the OGC Simple Feature Specification.

    The OGC document contains:

    enter image description here

    That's admittedly hard to parse, but that's the nature of technical specifications. Luckily, intersects is much simpler:

    a.Intersects(b) ⇔ ! a.Disjoint(b)

    DEMO

    I created a test program that

    Live On Compiler Explorer

    #include <boost/geometry.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <fstream>

namespace bg  = boost::geometry;
namespace bgm = bg::model;
namespace bgi = bg::index;

using CS    = bg::cs::geographic<bg::degree>;
using Point = bg::model::point<double, 2, CS>;
using Box   = bg::model::box<Point>;
using Tree  = std::pair<Box, uint64_t>;

constexpr static auto kRTreeMaxElementsPerNode = 4;
using RTreeAlgorithm = bgi::rstar<kRTreeMaxElementsPerNode>;
using RTree          = bgi::rtree<Tree, RTreeAlgorithm>;

#ifndef RELATION
#define RELATION bgi::intersects
#endif
#define RELATION_STR BOOST_PP_STRINGIZE(RELATION)

bool add(RTree& tree, uint64_t index, Box box) {
    if (std::string reason; !bg::is_valid(box, reason)) {
        std::cerr << "Trying to correct: " << bg::dsv(box) << " (" << reason << ")" << std::endl;
        bg::correct(box);

        assert(bg::is_valid(box));
    }

    for (auto it = qbegin(tree, RELATION(box)), e = qend(tree); it != e; ++it)
        return false;

    tree.insert(std::make_pair(box, index));

    return true;
}

int main() {
    std::cout << std::boolalpha;
    RTree tree;

    struct {
        Box         box;
        char const* name;
        char const* color;
    } const shapes[]{
        {{{1.00, 1.00}, {3.00, 3.00}}, "box1", "red"},
        {{{4.00, 1.00}, {9.00, 5.00}}, "box2", "green"},
        {{{1.00, 4.00}, {2.00, 9.00}}, "box3", "blue"},
        {{{1.00, 2.75}, {2.00, 9.00}}, "probe1", "orange"},
        {{{1.00, 4.00}, {3.50, 6.00}}, "probe2", "gray"},
    };

    for (auto const& s : shapes) {
        auto idx = (&s - shapes);
        auto added = add(tree, idx, s.box);
        std::cout << "Adding " << s.name << " as #" << idx << ": "
                  << (added ? "ACCEPT" : "REJECTED") << " " << bg::dsv(s.box)
                  << "\n";

        if (!added) {
            for (auto it = qbegin(tree, RELATION(s.box)), e = qend(tree); it != e; ++it) {
                std::cout << " - because " << s.name << " " << RELATION_STR
                          << " " << shapes[it->second].name << "\n";
            }
        }
    }

    {
        std::ofstream ofs("output.svg");
        bg::svg_mapper<Point> svg(ofs, 400, 400);

        auto style = [](std::string c) {
            return "fill-rule:nonzero;fill-opacity:0.25;fill:" + c +
                ";stroke:" + c + ";stroke-width:1;";
        };

        for (auto const& [b, name, color] : shapes) {
            svg.add(b);
        }

        for (auto const& [b, name, color] : shapes) {
            auto h  = b.max_corner().get<1>() - b.min_corner().get<1>();
            auto mc = b.max_corner();
            mc.set<1>(mc.get<1>() - h / 2);
            svg.text(mc,
                     name +
                         ("\n" + boost::lexical_cast<std::string>(bg::dsv(b))),
                     style(color));
        }

        for (auto const& [b, _, color] : shapes)
            svg.map(b, style(color));
    }
}

    With either -DRELATION=bgi::intersects or -DRELATION=!bgi::disjoint:

    Adding box1 as #0: ACCEPT ((1, 1), (3, 3))
Adding box2 as #1: ACCEPT ((4, 1), (9, 5))
Adding box3 as #2: ACCEPT ((1, 4), (2, 9))
Adding probe1 as #3: REJECTED ((1, 2.75), (2, 9))
 - because probe1 bgi::intersects box1
 - because probe1 bgi::intersects box3
Adding probe2 as #4: REJECTED ((1, 4), (3.5, 6))
 - because probe2 bgi::intersects box3

    (replacing bgi::intersects with !bgi::disjoint in the output obviously).

    Visualization:

    enter image description here

    Note how the output changes if you change probe1:

        {{{1.00, 2.75}, {1.50, 9.00}}, "probe1", "orange"},

    That looks like:

    enter image description here

    This would behave as you expected even with the bgi::overlaps. Overlaps is just too strict for the expectated output.