Segmentation fault in CGAL when intersecting polygons

Question

I'm writing a code where I need to calculate the overlapping area of two polygons very frequently (one always is a square, the other one is a simple but in general non-convex polygon). However, using CGAL for this, I occasionally encounter segmentation faults. The following code provides a minimal example. Note that as soon as one of the point coordinates is moved by 0.001, everything works fine.

#include <CGAL/Cartesian.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Polygon_with_holes_2.h>
#include <CGAL/Boolean_set_operations_2.h>
#include <CGAL/Polygon_2_algorithms.h>

typedef CGAL::Cartesian<double> K;
typedef K::Point_2 Point;
typedef CGAL::Polygon_2<K> Polygon;
typedef CGAL::Polygon_with_holes_2<K> PolygonH;

int main()
{    
     // Rectangle
     Point pointsA[] = { Point(0.4,0.35), Point(0.4,0.4), Point(0.45,0.4), Point(0.5,0.35) };

     // Triangle
     Point pointsB[] = { Point(0.428,0.412), Point(0.387,0.389), Point(0.359,0.393) };

     // Convert to CGAL polygon
     Polygon polyA(pointsA, pointsA + 4);
     Polygon polyB(pointsB, pointsB + 3);

     // Intersection
     std::list<PolygonH> polyAB;
     CGAL::intersection(polyA, polyB, std::back_inserter(polyAB));
}

Answer 1

If the first polygon is always convex (as most squares are) and you only care for the area of overlap, you can roll your own solution easily.

Consider the support line of one of the sides of the first polygon. You will implement clipping against a half plane, let us assume bounded by x=0, WLOG.

Create an empty output polygon. Then try all edges of the second polygon in turn. Assign every vertex the boolean condition x>=0. If the condition is different for both endpoints, append the intersection of the side with the axis x=0 to the new polygon; then if the condition is true for the final endpoint of the side, append it as well.

In the end, you'll get a new polygon which is the clipped version of the original. In case the clipped polygon should be disconnected, the algorithm will produce a connected polygon with extra connecting sides along x=0. But as regards the computation of the area, this doesn't matter at all.

Repeat with the four clipping sides and compute the area using the shoelace formula.