c++ - 像 union_ 这样的 Boost.Geometry 操作如何处理浮点类型的基本不精确性？

Question

我正在尝试判断是否/如何使 Boost.Geometry 适用于特定用例。但是，我无法在任何地方找到有关该库如何处理浮点类型的文档。
如果你在官方文档中搜索“epsilon”这个词，就我所知，你会得到零命中；然而，从库的行为中可以清楚地看出，在进行比较时，它隐式地使用了处理浮点数的典型方式的某种版本，因为例如， union_ 操作将 union 两个彼此靠近但不重叠的多边形，如果它们是够近了。
例如，考虑以下代码，该代码执行二分搜索以确定两个单位正方形需要在合并时被视为相邻的阈值距离:

namespace bg = boost::geometry;

using point = bg::model::d2::point_xy<double>;
using polygon = bg::model::polygon<point, false>;

polygon create_poly(std::vector<std::tuple<double, double>> pts) {
    polygon poly;
    for (const auto& [x, y] : pts)
        bg::append(poly, bg::make<point>(x, y));
    auto [x_1, y_1] = pts[0];
    bg::append(poly, bg::make<point>(x_1, y_1));
    return poly;
}

bool perform_simple_union(const polygon& p1, const polygon& p2) {
    std::vector<polygon> output; 
    bg::union_(p1, p2, output);
    return output.size() == 1;
}

double find_epsilon(double left, double right) {

    if (right - left < std::numeric_limits<double>::epsilon())
        return left;
    double eps = (left + right) / 2;

    polygon a = create_poly(
        std::vector<std::tuple<double, double>>{
            {1.0, 1.0}, { 2.0,1.0 }, { 2.0, 2.0 }, { 1.0,2.0 }
        }
    );

    polygon b = create_poly(
        std::vector<std::tuple<double, double>>{
            {2.0 + eps, 1.0}, { 3.0 + eps, 1.0 }, { 3.0 + eps, 2.0 }, { 2.0 + eps,2.0 }
        }
    );

    if ( perform_simple_union(a, b) ) {
        return find_epsilon(eps, right);
    } else {
        return find_epsilon(left, eps);
    }
}

int main()
{
    auto eps = find_epsilon(0.0, 1.0);
    std::cout << "eps == " << eps << "\n";
}

当我用 Visual Studio 编译并运行上面的代码时，我得到了输出

eps == 1e-07

这是关于单精度浮点数的数字限制 epsilon。因此，如果 double 坐标在彼此之间的单精度 epsilon 内，则它是否将 double 坐标视为等效？
基本上我只想知道默认行为是什么，所以我可以决定它是否适合我。

Answer 1

在[介绍][1]中，它指出:

The library supports high precision arithmetic numbers, such as ttmath.[1]: https://www.boost.org/doc/libs/1_70_0/libs/geometry/doc/html/geometry/introduction.html

库 design rationale 对此有更多的了解:

[...], it would be too long, and it is not related to geometry. We just assume that there is a meta-function select_most_precise selecting the best type.

它们还沿 OGC 简单功能规范实现，这可能意味着您可以找到 more algorithmic robustness guarantees there 。
我从阅读代码中知道，有某些算法会考虑边缘情况，在这些情况下，结果可以变得更加健壮(通过按特定顺序执行操作或注意特征非常接近，IIRC)。一个简单的 grep 例如 robust 可能会向您展示一些道路:

policies/robustness/robust_point_type.hpp:// Meta-function to typedef a robust point type for a poli

algorithms/detail/overlay/get_turn_info_helpers.hpp: // Used ranges - owned by get_turns or (for

algorithms/detail/overlay/get_turn_info_helpers.hpp:// Version with rescaling, having robust points

algorithms/detail/overlay/append_no_dups_or_spikes.hpp: // Try using specified robust policy

我只是浅尝辄止，我并不声称了解那里的大部分内容。
使用任意精度或小数
精度是一个维度，输入为十进制形式时的源保真度是另一个维度。如果不使用 MPFR/GMP/ttmath(如上所述)，您可以轻松地使用 Boost Multiprecision。这为您提供了快速的概念验证，因为它附带了 boost，并且还允许您透明地切换到 GMP 或 MPFR 后端。
也可以看看:

Boost Geometry and exact point types

How to use Boost::Geometry _union with integers

Boost geometry intersection does not output correctly

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <iostream>
namespace mp = boost::multiprecision;
namespace bg = boost::geometry;

//// Note, cpp_dec_float<0> is variable-precision!
// using Number = mp::number<mp::cpp_dec_float<0>, mp::et_off>;

// Fixed precision, avoids allocating and populates std::numeric_limits<>
// with concrete data
using Number = mp::number<mp::cpp_dec_float<50>, mp::et_off>;

using point = boost::geometry::model::d2::point_xy<Number>;
using polygon = bg::model::polygon<point, false>;

polygon create_poly(std::vector<std::tuple<Number, Number>> pts) {
    polygon poly;
    for (const auto& [x, y] : pts)
        bg::append(poly, bg::make<point>(x, y));
    auto [x_1, y_1] = pts[0];
    bg::append(poly, bg::make<point>(x_1, y_1));
    return poly;
}

bool perform_simple_union(const polygon& p1, const polygon& p2) {
    std::vector<polygon> output; 
    bg::union_(p1, p2, output);
    return output.size() == 1;
}

Number find_epsilon(Number left, Number right) {

    Number eps = (left + right) / 2;
    if (right - left < std::numeric_limits<Number>::epsilon())
        return left;

    polygon a = create_poly(
        std::vector<std::tuple<Number, Number>>{
            {1.0, 1.0}, { 2.0,1.0 }, { 2.0, 2.0 }, { 1.0,2.0 }
        }
    );

    polygon b = create_poly(
        std::vector<std::tuple<Number, Number>>{
            {2.0 + eps, 1.0}, { 3.0 + eps, 1.0 }, { 3.0 + eps, 2.0 }, { 2.0 + eps,2.0 }
        }
    );

    if ( perform_simple_union(a, b) ) {
        return find_epsilon(eps, right);
    } else {
        return find_epsilon(left, eps);
    }
}

int main()
{
    std::cout << "nextafter(0, 1):  " << nextafter(Number(0), Number(1)) << "\n";
    std::cout << "Number: eps()     " << std::numeric_limits<Number>::epsilon()      << "\n";
    std::cout << "Number: min_exp() " << std::numeric_limits<Number>::min_exponent10 << "\n";
    std::cout << "Number: max_exp() " << std::numeric_limits<Number>::max_exponent10 << "\n";
    std::cout << "Number: min()     " << std::numeric_limits<Number>::min()          << "\n";
    std::cout << "Number: max()     " << std::numeric_limits<Number>::max()          << "\n";

    auto eps = find_epsilon(0.0, 1.0);

    std::cout << std::setprecision(180);
    std::cout << "eps == " << eps << "\n";

    std::cout << std::boolalpha;
    std::cout << "zero? " << (eps == 0) << "\n";
}

打印

nextafter(0, 1):  1e-67108864
Number: eps()     1e-49
Number: min_exp() -67108864
Number: max_exp() 67108864
Number: min()     1e-67108864
Number: max()     1e+67108864
eps == 0
zero? true

对于 cpp_dec_float<0> 它打印(注意在可变精度情况下的“奇怪的” numeric_limits<>::eps`):
Live On Coliru

nextafter(0, 1):  1e-67108864
Number: eps()     1e-08
Number: min_exp() -67108864
Number: max_exp() 67108864
Number: min()     1e-67108864
Number: max()     1e+67108864
eps == 0
zero? true