IIGeometricFunctors.h

 
#pragma once
 
#if defined(IIGeometricFunctors_RECURSES)
#error Recursive header files inclusion detected in IIGeometricFunctors.h
#else // defined(IIGeometricFunctors_RECURSES)
#define IIGeometricFunctors_RECURSES
 
#if !defined IIGeometricFunctors_h
#define IIGeometricFunctors_h
 
// Inclusions
#include <iostream>
#include <tuple>
#include "DGtal/base/Common.h"
#include "DGtal/math/linalg/EigenDecomposition.h"
#include "DGtal/math/linalg/CMatrix.h"
 
// @since 0.8 In DGtal::functors
namespace DGtal {
  namespace functors {
 
    // template class IINormalDirectionFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    class IINormalDirectionFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IINormalDirectionFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef RealVector Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
 
      IINormalDirectionFunctor() {}
      IINormalDirectionFunctor( const Self& /* other */ ) {}
      Self& operator=( const Self& /* other */ ) { return *this; }
      Value operator()( const Argument& arg ) const
      {
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
#ifdef DEBUG
        for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
        {
          ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
        }
#endif
 
        return eigenVectors.column( 0 ); // normal vector is associated to smallest eigenvalue.
      }
 
      void init( Component /* h */, Component /* r */ ) {}
 
    private:
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IINormalDirectionFunctor
 
 
    // template class IITangentDirectionFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    class IITangentDirectionFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IITangentDirectionFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef RealVector Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 2 ));
 
      IITangentDirectionFunctor() {}
      IITangentDirectionFunctor( const Self& /* other */ ) {}
      Self& operator=( const Self& /* other */ ) { return *this; }
      Value operator()( const Argument& arg ) const
      {
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
#ifdef DEBUG
        for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
        {
          ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
        }
#endif
        return eigenVectors.column( 1 ); // tangent vector is associated to greatest eigenvalue.
      }
 
      void init( Component /* h */, Component /* r */ ) {}
 
    private:
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IITangentDirectionFunctor
 
 
    // template class IIFirstPrincipalDirectionFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    class IIFirstPrincipalDirectionFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIFirstPrincipalDirectionFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef RealVector Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension >= 2 ));
 
      IIFirstPrincipalDirectionFunctor() {}
      IIFirstPrincipalDirectionFunctor( const Self& /* other */ ) {}
      Self& operator=( const Self& /* other */ ) { return *this; }
      Value operator()( const Argument& arg ) const
      {
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
#ifdef DEBUG
        for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
        {
          ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
        }
#endif
 
        return eigenVectors.column( Space::dimension - 1 ); // first principal curvature direction is associated to greatest eigenvalue.
      }
 
      void init( Component /* h */, Component /* r */ ) {}
 
    private:
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IIFirstPrincipalDirectionFunctor
 
 
    // template class IISecondPrincipalDirectionFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    class IISecondPrincipalDirectionFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IISecondPrincipalDirectionFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef RealVector Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension >= 3 ));
 
      IISecondPrincipalDirectionFunctor() {}
      IISecondPrincipalDirectionFunctor( const Self& /* other */ ) {}
      Self& operator=( const Self& /* other */ ) { return *this; }
      Value operator()( const Argument& arg ) const
      {
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
#ifdef DEBUG
        for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
        {
          ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
        }
#endif
 
        return eigenVectors.column( Space::dimension - 2 ); // second principal curvature direction is associated to greatest eigenvalue.
      }
 
      void init( Component /* h */, Component /* r */ ) {}
 
    private:
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IISecondPrincipalDirectionFunctor
 
    // template class IIPrincipalDirectionsFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    class IIPrincipalDirectionsFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIPrincipalDirectionsFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef std::pair<RealVector,RealVector> Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension >= 3 ));
 
      IIPrincipalDirectionsFunctor() {}
      IIPrincipalDirectionsFunctor( const Self& /* other */ ) {}
      Self& operator=( const Self& /* other */ ) { return *this; }
      Value operator()( const Argument& arg ) const
      {
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
#ifdef DEBUG
        for( Dimension i_dim = 1; i_dim < Space::dimension; ++i_dim )
        {
          ASSERT ( std::abs(eigenValues[i_dim - 1]) <= std::abs(eigenValues[i_dim]) );
        }
#endif
 
        return Value(
                eigenVectors.column( Space::dimension - 1 ),
                eigenVectors.column( Space::dimension - 2 )
                );
      }
 
