CGAL 计算两个凸多边形相交的面积

发布于

我正在使用 CGAL 计算两个凸多边形相交的面积。在对 this 的接受答案中发布了执行此操作的简短演示代码。问题。但是，当我修改该代码以使用我感兴趣的多边形时，CGAL 从 CGAL::intersection() 例程的深处抛出运行时异常。

这是一个简短的示例代码，它是从上面链接的 SO 问题中复制粘贴的，除了它使用我自己的多边形并打印一些关于每个多边形的诊断信息以表明它们是凸面的并使用 CCW 绕组订单。

代码语言：javascript
复制
#include <iostream>
#include <CGAL/Simple_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::Simple_cartesian<double> K;

typedef K::Point_2 Point;

typedef CGAL::Polygon_2<K> Polygon_2;

typedef CGAL::Polygon_with_holes_2<K> Polygon_with_holes_2;

using std::cout; using std::endl;


int main(){

bool IsSimple, IsConvex, IsClockwise;

double Area;


// Define the polygons. P is long and skinny, Q is 20 equally
// spaced points around a circle of radius 0.025. The top-right
// corner of P protudes into Q.
int np=7;
int nq=20;
Point points[]={
    Point( 2.67905949256324049657e-02,-2.36307305336680428809e-02),
    Point( 2.19804738309115066386e-02, 1.25954939609798088895e-02),
    Point( 3.24142504626934169210e-03, 8.27109808759729503436e-03),
    Point(-7.07297020897392769712e-02,-1.31780902623212625713e-01),
    Point(-1.72945875771868884385e+00,-4.33625989924519217311e+00),
    Point(-5.64696462817379796206e+00,-2.91346747524371210147e+01),
    Point(-5.66016394049875870564e+00,-4.95259849820236013329e+01)};
Point points2[]={
    Point( 2.50000000000000013878e-02, 0.00000000000000000000e+00),
    Point( 2.37764129073788424429e-02, 7.72542485937368662852e-03),
    Point( 2.02254248593736890571e-02, 1.46946313073118266929e-02),
    Point( 1.46946313073118284276e-02, 2.02254248593736890571e-02),
    Point( 7.72542485937368836324e-03, 2.37764129073788424429e-02),
    Point( 1.53080849893419158600e-18, 2.49999999999999979183e-02),
    Point(-7.72542485937368489379e-03, 2.37764129073788424429e-02),
    Point(-1.46946313073118232234e-02, 2.02254248593736890571e-02),
    Point(-2.02254248593736890571e-02, 1.46946313073118284276e-02),
    Point(-2.37764129073788424429e-02, 7.72542485937368923060e-03),
    Point(-2.50000000000000013878e-02, 3.06161699786838317201e-18),
    Point(-2.37764129073788424429e-02,-7.72542485937368402643e-03),
    Point(-2.02254248593736925266e-02,-1.46946313073118232234e-02),
    Point(-1.46946313073118318970e-02,-2.02254248593736890571e-02),
    Point(-7.72542485937369096533e-03,-2.37764129073788389734e-02),
    Point(-4.59242549680257437283e-18,-2.50000000000000013878e-02),
    Point( 7.72542485937368229171e-03,-2.37764129073788424429e-02),
    Point( 1.46946313073118214887e-02,-2.02254248593736925266e-02),
    Point( 2.02254248593736855877e-02,-1.46946313073118318970e-02),
    Point( 2.37764129073788389734e-02,-7.72542485937369183269e-03)};
Polygon_2 polyP(points, points+np);
Polygon_2 polyQ(points2, points2+nq);

// Print some properties of polygon P. It is simple, convex, and CCW
// oriented
CGAL::set_pretty_mode(cout);
cout &lt;&lt; &#34;created the polygon P:&#34; &lt;&lt; endl;
cout &lt;&lt; polyP &lt;&lt; endl;
cout &lt;&lt; endl;
IsSimple    = polyP.is_simple();
IsConvex    = polyP.is_convex();
IsClockwise = (polyP.orientation() == CGAL::CLOCKWISE);
Area      = polyP.area();
cout &lt;&lt; &#34;polygon P is&#34;;
if (!IsSimple) cout &lt;&lt; &#34; not&#34;;
cout &lt;&lt; &#34; simple.&#34; &lt;&lt; endl;
cout &lt;&lt; &#34;polygon P is&#34;;
if (!IsConvex) cout &lt;&lt; &#34; not&#34;;
cout &lt;&lt; &#34; convex.&#34; &lt;&lt; endl;
cout &lt;&lt; &#34;polygon P is&#34;;
if (!IsClockwise) cout &lt;&lt; &#34; not&#34;;
cout &lt;&lt; &#34; clockwise oriented.&#34; &lt;&lt; endl;
cout &lt;&lt; &#34;the area of polygon P is &#34; &lt;&lt; Area &lt;&lt; endl;
cout &lt;&lt; endl;

// Print some properties of polygon Q. It is simple, convex, and CCW
// oriented
cout &lt;&lt; &#34;created the polygon Q:&#34; &lt;&lt; endl;
cout &lt;&lt; polyQ &lt;&lt; endl;
cout &lt;&lt; endl;
IsSimple    = polyQ.is_simple();
IsConvex    = polyQ.is_convex();
IsClockwise = (polyQ.orientation() == CGAL::CLOCKWISE);
Area      = polyQ.area();
cout &lt;&lt; &#34;polygon Q is&#34;;
if (!IsSimple) cout &lt;&lt; &#34; not&#34;;
cout &lt;&lt; &#34; simple.&#34; &lt;&lt; endl;
cout &lt;&lt; &#34;polygon Q is&#34;;
if (!IsConvex) cout &lt;&lt; &#34; not&#34;;
cout &lt;&lt; &#34; convex.&#34; &lt;&lt; endl;
cout &lt;&lt; &#34;polygon Q is&#34;;
if (!IsClockwise) cout &lt;&lt; &#34; not&#34;;
cout &lt;&lt; &#34; clockwise oriented.&#34; &lt;&lt; endl;
cout &lt;&lt; &#34;the area of polygon Q is &#34; &lt;&lt; Area &lt;&lt; endl;
cout &lt;&lt; endl;

// Compute the intersection
std::list&lt;Polygon_with_holes_2&gt; polyI;
CGAL::intersection(polyP, polyQ, std::back_inserter(polyI));

double totalArea = 0;
typedef std::list&lt;Polygon_with_holes_2&gt;::iterator LIT;
for(LIT lit = polyI.begin(); lit!=polyI.end(); lit++){
    totalArea+=lit-&gt;outer_boundary().area();
}
cout &lt;&lt; &#34;TotalArea::&#34; &lt;&lt; totalArea;

return 0;



}

