boost::geometry::intersection fails

[price@…]

I expect this has something to do with floating-point tolerances. This bug appears to be similar to #6958 (from which I adapted the below code), but must not be identical since that was fixed in 1.55.0. Hopefully the same fix can be applied easily.

$ cat polygonBug.cc
#include <vector>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/geometries/register/ring.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>

typedef boost::geometry::model::d2::point_xy<double> pt;
typedef boost::geometry::model::polygon<pt> polygon;
typedef boost::geometry::model::multi_polygon<polygon> multi_polygon;

int main()
{
    polygon green, blue;
    multi_polygon output;

    boost::geometry::read_wkt("POLYGON(("
                              " 14853.584999999999 20481.394 ,"
                              " 16484.061000000002 20466.799999999999 ,"
                              " 14853.584000000001 20481.392 ,"
                              " 14853.584999999999 20481.394))", green);

    boost::geometry::read_wkt("POLYGON(("
                              " 16483.432000000001 20553.794000000002 ,"
                              " 16483.432000000001 20553.793000000001 ,"
                              " 16484.060000000001 20466.799999999999 ,"
                              " 16484.061000000002 20466.799999999999 ,"
                              " 14853.584999999999 20481.394 ,"
                              " 14854.127 20520.719000000001 ,"
                              " 16483.432000000001 20553.794000000002))", blue);

    boost::geometry::correct(green);
    boost::geometry::correct(blue);
    boost::geometry::intersection(green, blue, output);

    return 0;
}

$ clang++ -o polygonBug polygonBug.cc -I$BOOST_DIR/include
$ ./polygonBug 
terminate called throwing an exceptionAbort trap: 6

Adding the following before using boost::geometry::intersection allows the call to succeed:

    boost::geometry::simplify(green, green, 0.001);
    boost::geometry::simplify(blue, blue, 0.001);

Relevant entries in the lldb backtrace:

    frame #8: 0x0000000100002941 polygonBug`std::back_insert_iterator<boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator> > boost::geometry::detail::overlay::overlay<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, false, false, false, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, (boost::geometry::overlay_type)1>::apply<std::back_insert_iterator<boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator> >, boost::geometry::strategy_intersection<boost::geometry::cartesian_tag, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, void> >(boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator> const&, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator> const&, std::back_insert_iterator<boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator> >, boost::geometry::strategy_intersection<boost::geometry::cartesian_tag, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, void> const&) + 385 at overlay.hpp:230
    frame #9: 0x000000010000279d polygonBug`bool boost::geometry::dispatch::intersection<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::polygon_tag, boost::geometry::polygon_tag, false>::apply<boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator>, boost::geometry::strategy_intersection<boost::geometry::cartesian_tag, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, void> >(boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator> const&, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator> const&, boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator>&, boost::geometry::strategy_intersection<boost::geometry::cartesian_tag, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, void> const&) + 77 at intersection.hpp:103
    frame #10: 0x00000001000026e3 polygonBug`bool boost::geometry::intersection<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator> >(boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator> const&, boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator> const&, boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true, std::vector, std::vector, std::allocator, std::allocator>, std::vector, std::allocator>&) + 51 at intersection.hpp:198

[Barend Gehrels]

The second polygon (blue) is invalid. Boost.Geometry raises an exception for that.

It is also reported as invalid by both SQL Server and PostGIS so we can conclude it is really invalid.

"Self-intersection[16484.1 20466.8]"

I close the ticket therefore as "invalid" too.