weighted periodic 3d alpha shapes crash

[mglisse]

Issue Details

$ ./a.out Periodic Delaunay computed zsh: floating point exception ./a.out

Source Code

(based on ex_weighted_periodic_alpha_shapes_3.cpp but using the more useful GENERAL mode)

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Alpha_shape_3.h>
#include <CGAL/Alpha_shape_cell_base_3.h>
#include <CGAL/Alpha_shape_vertex_base_3.h>
#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
#include <CGAL/Periodic_3_regular_triangulation_3.h>
#include <vector>

using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using PK = CGAL::Periodic_3_regular_triangulation_traits_3<Kernel>;
using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>;
using Vb = CGAL::Regular_triangulation_vertex_base_3<PK,DsVb>;
using AsVb = CGAL::Alpha_shape_vertex_base_3<PK,Vb>;
using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>;
using Cb = CGAL::Regular_triangulation_cell_base_3<PK,DsCb>;
using AsCb = CGAL::Alpha_shape_cell_base_3<PK,Cb>;
using Tds = CGAL::Triangulation_data_structure_3<AsVb,AsCb>;
using P3RT3 = CGAL::Periodic_3_regular_triangulation_3<PK,Tds>;
using Alpha_shape_3 = CGAL::Alpha_shape_3<P3RT3>;

using Point_3 = P3RT3::Bare_point;
using Weighted_point_3 = P3RT3::Weighted_point;

int main() {
  std::vector<Weighted_point_3> wp = {
    Weighted_point_3(Point_3(-0.905233,-0.424915,9.11647e-08),0),
    Weighted_point_3(Point_3(0.50771,-0.861528,5.5397e-08),0),
    Weighted_point_3(Point_3(0.604646,-0.796495,9.5223e-08),0),
    Weighted_point_3(Point_3(0.821595,0.570072,6.06969e-08),0),
    Weighted_point_3(Point_3(-0.0669612,-0.997756,1.56679e-08),0),
    Weighted_point_3(Point_3(-0.951394,-0.307977,2.18257e-08),0),
    Weighted_point_3(Point_3(0.521555,0.853218,6.37552e-08),0),
    Weighted_point_3(Point_3(0.728106,-0.685465,7.71358e-08),0),
    Weighted_point_3(Point_3(-0.564656,-0.825326,2.83315e-08),0),
    Weighted_point_3(Point_3(0.399941,-0.916541,9.49327e-08),0),
    Weighted_point_3(Point_3(0,0,0),0) };

  P3RT3 prt(PK::Iso_cuboid_3(-1,-1,-1,1,1,1));
  prt.insert(wp.begin(), wp.end(), true);
  if (prt.is_triangulation_in_1_sheet()) prt.convert_to_1_sheeted_covering();
  std::cerr << "Periodic Delaunay computed\n";
  Alpha_shape_3 as(prt, 0, Alpha_shape_3::GENERAL);
}

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): linux 64
• Compiler: gcc-7 or clang-6
• Release or debug mode: both
• CGAL version: 4.11 or master
• Boost version: 1.62

[MaelRL]

The problem is that the Is_gabriel functions were badly getting offsets when using 27-sheeted triangulations, resulting in the same point being passed twice to the power_side_of_bounded_power_sphere predicate. Same issue in P3DT3.

Note that you shouldn't use true in the constructor of the periodic triangulation when you only have a handful of points, since you'll spend an absurd amount of time adding and removing dummies.

[mglisse]

Thank you for looking at it.

Note that you shouldn't use true in the constructor of the periodic triangulation when you only have a handful of points, since you'll spend an absurd amount of time adding and removing dummies.

