mercredi 24 août 2016

Precision problems with `boost::geometry::difference`

Most of the time when I use boost::geometry::difference it does what I'd expect. However, when I have two polygons whose edges are almost coincident, I get some weird behavior. Consider this example:

#include <iostream>
#include <fstream>
#include <list>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/multi_polygon.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 printValidity(const TPolygon& polygons)
{
    boost::geometry::validity_failure_type failure;
    bool valid = boost::geometry::is_valid(polygons, failure);

    if(!valid) {
        std::cout << "not valid: " << failure << std::endl;
    }
    else {
        std::cout << "valid." << std::endl;
    }
}

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(size_t polygonID = 0; polygonID < polygons.size(); ++polygonID) {
        mapper.add(polygons[polygonID]);
        int redValue = 50 * polygonID; // NOTE: This will break with more than 5 polygons
        mapper.map(polygons[polygonID], "fill:rgb(" + std::to_string(redValue) + ",128,0);stroke:rgb(0,0,100);stroke-width:1");
    }
}

int main()
{
    MultiPolygonType polygon1;
    boost::geometry::read_wkt("MULTIPOLYGON(((-23.8915 -12.2238,-23.7604 -10.2739,-23.1774 -7.83411,-22.196 -5.52561,-20.8436 -3.41293,-19.7976 -2.26009,-19.8108 -2.00306,24.8732 2.51519,26.0802 -9.42191,-23.6662 -14.452,-23.8915 -12.2238)))", polygon1);
    WritePolygonsToSVG(polygon1, "polygon1.svg");

    MultiPolygonType polygon2;
    boost::geometry::read_wkt("MULTIPOLYGON(((-8.85138 -12.954,-8.89712 -12.7115,-8.9307 -12.5279,-3.35773 -12.078,-3.42937 -11.5832,-3.50013 -11.0882,-3.57007 -10.5932,-3.63925 -10.098,-3.70775 -9.60273,-3.77561 -9.10737,-3.84289 -8.61192,-3.90967 -8.1164,-3.976 -7.62081,-4.04194 -7.12517,-4.10756 -6.62948,-4.17291 -6.13375,-4.23805 -5.63799,-4.30305 -5.14222,-4.36797 -4.64643,-4.43287 -4.15064,-4.49781 -3.65485,-4.56285 -3.15909,-4.62806 -2.66331,-4.69349 -2.16756,-4.75921 -1.67185,-4.82528 -1.17619,-4.89175 -0.680597,-4.91655 -0.497015,17.8166 1.80166,17.8143 1.61653,17.8078 1.11656,17.8009 0.61658,17.7937 0.116605,17.7863 -0.383369,17.7786 -0.883343,17.7707 -1.3833,17.7627 -1.88324,17.7545 -2.38318,17.7463 -2.88312,17.7381 -3.38307,17.7299 -3.88301,17.7218 -4.38295,17.7139 -4.8829,17.706 -5.38285,17.6984 -5.88279,17.6911 -6.38274,17.684 -6.88269,17.6772 -7.38265,17.6709 -7.8826,17.6649 -8.38256,17.6595 -8.88251,17.6545 -9.38247,17.6501 -9.88244,17.6471 -10.2746,-8.85138 -12.954)))", polygon2);
    WritePolygonsToSVG(polygon2, "polygon2.svg");

    MultiPolygonType differencePolygon;

    boost::geometry::difference(polygon1, polygon2, differencePolygon);
    WritePolygonsToSVG(differencePolygon, "difference.svg");

    printValidity(differencePolygon);

    MultiPolygonType realDifference;
    boost::geometry::read_wkt("MULTIPOLYGON(((-23.8915 -12.2238,-23.7604 -10.2739,-23.1774 -7.83411,-22.196 -5.52561,-20.8436 -3.41293,-19.7976 -2.26009,-19.8108 -2.00306,24.8732 2.51519,26.0802 -9.42191,-23.6662 -14.452,-23.8915 -12.2238),(-8.85138 -12.954,17.6471 -10.2746,17.6501 -9.88244,17.6545 -9.38247,17.6595 -8.88251,17.6649 -8.38256,17.6709 -7.8826,17.6772 -7.38265,17.684 -6.88269,17.6911 -6.38274,17.6984 -5.88279,17.706 -5.38285,17.7139 -4.8829,17.7218 -4.38295,17.7299 -3.88301,17.7381 -3.38307,17.7463 -2.88312,17.7545 -2.38318,17.7627 -1.88324,17.7707 -1.3833,17.7786 -0.883343,17.7863 -0.383369,17.7937 0.116605,17.8009 0.61658,17.8078 1.11656,17.8143 1.61653,17.8166 1.80166,-4.91655 -0.497015,-4.89175 -0.680597,-4.82528 -1.17619,-4.75921 -1.67185,-4.69349 -2.16756,-4.62806 -2.66331,-4.56285 -3.15909,-4.49781 -3.65485,-4.43287 -4.15064,-4.36797 -4.64643,-4.30305 -5.14222,-4.23805 -5.63799,-4.17291 -6.13375,-4.10756 -6.62948,-4.04194 -7.12517,-3.976 -7.62081,-3.90967 -8.1164,-3.84289 -8.61192,-3.77561 -9.10737,-3.70775 -9.60273,-3.63925 -10.098,-3.57007 -10.5932,-3.50013 -11.0882,-3.42937 -11.5832,-3.35773 -12.078,-8.9307 -12.5279,-8.89712 -12.7115,-8.85138 -12.954)))", realDifference);
    WritePolygonsToSVG(realDifference, "realDifference.svg");
    printValidity(realDifference);

    return 0;
}

The 'differencePolygon' computed using boost::geometry::difference in this code is "valid", though it has an infinitesimally thin region: enter image description here

The 'realDifference' polygon is what I output from my "real code". Apparently there was some precision loss with the wkt writer in trying to construct this example, but you can see this polygon has two infinitesimally small regions, and actually returns not valid (21: failure_self_intersections): enter image description here

Is there some kind of tolerance parameter than can be set to avoid these infinitesimally small regions during the difference operation (which apparently sometimes break the polygon)? Or is there a function to "correct" the polygon by removing these self intersections? (I know the correct() function only fixes ring ordering).

Aucun commentaire:

Enregistrer un commentaire