boost read_wkt produced invalid polygon

1
votes

In the following code I use read_wkt to initialize a polygon. The polygon two holes.

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

After read_wkt the polygon is reported as not valid. 

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

I checked in debugger that in the internal representation, the points 9064.333999999999 177.8112 from the first hole and 9064.334000000001 177.8112 are indeed different.

(const boost::geometry::model::multi_polygon<boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian>, true, true> >) $0 = {
  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)
              }
            }
          }
        }
      }
    }
  }
}

Any suggestions what is going on here?

1
c++boostboost-geometry

1 Answers

2
votes

Printing the reason helps:

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

Prints

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}

Now, correcting it:

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

Does not remove the self intersection:

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

Bonus:

Visualization: enter image description here

Reducing stroke-width to 0:enter image description here

Conclusion

The inner rings are too close for the accuracy of the point type chosen. Fix it, by either making the resolution less dense or introducing a higher accuracy point type:

typedef boost::multiprecision::cpp_dec_float_100 base_type;

This works:

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

Prints

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