Alpha shape filtration doesn't produce a simplicial complex

[mrzv]

Issue Details

Alpha shape filtration, CGAL::Alpha_shape_3<Delaunay>::filtration() doesn't produce a simplicial complex. A full example, including input data and output, lives in a gist. If you compile the actual code, it produces output, which contains triangles

114 163 166
114 117 166
51 114 166

but not the edge 114 166.

Source Code

If you clone the gist, it has a Makefile that should compile a runnable example.

Environment

• Operating system: Mac, 64 bits
• Compiler: Apple LLVM version 10.0.1 (clang-1001.0.46.4)
• Release or debug mode: -O3
• Specific flags used (if any): -DCGAL_USE_GMP -DCGAL_USE_MPFR, see the Makefile for details.
• CGAL version: 4.14
• Boost version: 1.69.0

[mrzv]

I should add that most of the time, filtration() produces the right thing. But on this particular input (which happens to be very degenerate), it runs into this problem.

[sloriot]

If you look at the documentation page, you will see the ExactAlphaComparisonTag. Try setting it to CGAL::Tag_true and it should work correctly.

[lrineau]

@sloriot commented on May 14, 2019, 6:41 AM GMT+2:

If you look at the documentation page, you will see the ExactAlphaComparisonTag. Try setting it to CGAL::Tag_true and it should work correctly.

Your link seems the wrong one.

[mrzv]

When I add Tag_true to the Alpha_shape_3, I get the following compilation error:

