| 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 | #ifndef GEOGRAM_MESH_MESH_SAMPLING | ||
| 41 | #define GEOGRAM_MESH_MESH_SAMPLING | ||
| 42 | |||
| 43 | #include <geogram/basic/common.h> | ||
| 44 | #include <geogram/mesh/mesh.h> | ||
| 45 | #include <geogram/mesh/mesh_geometry.h> | ||
| 46 | #include <geogram/basic/geometry_nd.h> | ||
| 47 | #include <geogram/basic/logger.h> | ||
| 48 | #include <algorithm> | ||
| 49 | |||
| 50 | /** | ||
| 51 | * \file geogram/mesh/mesh_sampling.h | ||
| 52 | * \brief Functions to generate random samples | ||
| 53 | * on surfacic and in volumetric meshes | ||
| 54 | */ | ||
| 55 | |||
| 56 | namespace GEO { | ||
| 57 | |||
| 58 | /** | ||
| 59 | * \brief Computes the mass of a mesh facet. | ||
| 60 | * \details The function can optionally take into account the vertex | ||
| 61 | * weights. | ||
| 62 | * \param[in] mesh the surface mesh | ||
| 63 | * \param[in] f a facet index in \p mesh | ||
| 64 | * \param[in] vertex_weight a reference to a vertex attribute. If | ||
| 65 | * it is bound, it is used to weight the vertices. | ||
| 66 | * \return the mass of facet \p f in \p mesh | ||
| 67 | */ | ||
| 68 | template <index_t DIM> | ||
| 69 |
1/2✓ Branch 0 taken 156112 times.
✗ Branch 1 not taken.
|
171708 | inline double mesh_facet_mass( |
| 70 | const Mesh& mesh, | ||
| 71 | index_t f, | ||
| 72 | Attribute<double>& vertex_weight | ||
| 73 | ) { | ||
| 74 | geo_debug_assert(mesh.facets.are_simplices()); | ||
| 75 | geo_debug_assert(mesh.vertices.dimension() >= DIM); | ||
| 76 | |||
| 77 | index_t v1 = mesh.facets.vertex(f,0); | ||
| 78 | index_t v2 = mesh.facets.vertex(f,1); | ||
| 79 | index_t v3 = mesh.facets.vertex(f,2); | ||
| 80 | |||
| 81 | if(vertex_weight.is_bound()) { | ||
| 82 | ✗ | return Geom::triangle_mass( | |
| 83 | mesh.vertices.point<DIM>(v1), | ||
| 84 | mesh.vertices.point<DIM>(v2), | ||
| 85 | mesh.vertices.point<DIM>(v3), | ||
| 86 | vertex_weight[v1], | ||
| 87 | vertex_weight[v2], | ||
| 88 | vertex_weight[v3] | ||
| 89 | ✗ | ); | |
| 90 | } | ||
| 91 | 29360 | return Geom::triangle_area( | |
| 92 | mesh.vertices.point<DIM>(v1), | ||
| 93 | mesh.vertices.point<DIM>(v2), | ||
| 94 | mesh.vertices.point<DIM>(v3) | ||
| 95 | 171708 | ); | |
| 96 | } | ||
| 97 | |||
| 98 | /** | ||
| 99 | * \brief Generates a set of random samples over a surfacic mesh. | ||
| 100 | * \param[in] mesh the mesh | ||
| 101 | * \param[out] p pointer to an array of generated samples, of size | ||
| 102 | * \p nb_points times DIM. To be allocated by the caller. | ||
| 103 | * \param[in] nb_points number of points to generate | ||
| 104 | * \param[in] weight a reference to a vertex attribute. If bound, it | ||
| 105 | * is taken into account. | ||
| 106 | * \param[in] facets_begin_in if specified, first index of the facet | ||
| 107 | * sequence in which points should be generated. If left unspecified (-1), | ||
| 108 | * points are generated over all the facets of the mesh. | ||
| 109 | * \param[in] facets_end_in if specified, one position past the last | ||
| 110 | * index of the facet sequence in which points should be generated. | ||
| 111 | * If left unspecified (-1), points are generated over all the facets | ||
| 112 | * of the mesh. | ||
| 113 | * \tparam DIM dimension of the points, specified as a template argument | ||
| 114 | * for efficiency reasons | ||
| 115 | * \return true if everything went OK, false otherwise. Whenever all the | ||
| 116 | * points land in the same facet, the function returns false to notify | ||
| 117 | * a potential numerical problem. | ||
| 118 | */ | ||
| 119 | template <index_t DIM> | ||
| 120 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
|
28 | inline bool mesh_generate_random_samples_on_surface( |
| 121 | const Mesh& mesh, | ||
| 122 | double* p, | ||
| 123 | index_t nb_points, | ||
| 124 | Attribute<double>& weight, | ||
| 125 | index_t facets_begin_in = NO_INDEX, | ||
| 126 | index_t facets_end_in = NO_INDEX | ||
| 127 | ) { | ||
| 128 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
28 | geo_assert(mesh.facets.are_simplices()); |
| 129 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
28 | geo_assert(mesh.vertices.dimension() >= DIM); |
| 130 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
28 | geo_assert(mesh.facets.nb() > 0); |
| 131 | |||
| 132 | index_t facets_begin = 0; | ||
| 133 | index_t facets_end = mesh.facets.nb(); | ||
| 134 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
|
28 | if(facets_begin_in != NO_INDEX) { |
| 135 | facets_begin = facets_begin_in; | ||
| 136 | } | ||
| 137 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
|
28 | if(facets_end_in != NO_INDEX) { |
| 138 | facets_end = facets_end_in; | ||
| 139 | } | ||
| 140 | |||
| 141 | typedef vecng<DIM, double> Point; | ||
| 142 | |||
| 143 | // To ensure reproducibility accros successive | ||
| 144 | // runs, reset the random number generator. | ||
| 145 | 28 | Numeric::random_reset(); | |
| 146 | |||
| 147 | vector<double> s(nb_points); | ||
| 148 |
2/2✓ Branch 0 taken 795359 times.
✓ Branch 1 taken 22 times.
|
825387 | for(index_t i = 0; i < nb_points; i++) { |
| 149 |
1/2✓ Branch 1 taken 795359 times.
✗ Branch 2 not taken.
|
825359 | s[i] = Numeric::random_float64(); |
| 150 | } | ||
| 151 | 28 | std::sort(s.begin(), s.end()); | |
| 152 | |||
| 153 | double Atot = 0.0; | ||
| 154 |
2/2✓ Branch 0 taken 78056 times.
✓ Branch 1 taken 22 times.
|
85882 | for(index_t t = facets_begin; t < facets_end; ++t) { |
| 155 | 85854 | double At = mesh_facet_mass<DIM>(mesh, t, weight); | |
| 156 | 85854 | Atot += At; | |
| 157 | } | ||
| 158 | |||
| 159 | index_t first_t = NO_INDEX; | ||
| 160 | 28 | index_t last_t = 0; | |
| 161 | |||
| 162 | 28 | index_t cur_t = facets_begin; | |
| 163 | 28 | double cur_s = | |
| 164 | 28 | mesh_facet_mass<DIM>(mesh, facets_begin, weight) / Atot; | |
| 165 |
2/2✓ Branch 0 taken 795359 times.
✓ Branch 1 taken 22 times.
|
825387 | for(index_t i = 0; i < nb_points; i++) { |
| 166 | geo_debug_assert(i < s.size()); | ||
| 167 |
3/4✓ Branch 0 taken 78034 times.
✓ Branch 1 taken 795359 times.
✓ Branch 2 taken 78034 times.
✗ Branch 3 not taken.
|
911185 | while(s[i] > cur_s && cur_t < facets_end - 1) { |
| 168 | 85826 | cur_t++; | |
| 169 | geo_debug_assert(cur_t < facets_end); | ||
| 170 | 85826 | cur_s += mesh_facet_mass<DIM>(mesh, cur_t, weight) / Atot; | |
| 171 | } | ||
| 172 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 795337 times.
|
825359 | if(first_t == NO_INDEX) { |
| 173 | 28 | first_t = cur_t; | |
| 174 | } | ||
| 175 |
1/2✓ Branch 0 taken 795359 times.
✗ Branch 1 not taken.
|
825359 | last_t = std::max(last_t, cur_t); |
| 176 | |||
| 177 | // TODO: take weights into account | ||
| 178 | // with a new random_point_in_triangle_weighted() | ||
| 179 | // function. | ||
| 180 | index_t v1 = mesh.facets.vertex(cur_t,0); | ||
| 181 | index_t v2 = mesh.facets.vertex(cur_t,1); | ||
| 182 | index_t v3 = mesh.facets.vertex(cur_t,2); | ||
| 183 |
1/2✓ Branch 1 taken 795359 times.
✗ Branch 2 not taken.
|
825359 | Point cur_p = Geom::random_point_in_triangle( |
| 184 | mesh.vertices.point<DIM>(v1), | ||
| 185 | mesh.vertices.point<DIM>(v2), | ||
| 186 | mesh.vertices.point<DIM>(v3) | ||
| 187 | ); | ||
| 188 |
2/2✓ Branch 0 taken 2446077 times.
✓ Branch 1 taken 795359 times.
|
3421436 | for(coord_index_t coord = 0; coord < DIM; coord++) { |
| 189 | 2596077 | p[i * DIM + coord] = cur_p[coord]; | |
| 190 | } | ||
| 191 | } | ||
| 192 |
2/4✓ Branch 0 taken 22 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 22 times.
|
28 | if(mesh.facets.nb() > 1 && last_t == first_t) { |
| 193 | ✗ | Logger::warn("Sampler") | |
| 194 | << "Did put all the points in the same triangle" | ||
| 195 | << std::endl; | ||
| 196 | ✗ | return false; | |
| 197 | } | ||
| 198 | return true; | ||
| 199 | } | ||
| 200 | |||
| 201 | /************************************************************************/ | ||
| 202 | |||
| 203 | /** | ||
| 204 | * \brief Computes the mass of a mesh tetrahedron. | ||
| 205 | * \details The function can optionally take into account the vertex | ||
| 206 | * weights. | ||
| 207 | * \param[in] mesh the surface mesh | ||
| 208 | * \param[in] t a tetrahedron index in \p mesh | ||
| 209 | * \return the mass of tetrahedron \p t in \p mesh | ||
| 210 | */ | ||
| 211 | template <index_t DIM> | ||
| 212 |
1/2✓ Branch 0 taken 3072 times.
✗ Branch 1 not taken.
|
6144 | inline double mesh_tetra_mass( |
| 213 | const Mesh& mesh, | ||
| 214 | index_t t | ||
| 215 | ) { | ||
| 216 | geo_debug_assert(mesh.vertices.dimension() >= DIM); | ||
| 217 | |||
| 218 | 4608 | double result = Geom::tetra_volume( | |
| 219 | mesh.cells.point<DIM>(t,0), | ||
| 220 | mesh.cells.point<DIM>(t,1), | ||
| 221 | mesh.cells.point<DIM>(t,2), | ||
| 222 | mesh.cells.point<DIM>(t,3) | ||
| 223 | ); | ||
| 224 | |||
| 225 | 6144 | return result; | |
| 226 | } | ||
| 227 | |||
| 228 | /** | ||
| 229 | * \brief Computes the mass of a mesh tetrahedron. | ||
| 230 | * \details The function can optionally take into account the vertex | ||
| 231 | * weights. | ||
| 232 | * \param[in] mesh the surface mesh | ||
| 233 | * \param[in] t a tetrahedron index in \p mesh | ||
| 234 | * \param[in] weight a reference to a vertex weight attribute. If it | ||
| 235 | * is bound, it is taken into account in mass computation | ||
| 236 | * \return the mass of tetrahedron \p t in \p mesh | ||
| 237 | */ | ||
| 238 | template <index_t DIM> | ||
| 239 | 6144 | inline double mesh_tetra_mass( | |
| 240 | const Mesh& mesh, | ||
| 241 | index_t t, | ||
| 242 | const Attribute<double>& weight | ||
| 243 | ) { | ||
| 244 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 768 times.
|
6144 | double result = mesh_tetra_mass<DIM>(mesh, t); |
| 245 | |||
| 246 | if(weight.is_bound()) { | ||
| 247 | index_t v0 = mesh.cells.vertex(t, 0); | ||
| 248 | index_t v1 = mesh.cells.vertex(t, 1); | ||
| 249 | index_t v2 = mesh.cells.vertex(t, 2); | ||
| 250 | index_t v3 = mesh.cells.vertex(t, 3); | ||
| 251 | ✗ | result *= ( | |
| 252 | ✗ | weight[v0] + weight[v1] + weight[v2] + weight[v3] | |
| 253 | ✗ | ) / 4.0; | |
| 254 | } | ||
| 255 | |||
| 256 | 6144 | return result; | |
| 257 | } | ||
| 258 | |||
| 259 | /** | ||
| 260 | * \brief Generates a set of random samples in a volumetric mesh. | ||
| 261 | * \param[in] mesh the mesh | ||
| 262 | * \param[out] p pointer to an array of generated samples, of size | ||
| 263 | * \p nb_points times DIM. To be allocated by the caller. | ||
| 264 | * \param[in] nb_points number of points to generate | ||
| 265 | * \param[in] vertex_weight if bound, vertex weights are taken into account | ||
| 266 | * \param[in] tets_begin_in if specified, first index of the tetrahedron | ||
| 267 | * sequence in which points should be generated. If left unspecified (-1), | ||
| 268 | * points are generated over all the tetrahedra of the mesh. | ||
| 269 | * \param[in] tets_end_in if specified, one position past the last | ||
| 270 | * index of the tetrahedron sequence in which points should be generated. | ||
| 271 | * If left unspecified (-1), points are generated over all the tetrahedra | ||
| 272 | * of the mesh. | ||
| 273 | * \tparam DIM dimension of the points, specified as a template argument | ||
| 274 | * for efficiency reasons | ||
| 275 | * \return true if everything went OK, false otherwise. Whenever all the | ||
| 276 | * points land in the same tetrahedron, the function returns false | ||
| 277 | * to notify potential numerical problem. | ||
| 278 | */ | ||
| 279 | template <index_t DIM> | ||
| 280 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
8 | inline bool mesh_generate_random_samples_in_volume( |
| 281 | const Mesh& mesh, | ||
| 282 | double* p, | ||
| 283 | index_t nb_points, | ||
| 284 | Attribute<double>& vertex_weight, | ||
| 285 | index_t tets_begin_in = NO_INDEX, | ||
| 286 | index_t tets_end_in = NO_INDEX | ||
| 287 | ) { | ||
| 288 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8 | geo_assert(mesh.vertices.dimension() >= DIM); |
| 289 |
1/6✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
8 | geo_assert(mesh.cells.nb() > 0); |
| 290 | |||
| 291 | index_t tets_begin = 0; | ||
| 292 | index_t tets_end = mesh.cells.nb(); | ||
| 293 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
8 | if(tets_begin_in != NO_INDEX) { |
| 294 | tets_begin = tets_begin_in; | ||
| 295 | } | ||
| 296 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
8 | if(tets_end_in != NO_INDEX) { |
| 297 | tets_end = tets_end_in; | ||
| 298 | } | ||
| 299 | |||
| 300 | typedef vecng<DIM, double> Point; | ||
| 301 | |||
| 302 | // To ensure reproducibility accros successive | ||
| 303 | // runs, reset the random number generator. | ||
| 304 | 8 | Numeric::random_reset(); | |
| 305 | |||
| 306 | vector<double> s(nb_points); | ||
| 307 |
2/2✓ Branch 0 taken 1200 times.
✓ Branch 1 taken 4 times.
|
2408 | for(index_t i = 0; i < nb_points; i++) { |
| 308 |
1/2✓ Branch 1 taken 1200 times.
✗ Branch 2 not taken.
|
2400 | s[i] = Numeric::random_float64(); |
| 309 | } | ||
| 310 | 8 | std::sort(s.begin(), s.end()); | |
| 311 | |||
| 312 | double Vtot = 0.0; | ||
| 313 |
2/2✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 4 times.
|
3080 | for(index_t t = tets_begin; t < tets_end; ++t) { |
| 314 | 3072 | double Vt = mesh_tetra_mass<DIM>(mesh, t, vertex_weight); | |
| 315 | 3072 | Vtot += Vt; | |
| 316 | } | ||
| 317 | |||
| 318 | index_t first_t = NO_INDEX; | ||
| 319 | 8 | index_t last_t = 0; | |
| 320 | |||
| 321 | 8 | index_t cur_t = tets_begin; | |
| 322 | 8 | double cur_s = | |
| 323 | 8 | mesh_tetra_mass<DIM>(mesh, tets_begin, vertex_weight) / Vtot; | |
| 324 |
2/2✓ Branch 0 taken 1200 times.
✓ Branch 1 taken 4 times.
|
2408 | for(index_t i = 0; i < nb_points; i++) { |
| 325 | geo_debug_assert(i < s.size()); | ||
| 326 |
3/4✓ Branch 0 taken 1532 times.
✓ Branch 1 taken 1200 times.
✓ Branch 2 taken 1532 times.
✗ Branch 3 not taken.
|
5464 | while(s[i] > cur_s && cur_t < tets_end - 1) { |
| 327 | 3064 | cur_t++; | |
| 328 | geo_debug_assert(cur_t < tets_end); | ||
| 329 | 3064 | cur_s += mesh_tetra_mass<DIM>( | |
| 330 | mesh, cur_t, vertex_weight | ||
| 331 | 3064 | ) / Vtot; | |
| 332 | } | ||
| 333 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1196 times.
|
2400 | if(first_t == NO_INDEX) { |
| 334 | 8 | first_t = cur_t; | |
| 335 | } | ||
| 336 |
1/2✓ Branch 0 taken 1200 times.
✗ Branch 1 not taken.
|
2400 | last_t = std::max(last_t, cur_t); |
| 337 | |||
| 338 | index_t v0 = mesh.cells.vertex(cur_t, 0); | ||
| 339 | index_t v1 = mesh.cells.vertex(cur_t, 1); | ||
| 340 | index_t v2 = mesh.cells.vertex(cur_t, 2); | ||
| 341 | index_t v3 = mesh.cells.vertex(cur_t, 3); | ||
| 342 | |||
| 343 | // TODO: take weights into account | ||
| 344 | // with a new random_point_in_tetra_weighted() | ||
| 345 | // function. | ||
| 346 |
1/2✓ Branch 1 taken 1200 times.
✗ Branch 2 not taken.
|
2400 | Point cur_p = Geom::random_point_in_tetra( |
| 347 | mesh.vertices.point<DIM>(v0), | ||
| 348 | mesh.vertices.point<DIM>(v1), | ||
| 349 | mesh.vertices.point<DIM>(v2), | ||
| 350 | mesh.vertices.point<DIM>(v3) | ||
| 351 | ); | ||
| 352 |
2/2✓ Branch 0 taken 6300 times.
✓ Branch 1 taken 1200 times.
|
15000 | for(coord_index_t coord = 0; coord < DIM; coord++) { |
| 353 | 12600 | p[i * DIM + coord] = cur_p[coord]; | |
| 354 | } | ||
| 355 | } | ||
| 356 |
2/4✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
8 | if(mesh.cells.nb() > 1 && last_t == first_t) { |
| 357 | ✗ | Logger::warn("Sampler") | |
| 358 | << "Did put all the points in the same tetrahedron" | ||
| 359 | << std::endl; | ||
| 360 | ✗ | return false; | |
| 361 | } | ||
| 362 | return true; | ||
| 363 | } | ||
| 364 | } | ||
| 365 | |||
| 366 | #endif | ||
| 367 |