Boost point_circle 以奇怪的形状出现

Question

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

这是tutorial。

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

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

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

但是将其缩小到 10m 会产生意想不到的结果：

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

这似乎是一个浮点舍入错误，但这里的一切都应该使用双精度浮点数，即more than enough to represent GPS coordinates。计算似乎出了点问题。

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

这是怎么回事，我应该手动计算圆吗？

编辑 1

如果你使用bg::area 来计算面积，那就更奇怪了。我尝试在POINT(4.9 52.1) 周围绘制一个'10m' 半径的圆，并得到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;

Answer 1

据我了解，有两个不准确的来源，区域算法和点上的地理缓冲区算法。

关于前者https://github.com/boostorg/geometry/pull/801 提出了一个修复方案。使用此修复程序，上述误差函数 (godbolt.org/z/sTGqKK) 返回小于 1% 的相对误差。以下代码使用策略扩展了这一点。

#include <boost/geometry.hpp>
#include <cmath>
#include <iostream>
 
template <typename CT>
void error_function(CT area, CT theoreticalArea)
{
    std::cout << "area: " <<  area << " m², ";
    std::cout << "error: " <<  area - theoreticalArea << " m²,\t";
    std::cout << "normalised error: " <<  fabs(100 * (area - theoreticalArea)
        / theoreticalArea) << "%" << std::endl;
}

template <typename F, typename DegreeOrRadian>
void do_test(int n, F radius, F offset = {}) {
    namespace bg = boost::geometry;

    std::cout
        << "----- " << __PRETTY_FUNCTION__
        << " n:" << n << " radius:" << radius << " 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 (ten metres, around the point, on Earth)
    bg::strategy::buffer::distance_symmetric<F> distance_strategy(radius);

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

    auto area = bg::area(result);
    auto areat = bg::area(result,bg::strategy::area::geographic<bg::strategy::thomas>());
    auto areav = bg::area(result,bg::strategy::area::geographic<bg::strategy::vincenty>());
    auto areak = bg::area(result,bg::strategy::area::geographic<bg::strategy::karney>());

    // Assumes that the Earth is flat, which it clearly is.
    // A = n/2 * R^2 * sin(2*pi/n) where R is the circumradius
    auto theoreticalArea = n * radius * radius * std::sin(2.0 * 3.142 / n) / 2.0;

    std::cout << "reference: " << bg::dsv(Amsterdam)  << std::endl;
    std::cout << "point: " << bg::dsv(p)  << std::endl;
    std::cout << "radius: " <<  radius << " m" << std::endl;
    error_function(area, theoreticalArea);
    error_function(areat, theoreticalArea);
    error_function(areav, theoreticalArea);
    error_function(areak, theoreticalArea);
}

int main() {
    double offset = 1e-7;
    int n = 8;

    do_test<double,      boost::geometry::degree>(n, 10.);
    do_test<long double, boost::geometry::degree>(n, 10.);

    do_test<double,      boost::geometry::radian>(n, 10.);
    do_test<long double, boost::geometry::radian>(n, 10.);

    do_test<double,      boost::geometry::degree>(n, 10., offset);
    do_test<long double, boost::geometry::degree>(n, 10., offset);

    do_test<double,      boost::geometry::degree>(n, 1000.);
    do_test<double,      boost::geometry::degree>(n, 1000., offset);

    do_test<double,      boost::geometry::degree>(n, 1.);
    do_test<long double, boost::geometry::degree>(n, 1.);
}

返回：

----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 281.272 m²,   error: -1.59991 m²,     normalised error: 0.565596%
area: 282.843 m²,   error: -0.0284134 m²,   normalised error: 0.0100446%
area: 281.749 m²,   error: -1.12206 m²,     normalised error: 0.396666%
area: 282.843 m²,   error: -0.028415 m²,    normalised error: 0.0100452%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 283.57 m²,    error: 0.698736 m²,     normalised error: 0.247015%
area: 282.843 m²,   error: -0.0284201 m²,   normalised error: 0.010047%
area: 283.568 m²,   error: 0.696594 m²,     normalised error: 0.246258%
area: 282.843 m²,   error: -0.0284255 m²,   normalised error: 0.0100489%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::radian] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 282.715 m²,   error: -0.156633 m²,    normalised error: 0.0553726%
area: 282.843 m²,   error: -0.0286857 m²,   normalised error: 0.0101409%
area: 280.578 m²,   error: -2.29311 m²,     normalised error: 0.810656%
area: 282.843 m²,   error: -0.0286896 m²,   normalised error: 0.0101423%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::radian] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 283.135 m²,   error: 0.263058 m²,     normalised error: 0.0929955%
area: 282.843 m²,   error: -0.0287086 m²,   normalised error: 0.010149%
area: 283.164 m²,   error: 0.292786 m²,     normalised error: 0.103505%
area: 282.843 m²,   error: -0.0287018 m²,   normalised error: 0.0101466%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:1e-07 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 281.749 m²,   error: -1.12206 m²,     normalised error: 0.396666%
area: 282.843 m²,   error: -0.0283973 m²,   normalised error: 0.010039%
area: 281.749 m²,   error: -1.12206 m²,     normalised error: 0.396666%
area: 282.843 m²,   error: -0.0284534 m²,   normalised error: 0.0100588%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:1e-07 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 283.569 m²,   error: 0.697826 m²,     normalised error: 0.246694%
area: 282.843 m²,   error: -0.0284078 m²,   normalised error: 0.0100427%
area: 283.568 m²,   error: 0.696529 m²,     normalised error: 0.246235%
area: 282.843 m²,   error: -0.0283946 m²,   normalised error: 0.010038%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1000 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1000 m
area: 2.82843e+06 m²,   error: -284.28 m²,  normalised error: 0.0100498%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
area: 2.82843e+06 m²,   error: -284.259 m², normalised error: 0.0100491%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1000 offset:1e-07 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1000 m
area: 2.82843e+06 m²,   error: -284.372 m², normalised error: 0.010053%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
area: 2.82843e+06 m²,   error: -284.282 m², normalised error: 0.0100499%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1 m
area: 2.81749 m²,   error: -0.0112205 m²,   normalised error: 0.396663%
area: 2.8285 m²,    error: -0.000219998 m², normalised error: 0.0077773%
area: 2.83391 m²,   error: 0.0051987 m²,    normalised error: 0.183783%
area: 2.82848 m²,   error: -0.000234082 m², normalised error: 0.00827521%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1 m
area: 2.83535 m²,   error: 0.00663946 m²,   normalised error: 0.234717%
area: 2.82844 m²,   error: -0.000278463 m², normalised error: 0.00984417%
area: 2.83392 m²,   error: 0.005205 m²,     normalised error: 0.184006%
area: 2.82842 m²,   error: -0.000294424 m², normalised error: 0.0104084%

一些cmets：

• 使用不同的策略（即在 boost 几何中执行地理计算的算法）控制着算法的准确性和性能。
• 地理缓冲区仍有问题，请随时在 github 上提交问题以保持跟踪
• “theoreticalArea”仅适用于小区域，随着区域的增长，预计提升几何算法会比该区域更准确。
• 地球不是平的；）

Answer 2

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

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²

在我的机器上

---

¹超过处理时间