Boost point_circle coming out in weird shapes

3
votes

I'm attempting to create a 10m radius polygon on Earth using Boost's Geometry library.

Here's the tutorial.

To compile this example, I used Wandbox with the latest Clang and Boost 1.73.0.

I first discovered the issue in my production environment, which is Clang 12 and Boost 1.71.0.

Using a 1000m radius circle with 32 points yields expected results:

enter image description here

Shrinking that down to 10m however has unexpected results:

enter image description here

I used a WKT playground to display the results, and have confirmed that results are the same in other visualisation tools.

It seems to be a floating point rounding error, but everything here should be using double-precision floats which are more than enough to represent GPS coordinates. Something seems to be going wrong with the calculation.

The same thing happens with boost::geometry::point_circle using a radius of 0.0001.

What's going on, and should I just calculate the circle manually instead?

Edit 1

It gets even weirder if you use bg::area to calculate the area. I tried on a '10m' radius circle drawn around POINT(4.9 52.1) and got an area of 25984.4m. I tried the same at POINT(4.9 52.1000001) and got -1122.14.

See the following playground: https://godbolt.org/z/sTGqKK

Edit 2

I discovered that the issue with displaying the polygon is separate to the issue of the calculated area being incorrect. In fact, the display issue is as a result of rounding when printing to stdout. By increasing the precision of decimals, or using std::fixed, the display issue is resolved.

std::cout << std::fixed << bg::wkt(result) << std::endl;
2
c++boostgisboost-geometry
Further details on how boost::geometry::buffer works: boost.org/doc/libs/master/libs/geometry/doc/html/geometry/… – Chris Watts
If I read this I don't think I see actual numbers for double precision, but the single precision accuracies list worst case 1.64-2.37m accuracy. If you divide a circle of radius 10m to that kind of precision, it's like drawing a circle r=5 on an integer grid (10m/~2m = 5). That would explain a shape like you show. That the "resolution is enough for GPS" is hardly relevant because you couldn't use GPS to accurately outline a 10m radius circle either. – sehe
@sehe Yes, what we see is almost as if it's displaying a float. The machine epsilon for double is 2.22e-16, which according to this table ought to have a scale precision of much less than 1mm (2.471304e-11 metres to be precise). – Chris Watts

2 Answers

2
votes

To my understanding there are two sources of inaccuracy, the area algorithm and the geographic buffer over a point algorithm.

Regarding the former https://github.com/boostorg/geometry/pull/801 proposes a fix. Using this fix the above error function (godbolt.org/z/sTGqKK) returns less than 1% relative error. The following piece of code extends this by using strategies.

#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.);
}

which 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%

Some comments:

  • the use of different strategies (i.e. algorithms that perform the geographic computations in boost geometry) controls the accuracy but also the performance of algorithms.
  • there is still an issue in geographic buffer, feel free to file an issue on github to keep it in track
  • the "theoreticalArea" is accurate only for small areas, as the area grows the boost geometry algorithms are expected to be more accurate than that area.
  • the Earth is not flat ;)
2
votes

It does seem that there are accuracy issues. I tried to work around things but didn't get as far as I'd like.

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.

Here's a testbed that might help pinpoint/debug/trace things a bit further. Note that

  • for float the accuracy is such that the library is_valid will report invalid due to spikes.
  • long double seems to do reasonable

The over-arching problem (lack of control/predictability) however remains.

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>();
        }
    }
}

Prints

----- 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²

On my machine

¹ exceeds processing time