boost read_wkt 产生了无效的多边形

boost read_wkt produced invalid polygon

在下面的代码中我使用read_wkt来初始化一个polygon。多边形两个洞

#include <iostream>
#define BOOST_GEOMETRY_TEST_DEBUG

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

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

typedef double base_type;
typedef bgm::d2::point_xy<base_type> point_type;
typedef bgm::polygon<point_type> polygon_type;
typedef bgm::multi_polygon<polygon_type> multipolygon_type;

int main() {
    polygon_type in;
    bg::read_wkt("POLYGON ((0 0, 0 15998.49, 12798.76 15998.49, 12798.76 0, 0 0), "
                     "(3921.294 177.8112, 9064.333999999999 177.8112, 9064.333999999999 2951.2112, 3921.294 2951.2112, 3921.294 177.8112), "
                     "(9064.334000000001 177.8112, 12765.034 177.8112, 12765.034 5192.0872, 12743.139 5192.0872, 12743.139 6685.701000000001, 11439.19 6685.701000000001, 11439.19 5192.0872, 11438.834 5192.0872, 11438.834 2951.2112, 9064.334000000001 2951.2112, 9064.334000000001 177.8112), )", in);
    std::cout << (bg::is_valid(in)?"valid":"invalid") << std::endl;
    return 0;
}

read_wkt 之后,多边形被报告为无效。

checking exterior ring...
checking interior rings...
computing and analyzing turns...
turns: [t,x/i {-1, -1} {0, 1} {0, 9} (9064.33, 177.811)]

invalid

我在调试器里查了一下,在内部表示中,第一个孔的点9064.333999999999 177.81129064.334000000001 177.8112确实不一样。

(const boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true> >) [=13=] = {
  std::__1::vector<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true>, std::__1::allocator<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true> > > = size=1 {
    [0] = {
      m_outer = {
        std::__1::vector<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, std::__1::allocator<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian> > > = size=5 {
          [0] = {
            boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
              m_values = ([0] = 0, [1] = 0)
            }
          }
          [1] = {
            boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
              m_values = ([0] = 0, [1] = 15998.49)
            }
          }
          [2] = {
            boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
              m_values = ([0] = 12798.76, [1] = 15998.49)
        }
      }
      [3] = {
        boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
          m_values = ([0] = 12798.76, [1] = 0)
        }
      }
      [4] = {
        boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
          m_values = ([0] = 0, [1] = 0)
        }
      }
    }
  }
  m_inners = size=2 {
    [0] = {
      std::__1::vector<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, std::__1::allocator<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian> > > = size=5 {
        [0] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 3921.2939999999999, [1] = 177.81120000000001)
          }
        }
        [1] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 9064.3339999999989, [1] = 177.81120000000001)
          }
        }
        [2] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 9064.3339999999989, [1] = 2951.2112000000002)
          }
        }
        [3] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 3921.2939999999999, [1] = 2951.2112000000002)
          }
        }
        [4] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 3921.2939999999999, [1] = 177.81120000000001)
          }
        }
      }
    }
    [1] = {
      std::__1::vector<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, std::__1::allocator<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian> > > = size=11 {
        [0] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 9064.3340000000007, [1] = 177.81120000000001)
          }
        }
        [1] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 12765.034, [1] = 177.81120000000001)
          }
        }
        [2] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 12765.034, [1] = 5192.0871999999999)
          }
        }
        [3] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 12743.138999999999, [1] = 5192.0871999999999)
          }
        }
        [4] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 12743.138999999999, [1] = 6685.7010000000009)
          }
        }
        [5] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 11439.190000000001, [1] = 6685.7010000000009)
          }
        }
        [6] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 11439.190000000001, [1] = 5192.0871999999999)
          }
        }
        [7] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 11438.834000000001, [1] = 5192.0871999999999)
          }
        }
        [8] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 11438.834000000001, [1] = 2951.2112000000002)
          }
        }
        [9] = {
          boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
            m_values = ([0] = 9064.3340000000007, [1] = 2951.2112000000002)
          }
        }
            [10] = {
              boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> = {
                m_values = ([0] = 9064.3340000000007, [1] = 177.81120000000001)
              }
            }
          }
        }
      }
    }
  }
}

有什么建议吗?

打印原因有帮助:

std::string reason;
std::cout << (bg::is_valid(in, reason)?"valid ":"invalid ") << reason << std::endl;

版画

invalid Geometry has invalid self-intersections. A self-intersection point was found at (9064.33, 177.811); method: t; operations: x/i; segment IDs {source, multi, ring, segment}: {0, -1, 0, 0}/{0, -1, 1, 9}

现在更正:

bg::correct(in);
std::cout << std::fixed << std::setprecision(3) << bg::wkt(in) << "\n";

不删除自交集:

POLYGON((0.000 0.000,0.000 15998.490,12798.760 15998.490,12798.760 0.000,0.000 0.000),(3921.294 177.811,9064.334 177.811,9064.334 2951.211,3921.294 2951.211,3921.294 177.811),(9064.334 177.811,12765.034 177.811,12765.034 5192.087,12743.139 5192.087,12743.139 6685.701,11439.190 6685.701,11439.190 5192.087,11438.834 5192.087,11438.834 2951.211,9064.334 2951.211,9064.334 177.811))

奖金:

可视化:

将笔划宽度减小为 0:

结论

对于所选点类型的精度而言,内圈太近了。通过降低分辨率或引入更高精度的点类型来修复它:

typedef boost::multiprecision::cpp_dec_float_100 base_type;

这个有效:

Live On Coliru

#include <iostream>
#include <fstream>
//#define BOOST_GEOMETRY_TEST_DEBUG


#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/point_xy.hpp>

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

typedef boost::multiprecision::cpp_dec_float_100 base_type;
typedef bgm::d2::point_xy<base_type> point_type;
typedef bgm::polygon<point_type> polygon_type;
typedef bgm::multi_polygon<polygon_type> multipolygon_type;

int main() {
    polygon_type in;
    bg::read_wkt("POLYGON ((0 0, 0 15998.49, 12798.76 15998.49, 12798.76 0, 0 0), "
                     "(3921.294 177.8112, 9064.333999999999 177.8112, 9064.333999999999 2951.2112, 3921.294 2951.2112, 3921.294 177.8112), "
                     "(9064.334000000001 177.8112, 12765.034 177.8112, 12765.034 5192.0872, 12743.139 5192.0872, 12743.139 6685.701000000001, 11439.19 6685.701000000001, 11439.19 5192.0872, 11438.834 5192.0872, 11438.834 2951.2112, 9064.334000000001 2951.2112, 9064.334000000001 177.8112), )", in);
;
    {
        std::ofstream svg("svg.svg");
        boost::geometry::svg_mapper<point_type> mapper(svg, 400, 400);
        mapper.add(in);

        mapper.map(in, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:0");
    }

    std::string reason;
    std::cout << (bg::is_valid(in, reason)?"valid ":"invalid ") << reason << std::endl;
    std::cout << bg::wkt(in) << "\n";
}

版画

valid Geometry is valid
POLYGON((0 0,0 15998.5,12798.8 15998.5,12798.8 0,0 0),(3921.29 177.811,9064.33 177.811,9064.33 2951.21,3921.29 2951.21,3921.29 177.811),(9064.33 177.811,12765 177.811,12765 5192.09,12743.1 5192.09,12743.1 6685.7,11439.2 6685.7,11439.2 5192.09,11438.8 5192.09,11438.8 2951.21,9064.33 2951.21,9064.33 177.811))