Geometry, intersection

Question

I have a problem. I would like to intersect a quad with a quad.

int main(){
  typedef boost::geometry::model::point_xy<double> TBoostPoint;
  typedef boost::geometry::model::polygon<TBoostPoint> TBoostPoly;
  TBoostPoint point;
  TBoostPoly firstPoly, secondPoly;
  boost::geometry::read_wkt("POLYGON(
                (1.504477611940313, 3.761194029850755), 
                (1.504477611940305, 3.573134328358203),
                (1.316417910447765, 3.573134328358206),
                (1.316417910447769, 3.761194029850752))", firstPoly);
  boost::geometry::read_wkt("POLYGON(
                (1.504477611940313, 3.761194029850755), 
                (1.504477611940305, 3.573134328358203),
                (1.316417910447765, 3.573134328358206),
                (1.316417910447751, 3.761194029850769))", secondPoly);
 std::vector<TBoostPoly> outPoly;
 boost::geometry::intersection(firstPoly,secondPoly,outPoly);
}

outPoly - is empty, but it not so.

Answer 1

There were 2 main issues.

The output is undefined because the input is invalid.

1. The input WKT specifies a lot of invalid inner rings (consisting of single points), instead of what you expected, a single outer ring of 5 points (excl. closing point). Fix it:

   bg::read_wkt("POLYGON(( 1.504477611940313 3.761194029850755, 1.504477611940305 3.573134328358203, 1.316417910447765 3.573134328358206, 1.316417910447769 3.761194029850752))", first);
   bg::read_wkt("POLYGON(( 1.504477611940313 3.761194029850755, 1.504477611940305 3.573134328358203, 1.316417910447765 3.573134328358206, 1.316417910447751 3.761194029850769))", second); 
   

2. Boost Geometry assumes throughout that you never make errors against the documented preconditions. If you read in the polygon concept page and preconditions for intersection you'll see thew full list¹.

   If you don't, you get no friendly errors, just silent failure, corruption or just wrong answers. Yeah. That's bad.

   What's worse, BGeo didn't even have a is_valid facility to test the bulk of requirements until Boost 1_57 (IIRC). The good news is, if you upgrade to this version or later your life will be much simpler.

In this case you would have learned that the polygons weren't properly closed:

Live On Coliru

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/algorithms/intersection.hpp>
#include <boost/geometry/algorithms/correct.hpp>
#include <boost/geometry/algorithms/is_valid.hpp>

namespace bg = boost::geometry;

int main(){
    typedef bg::model::d2::point_xy<double> TPoint;
    typedef bg::model::polygon<TPoint>      TPoly;
    TPoly first, second;

    bg::read_wkt("POLYGON(( 1.504477611940313 3.761194029850755, 1.504477611940305 3.573134328358203, 1.316417910447765 3.573134328358206, 1.316417910447769 3.761194029850752))", first);
    bg::read_wkt("POLYGON(( 1.504477611940313 3.761194029850755, 1.504477611940305 3.573134328358203, 1.316417910447765 3.573134328358206, 1.316417910447751 3.761194029850769))", second); 

    std::string reason;
    // polys not closed!
    if (!bg::is_valid(first, reason))  std::cout << "First polygon not valid: "  << reason << "\n";
    if (!bg::is_valid(second, reason)) std::cout << "Second polygon not valid: " << reason << "\n";

    bg::correct(first);
    bg::correct(second);

    // no more output!
    if (!bg::is_valid(first, reason))  std::cout << "First polygon not valid: "  << reason << "\n";
    if (!bg::is_valid(second, reason)) std::cout << "Second polygon not valid: " << reason << "\n";

    std::vector<TPoly> out;
    bg::intersection(first, second, out);

    for (auto& g : out)
        std::cout << "\nresult: " << bg::wkt(g) << "\n";
}

Prints:

First polygon not valid: Geometry is defined as closed but is open
Second polygon not valid: Geometry is defined as closed but is open

Oops. The geos weren't closed! correct(poly) fixes this for us on auto-pilot:

result: POLYGON((1.50448 3.57313,1.31642 3.57313,1.31642 3.76119,1.50448 3.76119,1.50448 3.57313))

---

¹ outer ring must be counter clockwise, inner cw, polys must be closed... stuff like that.