      void init( Component /* h */, Component /* r */ ) {}
 
    private:
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IIPrincipalDirectionsFunctor
 
 
    // template class IIPrincipalCurvaturesAndDirectionsFunctor
    template  <typename TSpace,
               typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension,
               TSpace::dimension> >
    class IIPrincipalCurvaturesAndDirectionsFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIPrincipalCurvaturesAndDirectionsFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef std::tuple<double, double, RealVector, RealVector> Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
 
      Value operator()( const Argument& arg ) const
      {
        Argument cp_arg = arg;
        cp_arg *= dh5;
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( cp_arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
        ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
                && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
 
        Quantity res(d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r,
                     d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r,
                     eigenVectors.column( 1 ),
                     eigenVectors.column( 2 ));
        return res;
      }
 
      void init( Component h , Component  r )
      {
        double r3 = r * r * r;
        double r6 = r3 * r3;
        d6_PIr6 = 6.0 / ( M_PI * r6 );
        d8_5r = 8.0 / ( 5.0 * r );
        double h2 = h * h;
        dh5 = h2 * h2 * h;
       }
 
    private:
      double dh5;
      double d6_PIr6;
      double d8_5r;
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IIPrincipalCurvaturesAndDirectionsFunctor
 
    // template class IICurvatureFunctor
    template  <typename TSpace>
    struct IICurvatureFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IICurvatureFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef Component Argument;
      typedef Component Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 2 ));
 
      Value operator()( const Argument& arg ) const
      {
        Quantity cp_quantity = arg;
        cp_quantity *= dh2;
        return d3_r * ( dPI_2 - d1_r2 * cp_quantity );
      }
 
      void init( Component h, Component r )
      {
        d1_r2 = 1.0 / ( r * r );
        dPI_2 = M_PI / 2.0;
        d3_r = 3.0 / r;
        dh2 = h * h;
      }
 
