Search code examples
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

    • allows you to override the RELATION (e.g. bgi::intersects, bgi::overlaps, !bgi::disjoint)
    • also writes a SVG vizualtion of the shapes involved.

    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.