David Doria David Doria - 20 days ago 7
C++ Question

boost::geometry::union_ produces self intersections

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:
enter image description here

Polygon2:
enter image description here

Both input polygons:
enter image description here

Union polygon (self intersections):
enter image description here

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

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)))
Comments