VCMGeometricFunctors.h

#pragma once

#if defined(VCMGeometricFunctors_RECURSES)
#error Recursive header files inclusion detected in VCMGeometricFunctors.h
#else // defined(VCMGeometricFunctors_RECURSES)

#define VCMGeometricFunctors_RECURSES

#if !defined VCMGeometricFunctors_h

#define VCMGeometricFunctors_h

// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"

// @since 0.8 In DGtal::functors
namespace DGtal {
  namespace functors {

    template <typename TVCMOnDigitalSurface>
    struct VCMNormalVectorFunctor {
      typedef TVCMOnDigitalSurface VCMOnDigitalSurface;
      typedef typename VCMOnDigitalSurface::KSpace KSpace;
      typedef typename VCMOnDigitalSurface::Surfel Surfel;
      typedef typename VCMOnDigitalSurface::VectorN RealVector;
      typedef typename RealVector::Component Scalar;
      typedef Surfel Argument;
      typedef RealVector Quantity;
      
      VCMNormalVectorFunctor( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface = 0 ) 
        : myVCMOnDigitalSurface( aVCMOnDigitalSurface ) {}
      
      void attach( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface )
      {
        myVCMOnDigitalSurface = aVCMOnDigitalSurface;
      }

      Quantity operator()( const Surfel & s ) const
      {
        typedef typename VCMOnDigitalSurface::Surfel2Normals Surfel2Normals;
        ASSERT( myVCMOnDigitalSurface != 0 );
        typename Surfel2Normals::const_iterator itSN = myVCMOnDigitalSurface->mapSurfel2Normals().find( s );
        ASSERT( itSN != myVCMOnDigitalSurface->mapSurfel2Normals().end() );
        return - itSN->second.vcmNormal;
      }

    private:
      CountedConstPtrOrConstPtr<VCMOnDigitalSurface> myVCMOnDigitalSurface;
    };
    

    template <typename TVCMOnDigitalSurface>
    struct VCMAbsoluteCurvatureFunctor {
      typedef TVCMOnDigitalSurface VCMOnDigitalSurface;
      typedef typename VCMOnDigitalSurface::KSpace KSpace;
      typedef typename VCMOnDigitalSurface::Surfel Surfel;
      typedef typename VCMOnDigitalSurface::VectorN RealVector;
      typedef typename RealVector::Component Scalar;
      typedef Surfel Argument;
      typedef Scalar Quantity;
      
      VCMAbsoluteCurvatureFunctor( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface = 0 ) 
        : myVCMOnDigitalSurface( aVCMOnDigitalSurface ) 
      {
        BOOST_STATIC_ASSERT(( KSpace::dimension == 2 ));
      }
      
      void attach( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface )
      {
        myVCMOnDigitalSurface = aVCMOnDigitalSurface;
      }

      Quantity operator()( const Surfel & s ) const
      {
        ASSERT( myVCMOnDigitalSurface != 0 );
        RealVector lambda;
        bool ok = myVCMOnDigitalSurface->getChiVCMEigenvalues( lambda, s );
        ASSERT( ok ); boost::ignore_unused_variable_warning( ok );

        // The last eigenvalue l1 is approximately the mixed "area" 8pi R^3 r / 3
        // The curvature is related to the first eigenvalue l0.
        // k1^2 = 3*l0 / (l1*r^2)
        return sqrt( 3.0 * lambda[0] / lambda[1] ) / myVCMOnDigitalSurface->r();
      }

    private:
      CountedConstPtrOrConstPtr<VCMOnDigitalSurface> myVCMOnDigitalSurface;
    };

    template <typename TVCMOnDigitalSurface>
    struct VCMFirstPrincipalAbsoluteCurvatureFunctor {
      typedef TVCMOnDigitalSurface VCMOnDigitalSurface;
      typedef typename VCMOnDigitalSurface::KSpace KSpace;
      typedef typename VCMOnDigitalSurface::Surfel Surfel;
      typedef typename VCMOnDigitalSurface::VectorN RealVector;
      typedef typename RealVector::Component Scalar;
      typedef Surfel Argument;
      typedef Scalar Quantity;
      
      VCMFirstPrincipalAbsoluteCurvatureFunctor( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface = 0 ) 
        : myVCMOnDigitalSurface( aVCMOnDigitalSurface ) 
      {
        BOOST_STATIC_ASSERT(( KSpace::dimension == 3 ));
      }
      
      void attach( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface )
      {
        myVCMOnDigitalSurface = aVCMOnDigitalSurface;
      }

