Calculate polygon intersection area with GPC or Clipper

Question

I have a 2D particle system, where the particles are represented as ellipses. I need to calculate ellipse - ellipse overlap areas, but this is a hard analytical problem Ellipse-Ellipse Overlap. I now represent my ellipses as 20-gons, so that they are "polygonized" and I am using Boost.Geometry to make the necessary calculations.

However, very often, I get an exception from Boost.Geometry : boost.geometry overlay invalid input exception. I searched and found that this is a known bug of boost.Geometry and there is no fix as of version 1.53. Even the documentation on upcoming v1.54 does not say anything about fixing this problem.

I have stumbled upon Clipper and GPC - General Polygon Clipper Library. They seem to do what I want, but they only output boolean results. Does anyone know if there is a way to output the area of the calculated intersections with these libraries? I guess that since the intersection is stored in memory as some kind of polygon, I could use triangulation or some other method to calculate the area. Any pointers would be really appreciated!

The Boost overlay exception is consistent in MSVC 2010 and 2012 under Win7 x64, MinGW and Qt 4.8.1 under Linux Mint 14.

Answer 1

It seems you're wrong about Clipper and GPC "they only output boolean results". Both libraries calculate intersection polygon - for example, look at the code snippet with picture on Clipper page.

Answer 2

Can you submit a ticket into the Boost ticket system? Including a sample of a geometry pair giving invalid input exceptions? The error can have several causes, either it is indeed invalid, or it is a bug in the library. Some of these issues are fixed per 1.54.

Answer 3

In Angus' benchmark of Clipper with the other libraries, he has a method for polygon Area calculation. I used that and modified his ellipse method. The result is the following:

void Ellipse2Poly(double theta, double A1, double B1, double H1, double K1, Poly& p)
{ 
    const double pi = 3.1415926535898, tolerance = 0.125;
    const int n = 30;
    double step = pi/(double)n;
    double a = A1; // Long semi-axis length
    double b = B1; // short semi-axis length
    double xc = H1; // current X position
    double yc = K1; // current Y position
    double sintheta = sin(theta);
    double costheta = cos(theta);
    double t = 0;
    p.resize(2*n+1);
    for (int i = 0; i < 2*n+1; i++)
    {
        p[i].x = xc + a*cos(t)*costheta - b*sin(t)*sintheta;
        p[i].y = yc + a*cos(t)*sintheta + b*sin(t)*costheta;
        t += step;
    }
}

And the area calculation is:

double OverlapArea(Poly poly1, Poly poly2)
{
    Poly clip;
    Polys polys_poly1, polys_poly2, polys_clip;
    polys_poly1.resize(1);
    polys_poly2.resize(1);
    polys_clip.resize(1);
    polys_poly1[0] = poly1;
    polys_poly2[0] = poly2;
    polys_clip[0] = clip;
    Polygons clipper_polys_poly1, clipper_polys_poly2, clipper_polys_clip;
    LoadClipper(clipper_polys_poly1,polys_poly1);
    LoadClipper(clipper_polys_poly2,polys_poly2);
    ClipType op = ctIntersection;
    Clipper cp;
    cp.AddPolygons(clipper_polys_poly1,ptSubject);
    cp.AddPolygons(clipper_polys_poly2,ptClip);
    cp.Execute(op,clipper_polys_clip,pftEvenOdd,pftEvenOdd);
    if(clipper_polys_clip.size() != 0)
    {
        UnloadClipper(polys_clip,clipper_polys_clip);
        return Area(polys_clip[0]);
    }
    else
        return 0.0;
}

It is slower than Boost.Geometry but very stable. I run it for 2*10^5 time steps overnight and it didn't crash. It took 6h 12min, while Boost would take about 4h, but I am happy.

EDIT:

I had rewritten some of the other functions, so the comparison was unfair to Clipper! For the same configuration and 100K time steps, Clipper completed the operation in 2h 36min which I find very nice! Boost would complete the same operation in approximately 2h 15min, so I find them very comparable for these long runs!