ShortcutsGeometry.h

 
#pragma once
 
#if defined(ShortcutsGeometry_RECURSES)
#error Recursive header files inclusion detected in ShortcutsGeometry.h
#else // defined(ShortcutsGeometry_RECURSES)
#define ShortcutsGeometry_RECURSES
 
#if !defined ShortcutsGeometry_h
#define ShortcutsGeometry_h
 
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/kernel/BasicPointPredicates.h"
#include "DGtal/geometry/volumes/distance/LpMetric.h"
#include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
#include "DGtal/geometry/volumes/distance/VoronoiMap.h"
#include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
#include "DGtal/geometry/surfaces/estimation/TrueDigitalSurfaceLocalEstimator.h"
#include "DGtal/geometry/surfaces/estimation/VoronoiCovarianceMeasureOnDigitalSurface.h"
#include "DGtal/geometry/surfaces/estimation/VCMDigitalSurfaceLocalEstimator.h"
#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
#include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
#include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
#include "DGtal/geometry/meshes/CorrectedNormalCurrentComputer.h"
 
#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
#include "DGtal/dec/ATSolver2D.h"
 
namespace DGtal
{
 
  namespace sgf = ::DGtal::functors::ShapeGeometricFunctors;
 
  // template class ShortcutsGeometry
  template  < typename TKSpace >
    class ShortcutsGeometry : public Shortcuts< TKSpace >
    {
      BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
    public:
      typedef Shortcuts< TKSpace >                     Base;
      typedef ShortcutsGeometry< TKSpace >             Self;
      using Base::parametersKSpace;
      using Base::getKSpace;
      using Base::parametersDigitizedImplicitShape3D;
 
      // ----------------------- Usual space types --------------------------------------
    public:
 
      typedef TKSpace                                  KSpace;
      typedef typename KSpace::Space                   Space;
      typedef typename Space::Integer                  Integer;
      typedef typename Space::Point                    Point;
      typedef typename Space::Vector                   Vector;
      typedef typename Space::RealVector               RealVector;
      typedef typename Space::RealPoint                RealPoint;
      typedef typename RealVector::Component           Scalar;
      typedef HyperRectDomain<Space>                   Domain;
      typedef unsigned char                            GrayScale;
 
      // ----------------------- Shortcut types --------------------------------------
    public:
      typedef MPolynomial< Space::dimension, Scalar >      ScalarPolynomial;
      typedef ImplicitPolynomial3Shape<Space>              ImplicitShape3D;
      typedef GaussDigitizer< Space, ImplicitShape3D >     DigitizedImplicitShape3D;
      typedef ImageContainerBySTLVector<Domain, bool>      BinaryImage;
      typedef ImageContainerBySTLVector<Domain, GrayScale> GrayScaleImage;
      typedef ImageContainerBySTLVector<Domain, float>     FloatImage;
      typedef ImageContainerBySTLVector<Domain, double>    DoubleImage;
      typedef typename KSpace::SurfelSet                   SurfelSet;
      typedef LightImplicitDigitalSurface< KSpace, BinaryImage >  LightSurfaceContainer;
      typedef ::DGtal::DigitalSurface< LightSurfaceContainer >    LightDigitalSurface;
      typedef SetOfSurfels< KSpace, SurfelSet >                   ExplicitSurfaceContainer;
      typedef ::DGtal::DigitalSurface< ExplicitSurfaceContainer > DigitalSurface;
      typedef IndexedDigitalSurface< ExplicitSurfaceContainer >   IdxDigitalSurface;
      typedef typename LightDigitalSurface::Surfel                Surfel;
      typedef typename LightDigitalSurface::Cell                  Cell;
      typedef typename LightDigitalSurface::SCell                 SCell;
      typedef typename LightDigitalSurface::Vertex                Vertex;
      typedef typename LightDigitalSurface::Arc                   Arc;
      typedef typename LightDigitalSurface::Face                  Face;
      typedef typename LightDigitalSurface::ArcRange              ArcRange;
      typedef typename IdxDigitalSurface::Vertex                  IdxSurfel;
      typedef typename IdxDigitalSurface::Vertex                  IdxVertex;
      typedef typename IdxDigitalSurface::Arc                     IdxArc;
      typedef typename IdxDigitalSurface::ArcRange                IdxArcRange;
      typedef std::set< IdxSurfel >                               IdxSurfelSet;
      typedef std::vector< Surfel >                               SurfelRange;
      typedef std::vector< Cell >                                 CellRange;
      typedef std::vector< IdxSurfel >                            IdxSurfelRange;
      typedef std::vector< Scalar >                               Scalars;
      typedef std::vector< RealVector >                           RealVectors;
      typedef std::vector< RealPoint >                            RealPoints;
 
      typedef ::DGtal::Statistic<Scalar>                          ScalarStatistic;
 
      typedef sgf::ShapePositionFunctor<ImplicitShape3D>          PositionFunctor;
      typedef sgf::ShapeNormalVectorFunctor<ImplicitShape3D>      NormalFunctor;
      typedef sgf::ShapeMeanCurvatureFunctor<ImplicitShape3D>     MeanCurvatureFunctor;
      typedef sgf::ShapeGaussianCurvatureFunctor<ImplicitShape3D> GaussianCurvatureFunctor;
      typedef sgf::ShapeFirstPrincipalCurvatureFunctor<ImplicitShape3D> FirstPrincipalCurvatureFunctor;
      typedef sgf::ShapeSecondPrincipalCurvatureFunctor<ImplicitShape3D> SecondPrincipalCurvatureFunctor;
      typedef sgf::ShapeFirstPrincipalDirectionFunctor<ImplicitShape3D> FirstPrincipalDirectionFunctor;
      typedef sgf::ShapeSecondPrincipalDirectionFunctor<ImplicitShape3D> SecondPrincipalDirectionFunctor;
      typedef sgf::ShapePrincipalCurvaturesAndDirectionsFunctor<ImplicitShape3D> PrincipalCurvaturesAndDirectionsFunctor;
 
      typedef typename functors::IIPrincipalCurvaturesAndDirectionsFunctor<Space>::Quantity   CurvatureTensorQuantity;
      typedef std::vector< CurvatureTensorQuantity >              CurvatureTensorQuantities;
 
      typedef CorrectedNormalCurrentComputer<RealPoint, RealVector> CNCComputer;
 
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, PositionFunctor >                TruePositionEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, NormalFunctor >                  TrueNormalEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, MeanCurvatureFunctor >           TrueMeanCurvatureEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, GaussianCurvatureFunctor >       TrueGaussianCurvatureEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, FirstPrincipalCurvatureFunctor > TrueFirstPrincipalCurvatureEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, SecondPrincipalCurvatureFunctor > TrueSecondPrincipalCurvatureEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, FirstPrincipalDirectionFunctor > TrueFirstPrincipalDirectionEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, SecondPrincipalDirectionFunctor > TrueSecondPrincipalDirectionEstimator;
      typedef TrueDigitalSurfaceLocalEstimator
        < KSpace, ImplicitShape3D, PrincipalCurvaturesAndDirectionsFunctor > TruePrincipalCurvaturesAndDirectionsEstimator;
 