      Quantity operator()( const Surfel & s ) const
      {
        ASSERT( myVCMOnDigitalSurface != 0 );
        RealVector lambda;
        bool ok = myVCMOnDigitalSurface->getChiVCMEigenvalues( lambda, s );
        ASSERT( ok ); boost::ignore_unused_variable_warning( ok );
        
        // The last eigenvalue l2 is approximately the mixed "area" 2pi R^3 r^2 / 3
        // The greatest principal curvature is related to the second eigenvalue l1.
        // k1^2 = 4*l1 / (l2*r^2)
        return 2.0 * sqrt( lambda[1] / lambda[2] ) / myVCMOnDigitalSurface->r();
      }

    private:
      CountedConstPtrOrConstPtr<VCMOnDigitalSurface> myVCMOnDigitalSurface;
    };

    template <typename TVCMOnDigitalSurface>
    struct VCMSecondPrincipalAbsoluteCurvatureFunctor {
      typedef TVCMOnDigitalSurface VCMOnDigitalSurface;
      typedef typename VCMOnDigitalSurface::KSpace KSpace;
      typedef typename VCMOnDigitalSurface::Surfel Surfel;
      typedef typename VCMOnDigitalSurface::VectorN RealVector;
      typedef typename RealVector::Component Scalar;
      typedef Surfel Argument;
      typedef Scalar Quantity;
      
      VCMSecondPrincipalAbsoluteCurvatureFunctor( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface = 0 ) 
        : myVCMOnDigitalSurface( aVCMOnDigitalSurface ) 
      {
        BOOST_STATIC_ASSERT(( KSpace::dimension == 3 ));
      }
      
      void attach( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface )
      {
        myVCMOnDigitalSurface = aVCMOnDigitalSurface;
      }

      Quantity operator()( const Surfel & s ) const
      {
        ASSERT( myVCMOnDigitalSurface != 0 );
        RealVector lambda;
        bool ok = myVCMOnDigitalSurface->getChiVCMEigenvalues( lambda, s );
        ASSERT( ok ); boost::ignore_unused_variable_warning( ok );
        
        // The last eigenvalue l2 is approximately the mixed "area" 2pi R^3 r^2 / 3
        // The smallest principal curvature is related to the first eigenvalue l0.
        // k2^2 = 4*l0 / (l2*r^2)
        return 2.0 * sqrt( lambda[0] / lambda[2] ) / myVCMOnDigitalSurface->r();
      }

    private:
      CountedConstPtrOrConstPtr<VCMOnDigitalSurface> myVCMOnDigitalSurface;
    };


    template <typename TVCMOnDigitalSurface>
    struct VCMMeanAbsoluteCurvatures3DFunctor {
      typedef TVCMOnDigitalSurface VCMOnDigitalSurface;
      typedef typename VCMOnDigitalSurface::KSpace KSpace;
      typedef typename VCMOnDigitalSurface::Surfel Surfel;
      typedef typename VCMOnDigitalSurface::VectorN RealVector;
      typedef typename RealVector::Component Scalar;
      typedef Surfel Argument;
      typedef Scalar Quantity;
      
      VCMMeanAbsoluteCurvatures3DFunctor( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface = 0 ) 
        : myVCMOnDigitalSurface( aVCMOnDigitalSurface ) 
      {
        BOOST_STATIC_ASSERT(( KSpace::dimension == 3 ));
      }
      
      void attach( ConstAlias<VCMOnDigitalSurface> aVCMOnDigitalSurface )
      {
        myVCMOnDigitalSurface = aVCMOnDigitalSurface;
      }

      Quantity operator()( const Surfel & s ) const
      {
        ASSERT( myVCMOnDigitalSurface != 0 );
        RealVector lambda;
        bool ok = myVCMOnDigitalSurface->getChiVCMEigenvalues( lambda, s );
        ASSERT( ok ); boost::ignore_unused_variable_warning( ok );

        // The last eigenvalue l2 is approximately the mixed "area" 2pi R^3 r^2 / 3
        // The greatest principal curvature is related to the second eigenvalue l1.
        // k1^2 = 4*l1 / (l2*r^2)
        Quantity k1 = 2.0 * sqrt( lambda[1] / lambda[2] ) / myVCMOnDigitalSurface->r();

       // The last eigenvalue l2 is approximately the mixed "area" 2pi R^3 r^2 / 3
        // The smallest principal curvature is related to the first eigenvalue l0.
        // k2^2 = 4*l0 / (l2*r^2)
        Quantity k2 = 2.0 * sqrt( lambda[0] / lambda[2] ) / myVCMOnDigitalSurface->r();
        return (k1+k2) * 0.5;
      }

    private:
      CountedConstPtrOrConstPtr<VCMOnDigitalSurface> myVCMOnDigitalSurface;
    };

  } //namespace functors
} // namespace DGtal


// Includes inline functions.

//                                                                           //

#endif // !defined VCMGeometricFunctors_h

#undef VCMGeometricFunctors_RECURSES
#endif // else defined(VCMGeometricFunctors_RECURSES)