| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * Copyright (c) 2000-2022 Inria | ||
| 3 | * All rights reserved. | ||
| 4 | * | ||
| 5 | * Redistribution and use in source and binary forms, with or without | ||
| 6 | * modification, are permitted provided that the following conditions are met: | ||
| 7 | * | ||
| 8 | * * Redistributions of source code must retain the above copyright notice, | ||
| 9 | * this list of conditions and the following disclaimer. | ||
| 10 | * * Redistributions in binary form must reproduce the above copyright notice, | ||
| 11 | * this list of conditions and the following disclaimer in the documentation | ||
| 12 | * and/or other materials provided with the distribution. | ||
| 13 | * * Neither the name of the ALICE Project-Team nor the names of its | ||
| 14 | * contributors may be used to endorse or promote products derived from this | ||
| 15 | * software without specific prior written permission. | ||
| 16 | * | ||
| 17 | * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" | ||
| 18 | * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | ||
| 19 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | ||
| 20 | * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE | ||
| 21 | * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR | ||
| 22 | * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF | ||
| 23 | * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS | ||
| 24 | * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN | ||
| 25 | * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) | ||
| 26 | * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE | ||
| 27 | * POSSIBILITY OF SUCH DAMAGE. | ||
| 28 | * | ||
| 29 | * Contact: Bruno Levy | ||
| 30 | * | ||
| 31 | * https://www.inria.fr/fr/bruno-levy | ||
| 32 | * | ||
| 33 | * Inria, | ||
| 34 | * Domaine de Voluceau, | ||
| 35 | * 78150 Le Chesnay - Rocquencourt | ||
| 36 | * FRANCE | ||
| 37 | * | ||
| 38 | */ | ||
| 39 | |||
| 40 | #include <geogram/mesh/mesh_geometry.h> | ||
| 41 | #include <geogram/delaunay/LFS.h> | ||
| 42 | #include <geogram/voronoi/CVT.h> | ||
| 43 | #include <geogram/basic/attributes.h> | ||
| 44 | #include <geogram/basic/geometry.h> | ||
| 45 | #include <geogram/basic/logger.h> | ||
| 46 | |||
| 47 | namespace { | ||
| 48 | |||
| 49 | using namespace GEO; | ||
| 50 | |||
| 51 | /** | ||
| 52 | * \brief Computes a sizing field using local feature size | ||
| 53 | * \details The sizing field is stored into the vertex weights | ||
| 54 | * of the mesh | ||
| 55 | * \param[in] M the mesh | ||
| 56 | * \param[in] LFS the local feature size | ||
| 57 | * \param[in] gradation power to be applied to the sizing field | ||
| 58 | */ | ||
| 59 | ✗ | void compute_sizing_field_lfs( | |
| 60 | Mesh& M, const LocalFeatureSize& LFS, double gradation | ||
| 61 | ) { | ||
| 62 | // Avoid LFS points that are too close to the surface | ||
| 63 | ✗ | double min_distance2 = 0.1*surface_average_edge_length(M); | |
| 64 | ✗ | min_distance2 = min_distance2 * min_distance2; | |
| 65 | |||
| 66 | ✗ | Attribute<double> weight(M.vertices.attributes(),"weight"); | |
| 67 | ✗ | for(index_t v: M.vertices) { | |
| 68 | ✗ | double lfs2 = LFS.squared_lfs(M.vertices.point_ptr(v)); | |
| 69 | lfs2 = std::max(lfs2, min_distance2); | ||
| 70 | ✗ | double w = pow(lfs2, -2.0 * gradation); | |
| 71 | ✗ | weight[v] = w; | |
| 72 | } | ||
| 73 | ✗ | } | |
| 74 | } | ||
| 75 | |||
| 76 | /****************************************************************************/ | ||
| 77 | |||
| 78 | namespace GEO { | ||
| 79 | |||
| 80 | namespace Geom { | ||
| 81 | |||
| 82 | 130426 | vec3 mesh_facet_normal(const Mesh& M, index_t f) { | |
| 83 | vec3 result(0.0, 0.0, 0.0); | ||
| 84 |
2/2✓ Branch 0 taken 130426 times.
✓ Branch 1 taken 130426 times.
|
260852 | for(auto [p1, p2, p3]: M.facets.triangle_points(f)) { |
| 85 | result += cross(p2 - p1, p3 - p1); | ||
| 86 | } | ||
| 87 | 130426 | return result; | |
| 88 | } | ||
| 89 | |||
| 90 | ✗ | double mesh_unsigned_normal_angle( | |
| 91 | const Mesh& M, index_t f1, index_t f2 | ||
| 92 | ) { | ||
| 93 | ✗ | vec3 n1 = mesh_facet_normal(M,f1); | |
| 94 | ✗ | vec3 n2 = mesh_facet_normal(M,f2); | |
| 95 | ✗ | double l = length(n1)*length(n2); | |
| 96 | ✗ | double cos_angle = (l > 1e-30) ? dot(n1, n2)/l : 1.0; | |
| 97 | // Numerical precision problem may occur, and generate | ||
| 98 | // normalized dot products that are outside the valid | ||
| 99 | // range of acos. | ||
| 100 | geo_clamp(cos_angle, -1.0, 1.0); | ||
| 101 | ✗ | return acos(cos_angle); | |
| 102 | } | ||
| 103 | |||
| 104 | ✗ | double mesh_normal_angle(const Mesh& M, index_t c) { | |
| 105 | geo_debug_assert(M.facets.are_simplices()); | ||
| 106 | ✗ | index_t f1 = c/3; | |
| 107 | index_t f2 = M.facet_corners.adjacent_facet(c); | ||
| 108 | geo_debug_assert(f2 != NO_FACET); | ||
| 109 | ✗ | vec3 n1 = mesh_facet_normal(M,f1); | |
| 110 | ✗ | vec3 n2 = mesh_facet_normal(M,f2); | |
| 111 | ✗ | double l = length(n1)*length(n2); | |
| 112 | double sign = 1.0; | ||
| 113 | ✗ | if(dot(cross(n1,n2),mesh_corner_vector(M,c)) > 0.0) { | |
| 114 | sign = -1.0; | ||
| 115 | } | ||
| 116 | ✗ | double cos_angle = (l > 1e-30) ? dot(n1, n2)/l : 1.0; | |
| 117 | // Numerical precision problem may occur, and generate | ||
| 118 | // normalized dot products that are outside the valid | ||
| 119 | // range of acos. | ||
| 120 | geo_clamp(cos_angle, -1.0, 1.0); | ||
| 121 | ✗ | return sign*acos(cos_angle); | |
| 122 | } | ||
| 123 | |||
| 124 | 60 | double mesh_area(const Mesh& M, index_t dim) { | |
| 125 | double result = 0.0; | ||
| 126 |
2/2✓ Branch 0 taken 402841 times.
✓ Branch 1 taken 60 times.
|
402901 | for(index_t f: M.facets) { |
| 127 | 402841 | result += mesh_facet_area(M, f, dim); | |
| 128 | } | ||
| 129 | 60 | return result; | |
| 130 | } | ||
| 131 | |||
| 132 | ✗ | double GEOGRAM_API mesh_enclosed_volume(const Mesh& M) { | |
| 133 | double result = 0.0; | ||
| 134 | ✗ | for(index_t f: M.facets) { | |
| 135 | ✗ | for(auto [ p1, p2, p3 ] : M.facets.triangle_points(f)) { | |
| 136 | ✗ | result += dot(p1,cross(p2,p3)) / 6.0; | |
| 137 | } | ||
| 138 | } | ||
| 139 | ✗ | return ::fabs(result); | |
| 140 | } | ||
| 141 | |||
| 142 | } | ||
| 143 | |||
| 144 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | void compute_normals(Mesh& M) { |
| 145 |
1/2✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
|
5 | if(M.vertices.dimension() < 6) { |
| 146 | 5 | M.vertices.set_dimension(6); | |
| 147 | } else { | ||
| 148 | ✗ | for(index_t i: M.vertices) { | |
| 149 | ✗ | Geom::mesh_vertex_normal_ref(M, i) = vec3(0.0, 0.0, 0.0); | |
| 150 | } | ||
| 151 | } | ||
| 152 |
2/2✓ Branch 0 taken 7360 times.
✓ Branch 1 taken 5 times.
|
7365 | for(index_t f: M.facets) { |
| 153 | 7360 | vec3 N = Geom::mesh_facet_normal(M, f); | |
| 154 |
2/2✓ Branch 0 taken 22080 times.
✓ Branch 1 taken 7360 times.
|
29440 | for(index_t corner: M.facets.corners(f)) { |
| 155 | index_t v = M.facet_corners.vertex(corner); | ||
| 156 | Geom::mesh_vertex_normal_ref(M, v) += N; | ||
| 157 | } | ||
| 158 | } | ||
| 159 |
2/2✓ Branch 0 taken 3772 times.
✓ Branch 1 taken 5 times.
|
3777 | for(index_t i: M.vertices) { |
| 160 | 3772 | Geom::mesh_vertex_normal_ref(M, i) = normalize( | |
| 161 | Geom::mesh_vertex_normal(M, i) | ||
| 162 | ); | ||
| 163 | } | ||
| 164 | 5 | } | |
| 165 | |||
| 166 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | void simple_Laplacian_smooth(Mesh& M, index_t nb_iter, bool normals_only) { |
| 167 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
5 | geo_assert(M.vertices.dimension() >= 6); |
| 168 |
1/2✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
|
5 | std::vector<vec3> p(M.vertices.nb()); |
| 169 |
1/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
5 | std::vector<double> c(M.vertices.nb()); |
| 170 | |||
| 171 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 5 times.
|
20 | for(index_t k = 0; k < nb_iter; k++) { |
| 172 |
1/2✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
|
15 | p.assign(M.vertices.nb(), vec3(0.0, 0.0, 0.0)); |
| 173 |
1/4✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
15 | c.assign(M.vertices.nb(), 0); |
| 174 |
2/2✓ Branch 0 taken 22080 times.
✓ Branch 1 taken 15 times.
|
22095 | for(index_t f: M.facets) { |
| 175 | index_t b = M.facets.corners_begin(f); | ||
| 176 | index_t e = M.facets.corners_end(f); | ||
| 177 |
2/2✓ Branch 0 taken 66240 times.
✓ Branch 1 taken 22080 times.
|
88320 | for(index_t c1 = b; c1 != e; c1++) { |
| 178 |
2/2✓ Branch 0 taken 44160 times.
✓ Branch 1 taken 22080 times.
|
66240 | index_t c2 = (c1 == e - 1) ? b : c1 + 1; |
| 179 | index_t v1 = M.facet_corners.vertex(c1); | ||
| 180 | index_t v2 = M.facet_corners.vertex(c2); | ||
| 181 |
2/2✓ Branch 0 taken 33120 times.
✓ Branch 1 taken 33120 times.
|
66240 | if(v1 < v2) { |
| 182 | double a = 1.0; | ||
| 183 |
1/2✓ Branch 0 taken 33120 times.
✗ Branch 1 not taken.
|
33120 | c[v1] += a; |
| 184 | 33120 | c[v2] += a; | |
| 185 |
1/2✓ Branch 0 taken 33120 times.
✗ Branch 1 not taken.
|
33120 | if(normals_only) { |
| 186 | p[v1] += a * Geom::mesh_vertex_normal(M, v2); | ||
| 187 | p[v2] += a * Geom::mesh_vertex_normal(M, v1); | ||
| 188 | } else { | ||
| 189 | p[v1] += a * M.vertices.point(v2); | ||
| 190 | p[v2] += a * M.vertices.point(v1); | ||
| 191 | } | ||
| 192 | } | ||
| 193 | } | ||
| 194 | } | ||
| 195 |
2/2✓ Branch 0 taken 11316 times.
✓ Branch 1 taken 15 times.
|
11331 | for(index_t v: M.vertices) { |
| 196 |
1/2✓ Branch 0 taken 11316 times.
✗ Branch 1 not taken.
|
11316 | if(normals_only) { |
| 197 |
1/2✓ Branch 0 taken 11316 times.
✗ Branch 1 not taken.
|
11316 | double l = length(p[v]); |
| 198 |
1/2✓ Branch 0 taken 11316 times.
✗ Branch 1 not taken.
|
11316 | if(l > 1e-30) { |
| 199 | 11316 | Geom::mesh_vertex_normal_ref(M,v) = (1.0 / l) * p[v]; | |
| 200 | } | ||
| 201 | } else { | ||
| 202 | ✗ | M.vertices.point(v) = 1.0 / c[v] * p[v]; | |
| 203 | } | ||
| 204 | } | ||
| 205 | } | ||
| 206 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(!normals_only) { |
| 207 | ✗ | compute_normals(M); | |
| 208 | } | ||
| 209 | 5 | } | |
| 210 | |||
| 211 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | void get_bbox(const Mesh& M, double* xyzmin, double* xyzmax) { |
| 212 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
78 | geo_assert(M.vertices.dimension() >= 3); |
| 213 |
2/2✓ Branch 0 taken 234 times.
✓ Branch 1 taken 78 times.
|
312 | for(index_t c = 0; c < 3; c++) { |
| 214 | 234 | xyzmin[c] = Numeric::max_float64(); | |
| 215 | 234 | xyzmax[c] = Numeric::min_float64(); | |
| 216 | } | ||
| 217 |
2/2✓ Branch 0 taken 233886 times.
✓ Branch 1 taken 78 times.
|
233964 | for(const vec3& p: M.vertices.points()) { |
| 218 |
2/2✓ Branch 0 taken 701658 times.
✓ Branch 1 taken 233886 times.
|
935544 | for(index_t c = 0; c < 3; c++) { |
| 219 |
2/2✓ Branch 0 taken 7835 times.
✓ Branch 1 taken 693823 times.
|
701658 | xyzmin[c] = std::min(xyzmin[c], p[c]); |
| 220 |
2/2✓ Branch 0 taken 3878 times.
✓ Branch 1 taken 697780 times.
|
705536 | xyzmax[c] = std::max(xyzmax[c], p[c]); |
| 221 | } | ||
| 222 | } | ||
| 223 | 78 | } | |
| 224 | |||
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
|
78 | double bbox_diagonal(const Mesh& M) { |
| 226 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
78 | geo_assert(M.vertices.dimension() >= 3); |
| 227 | double xyzmin[3]; | ||
| 228 | double xyzmax[3]; | ||
| 229 | 78 | get_bbox(M, xyzmin, xyzmax); | |
| 230 | 78 | return ::sqrt( | |
| 231 | 78 | geo_sqr(xyzmax[0] - xyzmin[0]) + | |
| 232 | 78 | geo_sqr(xyzmax[1] - xyzmin[1]) + | |
| 233 | 78 | geo_sqr(xyzmax[2] - xyzmin[2]) | |
| 234 | 78 | ); | |
| 235 | } | ||
| 236 | |||
| 237 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | void set_anisotropy(Mesh& M, double s) { |
| 238 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(M.vertices.dimension() < 6) { |
| 239 | ✗ | compute_normals(M); | |
| 240 | } | ||
| 241 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
|
5 | if(s == 0.0) { |
| 242 | ✗ | unset_anisotropy(M); | |
| 243 | ✗ | return; | |
| 244 | } | ||
| 245 | 5 | s *= bbox_diagonal(M); | |
| 246 |
2/2✓ Branch 0 taken 3772 times.
✓ Branch 1 taken 5 times.
|
3777 | for(index_t i: M.vertices) { |
| 247 | 3772 | Geom::mesh_vertex_normal_ref(M, i) = | |
| 248 | 3772 | s * normalize(Geom::mesh_vertex_normal(M, i)); | |
| 249 | } | ||
| 250 | } | ||
| 251 | |||
| 252 | ✗ | void unset_anisotropy(Mesh& M) { | |
| 253 | ✗ | if(M.vertices.dimension() < 6) { | |
| 254 | return; | ||
| 255 | } | ||
| 256 | ✗ | for(index_t i: M.vertices) { | |
| 257 | ✗ | Geom::mesh_vertex_normal_ref(M, i) = normalize( | |
| 258 | Geom::mesh_vertex_normal(M, i) | ||
| 259 | ); | ||
| 260 | } | ||
| 261 | } | ||
| 262 | |||
| 263 | ✗ | void compute_sizing_field( | |
| 264 | Mesh& M, double gradation, index_t nb_lfs_samples | ||
| 265 | ) { | ||
| 266 | ✗ | if(nb_lfs_samples != 0) { | |
| 267 | ✗ | Logger::out("LFS") << "Sampling surface" << std::endl; | |
| 268 | ✗ | CentroidalVoronoiTesselation CVT(&M, 3); | |
| 269 | ✗ | CVT.compute_initial_sampling(nb_lfs_samples); | |
| 270 | ✗ | Logger::out("LFS") << "Optimizing sampling (Lloyd)" << std::endl; | |
| 271 | ✗ | CVT.Lloyd_iterations(5); | |
| 272 | ✗ | Logger::out("LFS") << "Optimizing sampling (Newton)" << std::endl; | |
| 273 | ✗ | CVT.Newton_iterations(10); | |
| 274 | ✗ | Logger::out("LFS") << "Computing medial axis" << std::endl; | |
| 275 | ✗ | LocalFeatureSize LFS(CVT.nb_points(), CVT.embedding(0)); | |
| 276 | ✗ | Logger::out("LFS") << "Computing sizing field" << std::endl; | |
| 277 | ✗ | compute_sizing_field_lfs(M, LFS, gradation); | |
| 278 | ✗ | } else { | |
| 279 | ✗ | if(M.vertices.dimension() == 3) { | |
| 280 | ✗ | LocalFeatureSize LFS(M.vertices.nb(), M.vertices.point_ptr(0)); | |
| 281 | ✗ | compute_sizing_field_lfs(M, LFS, gradation); | |
| 282 | ✗ | } else { | |
| 283 | std::vector<double> pts; | ||
| 284 | ✗ | pts.reserve(M.vertices.nb() * 3); | |
| 285 | ✗ | for(index_t v: M.vertices) { | |
| 286 | pts.push_back(M.vertices.point_ptr(v)[0]); | ||
| 287 | ✗ | pts.push_back(M.vertices.point_ptr(v)[1]); | |
| 288 | ✗ | pts.push_back(M.vertices.point_ptr(v)[2]); | |
| 289 | } | ||
| 290 | ✗ | LocalFeatureSize LFS(M.vertices.nb(), pts.data()); | |
| 291 | ✗ | compute_sizing_field_lfs(M, LFS, gradation); | |
| 292 | ✗ | } | |
| 293 | } | ||
| 294 | ✗ | } | |
| 295 | |||
| 296 | ✗ | void normalize_embedding_area(Mesh& M) { | |
| 297 | ✗ | if(M.vertices.dimension() == 3) { | |
| 298 | ✗ | if(M.vertices.attributes().is_defined("weight")) { | |
| 299 | ✗ | M.vertices.attributes().delete_attribute_store("weight"); | |
| 300 | } | ||
| 301 | ✗ | return; | |
| 302 | } | ||
| 303 | ✗ | Attribute<double> weight(M.vertices.attributes(), "weight"); | |
| 304 | ✗ | std::vector<double> area3d(M.vertices.nb(), 0.0); | |
| 305 | ✗ | std::vector<double> areaNd(M.vertices.nb(), 0.0); | |
| 306 | ✗ | for(index_t f: M.facets) { | |
| 307 | ✗ | double A3d = Geom::mesh_facet_area(M, f, 3); | |
| 308 | ✗ | double ANd = Geom::mesh_facet_area(M, f); | |
| 309 | ✗ | for(index_t c: M.facets.corners(f)) { | |
| 310 | index_t v = M.facet_corners.vertex(c); | ||
| 311 | ✗ | area3d[v] += A3d; | |
| 312 | ✗ | areaNd[v] += ANd; | |
| 313 | } | ||
| 314 | } | ||
| 315 | ✗ | for(index_t v: M.vertices) { | |
| 316 | ✗ | double A3d = area3d[v]; | |
| 317 | ✗ | double ANd = areaNd[v]; | |
| 318 | ANd = std::max(ANd, 1e-6); | ||
| 319 | ✗ | double w = ::pow(A3d / ANd, 2.0); | |
| 320 | ✗ | weight[v] = w; | |
| 321 | } | ||
| 322 | } | ||
| 323 | |||
| 324 | ✗ | double mesh_cell_volume( | |
| 325 | const Mesh& M, index_t c | ||
| 326 | ) { | ||
| 327 | ✗ | geo_assert(M.vertices.dimension() >= 3); | |
| 328 | |||
| 329 | // Connectors are virtual cells, they | ||
| 330 | // do not have a geometry. | ||
| 331 | ✗ | if(M.cells.type(c) == MESH_CONNECTOR) { | |
| 332 | return 0.0; | ||
| 333 | } | ||
| 334 | |||
| 335 | // Easy case: tetrahedra. | ||
| 336 | ✗ | if(M.cells.type(c) == MESH_TET) { | |
| 337 | const vec3& p0 = M.cells.point(c,0); | ||
| 338 | const vec3& p1 = M.cells.point(c,1); | ||
| 339 | const vec3& p2 = M.cells.point(c,2); | ||
| 340 | const vec3& p3 = M.cells.point(c,3); | ||
| 341 | ✗ | return ::fabs(Geom::tetra_signed_volume(p0,p1,p2,p3)); | |
| 342 | } | ||
| 343 | |||
| 344 | // Arbitrary cells are decomposed into tetrahedra, with one | ||
| 345 | // vertex in the center of the cell, and one vertex in the | ||
| 346 | // center of each face. This ensures that two adjacent cells | ||
| 347 | // are not overlapping (if simply triangulating the faces, | ||
| 348 | // there would be tiny overlaps / tiny gaps that would introduce | ||
| 349 | // errors in the total volume of the mesh). | ||
| 350 | // | ||
| 351 | // Note that the center point may fall outside the cell in some | ||
| 352 | // degenerate configurations. It is not a problem since we compute | ||
| 353 | // signed tetrahedra volumes, in such a way that overlapping volumes | ||
| 354 | // will cancel-out in such a configuration. | ||
| 355 | // Therefore, we could take an arbitrary point as the enter point, | ||
| 356 | // including the origin, but taking the center probably makes | ||
| 357 | // computations more stable, by cancelling the translations. | ||
| 358 | |||
| 359 | double result = 0.0; | ||
| 360 | ✗ | vec3 center{0.0, 0.0, 0.0}; | |
| 361 | index_t nbcv = M.cells.nb_vertices(c); | ||
| 362 | ✗ | for(const vec3& p: M.cells.points(c)) { | |
| 363 | center += p; | ||
| 364 | } | ||
| 365 | ✗ | center /= double(nbcv); | |
| 366 | ✗ | for(index_t lf=0; lf<M.cells.nb_facets(c); ++lf) { | |
| 367 | index_t nbcfv = M.cells.facet_nb_vertices(c,lf); | ||
| 368 | ✗ | if(nbcfv == 3) { | |
| 369 | ✗ | const vec3& p1 = M.vertices.point(M.cells.facet_vertex(c,lf,0)); | |
| 370 | ✗ | const vec3& p2 = M.vertices.point(M.cells.facet_vertex(c,lf,1)); | |
| 371 | ✗ | const vec3& p3 = M.vertices.point(M.cells.facet_vertex(c,lf,2)); | |
| 372 | ✗ | result += Geom::tetra_signed_volume(center, p1, p2, p3); | |
| 373 | } else { | ||
| 374 | ✗ | vec3 facet_center {0.0, 0.0, 0.0}; | |
| 375 | ✗ | for(index_t lfv=0; lfv<nbcfv; ++lfv) { | |
| 376 | facet_center += | ||
| 377 | ✗ | M.vertices.point(M.cells.facet_vertex(c,lf,lfv)); | |
| 378 | } | ||
| 379 | ✗ | facet_center /= double(nbcfv); | |
| 380 | ✗ | for(index_t lfv1=0; lfv1<nbcfv; ++lfv1) { | |
| 381 | ✗ | index_t lfv2 = (lfv1 + 1) % nbcfv; | |
| 382 | const vec3& p1 = | ||
| 383 | ✗ | M.vertices.point(M.cells.facet_vertex(c,lf,lfv1)); | |
| 384 | const vec3& p2 = | ||
| 385 | ✗ | M.vertices.point(M.cells.facet_vertex(c,lf,lfv2)); | |
| 386 | ✗ | result += | |
| 387 | ✗ | Geom::tetra_signed_volume(center, facet_center, p1, p2); | |
| 388 | } | ||
| 389 | } | ||
| 390 | } | ||
| 391 | ✗ | return ::fabs(result); | |
| 392 | } | ||
| 393 | |||
| 394 | |||
| 395 | ✗ | double mesh_cells_volume(const Mesh& M) { | |
| 396 | double result = 0.0; | ||
| 397 | ✗ | for(index_t c: M.cells) { | |
| 398 | ✗ | result += mesh_cell_volume(M,c); | |
| 399 | } | ||
| 400 | ✗ | return result; | |
| 401 | } | ||
| 402 | |||
| 403 | ✗ | vec3 GEOGRAM_API mesh_cell_facet_normal( | |
| 404 | const Mesh& M, index_t c, index_t lf | ||
| 405 | ) { | ||
| 406 | geo_debug_assert(M.vertices.dimension() >= 3); | ||
| 407 | geo_debug_assert(c < M.cells.nb()); | ||
| 408 | geo_debug_assert(lf < M.cells.nb_facets(c)); | ||
| 409 | |||
| 410 | ✗ | index_t v1 = M.cells.facet_vertex(c,lf,0); | |
| 411 | ✗ | index_t v2 = M.cells.facet_vertex(c,lf,1); | |
| 412 | ✗ | index_t v3 = M.cells.facet_vertex(c,lf,2); | |
| 413 | |||
| 414 | const vec3& p1 = M.vertices.point(v1); | ||
| 415 | const vec3& p2 = M.vertices.point(v2); | ||
| 416 | const vec3& p3 = M.vertices.point(v3); | ||
| 417 | |||
| 418 | ✗ | return cross(p2 - p1, p3 - p1); | |
| 419 | } | ||
| 420 | |||
| 421 | |||
| 422 | |||
| 423 | 1 | double surface_average_edge_length(const Mesh& M) { | |
| 424 | double result = 0.0; | ||
| 425 | index_t count = 0; | ||
| 426 |
2/2✓ Branch 0 taken 9804 times.
✓ Branch 1 taken 1 times.
|
9805 | for(index_t f: M.facets) { |
| 427 |
2/2✓ Branch 0 taken 29412 times.
✓ Branch 1 taken 9804 times.
|
39216 | for(index_t c1: M.facets.corners(f)) { |
| 428 | index_t c2 = M.facets.next_corner_around_facet(f,c1); | ||
| 429 | index_t v1 = M.facet_corners.vertex(c1); | ||
| 430 | index_t v2 = M.facet_corners.vertex(c2); | ||
| 431 | 29412 | result += Geom::distance( | |
| 432 | M.vertices.point_ptr(v1), | ||
| 433 | M.vertices.point_ptr(v2), | ||
| 434 | coord_index_t(M.vertices.dimension()) | ||
| 435 | ); | ||
| 436 | 29412 | ++count; | |
| 437 | } | ||
| 438 | } | ||
| 439 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | if(count != 0) { |
| 440 | 1 | result /= double(count); | |
| 441 | } | ||
| 442 | 1 | return result; | |
| 443 | } | ||
| 444 | } | ||
| 445 |