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