MongeJetFittingMeanCurvatureEstimator.h

 
#pragma once

 
#if defined(MongeJetFittingMeanCurvatureEstimator_RECURSES)
#error Recursive header files inclusion detected in MongeJetFittingMeanCurvatureEstimator.h
#else // defined(MongeJetFittingMeanCurvatureEstimator_RECURSES)
#define MongeJetFittingMeanCurvatureEstimator_RECURSES
 
#if !defined MongeJetFittingMeanCurvatureEstimator_h
#define MongeJetFittingMeanCurvatureEstimator_h

// Inclusions
#include <iostream>
#include <DGtal/base/Common.h>
#include <DGtal/topology/SCellsFunctors.h>
 
#ifndef  DGTAL_WITH_CGAL
#error You need to have activated CGAL  (DGTAL_WITH_CGAL) to include this file.
#endif
 
//CGAL
#include <CGAL/Cartesian.h>
#include <CGAL/Monge_via_jet_fitting.h>
#include <vector>

 
namespace DGtal
{
  namespace functors
  {
  // template class MongeJetFittingMeanCurvatureEstimator
  template <typename TSurfel, typename TEmbedder>
  class MongeJetFittingMeanCurvatureEstimator
  {
  public:
 
    typedef TSurfel Surfel;
    typedef TEmbedder SCellEmbedder;
    typedef double Quantity;
    typedef typename SCellEmbedder::RealPoint RealPoint;
 
    typedef CGAL::Cartesian<double> CGALKernel;
    typedef CGALKernel::Point_3  CGALPoint;
    typedef CGAL::Monge_via_jet_fitting<CGALKernel>  CGALMongeViaJet;
    typedef CGALMongeViaJet::Monge_form CGALMongeForm;

    MongeJetFittingMeanCurvatureEstimator(ConstAlias<SCellEmbedder> anEmbedder,
                                          const double h, unsigned int d = 4):
      myEmbedder(&anEmbedder), myH(h), myD(d)
    {
      FATAL_ERROR_MSG(d>=2,"Polynomial surface degree must be greater than 2");
    }

    void pushSurfel(const Surfel & aSurf,
                    const double aDistance)
    {
      BOOST_VERIFY(aDistance==aDistance);
 
      RealPoint p = myEmbedder->operator()(aSurf);
      CGALPoint pp(p[0]*myH,p[1]*myH,p[2]*myH);
      myPoints.push_back(pp);
    }

    Quantity eval()
    {
      CGALMongeForm monge_form;
      CGALMongeViaJet monge_fit;
 
      monge_form = monge_fit(myPoints.begin() , myPoints.end(), myD, (2<myD)? myD : 2);
 
      double k1 = monge_form.principal_curvatures ( 0 );
      double k2 = monge_form.principal_curvatures ( 1 );
      return 0.5*(k1+k2);
    }

    void reset()
    {
      myPoints.clear();
    }
 
 
  private:

    const SCellEmbedder * myEmbedder;

    std::vector<CGALPoint> myPoints;

    double myH;

    unsigned int myD;
 
  }; // end of class MongeJetFittingMeanCurvatureEstimator
  }
} // namespace DGtal
 
 
//                                                                           //
 
#endif // !defined MongeJetFittingMeanCurvatureEstimator_h
 
#undef MongeJetFittingMeanCurvatureEstimator_RECURSES
#endif // else defined(MongeJetFittingMeanCurvatureEstimator_RECURSES)