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

【问题标题】：Calculate polygon intersection area with GPC or Clipper使用 GPC 或 Clipper 计算多边形相交面积
【发布时间】：2013-05-15 13:20:54
【问题描述】：

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

但是，我经常收到来自Boost.Geometry 的异常：boost.geometry overlay invalid input exception。我搜了一下发现这是boost.Geometry的一个已知bug，从1.53版开始没有修复。甚至即将发布的 v1.54 的文档也没有说明解决此问题的任何内容。

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

Boost overlay 异常在 Win7 x64 下的 MSVC 2010 和 2012、Linux Mint 14 下的 MinGW 和 Qt 4.8.1 中是一致的。

【问题讨论】：

标签： c++ geometry polygons boost-geometry

【解决方案1】：

看来您对 Clipper 和 GPC 的看法是错误的“它们只输出布尔结果”。两个库都计算交点多边形 - 例如，在 Clipper 页面上查看带有图片的代码 sn-p。

【讨论】：

实际上，您说的有一部分是正确的。在代码 sn-p 中，最后可以看到对于 Clipper c 对象，Angus 使用了这个：c.Execute(ctIntersection, solution, pftNonZero, pftNonZero);，这是一个布尔表达式。但是在他的基准测试中，他正在计算面积，我找到了一种根据我的需要修改他的代码的方法。我会将您的答案标记为正确的答案，因为它让我思考并让我走上了正确的轨道！谢谢！

【解决方案2】：

您可以向 Boost 票证系统提交票证吗？包括给出无效输入异常的几何对样本？该错误可能有多种原因，要么确实无效，要么是库中的错误。其中一些问题已在 1.54 中得到修复。

【讨论】：

感谢您的回复！我会提交一张票。很高兴听到这可能会在 1.54 中修复！

【解决方案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 在 2h 36min 内完成了操作，我觉得非常好！ Boost 将在大约 2 小时 15 分钟内完成相同的操作，因此我发现它们对于这些长时间运行非常具有可比性！

【讨论】：

您提到的论文的作者发布了他们算法的实现，并出于验证目的发布了聚合方法（您刚刚使用的代码）。分析方法比数值方法快约 10 倍。有什么理由让您更喜欢多边形方法而不是解析解决方案？