boost::within 的问题：测试一个多边形是否完全在另一个多边形内

Question

我正在做一个需要一些简单几何函数的项目。基本上，目标是测试一个多边形是否完全在另一个多边形内。按照boost::within的说法，应该可以搞定的。

代码如下：

#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;
}

这是可视化：

显然，poly2(red) 不应该在 poly1(green) 之内/覆盖。但是函数内部/覆盖都返回true。我知道有时浮点问题可能很烦人，但这种特殊情况太简单了，结果不应该受到浮点的影响。我认为 boost 足够强大，这是一个错误，还是我没有以正确的方式使用该函数。

任何帮助将不胜感激，在此先感谢

Answer 1

几何图形确实重叠：

缩放：

确实，提高精度会改变判断，即使仍然会有（这次是很小的）重叠：

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");
}

有趣的是，事实证明降低浮点精度是让样本按预期运行的唯一快速方法：

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

这是一个带有生成的 svg 的 zip：http://stackoverflow-sehe.s3.amazonaws.com/e7953654-8618-4002-bb00-90316ee11fb3/svgs.zip

这是浮动案例的缩放渲染（为了更好的视图，不透明度和线宽减小了）：

如您所见，它似乎仍然是重叠的，但它在算法上显然只是在边缘情况的安全方面。