NormalCycleFormula.h

 
#pragma once

 
#if defined(NormalCycleFormula_RECURSES)
#error Recursive header files inclusion detected in NormalCycleFormula.h
#else // defined(NormalCycleFormula_RECURSES)
#define NormalCycleFormula_RECURSES
 
#if !defined NormalCycleFormula_h
#define NormalCycleFormula_h

// Inclusions
#include <iostream>
#include <map>
#include "DGtal/base/Common.h"
#include "DGtal/math/linalg/SimpleMatrix.h"

 
namespace DGtal
{

  // class NormalCycleFormula
  template < typename TRealPoint, typename TRealVector >
  struct NormalCycleFormula
  {
    typedef TRealPoint                     RealPoint;
    typedef TRealVector                    RealVector;
    typedef typename RealVector::Component Scalar;
    typedef std::vector< RealPoint >       RealPoints;
    typedef std::vector< RealVector >      RealVectors;
    typedef SimpleMatrix< Scalar, 3, 3 >   RealTensor;
    typedef std::size_t                    Size;
    typedef std::size_t                    Index;
    static const Dimension dimension = RealPoint::dimension;
 
    //-------------------------------------------------------------------------
  public:

    static
    Scalar area( const RealPoints& pts )
    {
      if ( pts.size() <  3 ) return 0.0;
      if ( pts.size() == 3 )
        return area( pts[ 0 ], pts[ 1 ], pts[ 2 ] );
      const RealPoint b = barycenter( pts );
      Scalar          a = 0.0;
      for ( Index i = 0; i < pts.size(); i++ )
        a += area( b, pts[ i ], pts[ (i+1)%pts.size() ] );
      return a;
    }

    static
    Scalar twiceMeanCurvature
    ( const RealPoint& a, const RealPoint& b,
      const RealVector& right, const RealVector& left )
    {
      const RealVector diedre = right.crossProduct( left );
      const Scalar n = std::min( 1.0, std::max( diedre.norm(), 0.0 ) );
      const Scalar angle = ( diedre.dot( b - a) < 0.0 )
        ? asin( n ) : - asin( n );
      return ( b - a ).norm() * angle;
    }

    static
    Scalar gaussianCurvature
    ( const RealPoint& a, const RealPoints& vtcs )
    {
      Scalar angle_sum = 0.0;
      for ( Size i = 0; i < vtcs.size(); i++ )
        angle_sum += acos( (vtcs[i] - a).getNormalized()
                           .dot( ( vtcs[(i+1)%vtcs.size()] - a ).getNormalized() ) );
      return 2.0 * M_PI - angle_sum;
    }

    static
    Scalar gaussianCurvatureWithPairs
    ( const RealPoint& a, const RealPoints& pairs )
    {
      Scalar angle_sum = 0.0;
      for ( Size i = 0; i < pairs.size(); i += 2 )
        angle_sum += acos( ( pairs[i] - a ).getNormalized()
                           .dot( ( pairs[i+1] - a ).getNormalized() ) );
      return 2.0 * M_PI - angle_sum;
    }
    
    static
    RealTensor anisotropicCurvatureH1
    ( const RealPoint& a, const RealPoint& b,
      const RealVector& right, const RealVector& left )
    {
      const RealVector diedre = right.crossProduct( left );
      const Scalar     length = std::max( 0.0, std::min( 1.0, diedre.norm() ) );
      const Scalar      angle = ( diedre.dot( b - a) > 0.0 )
        ? asin( length ) : - asin( length );
      RealVector          e_p = right + left;
      RealVector          e_m = right - left;
      const Scalar norm_e_p   = e_p.norm();
      const Scalar norm_e_m   = e_m.norm();
      e_p = norm_e_p > 1e-10 ? e_p / norm_e_p : RealVector::zero;
      e_m = norm_e_m > 1e-10 ? e_m / norm_e_m : RealVector::zero;
      const RealTensor T_p  =
        { e_p[ 0 ] * e_p[ 0 ], e_p[ 0 ] * e_p[ 1 ], e_p[ 0 ] * e_p[ 2 ],
          e_p[ 1 ] * e_p[ 0 ], e_p[ 1 ] * e_p[ 1 ], e_p[ 1 ] * e_p[ 2 ],
          e_p[ 2 ] * e_p[ 0 ], e_p[ 2 ] * e_p[ 1 ], e_p[ 2 ] * e_p[ 2 ] };
      const RealTensor T_m  =
        { e_m[ 0 ] * e_m[ 0 ], e_m[ 0 ] * e_m[ 1 ], e_m[ 0 ] * e_m[ 2 ],
          e_m[ 1 ] * e_m[ 0 ], e_m[ 1 ] * e_m[ 1 ], e_m[ 1 ] * e_m[ 2 ],
          e_m[ 2 ] * e_m[ 0 ], e_m[ 2 ] * e_m[ 1 ], e_m[ 2 ] * e_m[ 2 ] };
      return 0.5 * ( b - a ).norm()
        * ( ( angle - sin( angle ) ) * T_p + ( angle + sin( angle ) ) * T_m );
    }

    static
    RealTensor anisotropicCurvatureH2
    ( const RealPoint& a, const RealPoint& b,
      const RealVector& right, const RealVector& left )
    {
      const RealVector diedre = right.crossProduct( left );
      const Scalar     length = std::max( 0.0, std::min( 1.0, diedre.norm() ) );
      const Scalar      angle = ( diedre.dot( b - a) > 0.0 )
        ? asin( length ) : - asin( length );
      const Scalar norm_ab = (b - a).norm();
      const RealVector   e = norm_ab > 1e-10 ? (b - a) / norm_ab : RealVector::zero;
      const RealTensor   T =
        { e[ 0 ] * e[ 0 ], e[ 0 ] * e[ 1 ], e[ 0 ] * e[ 2 ],
          e[ 1 ] * e[ 0 ], e[ 1 ] * e[ 1 ], e[ 1 ] * e[ 2 ],
          e[ 2 ] * e[ 0 ], e[ 2 ] * e[ 1 ], e[ 2 ] * e[ 2 ] };
      return ( 0.5 * norm_ab * angle ) * T; // JOL * 0.5
    }
 
    
    //-------------------------------------------------------------------------
  public:
    
    static 
    RealPoint barycenter( const RealPoints& pts )
    {
      RealPoint b;
      for ( auto p : pts ) b += p;
      b /= pts.size();
      return b;
    }

    static
    RealVector normal( const RealPoint& a, const RealPoint& b, const RealPoint& c )
    {
      return ( ( b - a ).crossProduct( c - a ) ).getNormalized();
    }    

    static
    Scalar area( const RealPoint& a, const RealPoint& b, const RealPoint& c )
    {
      return 0.5 * ( ( b - a ).crossProduct( c - a ) ).norm();
    }    

    static 
    RealVector averageUnitVector( const RealVectors& vecs )
    {
      RealVector avg;
      for ( auto v : vecs ) avg += v;
      auto avg_norm = avg.norm();
      return avg_norm != 0.0 ? avg / avg_norm : avg;
    }
    
    
    
  };
 
} // namespace DGtal
 

// Includes inline functions.
//#include "NormalCycleFormula.ih"
 
//                                                                           //
 
#endif // !defined NormalCycleFormula_h
 
#undef NormalCycleFormula_RECURSES
#endif // else defined(NormalCycleFormula_RECURSES)