1 // Copyright (c) 2007-2009 INRIA Sophia-Antipolis (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Mesh_3/include/CGAL/Mesh_3/min_dihedral_angle.h $
7 // $Id: min_dihedral_angle.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Laurent RINEAU, Stephane Tayeb
12
13 #ifndef CGAL_MESH_3_MIN_DIHEDRAL_ANGLE_H
14 #define CGAL_MESH_3_MIN_DIHEDRAL_ANGLE_H
15
16 #include <CGAL/license/Mesh_3.h>
17
18
19 #include <cmath>
20 #include <CGAL/Kernel_traits.h>
21 #include <CGAL/utils.h>
22 #include <CGAL/number_utils.h>
23
24 namespace CGAL {
25
26 namespace Mesh_3 {
27
28 #ifdef CGAL_MESH_3_OLD_MINIMUM_DIHEDRAL_ANGLE
29
30 template <typename K>
31 typename K::FT
32 minimum_dihedral_angle(
33 const typename K::Point_3& p0,
34 const typename K::Point_3& p1,
35 const typename K::Point_3& p2,
36 const typename K::Point_3& p3,
37 K k = K())
38 {
39 typedef typename K::FT FT;
40
41 typename K::Compute_squared_distance_3 sq_distance =
42 k.compute_squared_distance_3_object();
43 typename K::Compute_volume_3 volume =
44 k.compute_volume_3_object();
45 typename K::Compute_area_3 area =
46 k.compute_area_3_object();
47
48 FT a_012 = area(p0,p1,p2);
49 FT a_013 = area(p0,p1,p3);
50 FT a_123 = area(p1,p2,p3);
51 FT a_023 = area(p0,p2,p3);
52
53 FT min_quotient =
54 CGAL::sqrt(sq_distance(p0, p1)) / a_012 / a_013;
55
56 min_quotient =
57 (CGAL::min)(min_quotient,
58 CGAL::sqrt(sq_distance(p0, p2)) / a_012 / a_023);
59 min_quotient =
60 (CGAL::min)(min_quotient,
61 CGAL::sqrt(sq_distance(p0, p3)) / a_013 / a_023);
62 min_quotient =
63 (CGAL::min)(min_quotient,
64 CGAL::sqrt(sq_distance(p1, p2)) / a_012 / a_123);
65 min_quotient =
66 (CGAL::min)(min_quotient,
67 CGAL::sqrt(sq_distance(p1, p3)) / a_013 / a_123);
68 min_quotient =
69 (CGAL::min)(min_quotient,
70 CGAL::sqrt(sq_distance(p2, p3)) / a_023 / a_123);
71
72 const FT result (std::asin( FT(1.5) * volume(p0, p1, p2, p3) * min_quotient )
73 * FT(180) / FT(CGAL_PI));
74
75 return CGAL::abs(result);
76 }
77
78 #else // not CGAL_MESH_3_OLD_MINIMUM_DIHEDRAL_ANGLE
79
80 template <typename K>
81 typename K::FT
82 minimum_dihedral_angle(
83 const typename K::Point_3& p0,
84 const typename K::Point_3& p1,
85 const typename K::Point_3& p2,
86 const typename K::Point_3& p3,
87 K k = K())
88 {
89 typedef typename K::FT FT;
90
91 typename K::Compute_determinant_3 determinant =
92 k.compute_determinant_3_object();
93 typename K::Construct_cross_product_vector_3 cp =
94 k.construct_cross_product_vector_3_object();
95 typename K::Compute_scalar_product_3 sp =
96 k.compute_scalar_product_3_object();
97 typename K::Construct_vector_3 cv =
98 k.construct_vector_3_object();
99
100 typename K::Vector_3 v01 = cv(p0,p1);
101 typename K::Vector_3 v02 = cv(p0,p2);
102 typename K::Vector_3 v03 = cv(p0,p3);
103 typename K::Vector_3 v12 = cv(p1,p2);
104 typename K::Vector_3 v13 = cv(p1,p3);
105 typename K::Vector_3 v23 = cv(p2,p3);
106
107 typename K::Vector_3 v_01_02 = cp(v01,v02);
108 FT a_012 = v_01_02*v_01_02;
109
110
111 typename K::Vector_3 v_01_03 = cp(v01,v03);
112 FT a_013 = v_01_03*v_01_03;
113
114
115 typename K::Vector_3 v_12_13 = cp(v12,v13);
116 FT a_123 = v_12_13*v_12_13;
117
118 typename K::Vector_3 v_02_03 = cp(v02,v03);
119 FT a_023 = v_02_03*v_02_03;
120
121 FT min_quotient = sp(v01,v01) / (a_012 * a_013);
122 min_quotient = (CGAL::min)(min_quotient,
123 sp(v02,v02) / (a_012 * a_023));
124 min_quotient = (CGAL::min)(min_quotient,
125 (v03*v03) / (a_013 * a_023));
126 min_quotient = (CGAL::min)(min_quotient,
127 sp(v12,v12) / (a_012 * a_123));
128 min_quotient = (CGAL::min)(min_quotient,
129 sp(v13,v13) / (a_013 * a_123));
130 min_quotient = (CGAL::min)(min_quotient,
131 sp(v23,v23) / (a_023 * a_123));
132 min_quotient = sqrt(min_quotient);
133
134 return CGAL::abs(std::asin(determinant(v01, v02, v03) * min_quotient)
135 * FT(180) / FT(CGAL_PI));
136 }
137
138 #endif // CGAL_MESH_3_OLD_MINIMUM_DIHEDRAL_ANGLE
139
140 template <typename K>
141 typename K::FT
142 minimum_dihedral_angle(const typename K::Tetrahedron_3& t, K k = K() )
143 {
144 return minimum_dihedral_angle(t[0],t[1],t[2],t[3],k);
145 }
146
147 template <typename Tetrahedron_3>
148 typename Kernel_traits<Tetrahedron_3>::Kernel::FT
minimum_dihedral_angle(const Tetrahedron_3 & t)149 minimum_dihedral_angle(const Tetrahedron_3& t )
150 {
151 return minimum_dihedral_angle(t, typename Kernel_traits<Tetrahedron_3>::Kernel() );
152 }
153
154 template<typename C3T3>
155 typename C3T3::Triangulation::Geom_traits::FT
minimum_dihedral_angle_in_c3t3(const C3T3 & c3t3)156 minimum_dihedral_angle_in_c3t3(const C3T3& c3t3)
157 {
158 typedef typename C3T3::Triangulation::Geom_traits K;
159 typename K::FT min_angle = (typename K::FT)(90.);
160
161 typename C3T3::Cells_in_complex_iterator cit;
162 for(cit = c3t3.cells_in_complex_begin();
163 cit != c3t3.cells_in_complex_end();
164 ++cit)
165 {
166 min_angle = (std::min)(min_angle,
167 minimum_dihedral_angle(c3t3.triangulation().tetrahedron(cit)));
168 }
169 return min_angle;
170 }
171
172
173 } // end namespace Mesh_3
174 } // end namespace CGAL
175
176 #endif // CGAL_MESH_3_MIN_DIHEDRAL_ANGLE_H
177