Boost point_circle 以奇怪的形状出现

编程问答 2022-06-01

问题描述

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

这是 tutorial。

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

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

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

enter image description here

将其缩小到 10m 却产生了意想不到的结果：

enter image description here

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

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

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

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

编辑 1

如果您使用 bg::area 来计算面积，它会变得更加奇怪。我尝试在 POINT(4.9 52.1) 周围绘制一个半径为 10 米的圆，得到的面积为 25984.4 米。我在 POINT(4.9 52.1000001) 尝试了同样的方法，结果是 -1122.14。

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

编辑 2

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

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

解决方法

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

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,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,offset);
            do_test<long double,offset);

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

            // not working yet
            //do_test<OctFloat,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)
area: 25984.4 m²
----- void do_test(int,F) [F = long double,52.1)
area: 301.264 m²
----- void do_test(int,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)
area: 1.85308e+08 m²
----- void do_test(int,DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9,52.1)
area: 6399.41 m²
----- void do_test(int,52.1)
area: 302.318 m²

在我的机器上

¹ 超过处理时间

,

据我所知，不准确的原因有两个，区域算法和点算法上的地理缓冲区。

关于前一个 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,F radius,F offset = {}) {
    namespace bg = boost::geometry;

    std::cout
        << "----- " << __PRETTY_FUNCTION__
        << " n:" << n << " radius:" << radius << " offset:" << offset
        << " ----\n";

    bg::model::point<F,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,on Earth)
    bg::strategy::buffer::distance_symmetric<F> distance_strategy(radius);

    // Declare other necessary strategies,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,10.);
    do_test<long double,10.);

    do_test<double,10.,offset);
    do_test<long double,offset);

    do_test<double,1000.);
    do_test<double,1000.,1.);
    do_test<long double,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)
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) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:0 ----
reference: (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) [with F = double; DegreeOrRadian = boost::geometry::radian] n:8 radius:10 offset:0 ----
reference: (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) [with F = long double; DegreeOrRadian = boost::geometry::radian] n:8 radius:10 offset:0 ----
reference: (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) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:1e-07 ----
reference: (4.9,52.1)
radius: 10 m
area: 281.749 m²,error: -0.0283973 m²,normalised error: 0.010039%
area: 281.749 m²,error: -0.0284534 m²,normalised error: 0.0100588%
----- void do_test(int,F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:1e-07 ----
reference: (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) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1000 offset:0 ----
reference: (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²,normalised error: 0.0100494%
----- void do_test(int,F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1000 offset:1e-07 ----
reference: (4.9,error: -284.372 m²,normalised error: 0.010053%
area: 2.82843e+06 m²,error: -284.282 m²,normalised error: 0.0100499%
area: 2.82843e+06 m²,F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1 offset:0 ----
reference: (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) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1 offset:0 ----
reference: (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%

一些评论：

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

相关问答

matplotlib报错：AttributeError: module 'backend_interagg' has no attribute 'FigureCanvas'. Did you mean: 'FigureCanvasAgg'?
使用本地python环境可以成功执行 import pandas as pd impor...
gitlab登录失败，报错：This challenge page was accidentally cached by an intermediary and is no longer available.
设置时间 控制面板
后端开发常见错误
错误1：Request method ‘DELETE‘ not supported 错误还原：...
docker常见错误
错误1：启动docker镜像时报错：Error response from daemon:...
idea常见错误
错误1：private field ‘xxx‘ is never assigned 按Alt...
pip安装依赖失败
报错如下，通过源不能下载，最后警告pip需升级版本 Requirem...