Issue with boost::within: test whether a polygon is completely inside the another polygon

user_12023517 0
Viewed 114

I am working on a project which requires some simple geometry functions. Basically, the objective is to test whether a polygon is completely inside another polygon. According to boost::within, it should be able to handle this.

Here is the code:

#include <iostream>
#include <list>

#include <string>

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


int main()
{
    typedef boost::geometry::model::d2::point_xy<double> point_type;
    typedef boost::geometry::model::polygon<point_type> polygon_type;


    polygon_type poly;

    boost::geometry::read_wkt(
            "POLYGON((76.00000000 430.00000000,63.00000000 430.00000000,63.00000000 441.00000000,64.00000000 441.00000000,64.00000000 448.00000000,65.00000000 448.00000000,65.00000000 462.00000000,64.00000000 464.90909091,64.00000000 466.00000000,46.66666667 518.00000000,47.00000000 518.00000000,47.00000000 519.00000000,48.00000000 519.00000000,48.00000000 520.00000000,49.00000000 520.00000000,49.00000000 521.00000000,50.00000000 521.00000000,50.00000000 522.00000000,51.00000000 522.00000000,51.00000000 524.00000000,52.00000000 524.00000000,52.00000000 524.80000000,56.00000000 509.00000000,56.75641026 509.00000000,57.00000000 508.00000000,58.66666667 508.00000000,74.00000000 439.00000000,74.50000000 439.00000000,75.00000000 436.00000000,76.00000000 436.00000000,76.00000000 430.00000000))"
            ,poly);
    polygon_type poly2;

    boost::geometry::read_wkt(
            "POLYGON((72.00000000 448.00000000,72.00000000 456.00000000,80.00000000 456.00000000,80.00000000 448.00000000,72.00000000 448.00000000))",
            poly2);


    std::cout << "within: " << (boost::geometry::within(poly2, poly) ? "yes" : "no") << std::endl;
    std::cout << "covered by: " << (boost::geometry::covered_by(poly2, poly) ? "yes" : "no") << std::endl;

    return 0;
}

Here is the visualization:

vis

Apparently, the poly2(red) should not be within/covered by the poly1(green). But both within/covered by function return true. I understand sometimes floating-point issues could be annoying, but this special case is way too simple and the result should not be affected by the floating-point. I thought boost is robust enough, is this a bug, or I didn't use the function in the right way.

Any helps would be appreciated, thanks in advance

boost-geometry polygon geometry boost
1 Answers

The geometries do overlap:

enter image description here

Zoom:

enter image description here

Indeed upping the precision will change the verdict, even though there will still be a (tiny, this time) overlap:

Live On Compiler Explorer

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <iostream>
#include <fstream>

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