In file included from compute-filtration.cpp:11:
/usr/local/include/CGAL/Alpha_shape_3.h:124:3: error: static_assert failed due to requirement 'boost::is_same<NT, typename Delaunay_triangulation_3<Epick,
      Triangulation_data_structure_3<Alpha_shape_vertex_base_3<Epick, Default, Boolean_tag<false>, Boolean_tag<false> >, Alpha_shape_cell_base_3<Epick, Default, Boolean_tag<false>, Boolean_tag<false> >,
      Sequential_tag>, Location_policy<Fast>, Default>::Cell::NT>::value' "(boost::is_same<NT,typename Dt::Cell::NT>::value)"
  CGAL_static_assertion( (boost::is_same<NT,typename Dt::Cell::NT>::value) );
  ^                       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/local/include/CGAL/assertions.h:130:6: note: expanded from macro 'CGAL_static_assertion'
     static_assert(EX, #EX)
     ^             ~~
compute-filtration.cpp:22:20: note: in instantiation of template class 'CGAL::Alpha_shape_3<CGAL::Delaunay_triangulation_3<CGAL::Epick,
      CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick,
      CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL::Location_policy<CGAL::Fast>, CGAL::Default>, CGAL::Boolean_tag<true> >' requested here
using FT         = AlphaShape::FT;
                   ^

What does it mean, and how do I fix it?

My second question: the warning in the documentation says that Alpha_shape_3 won't work with periodic triangulations, with Tag_true, but I need the periodic case. Is there anything that can be done in that situation?

[mrzv]

Another error further down is also concerning:

compute-filtration.cpp:100:72: error: no member named 'GENERAL' in 'CGAL::Alpha_shape_3<CGAL::Delaunay_triangulation_3<CGAL::Epick,
      CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick,
      CGAL::Default, CGAL::Boolean_tag<false>, CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL::Location_policy<CGAL::Fast>, CGAL::Default>, CGAL::Boolean_tag<true> >'
    AlphaShape as(dt, std::numeric_limits<FT>::infinity(), AlphaShape::GENERAL);

[MaelRL]

Note that the tag ExactAlphaComparisonTag has to be set for the class Alpha_shape_3 but also for the vertex and cell base classes, see this example.

[mrzv]

Ah, that fixed it. Thanks, @MaelRL.

My question about periodic triangulations remains. How do I solve this problem, if I run into it there?

[MaelRL]

To quote 63d293765fe88e6c5756b4a250e2145bf0b6d489:

When we set the exact comparison tag to true, we use a lazy evaluation
of predicates and store pointers to the arguments of the predicates (that is,
pointer to points). However, in the case of a periodic setting, the points 
are only temporary objects constructed on the fly, and it is thus dangerous
to take pointers to those temporary values because the pointers quickly become invalid.

Thus, periodic triangulations are not allowed to use the exact tag.

A possible way to make it work is to define a small class:
Periodic_triangulation_with_stored_real_points that would inherit the base
(periodic) triangulation but store all the points of its simplices.
Then, the function point(face_handle, int) would return the correct point from this
container, of which it would be safe to take a pointer.

Thus your steps will be:

1. Create the periodic triangulation.
2. Initialize your permanent point container. Note that your point container will be a map[Face_handle fh, int i] = Point := your_base_periodic_triangulation.point(fh, i). I don't think you need point(vertex_handle).
3. Run the alpha shape.

You also have to remove the static assertions that will force compilation failure when the triangulation has the periodic tag in the Alpha_shape_3 class.

I haven't created such class not because it's necessarily complicated to do, but because I wasn't sure anybody would ever need it, and it's a bit awkward. Let me know if everything is clear and if you get stuck.

[mrzv]

That sounds complicated, and unmanageable for downstream users because of having to modify CGAL assertions.

Can I ask what the pointer (to points) lifetime is? I only need them very briefly to look up the index in a map, and thus get the vertex id.

[mrzv]

Also, unrelated to periodic triangulations, when I run filtration_with_alpha_values() and get the FT argument, which now seems to be 'CGAL::internal::Lazy_alpha_nt_3<CGAL::Epick, true, CGAL::Boolean_tag<false> >, how do I convert it to a float (or my input coordinate type)? Simple assignment doesn't work.

[MaelRL]

Since there's apparently a need, I can modify the code to add this wrapping class officially in an upcoming CGAL version (and a branch in the meantime). I was only explaining what you'd have to do to make it work in case you wanted it right now.

I don't understand your question about the pointer lifetime, the alpha shape algorithm needs to keep pointers to the points of your triangulation when it is building the alpha shapes.

You can switch back to float values using CGAL::to_double().

[mrzv]

Hm, to_double() seems to say that I'm not guaranteed to get the same number, if I run it on the same input. This can be a problem for the alpha shape, if a simplex and its face have the same value, but to_double() gives a higher value for the face than for a simplex. Is there a way around this?

As for the pointer lifetime, I'm referring to the first paragraph of your quote from a commit. It says the points are constructed on the fly, and thus it's dangerous to use the pointers. My question is what is the lifetime of these points constructed on the fly? If I simply access them inside the iterator callback from the filtration(), are they guaranteed to be valid? If yes, then it would solve my problem, without the more complicated workaround.

As for incorporating the change into CGAL, I'd welcome that. I can certainly wait for the next version: this is a sufficiently rare problem that it can wait a little.

[MaelRL]

I'm not as sure for the first question (@sloriot will confirm -- or not), but I think there are two ways:

• You can call exact() before CGAL::to_double() to narrow the interval width.
• You can use set_relative_precision_of_to_double() to set a low-enough precision (below the double precision).

As for the second question, the point(fh,i) function is called somewhere within the alpha shape construction and thus the lifetime of the point is too short to even outlive the creation of the alpha shape. You will not be able to obtain any filtration because you cannot even construct the alpha shape properly in the first place.

I will try to give you a prototype in a branch in not too long.

[mrzv]

@MaelRL Thanks.

[mrzv]

@MaelRL Sorry, one more question. I'm trying to add this Tag_true to the weighted alpha shape, and I'm getting some problem that sounds like I'm forgetting to add Tag_true to one of the classes. Does anything jump out at you below?

    using K         = CGAL::Exact_predicates_inexact_constructions_kernel;

#if (CGAL_VERSION_MAJOR == 4 && CGAL_VERSION_MINOR <= 9) || (CGAL_VERSION_MAJOR < 4)
    using Gt        = CGAL::Regular_triangulation_euclidean_traits_3<K>;
    using Vb        = CGAL::Alpha_shape_vertex_base_3<Gt,CGAL::Default,CGAL::Tag_true>;
    using Fb        = CGAL::Alpha_shape_cell_base_3<Gt,CGAL::Default,CGAL::Tag_true>;
    using TDS       = CGAL::Triangulation_data_structure_3<Vb,Fb>;
    using Delaunay  = CGAL::Regular_triangulation_3<Gt,TDS>;
#else
    using Rvb       = CGAL::Regular_triangulation_vertex_base_3<K>;
    using Vb        = CGAL::Alpha_shape_vertex_base_3<K, Rvb, CGAL::Tag_true>;
    using Rcb       = CGAL::Regular_triangulation_cell_base_3<K>;
    using Cb        = CGAL::Alpha_shape_cell_base_3<K, Rcb, CGAL::Tag_true>;
    using TDS       = CGAL::Triangulation_data_structure_3<Vb,Cb>;
    using Delaunay  = CGAL::Regular_triangulation_3<K,TDS>;
#endif
    using AlphaShape    = CGAL::Alpha_shape_3<Delaunay,CGAL::Tag_true>;

The error I get is:

/usr/local/include/CGAL/Alpha_shape_3.h:1169:24: error: no viable conversion from 'Lazy_alpha_nt_3<[2 * ...], Boolean_tag<true>>' to 'const Lazy_alpha_nt_3<[2 * ...], Boolean_tag<false>>'
    cell_it->set_alpha(alpha);
                       ^~~~~
/usr/local/include/CGAL/Alpha_shape_3.h:314:32: note: in instantiation of member function 'CGAL::Alpha_shape_3<CGAL::Regular_triangulation_3<CGAL::Epick,
      CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Regular_triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<void> >, CGAL::Boolean_tag<true>,
      CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Regular_triangulation_cell_base_3<CGAL::Epick, CGAL::Triangulation_cell_base_3<CGAL::Epick, CGAL::Triangulation_ds_cell_base_3<void>
      >, CGAL::Hidden_points_memory_policy<CGAL::Boolean_tag<true> >, std::__1::list<CGAL::Weighted_point_3<CGAL::Epick>, std::__1::allocator<CGAL::Weighted_point_3<CGAL::Epick> > > >, CGAL::Boolean_tag<true>,
      CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL::Default>, CGAL::Boolean_tag<true> >::initialize_alpha_cell_map' requested here
    if (reinitialize == false) initialize_alpha_cell_map();
                               ^
/usr/local/include/CGAL/Alpha_shape_3.h:263:29: note: in instantiation of member function 'CGAL::Alpha_shape_3<CGAL::Regular_triangulation_3<CGAL::Epick,
      CGAL::Triangulation_data_structure_3<CGAL::Alpha_shape_vertex_base_3<CGAL::Epick, CGAL::Regular_triangulation_vertex_base_3<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_3<void> >, CGAL::Boolean_tag<true>,
      CGAL::Boolean_tag<false> >, CGAL::Alpha_shape_cell_base_3<CGAL::Epick, CGAL::Regular_triangulation_cell_base_3<CGAL::Epick, CGAL::Triangulation_cell_base_3<CGAL::Epick, CGAL::Triangulation_ds_cell_base_3<void>
      >, CGAL::Hidden_points_memory_policy<CGAL::Boolean_tag<true> >, std::__1::list<CGAL::Weighted_point_3<CGAL::Epick>, std::__1::allocator<CGAL::Weighted_point_3<CGAL::Epick> > > >, CGAL::Boolean_tag<true>,
      CGAL::Boolean_tag<false> >, CGAL::Sequential_tag>, CGAL::Default>, CGAL::Boolean_tag<true> >::initialize_alpha' requested here
      if (dimension() == 3) initialize_alpha();
                            ^

[MaelRL]

There is a weighted tag you need to also set in the vertex/cell base class since you're using a regular triangulation.

Note to us: these tag mechanisms really ought to be simplified: default weighted value should look into the triangulation's weighted tag and the exact tag should only be set once.

[mrzv]

Sorry, can you be more specific? I'm not following. This code was working before I added Tag_true. Can you give an example of the modified lines?

[mglisse]

Dmitriy, you could take some inspiration from https://github.com/GUDHI/gudhi-devel/blob/master/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h Don't bother with Tag_true, see #3460.

[mrzv]

@mglisse Thanks for the pointer. That's very helpful. I'm going to borrow liberally.

[mrzv]

@mglisse Is the point of the Gudhi example that if I want correctness, I should switch to Exact_predicates_exact_constructions_kernel? It does the trick on my data.

[mglisse]

It isn't an example, it is a documented interface in Gudhi. It shows how to get exact results, how to get correct approximate results, and how to get fast incorrect results. Yes, Epeck is at the heart of it, but it also shows how to use to_double, what to change for periodic/weighted, etc.

[sloriot]

Note that if you do not change alpha (at least too often) you can should use EPICK + Fixed_alpha_shape_3.

[mglisse]

if you do not change alpha

That's kind of the opposite of what a filtration is :wink:

[sloriot]

I forgot the filtration part indeed :)

