CGALでポリゴンを交差する場合の前提条件の例外

Question

2つの凸多角形の交点の面積を計算するためにCGALを使用しています。これを行うための短いデモコードが、this質問への受け入れられた回答に掲載されました。しかし、興味のあるポリゴンを使用するようにそのコードを変更すると、CGALはCGAL :: intersection（）ルーチンの深いところからランタイム例外をスローします。

ここでは、上にリンクされたSO質問からのコピーアンドペーストですが、それは自分のポリゴンを使用し、各ポリゴンに関するいくつかの診断をプリントして凸であることを示し、CCW巻線順序。

#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 << "created the polygon P:" << endl; 
    cout << polyP << endl; 
    cout << endl; 
    IsSimple = polyP.is_simple(); 
    IsConvex = polyP.is_convex(); 
    IsClockwise = (polyP.orientation() == CGAL::CLOCKWISE); 
    Area  = polyP.area(); 
    cout << "polygon P is"; 
    if (!IsSimple) cout << " not"; 
    cout << " simple." << endl; 
    cout << "polygon P is"; 
    if (!IsConvex) cout << " not"; 
    cout << " convex." << endl; 
    cout << "polygon P is"; 
    if (!IsClockwise) cout << " not"; 
    cout << " clockwise oriented." << endl; 
    cout << "the area of polygon P is " << Area << endl; 
    cout << endl; 

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

    // Compute the intersection 
    std::list<Polygon_with_holes_2> polyI; 
    CGAL::intersection(polyP, polyQ, std::back_inserter(polyI)); 

    double totalArea = 0; 
    typedef std::list<Polygon_with_holes_2>::iterator LIT; 
    for(LIT lit = polyI.begin(); lit!=polyI.end(); lit++){ 
     totalArea+=lit->outer_boundary().area(); 
    } 
    cout << "TotalArea::" << totalArea; 

    return 0; 
} 

ここにテストプログラムの出力があります。

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) 

エラーの原因を理解できません。どんな助けもありがとう。

Answer 1

私はこのエラーを再現できます（clang ++のMacOSではCGAL 4.9を使用しています）。私が理解する限り、このタイプのキャッチされない例外は起こりそうにない、つまり、あなたはCGALのバグを発見しました。 （？または恐らくは異なるバージョンの取得できませんでした） - エラーメッセージ に指示して、バグレポートを提出しますあなたが投稿していなかった一部：

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

私の知る限りとしてそのファイルから見ることができます。スローする関数は、分割点を指定してカーブを2つのサブカーブに分割します。前提条件は、スプリットポイントがカーブ上にあり、終点ではないことです。明らかに、この事前条件が満たされているかどうかは、別のCGALルーチンの呼び出し側の問題です。言い換えれば、入力に間違いはありません。間違ったのは、CGALの実装、またはより正確には、使用される数値表現の不正確さを処理する方法です。

---

次のsloriotのコメントを編集します。あなたがリンクする必要があり

私は（私は単に、それぞれ、autoとdecltype(Area)を使用）

typedef CGAL::Exact_predicates_exact_constructions_kernel K; 

で

typedef CGAL::Simple_cartesian<double> K; 

を交換し、AreaとtotalAreaのための適切なタイプを使用する場合は、コードのコンパイル（それはlibgmpとlibmpfrに対して）、実行中にクラッシュすることなく実行します。

totalArea::0.000876944