Я использую CGAL, чтобы вычислить площадь пересечения двух выпуклых многоугольников. Короткий демонстрационный код для этого был опубликован в принятом ответе на этот вопрос. Однако, когда я изменяю этот код, чтобы использовать интересующие меня полигоны, 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)
Я не понимаю, в чем причина ошибки. Любая помощь будет оценена.
Я могу воспроизвести эту ошибку (используя CGAL 4.9 на MacOS с clang ++). Насколько я понимаю, не должно произойти неисследованного исключения этого типа, другими словами, вы обнаружили ошибку в CGAL. Итак, подайте отчет об ошибке, как указано в сообщении об ошибке
–– часть, которую вы не опубликовали (или, возможно, не получили из-за другой версии?):
Объяснение: Обратитесь к инструкциям по сообщению об ошибках в http://www.cgal.org/bug_report.html
libc ++ abi.dylib: завершается с необработанным исключением типа CGAL :: Precondition_exception: CGAL ERROR: нарушение предусловия!
Пример: (m_traits.compare_y_at_x_2_object () (p, cv) == РАВНО) && compare_xy (cv.left (), p) == МАЛЕНЬКИЙ && Compare_xy (Cv.right (), P) == БОЛЬШОЙ
Файл: /usr/local/include/CGAL/Arr_segment_traits_2.h Строка: 706
Насколько я вижу из этого файла, функция throw разбивает кривую на две подогнуты с учетом точки разделения. Предварительным условием является то, что точка разделения лежит на кривой и не является конечной точкой. Понятно, соблюдается ли это предварительное условие или нет — это вопрос вызывающей стороны, еще одна процедура CGAL. Другими словами, в вашем входе нет ничего плохого. Что не так, так это реализация CGAL, или, точнее, ее способ справиться с неточностью используемого представления чисел.
Редактировать следующий комментарий sloriot.
Если я заменю
typedef CGAL::Simple_cartesian<double> K;
с
typedef CGAL::Exact_predicates_exact_constructions_kernel K;
и использовать соответствующий тип для Area
а также totalArea
(Я просто использовал auto
а также decltype(Area)
соответственно) код компилируется (необходимо связать его с libgmp
а также libmpfr
) и работает без сбоев, сообщая
totalArea::0.000876944
Других решений пока нет …