boost::geometry::union_ produces self intersections

Question

I have two valid polygons. When I take their union, I am getting an invalid polygon (there are self intersections). Is this a bug? I would expect that the union operation would always produce a valid polygon. I have provided an example below along with visualizations. Can anyone explain why this isn't a bug, or explain if there's a way to fix it?

#include <fstream>
#include <iostream>

#include <boost/geometry.hpp> // read_wkt
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <boost/geometry/algorithms/union.hpp>

using PointType = boost::geometry::model::d2::point_xy<double>;
using PolygonType = boost::geometry::model::polygon<PointType>;
using MultiPolygonType = boost::geometry::model::multi_polygon<PolygonType>;

template <typename TPolygon>
void WritePolygonsToSVG(const std::vector<TPolygon>& polygons, const std::string& filename)
{
    std::ofstream svg(filename);

    boost::geometry::svg_mapper<PointType> mapper(svg, 400, 400);

    for(unsigned int i = 0; i < polygons.size(); ++i) {
        mapper.add(polygons[i]);
        mapper.map(polygons[i], "fill:rgb(255,128,0);stroke:rgb(0,0,100);stroke-width:1");
    }

}

int main(int, char**)
{
    /// Create two polygons
    PolygonType singlePolygon1;
    PolygonType singlePolygon2;

    boost::geometry::read_wkt("POLYGON((-52.8018 -26.744,-57.5465 -27.9916,-62.2844 -29.2642,-63.19 -26.066,-57.564 -24.5243,-53.7273 -23.3394,-52.8018 -26.744))", singlePolygon1);
    boost::geometry::read_wkt("POLYGON((-56.9695 -33.6407,-58.009 -33.0132,-58.5119 -32.8011,-59.0182 -32.6562,-59.5392 -32.5779,-60.0858 -32.5658,"
                              "-61.3005 -32.7384,-62.2844 -29.2642,-59.997 -28.7264,-58.0701 -28.125,-56.7078 -27.7124,-55.3109 -27.4556,-55.0955 -27.4599,"
                              "-54.8762 -27.503,-54.6545 -27.5806,-54.4318 -27.6887,-53.9893 -27.98,-53.5601 -28.344,-52.7879 -29.1592,-52.207 -29.8722,-56.9695 -33.6407))", singlePolygon2);

    boost::geometry::correct(singlePolygon1);
    boost::geometry::correct(singlePolygon2);

    /// Run union and check validity (should fail check)
    MultiPolygonType unionResult;
    boost::geometry::union_(singlePolygon1, singlePolygon2, unionResult);

    boost::geometry::validity_failure_type failure_type;
    if(!boost::geometry::is_valid(unionResult, failure_type)) {
        std::cout << "Result of union operation is not valid! " << failure_type << std::endl; // failure_type is 21 == failure_self_intersections
    }

    // Put these into the same type so that they can be passed to WritePolygonsToSVG
    MultiPolygonType polygon1;
    polygon1.push_back(singlePolygon1);

    MultiPolygonType polygon2;
    polygon2.push_back(singlePolygon2);

    std::vector<MultiPolygonType> polygons = {polygon1, polygon2, unionResult};
    WritePolygonsToSVG(polygons, "allPolygons.svg");

    return EXIT_SUCCESS;
}

Polygon1:

Polygon2:

Both input polygons:

Union polygon (self intersections):

The actual problem that I'm running into is that if I then try to union the (invalid) output polygon with another polygon, it throws an exception:

terminate called after throwing an instance of 'boost::geometry::overlay_invalid_input_exception'
  what():  Boost.Geometry Overlay invalid input exception

This expanded example demonstrates the exception:

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

#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/algorithms/intersection.hpp>
#include <boost/geometry/algorithms/union.hpp>


