DGtal  1.5.beta
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 
50 namespace 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;
74  static const Dimension dimension = RealPoint::dimension;
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
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:136
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)
int max(int a, int b)
double angle(const DGtal::Z2i::RealPoint &point)