Shortcuts.h

 
#pragma once
 
#if defined(Shortcuts_RECURSES)
#error Recursive header files inclusion detected in Shortcuts.h
#else // defined(Shortcuts_RECURSES)
#define Shortcuts_RECURSES
 
#if !defined Shortcuts_h
#define Shortcuts_h
 
// Inclusions
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>
#include <tuple>
#include <iterator>
#include <string>
#include "DGtal/base/Common.h"
#include "DGtal/base/CountedPtr.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/kernel/RegularPointEmbedder.h"
#include "DGtal/math/MPolynomial.h"
#include "DGtal/math/Statistic.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/images/IntervalForegroundPredicate.h"
#include <DGtal/images/ImageLinearCellEmbedder.h>
#include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
#include "DGtal/shapes/GaussDigitizer.h"
#include "DGtal/shapes/ShapeGeometricFunctors.h"
#include "DGtal/shapes/MeshHelpers.h"
#include "DGtal/topology/CCellularGridSpaceND.h"
#include "DGtal/topology/LightImplicitDigitalSurface.h"
#include "DGtal/topology/SetOfSurfels.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/topology/IndexedDigitalSurface.h"
#include "DGtal/topology/SurfelAdjacency.h"
#include "DGtal/topology/CCellEmbedder.h"
#include "DGtal/topology/CanonicCellEmbedder.h"
#include "DGtal/topology/CanonicSCellEmbedder.h"
#include "DGtal/topology/helpers/Surfaces.h"
#include "DGtal/geometry/volumes/KanungoNoise.h"
#include "DGtal/shapes/Mesh.h"
#include "DGtal/io/Color.h"
#include "DGtal/io/colormaps/GradientColorMap.h"
#include "DGtal/io/colormaps/TickedColorMap.h"
#include "DGtal/io/readers/MPolynomialReader.h"
#include "DGtal/io/readers/GenericReader.h"
#include "DGtal/io/writers/GenericWriter.h"
#include "DGtal/io/writers/MeshWriter.h"
#include "DGtal/graph/BreadthFirstVisitor.h"
#include "DGtal/graph/DepthFirstVisitor.h"
#include "DGtal/graph/GraphVisitorRange.h"
#include "DGtal/helpers/Parameters.h"
 
namespace DGtal
{
 
