用GPC或Clipper计算多边形相交面积

Question

我有一个2D粒子系统，其中粒子表示为椭圆。我需要计算椭圆-椭圆重叠面积，但这是一个很难分析的问题Ellipse-Ellipse Overlap。我现在将椭圆表示为20个多边形，因此它们是“多边形化的”，并且我使用Boost.Geometry进行必要的计算。

但是，我经常从Boost.Geometry得到一个异常：boost.geometry overlay invalid input exception。我搜索了一下，发现这是boost.Geometry的一个已知错误，1.53版还没有修复。即使是即将发布的v1.54的文档也没有提到解决这个问题的任何内容。

我偶然发现了Clipper和GPC - General Polygon Clipper Library。它们似乎做了我想要做的事情，但它们只输出布尔结果。有没有人知道有没有办法输出这些库计算出的交叉点的面积？我猜，由于交点作为某种多边形存储在内存中，我可以使用三角剖分或其他方法来计算面积。如果有任何建议，我们将不胜感激！

在Win7 x64、MinGW和Linux Mint14下的Qt 4.8.1中，Boost overlay异常在MSVC2010和2012中是一致的。

Answer 1

看来你对Clipper和GPC的理解是错误的，“他们只输出布尔结果”。这两个库都计算相交多边形-例如，在Clipper页面上查看带有图片的代码片段。

Answer 2

你能向Boost工单系统提交工单吗？包括给出无效输入异常的几何图形对的样本？错误可能有几个原因，要么确实无效，要么是库中的错误。其中一些问题已在1.54版本中修复。

Answer 3

在Angus的Clipper基准测试和其他库中，他有一个计算多边形面积的方法。我用它修改了他的椭圆方法。结果如下：

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

面积的计算是：

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

它比Boost.Geometry慢，但非常稳定。我连夜运行了2*10^5个时间步长，它没有崩溃。它需要6小时12分钟，而Boost大约需要4小时，但我很高兴。

编辑：

我已经重写了一些其他函数，所以这种比较对Clipper是不公平的！对于相同的配置和100K的时间步长，Clipper在2小时36分钟内完成了操作，我发现这非常好！Boost可以在大约2小时15分钟内完成相同的操作，所以我发现它们在这些长跑中非常有可比性！