这是测试程序的输出。

代码语言：javascript
复制
created the polygon P:

Polygon_2(

PointC2(0.0267906, -0.0236307)

PointC2(0.0219805, 0.0125955)

PointC2(0.00324143, 0.0082711)

PointC2(-0.0707297, -0.131781)

PointC2(-1.72946, -4.33626)

PointC2(-5.64696, -29.1347)

PointC2(-5.66016, -49.526)

)


polygon P is simple.

polygon P is convex.

polygon P is not clockwise oriented.

the area of polygon P is 71.1027


created the polygon Q:

Polygon_2(

PointC2(0.025, 0)

PointC2(0.0237764, 0.00772542)

PointC2(0.0202254, 0.0146946)

PointC2(0.0146946, 0.0202254)

PointC2(0.00772542, 0.0237764)

PointC2(1.53081e-18, 0.025)

PointC2(-0.00772542, 0.0237764)

PointC2(-0.0146946, 0.0202254)

PointC2(-0.0202254, 0.0146946)

PointC2(-0.0237764, 0.00772542)

PointC2(-0.025, 3.06162e-18)

PointC2(-0.0237764, -0.00772542)

PointC2(-0.0202254, -0.0146946)

PointC2(-0.0146946, -0.0202254)

PointC2(-0.00772542, -0.0237764)

PointC2(-4.59243e-18, -0.025)

PointC2(0.00772542, -0.0237764)

PointC2(0.0146946, -0.0202254)

PointC2(0.0202254, -0.0146946)

PointC2(0.0237764, -0.00772542)

)


polygon Q is simple.

polygon Q is convex.

polygon Q is not clockwise oriented.

the area of polygon Q is 0.00193136


terminate called after throwing an instance of 'CGAL::Precondition_exception'

what():  CGAL ERROR: precondition violation!

Expr: (m_traits.compare_y_at_x_2_object()(p, cv) == EQUAL) && compare_xy(cv.left(), p) == SMALLER && compare_xy(cv.right(), p) == LARGER

File: /usr/local/include/CGAL/Arr_segment_traits_2.h

Line: 706

Aborted (core dumped)

我不明白是什么导致了错误。任何帮助将不胜感激。

最佳答案

我可以重现此错误(在带有 clang++ 的 MacOS 上使用 CGAL 4.9)。据我了解，这种类型的未捕获异常不应该发生，换句话说，您发现了 CGAL 中的错误。因此，请按照错误消息中的说明提交错误报告 –– 您没有发布的部分(或者可能因为版本不同而没有发布？):

代码语言：javascript
复制
Explanation: Refer to the bug-reporting instructions at http://www.cgal.org/bug_report.html


libc++abi.dylib: terminating with uncaught exception of type CGAL::Precondition_exception: CGAL ERROR: precondition violation!


Expr: (m_traits.compare_y_at_x_2_object()(p, cv) == EQUAL) && compare_xy(cv.left(), p) == SMALLER && compare_xy(cv.right(), p) == LARGER


File: /usr/local/include/CGAL/Arr_segment_traits_2.h Line: 706

据我从该文件中所见，函数 throwing 在给定分割点的情况下将一条曲线分成两条子曲线。前提是分割点在曲线上，而不是终点。显然，是否满足这个前提条件是调用者的问题，另一个 CGAL 例程。换句话说，您的输入没有任何问题。问题出在 CGAL 实现上，或者更准确地说，是它处理所用数字表示不精确的方式。

根据 sloriot 的评论进行编辑。

如果我更换

typedef CGAL::Simple_cartesian K;
与

typedef CGAL::Exact_predicates_exact_constructions_kernel K;
并为 Area 和 totalArea 使用适当的类型(我只是使用了 auto 和 decltype(Area) ，分别)，代码编译(你必须将它链接到 libgmp 和 libmpfr)并且运行没有崩溃，报告

#cgal#cout#point#polygon#异常