Boost, geometry

Question

I have a problem with boost::geomentry.

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <vector>

int main(){
typedef boost::geometry::model::d2::point_xy<double> TBoostPoint;
typedef boost::geometry::model::polygon<TBoostPoint> TBoostPoly;    
TBoostPoly square, square1;    
square.outer().push_back(TBoostPoint(0.5,4.25));
square.outer().push_back(TBoostPoint(0.5,4.5));
square.outer().push_back(TBoostPoint(1.0,4.5));
square.outer().push_back(TBoostPoint(1.0,4.25));
square.outer().push_back(TBoostPoint(0.5,4.25));    
const double eps[] = {1e-15,1e-15,2e-15,2e-15};    
square.outer().push_back(TBoostPoint(0.5,4.25 + eps[0]));
square.outer().push_back(TBoostPoint(0.5,4.5  + eps[1]));
square.outer().push_back(TBoostPoint(1.0,4.5  + eps[2]));
square.outer().push_back(TBoostPoint(1.0,4.25 + eps[3]));
square.outer().push_back(TBoostPoint(0.5,4.25 + eps[0]));    
boost::geometry::correct(square);
boost::geometry::correct(square1);    
std::vector<TBoostPoly> output;    
boost::geometry::intersection(square,square1,output);    
return 0;    
}

If I use Boost(1_58) output - is not correct, but if I use Boost(1_55 or 1_52), output - is correct.

Expected

{(0.5,4.25),(0.5,4.5),(1.0,4.25),(1.0,4.25),(0.5,4.25)}

Result (boost 1_58)

{(0.5,4.25),(0.5,4.5),(1.0,4.25),(1.0,4.25 + 5e-8),(0.5,4.25)}.

Might care to add what you'd expect (and what you get). Not all of us have "ancient" Boost installations — WorldSEnder, Jul 08 '15 at 16:11
You should use integral values. Have a look at http://stackoverflow.com/questions/20080868/segmentation-fault-with-boostpolygon/20081230#20081230 — , Jul 08 '15 at 16:24

score 1 · Answer 1 · answered Jul 08 '15 at 17:16

You have to use integral coordinates.

From the documentation: http://www.boost.org/doc/libs/1_58_0/libs/polygon/doc/index.htm

The coordinate data type is a template parameter of all data types and algorithms provided by the library, and is expected to be integral. Floating point coordinate data types are not supported by the algorithms implemented in the library due to the fact that the achieving floating point robustness implies a different set of algorithms and generally platform specific assumptions about floating point representations.

Same applies to earlier versions.

In your case the output of Boost(1_55 or 1_52) is correct (by accident).

＋1 for addressing the essence. My answer was just curiosity. This is the answer to be accepted — sehe, Jul 08 '15 at 22:52

sehe · Answer 2 · 2015-07-08T22:52:00.163

The output may seem correct the first around, but in fact it isn't if you look closer:

A slightly refactored sample: Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/array.hpp>
#include <vector>
#include <iomanip>

namespace bg = boost::geometry;

template <typename C, typename T = typename C::value_type>
    void append(C& container, std::initializer_list<T> init) {
        container.insert(container.end(), init);
    }

int main() {
    typedef bg::model::d2::point_xy<double> TBoostPoint;
    typedef bg::model::polygon<TBoostPoint> TBoostPoly;

    std::vector<TBoostPoly> squares;

    using Eps = boost::array<double, 4>;

    for (auto const& eps : { 
                Eps {{     0,     0,     0,     0 }},
                Eps {{ 1e-15, 1e-15, 2e-15, 2e-15 }},
            })
    {
        TBoostPoly square;
        append(square.outer(), {
             { 0.5,  4.25 + eps[0] },
             { 0.5,   4.5 + eps[1] },
             { 1.0,   4.5 + eps[2] },
             { 1.0,  4.25 + eps[3] },
             { 0.5,  4.25 + eps[0] }
            });

        squares.push_back(std::move(square));
    }

    for (auto& p : squares)
        bg::correct(p);

    std::vector<TBoostPoly> output;
    bg::intersection(squares[0], squares[1], output);

    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
    std::cout << std::fixed << std::setprecision(std::numeric_limits<bg::coordinate_type<TBoostPoint>::type >::max_digits10);
    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
}

Which prints

Output: POLYGON((0.5 4.5,1 4.5,1 4.25,0.5 4.25,0.5 4.5))
Output: POLYGON((0.50000000000000000 4.50000000000000000,1.00000000000000000 4.50000000000000000,1.00000000000000000 4.25000000000000178,0.50000000000000000 4.25000004999999970,0.50000000000000000 4.50000000000000000))

As you can see, the naive, natural output may seem to be 4.25 at some point, but the actual value stored is 4.25000000000000178 at that exact moment.

Depending on what requirements you have you might be happier with some arbitrary precision decimal representation types. As a proof of concept, here's the same program parameterized to use 50-digit decimal floats:

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/array.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <vector>
#include <iomanip>

namespace bg  = boost::geometry;
namespace bmp = boost::multiprecision;

template <typename C, typename T = typename C::value_type>
    void append(C& container, std::initializer_list<T> init) {
        container.insert(container.end(), init);
    }

int main() {
    typedef bmp::number<bmp::cpp_dec_float<50>, bmp::et_off> Decimal;
    typedef bg::model::d2::point_xy<Decimal> TBoostPoint;
    typedef bg::model::polygon<TBoostPoint> TBoostPoly;

    std::vector<TBoostPoly> squares;

    using Eps = boost::array<Decimal, 4>;

    for (auto const& eps : { 
                Eps {{     0,     0,     0,     0 }},
                Eps {{ 1e-15, 1e-15, 2e-15, 2e-15 }},
            })
    {
        TBoostPoly square;
        append(square.outer(), {
             { 0.5,  4.25 + eps[0] },
             { 0.5,   4.5 + eps[1] },
             { 1.0,   4.5 + eps[2] },
             { 1.0,  4.25 + eps[3] },
             { 0.5,  4.25 + eps[0] }
            });

        squares.push_back(std::move(square));
    }

    for (auto& p : squares)
        bg::correct(p);

    std::vector<TBoostPoly> output;
    bg::intersection(squares[0], squares[1], output);

    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
    std::cout << std::fixed << std::setprecision(std::numeric_limits<bg::coordinate_type<TBoostPoint>::type >::max_digits10);
    for (auto& p : output) std::cout << "Output: " << bg::wkt(p) << "\n";
}

Which prints:

Output: POLYGON((0.5 4.5,1 4.5,1 4.25,0.5 4.25,0.5 4.5))
Output: POLYGON((0.50000000000000000000000000000000000000000000000000000000000000000000000000000000 4.50000000000000000000000000000000000000000000000000000000000000000000000000000000,1.00000000000000000000000000000000000000000000000000000000000000000000000000000000 4.50000000000000000000000000000000000000000000000000000000000000000000000000000000,1.00000000000000000000000000000000000000000000000000000000000000000000000000000000 4.25000000000000200000000000000015541079975332215847661437120239003029098500000000,0.50000000000000000000000000000000000000000000000000000000000000000000000000000000 4.25000000000000100000000000000007770539987666107923830718560119501514549200000000,0.50000000000000000000000000000000000000000000000000000000000000000000000000000000 4.50000000000000000000000000000000000000000000000000000000000000000000000000000000))

I realize this answer doesn't add essential information beyond the existing answer. For the curious the Multiprecision example might be informative, though. — sehe, Jul 08 '15 at 22:01

Boost, geometry

2 Answers2