      typedef ::DGtal::Mesh<RealPoint>                            Mesh;
      typedef ::DGtal::TriangulatedSurface<RealPoint>             TriangulatedSurface;
      typedef ::DGtal::PolygonalSurface<RealPoint>                PolygonalSurface;
      typedef std::map<Surfel, IdxSurfel>                         Surfel2Index;
      typedef std::map<Cell,   IdxVertex>                         Cell2Index;
 
      typedef DigitalSetByAssociativeContainer<Domain, std::unordered_set<typename Domain::Point>> DigitalSet;
      typedef functors::NotPointPredicate<DigitalSet> VoronoiPointPredicate;
 
      // ----------------------- Static services --------------------------------------
    public:
 
      // ----------------------- Exact geometry services ------------------------------
    public:
 
      static Parameters defaultParameters()
      {
        return parametersShapeGeometry()
          | parametersGeometryEstimation()
          | parametersATApproximation()
          | parametersVoronoiMap();
      }
 
      static Parameters parametersShapeGeometry()
      {
        return Parameters
          ( "projectionMaxIter",  20 )
          ( "projectionAccuracy", 0.0001 )
          ( "projectionGamma",    0.5 )
          ( "gridstep",           1.0 );
      }
 
      static RealPoints
        getPositions
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        RealVectors         n_true_estimations;
        TruePositionEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, PositionFunctor(), maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static RealPoints
        getPositions
        ( CountedPtr<ImplicitShape3D> shape,
          const RealPoints&           points,
          const Parameters&           params = parametersShapeGeometry() )
      {
        RealPoints proj_points( points.size() );
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        for ( unsigned int i = 0; i < points.size(); ++i )
          proj_points[ i ] = shape->nearestPoint( points[ i ], accuracy,
                                                  maxIter, gamma );
        return proj_points;
      }
 
      static RealVectors
        getNormalVectors
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        RealVectors         n_true_estimations;
        TrueNormalEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, NormalFunctor(), maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static Scalars
        getMeanCurvatures
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        Scalars                n_true_estimations;
        TrueMeanCurvatureEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, MeanCurvatureFunctor(), maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
    static Scalars
      getCNCMeanCurvatures
      ( CountedPtr<typename Base::SurfaceMesh>  mesh, 
        const typename Base::SurfaceMesh::Faces faces, 
        const Parameters&                       params = parametersShapeGeometry() )
      {
        using Face = typename Base::SurfaceMesh::Face;
 
        bool unit_u = params["unit_u"].as<int>();
        double radius = params["r-radius"].as<double>();
        double alpha  = params["alpha"].as<double>();
        double h      = params["gridstep"].as<double>();
        if ( alpha != 1.0 ) radius *= pow( h, alpha-1.0 );
 
        CNCComputer computer(*mesh, unit_u);
        
        const auto& mu0 = computer.computeMu0();
        const auto& mu1 = computer.computeMu1();
        
        Scalars curvatures(faces.size());
        for (size_t i = 0; i < faces.size(); ++i)
        {
          const auto center = mesh->faceCentroid(faces[i]);
          const auto area = mu0.measure(center, radius, faces[i]);
          const auto lmu1 = mu1.measure(center, radius, faces[i]);
          curvatures[i] = CNCComputer::meanCurvature(area, lmu1);
        }
 
        return curvatures;
      }
 
      static Scalars getCNCMeanCurvatures(
      CountedPtr<typename Base::SurfaceMesh> mesh,
      const Parameters & params = parametersShapeGeometry() )
      {
        std::vector<typename Base::SurfaceMesh::Face> allFaces(mesh->nbFaces());
        std::iota(allFaces.begin(), allFaces.end(), 0);
 
        return getCNCMeanCurvatures(mesh, allFaces, params);
      }
 
      template <typename T>
      static Scalars getCNCMeanCurvatures(
      T & digitalObject, const Parameters & params = parametersShapeGeometry() )
      {
        CountedPtr<typename Base::SurfaceMesh> mesh = Base::makePrimalSurfaceMesh(digitalObject);
        return getCNCMeanCurvatures(mesh, params);
      }
 
      static Scalars
        getGaussianCurvatures
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        Scalars                n_true_estimations;
        TrueGaussianCurvatureEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, GaussianCurvatureFunctor(), maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static Scalars getCNCGaussianCurvatures(
      CountedPtr<typename Base::SurfaceMesh> mesh,
      const typename Base::SurfaceMesh::Faces & faces,
      const Parameters & params = parametersShapeGeometry() )
      {
        using Face = typename Base::SurfaceMesh::Face;
 
        bool unit_u = params["unit_u"].as<int>();
        double radius = params["r-radius"].as<double>();
        double alpha  = params["alpha"].as<double>();
        double h      = params["gridstep"].as<double>();
        if ( alpha != 1.0 ) radius *= pow( h, alpha-1.0 );
 
        CNCComputer computer(*mesh, unit_u);
        
        const auto& mu0 = computer.computeMu0();
        const auto& mu2 = computer.computeMu2();
 
        Scalars curvatures(faces.size());
        for (size_t i = 0; i < faces.size(); ++i)
        {
          const auto center = mesh->faceCentroid(faces[i]);
          const auto area = mu0.measure(center, radius, faces[i]);
          const auto lmu2 = mu2.measure(center, radius, faces[i]);
          curvatures[i] = CNCComputer::GaussianCurvature(area, lmu2);
        }
        
        return curvatures;
      }
 
      static Scalars getCNCGaussianCurvatures(
      CountedPtr<typename Base::SurfaceMesh> mesh,
      const Parameters & params = parametersShapeGeometry() )
      {
        std::vector<typename Base::SurfaceMesh::Face> allFaces(mesh->nbFaces());
        std::iota(allFaces.begin(), allFaces.end(), 0);
        
        return getCNCGaussianCurvatures(mesh, allFaces, params);
      }
 
      template <typename T>
      static Scalars getCNCGaussianCurvatures(
      T & digitalObject, const Parameters & params = parametersShapeGeometry() )
      {
        CountedPtr<typename Base::SurfaceMesh> mesh = Base::makePrimalSurfaceMesh(digitalObject);
        return getCNCGaussianCurvatures(mesh, params);
      }
 
