From 8076e20b71bc83407ec94d6e399d6e41f163c8de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 1 Aug 2023 13:07:55 +0200 Subject: [PATCH 1/6] Add debug code --- .../CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h | 14 ++++++++------ .../test_AW3_cavity_initializations.cpp | 1 + .../test/Alpha_wrap_3/test_AW3_manifoldness.cpp | 1 + .../test/Alpha_wrap_3/test_AW3_multiple_calls.cpp | 1 + .../test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp | 2 +- 5 files changed, 12 insertions(+), 7 deletions(-) diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index e88e00fafe0..1e9674276fa 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -1362,7 +1362,12 @@ private: for(Vertex_handle v : m_dt.finite_vertex_handles()) { if(is_non_manifold(v)) + { +#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS_PP + std::cout << v->point() << " is non-manifold" << std::endl; +#endif non_manifold_vertices.push(v); + } } // Some lambdas for the comparer @@ -1412,23 +1417,19 @@ private: squared_distance(m_dt.point(c, 2), m_dt.point(c, 3)) }); }; -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS +#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS_PP std::cout << non_manifold_vertices.size() << " initial NMV" << std::endl; #endif while(!non_manifold_vertices.empty()) { -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS +#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS_PP std::cout << non_manifold_vertices.size() << " NMV in queue" << std::endl; #endif Vertex_handle v = non_manifold_vertices.top(); non_manifold_vertices.pop(); -#ifdef CGAL_AW3_DEBUG_MANIFOLDNESS - std::cout << "ยท"; -#endif - if(!is_non_manifold(v)) continue; @@ -1503,6 +1504,7 @@ private: void check_queue_sanity() { std::cout << "Check queue sanity..." << std::endl; + std::vector queue_gates; Gate previous_top_gate = m_queue.top(); while(!m_queue.empty()) diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp index d7567bab205..c6591e000df 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_cavity_initializations.cpp @@ -1,5 +1,6 @@ #define CGAL_AW3_TIMER #define CGAL_AW3_DEBUG +#define CGAL_AW3_DEBUG_MANIFOLDNESS // #define CGAL_AW3_DEBUG_INITIALIZATION #include diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp index 928cade6495..8ac2e422d87 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_manifoldness.cpp @@ -1,5 +1,6 @@ #define CGAL_AW3_TIMER #define CGAL_AW3_DEBUG +#define CGAL_AW3_DEBUG_MANIFOLDNESS //#define CGAL_AW3_DEBUG_STEINER_COMPUTATION //#define CGAL_AW3_DEBUG_INITIALIZATION //#define CGAL_AW3_DEBUG_QUEUE diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp index e2abc6f1f12..eda64a8376a 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_AW3_multiple_calls.cpp @@ -1,5 +1,6 @@ #define CGAL_AW3_TIMER #define CGAL_AW3_DEBUG +#define CGAL_AW3_DEBUG_MANIFOLDNESS #include diff --git a/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp b/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp index 5091934b62f..302494d8c30 100644 --- a/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp +++ b/Alpha_wrap_3/test/Alpha_wrap_3/test_alpha_wrap_3_mesh.cpp @@ -1,6 +1,6 @@ #define CGAL_AW3_TIMER //#define CGAL_AW3_DEBUG -//#define CGAL_AW3_DEBUG_MANIFOLDNESS +#define CGAL_AW3_DEBUG_MANIFOLDNESS //#define CGAL_AW3_DEBUG_STEINER_COMPUTATION //#define CGAL_AW3_DEBUG_INITIALIZATION //#define CGAL_AW3_DEBUG_QUEUE From c7b9317a96117f874aa88b6ca9cd6b376f4d10cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Tue, 1 Aug 2023 13:16:18 +0200 Subject: [PATCH 2/6] Simplify choice of cells to un-carve while enforcing manifoldness This combinatorial choice seemed like a good idea, but it can have nasty cascading effects, adding very large tetrahedra. See this issue: https://github.com/CGAL/cgal/issues/7625 In the end, the only thing we care about is small volumes being added. I keep the artificial vertex for now, but I am not fully convinced these should be actually kept too. --- .../include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h | 5 ----- 1 file changed, 5 deletions(-) diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index 1e9674276fa..1b0a6968c09 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -1445,11 +1445,6 @@ private: if(has_artificial_vertex(r)) return true; - const int l_bf_count = count_boundary_facets(l, v); - const int r_bf_count = count_boundary_facets(r, v); - if(l_bf_count != r_bf_count) - return l_bf_count > r_bf_count; - return sq_longest_edge(l) < sq_longest_edge(r); }; From b4e207ab00237bffbd221bc7f4bd4b0bb2e69c47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Wed, 2 Aug 2023 10:00:34 +0200 Subject: [PATCH 3/6] Add some comments on AW3 manifoldness heuristics criteria --- .../CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h | 43 +++++++++++-------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index 1b0a6968c09..a21925cdc6d 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -1380,25 +1380,27 @@ private: return false; }; - auto is_on_boundary = [](Cell_handle c, int i) -> bool - { - return (c->info().is_outside != c->neighbor(i)->info().is_outside); - }; - - auto count_boundary_facets = [&](Cell_handle c, Vertex_handle v) -> int - { - const int v_index_in_c = c->index(v); - int boundary_facets = 0; - for(int i=0; i<3; ++i) // also take into account the opposite facet? - { - if(i == v_index_in_c) - continue; - - boundary_facets += is_on_boundary(c, i); - } - - return boundary_facets; - }; + // This seemed like a good idea, but in the end it can have strong cascading issues, + // whereas some cells with much lower volume would have solved the non-manifoldness. +// auto is_on_boundary = [](Cell_handle c, int i) -> bool +// { +// return (c->info().is_outside != c->neighbor(i)->info().is_outside); +// }; +// +// auto count_boundary_facets = [&](Cell_handle c, Vertex_handle v) -> int +// { +// const int v_index_in_c = c->index(v); +// int boundary_facets = 0; +// for(int i=0; i<3; ++i) // also take into account the opposite facet? +// { +// if(i == v_index_in_c) +// continue; +// +// boundary_facets += is_on_boundary(c, i); +// } +// +// return boundary_facets; +// }; // longest edge works better // auto sq_circumradius = [&](Cell_handle c) -> FT @@ -1407,6 +1409,9 @@ private: // return geom_traits().compute_squared_distance_3_object()(m_dt.point(c, 0), cc); // }; + // the reasonning behind using longest edge rather than volume is that we want to avoid + // spikes (which would have a small volume), and can often happen since we do not spend + // any care on the quality of tetrahedra. auto sq_longest_edge = [&](Cell_handle c) -> FT { return (std::max)({ squared_distance(m_dt.point(c, 0), m_dt.point(c, 1)), From 330ff2e65776a9be1f7ca757a813e5055e0a1994 Mon Sep 17 00:00:00 2001 From: Mael Date: Wed, 2 Aug 2023 10:23:12 +0200 Subject: [PATCH 4/6] Fix spelling Thanks @albert-github! --- Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index a21925cdc6d..6d544be058d 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -1409,7 +1409,7 @@ private: // return geom_traits().compute_squared_distance_3_object()(m_dt.point(c, 0), cc); // }; - // the reasonning behind using longest edge rather than volume is that we want to avoid + // the reasoning behind using longest edge rather than volume is that we want to avoid // spikes (which would have a small volume), and can often happen since we do not spend // any care on the quality of tetrahedra. auto sq_longest_edge = [&](Cell_handle c) -> FT From bc8351f15678b90955a9d5c2a5cc4b0389b49642 Mon Sep 17 00:00:00 2001 From: Mael Date: Wed, 27 Sep 2023 11:06:49 +0200 Subject: [PATCH 5/6] Fix typo --- Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index 6d544be058d..639c9a3f468 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -1418,7 +1418,7 @@ private: squared_distance(m_dt.point(c, 0), m_dt.point(c, 2)), squared_distance(m_dt.point(c, 0), m_dt.point(c, 3)), squared_distance(m_dt.point(c, 1), m_dt.point(c, 2)), - squared_distance(m_dt.point(c, 3), m_dt.point(c, 3)), + squared_distance(m_dt.point(c, 1), m_dt.point(c, 3)), squared_distance(m_dt.point(c, 2), m_dt.point(c, 3)) }); }; From 7de4f442e8ff9067d0ea3637fa504d9528d5bfcc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Mon, 9 Oct 2023 12:12:31 +0200 Subject: [PATCH 6/6] Remove obsolete sort at every iteration There was a need for sorting at every iteration when the sorting used criteria which were changing with every iteration. This is no longer the case after c7b9317. Also make it deterministic. --- .../CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h index 639c9a3f468..ba4635889b4 100644 --- a/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h +++ b/Alpha_wrap_3/include/CGAL/Alpha_wrap_3/internal/Alpha_wrap_3.h @@ -1440,15 +1440,17 @@ private: // Prioritize: // - cells without bbox vertices - // - cells that already have a large number of boundary facets // - small cells when equal number of boundary facets - // @todo give topmost priority to cells with > 1 non-manifold vertex? + // + // Note that these are properties that do not depend on where the surface is, so we can + // sort once. However, if a criterion such as the number of inside cells were added, + // one would need to sort again after each modification of is_outside() statuses. auto comparer = [&](Cell_handle l, Cell_handle r) -> bool { - if(has_artificial_vertex(l)) - return false; - if(has_artificial_vertex(r)) - return true; + CGAL_precondition(!m_dt.is_infinite(l) && !m_dt.is_infinite(r)); + + if(has_artificial_vertex(l) != has_artificial_vertex(r)) + return has_artificial_vertex(r); return sq_longest_edge(l) < sq_longest_edge(r); }; @@ -1457,17 +1459,13 @@ private: inc_cells.reserve(64); m_dt.finite_incident_cells(v, std::back_inserter(inc_cells)); -#define CGAL_AW3_USE_BRUTE_FORCE_MUTABLE_PRIORITY_QUEUE -#ifndef CGAL_AW3_USE_BRUTE_FORCE_MUTABLE_PRIORITY_QUEUE - std::sort(inc_cells.begin(), inc_cells.end(), comparer); // sort once -#endif + // 'std::stable_sort' to have determinism without having to write something like: + // if(longest_edge(l) == longest_edge(r)) return ... + // in the comparer. It's a small range, so the cost does not matter. + std::stable_sort(inc_cells.begin(), inc_cells.end(), comparer); for(auto cit=inc_cells.begin(), cend=inc_cells.end(); cit!=cend; ++cit) { -#ifdef CGAL_AW3_USE_BRUTE_FORCE_MUTABLE_PRIORITY_QUEUE - // sort at every iteration since the number of boundary facets evolves - std::sort(cit, cend, comparer); -#endif Cell_handle ic = *cit; CGAL_assertion(!m_dt.is_infinite(ic));