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

[qwenger]

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:

1. 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?
2. Is this to be considered a bug in the basic Epeck or rather an acceptable behavior given its properties?
3. 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: #4195 #4196

EDIT: it seems that #7976 (comment) 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]

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

[qwenger]

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

1. 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).
2. 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.