TensorVotingFeatureExtraction.h

 
#pragma once

 
#if defined(TensorVotingFeatureExtraction_RECURSES)
#error Recursive header files inclusion detected in TensorVotingFeatureExtraction.h
#else // defined(TensorVotingFeatureExtraction_RECURSES)
#define TensorVotingFeatureExtraction_RECURSES
 
#if !defined TensorVotingFeatureExtraction_h
#define TensorVotingFeatureExtraction_h

// Inclusions
#include <iostream>
#include <DGtal/base/Common.h>
#include <DGtal/topology/SCellsFunctors.h>
#include <vector>
#include "DGtal/math/linalg/SimpleMatrix.h"
#include "DGtal/math/linalg/EigenDecomposition.h"
 
namespace DGtal
{
  namespace functors
  {
    // template class TensorVotingFeatureExtraction
    template <typename TSurfel, typename TEmbedder>
    class TensorVotingFeatureExtraction
    {
    public:
 
      typedef TSurfel Surfel;
      typedef TEmbedder SCellEmbedder;
      typedef typename SCellEmbedder::RealPoint RealPoint;
      typedef double Quantity;

      TensorVotingFeatureExtraction(ConstAlias<SCellEmbedder> anEmbedder,
                                    const double h):
        myEmbedder(&anEmbedder), myH(h)
      {
        myId.identity();
        myArea = 0.0;
        myFirstSurfel = true;
        myAccum.constant(0.0);
      }

      void pushSurfel(const Surfel & aSurf,
                      const double aDistance)
      {
        if (myFirstSurfel)
          {
            myReceiver =  myEmbedder->operator()(aSurf);
            myFirstSurfel = false;
          }
        else
          {
            myArea+= aDistance;
            const RealPoint p = myEmbedder->operator()(aSurf);
            const RealPoint v = p - myReceiver;
 
            double maxcol = 0.0;
            double matnorm = 0.0;
 
            // I - vv^t/||vv^t||
            for(DGtal::Dimension i= 0; i <3; i++)
              {
                maxcol = 0.0;
                for(DGtal::Dimension j=0; j < 3; j++)
                  {
                    myVote.setComponent(i,j, v(i)*v(j));
                    if (std::abs(v(i)*v(j)) > maxcol)
                      maxcol = std::abs(v(i)*v(j));
                  }
                matnorm += maxcol;
              }
 
            myAccum += (myId - myVote/matnorm)*aDistance;
          }
      }

      Quantity eval( )
      {
        SimpleMatrix<double, 3, 3> eigenvectors;
        RealPoint eigenvalues;
 
        myAccum /= myArea;
        EigenDecomposition<3, double>::getEigenDecomposition( myAccum, eigenvectors, eigenvalues);
 
#ifdef DEBUG
        for( Dimension i_dim = 1; i_dim < 3; ++i_dim )
            ASSERT ( std::abs(eigenvalues[i_dim - 1]) <= std::abs(eigenvalues[i_dim]) );
#endif
 
        return ((eigenvalues[0] + eigenvalues[1])/(eigenvalues[2]));
      }

      void reset()
      {
        myArea = 0.0;
        myFirstSurfel = true;
        myAccum.constant(0.0);
      }
 
 
    private:

      bool myFirstSurfel;

      const SCellEmbedder * myEmbedder;
 
      //Grid step
      double myH;

      RealPoint myReceiver;

      SimpleMatrix<double, 3, 3> myAccum, myId, myVote;

      double myArea;
 
    }; // end of class TensorVotingFeatureExtraction
  }//namespace functors
}// namespace DGtal
 
 
//                                                                           //
 
#endif // !defined TensorVotingFeatureExtraction_h
 
#undef TensorVotingFeatureExtraction_RECURSES
#endif // else defined(TensorVotingFeatureExtraction_RECURSES)