DGtal 1.3.0
Loading...
Searching...
No Matches
NormalCycleFormula.h
1
17#pragma once
18
31#if defined(NormalCycleFormula_RECURSES)
32#error Recursive header files inclusion detected in NormalCycleFormula.h
33#else // defined(NormalCycleFormula_RECURSES)
35#define NormalCycleFormula_RECURSES
36
37#if !defined NormalCycleFormula_h
39#define NormalCycleFormula_h
40
42// Inclusions
43#include <iostream>
44#include <map>
45#include "DGtal/base/Common.h"
46#include "DGtal/math/linalg/SimpleMatrix.h"
47
49
50namespace DGtal
51{
52
54 // class NormalCycleFormula
63 template < typename TRealPoint, typename TRealVector >
65 {
66 typedef TRealPoint RealPoint;
67 typedef TRealVector RealVector;
68 typedef typename RealVector::Component Scalar;
69 typedef std::vector< RealPoint > RealPoints;
70 typedef std::vector< RealVector > RealVectors;
72 typedef std::size_t Size;
73 typedef std::size_t Index;
75
76 //-------------------------------------------------------------------------
77 public:
80
84 static
85 Scalar area( const RealPoints& pts )
86 {
87 if ( pts.size() < 3 ) return 0.0;
88 if ( pts.size() == 3 )
89 return area( pts[ 0 ], pts[ 1 ], pts[ 2 ] );
90 const RealPoint b = barycenter( pts );
91 Scalar a = 0.0;
92 for ( Index i = 0; i < pts.size(); i++ )
93 a += area( b, pts[ i ], pts[ (i+1)%pts.size() ] );
94 return a;
95 }
96
105 static
107 ( const RealPoint& a, const RealPoint& b,
108 const RealVector& right, const RealVector& left )
109 {
110 const RealVector diedre = right.crossProduct( left );
111 const Scalar n = std::min( 1.0, std::max( diedre.norm(), 0.0 ) );
112 const Scalar angle = ( diedre.dot( b - a) < 0.0 )
113 ? asin( n ) : - asin( n );
114 return ( b - a ).norm() * angle;
115 }
116
123 static
125 ( const RealPoint& a, const RealPoints& vtcs )
126 {
127 Scalar angle_sum = 0.0;
128 for ( Size i = 0; i < vtcs.size(); i++ )
129 angle_sum += acos( (vtcs[i] - a).getNormalized()
130 .dot( ( vtcs[(i+1)%vtcs.size()] - a ).getNormalized() ) );
131 return 2.0 * M_PI - angle_sum;
132 }
133
141 static
143 ( const RealPoint& a, const RealPoints& pairs )
144 {
145 Scalar angle_sum = 0.0;
146 for ( Size i = 0; i < pairs.size(); i += 2 )
147 angle_sum += acos( ( pairs[i] - a ).getNormalized()
148 .dot( ( pairs[i+1] - a ).getNormalized() ) );
149 return 2.0 * M_PI - angle_sum;
150 }
151
158 static
160 ( const RealPoint& a, const RealPoint& b,
161 const RealVector& right, const RealVector& left )
162 {
163 const RealVector diedre = right.crossProduct( left );
164 const Scalar length = std::max( 0.0, std::min( 1.0, diedre.norm() ) );
165 const Scalar angle = ( diedre.dot( b - a) > 0.0 )
166 ? asin( length ) : - asin( length );
167 RealVector e_p = right + left;
168 RealVector e_m = right - left;
169 const Scalar norm_e_p = e_p.norm();
170 const Scalar norm_e_m = e_m.norm();
171 e_p = norm_e_p > 1e-10 ? e_p / norm_e_p : RealVector::zero;
172 e_m = norm_e_m > 1e-10 ? e_m / norm_e_m : RealVector::zero;
173 const RealTensor T_p =
174 { e_p[ 0 ] * e_p[ 0 ], e_p[ 0 ] * e_p[ 1 ], e_p[ 0 ] * e_p[ 2 ],
175 e_p[ 1 ] * e_p[ 0 ], e_p[ 1 ] * e_p[ 1 ], e_p[ 1 ] * e_p[ 2 ],
176 e_p[ 2 ] * e_p[ 0 ], e_p[ 2 ] * e_p[ 1 ], e_p[ 2 ] * e_p[ 2 ] };
177 const RealTensor T_m =
178 { e_m[ 0 ] * e_m[ 0 ], e_m[ 0 ] * e_m[ 1 ], e_m[ 0 ] * e_m[ 2 ],
179 e_m[ 1 ] * e_m[ 0 ], e_m[ 1 ] * e_m[ 1 ], e_m[ 1 ] * e_m[ 2 ],
180 e_m[ 2 ] * e_m[ 0 ], e_m[ 2 ] * e_m[ 1 ], e_m[ 2 ] * e_m[ 2 ] };
181 return 0.5 * ( b - a ).norm()
182 * ( ( angle - sin( angle ) ) * T_p + ( angle + sin( angle ) ) * T_m );
183 }
184
191 static
193 ( const RealPoint& a, const RealPoint& b,
194 const RealVector& right, const RealVector& left )
195 {
196 const RealVector diedre = right.crossProduct( left );
197 const Scalar length = std::max( 0.0, std::min( 1.0, diedre.norm() ) );
198 const Scalar angle = ( diedre.dot( b - a) > 0.0 )
199 ? asin( length ) : - asin( length );
200 const Scalar norm_ab = (b - a).norm();
201 const RealVector e = norm_ab > 1e-10 ? (b - a) / norm_ab : RealVector::zero;
202 const RealTensor T =
203 { e[ 0 ] * e[ 0 ], e[ 0 ] * e[ 1 ], e[ 0 ] * e[ 2 ],
204 e[ 1 ] * e[ 0 ], e[ 1 ] * e[ 1 ], e[ 1 ] * e[ 2 ],
205 e[ 2 ] * e[ 0 ], e[ 2 ] * e[ 1 ], e[ 2 ] * e[ 2 ] };
206 return ( 0.5 * norm_ab * angle ) * T; // JOL * 0.5
207 }
208
209
210 //-------------------------------------------------------------------------
211 public:
214
218 static
220 {
221 RealPoint b;
222 for ( auto p : pts ) b += p;
223 b /= pts.size();
224 return b;
225 }
226
233 static
234 RealVector normal( const RealPoint& a, const RealPoint& b, const RealPoint& c )
235 {
236 return ( ( b - a ).crossProduct( c - a ) ).getNormalized();
237 }
238
244 static
245 Scalar area( const RealPoint& a, const RealPoint& b, const RealPoint& c )
246 {
247 return 0.5 * ( ( b - a ).crossProduct( c - a ) ).norm();
248 }
249
253 static
255 {
256 RealVector avg;
257 for ( auto v : vecs ) avg += v;
258 auto avg_norm = avg.norm();
259 return avg_norm != 0.0 ? avg / avg_norm : avg;
260 }
261
262
264
265 };
266
267} // namespace DGtal
268
269
271// Includes inline functions.
272//#include "NormalCycleFormula.ih"
273
274// //
276
277#endif // !defined NormalCycleFormula_h
278
279#undef NormalCycleFormula_RECURSES
280#endif // else defined(NormalCycleFormula_RECURSES)
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:614
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
static const Dimension dimension
Copy of the static dimension of the Point/Vector.
Definition: PointVector.h:626
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:76
DGtal is the top-level namespace which contains all DGtal functions and types.
auto crossProduct(PointVector< 3, LeftEuclideanRing, LeftContainer > const &lhs, PointVector< 3, RightEuclideanRing, RightContainer > const &rhs) -> decltype(DGtal::constructFromArithmeticConversion(lhs, rhs))
Cross product of two 3D Points/Vectors.
DGtal::uint32_t Dimension
Definition: Common.h:137
Aim: A helper class that provides static methods to compute normal cycle formulas of curvatures.
static Scalar area(const RealPoints &pts)
SimpleMatrix< Scalar, 3, 3 > RealTensor
std::vector< RealPoint > RealPoints
static RealTensor anisotropicCurvatureH2(const RealPoint &a, const RealPoint &b, const RealVector &right, const RealVector &left)
static Scalar area(const RealPoint &a, const RealPoint &b, const RealPoint &c)
static Scalar gaussianCurvature(const RealPoint &a, const RealPoints &vtcs)
static Scalar gaussianCurvatureWithPairs(const RealPoint &a, const RealPoints &pairs)
static Scalar twiceMeanCurvature(const RealPoint &a, const RealPoint &b, const RealVector &right, const RealVector &left)
static RealVector normal(const RealPoint &a, const RealPoint &b, const RealPoint &c)
RealVector::Component Scalar
static const Dimension dimension
static RealPoint barycenter(const RealPoints &pts)
std::vector< RealVector > RealVectors
static RealTensor anisotropicCurvatureH1(const RealPoint &a, const RealPoint &b, const RealVector &right, const RealVector &left)
static RealVector averageUnitVector(const RealVectors &vecs)