RINGMesh  Version 5.0.0
A programming library for geological model meshes
geometry.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2012-2017, Association Scientifique pour la Geologie et ses
3  * Applications (ASGA). 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  * * Redistributions of source code must retain the above copyright
8  * notice, this list of conditions and the following disclaimer.
9  * * Redistributions in binary form must reproduce the above copyright
10  * notice, this list of conditions and the following disclaimer in the
11  * documentation and/or other materials provided with the distribution.
12  * * Neither the name of ASGA nor the
13  * names of its contributors may be used to endorse or promote products
14  * derived from this software without specific prior written permission.
15  *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
18  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
19  * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ASGA BE LIABLE FOR ANY DIRECT,
20  * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
25  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26  *
27  * http://www.ring-team.org
28  *
29  * RING Project
30  * Ecole Nationale Superieure de Geologie - GeoRessources
31  * 2 Rue du Doyen Marcel Roubault - TSA 70605
32  * 54518 VANDOEUVRE-LES-NANCY
33  * FRANCE
34  */
35 
37 
38 #include <geogram/mesh/mesh.h>
39 
40 #include <geogram/numerics/predicates.h>
41 
49 namespace
50 {
51  using namespace RINGMesh;
52 
53  bool is_almost_zero( double value )
54  {
55  return value < global_epsilon && value > -global_epsilon;
56  }
57 } // namespace
58 
59 namespace RINGMesh
60 {
61  double dot_perp( const vec2& v0, const vec2& v1 )
62  {
63  return dot( v0, vec2( v1.y, -v1.x ) );
64  }
65 
66  double triangle_signed_area( const vec3& p0,
67  const vec3& p1,
68  const vec3& p2,
69  const vec3& triangle_normal )
70  {
71  double area{ GEO::Geom::triangle_area( p0, p1, p2 ) };
72  vec3 area_normal{ cross( p0 - p2, p1 - p2 ) };
73  if( dot( triangle_normal, area_normal ) < 0 )
74  {
75  area = -area;
76  }
77  return area;
78  }
79 
80  bool operator==( const vec3& u, const vec3& v )
81  {
82  return u.x == v.x && u.y == v.y && u.z == v.z;
83  }
84 
85  bool operator!=( const vec3& u, const vec3& v )
86  {
87  return u.x != v.x || u.y != v.y || u.z != v.z;
88  }
89 
90  std::tuple< bool, std::array< double, 4 > > tetra_barycentric_coordinates(
91  const vec3& p,
92  const vec3& p0,
93  const vec3& p1,
94  const vec3& p2,
95  const vec3& p3 )
96  {
97  double total_volume{ GEO::Geom::tetra_signed_volume( p0, p1, p2, p3 ) };
98  if( std::fabs( total_volume ) < global_epsilon_3 )
99  {
100  std::array< double, 4 > lambdas{ { 0., 0., 0., 0. } };
101  return std::make_tuple( false, lambdas );
102  }
103  double volume0{ GEO::Geom::tetra_signed_volume( p1, p3, p2, p ) };
104  double volume1{ GEO::Geom::tetra_signed_volume( p0, p2, p3, p ) };
105  double volume2{ GEO::Geom::tetra_signed_volume( p0, p3, p1, p ) };
106  double volume3{ GEO::Geom::tetra_signed_volume( p0, p1, p2, p ) };
107 
108  double lambda0{ volume0 / total_volume };
109  double lambda1{ volume1 / total_volume };
110  double lambda2{ volume2 / total_volume };
111  double lambda3{ volume3 / total_volume };
112  std::array< double, 4 > lambdas{ { lambda0, lambda1, lambda2,
113  lambda3 } };
114  return std::make_tuple( true, lambdas );
115  }
116 
117  std::tuple< bool, std::array< double, 3 > >
119  const vec3& p, const vec3& p0, const vec3& p1, const vec3& p2 )
120  {
121  double total_area{ GEO::Geom::triangle_area( p0, p1, p2 ) };
122  if( std::fabs( total_area ) < global_epsilon_sq )
123  {
124  std::array< double, 3 > lambdas{ { 0., 0., 0. } };
125  return std::make_tuple( false, lambdas );
126  }
127  vec3 triangle_normal{ cross( p2 - p0, p1 - p0 ) };
128  double area0{ triangle_signed_area( p2, p1, p, triangle_normal ) };
129  double area1{ triangle_signed_area( p0, p2, p, triangle_normal ) };
130  double area2{ triangle_signed_area( p1, p0, p, triangle_normal ) };
131 
132  double lambda0{ area0 / total_area };
133  double lambda1{ area1 / total_area };
134  double lambda2{ area2 / total_area };
135  std::array< double, 3 > lambdas{ { lambda0, lambda1, lambda2 } };
136  return std::make_tuple( true, lambdas );
137  }
138 
139  std::tuple< bool, std::array< double, 3 > >
141  const vec2& p, const vec2& p0, const vec2& p1, const vec2& p2 )
142  {
143  double total_area{ GEO::Geom::triangle_signed_area( p2, p1, p0 ) };
144  if( std::fabs( total_area ) < global_epsilon_sq )
145  {
146  std::array< double, 3 > lambdas{ { 0., 0., 0. } };
147  return std::make_tuple( false, lambdas );
148  }
149  double area0{ GEO::Geom::triangle_signed_area( p2, p1, p ) };
150  double area1{ GEO::Geom::triangle_signed_area( p0, p2, p ) };
151  double area2{ GEO::Geom::triangle_signed_area( p1, p0, p ) };
152 
153  double lambda0{ area0 / total_area };
154  double lambda1{ area1 / total_area };
155  double lambda2{ area2 / total_area };
156  std::array< double, 3 > lambdas{ { lambda0, lambda1, lambda2 } };
157  return std::make_tuple( true, lambdas );
158  }
159 
160  template < index_t DIMENSION >
161  std::tuple< bool, vecn< DIMENSION > > point_segment_projection(
162  const vecn< DIMENSION >& p,
163  const vecn< DIMENSION >& p0,
164  const vecn< DIMENSION >& p1 )
165  {
166  vecn< DIMENSION > center{ ( p0 + p1 ) * 0.5 };
167  vecn< DIMENSION > diff{ p - center };
168  vecn< DIMENSION > edge{ p1 - p0 };
169  double extent{ 0.5 * edge.length() };
170  edge = normalize( edge );
171  double d{ dot( edge, diff ) };
172 
173  if( std::fabs( d ) <= extent )
174  {
175  vecn< DIMENSION > new_p{ center + d * edge };
176  return std::make_tuple( true, new_p );
177  }
178  vecn< DIMENSION > empty_point;
179  return std::make_tuple( false, empty_point );
180  }
181 
182  std::tuple< double, vec3 > point_segment_distance(
183  const vec3& p, const vec3& p0, const vec3& p1 )
184  {
185  bool is_point_segment_projection_possible;
186  vec3 nearest_p;
187  std::tie( is_point_segment_projection_possible, nearest_p ) =
188  point_segment_projection( p, p0, p1 );
189  if( is_point_segment_projection_possible )
190  {
191  return std::make_tuple( length( nearest_p - p ), nearest_p );
192  }
193  double p0_distance_sq{ length2( p0 - p ) };
194  double p1_distance_sq{ length2( p1 - p ) };
195  if( p0_distance_sq < p1_distance_sq )
196  {
197  return std::make_tuple( std::sqrt( p0_distance_sq ), p0 );
198  }
199  return std::make_tuple( std::sqrt( p1_distance_sq ), p1 );
200  }
201 
202  GEO::Matrix< 4, double > rotation_matrix_about_arbitrary_axis(
203  const vec3& origin, const vec3& axis, double theta, bool degrees )
204  {
205  // Note: Rotation is impossible about an axis with null length.
206  ringmesh_assert( axis != vec3() );
207 
208  if( degrees )
209  {
210  theta = theta * M_PI / 180.;
211  }
212 
213  double axis_length{ axis.length() };
214  ringmesh_assert( axis_length > 0. );
215  double x1{ origin[0] };
216  double y1{ origin[1] };
217  double z1{ origin[2] };
218  double a{ axis[0] / axis_length };
219  double b{ axis[1] / axis_length };
220  double c{ axis[2] / axis_length };
221  double d{ std::sqrt( b * b + c * c ) };
222  double cos_angle{ std::cos( theta ) };
223  double sin_angle{ std::sin( theta ) };
224 
225  GEO::Matrix< 4, double > T;
226  T( 0, 0 ) = 1;
227  T( 0, 1 ) = 0;
228  T( 0, 2 ) = 0;
229  T( 0, 3 ) = -x1;
230  T( 1, 0 ) = 0;
231  T( 1, 1 ) = 1;
232  T( 1, 2 ) = 0;
233  T( 1, 3 ) = -y1;
234  T( 2, 0 ) = 0;
235  T( 2, 1 ) = 0;
236  T( 2, 2 ) = 1;
237  T( 2, 3 ) = -z1;
238  T( 3, 0 ) = 0;
239  T( 3, 1 ) = 0;
240  T( 3, 2 ) = 0;
241  T( 3, 3 ) = 1;
242 
243  GEO::Matrix< 4, double > inv_T;
244  inv_T( 0, 0 ) = 1.;
245  inv_T( 0, 1 ) = 0.;
246  inv_T( 0, 2 ) = 0.;
247  inv_T( 0, 3 ) = x1;
248  inv_T( 1, 0 ) = 0.;
249  inv_T( 1, 1 ) = 1.;
250  inv_T( 1, 2 ) = 0.;
251  inv_T( 1, 3 ) = y1;
252  inv_T( 2, 0 ) = 0.;
253  inv_T( 2, 1 ) = 0.;
254  inv_T( 2, 2 ) = 1.;
255  inv_T( 2, 3 ) = z1;
256  inv_T( 3, 0 ) = 0.;
257  inv_T( 3, 1 ) = 0.;
258  inv_T( 3, 2 ) = 0.;
259  inv_T( 3, 3 ) = 1.;
260 
261 #ifdef RINGMESH_DEBUG
262  GEO::Matrix< 4, double > computed_inv_T = T.inverse();
263 #endif
264  ringmesh_assert( inv_T( 0, 0 ) == computed_inv_T( 0, 0 ) );
265  ringmesh_assert( inv_T( 0, 1 ) == computed_inv_T( 0, 1 ) );
266  ringmesh_assert( inv_T( 0, 2 ) == computed_inv_T( 0, 2 ) );
267  ringmesh_assert( inv_T( 0, 3 ) == computed_inv_T( 0, 3 ) );
268  ringmesh_assert( inv_T( 1, 0 ) == computed_inv_T( 1, 0 ) );
269  ringmesh_assert( inv_T( 1, 1 ) == computed_inv_T( 1, 1 ) );
270  ringmesh_assert( inv_T( 1, 2 ) == computed_inv_T( 1, 2 ) );
271  ringmesh_assert( inv_T( 1, 3 ) == computed_inv_T( 1, 3 ) );
272  ringmesh_assert( inv_T( 2, 0 ) == computed_inv_T( 2, 0 ) );
273  ringmesh_assert( inv_T( 2, 1 ) == computed_inv_T( 2, 1 ) );
274  ringmesh_assert( inv_T( 2, 2 ) == computed_inv_T( 2, 2 ) );
275  ringmesh_assert( inv_T( 2, 3 ) == computed_inv_T( 2, 3 ) );
276  ringmesh_assert( inv_T( 3, 0 ) == computed_inv_T( 3, 0 ) );
277  ringmesh_assert( inv_T( 3, 1 ) == computed_inv_T( 3, 1 ) );
278  ringmesh_assert( inv_T( 3, 2 ) == computed_inv_T( 3, 2 ) );
279  ringmesh_assert( inv_T( 3, 3 ) == computed_inv_T( 3, 3 ) );
280 
281  // Note: If d = 0, so rotation is along x axis. So Rx = inv_Rx = Id
282  GEO::Matrix< 4, double > Rx;
283  Rx( 0, 0 ) = 1.;
284  Rx( 0, 1 ) = 0.;
285  Rx( 0, 2 ) = 0.;
286  Rx( 0, 3 ) = 0.;
287  Rx( 1, 0 ) = 0.;
288  Rx( 1, 3 ) = 0.;
289  Rx( 2, 0 ) = 0.;
290  Rx( 2, 3 ) = 0.;
291  Rx( 3, 0 ) = 0.;
292  Rx( 3, 1 ) = 0.;
293  Rx( 3, 2 ) = 0.;
294  Rx( 3, 3 ) = 1.;
295  if( d == 0. )
296  {
297  Rx( 1, 1 ) = 1.;
298  Rx( 1, 2 ) = 0.;
299  Rx( 2, 1 ) = 0.;
300  Rx( 2, 2 ) = 1.;
301  }
302  else
303  {
304  Rx( 1, 1 ) = c / d;
305  Rx( 1, 2 ) = -b / d;
306  Rx( 2, 1 ) = b / d;
307  Rx( 2, 2 ) = c / d;
308  }
309 
310  GEO::Matrix< 4, double > inv_Rx;
311  inv_Rx( 0, 0 ) = 1.;
312  inv_Rx( 0, 1 ) = 0.;
313  inv_Rx( 0, 2 ) = 0.;
314  inv_Rx( 0, 3 ) = 0.;
315  inv_Rx( 1, 0 ) = 0.;
316  inv_Rx( 1, 3 ) = 0.;
317  inv_Rx( 2, 0 ) = 0.;
318  inv_Rx( 2, 3 ) = 0.;
319  inv_Rx( 3, 0 ) = 0.;
320  inv_Rx( 3, 1 ) = 0.;
321  inv_Rx( 3, 2 ) = 0.;
322  inv_Rx( 3, 3 ) = 1.;
323  if( d == 0. )
324  {
325  inv_Rx( 1, 1 ) = 1.;
326  inv_Rx( 1, 2 ) = 0.;
327  inv_Rx( 2, 1 ) = 0.;
328  inv_Rx( 2, 2 ) = 1.;
329  }
330  else
331  {
332  inv_Rx( 1, 1 ) = c / d;
333  inv_Rx( 1, 2 ) = b / d;
334  inv_Rx( 2, 1 ) = -b / d;
335  inv_Rx( 2, 2 ) = c / d;
336  }
337 
338 #ifdef RINGMESH_DEBUG
339  GEO::Matrix< 4, double > computed_inv_Rx = Rx.inverse();
340 #endif
341  ringmesh_assert( inv_Rx( 0, 0 ) == computed_inv_Rx( 0, 0 ) );
342  ringmesh_assert( inv_Rx( 0, 1 ) == computed_inv_Rx( 0, 1 ) );
343  ringmesh_assert( inv_Rx( 0, 2 ) == computed_inv_Rx( 0, 2 ) );
344  ringmesh_assert( inv_Rx( 0, 3 ) == computed_inv_Rx( 0, 3 ) );
345  ringmesh_assert( inv_Rx( 1, 0 ) == computed_inv_Rx( 1, 0 ) );
346  ringmesh_assert( inv_Rx( 1, 1 ) == computed_inv_Rx( 1, 1 ) );
347  ringmesh_assert( inv_Rx( 1, 2 ) == computed_inv_Rx( 1, 2 ) );
348  ringmesh_assert( inv_Rx( 1, 3 ) == computed_inv_Rx( 1, 3 ) );
349  ringmesh_assert( inv_Rx( 2, 0 ) == computed_inv_Rx( 2, 0 ) );
350  ringmesh_assert( inv_Rx( 2, 1 ) == computed_inv_Rx( 2, 1 ) );
351  ringmesh_assert( inv_Rx( 2, 2 ) == computed_inv_Rx( 2, 2 ) );
352  ringmesh_assert( inv_Rx( 2, 3 ) == computed_inv_Rx( 2, 3 ) );
353  ringmesh_assert( inv_Rx( 3, 0 ) == computed_inv_Rx( 3, 0 ) );
354  ringmesh_assert( inv_Rx( 3, 1 ) == computed_inv_Rx( 3, 1 ) );
355  ringmesh_assert( inv_Rx( 3, 2 ) == computed_inv_Rx( 3, 2 ) );
356  ringmesh_assert( inv_Rx( 3, 3 ) == computed_inv_Rx( 3, 3 ) );
357 
358  GEO::Matrix< 4, double > Ry;
359  Ry( 0, 0 ) = d;
360  Ry( 0, 1 ) = 0.;
361  Ry( 0, 2 ) = -a;
362  Ry( 0, 3 ) = 0.;
363  Ry( 1, 0 ) = 0.;
364  Ry( 1, 1 ) = 1.;
365  Ry( 1, 2 ) = 0.;
366  Ry( 1, 3 ) = 0.;
367  Ry( 2, 0 ) = a;
368  Ry( 2, 1 ) = 0.;
369  Ry( 2, 2 ) = d;
370  Ry( 2, 3 ) = 0.;
371  Ry( 3, 0 ) = 0.;
372  Ry( 3, 1 ) = 0.;
373  Ry( 3, 2 ) = 0.;
374  Ry( 3, 3 ) = 1.;
375 
376  GEO::Matrix< 4, double > inv_Ry;
377  inv_Ry( 0, 0 ) = d;
378  inv_Ry( 0, 1 ) = 0.;
379  inv_Ry( 0, 2 ) = a;
380  inv_Ry( 0, 3 ) = 0.;
381  inv_Ry( 1, 0 ) = 0.;
382  inv_Ry( 1, 1 ) = 1.;
383  inv_Ry( 1, 2 ) = 0.;
384  inv_Ry( 1, 3 ) = 0.;
385  inv_Ry( 2, 0 ) = -a;
386  inv_Ry( 2, 1 ) = 0.;
387  inv_Ry( 2, 2 ) = d;
388  inv_Ry( 2, 3 ) = 0.;
389  inv_Ry( 3, 0 ) = 0.;
390  inv_Ry( 3, 1 ) = 0.;
391  inv_Ry( 3, 2 ) = 0.;
392  inv_Ry( 3, 3 ) = 1.;
393 
394 #ifdef RINGMESH_DEBUG
395  GEO::Matrix< 4, double > computed_inv_Ry = Ry.inverse();
396 #endif
397  ringmesh_assert( inv_Ry( 0, 0 ) == computed_inv_Ry( 0, 0 ) );
398  ringmesh_assert( inv_Ry( 0, 1 ) == computed_inv_Ry( 0, 1 ) );
399  ringmesh_assert( inv_Ry( 0, 2 ) == computed_inv_Ry( 0, 2 ) );
400  ringmesh_assert( inv_Ry( 0, 3 ) == computed_inv_Ry( 0, 3 ) );
401  ringmesh_assert( inv_Ry( 1, 0 ) == computed_inv_Ry( 1, 0 ) );
402  ringmesh_assert( inv_Ry( 1, 1 ) == computed_inv_Ry( 1, 1 ) );
403  ringmesh_assert( inv_Ry( 1, 2 ) == computed_inv_Ry( 1, 2 ) );
404  ringmesh_assert( inv_Ry( 1, 3 ) == computed_inv_Ry( 1, 3 ) );
405  ringmesh_assert( inv_Ry( 2, 0 ) == computed_inv_Ry( 2, 0 ) );
406  ringmesh_assert( inv_Ry( 2, 1 ) == computed_inv_Ry( 2, 1 ) );
407  ringmesh_assert( inv_Ry( 2, 2 ) == computed_inv_Ry( 2, 2 ) );
408  ringmesh_assert( inv_Ry( 2, 3 ) == computed_inv_Ry( 2, 3 ) );
409  ringmesh_assert( inv_Ry( 3, 0 ) == computed_inv_Ry( 3, 0 ) );
410  ringmesh_assert( inv_Ry( 3, 1 ) == computed_inv_Ry( 3, 1 ) );
411  ringmesh_assert( inv_Ry( 3, 2 ) == computed_inv_Ry( 3, 2 ) );
412  ringmesh_assert( inv_Ry( 3, 3 ) == computed_inv_Ry( 3, 3 ) );
413 
414  GEO::Matrix< 4, double > Rz;
415  Rz( 0, 0 ) = cos_angle;
416  Rz( 0, 1 ) = -sin_angle;
417  Rz( 0, 2 ) = 0.;
418  Rz( 0, 3 ) = 0.;
419  Rz( 1, 0 ) = sin_angle;
420  Rz( 1, 1 ) = cos_angle;
421  Rz( 1, 2 ) = 0.;
422  Rz( 1, 3 ) = 0.;
423  Rz( 2, 0 ) = 0.;
424  Rz( 2, 1 ) = 0.;
425  Rz( 2, 2 ) = 1.;
426  Rz( 2, 3 ) = 0.;
427  Rz( 3, 0 ) = 0.;
428  Rz( 3, 1 ) = 0.;
429  Rz( 3, 2 ) = 0.;
430  Rz( 3, 3 ) = 1.;
431 
432  return inv_T * inv_Rx * inv_Ry * Rz * Ry * Rx * T;
433  }
434 
435  template std::tuple< bool, vecn< 2 > >
436  RINGMESH_API point_segment_projection(
437  const vecn< 2 >&, const vecn< 2 >&, const vecn< 2 >& );
438 
439  template std::tuple< bool, vecn< 3 > >
440  RINGMESH_API point_segment_projection(
441  const vecn< 3 >&, const vecn< 3 >&, const vecn< 3 >& );
442 } // namespace RINGMesh
std::tuple< bool, vecn< DIMENSION > > point_segment_projection(const vecn< DIMENSION > &p, const vecn< DIMENSION > &p0, const vecn< DIMENSION > &p1)
Definition: geometry.cpp:161
std::tuple< double, vec3 > point_segment_distance(const vec3 &p, const vec3 &p0, const vec3 &p1)
Definition: geometry.cpp:182
GEO::vecng< DIMENSION, double > vecn
Definition: types.h:74
vecn< 3 > vec3
Definition: types.h:76
double RINGMESH_API triangle_signed_area(const vec3 &p0, const vec3 &p1, const vec3 &p2, const vec3 &triangle_normal)
Definition: geometry.cpp:66
double RINGMESH_API dot_perp(const vec2 &v0, const vec2 &v1)
Definition: geometry.cpp:61
std::tuple< bool, std::array< double, 4 > > RINGMESH_API tetra_barycentric_coordinates(const vec3 &p, const vec3 &p0, const vec3 &p1, const vec3 &p2, const vec3 &p3)
Definition: geometry.cpp:90
bool operator==(const vecn< DIMENSION > &u, const vecn< DIMENSION > &v)
Definition: geometry.h:50
bool operator!=(const vecn< DIMENSION > &u, const vecn< DIMENSION > &v)
Definition: geometry.h:61
std::tuple< bool, std::array< double, 3 > > RINGMESH_API triangle_barycentric_coordinates(const vec3 &p, const vec3 &p0, const vec3 &p1, const vec3 &p2)
Definition: geometry.cpp:118
vecn< 2 > vec2
Definition: types.h:78
GEO::Matrix< 4, double > RINGMESH_API rotation_matrix_about_arbitrary_axis(const vec3 &origin, const vec3 &axis, double theta, bool degrees)
Builds a rotational matrix about an arbitrary axis.
Definition: geometry.cpp:202
#define ringmesh_assert(x)
Classes to build GeoModel from various inputs.
Definition: algorithm.h:48