diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 2c162e7cc61..49217739942 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -3648,10 +3648,17 @@ incident_slivers(const Vertex_handle& v, std::vector incident_cells_; tr_.incident_cells(v, std::back_inserter(incident_cells_)); +#ifdef CGAL_CXX17 + std::remove_copy_if(incident_cells_.begin(), + incident_cells_.end(), + out, + std::not_fn(Is_sliver(c3t3_, criterion, sliver_bound))); +#else std::remove_copy_if(incident_cells_.begin(), incident_cells_.end(), out, std::not1(Is_sliver(c3t3_,criterion,sliver_bound))); +#endif return out; } 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 d10fa40fe45..7153fd93bfb 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -86,11 +87,12 @@ void update_c3t3_facets(C3t3& c3t3, } } -template +template Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge, C3t3& c3t3, const std::vector& vertices_around_edge, - const Flip_Criterion& criterion) + const Flip_Criterion& criterion, + IncCellsVectorMap& inc_cells) { typedef typename C3t3::Triangulation Tr; typedef typename C3t3::Facet Facet; @@ -186,19 +188,19 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge, if (criterion == MIN_ANGLE_BASED) { //Current worst dihedral angle - FT curr_min_dh = min_dihedral_angle(tr, ch0); - curr_min_dh = (std::min)(curr_min_dh, min_dihedral_angle(tr, ch1)); - curr_min_dh = (std::min)(curr_min_dh, min_dihedral_angle(tr, cell_to_remove)); + Dihedral_angle_cosine curr_max_cosdh = max_cos_dihedral_angle(tr, ch0); + curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, ch1)); + curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, cell_to_remove)); //Result worst dihedral angle - if (curr_min_dh > min_dihedral_angle(tr, vh2, + if (curr_max_cosdh < max_cos_dihedral_angle(tr, vh2, ch0->vertex(indices(vh0_id, 0)), ch0->vertex(indices(vh0_id, 1)), ch0->vertex(indices(vh0_id, 2))) - || curr_min_dh > min_dihedral_angle(tr, vh3, - ch1->vertex(indices(vh1_id, 0)), - ch1->vertex(indices(vh1_id, 1)), - ch1->vertex(indices(vh1_id, 2)))) + || curr_max_cosdh < max_cos_dihedral_angle(tr, vh3, + ch1->vertex(indices(vh1_id, 0)), + ch1->vertex(indices(vh1_id, 1)), + ch1->vertex(indices(vh1_id, 2)))) return NO_BEST_CONFIGURATION; } else if (criterion == AVERAGE_ANGLE_BASED) @@ -214,7 +216,7 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge, (min_dihedral_angle(tr, vh2, ch0->vertex(indices(vh0_id, 0)), ch0->vertex(indices(vh0_id, 1)), ch0->vertex(indices(vh0_id, 2))) - + min_dihedral_angle(tr, vh3, ch1->vertex(indices(vh1_id, 0)), + + min_dihedral_angle(tr, vh3, ch1->vertex(indices(vh1_id, 0)), ch1->vertex(indices(vh1_id, 1)), ch1->vertex(indices(vh1_id, 2)))); //Result worst dihedral angle @@ -312,6 +314,9 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge, ch->set_neighbor(v, mirror_facet.first); } ch->vertex(v)->set_cell(ch); + + inc_cells[ch->vertex(v)] = boost::none; + ch->reset_cache_validity(); } } @@ -372,7 +377,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, typename C3t3::Vertex_handle vh2, typename C3t3::Vertex_handle vh3, CandidatesQueue& candidates, - double curr_min_dh, + const Dihedral_angle_cosine& curr_max_cos_dh, bool is_sliver_well_oriented = true, int e_id = 0) { @@ -382,8 +387,6 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, typedef typename C3t3::Facet Facet; typedef typename Tr::Facet_circulator Facet_circulator; typedef typename Tr::Cell_circulator Cell_circulator; - typedef typename Tr::Geom_traits Gt; - typedef typename Gt::FT FT; // std::cout << "find_best_flip_to_improve_dh boundary " << std::endl; Tr& tr = c3t3.triangulation(); @@ -451,12 +454,12 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, Cell_circulator cell_circulator = tr.incident_cells(edge); Cell_circulator done = cell_circulator; + boost::container::small_vector facets; for (std::size_t i = 0; i < opposite_vertices.size(); ++i) { Vertex_handle vh = opposite_vertices[i]; bool keep = true; - std::vector facets; do { //Store it if it do not have vh @@ -474,7 +477,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, } while (++cell_circulator != done); - FT min_flip_dihedral_angle = (std::numeric_limits::max)(); + Dihedral_angle_cosine max_flip_cos_dh(CGAL::NEGATIVE, 1., 1.); for (const Facet& fi : facets) { if (!tr.is_infinite(fi.first)) @@ -483,32 +486,75 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, fi.first->vertex(indices(fi.second, 1)), fi.first->vertex(indices(fi.second, 2)))) { - min_flip_dihedral_angle = (std::min)(min_flip_dihedral_angle, - min_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)), - fi.first->vertex(indices(fi.second, 1)), - fi.first->vertex(indices(fi.second, 2)))); + max_flip_cos_dh = (std::max)( + max_flip_cos_dh, + max_cos_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)), + fi.first->vertex(indices(fi.second, 1)), + fi.first->vertex(indices(fi.second, 2)))); } else { keep = false; break; } + + if (max_flip_cos_dh.is_one())//it will not get worse than 1. + { + keep = false; + break; + } } } + facets.clear(); - if (keep && (curr_min_dh < min_flip_dihedral_angle || !is_sliver_well_oriented)) + if (keep && (max_flip_cos_dh < curr_max_cos_dh || !is_sliver_well_oriented)) { - //std::cout << "vh " << vh->info() <<" old " << curr_min_dh << " min " << min_flip_dihedral_angle << std::endl; - candidates.push(std::make_pair(min_flip_dihedral_angle, std::make_pair(vh, e_id))); + //std::cout << "vh " << vh->info() <<" old " << curr_max_cos_dh << " min " << min_flip_tan_dh << std::endl; + candidates.push(std::make_pair(max_flip_cos_dh, std::make_pair(vh, e_id))); } } } -template +template +bool is_edge_uv(Vertex_handle u, + Vertex_handle v, + const CellVector& cells_incident_to_u, + Cell_handle& cell, + int& i, + int& j) +{ + if (u == v) + return false; + + for (typename CellVector::value_type c : cells_incident_to_u) + { + if (c->has_vertex(v, j)) + { + cell = c; + i = cell->index(u); + return true; + } + } + return false; +} + +template +bool is_edge_uv(Vertex_handle u, + Vertex_handle v, + const CellVector& cells_incident_to_u) +{ + typename CellVector::value_type c; + int i, j; + return is_edge_uv(u, v, cells_incident_to_u, c, i, j); +} + +template void find_best_flip_to_improve_dh(C3t3& c3t3, typename C3t3::Edge& edge, CandidatesQueue& candidates, - double curr_min_dh, + const Dihedral_angle_cosine& curr_max_cosdh, + IncCellsVectorMap& inc_cells, bool is_sliver_well_oriented = true, int e_id = 0) { @@ -518,8 +564,6 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, typedef typename C3t3::Facet Facet; typedef typename Tr::Facet_circulator Facet_circulator; typedef typename Tr::Cell_circulator Cell_circulator; - typedef typename Tr::Geom_traits Gt; - typedef typename Gt::FT FT; Tr& tr = c3t3.triangulation(); @@ -547,6 +591,17 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, } } + if(tr.is_infinite(vh)) + continue; + + boost::optional>& o_inc_vh = inc_cells[vh]; + if (o_inc_vh == boost::none) + { + boost::container::small_vector inc_vec; + tr.incident_cells(vh, std::back_inserter(inc_vec)); + o_inc_vh = inc_vec; + } + Facet_circulator facet_circulator = curr_fcirc; Facet_circulator facet_done = curr_fcirc; @@ -563,15 +618,16 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, indices(facet_circulator->second, i)); if (curr_vertex != vh0 && curr_vertex != vh1) { - Cell_handle ch; - int i0, i1; - if (tr.is_edge(curr_vertex, vh, ch, i0, i1)) + if (is_edge_uv(vh, curr_vertex, boost::get(o_inc_vh))) + { is_edge = true; + break; + } } } } while (++facet_circulator != facet_done); - if (!is_edge && !tr.is_infinite(vh)) + if (!is_edge) opposite_vertices.push_back(vh); nb_cells_around_edge++; @@ -589,12 +645,10 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, Cell_circulator cell_circulator = tr.incident_cells(edge); Cell_circulator done = cell_circulator; - for (std::size_t i = 0; i < opposite_vertices.size(); ++i) + boost::container::small_vector facets; + for (Vertex_handle vh : opposite_vertices) { - Vertex_handle vh = opposite_vertices[i]; bool keep = true; - - std::vector facets; do { //Store it if it do not have vh @@ -612,7 +666,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, } while (++cell_circulator != done); - FT min_flip_dihedral_angle = (std::numeric_limits::max)(); + Dihedral_angle_cosine max_flip_cos_dh(CGAL::NEGATIVE, 1., 1.); for (const Facet& fi : facets) { if (!tr.is_infinite(fi.first)) @@ -621,31 +675,39 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, fi.first->vertex(indices(fi.second, 1)), fi.first->vertex(indices(fi.second, 2)))) { - min_flip_dihedral_angle = (std::min)(min_flip_dihedral_angle, - min_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)), - fi.first->vertex(indices(fi.second, 1)), - fi.first->vertex(indices(fi.second, 2)))); + max_flip_cos_dh = (std::max)(max_flip_cos_dh, + max_cos_dihedral_angle(tr, vh, fi.first->vertex(indices(fi.second, 0)), + fi.first->vertex(indices(fi.second, 1)), + fi.first->vertex(indices(fi.second, 2)))); } else { keep = false; break; } + + if (max_flip_cos_dh.is_one())//it will not get worse than 1. + { + keep = false; + break; + } } } + facets.clear(); - if (keep && (curr_min_dh < min_flip_dihedral_angle || !is_sliver_well_oriented)) + if (keep && (max_flip_cos_dh < curr_max_cosdh || !is_sliver_well_oriented)) { - //std::cout << "vh " << vh->info() <<" old " << curr_min_dh << " min " << min_flip_dihedral_angle << std::endl; - candidates.push(std::make_pair(min_flip_dihedral_angle, std::make_pair(vh, e_id))); + //std::cout << "vh " << vh->info() <<" old " << curr_max_cosdh << " min " << min_flip_tan_dh << std::endl; + candidates.push(std::make_pair(max_flip_cos_dh, std::make_pair(vh, e_id))); } } } -template +template Sliver_removal_result flip_n_to_m(C3t3& c3t3, typename C3t3::Edge& edge, typename C3t3::Vertex_handle vh, + IncCellsVectorMap& inc_cells, Visitor& visitor, bool check_validity = false) { @@ -693,7 +755,14 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, facet_circulator++; facet_circulator++; - std::vector vertices_around_edge; + boost::optional>& o_inc_vh = inc_cells[vh]; + if (o_inc_vh == boost::none) + { + boost::container::small_vector inc_vec; + tr.incident_cells(vh, std::back_inserter(inc_vec)); + o_inc_vh = inc_vec; + } + do { //Get the ids of the opposite vertices @@ -703,18 +772,13 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, indices(facet_circulator->second, i)); if (curr_vertex != vh0 && curr_vertex != vh1) { - Cell_handle ch; - int i0, i1; - if (tr.is_edge(curr_vertex, vh, ch, i0, i1)) + if (is_edge_uv(vh, curr_vertex, boost::get(o_inc_vh))) return NOT_FLIPPABLE; - - vertices_around_edge.push_back(curr_vertex); } } } while (++facet_circulator != facet_done); - std::vector cells_around_edge; boost::container::small_vector to_remove; //Neighbors that will need to be updated after flip @@ -733,8 +797,6 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, Cell_circulator done = cell_circulator; do { - cells_around_edge.push_back(cell_circulator); - //Facets opposite to vh0 Facet facet_vh0(cell_circulator, cell_circulator->index(vh0)); neighbor_facets.insert(tr.mirror_facet(facet_vh0)); @@ -877,6 +939,9 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, ch->set_neighbor(v, facet.first); } ch->vertex(v)->set_cell(ch); + + inc_cells[ch->vertex(v)] = boost::none; + ch->reset_cache_validity(); } } @@ -943,54 +1008,61 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, } -template +template Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge, C3t3& c3t3, - std::vector& boundary_vertices, + const std::vector& boundary_vertices, const Flip_Criterion& criterion, + IncCellsVectorMap& inc_cells, Visitor& visitor) { typedef typename C3t3::Vertex_handle Vertex_handle; typedef typename C3t3::Triangulation::Cell_circulator Cell_circulator; - typedef typename C3t3::Triangulation::Geom_traits Gt; - typedef typename Gt::FT FT; typename C3t3::Triangulation& tr = c3t3.triangulation(); Sliver_removal_result result = NOT_FLIPPABLE; - typedef std::pair > Angle_and_vertex; + typedef std::pair > CosAngle_and_vertex; //std::cout << "n_to_m_flip " << boundary_vertices.size() << std::endl; if (criterion == MIN_ANGLE_BASED) { - std::priority_queue candidates; + std::priority_queue, + std::greater + > candidates; Cell_circulator circ = c3t3.triangulation().incident_cells(edge); Cell_circulator done = circ; - FT curr_min_dh = min_dihedral_angle(tr, circ++); + Dihedral_angle_cosine curr_max_cosdh = max_cos_dihedral_angle(tr, circ++); while (circ != done) { - curr_min_dh = (std::min)(curr_min_dh, min_dihedral_angle(tr, circ++)); + curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, circ++)); } if (boundary_vertices.size() == 2) find_best_flip_to_improve_dh(c3t3, edge, boundary_vertices[0], boundary_vertices[1], - candidates, curr_min_dh); + candidates, curr_max_cosdh); else - find_best_flip_to_improve_dh(c3t3, edge, candidates, curr_min_dh); + find_best_flip_to_improve_dh(c3t3, edge, candidates, curr_max_cosdh, + inc_cells); bool flip_performed = false; while (!candidates.empty() && !flip_performed) { - Angle_and_vertex curr_cost_vpair = candidates.top(); + CosAngle_and_vertex curr_cost_vpair = candidates.top(); candidates.pop(); - //std::cout << curr_min_dh << " old, current " << curr_cost_vpair.second.first->info() <<" and angle " << curr_cost_vpair.first << std::endl; +// std::cout << "\tcurrent cos = " << curr_max_cosdh +// << "\t angle = " << std::acos(curr_max_cosdh) * 180./CGAL_PI << std::endl; +// std::cout << "\tcandidate cos = " << curr_cost_vpair.first +// << "\t angle = " << std::acos(curr_cost_vpair.first) * 180./CGAL_PI << std::endl; +// std::cout << std::endl; - if (curr_min_dh >= curr_cost_vpair.first) + if (curr_max_cosdh <= curr_cost_vpair.first) return NO_BEST_CONFIGURATION; - result = flip_n_to_m(c3t3, edge, curr_cost_vpair.second.first, visitor); + result = flip_n_to_m(c3t3, edge, curr_cost_vpair.second.first, inc_cells, visitor); if (result != NOT_FLIPPABLE) flip_performed = true; @@ -1000,10 +1072,11 @@ Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge, return result; } -template +template Sliver_removal_result find_best_flip(typename C3t3::Edge& edge, C3t3& c3t3, const Flip_Criterion& criterion, + IncCellsVectorMap& inc_cells, Visitor& visitor) { typedef typename C3t3::Triangulation Tr; @@ -1066,7 +1139,7 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge, { std::vector vertices; vertices.insert(vertices.end(), vertices_around_edge.begin(), vertices_around_edge.end()); - res = flip_3_to_2(edge, c3t3, vertices, criterion); + res = flip_3_to_2(edge, c3t3, vertices, criterion, inc_cells); } } else @@ -1078,12 +1151,11 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge, { std::vector vertices; vertices.insert(vertices.end(), boundary_vertices.begin(), boundary_vertices.end()); - res = flip_n_to_m(edge, c3t3, vertices, criterion, visitor); + res = flip_n_to_m(edge, c3t3, vertices, criterion, inc_cells, visitor); //return n_to_m_flip(edge, boundary_vertices, flip_criterion); } } - return res; } @@ -1096,20 +1168,33 @@ std::size_t flip_all_edges(const std::vector& edges, { typedef typename C3t3::Triangulation Tr; typedef typename Tr::Cell_handle Cell_handle; + typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Edge Edge; Tr& tr = c3t3.triangulation(); + std::unordered_map > > inc_cells; + std::size_t count = 0; for (const VertexPair vp : edges) { + boost::optional>& + o_inc_vh = inc_cells[vp.first]; + if (o_inc_vh == boost::none) + { + boost::container::small_vector inc_vec; + tr.incident_cells(vp.first, std::back_inserter(inc_vec)); + o_inc_vh = inc_vec; + } + Cell_handle ch; int i0, i1; - if (tr.is_edge(vp.first, vp.second, ch, i0, i1)) + if (is_edge_uv(vp.first, vp.second, boost::get(o_inc_vh), ch, i0, i1)) { Edge edge(ch, i0, i1); - Sliver_removal_result res = find_best_flip(edge, c3t3, criterion, visitor); + Sliver_removal_result res = find_best_flip(edge, c3t3, criterion, inc_cells, visitor); if (res == INVALID_CELL || res == INVALID_VERTEX || res == INVALID_ORIENTATION) { std::cout << "FLIP PROBLEM!!!!" << std::endl; @@ -1126,6 +1211,10 @@ std::size_t flip_all_edges(const std::vector& edges, } } } + + for(Cell_handle c : c3t3.triangulation().finite_cell_handles()) + c->reset_cache_validity(); + return count; } 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 fc9cf93d1dc..f0e71dfefa1 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h @@ -328,7 +328,12 @@ private: bool check_inversion_and_move(const typename Tr::Vertex_handle v, const typename Tr::Point& final_pos, const CellRange& inc_cells, - const Tr& /* tr */) + const Tr& /* tr */, +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + double& total_move) +#else + double&) +#endif { const typename Tr::Point backup = v->point(); //backup v's position const typename Tr::Geom_traits::Point_3 pv = point(backup); @@ -336,6 +341,7 @@ private: bool valid_orientation = false; double frac = 1.0; typename Tr::Geom_traits::Vector_3 move(pv, point(final_pos)); + do { v->set_point(typename Tr::Point(pv + frac * move)); @@ -354,15 +360,17 @@ private: } } valid_orientation = valid_try; - -// std::cout << std::boolalpha << "valid orientation = " << valid_orientation -// << "\tfrac = " << frac << std::endl; } while(!valid_orientation && frac > 0.1); if (!valid_orientation) //move failed v->set_point(backup); +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + else + total_move += CGAL::approximate_sqrt(CGAL::squared_distance(pv, point(v->point()))); +#endif + return valid_orientation; } @@ -407,7 +415,10 @@ public: std::cout << "Smooth vertices..."; std::cout.flush(); #endif - std::size_t nb_done = 0; + std::size_t nb_done_3d = 0; + std::size_t nb_done_2d = 0; + std::size_t nb_done_1d = 0; + FT total_move = 0.; Tr& tr = c3t3.triangulation(); @@ -521,8 +532,8 @@ public: #endif // move vertex const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z()); - if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr)) - nb_done++; + if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) + nb_done_1d++; } else if (neighbors[vid] > 0) { @@ -555,8 +566,8 @@ public: #endif // move vertex const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z()); - if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr)) - nb_done++; + if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) + nb_done_1d++; } } } @@ -628,8 +639,8 @@ public: os_surf << "2 " << current_pos << " " << final_position << std::endl; #endif const typename Tr::Point new_pos(final_position.x(), final_position.y(), final_position.z()); - if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr)) - nb_done++; + if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) + nb_done_2d++; } else if (neighbors[vid] > 0) { @@ -641,8 +652,8 @@ public: if (boost::optional mls_projection = project(si, current_pos)) { const typename Tr::Point new_pos(CGAL::ORIGIN + *mls_projection); - if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr)) - nb_done++; + if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) + nb_done_2d++; #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG os_surf0 << "2 " << current_pos << " " << new_pos << std::endl; @@ -693,8 +704,8 @@ public: #endif const Vector_3 p = smoothed_positions[vid] / static_cast(neighbors[vid]); typename Tr::Point new_pos(p.x(), p.y(), p.z()); - if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr)) - nb_done++; + if(check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)) + nb_done_3d++; #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG os_vol << " " << point(v->point()) << std::endl; @@ -704,7 +715,11 @@ public: CGAL_assertion(CGAL::Tetrahedral_remeshing::debug::are_cell_orientations_valid(tr)); #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << " done (" << nb_done << " vertices smoothed)." << std::endl; + std::size_t nb_done = nb_done_3d + nb_done_2d + nb_done_1d; + std::cout << " done (" + << nb_done_3d << "/" << nb_done_2d << "/" << nb_done_1d << " vertices smoothed," + << " average move = " << (total_move / nb_done) + << ")." << std::endl; #endif #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h index fb942deddcf..9a76b54a9f8 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h @@ -368,8 +368,18 @@ public: CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "99-postprocess"); #endif #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << "(peeling removed " << nb_slivers_peel << " slivers)" << std::endl; - std::cout << "done." << std::endl; + mindh = 180.; + for (Cell_handle cit : tr().finite_cell_handles()) + { + if (m_c3t3.is_in_complex(cit)) + { + const double dh = min_dihedral_angle(tr(), cit); + mindh = (std::min)(dh, mindh); + } + } + std::cout << "Peeling done (removed " << nb_slivers_peel << " slivers, " + << "min dihedral angle = " << mindh << ")." << std::endl; + #endif return nb_slivers_peel; } 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 0e31c30cb12..cb83de5a777 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 @@ -98,7 +98,7 @@ typename Geom_traits::FT min_dihedral_angle(const Point& p, FT a = CGAL::abs(dihedral_angle(p, q, r, s, gt)); FT min_dh = a; - a = CGAL::abs(dihedral_angle(p, r, q, s, gt)); + a = CGAL::abs(dihedral_angle(p, r, s, q, gt)); min_dh = (std::min)(a, min_dh); a = CGAL::abs(dihedral_angle(p, s, q, r, gt)); @@ -107,7 +107,7 @@ typename Geom_traits::FT min_dihedral_angle(const Point& p, a = CGAL::abs(dihedral_angle(q, r, p, s, gt)); min_dh = (std::min)(a, min_dh); - a = CGAL::abs(dihedral_angle(q, s, p, r, gt)); + a = CGAL::abs(dihedral_angle(q, s, r, p, gt)); min_dh = (std::min)(a, min_dh); a = CGAL::abs(dihedral_angle(r, s, p, q, gt)); @@ -141,6 +141,179 @@ typename Tr::Geom_traits::FT min_dihedral_angle(const Tr& tr, c->vertex(3)); } +struct Dihedral_angle_cosine +{ + CGAL::Sign m_sgn; + double m_sq_num; + double m_sq_den; + + Dihedral_angle_cosine(const CGAL::Sign& sgn, const double& sq_num, const double& sq_den) + : m_sgn(sgn) + , m_sq_num(sq_num) + , m_sq_den(sq_den) + {} + + bool is_one() const + { + return m_sgn == CGAL::POSITIVE && m_sq_num == m_sq_den; + } + double signed_square_value() const + { + switch(m_sgn) + { + case CGAL::POSITIVE: + return m_sq_num / m_sq_den; + case CGAL::ZERO: + return 0.; + default: + CGAL_assertion(m_sgn == CGAL::NEGATIVE); + return -1. * m_sq_num / m_sq_den; + }; + } + + friend bool operator<(const Dihedral_angle_cosine& l, + const Dihedral_angle_cosine& r) + { + //if numerators have different signs + if (l.m_sgn == CGAL::NEGATIVE && r.m_sgn != CGAL::NEGATIVE) + return true; + + else if (l.m_sgn == CGAL::POSITIVE && r.m_sgn != CGAL::POSITIVE) + return false; + + else if (l.m_sgn == CGAL::ZERO) + return (r.m_sgn == CGAL::POSITIVE); + + //else numerators have the same sign + else if (l.m_sgn == CGAL::POSITIVE) //both angles are in [0; PI/2[ + { + CGAL_assertion(r.m_sgn == CGAL::POSITIVE); + + return (l.m_sq_num * r.m_sq_den < r.m_sq_num* l.m_sq_den); + } + else //both angles are in [PI/2; PI] + { + CGAL_assertion(l.m_sgn != CGAL::POSITIVE); + CGAL_assertion(r.m_sgn != CGAL::POSITIVE); + + return (l.m_sq_num * r.m_sq_den >= r.m_sq_num* l.m_sq_den); + } + } + + friend bool operator<=(const Dihedral_angle_cosine& l, + const Dihedral_angle_cosine& r) + { + if(l < r) + return true; + else + return l.m_sgn == r.m_sgn + && l.m_sq_num * r.m_sq_den == r.m_sq_num * l.m_sq_den; + } +}; + +template +Dihedral_angle_cosine cos_dihedral_angle(const typename Gt::Point_3& i, + const typename Gt::Point_3& j, + const typename Gt::Point_3& k, + const typename Gt::Point_3& l, + const Gt& gt) +{ + CGAL_assertion(CGAL::orientation(i, j, k, l) != CGAL::NEGATIVE); + + typename Gt::Construct_vector_3 vector = gt.construct_vector_3_object(); + typename Gt::Construct_cross_product_vector_3 cross_product = + gt.construct_cross_product_vector_3_object(); + typename Gt::Compute_scalar_product_3 scalar_product = + gt.compute_scalar_product_3_object(); + + typedef typename Gt::FT FT; + typedef typename Gt::Vector_3 Vector_3; + const Vector_3 ij = vector(i, j); + const Vector_3 ik = vector(i, k); + const Vector_3 il = vector(i, l); + + const Vector_3 ijik = cross_product(ij, ik); + if(CGAL::NULL_VECTOR == ijik) + return Dihedral_angle_cosine(CGAL::POSITIVE, 1.,1.); + + const Vector_3 ilij = cross_product(il, ij); + if (CGAL::NULL_VECTOR == ilij) + return Dihedral_angle_cosine(CGAL::POSITIVE, 1.,1.); + + const FT num = scalar_product(ijik, ilij); + if(num == 0.) + return Dihedral_angle_cosine(CGAL::ZERO, 0.,1.); + + const double sqden = CGAL::to_double( + scalar_product(ijik, ijik) * scalar_product(ilij, ilij)); + + return Dihedral_angle_cosine(CGAL::sign(num), CGAL::square(num), sqden); +} + +template +Dihedral_angle_cosine max_cos_dihedral_angle(const Point& p, + const Point& q, + const Point& r, + const Point& s, + const Geom_traits& gt) +{ + Dihedral_angle_cosine a = cos_dihedral_angle(p, q, r, s, gt); + Dihedral_angle_cosine max_cos_dh = a; + if(max_cos_dh.is_one()) return max_cos_dh; + + a = cos_dihedral_angle(p, r, s, q, gt); + if(max_cos_dh < a) max_cos_dh = a; + if (max_cos_dh.is_one()) return max_cos_dh; + + a = cos_dihedral_angle(p, s, q, r, gt); + if (max_cos_dh < a) max_cos_dh = a; + if (max_cos_dh.is_one()) return max_cos_dh; + + a = cos_dihedral_angle(q, r, p, s, gt); + if (max_cos_dh < a) max_cos_dh = a; + if (max_cos_dh.is_one()) return max_cos_dh; + + a = cos_dihedral_angle(q, s, r, p, gt); + if (max_cos_dh < a) max_cos_dh = a; + if (max_cos_dh.is_one()) return max_cos_dh; + + a = cos_dihedral_angle(r, s, p, q, gt); + if (max_cos_dh < a) max_cos_dh = a; + + return max_cos_dh; +} + +template +Dihedral_angle_cosine max_cos_dihedral_angle(const Tr& tr, + const typename Tr::Vertex_handle v0, + const typename Tr::Vertex_handle v1, + const typename Tr::Vertex_handle v2, + const typename Tr::Vertex_handle v3) +{ + return max_cos_dihedral_angle(point(v0->point()), + point(v1->point()), + point(v2->point()), + point(v3->point()), + tr.geom_traits()); +} + +template +Dihedral_angle_cosine max_cos_dihedral_angle(const Tr& tr, + const typename Tr::Cell_handle c) +{ + if (c->is_cache_valid()) + return Dihedral_angle_cosine(CGAL::sign(c->sliver_value()), + CGAL::abs(c->sliver_value()), 1.); + + Dihedral_angle_cosine cos_dh = max_cos_dihedral_angle(tr, + c->vertex(0), + c->vertex(1), + c->vertex(2), + c->vertex(3)); + c->set_sliver_value(cos_dh.signed_square_value()); + return cos_dh; +} + template bool is_peelable(const C3t3& c3t3, const typename C3t3::Cell_handle ch,