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

How do Boost.Geometry operations like union_ deal with the fundamental imprecision of floating point types?

我正在尝试判断是否/如何让 Boost.Geometry 适用于特定用例。但是,我无法在任何地方找到有关库如何处理浮点类型的文档。

如果你在官方文档中搜索“epsilon”这个词,据我所知,你得到的结果为零;然而,从库的行为中可以清楚地看出,在进行比较时,它隐含地使用了处理浮点数的典型方式的某种版本,因为例如,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。因此,如果双精度坐标彼此在单精度 epsilon 范围内,那么它会将双精度坐标视为等效的?

基本上我只想知道默认行为是什么,这样我就可以决定它是否适合我。

在[介绍][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