[mrzv]

@mglisse I guess I should have provided some context. The code in the gist is boiled down from DioDe, specifically, from diode.hpp. DioDe has the code to do weighted/periodic alpha shape filtrations.

What I was surprised to discover is that filtration() does not guarantee correctness (with Epick). I don't understand why this is. As far as I can tell, the combinatorics of the alpha shape filtration can be determined purely via the predicates. So Epick should be enough to get a correct filtration, but perhaps I'm missing something. On top of that, I can't even guess how simplices go missing. Figuring this out would require diving into CGAL's alpha shapes code, and it would be easier to go back to my old code in Dionysus 1 that constructed alpha shape filtration from CGAL's Delaunay triangulation directly.

Ok, that was an aside. Given that I have working code, what I'd like to understand is the minimum number of changes required to guarantee correctness. Looking at the Gudhi code you linked, the main thing I can tell that's different than what I had is it using either CGAL::Exact_predicates_inexact_constructions_kernel or CGAL::Exact_predicates_exact_constructions_kernel, depending on a template parameter. What I'm trying to understand is whether I'm missing something else, or whether it's really enough to just plug in Epeck, to guarantee correctness?

I understand that it's also showing the correct use of to_double(), and using vertex handles directly, rather than storing points in the map.

[mglisse]

As far as I can tell, the combinatorics of the alpha shape filtration can be determined purely via the predicates.

