1

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?

David R.
  • 855
  • 8
  • 17

1 Answers1

2

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))
sehe
  • 374,641
  • 47
  • 450
  • 633
  • I do not understand why is there an intersection in the first place? – David R. Jun 15 '17 at 23:35
  • I dunno. Have you tried visualizing? Are the inner rings correctly oriented (reverse direction from the outer rings)? – sehe Jun 15 '17 at 23:36
  • Yes. The polygon is definitely valid and the points do not overlap. They are very close. Is there a resolution issue somewhere? – David R. Jun 15 '17 at 23:38
  • Could be. Do you _have_ the visualisation to share? I'm not big on reading WKT – sehe Jun 15 '17 at 23:39
  • Added a visualization: https://i.stack.imgur.com/5yguN.png I'd say the two inner rings might be too close, indeed. Here's with the stroke-width to 0: https://i.stack.imgur.com/CEmMA.png – sehe Jun 15 '17 at 23:43
  • The rectangle ends at X coordinate `9064.333999999999` and the rectilinear polygon begins at X coordinate `9064.3340000000007`. While close, they do not overlap. This WKT was produced by Shapely and Python subtracting the union of the two inner polygons (when they were polygons and not holes) from the larger rectangle. I would have thought it is a representation issue, but then LLDB prints the points as truly having different representation. Thus my question if there is any resolution that's being enforced. – David R. Jun 15 '17 at 23:48
  • It's clearly not being enforced. However, floating point representations are inexact, and hence the calculations are lossy with a certain margin of error. – sehe Jun 15 '17 at 23:52
  • A fix: **[Live On Coliru](http://coliru.stacked-crooked.com/a/5f02799549ede8bd)** (see answer update) – sehe Jun 15 '17 at 23:59
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/146805/discussion-between-sehe-and-david-r). – sehe Jun 16 '17 at 00:21