Boost point_circle 以奇怪的形状出现
Boost point_circle coming out in weird shapes
我正在尝试使用 Boost 的几何库在地球上创建一个 10m 半径的多边形。
这是 tutorial。
为了编译此示例，我使用了 Wandbox 以及最新的 Clang 和 Boost 1.73.0。
我首先在我的生产环境中发现了这个问题，即 Clang 12 和 Boost 1.71.0。
使用具有 32 个点的 1000 米半径圆产生预期结果：
然而，将其缩小到 10m 会产生意想不到的结果：
我用了一个WKT playground来显示结果，并且已经确认在其他可视化工具中结果是一样的。
这似乎是一个浮点舍入错误，但这里的所有内容都应该使用 more than enough to represent GPS coordinates 双精度浮点数。计算好像出了点问题。
使用半径 0.0001 的 boost::geometry::point_circle 也会发生同样的事情。
这是怎么回事，我应该手动计算圆吗？
编辑 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;
似乎确实存在准确性问题。我试图解决问题，但没有达到我想要的程度。
BGL uses some hard-qualified std::abs and std::acos calls that make it hard to use multiprecision types. I tried patching some of these, but the rabit hole was too deep for an afternoon.
这是一个可能有助于 pinpoint/debug/trace 更进一步的测试平台。注意
- 对于
float 的准确性使得库 is_valid 将由于尖峰而报告无效。
long double好像做的合理
然而，总体问题（缺少 control/predictability）仍然存在。
#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²
在我的机器上
¹ 超过处理时间
据我所知，有两个不准确的来源，区域算法和点算法上的地理缓冲区。
关于前者 https://github.com/boostorg/geometry/pull/801 提出修复。使用此修复上述误差函数 (godbolt.org/z/sTGqKK) returns 小于 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.);
}
哪个returns:
----- 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%
一些评论：
- 使用不同的策略（即在增强几何中执行地理计算的算法）控制着算法的准确性和性能。
- 地理缓冲区中仍然存在问题，请随时在 github 上提交问题以保持跟踪
- “theoreticalArea”仅适用于小区域，随着区域的增长，boost 几何算法预计会比该区域更准确。
- 地球不平;)
我正在尝试使用 Boost 的几何库在地球上创建一个 10m 半径的多边形。
这是 tutorial。
为了编译此示例，我使用了 Wandbox 以及最新的 Clang 和 Boost 1.73.0。
我首先在我的生产环境中发现了这个问题，即 Clang 12 和 Boost 1.71.0。
使用具有 32 个点的 1000 米半径圆产生预期结果：
然而，将其缩小到 10m 会产生意想不到的结果：
我用了一个WKT playground来显示结果，并且已经确认在其他可视化工具中结果是一样的。
这似乎是一个浮点舍入错误，但这里的所有内容都应该使用 more than enough to represent GPS coordinates 双精度浮点数。计算好像出了点问题。
使用半径 0.0001 的 boost::geometry::point_circle 也会发生同样的事情。
这是怎么回事，我应该手动计算圆吗？
编辑 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;
似乎确实存在准确性问题。我试图解决问题，但没有达到我想要的程度。
BGL uses some hard-qualified
std::absandstd::acoscalls that make it hard to use multiprecision types. I tried patching some of these, but the rabit hole was too deep for an afternoon.
这是一个可能有助于 pinpoint/debug/trace 更进一步的测试平台。注意
- 对于
float的准确性使得库is_valid将由于尖峰而报告无效。 long double好像做的合理
然而，总体问题（缺少 control/predictability）仍然存在。
#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²
在我的机器上
¹ 超过处理时间
据我所知，有两个不准确的来源，区域算法和点算法上的地理缓冲区。
关于前者 https://github.com/boostorg/geometry/pull/801 提出修复。使用此修复上述误差函数 (godbolt.org/z/sTGqKK) returns 小于 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.);
}
哪个returns:
----- 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%
一些评论：
- 使用不同的策略（即在增强几何中执行地理计算的算法）控制着算法的准确性和性能。
- 地理缓冲区中仍然存在问题，请随时在 github 上提交问题以保持跟踪
- “theoreticalArea”仅适用于小区域，随着区域的增长，boost 几何算法预计会比该区域更准确。
- 地球不平;)