You would need a predicate that says whether one circumradius is larger than another one (taking up to 8 points as arguments). We don't have such a predicate. If you compute the circumradius and store it, you need to do it exactly or bad things may happen. The Tag_true thing is a lazy way to do this, very similar to Epeck, but local inside Alpha_shape_3. Even if you have the correct combinatorics, you still need to compute some consistent filtration values, they cannot be in a different order than the combinatorial one.

I can't even guess how simplices go missing

No idea, weird things happen when you make wrong assumptions.

go back to my old code in Dionysus 1

Were you doing anything special about numerical issues?

whether it's really enough to just plug in Epeck, to guarantee correctness?

I don't know for sure. If you keep using Epeck::FT as the type of filtration values, it should be fine. If you convert to double using to_double(x.exact()), it should also be fine. If you convert to double with to_double(x), you get approximate values. Maybe there have been enough comparisons so the values have been refined and everything is fine. I think the sorting used to produce the filtration guarantees that. But I would hesitate to promise it. In gudhi we make a pass where for each simplex we increase its filtration value to the max of its faces if necessary, to be safe.

[mrzv]

Hm, in Dionysus 1, I think I followed roughly the following. You can determine if a simplex is critical or regular via the in-sphere test. Once you know which simplices are regular, their alpha value is the minimum of the values of their cofaces. I don't see what numerical issue arises when you take such a minimum, but it must be related to your comment about the 8-point predicate. I understand that to find the minimum exactly you need a predicate, but if all you want is that the value the face gets is not greater than any of its cofaces, I don't see the problem.

I guess the other issue I can imagine is that both a face and a coface are critical, but somehow because of the numerics, the value calculated on the face is greater than the value on the coface. I don't recall if I ever ran into such a problem, but it can be "fudged" using the same minimum of the cofaces, as above. In a sense, it's the right answer, and for persistence calculation (which I'm obviously after), the fudging won't make a difference. Only the combinatorics matters (i.e., that for any alpha value we get a simplicial complex).

Thanks for the rest of it. I've made a change in DioDe (adding an exact argument).

[mrzv]

Just to correct the previous comment, I meant using CGAL::side_of_bounded_sphere() test for criticality.

[mglisse]

[reordered]

Hm, in Dionysus 1, I think I followed roughly the following. You can determine if a simplex is critical or regular via the in-sphere test. Once you know which simplices are regular, their alpha value is the minimum of the values of their cofaces.

That's essentially the same we are doing for alpha-complexes in arbitrary dimensions in gudhi.

I guess the other issue I can imagine is that both a face and a coface are critical, but somehow because of the numerics, the value calculated on the face is greater than the value on the coface.

Yes, that's exactly the problem I have in mind (note that I did not study what problem you hit in CGAL with your example). Filtration values are computed with double, not robustly. In 3D, with the nice formula used, problems are not very common. In dD, with a much worse naive formula, it is numerically quite unstable and I regularly get filtration values in the millions for scenes of diameter 1... which motivated me to write Epeck_d.

I don't recall if I ever ran into such a problem, but it can be "fudged" using the same minimum of the cofaces, as above.

We do the reverse, when we have an inconsistency, we increase the coface instead of decreasing the face (we reuse the code from lower-star filtrations). Maybe the second one would cause slightly less issues, but at some point, if the computation of the filtration values sometimes gives nonsensical values, you can't avoid trouble. You are assuming that the error is small and can be fudged away, but there is no guarantee of that, that's the whole point of using a (lazy) exact kernel.

I don't see what numerical issue arises when you take such a minimum, but it must be related to your comment about the 8-point predicate. I understand that to find the minimum exactly you need a predicate, but if all you want is that the value the face gets is not greater than any of its cofaces, I don't see the problem.

In a sense, it's the right answer, and for persistence calculation (which I'm obviously after), the fudging won't make a difference. Only the combinatorics matters (i.e., that for any alpha value we get a simplicial complex).

CGAL::Alpha_shape_3 wasn't written for persistence, the filtration was added as an afterthought, it isn't particularly tuned for this use. But in any case, I don't see how you can guarantee an ok result without some special care for numerical issues.