时间戳

首页 > 解决方案 > Boost point_circle 以奇怪的形状出现

c++ - Boost point_circle 以奇怪的形状出现

问题描述

我正在尝试使用 Boost 的几何库在地球上创建一个 10m 半径的多边形。

这是教程。

为了编译这个例子，我使用了最新的 Clang 和 Boost 1.73.0 的 Wandbox 。

我首先在我的生产环境中发现了这个问题，即 Clang 12 和 Boost 1.71.0。

使用具有 32 个点的 1000m 半径圆产生预期结果：

在此处输入图像描述

然而，将其缩小到 10m 会产生意想不到的结果：

在此处输入图像描述

我使用WKT 游乐场来显示结果，并确认结果在其他可视化工具中是相同的。

这似乎是一个浮点舍入误差，但这里的一切都应该使用双精度浮点数，足以表示 GPS 坐标。计算似乎出了点问题。

使用半径为 0.0001 的boost::geometry::point_circle 也会发生同样的事情。

发生了什么事，我应该手动计算圆吗？

编辑 1

bg::area如果你用它来计算面积，那就更奇怪了。我试着画了一个半径为 10m 的圆POINT(4.9 52.1)，得到了 25984.4m 的面积。我尝试了同样的POINT(4.9 52.1000001)方法并得到了-1122.14。

请参阅以下操场：https ://godbolt.org/z/sTGqKK

编辑 2

我发现显示多边形的问题与计算面积不正确的问题是分开的。事实上，显示问题是打印到标准输出时四舍五入的结果。通过增加小数的精度，或使用std::fixed，解决了显示问题。

std::cout << std::fixed << bg::wkt(result) << std::endl;

标签： c++boostgisboost-geometry

解决方案

似乎确实存在准确性问题。我试图解决问题，但没有达到我想要的程度。

BGL 使用了一些硬限定std::abs和std::acos调用，这使得使用多精度类型变得困难。我尝试修补其中的一些，但一个下午的兔子洞太深了。

这是一个测试平台，可能有助于进一步查明/调试/跟踪事物。注意

  • 因为float准确性使得图书馆is_valid会因为尖峰而报告无效。
  • long double似乎做的合理

然而，首要问题（缺乏控制/可预测性）仍然存在。

Live On Compiler Explorer ¹

#include <boost/geometry.hpp>
#include <iostream>

#ifdef TRY_BOOST_MULTIPRECISION
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
    namespace bmp = boost::multiprecision;
    using OctFloat    = bmp::cpp_bin_float_oct;
    using Decimal     = bmp::number<bmp::cpp_dec_float<50>,  bmp::et_off>;
    using LongDecimal = bmp::number<bmp::cpp_dec_float<100>, bmp::et_off>;

    namespace boost::multiprecision {
        inline auto mod(OctFloat    const& a, OctFloat    const& b) { return fmod(a, b); }
        inline auto mod(Decimal     const& a, Decimal     const& b) { return fmod(a, b); }
        inline auto mod(LongDecimal const& a, LongDecimal const& b) { return fmod(a, b); }
        inline auto abs(OctFloat    const& a) { return fabs(a); }
        inline auto abs(Decimal     const& a) { return fabs(a); }
        inline auto abs(LongDecimal const& a) { return fabs(a); }
    }

    namespace std { // sadly BG overqualifies std::abs in places
        inline auto abs(OctFloat    const& a) { return fabs(a); }
    }
#endif

template <typename F, typename DegreeOrRadian>
void do_test(int n, F offset = {}) {
    namespace bg = boost::geometry;
    std::cout << "----- " << __PRETTY_FUNCTION__ << " n:" << n << " offset: " << offset << " ----\n";
    bg::model::point<F, 2, bg::cs::geographic<bg::degree> > Amsterdam { 4.9, 52.1 + offset };
    typedef bg::model::point<F, 2, bg::cs::geographic<DegreeOrRadian> > point;

    // Declare the geographic_point_circle strategy (with n points)
    // Default template arguments (taking Andoyer strategy)
    bg::strategy::buffer::geographic_point_circle<> point_strategy(n);

    // Declare the distance strategy (one kilometer, around the point, on Earth)
    bg::strategy::buffer::distance_symmetric<F> distance_strategy(10.0);

    // Declare other necessary strategies, unused for point
    bg::strategy::buffer::join_round    join_strategy;
    bg::strategy::buffer::end_round     end_strategy;
    bg::strategy::buffer::side_straight side_strategy;

    // Declare/fill a point on Earth, near Amsterdam
    point p;
    bg::convert(Amsterdam, p);

    // Create the buffer of a point on the Earth
    bg::model::multi_polygon<bg::model::polygon<point> > result;
    bg::buffer(p, result,
                distance_strategy, side_strategy,
                join_strategy, end_strategy, point_strategy);

    std::string reason;
    is_valid(result, reason);
    //std::cout << "result: " << wkt(result) << "\n";
    std::cout << reason << "\n";
    std::cout << "result: " << (bg::is_simple(result)?"simple":"compound") << "\n";

    auto area = bg::area(result);

    std::cout << "reference: " << bg::dsv(Amsterdam)  << std::endl;
    std::cout << "point: " << bg::dsv(p)  << std::endl;
    std::cout << "area: " <<  area << " m²" << std::endl;
}

int main() {
    for (long double offset : { 0.l/*, 1e-7l*/ }) {
        for (int n : { 36 }) {
            do_test<float,       boost::geometry::degree>(n, offset);
            do_test<double,      boost::geometry::degree>(n, offset);
            do_test<long double, boost::geometry::degree>(n, offset);

            do_test<float,       boost::geometry::radian>(n, offset);
            do_test<double,      boost::geometry::radian>(n, offset);
            do_test<long double, boost::geometry::radian>(n, offset);

            // not working yet
            //do_test<OctFloat,    boost::geometry::radian>(n, offset);
            //do_test<Decimal,     boost::geometry::degree>();
            //do_test<LongDecimal, boost::geometry::degree>();
        }
    }
}

印刷

----- void do_test(int, F) [F = float, DegreeOrRadian = boost::geometry::degree] n:36 offset: 0 ----
Geometry has spikes. A spike point was found with apex at (4.9, 52.0975)
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: -1.37916e+07 m²
----- void do_test(int, F) [F = double, DegreeOrRadian = boost::geometry::degree] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 25984.4 m²
----- void do_test(int, F) [F = long double, DegreeOrRadian = boost::geometry::degree] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 301.264 m²
----- void do_test(int, F) [F = float, DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry has spikes. A spike point was found with apex at (-1.38318, -1.30708)
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 1.85308e+08 m²
----- void do_test(int, F) [F = double, DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 6399.41 m²
----- void do_test(int, F) [F = long double, DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 302.318 m²

在我的机器上

¹ 超过处理时间

推荐阅读