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 WITH_CGAL
#error You need to have activated CGAL (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)