boost::geometry无法识别三点共线导致difference操作失败

问题描述 #

使用boost::geometry处理复杂折线布尔运算时，遇到共线点识别异常：无法从仅含端点的linestring中减去包含共线中点的linestring。预期完全重叠的两个linestring相减结果为空，但实际输出原linestring，不符合预期。

测试代码 #

#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>
#include <algorithm>

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

    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;

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

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

    linestring_type ls2_proj{ p1, pc_proj, p2 };

    multilinestring_type output1;
    boost::geometry::difference(ls1, ls2_proj, output1);

    std::cout << "Difference ls1 - ls2_proj: " << std::endl;
    for (auto& p : output1)
        std::cout << bg::wkt(p) << "\n";

    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

核心疑问 #

投影点pc_proj与ls1的距离为0，但bg::covered_by判定其不在ls1上，且ls1与ls2_proj的difference结果不为空，与预期不符，该如何解决？

解决方案 #

1. 浮点精度容差设置 #

bg::distance返回0仅表示点在linestring所在的直线上，但bg::covered_by需要点严格落在linestring的线段范围内（包括端点），浮点计算的微小误差会导致判定失败。需显式设置合适的容差：

• 全局容差设置：在程序开头设置全局epsilon，适用于所有几何判定：

  bg::set_epsilon<Number>(1e-10); // 根据实际精度需求调整
  

• 带容差的策略调用：针对单个covered_by调用指定容差策略：

  bg::strategy::within::winding<double> within_strategy(1e-10);
  if (!bg::covered_by(pc_proj, ls1, within_strategy))
  {
      std::cout << "Point " << bg::wkt(pc_proj) << " is not on ls1, but distance is: " << d2 << std::endl;
  }
  

2. 修正投影点参数范围 #

closest_points返回的投影点可能因浮点误差，其在线段上的比例参数超出[0,1]范围，导致covered_by判定失败。手动计算投影点并限制参数范围：

// 替换原closest_points的投影计算
bg::model::segment<point_type> seg(p1, p2);
Number t = bg::projected_point_ratio(pc, seg);
t = std::clamp(t, Number(0), Number(1)); // 限制t在[0,1]之间
point_type pc_proj_clamped;
bg::linear_interpolate(seg, t, pc_proj_clamped);

使用修正后的pc_proj_clamped构建ls2_proj，即可让covered_by和difference运算得到预期结果。

3. 布尔运算的精度适配 #

设置全局epsilon后，difference运算会自动使用该容差识别共线点的拓扑关系，此时ls1 - ls2_proj的结果会为空。

内容的提问来源于stack exchange，提问作者Vladimir Shutow