      static Scalars
      getFirstPrincipalCurvatures
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        Scalars n_true_estimations;
        TrueFirstPrincipalCurvatureEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, FirstPrincipalCurvatureFunctor(),
                                  maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static Scalars
      getSecondPrincipalCurvatures
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        Scalars n_true_estimations;
        TrueSecondPrincipalCurvatureEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, SecondPrincipalCurvatureFunctor(),
                                  maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static RealVectors
      getFirstPrincipalDirections
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        RealVectors n_true_estimations;
        TrueFirstPrincipalDirectionEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, FirstPrincipalDirectionFunctor(),
                                  maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static RealVectors
      getSecondPrincipalDirections
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        RealVectors n_true_estimations;
        TrueSecondPrincipalDirectionEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, SecondPrincipalDirectionFunctor(),
                                  maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static CurvatureTensorQuantities
      getPrincipalCurvaturesAndDirections
        ( CountedPtr<ImplicitShape3D> shape,
          const KSpace&               K,
          const SurfelRange&          surfels,
          const Parameters&           params = parametersShapeGeometry() )
      {
        CurvatureTensorQuantities n_true_estimations;
        TruePrincipalCurvaturesAndDirectionsEstimator true_estimator;
        int     maxIter = params[ "projectionMaxIter"  ].as<int>();
        double accuracy = params[ "projectionAccuracy" ].as<double>();
        double    gamma = params[ "projectionGamma"    ].as<double>();
        Scalar gridstep = params[ "gridstep"           ].as<Scalar>();
        true_estimator.attach( *shape );
        true_estimator.setParams( K, PrincipalCurvaturesAndDirectionsFunctor(),
                                  maxIter, accuracy, gamma );
        true_estimator.init( gridstep, surfels.begin(), surfels.end() );
        true_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_true_estimations ) );
        return n_true_estimations;
      }
 
      static std::tuple<Scalars, Scalars, RealVectors, RealVectors>
        getCNCPrincipalCurvaturesAndDirections
        ( CountedPtr<typename Base::SurfaceMesh>   mesh, 
          const typename Base::SurfaceMesh::Faces& faces, 
          const Parameters&                        params = parametersShapeGeometry() )
        {
          using Face = typename Base::SurfaceMesh::Face;
 
          bool unit_u = params["unit_u"].as<int>();
          double radius = params["r-radius"].as<double>();
          double alpha  = params["alpha"].as<double>();
          double h      = params["gridstep"].as<double>();
          if ( alpha != 1.0 ) radius *= pow( h, alpha-1.0 );
 
          CNCComputer computer(*mesh, unit_u);
          
          const auto& mu0  = computer.computeMu0();
          const auto& muxy = computer.computeMuXY();
          
          if (mesh->faceNormals().size() == 0)
          {
            // Try to use vertex normals if any
            if (mesh->vertexNormals().size() == 0)
              mesh->computeFaceNormalsFromPositions();
            else 
              mesh->computeFaceNormalsFromVertexNormals();
          }
          
          const auto& normals = mesh->faceNormals();
 
          Scalars k1(faces.size()), k2(faces.size());
          RealVectors d1(faces.size()), d2(faces.size());
 
          for (size_t i = 0; i < faces.size(); ++i)
          {
            const auto center = mesh->faceCentroid(faces[i]);
            const auto area  = mu0 .measure(center, radius, faces[i]);
            const auto lmuxy = muxy.measure(center, radius, faces[i]);
            std::tie(k1[i], k2[i], d1[i], d2[i]) = 
                CNCComputer::principalCurvatures(area, lmuxy, normals[faces[i]]);
          }
 
          return std::make_tuple(k1, k2, d1, d2);
        }
 
        static std::tuple<Scalars, Scalars, RealVectors, RealVectors>
        getCNCPrincipalCurvaturesAndDirections(
        CountedPtr<typename Base::SurfaceMesh> mesh,
        const Parameters & params = parametersShapeGeometry() )
        {
          std::vector<typename Base::SurfaceMesh::Face> allFaces(mesh->nbFaces());
          std::iota(allFaces.begin(), allFaces.end(), 0);
 
          return getCNCPrincipalCurvaturesAndDirections(mesh, allFaces, params);
        }
 
        template <typename T>
        static std::tuple<Scalars, Scalars, RealVectors, RealVectors>
        getCNCPrincipalCurvaturesAndDirections(
        T & digitalObject,
        const Parameters & params = parametersShapeGeometry() )
        {
          CountedPtr<typename Base::SurfaceMesh> mesh = Base::makePrimalSurfaceMesh(digitalObject);
          return getCNCPrincipalCurvaturesAndDirections(mesh, params);
        }
 
      // --------------------------- geometry estimation ------------------------------
    public:
 
      static Parameters parametersGeometryEstimation()
      {
        return Parameters
          ( "verbose",           1 )
          ( "t-ring",          3.0 )
          ( "kernel",        "hat" )
          ( "R-radius",       10.0 )
          ( "r-radius",        3.0 )
          ( "alpha",          0.33 )
          ( "surfelEmbedding",   0 )
          ( "unit_u"         ,   0 );
      }
 
      static RealVectors
        getTrivialNormalVectors( const KSpace&      K,
                                 const SurfelRange& surfels )
      {
        std::vector< RealVector > result;
        for ( auto s : surfels )
          {
            Dimension  k = K.sOrthDir( s );
            bool  direct = K.sDirect( s, k );
            RealVector t = RealVector::zero;
            t[ k ]       = direct ? -1.0 : 1.0;
            result.push_back( t );
          }
        return result;
      }
 
      template <typename TAnyDigitalSurface>
        static RealVectors
        getCTrivialNormalVectors
        ( CountedPtr<TAnyDigitalSurface> surface,
          const SurfelRange&             surfels,
          const Parameters&              params = parametersGeometryEstimation() )
        {
          int    verbose = params[ "verbose"  ].as<int>();
          Scalar       t = params[ "t-ring"   ].as<double>();
          typedef typename TAnyDigitalSurface::DigitalSurfaceContainer  SurfaceContainer;
          typedef LpMetric<Space>                                       Metric;
          typedef functors::HatFunction<Scalar>                         Functor;
          typedef functors::ElementaryConvolutionNormalVectorEstimator
            < Surfel, CanonicSCellEmbedder<KSpace> >                    SurfelFunctor;
          typedef LocalEstimatorFromSurfelFunctorAdapter
            < SurfaceContainer, Metric, SurfelFunctor, Functor>         NormalEstimator;
          if ( verbose > 0 )
            trace.info() << "- CTrivial normal t-ring=" << t << " (discrete)" << std::endl;
          const Functor fct( 1.0, t );
          const KSpace &  K = surface->container().space();
          Metric    aMetric( 2.0 );
          CanonicSCellEmbedder<KSpace> canonic_embedder( K );
          std::vector< RealVector >    n_estimations;
          SurfelFunctor                surfelFct( canonic_embedder, 1.0 );
          NormalEstimator              estimator;
          estimator.attach( *surface);
          estimator.setParams( aMetric, surfelFct, fct, t );
          estimator.init( 1.0, surfels.begin(), surfels.end());
          estimator.eval( surfels.begin(), surfels.end(),
                          std::back_inserter( n_estimations ) );
          std::transform( n_estimations.cbegin(), n_estimations.cend(), n_estimations.begin(),
                          [] ( RealVector v ) { return -v; } );
          return n_estimations;
        }
 
      template <typename TAnyDigitalSurface>
        static RealVectors
        getVCMNormalVectors
        ( CountedPtr<TAnyDigitalSurface> surface,
          const SurfelRange&             surfels,
          const Parameters&              params = parametersGeometryEstimation() )
        {
          typedef ExactPredicateLpSeparableMetric<Space,2> Metric;
          typedef typename TAnyDigitalSurface::DigitalSurfaceContainer SurfaceContainer;
          RealVectors n_estimations;
          int        verbose = params[ "verbose"   ].as<int>();
          std::string kernel = params[ "kernel"    ].as<std::string>();
          Scalar      h      = params[ "gridstep"  ].as<Scalar>();
          Scalar      R      = params[ "R-radius"  ].as<Scalar>();
          Scalar      r      = params[ "r-radius"  ].as<Scalar>();
          Scalar      t      = params[ "t-ring"    ].as<Scalar>();
          Scalar      alpha  = params[ "alpha"     ].as<Scalar>();
          int      embedding = params[ "embedding" ].as<int>();
          // Adjust parameters according to gridstep if specified.
          if ( alpha != 1.0 ) R *= pow( h, alpha-1.0 );
          if ( alpha != 1.0 ) r *= pow( h, alpha-1.0 );
          Surfel2PointEmbedding embType = embedding == 0 ? Pointels :
            embedding == 1 ? InnerSpel : OuterSpel;
          if ( verbose > 0 )
            {
              trace.info() << "- VCM normal kernel=" << kernel << " emb=" << embedding
                           << " alpha=" << alpha << std::endl;
              trace.info() << "- VCM normal r=" << (r*h)  << " (continuous) "
                           << r << " (discrete)" << std::endl;
              trace.info() << "- VCM normal R=" << (R*h)  << " (continuous) "
                           << R << " (discrete)" << std::endl;
              trace.info() << "- VCM normal t=" << t << " (discrete)" << std::endl;
            }
          if ( kernel == "hat" )
            {
              typedef functors::HatPointFunction<Point,Scalar>             KernelFunction;
              typedef VoronoiCovarianceMeasureOnDigitalSurface
                < SurfaceContainer, Metric, KernelFunction >               VCMOnSurface;
              typedef functors::VCMNormalVectorFunctor<VCMOnSurface>       NormalVFunctor;
              typedef VCMDigitalSurfaceLocalEstimator
                < SurfaceContainer, Metric, KernelFunction, NormalVFunctor> VCMNormalEstimator;
              KernelFunction chi_r( 1.0, r );
              VCMNormalEstimator estimator;
              estimator.attach( *surface );
              estimator.setParams( embType, R, r, chi_r, t, Metric(), verbose > 0 );
              estimator.init( h, surfels.begin(), surfels.end() );
              estimator.eval( surfels.begin(), surfels.end(),
                              std::back_inserter( n_estimations ) );
            }
          else if ( kernel == "ball" )
            {
              typedef functors::BallConstantPointFunction<Point,Scalar>    KernelFunction;
              typedef VoronoiCovarianceMeasureOnDigitalSurface
                < SurfaceContainer, Metric, KernelFunction >               VCMOnSurface;
              typedef functors::VCMNormalVectorFunctor<VCMOnSurface>       NormalVFunctor;
              typedef VCMDigitalSurfaceLocalEstimator
                < SurfaceContainer, Metric, KernelFunction, NormalVFunctor> VCMNormalEstimator;
              KernelFunction chi_r( 1.0, r );
              VCMNormalEstimator estimator;
              estimator.attach( *surface );
              estimator.setParams( embType, R, r, chi_r, t, Metric(), verbose > 0 );
              estimator.init( h, surfels.begin(), surfels.end() );
              estimator.eval( surfels.begin(), surfels.end(),
                              std::back_inserter( n_estimations ) );
            }
          else
            {
              trace.warning() << "[ShortcutsGeometry::getVCMNormalVectors] Unknown kernel: "
                              << kernel << std::endl;
            }
          return n_estimations;
        }
 
      static RealVectors
        getIINormalVectors( CountedPtr<BinaryImage> bimage,
                            const SurfelRange&      surfels,
                            const Parameters&       params
                            = parametersGeometryEstimation()
                            | parametersKSpace() )
      {
        auto K =  getKSpace( bimage, params );
        return getIINormalVectors( *bimage, K, surfels, params );
      }
 
      static RealVectors
        getIINormalVectors( CountedPtr< DigitizedImplicitShape3D > dshape,
                            const SurfelRange&      surfels,
                            const Parameters&       params
                            = parametersGeometryEstimation()
                            | parametersKSpace()
                            | parametersDigitizedImplicitShape3D() )
      {
        auto K =  getKSpace( params );
        return getIINormalVectors( *dshape, K, surfels, params );
      }
 
      template <typename TPointPredicate>
        static RealVectors
        getIINormalVectors( const TPointPredicate&  shape,
                            const KSpace&           K,
                            const SurfelRange&      surfels,
                            const Parameters&       params
                            = parametersGeometryEstimation()
                            | parametersKSpace() )
        {
          typedef functors::IINormalDirectionFunctor<Space> IINormalFunctor;
          typedef IntegralInvariantCovarianceEstimator
            <KSpace, TPointPredicate, IINormalFunctor>          IINormalEstimator;
 
          RealVectors n_estimations;
          int        verbose = params[ "verbose"   ].as<int>();
          Scalar     h       = params[ "gridstep"  ].as<Scalar>();
          Scalar     r       = params[ "r-radius"  ].as<Scalar>();
          Scalar     alpha   = params[ "alpha"     ].as<Scalar>();
          if ( alpha != 1.0 ) r *= pow( h, alpha-1.0 );
          if ( verbose > 0 )
            {
              trace.info() << "- II normal alpha=" << alpha << std::endl;
              trace.info() << "- II normal r=" << (r*h)  << " (continuous) "
                           << r << " (discrete)" << std::endl;
            }
          IINormalFunctor     functor;
          functor.init( h, r*h );
          IINormalEstimator   ii_estimator( functor );
          ii_estimator.attach( K, shape );
          ii_estimator.setParams( r );
          ii_estimator.init( h, surfels.begin(), surfels.end() );
          ii_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( n_estimations ) );
          const RealVectors n_trivial = getTrivialNormalVectors( K, surfels );
          orientVectors( n_estimations, n_trivial );
          return n_estimations;
        }
 
 
      static Scalars
        getIIMeanCurvatures( CountedPtr<BinaryImage> bimage,
                             const SurfelRange&      surfels,
                             const Parameters&       params
                             = parametersGeometryEstimation()
                             | parametersKSpace() )
      {
        auto K =  getKSpace( bimage, params );
        return getIIMeanCurvatures( *bimage, K, surfels, params );
      }
 
      static Scalars
        getIIMeanCurvatures( CountedPtr< DigitizedImplicitShape3D > dshape,
                             const SurfelRange&      surfels,
                             const Parameters&       params
                             = parametersGeometryEstimation()
                             | parametersKSpace()
                             | parametersDigitizedImplicitShape3D() )
      {
        auto K =  getKSpace( params );
        return getIIMeanCurvatures( *dshape, K, surfels, params );
      }
 
 
      template <typename TPointPredicate>
        static Scalars
        getIIMeanCurvatures( const TPointPredicate&  shape,
                             const KSpace&           K,
                             const SurfelRange&      surfels,
                             const Parameters&       params
                             = parametersGeometryEstimation()
                             | parametersKSpace() )
        {
          typedef functors::IIMeanCurvature3DFunctor<Space> IIMeanCurvFunctor;
          typedef IntegralInvariantVolumeEstimator
            <KSpace, TPointPredicate, IIMeanCurvFunctor>    IIMeanCurvEstimator;
 
          Scalars  mc_estimations;
          int      verbose = params[ "verbose"   ].as<int>();
          Scalar   h       = params[ "gridstep"  ].as<Scalar>();
          Scalar   r       = params[ "r-radius"  ].as<Scalar>();
          Scalar   alpha   = params[ "alpha"     ].as<Scalar>();
          if ( alpha != 1.0 ) r *= pow( h, alpha-1.0 );
          if ( verbose > 0 )
            {
              trace.info() << "- II mean curvature alpha=" << alpha << std::endl;
              trace.info() << "- II mean curvature r=" << (r*h)  << " (continuous) "
                           << r << " (discrete)" << std::endl;
            }
          IIMeanCurvFunctor   functor;
          functor.init( h, r*h );
          IIMeanCurvEstimator ii_estimator( functor );
          ii_estimator.attach( K, shape );
          ii_estimator.setParams( r );
          ii_estimator.init( h, surfels.begin(), surfels.end() );
          ii_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( mc_estimations ) );
          return mc_estimations;
        }
 
      static Scalars
        getIIGaussianCurvatures( CountedPtr<BinaryImage> bimage,
                                 const SurfelRange&      surfels,
                                 const Parameters&       params
                                 = parametersGeometryEstimation()
                                 | parametersKSpace() )
      {
        auto K =  getKSpace( bimage, params );
        return getIIGaussianCurvatures( *bimage, K, surfels, params );
      }
 
      static Scalars
        getIIGaussianCurvatures( CountedPtr< DigitizedImplicitShape3D > dshape,
                                 const SurfelRange&      surfels,
                                 const Parameters&       params
                                 = parametersGeometryEstimation()
                                 | parametersKSpace()
                                 | parametersDigitizedImplicitShape3D() )
      {
        auto K =  getKSpace( params );
        return getIIGaussianCurvatures( *dshape, K, surfels, params );
      }
 
 
      template <typename TPointPredicate>
        static Scalars
        getIIGaussianCurvatures( const TPointPredicate&  shape,
                                 const KSpace&           K,
                                 const SurfelRange&      surfels,
                                 const Parameters&       params
                                 = parametersGeometryEstimation()
                                 | parametersKSpace() )
        {
          typedef functors::IIGaussianCurvature3DFunctor<Space> IIGaussianCurvFunctor;
          typedef IntegralInvariantCovarianceEstimator
            <KSpace, TPointPredicate, IIGaussianCurvFunctor>    IIGaussianCurvEstimator;
 
          Scalars  mc_estimations;
          int      verbose = params[ "verbose"   ].as<int>();
          Scalar   h       = params[ "gridstep"  ].as<Scalar>();
          Scalar   r       = params[ "r-radius"  ].as<Scalar>();
          Scalar   alpha   = params[ "alpha"     ].as<Scalar>();
          if ( alpha != 1.0 ) r *= pow( h, alpha-1.0 );
          if ( verbose > 0 )
            {
              trace.info() << "- II Gaussian curvature alpha=" << alpha << std::endl;
              trace.info() << "- II Gaussian curvature r=" << (r*h)  << " (continuous) "
                           << r << " (discrete)" << std::endl;
            }
          IIGaussianCurvFunctor   functor;
          functor.init( h, r*h );
          IIGaussianCurvEstimator ii_estimator( functor );
          ii_estimator.attach( K, shape );
          ii_estimator.setParams( r );
          ii_estimator.init( h, surfels.begin(), surfels.end() );
          ii_estimator.eval( surfels.begin(), surfels.end(),
                             std::back_inserter( mc_estimations ) );
          return mc_estimations;
        }
 
      static CurvatureTensorQuantities
      getIIPrincipalCurvaturesAndDirections( CountedPtr<BinaryImage> bimage,
                                          const SurfelRange&      surfels,
                                          const Parameters&       params
                                            = parametersGeometryEstimation()
                                            | parametersKSpace() )
      {
        auto K =  getKSpace( bimage, params );
        return getIIPrincipalCurvaturesAndDirections( *bimage, K, surfels, params );
      }
 
      static CurvatureTensorQuantities
      getIIPrincipalCurvaturesAndDirections( CountedPtr< DigitizedImplicitShape3D > dshape,
                                          const SurfelRange&      surfels,
                                          const Parameters&       params
                                            = parametersGeometryEstimation()
                                            | parametersKSpace()
                                            | parametersDigitizedImplicitShape3D() )
      {
        auto K =  getKSpace( params );
        return getIIPrincipalCurvaturesAndDirections( *dshape, K, surfels, params );
      }
 
 
      template <typename TPointPredicate>
      static CurvatureTensorQuantities
      getIIPrincipalCurvaturesAndDirections( const TPointPredicate&  shape,
                                          const KSpace&           K,
                                          const SurfelRange&      surfels,
                                          const Parameters&       params
                                            = parametersGeometryEstimation()
                                            | parametersKSpace() )
      {
        typedef functors::IIPrincipalCurvaturesAndDirectionsFunctor<Space> IICurvFunctor;
        typedef IntegralInvariantCovarianceEstimator<KSpace, TPointPredicate, IICurvFunctor>    IICurvEstimator;
 
        CurvatureTensorQuantities  mc_estimations;
        int      verbose = params[ "verbose"   ].as<int>();
        Scalar   h       = params[ "gridstep"  ].as<Scalar>();
        Scalar   r       = params[ "r-radius"  ].as<Scalar>();
        Scalar   alpha   = params[ "alpha"     ].as<Scalar>();
        if ( alpha != 1.0 ) r *= pow( h, alpha-1.0 );
        if ( verbose > 0 )
        {
          trace.info() << "- II principal curvatures and directions alpha=" << alpha << std::endl;
          trace.info() << "- II principal curvatures and directions r=" << (r*h)  << " (continuous) "
          << r << " (discrete)" << std::endl;
        }
        IICurvFunctor   functor;
        functor.init( h, r*h );
        IICurvEstimator ii_estimator( functor );
        ii_estimator.attach( K, shape );
        ii_estimator.setParams( r );
        ii_estimator.init( h, surfels.begin(), surfels.end() );
        ii_estimator.eval( surfels.begin(), surfels.end(),
                          std::back_inserter( mc_estimations ) );
        return mc_estimations;
      }
 
 
      // --------------------------- AT approximation ------------------------------
    public:
 
      static Parameters parametersATApproximation()
      {
        return Parameters
          ( "at-enabled",        1 )
          ( "at-alpha",          0.1 )
          ( "at-lambda",         0.025 )
          ( "at-epsilon",        0.25 )
          ( "at-epsilon-start",  2.0 )
          ( "at-epsilon-ratio",  2.0 )
          ( "at-max-iter",      10 )
          ( "at-diff-v-max",     0.0001 )
          ( "at-v-policy",   "Maximum" );
      }
 
      template <typename TAnyDigitalSurface,
                typename VectorFieldInput>
      static
      VectorFieldInput
      getATVectorFieldApproximation( CountedPtr<TAnyDigitalSurface> surface,
                                     const SurfelRange&             surfels,
                                     const VectorFieldInput&        input,
                                     const Parameters&              params
                                     = parametersATApproximation() | parametersGeometryEstimation() )
      {
        (void)surface; //param not used. FIXME: JOL
 
        int      verbose   = params[ "verbose"          ].as<int>();
        Scalar   alpha_at  = params[ "at-alpha"         ].as<Scalar>();
        Scalar   lambda_at = params[ "at-lambda"        ].as<Scalar>();
        Scalar   epsilon1  = params[ "at-epsilon-start" ].as<Scalar>();
        Scalar   epsilon2  = params[ "at-epsilon"       ].as<Scalar>();
        Scalar   epsilonr  = params[ "at-epsilon-ratio" ].as<Scalar>();
        int      max_iter  = params[ "at-max-iter"      ].as<int>();
        Scalar   diff_v_max= params[ "at-diff-v-max"    ].as<Scalar>();
        typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;
        const auto calculus = CalculusFactory::createFromNSCells<2>( surfels.cbegin(), surfels.cend() );
        ATSolver2D< KSpace > at_solver( calculus, verbose );
        at_solver.initInputVectorFieldU2( input, surfels.cbegin(), surfels.cend() );
        at_solver.setUp( alpha_at, lambda_at );
        at_solver.solveGammaConvergence( epsilon1, epsilon2, epsilonr, false, diff_v_max, max_iter );
        auto output = input;
        at_solver.getOutputVectorFieldU2( output, surfels.cbegin(), surfels.cend() );
        return output;
      }
 
      template <typename TAnyDigitalSurface,
                typename VectorFieldInput,
                typename CellRangeConstIterator>
      static
      VectorFieldInput
      getATVectorFieldApproximation( Scalars&                       features,
                                     CellRangeConstIterator         itB,
                                     CellRangeConstIterator         itE,
                                     CountedPtr<TAnyDigitalSurface> surface,
                                     const SurfelRange&             surfels,
                                     const VectorFieldInput&        input,
                                     const Parameters&              params
                                     = parametersATApproximation() | parametersGeometryEstimation() )
      {
        (void)surface; //param not used FIXME: JOL
        
        int      verbose   = params[ "verbose"          ].as<int>();
        Scalar   alpha_at  = params[ "at-alpha"         ].as<Scalar>();
        Scalar   lambda_at = params[ "at-lambda"        ].as<Scalar>();
        Scalar   epsilon1  = params[ "at-epsilon-start" ].as<Scalar>();
        Scalar   epsilon2  = params[ "at-epsilon"       ].as<Scalar>();
        Scalar   epsilonr  = params[ "at-epsilon-ratio" ].as<Scalar>();
        int      max_iter  = params[ "at-max-iter"      ].as<int>();
        Scalar   diff_v_max= params[ "at-diff-v-max"    ].as<Scalar>();
        std::string policy = params[ "at-v-policy"      ].as<std::string>();
        typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;
        const auto calculus = CalculusFactory::createFromNSCells<2>( surfels.cbegin(), surfels.cend() );
        ATSolver2D< KSpace > at_solver( calculus, verbose );
        at_solver.initInputVectorFieldU2( input, surfels.cbegin(), surfels.cend() );
        at_solver.setUp( alpha_at, lambda_at );
        at_solver.solveGammaConvergence( epsilon1, epsilon2, epsilonr, false, diff_v_max, max_iter );
        auto output = input;
        at_solver.getOutputVectorFieldU2( output, surfels.cbegin(), surfels.cend() );
        auto p = ( policy == "Average" ) ? at_solver.Average
          :      ( policy == "Minimum" ) ? at_solver.Minimum
          :                                at_solver.Maximum;
        at_solver.getOutputScalarFieldV0( features, itB, itE, p );
        return output;
      }
 
      template <typename TAnyDigitalSurface>
      static
      Scalars
      getATScalarFieldApproximation( CountedPtr<TAnyDigitalSurface> surface,
                                     const SurfelRange&             surfels,
                                     const Scalars&                 input,
                                     const Parameters&              params
                                     = parametersATApproximation() | parametersGeometryEstimation() )
      {
        (void)surface; //param not used FIXME: JOL
 
        int      verbose   = params[ "verbose"          ].as<int>();
        Scalar   alpha_at  = params[ "at-alpha"         ].as<Scalar>();
        Scalar   lambda_at = params[ "at-lambda"        ].as<Scalar>();
        Scalar   epsilon1  = params[ "at-epsilon-start" ].as<Scalar>();
        Scalar   epsilon2  = params[ "at-epsilon"       ].as<Scalar>();
        Scalar   epsilonr  = params[ "at-epsilon-ratio" ].as<Scalar>();
        int      max_iter  = params[ "at-max-iter"      ].as<int>();
        Scalar   diff_v_max= params[ "at-diff-v-max"    ].as<Scalar>();
        typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;
        const auto calculus = CalculusFactory::createFromNSCells<2>( surfels.cbegin(), surfels.cend() );
        ATSolver2D< KSpace > at_solver( calculus, verbose );
        at_solver.initInputScalarFieldU2( input, surfels.cbegin(), surfels.cend() );
        at_solver.setUp( alpha_at, lambda_at );
        at_solver.solveGammaConvergence( epsilon1, epsilon2, epsilonr, false, diff_v_max, max_iter );
        auto output = input;
        at_solver.getOutputScalarFieldU2( output, surfels.cbegin(), surfels.cend() );
        return output;
      }
 
      template <typename TAnyDigitalSurface,
                typename CellRangeConstIterator>
      static
      Scalars
      getATScalarFieldApproximation( Scalars&                       features,
                                     CellRangeConstIterator         itB,
                                     CellRangeConstIterator         itE,
                                     CountedPtr<TAnyDigitalSurface> surface,
                                     const SurfelRange&             surfels,
                                     const Scalars&                 input,
                                     const Parameters&              params
                                     = parametersATApproximation() | parametersGeometryEstimation() )
      {
        (void)surface; //param not used FIXME: JOL
        
        int      verbose   = params[ "verbose"          ].as<int>();
        Scalar   alpha_at  = params[ "at-alpha"         ].as<Scalar>();
        Scalar   lambda_at = params[ "at-lambda"        ].as<Scalar>();
        Scalar   epsilon1  = params[ "at-epsilon-start" ].as<Scalar>();
        Scalar   epsilon2  = params[ "at-epsilon"       ].as<Scalar>();
        Scalar   epsilonr  = params[ "at-epsilon-ratio" ].as<Scalar>();
        int      max_iter  = params[ "at-max-iter"      ].as<int>();
        Scalar   diff_v_max= params[ "at-diff-v-max"    ].as<Scalar>();
        std::string policy = params[ "at-v-policy"      ].as<std::string>();
        typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;
        const auto calculus = CalculusFactory::createFromNSCells<2>( surfels.cbegin(), surfels.cend() );
        ATSolver2D< KSpace > at_solver( calculus, verbose );
        at_solver.initInputScalarFieldU2( input, surfels.cbegin(), surfels.cend() );
        at_solver.setUp( alpha_at, lambda_at );
        at_solver.solveGammaConvergence( epsilon1, epsilon2, epsilonr, false, diff_v_max, max_iter );
        auto output = input;
        at_solver.getOutputScalarFieldU2( output, surfels.cbegin(), surfels.cend() );
        auto p = ( policy == "Average" ) ? at_solver.Average
          :      ( policy == "Minimum" ) ? at_solver.Minimum
          :                                at_solver.Maximum;
        at_solver.getOutputScalarFieldV0( features, itB, itE, p );
        return output;
      }
 
 
      // ------------------------- Error measures services -------------------------
    public:
 
      static void
        orientVectors( RealVectors&       v,
                       const RealVectors& ref_v )
      {
        std::transform( ref_v.cbegin(), ref_v.cend(), v.cbegin(), v.begin(),
                        [] ( RealVector rw, RealVector w )
                        { return rw.dot( w ) >= 0.0 ? w : -w; } );
      }
 
      static ScalarStatistic
        getStatistic( const Scalars& v )
      {
        ScalarStatistic stat;
        stat.addValues( v.begin(), v.end() );
        stat.terminate();
        return stat;
      }
 
      static Scalars
        getVectorsAngleDeviation( const RealVectors& v1,
                                  const RealVectors& v2 )
      {
        Scalars v( v1.size() );
        if ( v1.size() == v2.size() )
          {
            auto outIt = v.begin();
            for ( auto it1 = v1.cbegin(), it2 = v2.cbegin(), itE1 = v1.cend();
                  it1 != itE1; ++it1, ++it2 )
              {
                Scalar angle_error = acos( (*it1).dot( *it2 ) );
                *outIt++ = angle_error;
              }
          }
        else
          {
            trace.warning() << "[ShortcutsGeometry::getVectorsAngleDeviation]"
                            << " v1.size()=" << v1.size() << " should be equal to "
                            << " v2.size()=" << v2.size() << std::endl;
          }
        return v;
      }
 
      static Scalars
        getScalarsAbsoluteDifference( const Scalars & v1,
                                      const Scalars & v2 )
      {
        Scalars result( v1.size() );
        std::transform( v2.cbegin(), v2.cend(), v1.cbegin(), result.begin(),
                        [] ( Scalar val1, Scalar val2 )
                        { return fabs( val1 - val2 ); } );
        return result;
      }
 
      static Scalar
        getScalarsNormL2( const Scalars & v1,
                          const Scalars & v2 )
      {
        Scalar sum = 0;
        for ( unsigned int i = 0; i < v1.size(); i++ )
          sum += ( v1[ i ] - v2[ i ] ) * ( v1[ i ] - v2[ i ] );
        return sqrt( sum / v1.size() );
      }
 
      static Scalar
        getScalarsNormL1( const Scalars & v1,
                          const Scalars & v2 )
      {
        Scalar sum = 0;
        for ( unsigned int i = 0; i < v1.size(); i++ )
          sum += fabs( v1[ i ] - v2[ i ] );
        return sum / v1.size();
      }
 
      static Scalar
        getScalarsNormLoo( const Scalars & v1,
                           const Scalars & v2 )
      {
        Scalar loo = 0;
        for ( unsigned int i = 0; i < v1.size(); i++ )
          loo = std::max( loo, fabs( v1[ i ] - v2[ i ] ) );
        return loo;
      }
      
      // ----------------------- VoronoiMap services ------------------------------
    public:
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      static Parameters parametersVoronoiMap() {
        return Parameters
          // Toricity might be moved elsewhere as this is quite a general parameter
          ( "toroidal-x" , false )
          ( "toroidal-y" , false )
          ( "toroidal-z" , false );
      }
 
 
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      template<uint32_t p, typename PointRange>
      static VoronoiMap<Space, VoronoiPointPredicate, ExactPredicateLpSeparableMetric<Space, p>>
        getVoronoiMap(Domain domain, 
                      const PointRange& sites,
                      const Parameters& params = parametersVoronoiMap())
      {
        using Metric = ExactPredicateLpSeparableMetric<Space, p>;
        using Map = VoronoiMap<Space, VoronoiPointPredicate, Metric>;
        DigitalSet set(domain); set.insert(sites.begin(), sites.end());
        VoronoiPointPredicate predicate(set);
        Metric metric;
 
        typename Map::PeriodicitySpec specs = {false, false, false};
        if (params["toroidal-x"].as<int>()) specs[0] = true;
        if (params["toroidal-y"].as<int>()) specs[1] = true;
        if (params["toroidal-z"].as<int>()) specs[2] = true;
 
 
        // Do not return a pointer here for two reasons:
        //  - The distance transform will not be passed anywhere else
        //  - The operator() is less accessible with pointers.
        return Map(domain, predicate, metric, specs);
      }
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      template<uint32_t p, typename PointRange>
      static VoronoiMap<Space, VoronoiPointPredicate, ExactPredicateLpSeparableMetric<Space, p>>
        getVoronoiMap(CountedPtr<Domain> domain, 
                      const PointRange& sites,
                      const Parameters& params = parametersVoronoiMap())
      {
        using Metric = ExactPredicateLpSeparableMetric<Space, p>;
        using Map = VoronoiMap<Space, VoronoiPointPredicate, Metric>;
        DigitalSet set(*domain); set.insert(sites.begin(), sites.end());
        VoronoiPointPredicate predicate(set);
        Metric metric;
 
        typename Map::PeriodicitySpec specs = {false, false, false};
        if (params["toroidal-x"].as<int>()) specs[0] = true;
        if (params["toroidal-y"].as<int>()) specs[1] = true;
        if (params["toroidal-z"].as<int>()) specs[2] = true;
 
        // Do not return a pointer here for two reasons:
        //  - The distance transform will not be passed anywhere else
        //  - The operator() is less accessible with pointers.
        return Map(*domain, predicate, metric, specs);
      }
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      template<uint32_t p, typename PointRange>
      static DistanceTransformation<Space, VoronoiPointPredicate, ExactPredicateLpSeparableMetric<Space, p>>
        getDistanceTransformation(Domain domain, 
                                   const PointRange& sites,
                                   const Parameters& params = parametersVoronoiMap())
      {
        using Metric = ExactPredicateLpSeparableMetric<Space, p>;
        using Map = DistanceTransformation<Space, VoronoiPointPredicate, Metric>;
        DigitalSet set(domain); set.insert(sites.begin(), sites.end());
        VoronoiPointPredicate predicate(set);
        Metric metric;
 
        typename Map::PeriodicitySpec specs = {false, false, false};
        if (params["toroidal-x"].as<int>()) specs[0] = true;
        if (params["toroidal-y"].as<int>()) specs[1] = true;
        if (params["toroidal-z"].as<int>()) specs[2] = true;
 
        // Do not return a pointer here for two reasons:
        //  - The distance transform will not be passed anywhere else
        //  - The operator() is less accessible with pointers.
        return Map(domain, predicate, metric, specs);
      }
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      template<uint32_t p, typename PointRange>
      static DistanceTransformation<Space, VoronoiPointPredicate, ExactPredicateLpSeparableMetric<Space, p>>
        getDistanceTransformation(CountedPtr<Domain> domain, 
                              const PointRange& sites,
                              const Parameters& params = parametersVoronoiMap())
      {
        using Metric = ExactPredicateLpSeparableMetric<Space, p>;
        using Map = DistanceTransformation<Space, VoronoiPointPredicate, Metric>;
        DigitalSet set(*domain); set.insert(sites.begin(), sites.end());
        VoronoiPointPredicate predicate(set);
        Metric metric;
 
        typename Map::PeriodicitySpec specs = {false, false, false};
        if (params["toroidal-x"].as<int>()) specs[0] = true;
        if (params["toroidal-y"].as<int>()) specs[1] = true;
        if (params["toroidal-z"].as<int>()) specs[2] = true;
        
        // Do not return a pointer here for two reasons:
        //  - The distance transform will not be passed anywhere else
        //  - The operator() is less accessible with pointers.
        return Map(*domain, predicate, metric, specs);
      }
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      template<uint32_t p, typename PointRangeSites, typename PointRange>
      static std::vector<Vector> getDirectionToClosestSite(
        const PointRange& points, 
        const PointRangeSites& sites,
        const Parameters& params = parametersVoronoiMap())
      {
        using Metric = ExactPredicateLpSeparableMetric<Space, p>;
        using Map = VoronoiMap<Space, VoronoiPointPredicate, Metric>;
 
        // Compute domain of points
        Point pmin = *points.begin();
        Point pmax = pmin;
 
        size_t pCount = 0;
        for (auto it = points.begin(); it != points.end(); ++it) 
        {
          pCount ++;
          for (size_t i = 0; i < Space::dimension; ++i)
          {
            pmin[i] = std::min(pmin[i], (*it)[i] - 1);
            pmax[i] = std::max(pmax[i], (*it)[i] + 1);
          }
        }
 
        for (auto it = sites.begin(); it != sites.end(); ++it) 
        {
          for (size_t i = 0; i < Space::dimension; ++i)
          {
            pmin[i] = std::min(pmin[i], (*it)[i] - 1);
            pmax[i] = std::max(pmax[i], (*it)[i] + 1);
          }
        }
 
        Domain domain(pmin, pmax);
 
        DigitalSet set(domain); set.insert(sites.begin(), sites.end());
        VoronoiPointPredicate predicate(set);
        Metric metric;
 
 
        typename Map::PeriodicitySpec specs = {false, false, false};
        if (params["toroidal-x"].as<int>()) specs[0] = true;
        if (params["toroidal-y"].as<int>()) specs[1] = true;
        if (params["toroidal-z"].as<int>()) specs[2] = true;
 
        auto map = Map(domain, predicate, metric, specs);
 
        std::vector<Vector> directions(pCount);
        size_t i = 0;
        for (auto it = points.begin(); it != points.end(); ++it)
        {
          directions[i++] = map(*it);
        }
        return directions;
      }
 
      //    - toroidal-x [false]: If the domain is toroidal in the first  dimension
      //    - toroidal-y [false]: If the domain is toroidal in the second dimension
      //    - toroidal-z [false]: If the domain is toroidal in the third  dimension
      template<uint32_t p, typename PointRangeSites, typename PointRange>
      static std::vector<typename ExactPredicateLpSeparableMetric<Space, p>::Value> getDistanceToClosestSite(
        const PointRange& points, 
        const PointRangeSites& sites,
        const Parameters& params = parametersVoronoiMap())
      {
        using Metric = ExactPredicateLpSeparableMetric<Space, p>;
        using Map = DistanceTransformation<Space, VoronoiPointPredicate, Metric>;
 
        // Compute domain of points
        Point pmin = *points.begin();
        Point pmax = pmin;
 
        size_t pCount = 0;
        for (auto it = points.begin(); it != points.end(); ++it) 
        {
          pCount ++;
          for (size_t i = 0; i < Space::dimension; ++i)
          {
            pmin[i] = std::min(pmin[i], (*it)[i] - 1);
            pmax[i] = std::max(pmax[i], (*it)[i] + 1);
          }
        }
 
        for (auto it = sites.begin(); it != sites.end(); ++it) 
        {
          for (size_t i = 0; i < Space::dimension; ++i)
          {
            pmin[i] = std::min(pmin[i], (*it)[i] - 1);
            pmax[i] = std::max(pmax[i], (*it)[i] + 1);
          }
        }
 
        Domain domain(pmin, pmax);
 
        DigitalSet set(domain); set.insert(sites.begin(), sites.end());
        VoronoiPointPredicate predicate(set);
        Metric metric;
 
        typename Map::PeriodicitySpec specs = {false, false, false};
        if (params["toroidal-x"].as<int>()) specs[0] = true;
        if (params["toroidal-y"].as<int>()) specs[1] = true;
        if (params["toroidal-z"].as<int>()) specs[2] = true;
 
        auto map = Map(domain, predicate, metric, specs);
 
        std::vector<typename Metric::Value> directions(pCount);
        size_t i = 0;
        for (auto it = points.begin(); it != points.end(); ++it)
        {
          directions[i++] = map(*it);
        }
        return directions;
      }
 
 
 
      // ----------------------- Standard services ------------------------------
    public:
 
      ShortcutsGeometry() = delete;
 
      ~ShortcutsGeometry() = delete;
 
      ShortcutsGeometry ( const ShortcutsGeometry & other ) = delete;
 
      ShortcutsGeometry ( ShortcutsGeometry && other ) = delete;
 
      ShortcutsGeometry & operator= ( const ShortcutsGeometry & other ) = delete;
 
      ShortcutsGeometry & operator= ( ShortcutsGeometry && other ) = delete;
 
 
      // ----------------------- Interface --------------------------------------
    public:
 
      // ------------------------- Protected Datas ------------------------------
    protected:
 
      // ------------------------- Private Datas --------------------------------
    private:
 
      // ------------------------- Hidden services ------------------------------
    protected:
 
      // ------------------------- Internals ------------------------------------
    private:
 
    }; // end of class ShortcutsGeometry
 
 
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined ShortcutsGeometry_h
 
#undef ShortcutsGeometry_RECURSES
#endif // else defined(ShortcutsGeometry_RECURSES)