NormalCycleComputer.h

 
#pragma once

 
#if defined(NormalCycleComputer_RECURSES)
#error Recursive header files inclusion detected in NormalCycleComputer.h
#else // defined(NormalCycleComputer_RECURSES)
#define NormalCycleComputer_RECURSES
 
#if !defined NormalCycleComputer_h
#define NormalCycleComputer_h

// Inclusions
#include "DGtal/base/Common.h"
#include "DGtal/math/linalg/EigenDecomposition.h"
#include "DGtal/geometry/meshes/SurfaceMeshMeasure.h"
#include "DGtal/geometry/meshes/NormalCycleFormula.h"
#include "DGtal/shapes/SurfaceMesh.h"
 
 
namespace DGtal
{
  // template class NormalCycleComputer
  template < typename TRealPoint, typename TRealVector >
  struct NormalCycleComputer
  {
    typedef TRealPoint                              RealPoint;
    typedef TRealVector                             RealVector;
    typedef NormalCycleComputer< RealPoint, RealVector > Self;
    static const Dimension dimension = RealPoint::dimension;
    BOOST_STATIC_ASSERT( ( dimension == 3 ) );
    typedef DGtal::SurfaceMesh< RealPoint, RealVector >  SurfaceMesh;
    typedef typename RealVector::Component               Scalar;
    typedef SimpleMatrix< Scalar, dimension, dimension > RealTensor;
    typedef NormalCycleFormula< RealPoint, RealVector >  Formula;
    typedef SurfaceMeshMeasure< RealPoint, RealVector, Scalar >     ScalarMeasure;
    typedef SurfaceMeshMeasure< RealPoint, RealVector, RealTensor > TensorMeasure;
    typedef std::vector< Scalar >                        Scalars;
    typedef std::vector< RealPoint >                     RealPoints;
    typedef std::vector< RealVector >                    RealVectors;
    typedef std::vector< RealTensor >                    RealTensors;
    typedef typename SurfaceMesh::Size                   Size; 
    typedef typename SurfaceMesh::Index                  Index;
    typedef typename SurfaceMesh::Vertex                 Vertex;
    typedef typename SurfaceMesh::Face                   Face;
 
    
    NormalCycleComputer( ConstAlias< SurfaceMesh > aMesh );

    ScalarMeasure computeMu0() const;
    ScalarMeasure computeMu1() const;
    ScalarMeasure computeMu2() const;
    TensorMeasure computeMuXY() const;
    TensorMeasure computeMuXYs() const;
    
    //-------------------------------------------------------------------------
  public:

    static
    Scalar meanCurvature( Scalar mu0, Scalar mu1 )
    {
      return ( mu0 != 0.0 ) ? mu1 / ( 2.0 * mu0 ) : 0.0;
    }

    static
    Scalar GaussianCurvature( Scalar mu0, Scalar mu2 )
    {
      return ( mu0 != 0.0 ) ? mu2 / mu0 : 0.0;
    }

    static
    std::tuple< Scalar, Scalar, RealVector, RealVector >
    principalCurvatures( Scalar mu0, RealTensor muXY, const RealVector& N ) 
    {
      muXY += muXY.transpose();
      muXY *= 0.5;
      const double   coef_N = 1000.0 * mu0;
      // Trick to force orthogonality to normal vector.
      for ( int j = 0; j < 3; j++ )
        for ( int k = 0; k < 3; k++ )
          muXY( j, k ) += coef_N * N[ j ] * N[ k ];
      RealTensor V;
      RealVector L;
      EigenDecomposition< 3, double>::getEigenDecomposition( muXY, V, L );
      return std::make_tuple( ( mu0 != 0.0 ) ? -L[ 1 ] / mu0 : 0.0,
                              ( mu0 != 0.0 ) ? -L[ 0 ] / mu0 : 0.0,
                              V.column( 1 ),
                              V.column( 0 ) );
    }    

 
    // ------------------------- Protected Datas ------------------------------
  protected:
    
    const SurfaceMesh& myMesh;
    
 
    
  }; // end of class NormalCycleComputer
    
} // namespace DGtal

// Includes inline functions.
#include "NormalCycleComputer.ih"
//                                                                           //
 
#endif // !defined NormalCycleComputer_h
 
#undef NormalCycleComputer_RECURSES
#endif // else defined(NormalCycleComputer_RECURSES)