int main(int argc, char** argv) {

    using PointType = boost::geometry::model::d2::point_xy<double>;
    using PolygonType = boost::geometry::model::polygon<PointType>;
    using MultiPolygonType = boost::geometry::model::multi_polygon<PolygonType>;

    /// Read in two polygons
    PolygonType polygon1;
    PolygonType polygon2;

    boost::geometry::read_wkt("POLYGON((-52.8018 -26.744,-57.5465 -27.9916,-62.2844 -29.2642,-63.19 -26.066,-57.564 -24.5243,-53.7273 -23.3394,-52.8018 -26.744))", polygon1);
    boost::geometry::read_wkt("POLYGON((-56.9695 -33.6407,-58.009 -33.0132,-58.5119 -32.8011,-59.0182 -32.6562,-59.5392 -32.5779,-60.0858 -32.5658,"
                              "-61.3005 -32.7384,-62.2844 -29.2642,-59.997 -28.7264,-58.0701 -28.125,-56.7078 -27.7124,-55.3109 -27.4556,-55.0955 -27.4599,"
                              "-54.8762 -27.503,-54.6545 -27.5806,-54.4318 -27.6887,-53.9893 -27.98,-53.5601 -28.344,-52.7879 -29.1592,-52.207 -29.8722,-56.9695 -33.6407))", polygon2);

    /// Run union and check validity (should fail check)
    MultiPolygonType unionResult1;
    boost::geometry::union_(polygon1, polygon2, unionResult1);

    boost::geometry::validity_failure_type failure_type;
    if(!boost::geometry::is_valid(unionResult1, failure_type)) {
        std::cout << "Result of union operation is not valid!" << std::endl;
    }

    /// Perform second union with third polygon (should throw)
    MultiPolygonType polygon3;
    boost::geometry::read_wkt("MULTIPOLYGON(((-36.5181 -22.1798,-39.072 -22.991,-41.6641 -23.6833,-52.8018 -26.744,-53.7273 -23.3394,-45.7888 -21.0022,"
                              "-41.8671 -19.6963,-37.9783 -18.2852,-36.5181 -22.1798)),((-52.207 -29.8722,-51.6368 -30.7643,-51.1556 -31.8424,-50.7568 -33.0527,"
                              "-50.4336 -34.3414,-49.9875 -36.9394,-49.8512 -38.1412,-49.7639 -39.2066,-54.0482 -39.4903,-54.1623 -38.8728,-54.4115 -38.1178,"
                              "-55.1907 -36.4023,-56.1365 -34.758,-56.594 -34.092,-57.0005 -33.5988,-57.0016 -33.5993,-57.0036 -33.6033,-57.0033 -33.6044,"
                              "-57.0029 -33.6055,-57.0014 -33.6079,-56.9995 -33.6106,-56.9941 -33.6165,-56.9695 -33.6407,-52.207 -29.8722)),((-61.3005 -32.7384,"
                              "-65.9902 -34.1028,-70.659 -35.5553,-75.2871 -37.1219,-79.8551 -38.829,-84.6551 -40.8483,-87.0116 -41.984,-89.2911 -43.2272,-92.5791 -37.6053,"
                              "-89.0742 -35.6754,-85.5037 -33.9031,-78.2194 -30.6022,-74.8627 -29.2273,-71.9907 -28.2985,-69.0648 -27.5347,-63.19 -26.066,-61.3005 -32.7384)),"
                              "((75.3002 -38.7765,72.918 -39.2671,70.5558 -39.533,68.2188 -39.5707,65.912 -39.3766,63.6406 -38.9473,61.4095 -38.2792,59.2239 -37.3688,"
                              "57.089 -36.2125,55.5367 -35.1756,54.0721 -34.0078,52.7036 -32.7209,51.4398 -31.3265,50.2893 -29.8362,49.2607 -28.2616,48.3624 -26.6144,"
                              "47.6032 -24.9062,47.158 -23.7161,46.803 -22.6211,46.3142 -20.5856,46.0379 -18.5376,45.9485 -17.427,45.8751 -16.215,45.5849 -10.2471,"
                              "45.4467 -7.25562,45.2962 -4.26324,45.0667 0.0336182,44.8334 4.32845,44.3738 10.6768,42.9884 29.2827,42.3055 38.5987,41.6559 47.9296,"
                              "45.1275 48.1622,45.7738 38.8379,46.4466 29.5288,47.7834 10.934,48.2323 4.58677,48.587 0.296648,48.9992 -4.00063,49.299 -6.95733,"
                              "49.6132 -9.93119,50.1835 -15.8585,50.3639 -17.8241,50.6215 -19.6471,51.0408 -21.4515,51.7068 -23.3609,52.3326 -24.729,53.0701 -26.047,"
                              "53.9117 -27.3065,54.8496 -28.4987,55.8759 -29.6152,56.9829 -30.6474,58.1629 -31.5867,59.408 -32.4246,60.9346 -33.2618,62.5325 -33.9455,"
                              "64.1878 -34.476,65.8864 -34.8539,67.3309 -35.051,68.7878 -35.1355,70.2399 -35.1053,72.7143 -34.7658,74.0086 -34.5489,75.3002 -38.7765)))", polygon3);

    MultiPolygonType unionResult2;
    boost::geometry::union_(unionResult1, polygon3, unionResult2); // throws

    std::cout << "Result: " << boost::geometry::wkt(unionResult2) << std::endl;

    return EXIT_SUCCESS;

}

Answer 1

I think whatever was the bug, was fixed.

Here's the result of the extended sample on my box with GCC 5.4 and Boost 1.62.0:

Result: MULTIPOLYGON(((75.3002 -38.7765,72.918 -39.2671,70.5558 -39.533,68.2188
-39.5707,65.912 -39.3766,63.6406 -38.9473,61.4095 -38.2792,59.2239
-37.3688,57.089 -36.2125,55.5367 -35.1756,54.0721 -34.0078,52.7036
-32.7209,51.4398 -31.3265,50.2893 -29.8362,49.2607 -28.2616,48.3624
-26.6144,47.6032 -24.9062,47.158 -23.7161,46.803 -22.6211,46.3142
-20.5856,46.0379 -18.5376,45.9485 -17.427,45.8751 -16.215,45.5849
-10.2471,45.4467 -7.25562,45.2962 -4.26324,45.0667 0.0336182,44.8334
4.32845,44.3738 10.6768,42.9884 29.2827,42.3055 38.5987,41.6559 47.9296,45.1275
48.1622,45.7738 38.8379,46.4466 29.5288,47.7834 10.934,48.2323 4.58677,48.587
0.296648,48.9992 -4.00063,49.299 -6.95733,49.6132 -9.93119,50.1835
-15.8585,50.3639 -17.8241,50.6215 -19.6471,51.0408 -21.4515,51.7068
-23.3609,52.3326 -24.729,53.0701 -26.047,53.9117 -27.3065,54.8496
-28.4987,55.8759 -29.6152,56.9829 -30.6474,58.1629 -31.5867,59.408
-32.4246,60.9346 -33.2618,62.5325 -33.9455,64.1878 -34.476,65.8864
-34.8539,67.3309 -35.051,68.7878 -35.1355,70.2399 -35.1053,72.7143
-34.7658,74.0086 -34.5489,75.3002 -38.7765),(-62.2843 -29.2642,-59.997
-28.7264,-58.2365 -28.1769,-62.2843 -29.2642)))