boost::geometry 无法识别三个点在一条线上（boost::geometry::difference 失败）

Question

我使用 boost::geometry 进行复杂的折线处理，我需要对它们进行布尔运算。我遇到了一个奇怪的情况，当 boost::geometry 无法看到 3 个点都在一条线上，并且无法从用 2 个点制成的线串中减去由 3 个共线点制成的线串（根据 boost::geometry 本身）时相同的端点（没有中间端点）。

如果 2 个线串重叠在一起，相减的结果不应该为空吗？

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

输出是

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

文档指出行为的定义是：

案例	行为
面积（例如多边形）	所有组合：长方体、圆环、多边形、多多边形
线性（例如线串）/面积（例如多边形）	（多）线串与（多）多边形的组合会产生线串集合
线性（例如线串）	所有组合：linestring、multi_linestring；结果是线串的集合
点状（例如点）	所有组合：point、multi_point；结果是点的集合
其他几何形状	此版本尚不支持
球形	此版本尚不支持
三维	此版本尚不支持

现在，稍微修改一下代码让我怀疑浮点不准确：Live On Coliru

这让我想到了其他方法。最简单的似乎是

simplify

，它看起来效果很好，至少对于许多坐标类型来说：

住在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>();
}

打印

 --- 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)

另一种选择可能是在查询 buffer

 关系之前应用 

covered_by 来允许公差“区域”。我不确定我期望哪种方法更优化，但我怀疑这将是

simplify

方法。

---

注意，根据文档，

simplify

可能会导致自相交。你 可能会从我的

first 找到 

check 帮手 列表方便：

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