1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkTriangleHelper_hxx
19 #define itkTriangleHelper_hxx
20
21 #include "itkTriangleHelper.h"
22 #include "itkMath.h"
23
24 namespace itk
25 {
26 template< typename TPoint >
IsObtuse(const PointType & iA,const PointType & iB,const PointType & iC)27 bool TriangleHelper< TPoint >::IsObtuse(const PointType & iA, const PointType & iB, const PointType & iC)
28 {
29 VectorType v01 = iB - iA;
30 VectorType v02 = iC - iA;
31 VectorType v12 = iC - iB;
32
33 if ( v01 * v02 < 0.0 )
34 {
35 return true;
36 }
37 else
38 {
39 if ( v02 * v12 < 0.0 )
40 {
41 return true;
42 }
43 else
44 {
45 if ( v01 * -v12 < 0.0 )
46 {
47 return true;
48 }
49 else
50 {
51 return false;
52 }
53 }
54 }
55 }
56
57 template< typename TPoint >
58 typename TriangleHelper< TPoint >::VectorType
ComputeNormal(const PointType & iA,const PointType & iB,const PointType & iC)59 TriangleHelper< TPoint >::ComputeNormal(const PointType & iA,
60 const PointType & iB,
61 const PointType & iC)
62 {
63 CrossVectorType cross;
64 VectorType w = cross (iB - iA, iC - iA);
65 CoordRepType l2 = w.GetSquaredNorm();
66
67 if ( l2 != 0.0 )
68 {
69 w /= std::sqrt(l2);
70 }
71
72 return w;
73 }
74
75 template< typename TPoint >
76 typename TriangleHelper< TPoint >::CoordRepType
Cotangent(const PointType & iA,const PointType & iB,const PointType & iC)77 TriangleHelper< TPoint >::Cotangent(const PointType & iA,
78 const PointType & iB,
79 const PointType & iC)
80 {
81 VectorType v21 = iA - iB;
82
83 CoordRepType v21_l2 = v21.GetSquaredNorm();
84
85 if ( Math::NotAlmostEquals( v21_l2, NumericTraits< CoordRepType >::ZeroValue() ) )
86 {
87 v21 /= std::sqrt(v21_l2);
88 }
89
90 VectorType v23 = iC - iB;
91 CoordRepType v23_l2 = v23.GetSquaredNorm();
92 if ( Math::NotAlmostEquals( v23_l2, NumericTraits< CoordRepType >::ZeroValue() ) )
93 {
94 v23 /= std::sqrt(v23_l2);
95 }
96
97 CoordRepType bound(0.999999);
98
99 CoordRepType cos_theta = std::max( -bound,
100 std::min(bound, v21 * v23) );
101
102 return 1.0 / std::tan( std::acos(cos_theta) );
103 }
104
105 template< typename TPoint >
106 typename TriangleHelper< TPoint >::PointType
ComputeBarycenter(const CoordRepType & iA1,const PointType & iP1,const CoordRepType & iA2,const PointType & iP2,const CoordRepType & iA3,const PointType & iP3)107 TriangleHelper< TPoint >::ComputeBarycenter(
108 const CoordRepType & iA1, const PointType & iP1,
109 const CoordRepType & iA2, const PointType & iP2,
110 const CoordRepType & iA3, const PointType & iP3)
111 {
112 PointType oPt;
113
114 CoordRepType total = iA1 + iA2 + iA3;
115
116 if ( Math::AlmostEquals( total, NumericTraits< CoordRepType >::ZeroValue() ) )
117 {
118 //in such case there is no barycenter;
119 oPt.Fill(0.);
120 return oPt;
121 }
122
123 CoordRepType inv_total = 1. / total;
124 CoordRepType a1 = iA1 * inv_total;
125 CoordRepType a2 = iA2 * inv_total;
126 CoordRepType a3 = iA3 * inv_total;
127
128 for ( unsigned int dim = 0; dim < PointDimension; ++dim )
129 {
130 oPt[dim] = a1 * iP1[dim] + a2 * iP2[dim] + a3 * iP3[dim];
131 }
132
133 return oPt;
134 }
135
136 template< typename TPoint >
137 typename TriangleHelper< TPoint >::CoordRepType
ComputeAngle(const PointType & iP1,const PointType & iP2,const PointType & iP3)138 TriangleHelper< TPoint >::ComputeAngle(const PointType & iP1,
139 const PointType & iP2,
140 const PointType & iP3)
141 {
142 VectorType v21 = iP1 - iP2;
143 VectorType v23 = iP3 - iP2;
144
145 CoordRepType v21_l2 = v21.GetSquaredNorm();
146 CoordRepType v23_l2 = v23.GetSquaredNorm();
147
148 if ( v21_l2 != 0.0 )
149 {
150 v21 /= std::sqrt(v21_l2);
151 }
152 if ( v23_l2 != 0.0 )
153 {
154 v23 /= std::sqrt(v23_l2);
155 }
156
157 CoordRepType bound(0.999999);
158
159 CoordRepType cos_theta = std::max( -bound,
160 std::min(bound, v21 * v23) );
161
162 return std::acos(cos_theta);
163 }
164
165 template< typename TPoint >
166 typename TriangleHelper< TPoint >::PointType
ComputeGravityCenter(const PointType & iP1,const PointType & iP2,const PointType & iP3)167 TriangleHelper< TPoint >::ComputeGravityCenter(
168 const PointType & iP1,
169 const PointType & iP2,
170 const PointType & iP3)
171 {
172 return ComputeBarycenter(1., iP1, 1., iP2, 1., iP3);
173 }
174
175 template< typename TPoint >
176 typename TriangleHelper< TPoint >::PointType
ComputeCircumCenter(const PointType & iP1,const PointType & iP2,const PointType & iP3)177 TriangleHelper< TPoint >::ComputeCircumCenter(
178 const PointType & iP1,
179 const PointType & iP2,
180 const PointType & iP3)
181 {
182 PointType oPt;
183
184 oPt.Fill (0.0);
185
186 CoordRepType a = iP2.SquaredEuclideanDistanceTo (iP3);
187 CoordRepType b = iP1.SquaredEuclideanDistanceTo (iP3);
188 CoordRepType c = iP2.SquaredEuclideanDistanceTo (iP1);
189
190 CoordRepType Weight[3];
191 Weight[0] = a * ( b + c - a );
192 Weight[1] = b * ( c + a - b );
193 Weight[2] = c * ( a + b - c );
194
195 return ComputeBarycenter(Weight[0], iP1, Weight[1], iP2, Weight[2], iP3);
196 }
197
198 template< typename TPoint >
199 typename TriangleHelper< TPoint >::PointType
ComputeConstrainedCircumCenter(const PointType & iP1,const PointType & iP2,const PointType & iP3)200 TriangleHelper< TPoint >::ComputeConstrainedCircumCenter(const PointType & iP1,
201 const PointType & iP2, const PointType & iP3)
202 {
203 const CoordRepType a = iP2.SquaredEuclideanDistanceTo (iP3);
204 const CoordRepType b = iP1.SquaredEuclideanDistanceTo (iP3);
205 const CoordRepType c = iP2.SquaredEuclideanDistanceTo (iP1);
206
207 CoordRepType Weight[3] = {
208 a * ( b + c - a ),
209 b * ( c + a - b ),
210 c * ( a + b - c )
211 };
212
213 for (auto & i : Weight)
214 {
215 if ( i < 0.0 )
216 {
217 i = 0.;
218 }
219 }
220
221 return ComputeBarycenter(Weight[0], iP1, Weight[1], iP2, Weight[2], iP3);
222 }
223
224 template< typename TPoint >
225 typename TriangleHelper< TPoint >::CoordRepType
ComputeArea(const PointType & iP1,const PointType & iP2,const PointType & iP3)226 TriangleHelper< TPoint >::ComputeArea(const PointType & iP1,
227 const PointType & iP2,
228 const PointType & iP3)
229 {
230 CoordRepType a = iP2.EuclideanDistanceTo (iP3);
231 CoordRepType b = iP1.EuclideanDistanceTo (iP3);
232 CoordRepType c = iP2.EuclideanDistanceTo (iP1);
233
234 CoordRepType s = 0.5 * ( a + b + c );
235
236 return static_cast< CoordRepType >( std::sqrt ( s * ( s - a ) * ( s - b ) * ( s - c ) ) );
237 }
238
239 template< typename TPoint >
240 typename TriangleHelper< TPoint >::CoordRepType
ComputeMixedArea(const PointType & iP1,const PointType & iP2,const PointType & iP3)241 TriangleHelper< TPoint >::ComputeMixedArea(const PointType & iP1,
242 const PointType & iP2,
243 const PointType & iP3)
244 {
245 using TriangleType = TriangleHelper< TPoint >;
246
247 if ( !TriangleType::IsObtuse(iP1, iP2, iP3) )
248 {
249 auto sq_d01 = static_cast< CoordRepType >( iP1.SquaredEuclideanDistanceTo(iP2) );
250 auto sq_d02 = static_cast< CoordRepType >( iP1.SquaredEuclideanDistanceTo(iP3) );
251
252 CoordRepType cot_theta_210 = TriangleType::Cotangent(iP3, iP2, iP1);
253 CoordRepType cot_theta_021 = TriangleType::Cotangent(iP1, iP3, iP2);
254
255 return 0.125 * ( sq_d02 * cot_theta_210 + sq_d01 * cot_theta_021 );
256 }
257 else
258 {
259 auto area = static_cast< CoordRepType >( TriangleType::ComputeArea(iP1, iP2, iP3) );
260
261 if ( ( iP2 - iP1 ) * ( iP3 - iP1 ) < NumericTraits< CoordRepType >::ZeroValue() )
262 {
263 return 0.5 * area;
264 }
265 else
266 {
267 return 0.25 * area;
268 }
269 }
270 }
271 }
272
273 #endif
274