exampleEstimatorFromSurfelFunctors.cpp

 
#include <iostream>
#include "DGtal/base/Common.h"
 
// Shape constructors
#include "DGtal/io/readers/VolReader.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
#include "DGtal/topology/SurfelAdjacency.h"
#include "DGtal/topology/helpers/Surfaces.h"
#include "DGtal/topology/LightImplicitDigitalSurface.h"
#include "DGtal/images/ImageHelper.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/graph/DepthFirstVisitor.h"
#include "DGtal/graph/GraphVisitorRange.h"
#include "DGtal/shapes/Shapes.h"
 
#include "DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h"
#include "DGtal/geometry/volumes/distance/LpMetric.h"
#ifdef WITH_CGAL
#include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingGaussianCurvatureEstimator.h"
#include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingMeanCurvatureEstimator.h"
#include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingNormalVectorEstimator.h"
#include "DGtal/geometry/surfaces/estimation/estimationFunctors/LinearLeastSquareFittingNormalVectorEstimator.h"
#endif
#include "DGtal/geometry/surfaces/estimation/estimationFunctors/ElementaryConvolutionNormalVectorEstimator.h"
#include "DGtal/base/BasicFunctors.h"
 
using namespace std;
using namespace DGtal;
 
 
template <typename TPoint3>
struct ImplicitDigitalEllipse3 {
  typedef TPoint3 Point;
  inline
  ImplicitDigitalEllipse3( double a, double b, double c )
  : myA( a ), myB( b ), myC( c )
  {}
  inline
  bool operator()( const TPoint3 & p ) const
  {
    double x = ( (double) p[ 0 ] / myA );
    double y = ( (double) p[ 1 ] / myB );
    double z = ( (double) p[ 2 ] / myC );
    return ( x*x + y*y + z*z ) <= 1.0;
  }
  double myA, myB, myC;
};
 
 
int main(  )
{
  using namespace Z3i;
  typedef ImplicitDigitalEllipse3<Point> ImplicitDigitalEllipse;
  typedef LightImplicitDigitalSurface<KSpace,ImplicitDigitalEllipse> SurfaceContainer;
  typedef DigitalSurface< SurfaceContainer > Surface;
  typedef SurfaceContainer::Surfel Surfel;
 
  Point p1( -10, -10, -10 );
  Point p2( 10, 10, 10 );
  KSpace K;
  K.init( p1, p2, true );
  ImplicitDigitalEllipse ellipse( 6.0, 4.5, 3.4 );
  Surfel bel = Surfaces<KSpace>::findABel( K, ellipse, 10000 );
  SurfaceContainer* surfaceContainer = new SurfaceContainer
    ( K, ellipse, SurfelAdjacency<KSpace::dimension>( true ), bel );
  Surface surface( surfaceContainer ); // acquired
 
 
#ifdef WITH_CGAL
  typedef DGtal::functors::MongeJetFittingGaussianCurvatureEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorGaussian;
  typedef functors::MongeJetFittingMeanCurvatureEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorMean;
  typedef functors::MongeJetFittingNormalVectorEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorNormal;
  typedef functors::LinearLeastSquareFittingNormalVectorEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorNormalLeast;
 
  //constant convolution functor
  typedef functors::ConstValue<double> ConstFunctor;
 
  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorGaussian, ConstFunctor> ReporterK;
  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorMean, ConstFunctor> ReporterH;
  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorNormal, ConstFunctor> ReporterNormal;
  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>, FunctorNormalLeast, ConstFunctor> ReporterNormalLeast;
#endif
 
  typedef functors::ElementaryConvolutionNormalVectorEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorNormalElem;
  typedef LocalEstimatorFromSurfelFunctorAdapter<SurfaceContainer, LpMetric<Z3i::Space>,
                                                 FunctorNormalElem, DGtal::functors::GaussianKernel> ReporterNormalElem;
 
 
  CanonicSCellEmbedder<KSpace> embedder(surface.container().space());
#ifdef WITH_CGAL
  FunctorGaussian estimatorK(embedder,1.0);
  FunctorMean estimatorH(embedder, 1.0);
  FunctorNormal estimatorN(embedder,1.0);
  FunctorNormalLeast estimatorL(embedder,1.0);
 
  ConstFunctor constFunctor(1.0);
 
  ReporterK reporterK;
  ReporterH reporterH;
  ReporterNormal reporterN;
  ReporterNormalLeast reporterL;
#endif
 
  LpMetric<Z3i::Space> l2(2.0);
  FunctorNormalElem estimatorNormalElem(embedder,1.0);
  DGtal::functors::GaussianKernel gaussian(2.0);
  ReporterNormalElem reporterElem(surface, l2,
                                  estimatorNormalElem, gaussian);
 
#ifdef WITH_CGAL
  reporterK.attach(surface);
  reporterH.attach(surface);
  reporterN.attach(surface);
  reporterL.attach(surface);
 
  reporterK.init(1, surface.begin(), surface.end());
  reporterH.init(1, surface.begin(), surface.end());
  reporterN.init(1, surface.begin(), surface.end());
  reporterL.init(1, surface.begin(), surface.end());
 
  reporterK.setParams(l2, estimatorK, constFunctor, 5.0);
  reporterH.setParams(l2, estimatorH, constFunctor, 5.0);
  reporterN.setParams(l2, estimatorN, constFunctor, 5.0);
  reporterL.setParams(l2, estimatorL, constFunctor, 5.0);
 
  FunctorGaussian::Quantity valK = reporterK.eval( surface.begin());
  FunctorMean::Quantity valH = reporterH.eval( surface.begin());
  FunctorNormal::Quantity valN = reporterN.eval( surface.begin());
  FunctorNormalLeast::Quantity valL = reporterL.eval( surface.begin());
#endif
  
  reporterElem.attach(surface);
  reporterElem.setParams(l2,
                         estimatorNormalElem, gaussian, 5.0);
  reporterElem.init(1.0, surface.begin(), surface.end());
  FunctorNormalElem::Quantity valElem = reporterElem.eval( surface.begin());
 
#ifdef WITH_CGAL
  trace.info() << "Gaussian = "<<valK <<std::endl;
  trace.info() << "Mean = "<<valH<< std::endl;
  trace.info() << "Normal Vector (from Monge form) = "<<valN<< std::endl;
  trace.info() << "Normal Vector (linear least square) = "<<valL<< std::endl;
#endif
  trace.info() << "Normal Vector (Elementary conv) = "<<valElem<< std::endl;
 
  return 0;
}
//                                                                           //