DGtal  1.5.beta
CorrectedNormalCurrentFormula.h
1 
17 #pragma once
18 
31 #if defined(CorrectedNormalCurrentFormula_RECURSES)
32 #error Recursive header files inclusion detected in CorrectedNormalCurrentFormula.h
33 #else // defined(CorrectedNormalCurrentFormula_RECURSES)
35 #define CorrectedNormalCurrentFormula_RECURSES
36 
37 #if !defined CorrectedNormalCurrentFormula_h
39 #define CorrectedNormalCurrentFormula_h
40 
42 // Inclusions
43 #include <iostream>
44 #include <map>
45 #include "DGtal/base/Common.h"
46 #include "DGtal/math/linalg/SimpleMatrix.h"
47 #include "DGtal/geometry/tools/SphericalTriangle.h"
48 
50 
51 namespace DGtal
52 {
53 
55  // class CorrectedNormalCurrentFormula
98  template < typename TRealPoint, typename TRealVector >
99  struct CorrectedNormalCurrentFormula
100  {
101  typedef TRealPoint RealPoint;
102  typedef TRealVector RealVector;
103  typedef typename RealVector::Component Scalar;
104  typedef std::vector< RealPoint > RealPoints;
105  typedef std::vector< RealVector > RealVectors;
106  typedef SimpleMatrix< Scalar, 3, 3 > RealTensor;
107  typedef std::size_t Index;
108  static const Dimension dimension = RealPoint::dimension;
109 
110  //-------------------------------------------------------------------------
111  public:
114 
122  static
123  Scalar mu0ConstantU
124  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
125  const RealVector& u )
126  {
127  return 0.5 * ( b - a ).crossProduct( c - a ).dot( u );
128  }
129 
142  static
143  Scalar mu0InterpolatedU
144  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
145  const RealVector& ua, const RealVector& ub, const RealVector& uc,
146  bool unit_u = false )
147  {
148  // MU0=1/2*det( uM, B-A, C-A )
149  // = 1/2 < ( (u_A + u_B + u_C)/3.0 ) | (AB x AC ) >
150  RealVector uM = ( ua+ub+uc ) / 3.0;
151  if ( unit_u )
152  {
153  auto uM_norm = uM.norm();
154  uM = uM_norm == 0.0 ? uM : uM / uM_norm;
155  }
156  return 0.5 * ( b - a ).crossProduct( c - a ).dot( uM );
157  }
158 
164  static
165  Scalar mu0ConstantU( const RealPoints& pts, const RealVector& u )
166  {
167  if ( pts.size() < 3 ) return 0.0;
168  if ( pts.size() == 3 )
169  return mu0ConstantU( pts[ 0 ], pts[ 1 ], pts[ 2 ], u );
170  const RealPoint b = barycenter( pts );
171  Scalar a = 0.0;
172  for ( Index i = 0; i < pts.size(); i++ )
173  a += mu0ConstantU( b, pts[ i ], pts[ (i+1)%pts.size() ], u );
174  return a;
175  }
176 
185  static
186  Scalar mu0InterpolatedU( const RealPoints& pts, const RealVectors& u,
187  bool unit_u = false )
188  {
189  ASSERT( pts.size() == u.size() );
190  if ( pts.size() < 3 ) return 0.0;
191  if ( pts.size() == 3 )
192  return mu0InterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
193  u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
194  const RealPoint b = barycenter( pts );
195  const RealVector ub = averageUnitVector( u );
196  Scalar a = 0.0;
197  for ( Index i = 0; i < pts.size(); i++ )
198  a += mu0InterpolatedU( b, pts[ i ], pts[ (i+1)%pts.size() ],
199  ub, u[ i ], u[ (i+1)%pts.size() ], unit_u );
200  return a;
201  }
202 
204 
205  //-------------------------------------------------------------------------
206  public:
209 
225  static
226  Scalar mu1ConstantUAtEdge
227  ( const RealPoint& a, const RealPoint& b,
228  const RealVector& ur, const RealVector& ul )
229  {
230  RealVector e = b - a;
231  RealVector e1 = ur.crossProduct( ul );
232  const Scalar l1 = std::min( e1.norm(), 1.0 );
233  if ( l1 < 1e-10 ) return 0.0;
234  e1 /= l1;
235  const Scalar Psi = asin( l1 );
236  return -Psi * e.dot( e1 );
237  }
238 
247  static
248  Scalar mu1ConstantU
249  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
250  const RealVector& u )
251  {
252  (void)a; (void)b; (void)c; (void)u;
253  return 0.0;
254  }
255 
270  static
271  Scalar mu1InterpolatedU
272  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
273  const RealVector& ua, const RealVector& ub, const RealVector& uc,
274  bool unit_u = false )
275  {
276  // MU1=1/2( | uM u_C-u_B A | + | uM u_A-u_C B | + | uM u_B-u_A C |
277  RealVector uM = ( ua+ub+uc ) / 3.0;
278  if ( unit_u ) uM /= uM.norm();
279  return 0.5 * ( uM.crossProduct( uc - ub ).dot( a )
280  + uM.crossProduct( ua - uc ).dot( b )
281  + uM.crossProduct( ub - ua ).dot( c ) );
282  }
283 
290  static
291  Scalar mu1ConstantU( const RealPoints& pts, const RealVector& u )
292  {
293  (void)pts; //not used parameter
294  (void)u;
295  return 0.0;
296  }
297 
306  static
307  Scalar mu1InterpolatedU( const RealPoints& pts, const RealVectors& u,
308  bool unit_u = false )
309  {
310  ASSERT( pts.size() == u.size() );
311  if ( pts.size() < 3 ) return 0.0;
312  if ( pts.size() == 3 )
313  return mu1InterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
314  u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
315  const RealPoint b = barycenter( pts );
316  const RealVector ub = averageUnitVector( u );
317  Scalar a = 0.0;
318  for ( Index i = 0; i < pts.size(); i++ )
319  a += mu1InterpolatedU( b, pts[ i ], pts[ (i+1)%pts.size() ],
320  ub, u[ i ], u[ (i+1)%pts.size() ], unit_u );
321  return a;
322  }
323 
325 
326  //-------------------------------------------------------------------------
327  public:
330 
338  static
339  Scalar mu2ConstantU
340  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
341  const RealVector& u )
342  {
343  (void)a; (void)b; (void)c; (void)u;
344  return 0.0;
345  }
346 
359  static
360  Scalar mu2ConstantUAtVertex
361  ( const RealPoint& a, const RealVectors& vu )
362  {
363  typedef SpaceND< dimension > Space;
364  typedef SphericalTriangle<Space> ST;
365  RealVector avg_u;
366  for ( const auto& u : vu ) avg_u += u;
367  const Scalar l = avg_u.norm();
368  if ( l < 1e-10 )
369  {
370  trace.warning() << "[CorrectedNormalCurrentFormula::mu2ConstantUAtVertex]"
371  << " Invalid surrounding normal vectors at vertex "
372  << a << "."
373  << " Unit normal vectors should lie strictly in the same"
374  << " hemisphere." << std::endl;
375  return 0.0;
376  }
377  avg_u /= l;
378  Scalar mu2 = 0.0;
379  const auto n = vu.size();
380  for ( size_t i = 0; i < n; ++i )
381  {
382  ST st( avg_u, vu[ i ], vu[ (i+1) % n ] );
383  mu2 += st.algebraicArea();
384  }
385  return mu2;
386  }
387 
402  static
403  Scalar mu2InterpolatedU
404  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
405  const RealVector& ua, const RealVector& ub, const RealVector& uc,
406  bool unit_u = false )
407  {
408  (void)a; //not used
409  (void)b;
410  (void)c;
411 
412  // Using non unitary interpolated normals give
413  // MU2=1/2*det( uA, uB, uC )
414  // When normals are unitary, it is the area of a spherical triangle.
415  if ( unit_u )
416  {
417  typedef SpaceND< dimension > Space;
418  SphericalTriangle<Space> ST( ua, ub, uc ); // check order.
419  return ST.algebraicArea();
420  }
421  else
422  return 0.5 * ( ua.crossProduct( ub ).dot( uc ) );
423  }
424 
430  static
431  Scalar mu2ConstantU( const RealPoints& pts, const RealVector& u )
432  {
433  (void) pts; (void) u;
434  return 0.0;
435  }
436 
445  static
446  Scalar mu2InterpolatedU( const RealPoints& pts, const RealVectors& u,
447  bool unit_u = false )
448  {
449  ASSERT( pts.size() == u.size() );
450  if ( pts.size() < 3 ) return 0.0;
451  if ( pts.size() == 3 )
452  return mu2InterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
453  u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
454  const RealPoint b = barycenter( pts );
455  const RealVector ub = averageUnitVector( u );
456  Scalar a = 0.0;
457  for ( Index i = 0; i < pts.size(); i++ )
458  a += mu2InterpolatedU( b, pts[ i ], pts[ (i+1)%pts.size() ],
459  ub, u[ i ], u[ (i+1)%pts.size() ], unit_u );
460  return a;
461  }
462 
464 
465  //-------------------------------------------------------------------------
466  public:
469 
477  static
478  RealTensor muXYConstantU
479  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
480  const RealVector& u )
481  {
482  (void)a; (void)b; (void)c; (void)u;
483  return RealTensor { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
484  }
485 
501  static
502  RealTensor muXYConstantUAtEdge
503  ( const RealPoint& a, const RealPoint& b,
504  const RealVector& ur, const RealVector& ul )
505  {
506  RealTensor M { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
507  RealVector e = b - a;
508  RealVector e1 = ur.crossProduct( ul );
509  const Scalar l1 = std::min( e1.norm(), 1.0 );
510  if ( l1 < 1e-10 ) return M;
511  e1 /= l1;
512  const Scalar psi = asin( l1 );
513  const Scalar sin_psi = sin( psi );
514  const Scalar sin_2psi = sin( 2.0 * psi );
515  const RealVector e_x_ur = e. crossProduct( ur );
516  const RealVector e1_x_ur = e1.crossProduct( ur );
517  const RealVector e_x_e1_x_ur = e. crossProduct( e1_x_ur );
518  for ( Dimension i = 0; i < 3; i++ )
519  for ( Dimension j = 0; j < 3; j++ )
520  {
521  Scalar v = (-psi - 0.5*sin_2psi) * e_x_ur[ i ] * e1_x_ur[ j ]
522  + (sin_psi*sin_psi) * ( e_x_ur[ i ] * ur[ j ]
523  - e_x_e1_x_ur[ i ] * e1_x_ur[ j ] )
524  + (psi - 0.5*sin_2psi) * e_x_e1_x_ur[ i ] * ur[ j ];
525  M.setComponent( i, j, -0.5 * v );
526  }
527  return M;
528  }
529 
542  static
543  RealTensor muXYInterpolatedU
544  ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
545  const RealVector& ua, const RealVector& ub, const RealVector& uc,
546  bool unit_u = false )
547  {
548  RealTensor T;
549  // MUXY = 1/2 < uM | < Y | uc-ua > X x (b-a) - < Y | ub-ua > X x (c-a) >
550  // MUXY = 1/2 ( < Y | ub-ua > | X uM (c-a) | - < Y | uc-ua > | X uM (b-a) | )
551  RealVector uM = ( ua+ub+uc ) / 3.0;
552  if ( unit_u ) uM /= uM.norm();
553  const RealVector uac = uc - ua;
554  const RealVector uab = ub - ua;
555  const RealVector ab = b - a;
556  const RealVector ac = c - a;
557  for ( Dimension i = 0; i < 3; ++i ) {
558  RealVector X = RealVector::base( i, 1.0 );
559  for ( Dimension j = 0; j < 3; ++j ) {
560  // Since RealVector Y = RealVector::base( j, 1.0 );
561  // < Y | uac > = uac[ j ]
562  const Scalar tij =
563  0.5 * uM.dot( uac[ j ] * X.crossProduct( ab )
564  - uab[ j ] * X.crossProduct( ac ) );
565  T.setComponent( i, j, tij );
566  }
567  }
568  return T;
569  }
570 
576  static
577  RealTensor muXYConstantU( const RealPoints& pts, const RealVector& u )
578  {
579  (void)pts; (void)u;
580  return RealTensor { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
581  }
582 
591  static
592  RealTensor muXYInterpolatedU( const RealPoints& pts, const RealVectors& u,
593  bool unit_u = false )
594  {
595  ASSERT( pts.size() == u.size() );
596  if ( pts.size() < 3 ) return RealTensor { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
597  if ( pts.size() == 3 )
598  return muXYInterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
599  u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
600  const RealPoint b = barycenter( pts );
601  const RealVector ub = averageUnitVector( u );
602  RealTensor T = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
603  for ( Index i = 0; i < pts.size(); i++ )
604  T += muXYInterpolatedU( b, pts[ i ], pts[ (i+1)%pts.size() ],
605  ub, u[ i ], u[ (i+1)%pts.size() ], unit_u );
606  return T;
607  }
608 
610 
611 
612  //-------------------------------------------------------------------------
613  public:
616 
620  static
621  RealPoint barycenter( const RealPoints& pts )
622  {
623  RealPoint b;
624  for ( auto p : pts ) b += p;
625  b /= pts.size();
626  return b;
627  }
628 
632  static
633  RealVector averageUnitVector( const RealVectors& vecs )
634  {
635  RealVector avg;
636  for ( auto v : vecs ) avg += v;
637  auto avg_norm = avg.norm();
638  return avg_norm != 0.0 ? avg / avg_norm : avg;
639  }
640 
642 
643  };
644 
645 } // namespace DGtal
646 
647 
649 // Includes inline functions.
650 //#include "CorrectedNormalCurrentFormula.ih"
651 
652 // //
654 
655 #endif // !defined CorrectedNormalCurrentFormula_h
656 
657 #undef CorrectedNormalCurrentFormula_RECURSES
658 #endif // else defined(CorrectedNormalCurrentFormula_RECURSES)
659 
DGtal::SimpleMatrix
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:76
DGtal::SimpleMatrix::setComponent
void setComponent(const DGtal::Dimension i, const DGtal::Dimension j, const Component &aValue)
DGtal::SphericalTriangle
Aim: Represent a triangle drawn onto a sphere of radius 1.
Definition: SphericalTriangle.h:62
DGtal::Trace::warning
std::ostream & warning()
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::crossProduct
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::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:136
DGtal::trace
Trace trace
Definition: Common.h:153
DGtal::CorrectedNormalCurrentFormula
Aim: A helper class that provides static methods to compute corrected normal current formulas of curv...
Definition: CorrectedNormalCurrentFormula.h:100
DGtal::CorrectedNormalCurrentFormula::averageUnitVector
static RealVector averageUnitVector(const RealVectors &vecs)
Definition: CorrectedNormalCurrentFormula.h:633
DGtal::CorrectedNormalCurrentFormula::muXYInterpolatedU
static RealTensor muXYInterpolatedU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &ua, const RealVector &ub, const RealVector &uc, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:544
DGtal::CorrectedNormalCurrentFormula::mu2InterpolatedU
static Scalar mu2InterpolatedU(const RealPoints &pts, const RealVectors &u, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:446
DGtal::CorrectedNormalCurrentFormula::mu2InterpolatedU
static Scalar mu2InterpolatedU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &ua, const RealVector &ub, const RealVector &uc, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:404
DGtal::CorrectedNormalCurrentFormula::mu1InterpolatedU
static Scalar mu1InterpolatedU(const RealPoints &pts, const RealVectors &u, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:307
DGtal::CorrectedNormalCurrentFormula::mu1ConstantU
static Scalar mu1ConstantU(const RealPoints &pts, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:291
DGtal::CorrectedNormalCurrentFormula::mu2ConstantU
static Scalar mu2ConstantU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:340
DGtal::CorrectedNormalCurrentFormula::barycenter
static RealPoint barycenter(const RealPoints &pts)
Definition: CorrectedNormalCurrentFormula.h:621
DGtal::CorrectedNormalCurrentFormula::muXYInterpolatedU
static RealTensor muXYInterpolatedU(const RealPoints &pts, const RealVectors &u, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:592
DGtal::CorrectedNormalCurrentFormula::mu2ConstantU
static Scalar mu2ConstantU(const RealPoints &pts, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:431
DGtal::CorrectedNormalCurrentFormula::mu1ConstantUAtEdge
static Scalar mu1ConstantUAtEdge(const RealPoint &a, const RealPoint &b, const RealVector &ur, const RealVector &ul)
Definition: CorrectedNormalCurrentFormula.h:227
DGtal::CorrectedNormalCurrentFormula::mu1InterpolatedU
static Scalar mu1InterpolatedU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &ua, const RealVector &ub, const RealVector &uc, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:272
DGtal::CorrectedNormalCurrentFormula::muXYConstantUAtEdge
static RealTensor muXYConstantUAtEdge(const RealPoint &a, const RealPoint &b, const RealVector &ur, const RealVector &ul)
Definition: CorrectedNormalCurrentFormula.h:503
DGtal::CorrectedNormalCurrentFormula::muXYConstantU
static RealTensor muXYConstantU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:479
DGtal::CorrectedNormalCurrentFormula::mu2ConstantUAtVertex
static Scalar mu2ConstantUAtVertex(const RealPoint &a, const RealVectors &vu)
Definition: CorrectedNormalCurrentFormula.h:361
DGtal::CorrectedNormalCurrentFormula::mu0InterpolatedU
static Scalar mu0InterpolatedU(const RealPoints &pts, const RealVectors &u, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:186
DGtal::CorrectedNormalCurrentFormula::mu1ConstantU
static Scalar mu1ConstantU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:249
DGtal::CorrectedNormalCurrentFormula::mu0ConstantU
static Scalar mu0ConstantU(const RealPoints &pts, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:165
DGtal::CorrectedNormalCurrentFormula::muXYConstantU
static RealTensor muXYConstantU(const RealPoints &pts, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:577
DGtal::CorrectedNormalCurrentFormula::mu0InterpolatedU
static Scalar mu0InterpolatedU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &ua, const RealVector &ub, const RealVector &uc, bool unit_u=false)
Definition: CorrectedNormalCurrentFormula.h:144
DGtal::CorrectedNormalCurrentFormula::mu0ConstantU
static Scalar mu0ConstantU(const RealPoint &a, const RealPoint &b, const RealPoint &c, const RealVector &u)
Definition: CorrectedNormalCurrentFormula.h:124