Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
mesh_sampling.h
Go to the documentation of this file.
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>
47 #include <geogram/basic/logger.h>
48 #include <algorithm>
49 
56 namespace GEO {
57 
68  template <index_t DIM>
69  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  typedef vecng<DIM, double> Point;
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  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v1)),
84  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v2)),
85  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v3)),
86  vertex_weight[v1],
87  vertex_weight[v2],
88  vertex_weight[v3]
89  );
90  }
91  return Geom::triangle_area(
92  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v1)),
93  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v2)),
94  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v3))
95  );
96  }
97 
119  template <index_t DIM>
121  const Mesh& mesh,
122  double* p,
123  index_t nb_points,
124  Attribute<double>& weight,
125  signed_index_t facets_begin_in = -1,
126  signed_index_t facets_end_in = -1
127  ) {
128  geo_assert(mesh.facets.are_simplices());
129  geo_assert(mesh.vertices.dimension() >= DIM);
130  geo_assert(mesh.facets.nb() > 0);
131 
132  index_t facets_begin = 0;
133  index_t facets_end = mesh.facets.nb();
134  if(facets_begin_in != -1) {
135  facets_begin = index_t(facets_begin_in);
136  }
137  if(facets_end_in != -1) {
138  facets_end = index_t(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.
146 
147  vector<double> s(nb_points);
148  for(index_t i = 0; i < nb_points; i++) {
149  s[i] = Numeric::random_float64();
150  }
151  std::sort(s.begin(), s.end());
152 
153  double Atot = 0.0;
154  for(index_t t = facets_begin; t < facets_end; ++t) {
155  double At = mesh_facet_mass<DIM>(mesh, t, weight);
156  Atot += At;
157  }
158 
159  signed_index_t first_t = -1;
160  signed_index_t last_t = 0;
161 
162  index_t cur_t = facets_begin;
163  double cur_s =
164  mesh_facet_mass<DIM>(mesh, facets_begin, weight) / Atot;
165  for(index_t i = 0; i < nb_points; i++) {
166  geo_debug_assert(i < s.size());
167  while(s[i] > cur_s && cur_t < facets_end - 1) {
168  cur_t++;
169  geo_debug_assert(cur_t < facets_end);
170  cur_s += mesh_facet_mass<DIM>(mesh, cur_t, weight) / Atot;
171  }
172  if(first_t == -1) {
173  first_t = signed_index_t(cur_t);
174  }
175  last_t = std::max(last_t, signed_index_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  Point cur_p = Geom::random_point_in_triangle(
184  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v1)),
185  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v2)),
186  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v3))
187  );
188  for(coord_index_t coord = 0; coord < DIM; coord++) {
189  p[i * DIM + coord] = cur_p[coord];
190  }
191  }
192  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 
211  template <index_t DIM>
212  inline double mesh_tetra_mass(
213  const Mesh& mesh,
214  index_t t
215  ) {
216  geo_debug_assert(mesh.vertices.dimension() >= DIM);
217  typedef vecng<DIM, double> Point;
218 
219  index_t v0 = mesh.cells.vertex(t, 0);
220  index_t v1 = mesh.cells.vertex(t, 1);
221  index_t v2 = mesh.cells.vertex(t, 2);
222  index_t v3 = mesh.cells.vertex(t, 3);
223 
224  double result = Geom::tetra_volume(
225  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v0)),
226  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v1)),
227  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v2)),
228  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v3))
229  );
230 
231  return result;
232  }
233 
244  template <index_t DIM>
245  inline double mesh_tetra_mass(
246  const Mesh& mesh,
247  index_t t,
248  const Attribute<double>& weight
249  ) {
250  double result = mesh_tetra_mass<DIM>(mesh, t);
251 
252  if(weight.is_bound()) {
253  index_t v0 = mesh.cells.vertex(t, 0);
254  index_t v1 = mesh.cells.vertex(t, 1);
255  index_t v2 = mesh.cells.vertex(t, 2);
256  index_t v3 = mesh.cells.vertex(t, 3);
257  result *= (
258  weight[v0] + weight[v1] +
259  weight[v2] + weight[v3]
260  ) / 4.0;
261  }
262 
263  return result;
264  }
265 
286  template <index_t DIM>
288  const Mesh& mesh,
289  double* p,
290  index_t nb_points,
291  Attribute<double>& vertex_weight,
292  signed_index_t tets_begin_in = -1,
293  signed_index_t tets_end_in = -1
294  ) {
295  geo_assert(mesh.vertices.dimension() >= DIM);
296  geo_assert(mesh.cells.nb() > 0);
297 
298  index_t tets_begin = 0;
299  index_t tets_end = mesh.cells.nb();
300  if(tets_begin_in != -1) {
301  tets_begin = index_t(tets_begin_in);
302  }
303  if(tets_end_in != -1) {
304  tets_end = index_t(tets_end_in);
305  }
306 
307  typedef vecng<DIM, double> Point;
308 
309  // To ensure reproducibility accros successive
310  // runs, reset the random number generator.
312 
313  vector<double> s(nb_points);
314  for(index_t i = 0; i < nb_points; i++) {
315  s[i] = Numeric::random_float64();
316  }
317  std::sort(s.begin(), s.end());
318 
319  double Vtot = 0.0;
320  for(index_t t = tets_begin; t < tets_end; ++t) {
321  double Vt = mesh_tetra_mass<DIM>(mesh, t, vertex_weight);
322  Vtot += Vt;
323  }
324 
325  signed_index_t first_t = -1;
326  signed_index_t last_t = 0;
327 
328  index_t cur_t = tets_begin;
329  double cur_s =
330  mesh_tetra_mass<DIM>(mesh, tets_begin, vertex_weight) / Vtot;
331  for(index_t i = 0; i < nb_points; i++) {
332  geo_debug_assert(i < s.size());
333  while(s[i] > cur_s && cur_t < tets_end - 1) {
334  cur_t++;
335  geo_debug_assert(cur_t < tets_end);
336  cur_s += mesh_tetra_mass<DIM>(
337  mesh, cur_t, vertex_weight
338  ) / Vtot;
339  }
340  if(first_t == -1) {
341  first_t = signed_index_t(cur_t);
342  }
343  last_t = std::max(last_t, signed_index_t(cur_t));
344 
345  index_t v0 = mesh.cells.vertex(cur_t, 0);
346  index_t v1 = mesh.cells.vertex(cur_t, 1);
347  index_t v2 = mesh.cells.vertex(cur_t, 2);
348  index_t v3 = mesh.cells.vertex(cur_t, 3);
349 
350  // TODO: take weights into account
351  // with a new random_point_in_tetra_weighted()
352  // function.
353  Point cur_p = Geom::random_point_in_tetra(
354  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v0)),
355  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v1)),
356  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v2)),
357  *reinterpret_cast<const Point*>(mesh.vertices.point_ptr(v3))
358  );
359  for(coord_index_t coord = 0; coord < DIM; coord++) {
360  p[i * DIM + coord] = cur_p[coord];
361  }
362  }
363  if(mesh.cells.nb() > 1 && last_t == first_t) {
364  Logger::warn("Sampler")
365  << "Did put all the points in the same tetrahedron"
366  << std::endl;
367  return false;
368  }
369  return true;
370  }
371 }
372 
373 #endif
#define geo_assert(x)
Verifies that a condition is met.
Definition: assert.h:149
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
bool is_bound() const
Tests whether an Attribute is bound.
Definition: attributes.h:1139
static std::ostream & warn(const std::string &feature)
Gets the stream to send warning messages.
index_t vertex(index_t c, index_t lv) const
Gets a vertex of a cell by local vertex index.
Definition: mesh.h:1950
bool are_simplices() const
Tests whether all the facets are triangles.
Definition: mesh.h:809
index_t vertex(index_t f, index_t lv) const
Gets a vertex by facet and local vertex index.
Definition: mesh.h:1054
index_t nb() const
Gets the number of (sub-)elements.
Definition: mesh.h:89
const double * point_ptr(index_t v) const
Gets a point.
Definition: mesh.h:455
index_t dimension() const
Gets the dimension of the vertices.
Definition: mesh.h:427
Represents a mesh.
Definition: mesh.h:2693
Generic maths vector.
Definition: vecg.h:70
index_t size() const
Gets the number of elements.
Definition: memory.h:674
Common include file, providing basic definitions. Should be included before anything else by all head...
Geometric functions in arbitrary dimension.
Generic logging mechanism.
The class that represents a mesh.
Functions for accessing the geometry in a mesh.
vec3 random_point_in_triangle(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Generates a random point in a 3d triangle.
Definition: geometry.h:582
double tetra_volume(const vec3 &p1, const vec3 &p2, const vec3 &p3, const vec3 &p4)
Computes the volume of a 3d tetrahedron.
Definition: geometry.h:500
double triangle_mass(const vec3 &p, const vec3 &q, const vec3 &r, double a, double b, double c)
Computes the mass of a 3d triangle with weighted points.
Definition: geometry.h:563
double triangle_area(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Computes the area of a 3d triangle.
Definition: geometry.h:348
VEC random_point_in_tetra(const VEC &p1, const VEC &p2, const VEC &p3, const VEC &p4)
Generates a random point in a nd tetrahedron.
Definition: geometry_nd.h:364
float64 random_float64()
Returns a 64 bits float between 0 and 1.
void random_reset()
Resets the random number generator.
Global Vorpaline namespace.
Definition: basic.h:55
bool mesh_generate_random_samples_in_volume(const Mesh &mesh, double *p, index_t nb_points, Attribute< double > &vertex_weight, signed_index_t tets_begin_in=-1, signed_index_t tets_end_in=-1)
Generates a set of random samples in a volumetric mesh.
geo_signed_index_t signed_index_t
The type for storing and manipulating indices differences.
Definition: numeric.h:343
double mesh_tetra_mass(const Mesh &mesh, index_t t)
Computes the mass of a mesh tetrahedron.
bool mesh_generate_random_samples_on_surface(const Mesh &mesh, double *p, index_t nb_points, Attribute< double > &weight, signed_index_t facets_begin_in=-1, signed_index_t facets_end_in=-1)
Generates a set of random samples over a surfacic mesh.
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
void sort(const ITERATOR &begin, const ITERATOR &end)
Sorts elements in parallel.
Definition: algorithm.h:90
double mesh_facet_mass(const Mesh &mesh, index_t f, Attribute< double > &vertex_weight)
Computes the mass of a mesh facet.
Definition: mesh_sampling.h:69
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition: numeric.h:363