boost::geometry can't recognize that three points are in line (boost::geometry::difference fails)

Question

I have a complex polyline processing using boost::geometry where I need to do boolean operations with them. I got a strange situation, when boost::geometry can't see that 3 points are all in line and fails to subtract a linestring made with 3 points that are colinear (according to boost::geometry itself) from a linestring made with 2 same endpoints (without the middle one).

Shouldn't the result of subtraction of 2 linestrings be empty if they lay one over another?

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <iostream>

namespace bg = boost::geometry;

using Number = double;
typedef bg::model::d2::point_xy<Number> point_type;

typedef bg::model::linestring<point_type> linestring_type;
typedef bg::model::multi_linestring<linestring_type> multilinestring_type;
typedef bg::model::segment<point_type> segment_type;
    
int main() {
    std::cout << std::setprecision(17);
    
    point_type p1{ 41.746999534390177, 58.355029632348561 };
    point_type pc{ 41.803915890274112, 58.454652240833830 };
    point_type p2{ 41.856075653483821, 58.54593925181792 };

    linestring_type ls1{ p1, p2 };
    linestring_type ls2{ p1, pc, p2 }; // same as ls1 endpoints with a center point which is slightly off

    auto d1 = bg::distance(pc, ls1);
    std::cout << "Distance between Point pc " << bg::wkt(pc) << "and line ls1 " << bg::wkt(ls1) << " is " << d1 << std::endl;

    // now project pc onto ls1
    bg::model::segment<point_type> sout;
    bg::closest_points(pc, ls1, sout);
    point_type pc_proj = sout.second;

    // check if pc_proj is on ls1 (expect zero distance since it is the result of the projection)
    auto d2 = bg::distance(pc_proj, ls1); // should be 0
    std::cout << "Distance between Point pc_proj " << bg::wkt(pc_proj) << "and line ls1 " << bg::wkt(ls1) << " is " << d2 << std::endl;

    //now, create another linestring where all three points should be in line and aligned with ls1
    linestring_type ls2_proj{ p1, pc_proj, p2 }; // same as ls2 with the aligned with ls1 center point 

    // the idea is that since all the points of ls2_proj lie on ls1, the difference should be empty
    multilinestring_type output1;
    boost::geometry::difference(ls1, ls2_proj, output1);

    // oops, non-empty difference
    std::cout << "Difference ls1 - ls2_proj: " << std::endl;
    for (auto& p : output1)
        std::cout << bg::wkt(p) << "\n";

    // and if we ask whether pc_proj is covered by ls1, it is not. But why?
    if (!bg::covered_by(pc_proj, ls1))
    {
        std::cout << "Point " << bg::wkt(pc_proj) << " is not on ls1, but distance is: " << d2 << std::endl;
    }
}

Output is

Distance between Point pc POINT(41.803915890274112 58.45465224083383)and line ls1 LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792) is 2.5817835214104254e-06
Distance between Point pc_proj POINT(41.803918131966284 58.454650960044091)and line ls1 LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792) is 0
Difference ls1 - ls2_proj: 
LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792)
Point POINT(41.803918131966284 58.454650960044091) is not on ls1, but distance is: 0

score 2 · Answer 1 · answered Dec 09 '22 at 03:25

The documentation states that behaviour is defined for:

Now, playing around with the code a bit led me to suspect floating point inaccuracies: Live On Coliru

This made me think of alternative appraoches. The simplest appears to be simplify which appears to work well, at least for a number of coordinate types:

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/linestring.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <iostream>

template <typename Number, intmax_t Scale = 1> void do_test() {
    namespace bg = boost::geometry;

    using point_type      = bg::model::d2::point_xy<Number>;
    using linestring_type = bg::model::linestring<point_type>;

    point_type const p1(41.746999534390177 * Scale, 58.355029632348561 * Scale);
    point_type const pc(41.803915890274112 * Scale, 58.454652240833830 * Scale);
    point_type const p2(41.856075653483821 * Scale, 58.54593925181792 * Scale);

    linestring_type const ls1{p1, p2};
    linestring_type const ls2{p1, pc, p2}; // center point slightly off

    linestring_type out;
    bg::simplify(ls2, out, 1e-4 * Scale);

    std::cout << " --- simplification equals ls1? " << bg::equals(out, ls1) << "\n";
    std::cout << "          ls1: " << bg::wkt(ls1) << "\n";
    std::cout << "          out: " << bg::wkt(out) << "\n";
}

int main() {
    std::cout << std::setprecision(17) << std::boolalpha;

    do_test<intmax_t, 100'000'000>();
    do_test<long double>();
    do_test<double>();
}

Prints

 --- simplification equals ls1? true
          ls1: LINESTRING(4174699953 5835502963,4185607565 5854593925)
          out: LINESTRING(4174699953 5835502963,4185607565 5854593925)
 --- simplification equals ls1? true
          ls1: LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792)
          out: LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792)
 --- simplification equals ls1? true
          ls1: LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792)
          out: LINESTRING(41.746999534390177 58.355029632348561,41.856075653483821 58.54593925181792)

Another alternative might be to allow for a tolerance "zone" by applying buffer before querying the covered_by relation. I'm not sure which of the approaches I expect to be more optimal, but I suspect it will be the simplify approach.

Note, per the docs, simplify might lead to self-intersections. You might find the check helper from my first listing handy:

auto check = [](auto& geo) {
    if (std::string reason; !bg::is_valid(geo, reason)) {
        std::cout << "Trying to correct " << bg::wkt(geo) << " reason: " << reason << "\n";
        bg::correct(geo);
        if (!bg::is_valid(geo, reason))
            std::cout << " -- failed: " << reason << "\n";
        else
            std::cout << " -- result: " << bg::wkt(geo) << "\n";
    }
};

boost::geometry can't recognize that three points are in line (boost::geometry::difference fails)

1 Answers1