  // template class Shortcuts
  template  < typename TKSpace >
    class Shortcuts
    {
      BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< TKSpace > ));
 
      // ----------------------- 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< SCell >                                SCellRange;
      typedef std::vector< Cell >                                 CellRange;
      typedef CellRange                                           PointelRange;
      typedef SCellRange                                          SurfelRange;
      typedef std::vector< IdxSurfel >                            IdxSurfelRange;
      typedef std::vector< Scalar >                               Scalars;
      typedef std::vector< RealVector >                           RealVectors;
      typedef std::vector< RealPoint >                            RealPoints;
      typedef IdxVertex                                           Idx;
      typedef std::vector< IdxVertex >                            IdxRange;
 
      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 ::DGtal::Color                                      Color;
      typedef std::vector< Color >                                Colors;
      typedef GradientColorMap<Scalar>                            ColorMap;
      typedef TickedColorMap<Scalar,ColorMap>                     ZeroTickedColorMap;
    
      // ----------------------- Static services --------------------------------------
    public:
 
      // ----------------------- General static services ------------------------------
    public:
    
      static Parameters defaultParameters()
      {
        return parametersImplicitShape3D()
          | parametersKSpace()
          | parametersDigitizedImplicitShape3D()
          | parametersBinaryImage()
          | parametersGrayScaleImage()
          | parametersDigitalSurface()
          | parametersMesh()
          | parametersUtilities();
      }
 
      // ----------------------- ImplicitShape3D static services ------------------------
    public:
 
      static std::map< std::string, std::string >
        getPolynomialList()
        {
          std::vector< std::pair< std::string, std::string > >
            Ps = { { "sphere1", "x^2+y^2+z^2-1" },
                   { "sphere9", "x^2+y^2+z^2-81" },
                   { "ellipsoid", "3*x^2+2*y^2+z^2-90" },
                   { "cylinder", "x^2+2*z^2-90" },
                   { "torus",   "(x^2+y^2+z^2+6*6-2*2)^2-4*6*6*(x^2+y^2)" },
                   { "rcube",   "x^4+y^4+z^4-6561" },
                   { "goursat", "-1*(8-0.03*x^4-0.03*y^4-0.03*z^4+2*x^2+2*y^2+2*z^2)" },
                   { "goursat-hole", "x^4+y^4+z^4-2*4*(x^2+y^2+z^2)+2*4*4-2" },
                   { "distel",  "10000-(x^2+y^2+z^2+1000*(x^2+y^2)*(x^2+z^2)*(y^2+z^2))"},
                   { "leopold", "(x^2*y^2*z^2+4*x^2+4*y^2+3*z^2)-100" },
                   { "diabolo", "x^2-(y^2+z^2)^2" },
                   { "heart",   "-1*(x^2+2.25*y^2+z^2-1)^3+x^2*z^3+0.1125*y^2*z^3" },
                   { "crixxi",  "-0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3" } };
          std::map< std::string, std::string > L;
          for ( auto p : Ps )
            L[ p.first ] = p.second;
          return L;
        }
 
      static Parameters parametersImplicitShape3D()
      {
        return Parameters
          ( "polynomial", "sphere1" )
          ( "projectionMaxIter", 20 )
          ( "projectionAccuracy", 0.0001 )
          ( "projectionGamma",    0.5 );
      }
    
      static CountedPtr<ImplicitShape3D>
        makeImplicitShape3D( const Parameters& params = parametersImplicitShape3D() )
      {
        typedef MPolynomialReader< Space::dimension, Scalar> Polynomial3Reader;
        std::string poly_str = params[ "polynomial" ].as<std::string>();
        // Recognizes specific strings as polynomials.
        auto PL = getPolynomialList();
        if ( PL[ poly_str ] != "" ) poly_str = PL[ poly_str ];
        ScalarPolynomial poly;
        Polynomial3Reader reader;
        std::string::const_iterator iter
          = reader.read( poly, poly_str.begin(), poly_str.end() );
        if ( iter != poly_str.end() )
          {
            trace.error() << "[Shortcuts::makeImplicitShape3D]"
                          << " ERROR reading polynomial: I read only <"
                          << poly_str.substr( 0, iter - poly_str.begin() )
                          << ">, and I built P=" << poly << std::endl;
          }
        return CountedPtr<ImplicitShape3D>( new ImplicitShape3D( poly ) );
      }
 
      // ----------------------- KSpace static services ------------------------------
    public:
    
      static Parameters parametersKSpace()
      {
        return Parameters
          ( "closed",    1   )
          ( "gridsizex", 1.0 )
          ( "gridsizey", 1.0 )
          ( "gridsizez", 1.0 );
      }
 
      static KSpace getKSpace( const Point& low, const Point& up,
                               Parameters params = parametersKSpace() )
      {
        int closed  = params[ "closed"  ].as<int>();
        KSpace K;
        if ( ! K.init( low, up, closed ) )
          trace.error() << "[Shortcuts::getKSpace]"
                        << " Error building Khalimsky space K=" << K << std::endl;
        return K;
      }
 
      static KSpace getKSpace( CountedPtr<BinaryImage> bimage,
                               Parameters params = parametersKSpace() )
      {
        int closed  = params[ "closed"  ].as<int>();
        KSpace K;
        if ( ! K.init( bimage->domain().lowerBound(),
                       bimage->domain().upperBound(),
                       closed ) )
          trace.error() << "[Shortcuts::getKSpace]"
                        << " Error building Khalimsky space K=" << K << std::endl;
        return K;
      }
 
      static KSpace getKSpace( CountedPtr<GrayScaleImage> gimage,
                               Parameters params = parametersKSpace() )
      {
        int closed  = params[ "closed"  ].as<int>();
        KSpace K;
        if ( ! K.init( gimage->domain().lowerBound(),
                       gimage->domain().upperBound(),
                       closed ) )
          trace.error() << "[Shortcuts::getKSpace]"
                        << " Error building Khalimsky space K=" << K << std::endl;
        return K;
      }
 
      template <typename TDigitalSurfaceContainer>
      static KSpace
      getKSpace
      ( CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer> > surface )
      {
        return surface->container().space();
      }
 
      template <typename TDigitalSurfaceContainer>
      static KSpace
      getKSpace
      ( CountedPtr< ::DGtal::IndexedDigitalSurface< TDigitalSurfaceContainer> > surface )
      {
        return surface->container().space();
      }
      
      template <typename TDigitalSurfaceContainer>
      static const KSpace&
      refKSpace
      ( CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer> > surface )
      {
        return surface->container().space();
      }
      
      template <typename TDigitalSurfaceContainer>
      static const KSpace&
      refKSpace
      ( CountedPtr< ::DGtal::IndexedDigitalSurface< TDigitalSurfaceContainer> > surface )
      {
        return surface->container().space();
      }
 
      static CanonicCellEmbedder<KSpace>
        getCellEmbedder( const KSpace& K )
      {
        return CanonicCellEmbedder<KSpace>( K );
      }
 
      static CanonicSCellEmbedder<KSpace>
        getSCellEmbedder( const KSpace& K )
      {
        return CanonicSCellEmbedder<KSpace>( K );
      }
      
    
      // ----------------------- DigitizedImplicitShape3D static services --------------
    public:
    
      static Parameters parametersDigitizedImplicitShape3D()
      {
        return Parameters
          ( "minAABB",  -10.0 )
          ( "maxAABB",   10.0 )
          ( "gridstep",   1.0 )
          ( "offset",     5.0 );
      }
 
    
      static KSpace
        getKSpace( Parameters params =
                   parametersKSpace() | parametersDigitizedImplicitShape3D() )
      {
        Scalar min_x  = params[ "minAABB"  ].as<Scalar>();
        Scalar max_x  = params[ "maxAABB"  ].as<Scalar>();
        Scalar h      = params[ "gridstep" ].as<Scalar>();
        Scalar offset = params[ "offset"   ].as<Scalar>();
        bool   closed = params[ "closed"   ].as<int>();
        RealPoint p1( min_x - offset * h, min_x - offset * h, min_x - offset * h );
        RealPoint p2( max_x + offset * h, max_x + offset * h, max_x + offset * h );
        CountedPtr<DigitizedImplicitShape3D> dshape( new DigitizedImplicitShape3D() );
        dshape->init( p1, p2, h );
        Domain domain = dshape->getDomain();
        KSpace K;
        if ( ! K.init( domain.lowerBound(), domain.upperBound(), closed ) )
          trace.error() << "[Shortcuts::getKSpace]"
                        << " Error building Khalimsky space K=" << K << std::endl
                        << "Note: if you use decimal values, check your locale for decimal point '.' or ','."
                        << std::endl;
        return K;
      }
 
      static CountedPtr<DigitizedImplicitShape3D>
        makeDigitizedImplicitShape3D
        ( CountedPtr<ImplicitShape3D> shape,
          Parameters params = parametersDigitizedImplicitShape3D() )
      {
        Scalar min_x  = params[ "minAABB"  ].as<Scalar>();
        Scalar max_x  = params[ "maxAABB"  ].as<Scalar>();
        Scalar h      = params[ "gridstep" ].as<Scalar>();
        Scalar offset = params[ "offset"   ].as<Scalar>();
        RealPoint p1( min_x - offset * h, min_x - offset * h, min_x - offset * h );
        RealPoint p2( max_x + offset * h, max_x + offset * h, max_x + offset * h );
        CountedPtr<DigitizedImplicitShape3D> dshape( new DigitizedImplicitShape3D() );
        dshape->attach( shape );
        dshape->init( p1, p2, h );
        return dshape;
      }
 
 
      // ----------------------- BinaryImage static services --------------------------
    public:
    
      static Parameters parametersBinaryImage()
      {
        return Parameters
          ( "noise", 0.0 )
          ( "thresholdMin", 0 )
          ( "thresholdMax", 255 );
      }
    
      static CountedPtr<BinaryImage>
        makeBinaryImage( Domain shapeDomain )
      {
        return CountedPtr<BinaryImage>( new BinaryImage( shapeDomain ) );
      }
    
      static CountedPtr<BinaryImage>
        makeBinaryImage( CountedPtr<DigitizedImplicitShape3D> shape_digitization,
                         Parameters params = parametersBinaryImage() )
      {
        return makeBinaryImage( shape_digitization,
                                shape_digitization->getDomain(),
                                params );
      }
    
      static CountedPtr<BinaryImage>
        makeBinaryImage( CountedPtr<DigitizedImplicitShape3D> shape_digitization,
                         Domain shapeDomain,
                         Parameters params = parametersBinaryImage() )
      {
        const Scalar noise        = params[ "noise"  ].as<Scalar>();
        CountedPtr<BinaryImage> img ( new BinaryImage( shapeDomain ) );
        if ( noise <= 0.0 )
          {
            std::transform( shapeDomain.begin(), shapeDomain.end(),
                            img->begin(),
                            [&shape_digitization]
                            ( const Point& p ) { return (*shape_digitization)(p); } );
          }
        else
          {
            typedef KanungoNoise< DigitizedImplicitShape3D, Domain > KanungoPredicate;
            KanungoPredicate noisy_dshape( *shape_digitization, shapeDomain, noise );
            std::transform( shapeDomain.begin(), shapeDomain.end(),
                            img->begin(),
                            [&noisy_dshape] ( const Point& p ) { return noisy_dshape(p); } );
          }
        return img;
      }
 
      static CountedPtr<BinaryImage>
        makeBinaryImage( CountedPtr<BinaryImage> bimage,
                         Parameters params = parametersBinaryImage() )
      {
        const Scalar noise = params[ "noise"  ].as<Scalar>();
        if ( noise <= 0.0 ) return bimage;
        typedef KanungoNoise< BinaryImage, Domain > KanungoPredicate;
        const Domain shapeDomain    = bimage->domain();
        CountedPtr<BinaryImage> img ( new BinaryImage( shapeDomain ) );
        KanungoPredicate noisy_dshape( *bimage, shapeDomain, noise );
        std::transform( shapeDomain.begin(), shapeDomain.end(),
                        img->begin(),
                        [&noisy_dshape] ( const Point& p ) { return noisy_dshape(p); } );
        return img;
      }
 
      static CountedPtr<BinaryImage>
        makeBinaryImage
        ( std::string input,
          Parameters params = parametersBinaryImage() )
      {
        int     thresholdMin = params["thresholdMin"].as<int>();
        int     thresholdMax = params["thresholdMax"].as<int>();
        GrayScaleImage image = GenericReader<GrayScaleImage>::import( input );
        Domain        domain = image.domain();
        typedef functors::IntervalForegroundPredicate<GrayScaleImage> ThresholdedImage;
        ThresholdedImage tImage( image, thresholdMin, thresholdMax );
        CountedPtr<BinaryImage> img ( new BinaryImage( domain ) );
        std::transform( domain.begin(), domain.end(),
                        img->begin(),
                        [tImage] ( const Point& p ) { return tImage(p); } );
        return makeBinaryImage( img, params );
      }
 
      static CountedPtr<BinaryImage>
        makeBinaryImage
        ( CountedPtr<GrayScaleImage> gray_scale_image,
          Parameters params = parametersBinaryImage() )
      {
        int     thresholdMin = params["thresholdMin"].as<int>();
        int     thresholdMax = params["thresholdMax"].as<int>();
        Domain        domain = gray_scale_image->domain();
        typedef functors::IntervalForegroundPredicate<GrayScaleImage> ThresholdedImage;
        ThresholdedImage tImage( *gray_scale_image, thresholdMin, thresholdMax );
        CountedPtr<BinaryImage> img ( new BinaryImage( domain ) );
        std::transform( domain.begin(), domain.end(),
                        img->begin(),
                        [tImage] ( const Point& p ) { return tImage(p); } );
        return makeBinaryImage( img, params );
      }
 
    
      static bool
        saveBinaryImage
        ( CountedPtr<BinaryImage> bimage, std::string output )
      {
        auto gray_scale_image = makeGrayScaleImage( bimage );
        return saveGrayScaleImage( gray_scale_image, output );
      }
 
 
      // ----------------------- GrayScaleImage static services -------------------------
    public:
 
      static Parameters parametersGrayScaleImage()
      {
        return Parameters
          ( "qShift",   128.0 )
          ( "qSlope",     1.0 );
      }
    
      static CountedPtr<GrayScaleImage>
        makeGrayScaleImage( Domain aDomain )
      {
        return CountedPtr<GrayScaleImage>( new GrayScaleImage( aDomain ) );
      }
 
      static CountedPtr<GrayScaleImage>
        makeGrayScaleImage
        ( std::string input )
      {
        GrayScaleImage image = GenericReader<GrayScaleImage>::import( input );
        return CountedPtr<GrayScaleImage>( new GrayScaleImage( image ) );
      }
 
      static CountedPtr<GrayScaleImage>
        makeGrayScaleImage
        ( CountedPtr<BinaryImage> binary_image,
          std::function< GrayScale( bool ) > const & bool2grayscale
          = [] ( bool v ) { return v ? (unsigned char) 255 : (unsigned char) 0; }
          // JOL: (GrayScale) was not working with visual C++: error C2065
          )
      {
        const Domain domain = binary_image->domain(); 
        CountedPtr<GrayScaleImage> gray_scale_image( new GrayScaleImage( domain ) );
        std::transform( binary_image->begin(), binary_image->end(),
                        gray_scale_image->begin(),
                        bool2grayscale );
        return gray_scale_image;
      }
 
      static bool
        saveGrayScaleImage
        ( CountedPtr<GrayScaleImage> gray_scale_image, std::string output )
      {
        return GenericWriter< GrayScaleImage > 
          ::exportFile( output, *gray_scale_image );
      }
 
 
      static CountedPtr<GrayScaleImage>
        makeGrayScaleImage
        ( CountedPtr<FloatImage> fimage,
          Parameters params = parametersGrayScaleImage() )
      {
        float qShift = params[ "qShift"   ].as<float>();
        float qSlope = params[ "qSlope"   ].as<float>();
        std::function< unsigned char( float ) > f
          = [qShift,qSlope] (float v)
          { return (unsigned char) std::min( 255.0f, std::max( 0.0f, qSlope * v + qShift ) ); };
        Domain domain = fimage->domain();
        auto   gimage = makeGrayScaleImage( domain );
        auto       it = gimage->begin();
        for ( auto p : domain )
          {
            float val = (*fimage)( p );
            *it++      = f( val );
          }
        return gimage;
      }
 
      static CountedPtr<GrayScaleImage>
        makeGrayScaleImage
        ( CountedPtr<DoubleImage> fimage,
          Parameters params = parametersGrayScaleImage() )
      {
        double qShift = params[ "qShift"   ].as<double>();
        double qSlope = params[ "qSlope"   ].as<double>();
        std::function< unsigned char( double ) > f
          = [qShift,qSlope] (double v)
          { return (unsigned char) std::min( 255.0, std::max( 0.0, qSlope * v + qShift ) ); };
        Domain domain = fimage->domain();
        auto   gimage = makeGrayScaleImage( domain );
        auto       it = gimage->begin();
        for ( auto p : domain )
          {
            double val = (*fimage)( p );
            *it++      = f( val );
          }
        return gimage;
      }
 
      // ----------------------- FloatImage static services -------------------------
    public:
 
      static CountedPtr<FloatImage>
        makeFloatImage( Domain aDomain )
      {
        return CountedPtr<FloatImage>( new FloatImage( aDomain ) );
      }
 
      static CountedPtr<FloatImage>
        makeFloatImage
        ( std::string input )
      {
        FloatImage image = GenericReader<FloatImage>::import( input );
        return CountedPtr<FloatImage>( new FloatImage( image ) );
      }
 
      static CountedPtr<FloatImage>
        makeFloatImage
        ( CountedPtr<ImplicitShape3D> shape,
          Parameters params = parametersDigitizedImplicitShape3D() )
      {
        Scalar min_x  = params[ "minAABB"  ].as<Scalar>();
        Scalar max_x  = params[ "maxAABB"  ].as<Scalar>();
        Scalar h      = params[ "gridstep" ].as<Scalar>();
        Scalar offset = params[ "offset"   ].as<Scalar>();
        RealPoint p1( min_x - offset * h, min_x - offset * h, min_x - offset * h );
        RealPoint p2( max_x + offset * h, max_x + offset * h, max_x + offset * h );
        CountedPtr<DigitizedImplicitShape3D> dshape( new DigitizedImplicitShape3D() );
        dshape->attach( shape );
        dshape->init( p1, p2, h );
        Domain domain = dshape->getDomain();
        auto   fimage = makeFloatImage( domain );
        auto       it = fimage->begin();
        for ( auto p : domain )
          {
            float val = (float) (*shape)( p );
            *it++      = val;
          }
        return fimage;
      }
 
      // ----------------------- DoubleImage static services -------------------------
    public:
 
      static CountedPtr<DoubleImage>
        makeDoubleImage( Domain aDomain )
      {
        return CountedPtr<DoubleImage>( new DoubleImage( aDomain ) );
      }
 
      static CountedPtr<DoubleImage>
        makeDoubleImage
        ( std::string input )
      {
        DoubleImage image = GenericReader<DoubleImage>::import( input );
        return CountedPtr<DoubleImage>( new DoubleImage( image ) );
      }
 
      static CountedPtr<DoubleImage>
        makeDoubleImage
        ( CountedPtr<ImplicitShape3D> shape,
          Parameters params = parametersDigitizedImplicitShape3D() )
      {
        Scalar min_x  = params[ "minAABB"  ].as<Scalar>();
        Scalar max_x  = params[ "maxAABB"  ].as<Scalar>();
        Scalar h      = params[ "gridstep" ].as<Scalar>();
        Scalar offset = params[ "offset"   ].as<Scalar>();
        RealPoint p1( min_x - offset * h, min_x - offset * h, min_x - offset * h );
        RealPoint p2( max_x + offset * h, max_x + offset * h, max_x + offset * h );
        CountedPtr<DigitizedImplicitShape3D> dshape( new DigitizedImplicitShape3D() );
        dshape->attach( shape );
        dshape->init( p1, p2, h );
        Domain domain = dshape->getDomain();
        auto   fimage = makeDoubleImage( domain );
        auto       it = fimage->begin();
        for ( auto p : domain )
          {
            double val = (double) (*shape)( p );
            *it++      = val;
          }
        return fimage;
      }
    
      // ----------------------- DigitalSurface static services ------------------------
    public:
        
      static Parameters parametersDigitalSurface()
      {
        return Parameters
          ( "surfelAdjacency",   0 )
          ( "nbTriesToFindABel", 100000 )
          ( "surfaceComponents", "AnyBig" )
          ( "surfaceTraversal",  "Default" );
      }
 
      template <typename TDigitalSurfaceContainer>
        static CanonicCellEmbedder<KSpace>
        getCellEmbedder
        ( CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer> > surface )
        {
          return getCellEmbedder( refKSpace( surface ) );
        }
 
      template <typename TDigitalSurfaceContainer>
        static CanonicSCellEmbedder<KSpace>
        getSCellEmbedder
        ( CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer> > surface )
        {
          return getSCellEmbedder( refKSpace( surface ) );
        }
 
      template <typename TDigitalSurfaceContainer>
        static CanonicCellEmbedder<KSpace>
        getCellEmbedder
        ( CountedPtr< ::DGtal::IndexedDigitalSurface< TDigitalSurfaceContainer> > surface )
        {
          return getCellEmbedder( refKSpace( surface ) );
        }
 
      template <typename TDigitalSurfaceContainer>
        static CanonicSCellEmbedder<KSpace>
        getSCellEmbedder
        ( CountedPtr< ::DGtal::IndexedDigitalSurface< TDigitalSurfaceContainer> > surface )
        {
          return getSCellEmbedder( refKSpace( surface ) );
        }
    
      static CountedPtr<LightDigitalSurface>
        makeLightDigitalSurface
        ( CountedPtr<BinaryImage> bimage,
          const KSpace&           K,
          const Parameters&       params = parametersDigitalSurface() )
      {
        bool surfel_adjacency      = params[ "surfelAdjacency" ].as<int>();
        int nb_tries_to_find_a_bel = params[ "nbTriesToFindABel" ].as<int>();
        SurfelAdjacency< KSpace::dimension > surfAdj( surfel_adjacency );
 
        // We have to search for a surfel that belongs to a big connected component.
        CountedPtr<LightDigitalSurface> ptrSurface;
        Surfel       bel;
        Scalar       minsize    = bimage->extent().norm();
        unsigned int nb_surfels = 0;
        unsigned int tries      = 0;
        do
          {
            try { // Search initial bel
              bel = Surfaces<KSpace>::findABel( K, *bimage, nb_tries_to_find_a_bel );
            } catch (DGtal::InputException& e) {
              trace.error() << "[Shortcuts::makeLightDigitalSurface]"
                            << " ERROR Unable to find bel. " << e << std::endl;
              return ptrSurface;
            }
            // this pointer will be acquired by the surface.
            LightSurfaceContainer* surfContainer
              = new LightSurfaceContainer( K, *bimage, surfAdj, bel );
            ptrSurface = CountedPtr<LightDigitalSurface>
              ( new LightDigitalSurface( surfContainer ) ); // acquired
            nb_surfels = ptrSurface->size();
          }
        while ( ( nb_surfels < 2 * minsize ) && ( tries++ < 150 ) );
        if( tries >= 150 )
          trace.warning() << "[Shortcuts::makeLightDigitalSurface]"
                          << "ERROR cannot find a proper bel in a big enough component."
                          << std::endl;
        return ptrSurface;
      }
 
      static std::vector< CountedPtr<LightDigitalSurface> >
        makeLightDigitalSurfaces
        ( CountedPtr<BinaryImage> bimage,
          const KSpace&           K,
          const Parameters&       params = parametersDigitalSurface() )
      {
        SurfelRange surfel_reps;
        return makeLightDigitalSurfaces( surfel_reps, bimage, K, params );
      }
 
      static std::vector< CountedPtr<LightDigitalSurface> >
        makeLightDigitalSurfaces
        ( SurfelRange&            surfel_reps,
          CountedPtr<BinaryImage> bimage,
          const KSpace&           K,
          const Parameters&       params = parametersDigitalSurface() )
      {
        std::vector< CountedPtr<LightDigitalSurface> > result;
        std::string component      = params[ "surfaceComponents" ].as<std::string>();
        if ( component == "AnyBig" )
          {
            result.push_back( makeLightDigitalSurface( bimage, K, params ) );
            surfel_reps.push_back( *( result[ 0 ]->begin() ) );
            return result;
          }     
        bool surfel_adjacency      = params[ "surfelAdjacency" ].as<int>();
        SurfelAdjacency< KSpace::dimension > surfAdj( surfel_adjacency );
        // Extracts all boundary surfels
        SurfelSet all_surfels;
        Surfaces<KSpace>::sMakeBoundary( all_surfels, K, *bimage,
                                         K.lowerBound(), K.upperBound() );
        // Builds all connected components of surfels.
        SurfelSet marked_surfels;
        CountedPtr<LightDigitalSurface> ptrSurface;
        for ( auto bel : all_surfels )
          {
            if ( marked_surfels.count( bel ) != 0 ) continue;
            surfel_reps.push_back( bel );
            LightSurfaceContainer* surfContainer
              = new LightSurfaceContainer( K, *bimage, surfAdj, bel );
            ptrSurface = CountedPtr<LightDigitalSurface>
              ( new LightDigitalSurface( surfContainer ) ); // acquired
            // mark all surfels of the surface component.
            marked_surfels.insert( ptrSurface->begin(), ptrSurface->end() );
            // add surface component to result.
            result.push_back( ptrSurface );
          }
        return result;
      }
 
    
      template <typename TPointPredicate>
        static CountedPtr< DigitalSurface >
        makeDigitalSurface
        ( CountedPtr< TPointPredicate > bimage,
          const KSpace&           K,
          const Parameters&       params = parametersDigitalSurface() )
        {
          SurfelSet all_surfels;
          bool      surfel_adjacency = params[ "surfelAdjacency" ].as<int>();
          SurfelAdjacency< KSpace::dimension > surfAdj( surfel_adjacency );
          // Extracts all boundary surfels
          Surfaces<KSpace>::sMakeBoundary( all_surfels, K, *bimage,
                                           K.lowerBound(), K.upperBound() );
          ExplicitSurfaceContainer* surfContainer
            = new ExplicitSurfaceContainer( K, surfAdj, all_surfels );
          return CountedPtr< DigitalSurface >
            ( new DigitalSurface( surfContainer ) ); // acquired
        }
 
 
      static CountedPtr< DigitalSurface >
        makeDigitalSurface
        ( CountedPtr<IdxDigitalSurface> idx_surface,
          const Parameters&             params = parametersDigitalSurface() )
      {
        bool surfel_adjacency      = params[ "surfelAdjacency" ].as<int>();
        const KSpace& K = refKSpace( idx_surface );
        SurfelAdjacency< KSpace::dimension > surfAdj( surfel_adjacency );
        auto all_idx_surfels
          = getIdxSurfelRange( idx_surface, Parameters( "surfaceTraversal", "Default" ) );
        auto idx2surfel = idx_surface->surfels();
        SurfelSet all_surfels;
        for ( auto idx : all_idx_surfels ) all_surfels.insert( idx2surfel[ idx ] );
        ExplicitSurfaceContainer* surfContainer
          = new ExplicitSurfaceContainer( K, surfAdj, all_surfels );
        return CountedPtr<DigitalSurface>
          ( new DigitalSurface( surfContainer ) ); // acquired
      }
    
      static CountedPtr<IdxDigitalSurface>
        makeIdxDigitalSurface
        ( CountedPtr<BinaryImage> bimage,
          const KSpace&           K,
          const Parameters&       params = parametersDigitalSurface() )
      {
        std::string component      = params[ "surfaceComponents" ].as<std::string>();
        SurfelSet surfels;
        if ( component == "AnyBig" )
          {
            auto light_surface = makeLightDigitalSurface( bimage, K, params );
            surfels.insert( light_surface->begin(), light_surface->end() );
          }
        else if ( component == "All" )
          {
            Surfaces<KSpace>::sMakeBoundary( surfels, K, *bimage,
                                             K.lowerBound(), K.upperBound() );
          }
        return makeIdxDigitalSurface( surfels, K, params );
      }    
 
      template <typename TSurfelRange>
        static CountedPtr<IdxDigitalSurface>
        makeIdxDigitalSurface
        ( const TSurfelRange&  surfels,
          ConstAlias< KSpace > K,
          const Parameters&    params = parametersDigitalSurface() )
        {
          bool surfel_adjacency      = params[ "surfelAdjacency" ].as<int>();
          SurfelAdjacency< KSpace::dimension > surfAdj( surfel_adjacency );
          // Build indexed digital surface.
          CountedPtr<ExplicitSurfaceContainer> ptrSurfContainer
            ( new ExplicitSurfaceContainer( K, surfAdj, surfels ) );
          CountedPtr<IdxDigitalSurface> ptrSurface
            ( new IdxDigitalSurface() );
          bool ok = ptrSurface->build( ptrSurfContainer );
          if ( !ok )
            trace.warning() << "[Shortcuts::makeIdxDigitalSurface]"
                            << " Error building indexed digital surface." << std::endl;
          return ptrSurface;
        }
 
      template <typename TDigitalSurfaceContainer>
        static CountedPtr<IdxDigitalSurface>
        makeIdxDigitalSurface
        ( CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer> > surface,
          const Parameters&               params = parametersDigitalSurface() )
        {
          const KSpace& K = refKSpace( surface );
          SurfelSet     surfels;
          surfels.insert( surface->begin(), surface->end() );
          return makeIdxDigitalSurface( surfels, K, params );
        }    
 
      static CountedPtr<IdxDigitalSurface>
        makeIdxDigitalSurface
        ( const std::vector< CountedPtr<LightDigitalSurface> >& surfaces,
          const Parameters&       params = parametersDigitalSurface()  )
      {
        if ( surfaces.empty() ) return CountedPtr<IdxDigitalSurface>( 0 );
        const KSpace& K = surfaces[ 0 ]->container().space();
        SurfelSet     surfels;
        for ( std::size_t i = 0; i < surfaces.size(); ++i )
          {
            const KSpace& Ki = surfaces[ i ]->container().space();
            if ( ( Ki.lowerBound() != K.lowerBound() )
                 || ( Ki.upperBound() != K.upperBound() ) )
              trace.warning() << "[Shortcuts::makeIdxDigitalSurface]"
                              << " Incompatible digital spaces for surface " << i << std::endl;
            surfels.insert( surfaces[ i ]->begin(), surfaces[ i ]->end() );
          }
        return makeIdxDigitalSurface( surfels, K, params );
      }    
 
    
      template <typename TDigitalSurfaceContainer>
      static CellRange
      getCellRange
      ( Cell2Index& c2i,
        CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > surface,
        const Dimension k )
      {
        CellRange result;
        // Approximate number of pointels given the number of 2-cells (valid for 2d surfaces in nD).
        result.reserve( 2 * surface->size() + 100 );
        const KSpace& K = refKSpace( surface );
        Idx n = 0;
        for ( auto&& surfel : *surface )
          {
            CellRange primal_cells = getPrimalCells( K, surfel, k );
            for ( auto&& primal_cell : primal_cells )
              {
                if ( ! c2i.count( primal_cell ) )
                  {
                    result.push_back( primal_cell );
                    c2i[ primal_cell ] = n++;
                  }
              }
          }
        return result;
      }
 
      template <typename TDigitalSurfaceContainer>
      static PointelRange
      getCellRange
      (
        CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > surface,
        const Dimension k )
      {
        Cell2Index c2i;
        return getCellRange( c2i, surface, k );
      }
      
      template <typename TDigitalSurfaceContainer>
        static PointelRange
        getPointelRange
        ( Cell2Index& c2i,
          CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > surface )
        {
          PointelRange result;
          result.reserve( surface->size() );
          const KSpace& K = refKSpace( surface );
          Idx n = 0;
          for ( auto&& surfel : *surface )
            {
              CellRange primal_vtcs = getPointelRange( K, surfel );
              for ( auto&& primal_vtx : primal_vtcs )
                {
                  if ( ! c2i.count( primal_vtx ) )
                    {
                      result.push_back( primal_vtx );
                      c2i[ primal_vtx ] = n++;
                    }
                }
            }
          return result;
        }
 
      template <typename TDigitalSurfaceContainer>
      static PointelRange
      getPointelRange
      ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > surface )
      {
        Cell2Index c2i;
        return getPointelRange( c2i, surface );
      }
 
      static PointelRange
      getPointelRange
      ( const KSpace& K, const SCell& surfel )
      {
        return KSpace::dimension == 3
          ? getPrimalVertices( K, surfel, true )
          : getPrimalVertices( K, surfel );
      }
      
      template <typename TDigitalSurfaceContainer>
      static SurfelRange
      getSurfelRange
      ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > surface,
        const Parameters&   params = parametersDigitalSurface() )
      {
        return getSurfelRange( surface, *( surface->begin() ), params );
      }
      
      template <typename TDigitalSurfaceContainer>
      static SurfelRange
      getSurfelRange
      ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > surface,
        const Surfel&       start_surfel,
        const Parameters&   params = parametersDigitalSurface() )
      {
        typedef ::DGtal::DigitalSurface<TDigitalSurfaceContainer> AnyDigitalSurface;
        SurfelRange result;
        std::string traversal = params[ "surfaceTraversal" ].as<std::string>();
        if ( traversal == "DepthFirst" )
          {
            typedef DepthFirstVisitor< AnyDigitalSurface > Visitor;
            typedef GraphVisitorRange< Visitor > VisitorRange;
            VisitorRange range( new Visitor( *surface, start_surfel ) );
            std::for_each( range.begin(), range.end(),
                           [&result] ( Surfel s ) { result.push_back( s ); } );
          }
        else if ( traversal == "BreadthFirst" )
          {
            typedef BreadthFirstVisitor< AnyDigitalSurface > Visitor;
            typedef GraphVisitorRange< Visitor > VisitorRange;
            VisitorRange range( new Visitor( *surface, start_surfel ) );
            std::for_each( range.begin(), range.end(),
                           [&result] ( Surfel s ) { result.push_back( s ); } );
          }
        else
          {
            std::for_each( surface->begin(), surface->end(),
                           [&result] ( Surfel s ) { result.push_back( s ); } );
          }
        return result;
      }
 
      static IdxSurfelRange
      getIdxSurfelRange
      ( CountedPtr<IdxDigitalSurface> surface,
        const Parameters&   params = parametersDigitalSurface() )
      {
        return getIdxSurfelRange( surface, (IdxSurfel) 0, params );
      }
    
      static IdxSurfelRange
      getIdxSurfelRange
      ( CountedPtr<IdxDigitalSurface> surface,
        const IdxSurfel&    start_surfel,
        const Parameters&   params = parametersDigitalSurface() )
      {
        IdxSurfelRange result;
        std::string traversal = params[ "surfaceTraversal" ].as<std::string>();
        if ( traversal == "DepthFirst" )
          {
            typedef DepthFirstVisitor< IdxDigitalSurface > Visitor;
            typedef GraphVisitorRange< Visitor > VisitorRange;
            VisitorRange range( new Visitor( *surface, start_surfel ) );
            std::for_each( range.begin(), range.end(),
                           [&result] ( IdxSurfel s ) { result.push_back( s ); } );
          }
        else if ( traversal == "BreadthFirst" )
          {
            typedef BreadthFirstVisitor< IdxDigitalSurface > Visitor;
            typedef GraphVisitorRange< Visitor > VisitorRange;
            VisitorRange range( new Visitor( *surface, start_surfel ) );
            std::for_each( range.begin(), range.end(),
                           [&result] ( IdxSurfel s ) { result.push_back( s ); } );
          }
        else return surface->allVertices();
        return result;
      }
      template <typename TDigitalSurfaceContainer,
        typename TCellEmbedder>
        static bool
        saveOFF
        ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > digsurf,
          const TCellEmbedder&           embedder,
          std::string                    off_file,
          const Color&                   face_color=DGtal::Color::None)
 
        {
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
 
          std::ofstream output_off( off_file.c_str() );
          output_off << "OFF" << std::endl;
          output_off << "# Generated from  DGtal::Shortcuts from the DGTal library" << std::endl;
          Cell2Index      c2i;
          auto       pointels = getPointelRange( c2i, digsurf );
          output_off <<  pointels.size()  << " " << digsurf->size() << " " << 0 << " " << std::endl;
          
          
          // Number and output vertices.
          const KSpace&     K = refKSpace( digsurf );
          for ( auto&& pointel : pointels )
            {
              RealPoint p = embedder( pointel );
              output_off << p[ 0 ] << " " << p[ 1 ] << " " << p[ 2 ] << std::endl;
            }   
          
          // Taking care of faces
          for ( auto&& surfel : *digsurf )
            {              
              auto primal_vtcs = getPointelRange( K, surfel );
              output_off << primal_vtcs.size();
                {
                  for ( auto&& primal_vtx : primal_vtcs )
                    output_off << " " << (c2i[ primal_vtx ]);
                }
                if(face_color != DGtal::Color::None)
                {
                  output_off << " ";
                  output_off  << face_color.r() << " " << face_color.g()
                              << " " << face_color.b() << " " << face_color.a();
                }
                output_off << std::endl;
            }
          return output_off.good();
        }
    
 
      template <typename TDigitalSurfaceContainer,
        typename TCellEmbedder>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > digsurf,
          const TCellEmbedder&           embedder,
          const RealVectors&             normals,
          const Colors&                  diffuse_colors,
          std::string                    objfile,
          const Color&                   ambient_color  = Color( 32, 32, 32 ),
          const Color&                   diffuse_color  = Color( 200, 200, 255 ),
          const Color&                   specular_color = Color::White )
        {
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          std::string mtlfile;
          auto lastindex = objfile.find_last_of(".");
          if ( lastindex == std::string::npos )
            {
              mtlfile  = objfile + ".mtl";
              objfile  = objfile + ".obj";
            }
          else
            {
              mtlfile  = objfile.substr(0, lastindex) + ".mtl"; 
            }
 
          std::ofstream output_obj( objfile.c_str() );
          output_obj << "#  OBJ format" << std::endl;
          output_obj << "# DGtal::MeshHelpers::exportOBJwithFaceNormalAndColor" << std::endl;
          output_obj << "o anObject" << std::endl;
                  //remove directory to write material
          auto indexpath = objfile.find_last_of("/");
          output_obj << "mtllib " << mtlfile.substr(indexpath+1) << std::endl;
          std::ofstream output_mtl( mtlfile.c_str() );
          output_mtl << "#  MTL format"<< std::endl;
          output_mtl << "# generated from MeshWriter from the DGTal library"<< std::endl;
          // Number and output vertices.
          const KSpace&     K = refKSpace( digsurf );
          Cell2Index      c2i;
          auto       pointels = getPointelRange( c2i, digsurf );
          for ( auto&& pointel : pointels )
            {
              RealPoint p = embedder( pointel );
              output_obj << "v " << p[ 0 ] << " " << p[ 1 ] << " " << p[ 2 ] << std::endl;
            }   
          // Taking care of normals
          Idx nbfaces = digsurf->size();
          bool has_normals = ( nbfaces == normals.size() );
          if ( has_normals )
            {
              for ( Idx f = 0; f < nbfaces; ++f )
                {
                  const auto& p = normals[ f ];
                  output_obj << "vn " << p[ 0 ] << " " << p[ 1 ] << " " << p[ 2 ] << std::endl;
                }
            }
          // Taking care of materials
          bool has_material = ( nbfaces == diffuse_colors.size() );
          Idx   idxMaterial = 0;
          std::map<Color, unsigned int > mapMaterial;
          if ( has_material )
            {
              for ( Idx f = 0; f < nbfaces; ++f )
                {
                  Color c = diffuse_colors[ f ];
                  if ( mapMaterial.count( c ) == 0 )
                    {
                      MeshHelpers::exportMTLNewMaterial
                        ( output_mtl, idxMaterial, ambient_color, c, specular_color );
                      mapMaterial[ c ] = idxMaterial++;
                    }
                }
            }
          else
            {
              MeshHelpers::exportMTLNewMaterial
                ( output_mtl, idxMaterial, ambient_color, diffuse_color, specular_color );
            }
      
          // Taking care of faces
          Idx f = 0;
          for ( auto&& surfel : *digsurf )
            {
              output_obj << "usemtl material_"
                         << ( has_material ? mapMaterial[ diffuse_colors[ f ] ] : idxMaterial )
                         << std::endl; 
              output_obj << "f";
              auto primal_vtcs = getPointelRange( K, surfel );
              // The +1 in lines below is because indexing starts at 1 in OBJ file format.
              if ( has_normals )
                {
                  for ( auto&& primal_vtx : primal_vtcs )
                    output_obj << " " << (c2i[ primal_vtx ]+1) << "//" << (f+1);
                }
              else
                {
                  for ( auto&& primal_vtx : primal_vtcs )
                    output_obj << " " << (c2i[ primal_vtx ]+1);
                }
              output_obj << std::endl;
              f += 1;
            }
          output_mtl.close();
          return output_obj.good();
        }
    
      template <typename TDigitalSurfaceContainer>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > digsurf,
          const RealVectors&                            normals,
          const Colors&                                 diffuse_colors,
          std::string                                   objfile,
          const Color&                   ambient_color  = Color( 32, 32, 32 ),
          const Color&                   diffuse_color  = Color( 200, 200, 255 ),
          const Color&                   specular_color = Color::White )
        {
          auto embedder = getCellEmbedder( digsurf );
          return saveOBJ( digsurf, embedder, normals, diffuse_colors, objfile,
                          ambient_color, diffuse_color, specular_color );
        }
 
      template <typename TDigitalSurfaceContainer>
        static bool
        saveOFF
        ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > digsurf,
          std::string                                   off_file,
          const Color&                   face_color  = Color( 32, 32, 32 ))
          {
          auto embedder = getCellEmbedder( digsurf );
          return saveOFF( digsurf, embedder, off_file, face_color);
        }
      template <typename TDigitalSurfaceContainer>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::DigitalSurface<TDigitalSurfaceContainer> > digsurf,
          std::string                                   objfile,
          const Color&                   ambient_color  = Color( 32, 32, 32 ),
          const Color&                   diffuse_color  = Color( 200, 200, 255 ),
          const Color&                   specular_color = Color::White )
        {
          auto embedder = getCellEmbedder( digsurf );
          return saveOBJ( digsurf, embedder, RealVectors(), Colors(), objfile,
                          ambient_color, diffuse_color, specular_color );
        }
 
 
      static bool
        saveVectorFieldOBJ
        ( const RealPoints&  positions,
          const RealVectors& vf,
          double             thickness,
          const Colors&      diffuse_colors,
          std::string        objfile,
          const Color&       ambient_color  = Color( 32, 32, 32 ),
          const Color&       diffuse_color  = Color( 200, 200, 255 ),
          const Color&       specular_color = Color::White )
        {
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          std::string mtlfile;
          auto lastindex = objfile.find_last_of(".");
          if ( lastindex == std::string::npos )
            {
              mtlfile  = objfile + ".mtl";
              objfile  = objfile + ".obj";
            }
          else
            {
              mtlfile  = objfile.substr(0, lastindex) + ".mtl"; 
            }
          std::ofstream output_obj( objfile.c_str() );
          output_obj << "#  OBJ format" << std::endl;
          output_obj << "# DGtal::saveOBJ" << std::endl;
          output_obj << "o vectors" << std::endl;
          output_obj << "mtllib " << mtlfile << std::endl;
          std::ofstream output_mtl( mtlfile.c_str() );
          output_mtl << "#  MTL format"<< std::endl;
          output_mtl << "# generated from MeshWriter from the DGTal library"<< std::endl;
          // Output vertices
          unsigned int n = std::min( positions.size(), vf.size() );
          for ( unsigned int i = 0; i < n; ++i )
            {
              RealPoint    p0 = positions[ i ];
              RealPoint    p1 = p0 + vf[ i ];
              RealVector    v = vf[ i ];
              RealVector absv = RealVector( fabs( v[ 0 ] ), fabs( v[ 1 ] ), fabs( v[ 2 ] ) );
              unsigned int mc = std::max_element( absv.begin(), absv.end() ) - absv.begin();
              RealVector   b  =
                mc == 2 ? RealVector( 1, 0, 0 ) :
                mc == 1 ? RealVector( 0, 0, 1 ) : RealVector( 0, 1, 0 );
              RealVector   e0 = v.crossProduct( b ).getNormalized();
              RealVector   e1 = v.crossProduct( e0 ).getNormalized();
              RealPoint    t[4] = {  thickness * e0,  thickness * e1,
                                    -thickness * e0, -thickness * e1 };
              for ( unsigned int j = 0; j < 4; ++j ) {
                RealPoint pt0 = p0 + t[ j ];
                RealPoint pt1 = p1 + t[ j ];
                output_obj << "v " << pt0[ 0 ] << " " << pt0[ 1 ] << " " << pt0[ 2 ]
                           << std::endl;
                output_obj << "v " << pt1[ 0 ] << " " << pt1[ 1 ] << " " << pt1[ 2 ]
                           << std::endl;
              }
            }
          // Simplify materials (very useful for blender).
          std::map<Color,Idx> map_colors;
          {
            Idx j = 0;
            for ( auto && c : diffuse_colors )
              if ( ! map_colors.count( c ) )
                map_colors[ c ] = j++;
          }
          
          // Output materials
          bool has_material = ! diffuse_colors.empty();
          if ( has_material )
            for ( auto&& pair : map_colors )
              MeshHelpers::exportMTLNewMaterial
                ( output_mtl, pair.second, ambient_color, pair.first, specular_color);
          else
            MeshHelpers::exportMTLNewMaterial
              ( output_mtl, 0, ambient_color, diffuse_color, specular_color );
          // Output faces
          for ( Idx i = 0; i < n; ++i )
            {
              output_obj << "usemtl material_" // << ( has_material ? i : 0 )
                         << ( has_material ? map_colors[ diffuse_colors[ i ] ] : 0 )
                         << std::endl;
              Idx b = 8*i+1;
              for ( Idx j = 0; j < 8; j += 2 )
                output_obj << "f " << (b+j) << " " << (b+j+1)
                           << " " << (b+(j+3)%8) << " " << (b+(j+2)%8) << std::endl;
              output_obj << "f " << b << " " << (b+2)
                         << " " << (b+4) << " " << (b+6) << std::endl;
              output_obj << "f " << (b+1) << " " << (b+7)
                         << " " << (b+5) << " " << (b+3) << std::endl;
            }
          output_mtl.close();
          return output_obj.good();
        }
      
      // ----------------------- Mesh services ------------------------------
    public:
 
      static Parameters parametersMesh()
      {
        return Parameters
          ( "faceSubdivision", "Centroid" ); 
      }
      
      static CountedPtr< TriangulatedSurface >
        makeTriangulatedSurface( CountedPtr< Mesh > aMesh )
      {
        auto pTriSurf = CountedPtr<TriangulatedSurface>
          ( new TriangulatedSurface ); // acquired
        bool ok = MeshHelpers::mesh2TriangulatedSurface( *aMesh, *pTriSurf );
        return ok ? pTriSurf : CountedPtr< TriangulatedSurface >( nullptr );
      }
 
      static CountedPtr< Mesh >
        makeMesh( CountedPtr< TriangulatedSurface > triSurf,
                  const Color& aColor = Color::White )
      {
        auto pMesh = CountedPtr<Mesh>( new Mesh( aColor ) ); // acquired
        MeshHelpers::triangulatedSurface2Mesh( *triSurf, *pMesh );
        return pMesh;
      }
 
      static CountedPtr< Mesh >
        makeMesh( CountedPtr< PolygonalSurface > polySurf,
                  const Color& aColor = Color::White )
      {
        auto pMesh = CountedPtr<Mesh>( new Mesh( aColor ) ); // acquired
        MeshHelpers::polygonalSurface2Mesh( *polySurf, *pMesh );
        return pMesh;
      }
 
      template < typename TContainer >
        static CountedPtr< TriangulatedSurface >
        makeTriangulatedSurface( Surfel2Index& s2i,
                                 CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface )
        {
          auto embedder = getCellEmbedder( aSurface );
          auto pTriSurf = CountedPtr<TriangulatedSurface>
            ( new TriangulatedSurface ); // acquired
          MeshHelpers::digitalSurface2DualTriangulatedSurface
            ( *aSurface, embedder, *pTriSurf, s2i );
          return pTriSurf;
        }
 
      template < typename TContainer >
        static CountedPtr< TriangulatedSurface >
        makeTriangulatedSurface( CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface )
        {
          Surfel2Index s2i;
          return makeTriangulatedSurface( s2i, aSurface );
        }
 
      static CountedPtr< TriangulatedSurface >
        makeTriangulatedSurface( CountedPtr< PolygonalSurface > polySurf,
                                 const Parameters&   params = parametersMesh() )
      {
        std::string faceSubdivision = params[ "faceSubdivision" ].as<std::string>();
        bool centroid = ( faceSubdivision == "Centroid" );
        auto pTriSurf = CountedPtr<TriangulatedSurface>
          ( new TriangulatedSurface ); // acquired
        MeshHelpers::polygonalSurface2TriangulatedSurface( *polySurf, *pTriSurf, centroid );
        return pTriSurf;
      }
    
      static CountedPtr< PolygonalSurface >
        makePolygonalSurface( CountedPtr< Mesh > aMesh )
      {
        auto pPolySurf = CountedPtr<PolygonalSurface>
          ( new PolygonalSurface ); // acquired
        bool ok = MeshHelpers::mesh2PolygonalSurface( *aMesh, *pPolySurf );
        return ok ? pPolySurf : CountedPtr< PolygonalSurface >( nullptr );
      }
 
      static CountedPtr< PolygonalSurface >
        makePolygonalSurface( CountedPtr<GrayScaleImage> gray_scale_image,
                              const Parameters&          params =
                              parametersKSpace()
                              | parametersBinaryImage()
                              | parametersDigitalSurface() )
      {
        auto K       = getKSpace( gray_scale_image );
        auto bimage  = makeBinaryImage( gray_scale_image, params );
        auto digSurf = makeDigitalSurface( bimage, K, params );
        RealVector gh = { params[ "gridsizex" ].as<double>(),
                          params[ "gridsizey" ].as<double>(),
                          params[ "gridsizez" ].as<double>() };
        double threshold = params[ "thresholdMin" ].as<double>() + 0.5;
        typedef RegularPointEmbedder<Space>         PointEmbedder;
        typedef ImageLinearCellEmbedder
          < KSpace, GrayScaleImage, PointEmbedder > ImageCellEmbedder;
        PointEmbedder pembedder;
        pembedder.init( gh );
        ImageCellEmbedder cembedder;
        cembedder.init( K, *gray_scale_image, pembedder, threshold );
        auto pPolySurf = CountedPtr<PolygonalSurface>
          ( new PolygonalSurface ); // acquired
        Surfel2Index s2i;
        MeshHelpers::digitalSurface2DualPolygonalSurface
          ( *digSurf, cembedder, *pPolySurf, s2i );
        return pPolySurf;
      }
 
      static CountedPtr< TriangulatedSurface >
        makeTriangulatedSurface( CountedPtr<GrayScaleImage> gray_scale_image,
                                 const Parameters&          params =
                                 parametersKSpace()
                                 | parametersBinaryImage()
                                 | parametersDigitalSurface() )
      {
        auto K       = getKSpace( gray_scale_image );
        auto bimage  = makeBinaryImage( gray_scale_image, params );
        auto digSurf = makeDigitalSurface( bimage, K, params );
        RealVector gh = { params[ "gridsizex" ].as<double>(),
                          params[ "gridsizey" ].as<double>(),
                          params[ "gridsizez" ].as<double>() };
        double threshold = params[ "thresholdMin" ].as<double>() + 0.5;
        typedef RegularPointEmbedder<Space>         PointEmbedder;
        typedef ImageLinearCellEmbedder
          < KSpace, GrayScaleImage, PointEmbedder > ImageCellEmbedder;
        PointEmbedder pembedder;
        pembedder.init( gh );
        ImageCellEmbedder cembedder;
        cembedder.init( K, *gray_scale_image, pembedder, threshold );
        auto pPolySurf = CountedPtr<TriangulatedSurface>
          ( new TriangulatedSurface ); // acquired
        Surfel2Index s2i;
        MeshHelpers::digitalSurface2DualTriangulatedSurface
          ( *digSurf, cembedder, *pPolySurf, s2i );
        return pPolySurf;
      }
 
      template < typename TContainer >
        static CountedPtr< PolygonalSurface >
        makeDualPolygonalSurface( Surfel2Index& s2i,
                                  CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface )
        {
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          auto embedder = getCellEmbedder( aSurface );
          auto pPolySurf = CountedPtr<PolygonalSurface>
            ( new PolygonalSurface ); // acquired
          MeshHelpers::digitalSurface2DualPolygonalSurface
            ( *aSurface, embedder, *pPolySurf, s2i );
          return pPolySurf;
        }
 
      template < typename TContainer >
        static CountedPtr< PolygonalSurface >
        makeDualPolygonalSurface( CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface )
        {
          Surfel2Index s2i;
          return makeDualPolygonalSurface( s2i, aSurface );
        }
 
      template < typename TContainer >
        static CountedPtr< PolygonalSurface >
        makeDualPolygonalSurface
        ( CountedPtr< ::DGtal::IndexedDigitalSurface< TContainer > > aSurface )
        {
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          auto pPolySurf = CountedPtr<PolygonalSurface>
            ( new PolygonalSurface( aSurface->heds(),
                                    aSurface->positions().storage() ) );
          return pPolySurf;
        }
 
      template < typename TContainer >
        static CountedPtr< PolygonalSurface >
        makePrimalPolygonalSurface( Cell2Index& c2i,
                                    CountedPtr< ::DGtal::DigitalSurface<TContainer> > aSurface )
        {
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          auto embedder = getCellEmbedder( aSurface );
          auto pPolySurf = CountedPtr<PolygonalSurface>
            ( new PolygonalSurface ); // acquired
          bool ok = MeshHelpers::digitalSurface2PrimalPolygonalSurface
            ( *aSurface, embedder, *pPolySurf, c2i );
          return ok ? pPolySurf : CountedPtr< PolygonalSurface >( nullptr );
        }
 
      template < typename TContainer >
        static CountedPtr< PolygonalSurface >
        makePrimalPolygonalSurface( CountedPtr< ::DGtal::DigitalSurface<TContainer> > aSurface )
        {
          Cell2Index c2i;
          return makePrimalPolygonalSurface( c2i, aSurface );
        }
 
      template < typename TContainer >
        static CountedPtr< PolygonalSurface >
        makePrimalPolygonalSurface
        ( CountedPtr< ::DGtal::IndexedDigitalSurface<TContainer> > aSurface )
        {
          auto dsurf = makeDigitalSurface( aSurface );
          Cell2Index c2i;
          return makePrimalPolygonalSurface( c2i, dsurf );
        }
 
      template <typename TPoint>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::PolygonalSurface<TPoint> > polysurf,
          const std::string&                            objfile )
        {
          std::ofstream output( objfile.c_str() );
          bool ok = MeshHelpers::exportOBJ( output, *polysurf );
          output.close();
          return ok;
        }
 
      template <typename TPoint>
      static bool
      saveOFF ( CountedPtr< ::DGtal::PolygonalSurface<TPoint> > polysurf,
                  std::string                                   off_file,
                  const Color&                                 face_color  = DGtal::Color::None)
          
        {
          DGtal::Mesh< TPoint > m;
          MeshHelpers::polygonalSurface2Mesh(*polysurf,m);
          std::ofstream output( off_file.c_str() );
          for (unsigned int i=0; i< m.nbFaces(); i++)
          {
            m.setFaceColor(i, face_color);
          }
          bool ok = MeshWriter< TPoint >::export2OFF( output, m, true );
          output.close();
          return ok;
        }
 
      template <typename TPoint>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::TriangulatedSurface<TPoint> > trisurf,
          const std::string&                               objfile )
        {
          std::ofstream output( objfile.c_str() );
          bool ok = MeshHelpers::exportOBJ( output, *trisurf );
          output.close();
          return ok;
        }
 
      template <typename TPoint>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::PolygonalSurface<TPoint> > polysurf,
          const RealVectors&                            normals,
          const Colors&                                 diffuse_colors,
          std::string                                   objfile,
          const Color&                   ambient_color  = Color( 32, 32, 32 ),
          const Color&                   diffuse_color  = Color( 200, 200, 255 ),
          const Color&                   specular_color = Color::White )
        {
          std::string mtlfile;
          auto lastindex = objfile.find_last_of(".");
          if ( lastindex == std::string::npos )
            {
              mtlfile  = objfile + ".mtl";
              objfile  = objfile + ".obj";
            }
          else
            {
              mtlfile  = objfile.substr(0, lastindex) + ".mtl"; 
            }
          std::ofstream output( objfile.c_str() );
          bool ok = MeshHelpers::exportOBJwithFaceNormalAndColor
            ( output, mtlfile, *polysurf, normals, diffuse_colors,
              ambient_color, diffuse_color, specular_color );
          output.close();
          return ok;
        }
 
      template <typename TPoint>
        static bool
        saveOBJ
        ( CountedPtr< ::DGtal::TriangulatedSurface<TPoint> > trisurf,
          const RealVectors&                            normals,
          const Colors&                                 diffuse_colors,
          std::string                                   objfile,
          const Color&                   ambient_color  = Color( 32, 32, 32 ),
          const Color&                   diffuse_color  = Color( 200, 200, 255 ),
          const Color&                   specular_color = Color::White )
        {
          std::string mtlfile;
          auto lastindex = objfile.find_last_of(".");
          if ( lastindex == std::string::npos )
            {
              mtlfile  = objfile + ".mtl";
              objfile  = objfile + ".obj";
            }
          else
            {
              mtlfile  = objfile.substr(0, lastindex) + ".mtl"; 
            }
          std::ofstream output( objfile.c_str() );
          bool ok = MeshHelpers::exportOBJwithFaceNormalAndColor
            ( output, mtlfile, *trisurf, normals, diffuse_colors,
              ambient_color, diffuse_color, specular_color );
          output.close();
          return ok;
        }
 
      template <typename TPoint>
      static bool
      saveOFF ( CountedPtr< ::DGtal::TriangulatedSurface<TPoint> > trisurf,
                  std::string                                   off_file,
                  const Color&                   face_color  = DGtal::Color::None)
          
        {
          DGtal::Mesh< TPoint > m;
          MeshHelpers::triangulatedSurface2Mesh(*trisurf,m);
          std::ofstream output( off_file.c_str() );
          for (unsigned int i=0; i< m.nbFaces(); i++)
          {
            m.setFaceColor(i, face_color);
          }
          bool ok = MeshWriter< TPoint >::export2OFF( output, m, true );
          output.close();
          return ok;
        }
 
    
      // ------------------------------ utilities ------------------------------
    public:
 
      static Parameters parametersUtilities()
      {
        return Parameters
          ( "colormap", "Custom" )
          ( "zero-tic", 0.0 ); 
      }
 
      template <typename TValue>
        static IdxRange
        getRangeMatch( const std::vector< TValue >& s1, const std::vector< TValue >& s2,
                       bool perfect = false )
        {
          if ( perfect && ( s1.size() != s2.size() ) ) return IdxRange();
          std::map<TValue, Idx> M;
          Idx idx = 0;
          for ( auto val : s2 ) M[ val ] = idx++;
          IdxRange V( s1.size() );
          idx = 0;
          for ( auto val : s1 )
            {
              auto it = M.find( val );
              if ( it != M.end() ) V[ idx++ ] = it->second;
              else
                {
                  if ( perfect ) return IdxRange();
                  V[ idx++ ] = s2.size();
                }
            }
          return V;
        }
 
      template <typename TValue>
        static std::vector< TValue >
        getMatchedRange( const std::vector< TValue >& range, const IdxRange& match )
        {
          std::vector< TValue > result( match.size() );
          for ( Idx i = 0; i < result.size(); i++ )
            result[ i ] = range[ match[ i ] ]; 
          return result;
        }
    
    
      static ColorMap
        getColorMap( Scalar            min,
                     Scalar            max,
                     const Parameters& params = parametersUtilities() )
      {
        std::string cmap = params[ "colormap" ].as<std::string>();
        if      ( cmap == "Cool" )   return ColorMap( min, max, CMAP_COOL );
        else if ( cmap == "Copper" ) return ColorMap( min, max, CMAP_COPPER );
        else if ( cmap == "Hot" )    return ColorMap( min, max, CMAP_HOT );
        else if ( cmap == "Jet" )    return ColorMap( min, max, CMAP_JET );
        else if ( cmap == "Spring" ) return ColorMap( min, max, CMAP_SPRING );
        else if ( cmap == "Summer" ) return ColorMap( min, max, CMAP_SUMMER );
        else if ( cmap == "Autumn" ) return ColorMap( min, max, CMAP_AUTUMN );
        else if ( cmap == "Winter" ) return ColorMap( min, max, CMAP_WINTER );
        else if ( cmap == "Error" )
          {
            ColorMap gradcmap( min, max );
            gradcmap.addColor( Color( 255, 255, 255 ) );
            gradcmap.addColor( Color( 255,   0,   0 ) );
            gradcmap.addColor( Color( 0,   0,   0 ) );
            return gradcmap;
          }
        // Custom colormap 
        ColorMap gradcmap( min, max );
        gradcmap.addColor( Color( 0, 0, 255 ) );
        gradcmap.addColor( Color( 0, 255, 255 ) );
        gradcmap.addColor( Color( 255, 255, 255 ) );
        gradcmap.addColor( Color( 255, 255, 0 ) );
        gradcmap.addColor( Color( 255, 0, 0 ) );
        return gradcmap;
      }
 
      static ZeroTickedColorMap
        getZeroTickedColorMap( Scalar           min,
                               Scalar           max,
                               const Parameters& params = parametersUtilities() )
      {
        auto cmap = getColorMap( min, max, params );
        auto ztic = params[ "zero-tic" ].as<double>();
        ZeroTickedColorMap ztic_cmap( cmap, Color::Black );
        if ( ztic <= 0.0 ) return ztic_cmap;
        if ( min <= 0.0 && 0.0 <= max )
          ztic_cmap.addTick( 0.0, ztic );
        ztic_cmap.finalize();
        return ztic_cmap;
      }
 
    
      template <typename TCellEmbedder = CanonicCellEmbedder< KSpace > >
        static bool
        outputSurfelsAsObj
        ( std::ostream&              output,
          const SurfelRange&         surfels,
          const TCellEmbedder&       embedder )
        {
          typedef unsigned long Size;
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          BOOST_CONCEPT_ASSERT(( concepts::CCellEmbedder< TCellEmbedder > ));
          const KSpace& K = embedder.space();
          // Number and output vertices.
          std::map< Cell, Size > vtx_numbering;
          Size n = 1;  // OBJ vertex numbering start at 1 
          for ( auto&& s : surfels )
            {
              auto primal_vtcs = getPointelRange( K, s, true );
              for ( auto&& primal_vtx : primal_vtcs )
                {
                  if ( ! vtx_numbering.count( primal_vtx ) )
                    {
                      vtx_numbering[ primal_vtx ] = n++;
                      // Output vertex positions
                      RealPoint p = embedder( primal_vtx );
                      output << "v " << p[ 0 ] << " " << p[ 1 ] << " " << p[ 2 ] << std::endl;
                    }
                }
            }
          // Outputs all faces
          for ( auto&& s : surfels )
            {
              output << "f";
              auto primal_vtcs = getPointelRange( K, s, true );
              for ( auto&& primal_vtx : primal_vtcs )
                output << " " << vtx_numbering[ primal_vtx ];
              output << std::endl;
            }      
          return output.good();
        }
    
      template <typename TAnyDigitalSurface>
        static bool
        outputPrimalDigitalSurfaceAsObj
        ( std::ostream&              output,
          CountedPtr<TAnyDigitalSurface> surface )
        {
          auto embedder = getCellEmbedder( surface );
          return outputPrimalDigitalSurfaceAsObj( output, surface, embedder );
        }
    
      template < typename TAnyDigitalSurface,
        typename TCellEmbedder = CanonicCellEmbedder< KSpace > >
        static bool
        outputPrimalDigitalSurfaceAsObj
        ( std::ostream&              output,
          CountedPtr<TAnyDigitalSurface> surface,
          const TCellEmbedder&       embedder )
        {
          auto surfels = getSurfelRange( surface, Parameters( "Traversal", "Default" ) );
          return outputSurfelsAsObj( output, surfels, embedder );
        }
 
      static bool
        outputPrimalIdxDigitalSurfaceAsObj
        ( std::ostream&              output,
          CountedPtr<IdxDigitalSurface> surface )
      {
        auto embedder = getCellEmbedder( surface );
        return outputPrimalIdxDigitalSurfaceAsObj( output, surface, embedder );
      }
    
      template <typename TCellEmbedder = CanonicCellEmbedder< KSpace > >
        static bool
        outputPrimalIdxDigitalSurfaceAsObj
        ( std::ostream&                 output,
          CountedPtr<IdxDigitalSurface> surface,
          const TCellEmbedder&          embedder )
        {
          auto idxsurfels = getIdxSurfelRange( surface, Parameters( "Traversal", "Default" ) );
          auto  surfelmap = surface->surfels();
          SurfelRange surfels;
          for ( auto&& idx : idxsurfels )
            surfels.push_back( surfelmap[ idx ] );
          return outputSurfelsAsObj( output, surfels, embedder );
        }
 
      template <typename TDigitalSurfaceContainer>
        static bool
        outputDualDigitalSurfaceAsObj
        ( std::ostream&              output,
          CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer> > surface,
          const Parameters&   params = parametersMesh() )
        {
          auto embedder = getCellEmbedder( surface );
          return outputDualDigitalSurfaceAsObj( output, surface, embedder, params );
        }
    
      template < typename TDigitalSurfaceContainer,
        typename TCellEmbedder = CanonicCellEmbedder< KSpace > >
        static bool
        outputDualDigitalSurfaceAsObj
        ( std::ostream&              output,
          CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface,
          const TCellEmbedder&       embedder,
          const Parameters&   params = parametersMesh() )
        {
          typedef unsigned long Size;
          BOOST_STATIC_ASSERT (( KSpace::dimension == 3 ));
          BOOST_CONCEPT_ASSERT(( concepts::CCellEmbedder< TCellEmbedder > ));
          std::string dualFaceSubdivision = params[ "faceSubdivision" ].as<std::string>();
          const int   subdivide
          = dualFaceSubdivision == "Naive"    ? 1
          : dualFaceSubdivision == "Centroid" ? 2
          : 0;
          const KSpace& K = embedder.space();
          // Number and ouput vertices.
          std::map< Vertex, Size > vtx_numbering;
          std::map< Face,   Size > sub_numbering;
          Size n = 1;  // OBJ vertex numbering start at 1 
          for ( auto && s : *surface )
            {
              if ( ! vtx_numbering.count( s ) )
                {
                  vtx_numbering[ s ] = n++;
                  // Output vertex positions
                  RealPoint p = embedder( K.unsigns( s ) );
                  output << "v " << p[ 0 ] << " " << p[ 1 ] << " " << p[ 2 ] << std::endl;
                }
            }
          auto faces = surface->allClosedFaces();
          // Prepare centroids if necessary
          if ( subdivide == 2 )
            {
              for ( auto&& f : faces )
                {
                  auto vtcs = surface->verticesAroundFace( f );
                  Size   nv = vtcs.size();
                  if ( nv > 3 )
                    {
                      sub_numbering[ f ] = n++;
                      RealPoint p = RealPoint::zero;
                      for ( auto&& s : vtcs ) p += embedder( K.unsigns( s ) );
                      p /= nv;
                      output << "v " << p[ 0 ] << " " << p[ 1 ] << " " << p[ 2 ] << std::endl;
                    }
                }
            }
          // Outputs closed faces.
          if ( subdivide == 0 )
            { // No subdivision
              for ( auto&& f : faces )
                {
                  output << "f";
                  auto vtcs = surface->verticesAroundFace( f );
                  std::reverse( vtcs.begin(), vtcs.end() );
                  for ( auto&& s : vtcs )
                    output << " " << vtx_numbering[ s ];
                  output << std::endl;
                }
            }
          else if ( subdivide == 1 )
            { // naive subdivision
              for ( auto&& f : faces )
                {
                  auto vtcs = surface->verticesAroundFace( f );
                  Size   nv = vtcs.size();
                  for ( Size i = 1; i < nv - 1; ++i )
                    output << "f " << vtx_numbering[ vtcs[ 0 ] ]
                           << " "  << vtx_numbering[ vtcs[ i+1 ] ]
                           << " "  << vtx_numbering[ vtcs[ i ] ] << std::endl;
                }
            }
          else if ( subdivide == 2 )
            { // centroid subdivision
              for ( auto&& f : faces )
                {
                  auto vtcs = surface->verticesAroundFace( f );
                  Size   nv = vtcs.size();
                  if ( nv == 3 )
                    output << "f " << vtx_numbering[ vtcs[ 0 ] ]
                           << " "  << vtx_numbering[ vtcs[ 2 ] ]
                           << " "  << vtx_numbering[ vtcs[ 1 ] ] << std::endl;
                  else {
                    Size c = sub_numbering[ f ];
                    for ( Size i = 0; i < nv; ++i )
                      {
                        output << "f " << c
                               << " "  << vtx_numbering[ vtcs[ (i+1)%nv ] ]
                               << " "  << vtx_numbering[ vtcs[ i ] ] << std::endl;
                      }
                  }
                }
            }
          return output.good();
        }
    
 
      // -------------------- map I/O services ------------------------------------------
    public:
      struct CellWriter {
        void operator()( std::ostream& output, const KSpace& K, const Cell & cell )
        {
          for ( Dimension d = 0; d < KSpace::dimension; ++d )
            output << " " << K.sKCoord( cell, d );
        }
      };
      struct CellReader {
        Cell operator()( std::istream& input, const KSpace& K )
        {
          Point kp;
          for ( Dimension d = 0; d < KSpace::dimension; ++d )
            input >> kp[ d ];
          return K.uCell( kp );
        }
      };
      struct SCellWriter {
        void operator()( std::ostream& output, const KSpace& K, const SCell & scell )
        {
          CellWriter::operator()( output, K, K.unsigns( scell ) );
          output << " " << K.sSign( scell );
        }
      };
      struct SCellReader {
        SCell operator()( std::istream& input, const KSpace& K )
        {
          Point                 kp;
          typename KSpace::Sign s;
          for ( Dimension d = 0; d < KSpace::dimension; ++d )
            input >> kp[ d ];
          input >> s;
          return K.sCell( kp, s );
        }
      };
 
      template <typename Value>
        struct ValueWriter {
          void operator()( std::ostream& output, const Value& v )
          {
            output << " " << v;
          }
          void operator()( std::ostream& output, const std::vector<Value>& vv )
          {
            for ( auto&& v : vv ) output << " " << v;
          }
        };
 
      template <typename Value>
        struct ValueReader {
          bool operator()( std::istream& input, Value& value )
          {
            std::string str;
            std::getline( input, str );
            // construct a stream from the string
            std::stringstream strstr(str);
            // use stream iterators to copy the stream to the vector as whitespace separated strings
            std::istream_iterator<std::string> it(strstr);
            std::istream_iterator<std::string> end;
            std::vector<std::string> results(it, end);
            std::stringstream sstr( results[ 0 ] );
            sstr >> value;
            return ( results.size() == 1 ) && input.good();
          }
      
          bool operator()( std::istream& input, std::vector<Value>& values )
          {
            std::string str;
            std::getline( input, str );
            // construct a stream from the string
            std::stringstream strstr(str);
            // use stream iterators to copy the stream to the vector as whitespace separated strings
            std::istream_iterator<std::string> it(strstr);
            std::istream_iterator<std::string> end;
            std::vector<std::string> results(it, end);
            values.resize( results.size() );
            for ( unsigned int i = 0; i < results.size(); ++i ) {
              std::stringstream sstr( results[ i ] );
              sstr >> values[ i ];
            }
            return input.good();
          }
        };
 
      // Outputs in \a output a map \a anyMap: SCell -> Value given the
      // appropriate value \a writer.
      //
      // @tparam TSCellMap any model of map SCell -> Value., e.g. std::map<SCell,double>
      // @tparam TValueWriter any model of value writer, e.g. ValueWriter<double>
      //
      // @param[out] output the output stream
      // @param[in]  K the Khalimsky space where cells are defined.
      // @param[in]  anyMap the map associated a value to signed cells.
      // @param[in]  writer the writer that can write values on the ouput
      // stream, e.g. ValueWriter<double> to write double value or
      // vector<double> values.
      template <typename TSCellMap, typename TValueWriter>
        static
        bool outputSCellMapAsCSV
        ( std::ostream&       output,
          const KSpace&       K,
          const TSCellMap&     anyMap,
          const TValueWriter& writer )
      {
        SCellWriter w;
        for ( auto&& v : anyMap )
          {
            w( output, K, v.first );
            writer( output, v.second );
            output << std::endl;
          }
        return output.good();
      }
    
      // Outputs in \a output a map \a anyMap: SCell -> value given the
      // appropriate value \a writer.
      //
      // @tparam TCellMap any model of map Cell -> Value., e.g. std::map<Cell,double>
      // @tparam TValueWriter any model of value writer, e.g. ValueWriter<double>
      //
      // @param[out] output the output stream
      // @param[in]  K the Khalimsky space where cells are defined.
      // @param[in]  anyMap the map associated a value to signed cells.
      // @param[in]  writer the writer that can write values on the ouput
      // stream, e.g. ValueWriter<double> to write double value or
      // vector<double> values.
      template <typename TCellMap, typename TValueWriter>
        static
        bool outputCellMapAsCSV
        ( std::ostream&       output,
          const KSpace&       K,
          const TCellMap&      anyMap,
          const TValueWriter& writer )
      {
        CellWriter w;
        for ( auto&& v : anyMap )
          {
            w( output, K, v.first );
            writer( output, v.second );
            output << std::endl;
          }
        return output.good();
      }
    
      static
      CellRange getPrimalCells( const KSpace& K, const SCell& s, const Dimension k )
      {
        auto faces = K.uFaces( K.unsigns( s ) );
        CellRange primal_cells;
        for ( auto&& f : faces )
          {
            if ( K.uDim( f ) == k ) primal_cells.push_back( f );
          }
        return primal_cells;
      }
 
      static
      CellRange getPrimalVertices( const KSpace& K, const SCell& s )
      {
        return getPrimalCells( K, s, 0 );
      }
    
      static
      CellRange getPrimalVertices( const KSpace& K, const Surfel& s, bool ccw )
      {
        BOOST_STATIC_ASSERT(( KSpace::dimension == 3 ));
        CellRange vtcs = getPrimalVertices( K, s );
        std::swap( vtcs[ 2 ], vtcs[ 3 ] );
        auto  orth_dir = K.sOrthDir( s );
        auto    direct = K.sDirect( s, orth_dir ) ? ccw : ! ccw;
        Vector    s0s1 = K.uCoords( vtcs[ 1 ] ) - K.uCoords( vtcs[ 0 ] );
        Vector    s0s2 = K.uCoords( vtcs[ 2 ] ) - K.uCoords( vtcs[ 0 ] );
        Vector       t = s0s1.crossProduct( s0s2 );
        if ( ( ( t[ orth_dir ] > 0.0 ) && direct )
             || ( ( t[ orth_dir ] < 0.0 ) && ! direct ) )
          std::reverse( vtcs.begin(), vtcs.end() );
        return vtcs;
      }
 
    
      // ----------------------- Standard services ------------------------------
    public:
 
      Shortcuts() = delete;
 
      ~Shortcuts() = delete;
 
      Shortcuts ( const Shortcuts & other ) = delete;
 
      Shortcuts ( Shortcuts && other ) = delete;
 
      Shortcuts & operator= ( const Shortcuts & other ) = delete;
 
      Shortcuts & operator= ( Shortcuts && other ) = delete;
 
      // ----------------------- Interface --------------------------------------
    public:
 
      void selfDisplay ( std::ostream & out ) const
      {
        out << "[Shortcuts]";
      }
 
      bool isValid() const
      {
        return true;
      }
 
      // ------------------------- Protected Datas ------------------------------
    protected:
 
      // ------------------------- Private Datas --------------------------------
    private:
 
      // ------------------------- Hidden services ------------------------------
    protected:
 
      // ------------------------- Internals ------------------------------------
    private:
 
    }; // end of class Shortcuts
 
 
  template <typename T>
    std::ostream&
    operator<< ( std::ostream & out, const Shortcuts<T> & object );
 
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined Shortcuts_h
 
#undef Shortcuts_RECURSES
#endif // else defined(Shortcuts_RECURSES)