    private:
      Quantity dh2;
      Quantity d3_r;
      Quantity dPI_2;
      Quantity d1_r2;
    }; // end of class IICurvatureFunctor
 
 
    // template class IIMeanCurvature3DFunctor
    template  <typename TSpace>
    struct IIMeanCurvature3DFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIMeanCurvature3DFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef Component Argument;
      typedef Component Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
 
      Value operator()( const Argument& arg ) const
      {
        Quantity cp_quantity = arg;
        cp_quantity *= dh3;
        return d8_3r - d_4_PIr4 * cp_quantity;
      }
 
      void init( Component h, Component r )
      {
        d8_3r = 8.0 / ( 3.0 * r );
        double r2 = r * r;
        d_4_PIr4 = 4.0 / ( M_PI * r2 * r2 );
        dh3 = h * h * h;
      }
 
    private:
      Quantity dh3;
      Quantity d8_3r;
      Quantity d_4_PIr4;
    }; // end of class IIMeanCurvature3DFunctor
 
    // template class IIGaussianCurvature3DFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    struct IIGaussianCurvature3DFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIGaussianCurvature3DFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef Component Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
 
      Value operator()( const Argument& arg ) const
      {
        Argument cp_arg = arg;
        cp_arg *= dh5;
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( cp_arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
        ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
              && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
 
        Value k1 = d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r;
        Value k2 = d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r;
        return k1 * k2;
      }
 
      void init( Component h, Component r )
      {
        double r3 = r * r * r;
        double r6 = r3 * r3;
        d6_PIr6 = 6.0 / ( M_PI * r6 );
        d8_5r = 8.0 / ( 5.0 * r );
        double h2 = h * h;
        dh5 = h2 * h2 * h;
      }
 
    private:
      Quantity dh5;
      Quantity d6_PIr6;
      Quantity d8_5r;
 
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IIGaussianCurvature3DFunctor
 
    // template class IIFirstPrincipalCurvature3DFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    struct IIFirstPrincipalCurvature3DFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIFirstPrincipalCurvature3DFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef Component Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
 
      Value operator()( const Argument& arg ) const
      {
        Argument cp_arg = arg;
        cp_arg *= dh5;
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( cp_arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
        ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
              && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
 
 
        return d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r;
      }
 
      void init( Component h, Component r )
      {
        double r3 = r * r * r;
        double r6 = r3 * r3;
        d6_PIr6 = 6.0 / ( M_PI * r6 );
        d8_5r = 8.0 / ( 5.0 * r );
        double h2 = h * h;
        dh5 = h2 * h2 * h;
      }
 
    private:
      Quantity dh5;
      Quantity d6_PIr6;
      Quantity d8_5r;
 
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IIFirstPrincipalCurvature3DFunctor
 
    // template class IISecondPrincipalCurvature3DFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    struct IISecondPrincipalCurvature3DFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IISecondPrincipalCurvature3DFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef Component Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
 
      Value operator()( const Argument& arg ) const
      {
        Argument cp_arg = arg;
        cp_arg *= dh5;
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( cp_arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
        ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
              && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
 
        return d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r;
      }
 
      void init( Component h, Component r )
      {
        double r3 = r * r * r;
        double r6 = r3 * r3;
        d6_PIr6 = 6.0 / ( M_PI * r6 );
        d8_5r = 8.0 / ( 5.0 * r );
        double h2 = h * h;
        dh5 = h2 * h2 * h;
      }
 
    private:
      Quantity dh5;
      Quantity d6_PIr6;
      Quantity d8_5r;
 
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IISecondPrincipalCurvature3DFunctor
 
 
    // template class IIPrincipalCurvatures3DFunctor
    template  <typename TSpace, typename TMatrix=SimpleMatrix< typename TSpace::RealVector::Component, TSpace::dimension, TSpace::dimension> >
    struct IIPrincipalCurvatures3DFunctor
    {
      // ----------------------- Standard services ------------------------------
    public:
      typedef IIPrincipalCurvatures3DFunctor<TSpace> Self;
      typedef TSpace Space;
      typedef typename Space::RealVector RealVector;
      typedef typename RealVector::Component Component;
      typedef TMatrix Matrix;
      typedef Matrix Argument;
      typedef std::pair<Component, Component> Quantity;
      typedef Quantity Value;
 
      BOOST_CONCEPT_ASSERT(( concepts::CMatrix<Matrix> ));
      BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
      BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
 
      Value operator()( const Argument& arg ) const
      {
        Argument cp_arg = arg;
        cp_arg *= dh5;
        EigenDecomposition<Space::dimension, Component, Matrix>
          ::getEigenDecomposition( cp_arg, eigenVectors, eigenValues );
 
        ASSERT ( !std::isnan(eigenValues[0]) ); // NaN
        ASSERT ( (std::abs(eigenValues[0]) <= std::abs(eigenValues[1]))
              && (std::abs(eigenValues[1]) <= std::abs(eigenValues[2])) );
 
        return Value(
                d6_PIr6 * ( eigenValues[2] - ( 3.0 * eigenValues[1] )) + d8_5r,
                d6_PIr6 * ( eigenValues[1] - ( 3.0 * eigenValues[2] )) + d8_5r
               );
      }
 
      void init( Component h, Component r )
      {
        double r3 = r * r * r;
        double r6 = r3 * r3;
        d6_PIr6 = 6.0 / ( M_PI * r6 );
        d8_5r = 8.0 / ( 5.0 * r );
        double h2 = h * h;
        dh5 = h2 * h2 * h;
      }
 
    private:
      double dh5;
      double d6_PIr6;
      double d8_5r;
 
      mutable Matrix eigenVectors;
      mutable RealVector eigenValues;
    }; // end of class IIPrincipalCurvatures3DFunctor
 
} // namespace functors
 
} // namespace DGtal
 
 
//                                                                           //
 
#endif // !defined IIGeometricFunctors_h
 
#undef IIGeometricFunctors_RECURSES
#endif // else defined(IIGeometricFunctors_RECURSES)