diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h index 831e0e94eaf..54d950e0072 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h @@ -477,7 +477,8 @@ template bool collapse_keeps_orientations(const CellRange& star, const typename C3t3::Vertex_handle v_moved, const typename C3t3::Vertex_handle v_other, - const typename C3t3::Triangulation::Geom_traits::Point_3& new_pos) + const typename C3t3::Triangulation::Geom_traits::Point_3& new_pos, + const C3t3& c3t3) { typedef typename C3t3::Triangulation::Geom_traits::Point_3 Point; @@ -490,6 +491,11 @@ bool collapse_keeps_orientations(const CellRange& star, point(ch->vertex(1)->point()), point(ch->vertex(2)->point()), point(ch->vertex(3)->point()) }; + + if(c3t3.triangulation().may_have_badly_oriented_cells() && + CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE) + return true; //in this case we don't care about orientation change, locally + pts[ch->index(v_moved)] = new_pos; if (CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE) return false; @@ -497,102 +503,6 @@ bool collapse_keeps_orientations(const CellRange& star, return true; } -template -bool is_valid_collapse(const typename C3t3::Edge& edge, - const Collapse_type& collapse_type, - const typename C3t3::Triangulation::Point& new_pos, - const C3t3& c3t3) -{ - typedef typename C3t3::Vertex_handle Vertex_handle; - typedef typename C3t3::Cell_handle Cell_handle; - typedef typename C3t3::Triangulation::Point Point; - - const Vertex_handle v0 = edge.first->vertex(edge.second); - const Vertex_handle v1 = edge.first->vertex(edge.third); - -#ifdef CGAL_DEBUG_TET_REMESHING_IN_PLUGIN - const bool in_cx = c3t3.is_in_complex(edge); - if (in_cx) - { - if (collapse_type == TO_MIDPOINT) - nb_test_midpoint++; - else if (collapse_type == TO_V1) - nb_test_v1++; - else - nb_test_v0++; - } -#endif - - if (collapse_type == TO_V1 || collapse_type == TO_MIDPOINT) - { - std::vector cells_to_check; - c3t3.triangulation().finite_incident_cells(v0, - std::back_inserter(cells_to_check)); - - for (const Cell_handle& ch : cells_to_check) - { - if (!ch->has_vertex(v1)) - { - //check orientation - std::array pts = { ch->vertex(0)->point(), - ch->vertex(1)->point(), - ch->vertex(2)->point(), - ch->vertex(3)->point()}; - pts[ch->index(v0)] = new_pos; - if (CGAL::orientation(point(pts[0]), point(pts[1]), point(pts[2]), point(pts[3])) - != CGAL::POSITIVE) - { -#ifdef CGAL_DEBUG_TET_REMESHING_IN_PLUGIN - if (in_cx) - { - if (collapse_type == TO_MIDPOINT) - nb_orientation_midpoint++; - else - nb_orientation_v1++; - } -#endif - return false; - } - } - } - } - if (collapse_type == TO_V0 || collapse_type == TO_MIDPOINT) - { - std::vector cells_to_check; - c3t3.triangulation().finite_incident_cells(v1, - std::back_inserter(cells_to_check)); - - for (const Cell_handle& ch : cells_to_check) - { - if (!ch->has_vertex(v0)) - { - //check orientation - std::array pts = { ch->vertex(0)->point(), - ch->vertex(1)->point(), - ch->vertex(2)->point(), - ch->vertex(3)->point() }; - pts[ch->index(v1)] = new_pos; - if (CGAL::orientation(point(pts[0]), point(pts[1]), point(pts[2]), point(pts[3])) - != CGAL::POSITIVE) - { -#ifdef CGAL_DEBUG_TET_REMESHING_IN_PLUGIN - if (in_cx) - { - if (collapse_type == TO_MIDPOINT) - nb_orientation_midpoint++; - else - nb_orientation_v0++; - } -#endif - return false; - } - } - } - } - - return is_valid_collapse(edge, c3t3); -} - template bool facet_has_edge(const Facet& f, const Vh v0, const Vh v1) { @@ -1306,10 +1216,10 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge, const auto orientations_ok = [&](const Collapse_type ct, const Point& pos) { if ((ct == TO_V1 || ct == TO_MIDPOINT) - && !collapse_keeps_orientations(star_of_v0(), v0, v1, point(pos))) + && !collapse_keeps_orientations(star_of_v0(), v0, v1, point(pos), c3t3)) return false; if ((ct == TO_V0 || ct == TO_MIDPOINT) - && !collapse_keeps_orientations(star_of_v1(), v1, v0, point(pos))) + && !collapse_keeps_orientations(star_of_v1(), v1, v0, point(pos), c3t3)) return false; return true; }; diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h index ed33ce8a714..d79ef82abaf 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h @@ -1988,7 +1988,8 @@ void flip_edges(C3T3& c3t3, vertices_subdomain_indices); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - if(!debug::are_cell_orientations_valid(c3t3.triangulation())) + if(!c3t3.triangulation().may_have_badly_oriented_cells() + && !debug::are_cell_orientations_valid(c3t3.triangulation())) std::cerr << "ERROR in ORIENTATION" << std::endl; #endif diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h index 47512f50f59..adf26d26e79 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h @@ -496,8 +496,6 @@ class Tetrahedral_remeshing_smoother const typename Tr::Point backup = v->point(); //backup v's position const typename Tr::Geom_traits::Point_3 pv = point(backup); - bool valid_orientation = false; - bool angles_improved = true; double frac = 1.0; typename Tr::Geom_traits::Vector_3 move(pv, final_pos); @@ -505,7 +503,11 @@ class Tetrahedral_remeshing_smoother ? max_cosine(tr, inc_cells) : Dihedral_angle_cosine(CGAL::ZERO, 0., 1.);//Dummy unused value + const bool check_orientation = need_to_check_orientation_after_change(tr, inc_cells); + bool valid_try = true; + bool valid_orientation = true; + bool angles_improved = true; do { v->set_point(typename Tr::Point(pv + frac * move)); @@ -516,10 +518,7 @@ class Tetrahedral_remeshing_smoother for (const typename Tr::Cell_handle& ci : inc_cells) { - if (CGAL::POSITIVE != CGAL::orientation(point(ci->vertex(0)->point()), - point(ci->vertex(1)->point()), - point(ci->vertex(2)->point()), - point(ci->vertex(3)->point()))) + if(check_orientation && !is_well_oriented(tr, ci)) { frac = 0.5 * frac; valid_try = false; diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h index 8a53ed2ea9e..c9e360bce51 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h @@ -40,10 +40,10 @@ namespace internal { template -bool positive_orientation_after_edge_split(const typename C3t3::Edge& e, +bool valid_orientation_after_edge_split(const typename C3t3::Edge& e, const typename C3t3::Cell_handle circ, const typename C3t3::Triangulation::Geom_traits::Point_3& steiner, - const C3t3&) + const C3t3& c3t3) { using Point = typename C3t3::Triangulation::Geom_traits::Point_3; @@ -54,6 +54,13 @@ bool positive_orientation_after_edge_split(const typename C3t3::Edge& e, point(circ->vertex(1)->point()), point(circ->vertex(2)->point()), point(circ->vertex(3)->point())}; + + // if the orientation is already bad before split, do not bother checking + // orientation after split, and return true straight away + if(c3t3.triangulation().may_have_badly_oriented_cells() + && CGAL::orientation(pts[0], pts[1], pts[2], pts[3]) != CGAL::POSITIVE) + return true; + // 1st half-cell const int i1 = circ->index(v1); const Point p1 = pts[i1]; @@ -103,7 +110,7 @@ construct_steiner_point(const typename C3t3::Edge& e, do { Cell_handle c = circ; - if(!positive_orientation_after_edge_split(e, c, steiner, c3t3)) + if(!valid_orientation_after_edge_split(e, c, steiner, c3t3)) { steiner_successful = false; break; @@ -184,7 +191,7 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e, } const Cell_handle c = circ; - if(!positive_orientation_after_edge_split(e, c, m, c3t3)) + if(!valid_orientation_after_edge_split(e, c, m, c3t3)) { const std::optional steiner = construct_steiner_point(e, c3t3); if (steiner != std::nullopt) diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h index 79eada1496c..edb6622be33 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h @@ -287,8 +287,7 @@ Dihedral_angle_cosine cos_dihedral_angle(const typename Gt::Point_3& i, const typename Gt::Point_3& l, const Gt& gt) { - CGAL_expensive_assertion(CGAL::orientation(i, j, k, l) != CGAL::NEGATIVE); - + //valid however orientation(i,j,k,l) is positive or negative typename Gt::Construct_cross_product_vector_3 cross_product = gt.construct_cross_product_vector_3_object(); typename Gt::Compute_scalar_product_3 scalar_product = @@ -365,7 +364,8 @@ Dihedral_angle_cosine max_cos_dihedral_angle(const Point& p, const Vector_3 ps = vector(p, s); const Vector_3 pr = vector(p, r); - //compute normals pointing outside tetrahedron + // compute normals pointing outside tetrahedron if orientation(p,q,r,s) is POSITIVE + // and inside tetrahedron otherwise const Vector_3 n_pqr = cross(qp, qr); if (CGAL::NULL_VECTOR == n_pqr) return Dihedral_angle_cosine(CGAL::POSITIVE, 1., 1.); @@ -698,6 +698,28 @@ bool is_well_oriented(const Tr& tr, const typename Tr::Cell_handle ch) ch->vertex(3)); } +template +bool need_to_check_orientation_after_change(const Tr& tr, const CellRange& cells) +{ + for (const auto& cell : cells) + { + if (!is_well_oriented(tr, cell)) + { + // if a cell is badly oriented, but the triangulation is known + // to have badly oriented cells, then no need to re-check orientation + // after change. It will remain bad + if(tr.may_have_badly_oriented_cells()) + return false; + else + CGAL_assertion(false);//this function is called before the change, + //so this point should not be reached + } + } + + // all cells are positively oriented + return true; +} + template bool is_boundary(const C3T3& c3t3, const typename C3T3::Facet& f, @@ -2231,6 +2253,11 @@ void count_far_points(const C3t3& c3t3) template bool are_cell_orientations_valid(const Tr& tr) { + // skip the test if we know + // that the input triangulation already has inverted cells + if(tr.may_have_badly_oriented_cells()) + return true; + typedef typename Tr::Geom_traits::Point_3 Point_3; typedef typename Tr::Facet Facet; diff --git a/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_from_mesh_file.cpp b/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_from_mesh_file.cpp index 30c83e1c533..45fa0c39d21 100644 --- a/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_from_mesh_file.cpp +++ b/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_from_mesh_file.cpp @@ -7,26 +7,96 @@ #include #include +#include +#include + #include #include -#include - typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3 Remeshing_triangulation; +typedef Remeshing_triangulation Tr; -int main(int argc, char* argv[]) + +int test(const std::string& filename, const bool perturb_vertices_positions = false) { - const double target_edge_length = (argc > 1) ? atof(argv[1]) : 0.1; + const double target_edge_length = 0.05; - Remeshing_triangulation tr; + std::ifstream in(filename); - std::ifstream in("data/sphere.mesh"); - if (CGAL::SMDS_3::build_triangulation_from_file(in, tr)) + Remeshing_triangulation tr; + if(CGAL::SMDS_3::build_triangulation_from_file(in, tr, + false /*verbose*/, + false /*replace_domain_0*/, + false /*allow_non_manifold*/, + true /*allow_negative_orientation*/)) std::cout << "build triangulation ok" << std::endl; + else + return EXIT_FAILURE; + + if (perturb_vertices_positions) + { + std::cout << "CGAL Random seed = " << CGAL::get_default_random().get_seed() << std::endl; + using P = Remeshing_triangulation::Point; + const int percentage = 10;// 1/percentage = fraction of perturbed vertices + const double max_move = 0.1 * target_edge_length; + + CGAL::Random rng; + int vid = 0; + for (auto v : tr.finite_vertex_handles()) + { + if(++vid % percentage != 0) + continue; + + const P& pv = v->point(); + v->set_point(P{pv.x() + rng.get_double(-1., 1.) * max_move, + pv.y() + rng.get_double(-1., 1.) * max_move, + pv.z() + rng.get_double(-1., 1.) * max_move}); + } + } + + auto count_negative_cells = [&](const Tr& tr) -> unsigned int + { + return std::count_if(tr.finite_cells_begin(), + tr.finite_cells_end(), + [](const auto& c) + { + return CGAL::POSITIVE != CGAL::orientation(c.vertex(0)->point(), + c.vertex(1)->point(), + c.vertex(2)->point(), + c.vertex(3)->point()); + }); + }; + + unsigned int neg = count_negative_cells(tr); + if(neg > 0) + tr.may_have_badly_oriented_cells(true); + std::cout << neg << " finite cells are badly oriented before remeshing" << std::endl; + + CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length, + CGAL::parameters::number_of_iterations(5) + .remesh_boundaries(false)); + + unsigned int neg2 = count_negative_cells(tr); + if(neg2 > 0) + std::cout << neg2 << " finite cells are badly oriented after remeshing" << std::endl; + + in.close(); + return EXIT_SUCCESS; +} + +int main(int argc, char* argv[]) +{ + std::cout << "Run remeshing from elephant.mesh" << std::endl; + const std::string file1 = CGAL::data_file_path("meshes/elephant.mesh"); + if(test(file1) != EXIT_SUCCESS) + return EXIT_FAILURE; - CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length); + std::cout << "\n\nRun remeshing from perturbed elephant.mesh" << std::endl; + const std::string file2 = CGAL::data_file_path("meshes/elephant.mesh"); + if(test(file2, true) != EXIT_SUCCESS) + return EXIT_FAILURE; return EXIT_SUCCESS; }