template <typename coord_type>
void run_tests(std::string caption) {
    using point_type   = bgm::d2::point_xy<coord_type>;
    using polygon_type = bgm::polygon<point_type>;

    polygon_type poly1, poly2;

    bg::read_wkt(
        "POLYGON((76.00000000 430.00000000,63.00000000 "
        "430.00000000,63.00000000 441.00000000,64.00000000 "
        "441.00000000,64.00000000 448.00000000,65.00000000 "
        "448.00000000,65.00000000 462.00000000,64.00000000 "
        "464.90909091,64.00000000 466.00000000,46.66666667 "
        "518.00000000,47.00000000 518.00000000,47.00000000 "
        "519.00000000,48.00000000 519.00000000,48.00000000 "
        "520.00000000,49.00000000 520.00000000,49.00000000 "
        "521.00000000,50.00000000 521.00000000,50.00000000 "
        "522.00000000,51.00000000 522.00000000,51.00000000 "
        "524.00000000,52.00000000 524.00000000,52.00000000 "
        "524.80000000,56.00000000 509.00000000,56.75641026 "
        "509.00000000,57.00000000 508.00000000,58.66666667 "
        "508.00000000,74.00000000 439.00000000,74.50000000 "
        "439.00000000,75.00000000 436.00000000,76.00000000 "
        "436.00000000,76.00000000 430.00000000))",
        poly1);

    bg::read_wkt(
        "POLYGON((72.00000000 448.00000000,72.00000000 "
        "456.00000000,80.00000000 456.00000000,80.00000000 "
        "448.00000000,72.00000000 448.00000000))",
        poly2);

    std::string reason;
    if (!bg::is_valid(poly1, reason)) {
        std::cout << caption << ", Correcting poly1: " << reason << "\n";
        bg::correct(poly1);
    }
    if (!bg::is_valid(poly2, reason)) {
        std::cout << caption << ", Correcting poly2: " << reason << "\n";
        bg::correct(poly2);
    }

    std::cout << std::boolalpha;
    std::cout << caption << ", within: "     << bg::within(poly2, poly1)     << std::endl;
    std::cout << caption << ", covered by: " << bg::covered_by(poly2, poly1) << std::endl;

    bgm::multi_polygon<polygon_type> out;
    bg::intersection(poly1, poly2, out);
    std::cout << caption << ", intersection " << bg::wkt(out) << " with area "
        << bg::area(out) << "\n";

    {
        // Declare a stream and an SVG mapper
        std::ofstream svg("output-" + caption + ".svg");
        bg::svg_mapper<point_type> mapper(svg, 400, 400);

        // Add geometries such that all these geometries fit on the map
        mapper.add(poly1);
        mapper.add(poly2);
        mapper.add(out);

        // Draw the geometries on the SVG map, using a specific SVG style
        mapper.map(poly1, "fill-opacity:0.3;fill:rgb(51,0,0);stroke:rgb(51,0,0);stroke-width:1");
        mapper.map(poly2, "fill-opacity:0.3;fill:rgb(0,51,0);stroke:rgb(0,51,0);stroke-width:1");
        //mapper.map(c, "fill-opacity:0.3;fill:rgb(0,0,51);stroke:rgb(0,0,51);stroke-width:2");
    }
}

int main()
{
    run_tests<float>("float");
    run_tests<double>("double");
    run_tests<long double>("long double");
    run_tests<bmp::number<bmp::backends::cpp_dec_float<18>, bmp::et_off>>( "dec18");
    run_tests<bmp::number<bmp::backends::cpp_dec_float<50>, bmp::et_off>>( "dec50");
    run_tests<bmp::number<bmp::backends::cpp_dec_float<100>, bmp::et_off>>( "dec100");
    run_tests<bmp::number<bmp::backends::cpp_bin_float<100>, bmp::et_off>>( "bin100");
}

Interestingly, it turns out that reducing the floating point precision is the only quick way yo get your sample to behave as you expect:

float, within: false
float, covered by: false
float, intersection MULTIPOLYGON() with area 0
double, within: true
double, covered by: true
double, intersection MULTIPOLYGON(((72 448,72 456,80 456,80 448,72 448))) with area 64
long double, within: true
long double, covered by: true
long double, intersection MULTIPOLYGON(((72 448,72 456,80 456,80 448,72 448))) with area 64
dec18, within: false
dec18, covered by: false
dec18, intersection MULTIPOLYGON(((72 448,72 448,72 448,72 448))) with area 4.25331e-19
dec50, within: false
dec50, covered by: false
dec50, intersection MULTIPOLYGON(((72 448,72 448,72 448,72 448))) with area 4.25331e-19
dec100, within: false
dec100, covered by: false
dec100, intersection MULTIPOLYGON(((72 448,72 448,72 448,72 448))) with area 4.25331e-19
bin100, within: false
bin100, covered by: false
bin100, intersection MULTIPOLYGON(((72 448,72 448,72 448,72 448))) with area 4.25331e-19

Here's a zip with the generated svgs: http://stackoverflow-sehe.s3.amazonaws.com/e7953654-8618-4002-bb00-90316ee11fb3/svgs.zip

And here's a zoomed render of the float case (with opacities and line widths reduced for better view):

enter image description here

As you can see, still seems to be overlapping, but it algorithmically apparently falls just on the safe side of the edge-case.

user_853710
Related