Geogram  Version 1.9.1-rc
A programming library of geometric algorithms
exact_geometry.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2000-2023 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_NUMERICS_EXACT_GEOMETRY
41 #define GEOGRAM_NUMERICS_EXACT_GEOMETRY
42 
43 #include <geogram/basic/common.h>
44 #include <geogram/basic/geometry.h>
45 #include <geogram/basic/vechg.h>
47 #include <geogram/numerics/interval_nt.h>
48 
51 
52 #ifdef GEOGRAM_WITH_GEOGRAMPLUS
53 #include <geogram/geogramplus/numerics/exact_geometry.h>
54 #endif
55 
66 // If Tessael's geogramplus is available, use exact_nt coordinates,
67 // else use expansion_nt coordinates.
68 // exact_nt coordinates makes the algorithm 10x to 20x faster
69 // and have no risk of underflow / overflow.
70 #ifdef GEOGRAM_WITH_GEOGRAMPLUS
71 #define GEOGRAM_USE_EXACT_NT
72 #endif
73 
74 namespace GEO {
75 
81 
87 
94 
101 
109 
117 
125 
133 
134  /***********************************************************************/
135 
143  template <class VEC3 = vec3>
144  inline VEC3 make_vec3(const vec3& p1, const vec3& p2) {
145  typedef typename VEC3::value_type value_type;
146  return VEC3(
147  value_type(p2.x) - value_type(p1.x),
148  value_type(p2.y) - value_type(p1.y),
149  value_type(p2.z) - value_type(p1.z)
150  );
151  }
152 
160  template <class VEC2>
161  inline VEC2 make_vec2(
162  const vec2& p1, const vec2& p2
163  ) {
164  typedef typename VEC2::value_type value_type;
165  return VEC2(
166  value_type(p2.x) - value_type(p1.x),
167  value_type(p2.y) - value_type(p1.y)
168  );
169  }
170 
179  template <class VEC3>
180  inline VEC3 triangle_normal(
181  const vec3& p1, const vec3& p2, const vec3& p3
182  ) {
183  return cross(
184  make_vec3<VEC3>(p1,p2),
185  make_vec3<VEC3>(p1,p3)
186  );
187  }
188 
189  /***********************************************************************/
190 
191  namespace PCK {
192 
204  Sign GEOGRAM_API orient_2d(
205  const vec2HE& p0, const vec2HE& p1, const vec2HE& p2
206  );
207 
224  const vec3HE& p0, const vec3HE& p1, const vec3HE& p2,
225  coord_index_t axis
226  );
227 
239  Sign GEOGRAM_API orient_3d(
240  const vec3HE& p0, const vec3HE& p1,
241  const vec3HE& p2, const vec3HE& p3
242  );
243 
252  Sign GEOGRAM_API dot_2d(
253  const vec2HE& p0, const vec2HE& p1, const vec2HE& p2
254  );
255 
272  const vec2HE& p0, const vec2HE& p1,
273  const vec2HE& p2, const vec2HE& p3,
274  double l0, double l1, double l2, double l3
275  );
276 
293  const vec2HE& p0, const vec2HE& p1,
294  const vec2HE& p2, const vec2HE& p3,
295  double l0, double l1, double l2, double l3
296  );
297 
315  const vec2HE& p0, const vec2HE& p1,
316  const vec2HE& p2, const vec2HE& p3
317  ) {
318  double l0 = (geo_sqr(p0.x) + geo_sqr(p0.y)).estimate() /
319  geo_sqr(p0.w).estimate();
320  double l1 = (geo_sqr(p1.x) + geo_sqr(p1.y)).estimate() /
321  geo_sqr(p1.w).estimate();
322  double l2 = (geo_sqr(p2.x) + geo_sqr(p2.y)).estimate() /
323  geo_sqr(p2.w).estimate();
324  double l3 = (geo_sqr(p3.x) + geo_sqr(p3.y)).estimate() /
325  geo_sqr(p3.w).estimate();
326  return incircle_2d_SOS_with_lengths(p0,p1,p2,p3,l0,l1,l2,l3);
327  }
328 
340  const vec3& p1, const vec3& p2, const vec3& p3
341  );
342 
351  bool GEOGRAM_API aligned_3d(
352  const vec3HE& p0, const vec3HE& p1, const vec3HE& p2
353  );
354 
363  bool GEOGRAM_API on_segment_3d(
364  const vec3HE& p, const vec3HE& q1, const vec3HE& q2
365  );
366 
374  vec3 GEOGRAM_API approximate(const vec3HE& p);
375 
383  vec2 GEOGRAM_API approximate(const vec2HE& p);
384 
385  }
386 
387  /************************************************************************/
388 
392  template <>
393  inline vec2E make_vec2<vec2E>(const vec2& p1, const vec2& p2) {
394  return vec2E(
395  expansion_nt(expansion_nt::DIFF, p2.x, p1.x),
396  expansion_nt(expansion_nt::DIFF, p2.y, p1.y)
397  );
398  }
399 
403  template <>
404  inline vec3E make_vec3<vec3E>(const vec3& p1, const vec3& p2) {
405  return vec3E(
406  expansion_nt(expansion_nt::DIFF, p2.x, p1.x),
407  expansion_nt(expansion_nt::DIFF, p2.y, p1.y),
408  expansion_nt(expansion_nt::DIFF, p2.z, p1.z)
409  );
410  }
411 
412 // Under Linux we got 10 Mb of stack (!) Then some operations can be
413 // made faster by using the low-level expansion API (that allocates
414 // intermediary multiprecision values on stack rather than in the heap).
415 // These optimized functions are written as template specializations
416 // (used automatically).
417 
418 #ifdef GEO_HAS_BIG_STACK
419 
423  template<> expansion_nt GEOGRAM_API det(const vec2E& v1, const vec2E& v2);
424 
428  template<> expansion_nt GEOGRAM_API dot(const vec2E& v1, const vec2E& v2);
429 
433  template<> expansion_nt GEOGRAM_API dot(const vec3E& v1, const vec3E& v2);
434 
438  template<> vec2Hg<expansion_nt> GEOGRAM_API mix(
439  const rationalg<expansion_nt>& t,
440  const vecng<2,double>& p1, const vecng<2,double>& p2
441  );
442 
446  template<> vec3Hg<expansion_nt> GEOGRAM_API mix(
447  const rationalg<expansion_nt>& t,
448  const vecng<3,double>& p1, const vecng<3,double>& p2
449  );
450 
454  template <> GEOGRAM_API vec3E triangle_normal<vec3E>(
455  const vec3& p1, const vec3& p2, const vec3& p3
456  );
457 
458 #endif
459 
460  /************************************************************************/
461 
467  namespace exact {
468 #ifdef GEOGRAM_USE_EXACT_NT
469  typedef exact_nt scalar;
470 #else
472 #endif
480 
485 
490  }
491 }
492 
493 #endif
Expansion_nt (expansion Number Type) is used to compute the sign of polynoms exactly.
Definition: expansion_nt.h:70
rationalg (generic rational) is used to compute the sign of rational fractions exactly.
Definition: rationalg.h:59
2d vector with homogeneous coordinates
Definition: vechg.h:60
3d vector with homogeneous coordinates
Definition: vechg.h:185
Generic maths vector.
Definition: vecg.h:70
Exact predicates and constructs.
High-level interface to multi-precision arithmetics.
Common include file, providing basic definitions. Should be included before anything else by all head...
Geometric functions in 2d and 3d.
Sign orient_2d_projected(const vec3HE &p0, const vec3HE &p1, const vec3HE &p2, coord_index_t axis)
Computes the orientation predicate in 2d projected along an axis.
Sign orient_2d(const vec2HE &p0, const vec2HE &p1, const vec2HE &p2)
Computes the orientation predicate in 2d.
Sign incircle_2d_SOS_with_lengths(const vec2HE &p0, const vec2HE &p1, const vec2HE &p2, const vec2HE &p3, double l0, double l1, double l2, double l3)
Tests whether a point is in the circumscribed circle of three other points.
bool on_segment_3d(const vec3HE &p, const vec3HE &q1, const vec3HE &q2)
Tests whether a point is on a segment.
vec3 approximate(const vec3HE &p)
Gets a 3D floating-point approximation of a 3D point with exact coordinates.
coord_index_t triangle_normal_axis(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Gets the axis that is most normal to a triangle.
Sign dot_2d(const vec2HE &p0, const vec2HE &p1, const vec2HE &p2)
Computes the sign of the dot product between two vectors defined by three points.
bool aligned_3d(const vec3HE &p0, const vec3HE &p1, const vec3HE &p2)
Tests whether three 3d points are aligned.
Sign orient_3d(const vec3HE &p0, const vec3HE &p1, const vec3HE &p2, const vec3HE &p3)
Computes the orientation predicate in 3d.
Sign incircle_2d_SOS(const vec2HE &p0, const vec2HE &p1, const vec2HE &p2, const vec2HE &p3)
Tests whether a point is in the circumscribed circle of three other points.
vecng< 2, scalar > vec2
expansion_nt scalar
vec2Hg< scalar > vec2h
2d vector with exact homogeneous coordinates
vecng< 3, scalar > vec3
rationalg< scalar > rational
rational with exact numerator and denominator
vec3Hg< scalar > vec3h
3d vector with exact homogeneous coordinates
Global Vorpaline namespace.
Definition: basic.h:55
vecng< 3, interval_nt > vec3I
vec3 with coordinates as interval_nt
vec2Hg< expansion_nt > vec2HE
2D vector in homogeneous coordinates with coordinates as expansions
vecng< 2, expansion_nt > vec2E
vec2 with coordinates as expansions
T geo_sqr(T x)
Gets the square value of a value.
Definition: numeric.h:301
vec2Hg< interval_nt > vec2HI
2D vector in homogeneous coordinates with coordinates as intervals.
VEC3 make_vec3(const vec3 &p1, const vec3 &p2)
Creates a vector with coordinates of arbitrary type from two points with double coordinates.
vecng< 2, interval_nt > vec2I
vec2 with coordinates as interval_nt
vecng< 3, Numeric::float64 > vec3
Represents points and vectors in 3d.
Definition: geometry.h:65
vec3E make_vec3< vec3E >(const vec3 &p1, const vec3 &p2)
Specialization of make_vec3() for vec3E.
double det(const mat2 &M)
Computes the determinant of a 2x2 matrix.
Definition: geometry.h:142
vec3Hg< expansion_nt > vec3HE
3D vector in homogeneous coordinates with coordinates as expansions
vec2E make_vec2< vec2E >(const vec2 &p1, const vec2 &p2)
Specialization of make_vec2() for vec2E.
Sign
Integer constants that represent the sign of a value.
Definition: numeric.h:68
vec3Hg< interval_nt > vec3HI
3D vector in homogeneous coordinates with coordinates as intervals.
VEC2 make_vec2(const vec2 &p1, const vec2 &p2)
Creates a vector with coordinates of arbitrary type from two points with double coordinates.
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
VEC3 triangle_normal(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Computes the normal to a triangle from its three vertices.
vecng< 3, expansion_nt > vec3E
vec3 with coordinates as expansions
Filtered exact predicates for restricted Voronoi diagrams.
Generic implementation of geometric vectors in homogeneous coordinates.