1 // Copyright (c) 2020 GeometryFactory (France) and Telecom Paris (France). 2 // All rights reserved. 3 // 4 // This file is part of CGAL (www.cgal.org) 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/compute_c3t3_statistics.h $ 7 // $Id: compute_c3t3_statistics.h ab05dde 2020-06-12T08:08:56+02:00 Jane Tournois 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // 11 // Author(s) : Jane Tournois, Noura Faraj, Jean-Marc Thiery, Tamy Boubekeur 12 13 #ifndef CGAL_TR_INTERNAL_COMPUTE_C3T3_STATISTICS_H 14 #define CGAL_TR_INTERNAL_COMPUTE_C3T3_STATISTICS_H 15 16 #include <CGAL/license/Tetrahedral_remeshing.h> 17 18 #include <limits> 19 #include <vector> 20 #include <algorithm> 21 #include <fstream> 22 23 #include <boost/unordered_set.hpp> 24 25 #include <CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h> 26 27 namespace CGAL 28 { 29 namespace Tetrahedral_remeshing 30 { 31 namespace internal 32 { 33 template<typename Triangulation, typename CellSelector> 34 void compute_statistics(const Triangulation& tr, 35 CellSelector cell_selector, 36 const char* filename = "statistics_c3t3.txt") 37 { 38 typedef Triangulation Tr; 39 typedef typename Tr::Geom_traits Gt; 40 typedef typename Tr::Cell_handle Cell_handle; 41 typedef typename Tr::Vertex_handle Vertex_handle; 42 typedef typename Gt::Point_3 Point; 43 typedef typename Tr::Finite_facets_iterator Finite_facets_iterator; 44 typedef typename Tr::Finite_cells_iterator Finite_cells_iterator; 45 typedef typename Tr::Cell::Subdomain_index Subdomain_index; 46 47 std::size_t nb_edges = 0; 48 double total_edges = 0; 49 std::size_t nb_angle = 0; 50 double total_angle = 0; 51 52 double min_edges_length = (std::numeric_limits<double>::max)(); 53 double max_edges_length = 0.; 54 55 double smallest_edge_radius = (std::numeric_limits<double>::max)(); 56 double smallest_radius_radius = (std::numeric_limits<double>::max)(); 57 double biggest_v_sma_cube = 0.; 58 double max_dihedral_angle = 0.; 59 double min_dihedral_angle = 180.; 60 61 for (Finite_facets_iterator fit = tr.finite_facets_begin(); 62 fit != tr.finite_facets_end(); ++fit) 63 { 64 const Cell_handle cell = fit->first; 65 const int& index = fit->second; 66 if (!cell_selector(cell) || !cell_selector(cell->neighbor(index))) 67 continue; 68 69 const Point& pa = point(cell->vertex((index + 1) & 3)->point()); 70 const Point& pb = point(cell->vertex((index + 2) & 3)->point()); 71 const Point& pc = point(cell->vertex((index + 3) & 3)->point()); 72 73 double edges[3]; 74 edges[0] = (CGAL::sqrt(CGAL::squared_distance(pa, pb))); 75 edges[1] = (CGAL::sqrt(CGAL::squared_distance(pa, pc))); 76 edges[2] = (CGAL::sqrt(CGAL::squared_distance(pb, pc))); 77 for (int i = 0; i < 3; ++i) 78 { 79 if (edges[i] < min_edges_length){ min_edges_length = edges[i]; } 80 if (edges[i] > max_edges_length){ max_edges_length = edges[i]; } 81 total_edges += edges[i]; 82 ++nb_edges; 83 } 84 } 85 86 double mean_edges_length = total_edges / (double)nb_edges; 87 88 typename Gt::Compute_approximate_dihedral_angle_3 approx_dihedral_angle 89 = tr.geom_traits().compute_approximate_dihedral_angle_3_object(); 90 91 std::size_t nb_tets = 0; 92 boost::unordered_set<Vertex_handle> selected_vertices; 93 std::vector<Subdomain_index> sub_ids; 94 for (Finite_cells_iterator cit = tr.finite_cells_begin(); 95 cit != tr.finite_cells_end(); 96 ++cit) 97 { 98 const Subdomain_index& si = cit->subdomain_index(); 99 if (si == Subdomain_index() || !cell_selector(cit)) 100 continue; 101 102 ++nb_tets; 103 if (std::find(sub_ids.begin(), sub_ids.end(), si) == sub_ids.end()) 104 sub_ids.push_back(cit->subdomain_index()); 105 for (int i = 0; i < 4; ++i) 106 selected_vertices.insert(cit->vertex(i)); 107 108 const Point& p0 = point(cit->vertex(0)->point()); 109 const Point& p1 = point(cit->vertex(1)->point()); 110 const Point& p2 = point(cit->vertex(2)->point()); 111 const Point& p3 = point(cit->vertex(3)->point()); 112 double v = CGAL::abs(tr.tetrahedron(cit).volume()); 113 if (v == 0.) 114 { 115 std::cout << "degenerate cell :\n\t"; 116 std::cout << p0 << "\n\t" << p1 << "\n\t" << p2 << "\n\t" << p3 << std::endl; 117 } 118 double circumradius = (v == 0.) 119 ? CGAL::sqrt(CGAL::squared_radius(p0, p1, p2)) 120 : CGAL::sqrt(CGAL::squared_radius(p0, p1, p2, p3)); 121 122 //find shortest edge 123 double edges[6]; 124 edges[0] = CGAL::sqrt(CGAL::squared_distance(p0, p1)); 125 edges[1] = CGAL::sqrt(CGAL::squared_distance(p0, p2)); 126 edges[2] = CGAL::sqrt(CGAL::squared_distance(p0, p3)); 127 edges[3] = CGAL::sqrt(CGAL::squared_distance(p2, p1)); 128 edges[4] = CGAL::sqrt(CGAL::squared_distance(p2, p3)); 129 edges[5] = CGAL::sqrt(CGAL::squared_distance(p1, p3)); 130 131 double min_edge = edges[0]; 132 for (int i = 1; i < 6; ++i) 133 { 134 if (edges[i] < min_edge) 135 min_edge = edges[i]; 136 } 137 138 double sumar = CGAL::sqrt(CGAL::squared_area(p0, p1, p2)) 139 + CGAL::sqrt(CGAL::squared_area(p1, p2, p3)) 140 + CGAL::sqrt(CGAL::squared_area(p2, p3, p0)) 141 + CGAL::sqrt(CGAL::squared_area(p3, p1, p0)); 142 double inradius = 3. * v / sumar; 143 double smallest_edge_radius_ = min_edge / circumradius*CGAL::sqrt(6.) / 4.;//*sqrt(6)/4 so that the perfect tet ratio is 1 144 double smallest_radius_radius_ = inradius / circumradius * 3.; //*3 so that the perfect tet ratio is 1 instead of 1/3 145 double biggest_v_sma_cube_ = v / std::pow(min_edge, 3) * 6. * CGAL::sqrt(2.);//*6*sqrt(2) so that the perfect tet ratio is 1 instead 146 147 if (smallest_edge_radius_ < smallest_edge_radius) 148 smallest_edge_radius = smallest_edge_radius_; 149 150 if (smallest_radius_radius_ < smallest_radius_radius) 151 smallest_radius_radius = smallest_radius_radius_; 152 153 if (biggest_v_sma_cube_ > biggest_v_sma_cube) 154 biggest_v_sma_cube = biggest_v_sma_cube_; 155 156 double a = CGAL::to_double(CGAL::abs(approx_dihedral_angle(p0, p1, p2, p3))); 157 if (a < min_dihedral_angle) { min_dihedral_angle = a; } 158 if (a > max_dihedral_angle) { max_dihedral_angle = a; } 159 total_angle += a; 160 ++nb_angle; 161 a = CGAL::to_double(CGAL::abs(approx_dihedral_angle(p0, p2, p1, p3))); 162 if (a < min_dihedral_angle) { min_dihedral_angle = a; } 163 if (a > max_dihedral_angle) { max_dihedral_angle = a; } 164 total_angle += a; 165 ++nb_angle; 166 a = CGAL::to_double(CGAL::abs(approx_dihedral_angle(p0, p3, p1, p2))); 167 if (a < min_dihedral_angle) { min_dihedral_angle = a; } 168 if (a > max_dihedral_angle) { max_dihedral_angle = a; } 169 total_angle += a; 170 ++nb_angle; 171 a = CGAL::to_double(CGAL::abs(approx_dihedral_angle(p1, p2, p0, p3))); 172 if (a < min_dihedral_angle) { min_dihedral_angle = a; } 173 if (a > max_dihedral_angle) { max_dihedral_angle = a; } 174 total_angle += a; 175 ++nb_angle; 176 a = CGAL::to_double(CGAL::abs(approx_dihedral_angle(p1, p3, p0, p2))); 177 if (a < min_dihedral_angle) { min_dihedral_angle = a; } 178 if (a > max_dihedral_angle) { max_dihedral_angle = a; } 179 total_angle += a; 180 ++nb_angle; 181 a = CGAL::to_double(CGAL::abs(approx_dihedral_angle(p2, p3, p0, p1))); 182 if (a < min_dihedral_angle) { min_dihedral_angle = a; } 183 if (a > max_dihedral_angle) { max_dihedral_angle = a; } 184 total_angle += a; 185 ++nb_angle; 186 } 187 188 std::size_t nb_subdomains = sub_ids.size(); 189 //std::size_t nb_vertices = d->c3t3.number_of_vertices_in_complex(); 190 191 std::ofstream ofs(filename); 192 if (!ofs) 193 return; 194 195 ofs << "Nb subdomains : " << nb_subdomains << std::endl; 196 ofs << "Total number of vertices : " << tr.number_of_vertices() << std::endl; 197 ofs << "Number of selected cells : " << nb_tets << std::endl; 198 ofs << "Number of selected vertices : " << selected_vertices.size() << std::endl; 199 ofs << std::endl; 200 ofs << "Min dihedral angle : " << min_dihedral_angle << std::endl; 201 ofs << "Max dihedral angle : " << max_dihedral_angle << std::endl; 202 ofs << std::endl; 203 ofs << "Shortest edge : " << min_edges_length << std::endl; 204 ofs << "Longest edge : " << max_edges_length << std::endl; 205 ofs << "Average edge length : " << mean_edges_length << std::endl; 206 207 ofs.close(); 208 } 209 210 }//end namespace internal 211 }//end namespace Tetrahedral_remeshing 212 }//end namespace CGAL 213 214 #endif // CGAL_TR_INTERNAL_COMPUTE_C3T3_STATISTICS_H 215