Division by zero using Epeck to compute the position of vertices in SVD

qwenger commented 4 months ago

Issue Details

Hi and thank you very much for developing CGAL,

I use the library for computing segment Voronoi diagrams and I came across an issue (hereafter drastically simplified to a MWE).

The following code crashes at runtime with a precondition exception / division by zero error when computing the position of one of the vertices. From what I could investigate, during the building phase the two segments are found to be non-parallel, thus a vertex is created with both segments and the point as sites. However, on coordinate computation the value sqrt(D1D2) - a1a2 - b1b2 in compute_pll evaluates to zero and the code crashes.

After much tweaking and testing I came to the idea of swapping the kernel in favor of Exact_predicates_exact_constructions_kernel_with_sqrt, and then the code executes smoothly.

For me this raises a few questions:

Is Exact_predicates_exact_constructions_kernel_with_sqrt "guaranteed" to be stable in this use case or am I at risk of similar crashes for different inputs?
Is this to be considered a bug in the basic Epeck or rather an acceptable behavior given its properties?
In the latter case, would it make sense to write something about it in the docs? Reading around the exact computation paradigm and the kernel definitions, I (naively?) thought that Epeck would solve this kind of precision issues.

FWIW I also just came across following somewhat similar reports: https://github.com/CGAL/cgal/issues/4195 https://github.com/CGAL/cgal/issues/4196

EDIT: it seems that https://github.com/CGAL/cgal/pull/7976#issuecomment-1887194744 somewhat answers my point 2.

Thanks in advance!

Source Code

#include <iostream>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_policies_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Voronoi_diagram_2.h>

using K = CGAL::Exact_predicates_exact_constructions_kernel;
using Point = CGAL::Point_2<K>;
using Gt = CGAL::Segment_Delaunay_graph_traits_2<K>;
using SDG2 = CGAL::Segment_Delaunay_graph_2<Gt>;
using AT = CGAL::Segment_Delaunay_graph_adaptation_traits_2<SDG2>;
using AP = CGAL::Segment_Delaunay_graph_degeneracy_removal_policy_2<SDG2>;
using VoronoiDiagram = CGAL::Voronoi_diagram_2<SDG2, AT, AP>;
using Site = AT::Site_2;

int main() {
    VoronoiDiagram vd;
    vd.insert(Site::construct_site_2(Point(300.405065, 113.2793247), Point(310.0, 112.0)));
    vd.insert(Site::construct_site_2(Point(310.0,      112.0      ), Point(340.0, 108.0)));
    vd.insert(Site::construct_site_2(Point(  0.0,        0.0      )));
    for (auto vertex_iter = vd.vertices_begin(); vertex_iter != vd.vertices_end(); ++vertex_iter) {
        std::cout << vertex_iter->point() << std::endl;
    }
    return 0;
}

Environment

Operating system (Windows/Mac/Linux, 32/64 bits): Windows & Linux 64 bits
Compiler: MSVC & GCC
Release or debug mode: both (division by zero in release mode, precondition exception in debug)
Specific flags used (if any):
CGAL version: 5.6.1
Boost version: 1.83.0
Other libraries versions if used (Eigen, TBB, etc.):

afabri commented 4 months ago

Not a reply, but a question out of curiosity: What is the goal of https://gitlab.com/matpi/carina ?

qwenger commented 4 months ago

Oh, you found that already! ;)

The goal of the Carina package is to provide a Python library for computing straight skeletons and medial axes of polygons, based on CGAL. I created it because

AFAIK there are no up-to-date Python bindings for CGAL containing both (maybe I'm blind but I couldn't find access to straight skeleton generation in the official CGAL swig bindings, and other wrapping packages on PyPI are outdated).
My use case is splitting the vicinity of polygonal features in GIS by nearest neighbor, so I have special requirements on the outputs (A. medial axes are required, so some filtering of the raw Voronoi diagram is needed; B. parabola segments are discretized and the user can control the precision of that discretization; C. each vertex of the input polygon should generate exactly one single edge, even reflex vertices - unlike the standard definition of MAs). I didn't find built-in support for those requirements in the available packages, hence I created one.

Did I miss something? Is all of this already available in an official package? If so (or if that happen in the future) I can happily deprecate Carina as to not reinvent the wheel.

afabri commented 1 week ago

@MaelRL if I understand it correctly the vertex_iter evaluates to a face handle. We would still be in the realm of exact computing, if we had a function which returns a bool telling us if the dual is a finite point or not, wouldn't we?

CGAL / cgal