CGAL users discussion list

[Efi Fogel <>]

The choice whether to use Epic or Epec is not for you (the user) to make. If it had been, then all users would have chosen Epic over Epec.

If constructions take place while computing, and the implementation is sensitive to inexactness, then Epec must be used. Sometimes constructions occur only to produce intermediate results; even then, Epec is required. If the input points are far apart from each other this still cannot guarantee correct computation, because the computation can still reach  a stage of uncertainty, and yes, when the computation reaches a stage of uncertainty, anything can happen, including a segfault.

Let me reiterate, in order to avoid segfaults caused by limited precision arithmetic, we need to avoid reaching stages of uncertainty, and Epec does exactly that.

   ____  _        ____             _
  /_____/_) o    /__________  __  //
 (____ (   (    (    (_/ (_/-(-'_(/
                         _/

On Tue, Sep 13, 2016 at 11:29 PM, Christopher Mitchell <> wrote:

I've been running into a weird segfault with constrained Delauney triangulation of two perfectly innocent polygons with coordinates on the order of 10^2 units as well; are you telling me that the underlying constructions are expected to segfault in some way if there's a math error? That seems like bafflingly bugged behavior to me. I note that applying a scaling factor of 10^-2 to the coordinates in the example I have here made the code work properly; is there some sort of unstated bound on coordinate magnitude for inexact constructions that I should be obeying, like translating and scaling all of my operations to [-10, 10], for example?

I should clarify that my hesitation to switch to exact constructions is twofold: (1) the slowdown of using those constructions. I'm attempting to write very performant code that will be executed on millions of objects; (2) I see that exact constructions are not drag-and-drop; I'll have to sprinkle conversions between doubles and gmp numbers all over the place.

Thanks in advance for your elucidations!

On Tue, Sep 13, 2016 at 5:42 AM, Efi Fogel <> wrote:

When I run the program as is (using inexact construction, I get a core dump, which is most probably caused by the inexact constructions.

When I run the program with an exact construction kernel, I get the correct result.

You can simply:

void dumpPoly(CGALPolygon& poly) { std::cout << poly << std::endl; }

When using Exact_predicates_exact_constructions_kernel

   ____  _        ____             _
  /_____/_) o    /__________  __  //
 (____ (   (    (    (_/ (_/-(-'_(/
                         _/

On Mon, Sep 12, 2016 at 10:48 PM, Christopher Mitchell <> wrote:

Sebastien: Thanks for the advice. Unfortunately, for speed purposes I'd prefer to stick with inexact constructions. This polygon pair doesn't seem like it's really approaching the limits of numeric precision in terms of the huge overlap area.

Efi: Ask and you shall receive. This code is a little unwieldy for a minimum working example (102 lines), but I figure its longevity will be better here than in a Pastebin. To my great frustration, I tried removing points, and the problem remained, but when I tried subtracting 100.0 from every X coordinate, it indicated that the polygons intersected. Working with doubles, the intersection definitely shouldn't be hitting odd precision problems with coordinates that are only on the order of 10^2 units. This seems like it could possibly be a bug?

Thanks,

Christopher

// -----snip-----

#include <stdio.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Boolean_set_operations_2.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel CGAL_Kernel;
typedef CGAL_Kernel::Point_2 CGALPoint;
typedef CGAL::Polygon_2<CGAL_Kernel> CGALPolygon;

double pts_p1[][2] = {
        {138.434248, 54.124268},
        {156.251888, 53.447439},
        {159.268861, 55.635141},
        {168.750594, 82.497614},
        {165.721545, 83.640381},
        {159.550635, 86.191099},
        {149.385779, 90.040283},
        {130.971384, 97.216211},
        {119.770956, 101.497880},
        {112.218620, 104.432790},
        {109.481271, 93.899521},
        {108.693051, 87.351571},
        {108.600034, 80.589390},
        {109.019489, 73.333375},
        {109.398422, 60.659113},
        {109.614696, 54.943998},
        {128.132617, 54.352682}
};

double pts_p2[][2] = {
        {139.434399, 57.599339},
        {155.707971, 56.738629},
        {156.812714, 57.672514},
        {165.721545, 83.640381},
        {159.550635, 86.191099},
        {149.385779, 90.040283},
        {130.971384, 97.216211},
        {119.770956, 101.49788},
        {112.218620, 104.43279},
        {109.481271, 93.899521},
        {108.693051, 87.351571},
        {108.600034, 80.589390},
        {109.019489, 73.333375},
        {109.398422, 60.659113},
        {110.103104, 59.601234},
        {127.188178, 57.950231},
        {126.599293, 61.485162},
        {124.600343, 71.709878},
        {139.439351, 66.181265},
        {139.584575, 60.776208}
};

int nodes2Polygon(const double pts[][2], const size_t size, CGALPolygon& poly) {
        printf("Generating %zu-point poly\n", size);
        // Copy points into polygon. Omit the last point if it's a copy of the first point.
        for (size_t i = 0; i < size; i++) {
                if (i == 0 || i != (size - 1) || (pts[i][0] != pts[0][0] || pts[i][1] != pts[0][1])) {
                        CGALPoint point(pts[i][0], pts[i][1]);
                        // Make sure the point is not identical to the previous point in the poly.
                        if (!poly.size() || point != *std::prev(poly.vertices_end())) {
                                poly.push_back(point);
                        }
                }
        }

        if (!poly.is_simple()) {
                return -1;
        }

        // Correct orientation to make it an outline, not a hole
        if (poly.orientation() == CGAL::CLOCKWISE) {
                poly.reverse_orientation();
        }

        return 0;
}

void dumpPoly(CGALPolygon& poly) {
        for(auto it = poly.vertices_begin(); it != poly.vertices_end(); it++) {
                printf("%f %f\n", it->x(), it->y());
        }
        printf("\n");
}

int main(int argc, char* argv[]) {
        CGALPolygon p1;
        CGALPolygon p2;

        nodes2Polygon(pts_p1, sizeof(pts_p1) / sizeof(double) / 2, p1);
        dumpPoly(p1);

        nodes2Polygon(pts_p2, sizeof(pts_p2) / sizeof(double) / 2, p2);
        dumpPoly(p2);

        if (CGAL::do_intersect(p1, p2)) {
                printf("Polygons intersect\n");
        } else {
                printf("Polygons do not intersect\n");
        }

        if (CGAL::do_intersect(p2, p1)) {
                printf("Polygons intersect\n");
        } else {
                printf("Polygons do not intersect\n");
        }
        return 0;
}

On Mon, Sep 12, 2016 at 6:08 AM, Sebastien Loriot (GeometryFactory) <> wrote:

The kernel used does not have exact constructions.
Try again with CGAL::Exact_predicates_exact_constructions_kernel

Sebastien.

On 09/12/2016 12:00 PM, Efi Fogel wrote:

It would easier if you just send the code and the input polygon.
This way we can quickly try to reproduce the problem.

    ____  _        ____             _
   /_____/_) o    /__________  __  //
  (____ (   (    (    (_/ (_/-(-'_(/
                          _/



On Mon, Sep 12, 2016 at 1:03 AM, Christopher Mitchell <
<mailto:>> wrote:

    Thanks for the quick response. The code is sufficiently trivial (ie,
    a single if (CGAL::do_intersect(poly1, poly2)) { ... }) that I'm not
    sure it would help. I'm using a simple loop to convert a set of
    non-CGAL points into a Polygon_2, making sure to not duplicate the
    starting point, as you can see from the image attached to the first
    post. My polygons are CGAL::Polygon_2<CGAL_Kernel>, where
    CGAL_Kernel is CGAL::Exact_predicates_inexact_constructions_kernel.
    I hope this helps. Thanks!

    On Sat, Sep 10, 2016 at 6:25 PM, Efi Fogel <
    <mailto:>> wrote:

        What number type are you using?  Better yet, could you share the
        code?


        On Sep 10, 2016 18:31, "Christopher Mitchell" <
        <mailto:>> wrote:

            I'm attempting to determine if two Polygon_2s overlap in
            anyway, ie, share any portion of their interiors, including
            one being entirely contained by the other. In most cases
            this works, but for a few polygon pairs, do_intersect()
            returns false unexpectedly. I make sure (1) all polygons are
            simple (2) all polygons are not oriented CGAL::CLOCKWISE (I
            call .reverse_orientation() on the polys if they are
            CGAL::CLOCKWISE), and (3) that the starting vertex is not
            duplicated. I have attached an example of a polygon pair for
            which do_intersect() returns false. Thanks in advance for
            any insight; I'd be happy to provide additional information
            as necessary.

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss