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

Answer 1

The documentation states that behaviour is defined for:

Case	Behavior
areal (e.g. polygon)	All combinations of: box, ring, polygon, multi_polygon
linear (e.g. linestring) / areal (e.g. polygon)	A combinations of a (multi) linestring with a (multi) polygon results in a collection of linestrings
linear (e.g. linestring)	All combinations of: linestring, multi_linestring; results in a collection of linestrings
pointlike (e.g. point)	All combinations of: point, multi_point; results in a collection of points
Other geometries	Not yet supported in this version
Spherical	Not yet supported in this version
Three dimensional	Not yet supported in this version

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";
    }
};