Well, we have a program that reads user files, and it seemed natural to optimize it for the case with many points rather than the case with few points (or more simply we copy-pasted ex_weighted_periodic_alpha_shapes_3.cpp and trusted it was advertising the right options). If you think false is a better default (i.e. it isn't significantly slower for large point sets), we'll use that.

[MaelRL]

The tag has quite an influence: for your case, this is the difference between instantaneous and minutes. It's basically expensive to use true if your point set will switch back to 27-sheets when you remove the dummy points that are added by the heuristic.

Since you have a wrapper around the periodic triangulations, can't you decide on a per-case basis whether the tag should be true or false ? With a uniform distribution, about a couple hundreds points will create a safe 1-sheet triangulation and you can use the tag.

[mglisse]

Since you have a wrapper around the periodic triangulations, can't you decide on a per-case basis whether the tag should be true or false ? With a uniform distribution, about a couple hundreds points will create a safe 1-sheet triangulation and you can use the tag.

I don't really see how I could guess if the user's point set is going to be nice. The first point set our tester tried had all the points on a small circle very close to the center of the cuboid, so just counting the points is not a good enough heuristic. How much performance do I lose by using false for a nice, large point set? An additive constant (as soon as it is 1-sheeted, the speed on insertion becomes the same, so only the first few points are slower)? Then using false by default sounds good.

On the other hand, when we want persistence afterwards, we are essentially only interested in the 1-sheeted case. If it cannot fit on 1 sheet, we will error out, because the persistence computed on the multiple sheets is rather weird. Ideally we would subdivide the triangulation enough that we can produce a filtered simplicial complex with the right properties, but that's too much work for now. If using "true" is only slow in the case where we are going to produce an error, it is tempting not to care...

[MaelRL]

I don't really see how I could guess if the user's point set is going to be nice. The first point set our tester tried had all the points on a small circle very close to the center of the cuboid, so just counting the points is not a good enough heuristic.

I'd say it's better than nothing though.

How much performance do I lose by using false for a nice, large point set? An additive constant (as soon as it is 1-sheeted, the speed on insertion becomes the same, so only the first few points are slower)? Then using false by default sounds good.

Obviously it depends on how large the point set is, but it's also a significant cost to introduce manually points until you can switch to 1-cover when you know that you'll get a 1-cover with your point set (significant enough for this heuristic to be introduced in the first place :)). However, I guess that this is less significant than the other way (maybe @MoniqueTeillaud has done more testing on this ?).

On the other hand, when we want persistence afterwards, we are essentially only interested in the 1-sheeted case. If it cannot fit on 1 sheet, we will error out, because the persistence computed on the multiple sheets is rather weird. Ideally we would subdivide the triangulation enough that we can produce a filtered simplicial complex with the right properties, but that's too much work for now.

We debated recently returning "unique" elements in the alpha shapes regardless of the covering (if that's what you mean). I'm rather against that: the alpha shape should be consistent with the triangulation.

If using "true" is only slow in the case where we are going to produce an error, it is tempting not to care...

Comparing the two options, If you're anyway aborting if the triangulation is not a 1-cover, then I'd probably leave it to true.

[mglisse]

We debated recently returning "unique" elements in the alpha shapes regardless of the covering (if that's what you mean). I'm rather against that: the alpha shape should be consistent with the triangulation.

Since currently we want to build a simplicial complex, the "unique" thing you describe doesn't sound like it would help us. If a point is connected to one of its copies for instance, I cannot handle a self-loop. What I was talking about was subdividing the periodic triangulation (not maintaining the Delaunay property, but giving the new simplices filtration values consistent with the Delaunay triangulation) until the 1-cover becomes a simplicial complex. To go for overkill, 2 barycentric subdivisions are sufficient.

If we had the datastructures to handle more general cell complexes, this unique reporting might make sense, although we would need to make sure we transfer the connectivity properly... Trying to think, for the case of a single point in 2D, we would report 1 0-cell, 3 self-loops attached to it, and 2 2-cells, each attached to the 3 loops... rather inconvenient to work with, but it does seem meaningful.

I don't know if there are cases where getting the alpha-complex of the 27-cover is useful, I don't have one in mind at the moment, but that's a very weak argument. I guess there are cases where, ignoring the essential cycles, one could simply divide the multiplicity of persistence intervals by 27, but that's likely cases where we can actually get a single sheet.

[lrineau]

Fixed by #2613.