Geogram Version 1.8.5
A programming library of geometric algorithms
Loading...
Searching...
No Matches
geometry.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_BASIC_GEOMETRY
41#define GEOGRAM_BASIC_GEOMETRY
42
45
51namespace GEO {
52
53 /************************************************************************/
54
60
66
72
79
86
93
94
101
108
115
116
122
128
134
135 /************************************************************************/
136
140 namespace Geom {
141
148 inline vec3 barycenter(const vec3& p1, const vec3& p2) {
149 return vec3(
150 0.5 * (p1.x + p2.x),
151 0.5 * (p1.y + p2.y),
152 0.5 * (p1.z + p2.z)
153 );
154 }
155
162 inline vec2 barycenter(const vec2& p1, const vec2& p2) {
163 return vec2(
164 0.5 * (p1.x + p2.x),
165 0.5 * (p1.y + p2.y)
166 );
167 }
168
177 const vec3& p1, const vec3& p2, const vec3& p3
178 ) {
179 return vec3(
180 (p1.x + p2.x + p3.x) / 3.0,
181 (p1.y + p2.y + p3.y) / 3.0,
182 (p1.z + p2.z + p3.z) / 3.0
183 );
184 }
185
193 inline vec2 barycenter(
194 const vec2& p1, const vec2& p2, const vec2& p3
195 ) {
196 return vec2(
197 (p1.x + p2.x + p3.x) / 3.0,
198 (p1.y + p2.y + p3.y) / 3.0
199 );
200 }
201
208 inline double cos_angle(const vec3& a, const vec3& b) {
209 double lab = ::sqrt(length2(a)*length2(b));
210 double result = (lab > 1e-50) ? (dot(a, b) / lab) : 1.0;
211 // Numerical precision problem may occur, and generate
212 // normalized dot products that are outside the valid
213 // range of acos.
214 geo_clamp(result, -1.0, 1.0);
215 return result;
216 }
217
225 inline double angle(const vec3& a, const vec3& b) {
226 return ::acos(cos_angle(a, b));
227 }
228
235 inline double cos_angle(const vec2& a, const vec2& b) {
236 double lab = ::sqrt(length2(a)*length2(b));
237 double result = (lab > 1e-20) ? (dot(a, b) / lab) : 1.0;
238 // Numerical precision problem may occur, and generate
239 // normalized dot products that are outside the valid
240 // range of acos.
241 geo_clamp(result, -1.0, 1.0);
242 return result;
243 }
244
251 inline double det(const vec2& a, const vec2& b) {
252 return a.x * b.y - a.y * b.x;
253 }
254
262 inline double angle(const vec2& a, const vec2& b) {
263 return det(a, b) > 0 ?
264 ::acos(cos_angle(a, b)) :
265 -::acos(cos_angle(a, b));
266 }
267
275 const vec3& p1, const vec3& p2, const vec3& p3
276 ) {
277 return cross(p2 - p1, p3 - p1);
278 }
279
288 const vec3& p1, const vec3& p2, const vec3& p3
289 );
290
291
297 inline double triangle_area_3d(
298 const double* p1, const double* p2, const double* p3
299 ) {
300 double Ux = p2[0] - p1[0];
301 double Uy = p2[1] - p1[1];
302 double Uz = p2[2] - p1[2];
303
304 double Vx = p3[0] - p1[0];
305 double Vy = p3[1] - p1[1];
306 double Vz = p3[2] - p1[2];
307
308 double Nx = Uy*Vz - Uz*Vy;
309 double Ny = Uz*Vx - Ux*Vz;
310 double Nz = Ux*Vy - Uy*Vx;
311 return 0.5 * ::sqrt(Nx*Nx+Ny*Ny+Nz*Nz);
312 }
313
319 inline double triangle_area(
320 const vec3& p1, const vec3& p2, const vec3& p3
321 ) {
322 return triangle_area_3d(p1.data(), p2.data(), p3.data());
323 }
324
334 const double* p1, const double* p2, const double* p3
335 ) {
336 double a = p2[0]-p1[0];
337 double b = p3[0]-p1[0];
338 double c = p2[1]-p1[1];
339 double d = p3[1]-p1[1];
340 return 0.5*(a*d-b*c);
341 }
342
351 inline double triangle_signed_area(
352 const vec2& p1, const vec2& p2, const vec2& p3
353 ) {
354 return 0.5 * det(p2 - p1, p3 - p1);
355 }
356
364 inline double triangle_area(
365 const vec2& p1, const vec2& p2, const vec2& p3
366 ) {
367 return ::fabs(triangle_signed_area(p1, p2, p3));
368 }
369
377 inline double triangle_area_2d(
378 const double* p1, const double* p2, const double* p3
379 ) {
380 return ::fabs(triangle_signed_area_2d(p1,p2,p3));
381 }
382
392 const vec2& p1, const vec2& p2, const vec2& p3
393 );
394
400 inline bool has_nan(const vec3& v) {
401 return
402 Numeric::is_nan(v.x) ||
403 Numeric::is_nan(v.y) ||
404 Numeric::is_nan(v.z);
405 }
406
412 inline bool has_nan(const vec2& v) {
413 return
414 Numeric::is_nan(v.x) ||
415 Numeric::is_nan(v.y);
416 }
417
423 vec3 GEOGRAM_API perpendicular(const vec3& V);
424
434 inline double tetra_signed_volume(
435 const vec3& p1, const vec3& p2,
436 const vec3& p3, const vec3& p4
437 ) {
438 return dot(p2 - p1, cross(p3 - p1, p4 - p1)) / 6.0;
439 }
440
450 inline double tetra_signed_volume(
451 const double* p1, const double* p2,
452 const double* p3, const double* p4
453 ) {
454 return tetra_signed_volume(
455 *reinterpret_cast<const vec3*>(p1),
456 *reinterpret_cast<const vec3*>(p2),
457 *reinterpret_cast<const vec3*>(p3),
458 *reinterpret_cast<const vec3*>(p4)
459 );
460 }
461
471 inline double tetra_volume(
472 const vec3& p1, const vec3& p2,
473 const vec3& p3, const vec3& p4
474 ) {
475 return ::fabs(tetra_signed_volume(p1, p2, p3, p4));
476 }
477
489 const vec3& p1, const vec3& p2,
490 const vec3& p3, const vec3& p4
491 );
492
505 inline void triangle_centroid(
506 const vec3& p, const vec3& q, const vec3& r,
507 double a, double b, double c,
508 vec3& Vg, double& V
509 ) {
510 double abc = a + b + c;
511 double area = Geom::triangle_area(p, q, r);
512 V = area / 3.0 * abc;
513 double wp = a + abc;
514 double wq = b + abc;
515 double wr = c + abc;
516 double s = area / 12.0;
517 Vg.x = s * (wp * p.x + wq * q.x + wr * r.x);
518 Vg.y = s * (wp * p.y + wq * q.y + wr * r.y);
519 Vg.z = s * (wp * p.z + wq * q.z + wr * r.z);
520 }
521
534 inline double triangle_mass(
535 const vec3& p, const vec3& q, const vec3& r,
536 double a, double b, double c
537 ) {
538 return Geom::triangle_area(p, q, r) / 3.0 * (
539 sqrt(::fabs(a)) + sqrt(::fabs(b)) + sqrt(::fabs(c))
540 );
541 }
542
554 const vec3& p1,
555 const vec3& p2,
556 const vec3& p3
557 ) {
558 double s = Numeric::random_float64();
559 double t = Numeric::random_float64();
560 if(s + t > 1) {
561 s = 1.0 - s;
562 t = 1.0 - t;
563 }
564 double u = 1.0 - s - t;
565 return vec3(
566 u * p1.x + s * p2.x + t * p3.x,
567 u * p1.y + s * p2.y + t * p3.y,
568 u * p1.z + s * p2.z + t * p3.z
569 );
570 }
571 }
572
578 struct Plane {
579
586 Plane(const vec3& p1, const vec3& p2, const vec3& p3) {
587 vec3 n = cross(p2 - p1, p3 - p1);
588 a = n.x;
589 b = n.y;
590 c = n.z;
591 d = -(a * p1.x + b * p1.y + c * p1.z);
592 }
593
600 Plane(const vec3& p, const vec3& n) {
601 a = n.x;
602 b = n.y;
603 c = n.z;
604 d = -(a * p.x + b * p.y + c * p.z);
605 }
606
611 double a_in, double b_in, double c_in, double d_in
612 ) :
613 a(a_in),
614 b(b_in),
615 c(c_in),
616 d(d_in) {
617 }
618
623 }
624
628 vec3 normal() const {
629 return vec3(a, b, c);
630 }
631
632 double a, b, c, d;
633 };
634
635 /*******************************************************************/
636
640 class Box {
641 public:
642 double xyz_min[3];
643 double xyz_max[3];
644
650 bool contains(const vec3& b) const {
651 for(coord_index_t c = 0; c < 3; ++c) {
652 if(b[c] < xyz_min[c]) {
653 return false;
654 }
655 if(b[c] > xyz_max[c]) {
656 return false;
657 }
658 }
659 return true;
660 }
661 };
662
663 typedef Box Box3d;
664
672 inline bool bboxes_overlap(const Box& B1, const Box& B2) {
673 for(coord_index_t c = 0; c < 3; ++c) {
674 if(B1.xyz_max[c] < B2.xyz_min[c]) {
675 return false;
676 }
677 if(B1.xyz_min[c] > B2.xyz_max[c]) {
678 return false;
679 }
680 }
681 return true;
682 }
683
691 inline void bbox_union(Box& target, const Box& B1, const Box& B2) {
692 for(coord_index_t c = 0; c < 3; ++c) {
693 target.xyz_min[c] = std::min(B1.xyz_min[c], B2.xyz_min[c]);
694 target.xyz_max[c] = std::max(B1.xyz_max[c], B2.xyz_max[c]);
695 }
696 }
697
698 /*******************************************************************/
699
703 class Box2d {
704 public:
705 double xy_min[2];
706 double xy_max[2];
707
713 bool contains(const vec2& b) const {
714 for(coord_index_t c = 0; c < 2; ++c) {
715 if(b[c] < xy_min[c]) {
716 return false;
717 }
718 if(b[c] > xy_max[c]) {
719 return false;
720 }
721 }
722 return true;
723 }
724 };
725
726
734 inline bool bboxes_overlap(const Box2d& B1, const Box2d& B2) {
735 for(coord_index_t c = 0; c < 2; ++c) {
736 if(B1.xy_max[c] < B2.xy_min[c]) {
737 return false;
738 }
739 if(B1.xy_min[c] > B2.xy_max[c]) {
740 return false;
741 }
742 }
743 return true;
744 }
745
753 inline void bbox_union(Box2d& target, const Box2d& B1, const Box2d& B2) {
754 for(coord_index_t c = 0; c < 2; ++c) {
755 target.xy_min[c] = std::min(B1.xy_min[c], B2.xy_min[c]);
756 target.xy_max[c] = std::max(B1.xy_max[c], B2.xy_max[c]);
757 }
758 }
759
760 /*******************************************************************/
761
775 template <class FT> vecng<3,FT> transform_vector(
776 const vecng<3,FT>& v,
777 const Matrix<4,FT>& m
778 ){
779 index_t i,j ;
780 FT result[4] ;
781
782 for(i=0; i<4; i++) {
783 result[i] = 0 ;
784 }
785 for(i=0; i<4; i++) {
786 for(j=0; j<3; j++) {
787 result[i] += v[j] * m(j,i) ;
788 }
789 }
790
791 return vecng<3,FT>(
792 result[0], result[1], result[2]
793 ) ;
794 }
795
811 template <class FT> vecng<3,FT> transform_point(
812 const vecng<3,FT>& v,
813 const Matrix<4,FT>& m
814 ){
815 index_t i,j ;
816 FT result[4] ;
817
818 for(i=0; i<4; i++) {
819 result[i] = 0 ;
820 }
821 for(i=0; i<4; i++) {
822 for(j=0; j<3; j++) {
823 result[i] += v[j] * m(j,i) ;
824 }
825 result[i] += m(3,i);
826 }
827
828 return vecng<3,FT>(
829 result[0] / result[3],
830 result[1] / result[3],
831 result[2] / result[3]
832 ) ;
833 }
834
845 template <class FT> vecng<4,FT> transform_vector(
846 const vecng<4,FT>& v,
847 const Matrix<4,FT>& m
848 ) {
849 index_t i,j ;
850 FT res[4] = {FT(0), FT(0), FT(0), FT(0)};
851
852 for(i=0; i<4; i++) {
853 for(j=0; j<4; j++) {
854 res[i] += v[j] * m(j,i) ;
855 }
856 }
857
858 return vecng<4,FT>(res[0], res[1], res[2], res[3]) ;
859 }
860
861 /******************************************************************/
862
872 mat4 result;
873 result.load_identity();
874 result(3,0) = T.x;
875 result(3,1) = T.y;
876 result(3,2) = T.z;
877 return result;
878 }
879
888 inline mat4 create_scaling_matrix(double s) {
889 mat4 result;
890 result.load_identity();
891 result(0,0) = s;
892 result(1,1) = s;
893 result(2,2) = s;
894 return result;
895 }
896
897 /******************************************************************/
898
902 struct Ray {
908 Ray(vec3 O, vec3 D) : origin(O), direction(D) {
909 }
913 Ray() {
914 }
915 vec3 origin;
916 vec3 direction;
917 };
918
919 /******************************************************************/
920
921}
922
923#endif
924
Common include file, providing basic definitions. Should be included before anything else by all head...
Axis-aligned bounding box.
Definition geometry.h:703
bool contains(const vec2 &b) const
Tests whether a box contains a point.
Definition geometry.h:713
Axis-aligned bounding box.
Definition geometry.h:640
bool contains(const vec3 &b) const
Tests whether a box contains a point.
Definition geometry.h:650
A matrix type.
Definition matrix.h:65
void load_identity()
Sets the matrix to identity.
Definition matrix.h:125
Generic maths vector.
Definition vecg.h:69
T * data()
Gets modifiable vector data.
Definition vecg.h:147
Generic matrix type.
double det(const vec2 &a, const vec2 &b)
Computes the determinant of two vectors.
Definition geometry.h:251
double triangle_signed_area(const vec2 &p1, const vec2 &p2, const vec2 &p3)
Computes the area of a 2d triangle.
Definition geometry.h:351
vec3 tetra_circum_center(const vec3 &p1, const vec3 &p2, const vec3 &p3, const vec3 &p4)
Computes the center of the circumscribed sphere of 3d tetrahedron.
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:553
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:471
vec3 perpendicular(const vec3 &V)
Computes a 3d vector orthogonal to another one.
vec3 triangle_normal(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Computes the normal of a 3d triangle.
Definition geometry.h:274
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:534
double cos_angle(const vec3 &a, const vec3 &b)
Computes the cosine of the angle between two 3d vectors.
Definition geometry.h:208
double triangle_area(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Computes the area of a 3d triangle.
Definition geometry.h:319
double triangle_signed_area_2d(const double *p1, const double *p2, const double *p3)
Computes the area of a 2d triangle.
Definition geometry.h:333
void triangle_centroid(const vec3 &p, const vec3 &q, const vec3 &r, double a, double b, double c, vec3 &Vg, double &V)
Computes the centroid of a 3d triangle with weighted points.
Definition geometry.h:505
bool has_nan(const vec3 &v)
Tests whether a 3d vector has a NaN (not a number) coordinate.
Definition geometry.h:400
double triangle_area_2d(const double *p1, const double *p2, const double *p3)
Computes the area of a 2d triangle.
Definition geometry.h:377
double angle(const vec3 &a, const vec3 &b)
Computes the angle between two 3d vectors.
Definition geometry.h:225
double triangle_area_3d(const double *p1, const double *p2, const double *p3)
Computes the area of a 3d triangle.
Definition geometry.h:297
double tetra_signed_volume(const vec3 &p1, const vec3 &p2, const vec3 &p3, const vec3 &p4)
Computes the signed volume of a 3d tetrahedron.
Definition geometry.h:434
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.
vec3 barycenter(const vec3 &p1, const vec3 &p2)
Computes the barycenter of two points in 3d.
Definition geometry.h:148
vec2 triangle_circumcenter(const vec2 &p1, const vec2 &p2, const vec2 &p3)
Computes the center of the circumscribed circle of a 2d triangle.
float64 random_float64()
Returns a 64 bits float between 0 and 1.
bool is_nan(float32 x)
Checks whether a 32 bits float is "not a number".
Global Vorpaline namespace.
Definition algorithm.h:64
Matrix< 4, Numeric::float64 > mat4
Represents a 4x4 matrix.
Definition geometry.h:133
bool bboxes_overlap(const Box &B1, const Box &B2)
Tests whether two Boxes have a non-empty intersection.
Definition geometry.h:672
vecng< 3, Numeric::float32 > vec3f
Represents points and vectors in 3d with single-precision coordinates.
Definition geometry.h:85
T dot(const vecng< 3, T > &v1, const vecng< 3, T > &v2)
Computes the dot product of 2 vectors. vecng
Definition vecg.h:869
vecng< 4, Numeric::float32 > vec4f
Represents points and vectors in 4d with single-precision coordinates.
Definition geometry.h:92
void bbox_union(Box &target, const Box &B1, const Box &B2)
Computes the smallest Box that encloses two Boxes.
Definition geometry.h:691
vecng< 3, FT > transform_vector(const vecng< 3, FT > &v, const Matrix< 4, FT > &m)
Applies a 3d transform to a 3d vector.
Definition geometry.h:775
vecng< 3, Numeric::float64 > vec3
Represents points and vectors in 3d.
Definition geometry.h:65
vecng< 3, Numeric::int32 > vec3i
Represents points and vectors in 3d with integer coordinates.
Definition geometry.h:107
void geo_clamp(T &x, T min, T max)
Clamps a value to a range.
Definition numeric.h:272
mat4 create_scaling_matrix(double s)
Creates a scaling matrix.
Definition geometry.h:888
vecng< 3, FT > transform_point(const vecng< 3, FT > &v, const Matrix< 4, FT > &m)
Applies a 3d transform to a 3d point.
Definition geometry.h:811
geo_index_t index_t
The type for storing and manipulating indices.
Definition numeric.h:287
vecng< 4, Numeric::float64 > vec4
Represents points and vectors in 4d.
Definition geometry.h:71
mat4 create_translation_matrix(const vec3 &T)
Creates a translation matrix from a vector.
Definition geometry.h:871
vecng< 2, Numeric::int32 > vec2i
Represents points and vectors in 2d with integer coordinates.
Definition geometry.h:100
Matrix< 3, Numeric::float64 > mat3
Represents a 3x3 matrix.
Definition geometry.h:127
vecng< 2, Numeric::float32 > vec2f
Represents points and vectors in 2d with single-precision coordinates.
Definition geometry.h:78
vecng< 2, Numeric::float64 > vec2
Represents points and vectors in 2d.
Definition geometry.h:59
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition numeric.h:321
vecng< 4, Numeric::int32 > vec4i
Represents points and vectors in 4d with integer coordinates.
Definition geometry.h:114
Matrix< 2, Numeric::float64 > mat2
Represents a 2x2 matrix.
Definition geometry.h:121
A 3D Plane.
Definition geometry.h:578
Plane(const vec3 &p1, const vec3 &p2, const vec3 &p3)
Constructs the plane passing through three points.
Definition geometry.h:586
Plane(const vec3 &p, const vec3 &n)
Constructs a plane passign through a point and orthogonal to a vector.
Definition geometry.h:600
Plane(double a_in, double b_in, double c_in, double d_in)
Constructs a plane from the coefficients of its equation.
Definition geometry.h:610
Plane()
Constructs an uninitialized plane.
Definition geometry.h:622
vec3 normal() const
Gets the normal vector of the plane.
Definition geometry.h:628
A Ray, in parametric form.
Definition geometry.h:902
Ray()
Ray constructor.
Definition geometry.h:913
Ray(vec3 O, vec3 D)
Ray constructor.
Definition geometry.h:908