CGAL/cgal

Triangulation_full_cell segfault bug for degenerate cases

Opened this issue · 1 comments

Issue Details

A Delaunay triangulation of degenerate points has maximal cells that are of dimension lower than the dimension of the ambient space. The Triangulation_ds_full_cell class does not recognise this, and its vertices.end() method results in a segfault.

Source Code

#include <CGAL/Delaunay_triangulation.h>
#include <CGAL/Dimension.h>
#include <CGAL/Epick_d.h>
#include <CGAL/Gmpzf.h>
#include <CGAL/Triangulation.h>
#include <Eigen/Dense>
#include <iostream>

using Kernel_d                     = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
using Triangulation_data_structure = CGAL::Triangulation_data_structure<
	Kernel_d::Dimension,
	CGAL::Triangulation_vertex<Kernel_d, long>,
	CGAL::Triangulation_full_cell<Kernel_d>>;
using DelaunayTriangulation = CGAL::Delaunay_triangulation<Kernel_d, Triangulation_data_structure>;
using Point_d               = Kernel_d::Point_d;

using Eigen::MatrixXd, Eigen::Index;

using std::tuple, std::vector, std::cout;

// Convert a matrix of coordinate column vectors to a vector of CGAL points.
auto matrix_columns_to_points_vec(const MatrixXd& X) -> vector<Point_d> {
	vector<Point_d> points(X.cols());
	auto            cols_begin = X.colwise().cbegin();
	auto            cols_end   = X.colwise().cend();
	for (auto&& [i, column] = tuple<Index, decltype(cols_begin)>{0, cols_begin}; column != cols_end;
	     column++, i++) {
		points[i] = Point_d(column->cbegin(), column->cend());
	}
	return points;
}

auto main() -> int {
	MatrixXd X{
		{0.0, 1.0, 2.0},
		{0.0, 0.0, 0.0},
	};
	auto delX   = DelaunayTriangulation(static_cast<int>(X.rows()));
	auto points = matrix_columns_to_points_vec(X);
	delX.insert(points.begin(), points.end());

	std::cout << "delX.dimension() = " << delX.current_dimension() << '\n';
	for (auto cell_it = delX.finite_full_cells_begin(); cell_it != delX.finite_full_cells_end();
	     cell_it++) {
		int temp = 0;
		// Iterate over the vertices of the cell
		for (auto&& vert_it = cell_it->vertices_begin(); vert_it != cell_it->vertices_end();
		     vert_it++) {
			temp += 1;
			// The following line will segfault
			// *vert_it;
		}
		cout << "Finite cell vertex iterator had " << temp << " iterations.\n";
	}
}

Output

delX.dimension() = 1
Finite cell vertex iterator had 3 iterations.
Finite cell vertex iterator had 3 iterations.

Environment

  • Operating system (Windows/Mac/Linux, 32/64 bits): Ubuntu 24.04 x64
  • Compiler: GCC 13.0.3
  • Release or debug mode: Both
  • Specific flags used (if any): -frounding-math
  • CGAL version: 6.0.1
  • Boost version: 1.88
  • Other libraries versions if used (Eigen, TBB, etc.): Eigen 3.4
  • Inserting 1 point is enough to reproduce.
  • I see at least one place in CGAL that uses s->vertices_begin() + cur_dim + 1 as the "end" iterator.
  • A cell doesn't seem to know the current dimension (it is passed as argument to the functions that need it).
  • The doc does not seem to say anything about this (or I missed it)

Yes, it does look like we need to modify something in CGAL, either document things better, or move vertices_begin/end so they are called tr.vertices_end(s) instead of s->vertices_end(). And maybe audit the code in Triangulation to make sure that this behavior does not cause trouble anywhere.