CorrectedNormalCurrentFormula.h

 
#pragma once

 
#if defined(CorrectedNormalCurrentFormula_RECURSES)
#error Recursive header files inclusion detected in CorrectedNormalCurrentFormula.h
#else // defined(CorrectedNormalCurrentFormula_RECURSES)
#define CorrectedNormalCurrentFormula_RECURSES
 
#if !defined CorrectedNormalCurrentFormula_h
#define CorrectedNormalCurrentFormula_h

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

 
namespace DGtal
{

  // class CorrectedNormalCurrentFormula
  template < typename TRealPoint, typename TRealVector >
  struct CorrectedNormalCurrentFormula
  {
    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                    Index;
    static const Dimension dimension = RealPoint::dimension;
 
    //-------------------------------------------------------------------------
  public:
    
    static
    Scalar mu0ConstantU
    ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
      const RealVector& u )
    {
      return 0.5 * ( b - a ).crossProduct( c - a ).dot( u );
    }

    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 )
    {
      // MU0=1/2*det( uM, B-A, C-A )
      //    =  1/2 < ( (u_A + u_B + u_C)/3.0 ) | (AB x AC ) >
      RealVector uM = ( ua+ub+uc ) / 3.0;
      if ( unit_u )
        {
          auto uM_norm = uM.norm();
          uM = uM_norm == 0.0 ? uM : uM / uM_norm;
        }
      return 0.5 * ( b - a ).crossProduct( c - a ).dot( uM );
    }
    
    static
    Scalar mu0ConstantU( const RealPoints& pts, const RealVector& u )
    {
      if ( pts.size() <  3 ) return 0.0;
      if ( pts.size() == 3 )
        return mu0ConstantU( pts[ 0 ], pts[ 1 ], pts[ 2 ], u );
      const RealPoint b = barycenter( pts );
      Scalar          a = 0.0;
      for ( Index i = 0; i < pts.size(); i++ )
        a += mu0ConstantU( b, pts[ i ], pts[ (i+1)%pts.size() ], u );
      return a;
    }

    static
    Scalar mu0InterpolatedU( const RealPoints& pts, const RealVectors& u,
                             bool unit_u = false )
    {
      ASSERT( pts.size() == u.size() );
      if ( pts.size() <  3 ) return 0.0;
      if ( pts.size() == 3 )
        return mu0InterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
                                 u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
      const RealPoint   b = barycenter( pts );
      const RealVector ub = averageUnitVector( u );
      Scalar          a = 0.0;
      for ( Index i = 0; i < pts.size(); i++ )
        a += mu0InterpolatedU( b,  pts[ i ], pts[ (i+1)%pts.size() ],
                               ub,   u[ i ],   u[ (i+1)%pts.size() ], unit_u );
      return a;
    }

    
    //-------------------------------------------------------------------------
  public:

    static
    Scalar mu1ConstantUAtEdge
    ( const RealPoint& a, const RealPoint& b, 
      const RealVector& ur, const RealVector& ul )
    {
      RealVector   e  = b - a;
      RealVector   e1 = ur.crossProduct( ul );
      const Scalar l1 = std::min( e1.norm(), 1.0 );
      if ( l1 < 1e-10 ) return 0.0;
      e1 /= l1;
      const Scalar Psi = asin( l1 );
      return -Psi * e.dot( e1 );
    }

    static
    Scalar mu1ConstantU
    ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
      const RealVector& u )
    {
      (void)a; (void)b; (void)c; (void)u;
      return 0.0;
    }

    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 )
    {
      // MU1=1/2( | uM u_C-u_B A | + | uM u_A-u_C B | + | uM u_B-u_A C |
      RealVector uM = ( ua+ub+uc ) / 3.0;
      if ( unit_u ) uM /= uM.norm();
      return 0.5 * ( uM.crossProduct( uc - ub ).dot( a )
                     + uM.crossProduct( ua - uc ).dot( b )
                     + uM.crossProduct( ub - ua ).dot( c ) );
    }
    
    static
    Scalar mu1ConstantU( const RealPoints& pts, const RealVector& u )
    {
      (void)pts; //not used parameter
      (void)u;
      return 0.0;
    }

    static
    Scalar mu1InterpolatedU( const RealPoints& pts, const RealVectors& u,
                             bool unit_u = false )
    {
      ASSERT( pts.size() == u.size() );
      if ( pts.size() <  3 ) return 0.0;
      if ( pts.size() == 3 )
        return mu1InterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
                                 u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
      const RealPoint   b = barycenter( pts );
      const RealVector ub = averageUnitVector( u );
      Scalar          a = 0.0;
      for ( Index i = 0; i < pts.size(); i++ )
        a += mu1InterpolatedU( b,  pts[ i ], pts[ (i+1)%pts.size() ],
                               ub,   u[ i ],   u[ (i+1)%pts.size() ], unit_u );
      return a;
    }

 
    //-------------------------------------------------------------------------
  public:

    static
    Scalar mu2ConstantU
    ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
      const RealVector& u )
    {
      (void)a; (void)b; (void)c; (void)u;
      return 0.0;
    }

    static
    Scalar mu2ConstantUAtVertex
    ( const RealPoint& a, const RealVectors& vu )
    {
      typedef SpaceND< dimension >     Space;
      typedef SphericalTriangle<Space> ST;
      RealVector avg_u;
      for ( const auto& u : vu ) avg_u += u;
      const Scalar l = avg_u.norm();
      if ( l < 1e-10 )
        {
          trace.warning() << "[CorrectedNormalCurrentFormula::mu2ConstantUAtVertex]"
                          << " Invalid surrounding normal vectors at vertex "
                          << a << "."
                          << " Unit normal vectors should lie strictly in the same"
                          << " hemisphere." << std::endl;
          return 0.0;
        }
      avg_u /= l;
      Scalar mu2 = 0.0;
      const auto n = vu.size();
      for ( size_t i = 0; i < n; ++i )
        {
          ST st( avg_u, vu[ i ], vu[ (i+1) % n ] );
          mu2 += st.algebraicArea();
        }
      return mu2;
    }

    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 )
    {
      (void)a; //not used
      (void)b;
      (void)c;
      
      // Using non unitary interpolated normals give
      // MU2=1/2*det( uA, uB, uC )
      // When normals are unitary, it is the area of a spherical triangle.
      if ( unit_u )
        {
          typedef SpaceND< dimension > Space;
          SphericalTriangle<Space> ST( ua, ub, uc ); // check order.
          return ST.algebraicArea();
        }
      else
        return 0.5 * ( ua.crossProduct( ub ).dot( uc ) );
    }
    
    static
    Scalar mu2ConstantU( const RealPoints& pts, const RealVector& u )
    {
      (void) pts; (void) u;
      return 0.0;
    }
    
    static
    Scalar mu2InterpolatedU( const RealPoints& pts, const RealVectors& u,
                             bool unit_u = false )
    {
      ASSERT( pts.size() == u.size() );
      if ( pts.size() <  3 ) return 0.0;
      if ( pts.size() == 3 )
        return mu2InterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
                                 u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
      const RealPoint   b = barycenter( pts );
      const RealVector ub = averageUnitVector( u );
      Scalar          a = 0.0;
      for ( Index i = 0; i < pts.size(); i++ )
        a += mu2InterpolatedU( b,  pts[ i ], pts[ (i+1)%pts.size() ],
                               ub,   u[ i ],   u[ (i+1)%pts.size() ], unit_u );
      return a;
    }
    
    
    //-------------------------------------------------------------------------
  public:
    
    static
    RealTensor muXYConstantU
    ( const RealPoint& a, const RealPoint& b, const RealPoint& c,
      const RealVector& u )
    {
      (void)a; (void)b; (void)c; (void)u;
      return RealTensor { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    }

    static
    RealTensor muXYConstantUAtEdge
    ( const RealPoint& a, const RealPoint& b, 
      const RealVector& ur, const RealVector& ul )
    {
      RealTensor M { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
      RealVector   e  = b - a;
      RealVector   e1 = ur.crossProduct( ul );
      const Scalar l1 = std::min( e1.norm(), 1.0 );
      if ( l1 < 1e-10 ) return M;
      e1 /= l1;
      const Scalar             psi = asin( l1 );
      const Scalar         sin_psi = sin( psi );
      const Scalar        sin_2psi = sin( 2.0 * psi );
      const RealVector      e_x_ur = e. crossProduct( ur );
      const RealVector     e1_x_ur = e1.crossProduct( ur );
      const RealVector e_x_e1_x_ur = e. crossProduct( e1_x_ur );
      for ( Dimension i = 0; i < 3; i++ )
        for ( Dimension j = 0; j < 3; j++ )
          {
            Scalar v = (-psi - 0.5*sin_2psi) * e_x_ur[ i ] * e1_x_ur[ j ]
              + (sin_psi*sin_psi) * ( e_x_ur[ i ] * ur[ j ]
                                      - e_x_e1_x_ur[ i ] * e1_x_ur[ j ] )
              + (psi - 0.5*sin_2psi) * e_x_e1_x_ur[ i ] * ur[ j ];
            M.setComponent( i, j, -0.5 * v );
          }
      return M;
    }
    
    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 )
    {
      RealTensor T;
      //  MUXY = 1/2 < uM | < Y | uc-ua > X x (b-a) - < Y | ub-ua > X x (c-a) >
      //  MUXY = 1/2 ( < Y | ub-ua > | X uM (c-a) | - < Y | uc-ua > | X uM (b-a) | )
      RealVector uM = ( ua+ub+uc ) / 3.0;
      if ( unit_u ) uM /= uM.norm();
      const RealVector uac = uc - ua;
      const RealVector uab = ub - ua;
      const RealVector  ab = b - a;
      const RealVector  ac = c - a;
      for ( Dimension i = 0; i < 3; ++i ) {
        RealVector X = RealVector::base( i, 1.0 );
        for ( Dimension j = 0; j < 3; ++j ) {
          // Since RealVector Y = RealVector::base( j, 1.0 );
          // < Y | uac > = uac[ j ]
          const Scalar tij =
            0.5 * uM.dot( uac[ j ] * X.crossProduct( ab )
                          - uab[ j ] * X.crossProduct( ac ) );
          T.setComponent( i, j, tij );
        }
      }
      return T;
    }
    
    static
    RealTensor muXYConstantU( const RealPoints& pts, const RealVector& u )
    {
      (void)pts; (void)u;
      return RealTensor { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
    }

    static
    RealTensor muXYInterpolatedU( const RealPoints& pts, const RealVectors& u,
                                  bool unit_u = false )
    {
      ASSERT( pts.size() == u.size() );
      if ( pts.size() <  3 ) return RealTensor { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
      if ( pts.size() == 3 )
        return muXYInterpolatedU( pts[ 0 ], pts[ 1 ], pts[ 2 ],
                                 u[ 0 ], u[ 1 ], u[ 2 ], unit_u );
      const RealPoint   b = barycenter( pts );
      const RealVector ub = averageUnitVector( u );
      RealTensor        T = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
      for ( Index i = 0; i < pts.size(); i++ )
        T += muXYInterpolatedU( b,  pts[ i ], pts[ (i+1)%pts.size() ],
                                ub,   u[ i ],   u[ (i+1)%pts.size() ], unit_u );
      return T;
    }
    
 
    
    //-------------------------------------------------------------------------
  public:
    
    static 
    RealPoint barycenter( const RealPoints& pts )
    {
      RealPoint b;
      for ( auto p : pts ) b += p;
      b /= pts.size();
      return b;
    }

    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 "CorrectedNormalCurrentFormula.ih"
 
//                                                                           //
 
#endif // !defined CorrectedNormalCurrentFormula_h
 
#undef CorrectedNormalCurrentFormula_RECURSES
#endif // else defined(CorrectedNormalCurrentFormula_RECURSES)