DECImageHelpers.h

 
#pragma once
 
#if defined(DECImageHelpers_RECURSES)
#error Recursive header files inclusion detected in DECImageHelpers.h
#else // defined(DECImageHelpers_RECURSES)
#define DECImageHelpers_RECURSES
 
#if !defined DECImageHelpers_h
#define DECImageHelpers_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/images/CImage.h"
#include "DGtal/topology/CCellularGridSpaceND.h"
#include "DGtal/math/linalg/CLinearAlgebra.h"
#include "DGtal/dec/DiscreteExteriorCalculus.h"
 
 
namespace DGtal {
  namespace functions {
    namespace dec {
      
      template <typename Calculus, DGtal::Dimension dim, DGtal::Duality duality> 
      DGtal::LinearOperator<Calculus, dim, duality, dim, duality> 
      diagonal(const DGtal::KForm<Calculus, dim, duality>& kform) 
      { 
        typedef DGtal::LinearOperator<Calculus,dim, duality, dim, duality> Operator; 
        typedef typename Calculus::LinearAlgebraBackend::Triplet           Triplet; 
        typedef typename Calculus::Index                                   Index; 
        typedef std::vector<Triplet>                                       Triplets;
 
        Triplets triplets;
        for (Index index=0; index<kform.length(); index++) 
          triplets.push_back(Triplet(index, index, kform.myContainer(index)));
        
        Operator op( kform.myCalculus );
        op.myContainer.setFromTriplets( triplets.begin(), triplets.end() );
        return op;
      }
 
      template <typename Calculus, DGtal::Dimension dim, DGtal::Duality duality> 
      void
      squares(DGtal::KForm<Calculus, dim, duality>& kform) 
      {
        kform.myContainer.array() = kform.myContainer.array().square();
      }
      
      template <typename Calculus, DGtal::Dimension dim, DGtal::Duality duality> 
      DGtal::LinearOperator<Calculus, dim, duality, dim, duality> 
      squaredDiagonal(const DGtal::KForm<Calculus, dim, duality>& kform) 
      {
        auto v2 = kform;
        squares( v2 );
        return diagonal( v2 );
      }
 
      template <typename Image>
      void writePixel( Image& image, typename Image::Point pt, typename Image::Value val,
                       int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        typedef typename Image::Point      Point;
        typedef typename Point::Coordinate Coordinate;
        pt *= pixel_size;
        for ( int y = 0; y < pixel_size; y++ )
          for ( int x = 0; x < pixel_size; x++ )
            {
              Point q( (Coordinate) x, (Coordinate) y );
              image.setValue( pt + q, val );
            }
      }
 
      template <typename Image>
      void writePrimalLinel( Image& image, typename Image::Point pt, typename Image::Value val,
                       int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        typedef typename Image::Point      Point;
        typedef typename Point::Coordinate Coordinate;
        int pixel_size_x = NumberTraits<Coordinate>::even( pt[ 0 ] ) ? 1 : pixel_size;
        int pixel_size_y = NumberTraits<Coordinate>::even( pt[ 1 ] ) ? 1 : pixel_size;
        pt /= 2;
        pt *= pixel_size;
        for ( int y = 0; y < pixel_size_y; y++ )
          for ( int x = 0; x < pixel_size_x; x++ )
            {
              Point q( (Coordinate) x, (Coordinate) y );
              image.setValue( pt + q, val );
            }
      }
      
      template <typename Image>
      void writeDualLinel( Image& image, typename Image::Point pt, typename Image::Value val,
                           int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        typedef typename Image::Point      Point;
        typedef typename Point::Coordinate Coordinate;
        int pixel_size_x = NumberTraits<Coordinate>::even( pt[ 0 ] ) ? 0 : pixel_size-1;
        int pixel_size_y = NumberTraits<Coordinate>::even( pt[ 1 ] ) ? 0 : pixel_size-1;
        pt /= 2;
        pt *= pixel_size;
        for ( int y = pixel_size_y; y < pixel_size; y++ )
          for ( int x = pixel_size_x; x < pixel_size; x++ )
            {
              Point q( (Coordinate) x, (Coordinate) y );
              image.setValue( pt + q, val );
            }
      }
      
      template <typename Calculus, typename AnyForm2, typename Image, typename Function>
      void form2ToImage
      ( const Calculus& calculus, 
        const AnyForm2& u, 
        Image& image,
        const Function& functor,
        double cut_low = 0.0, double cut_up = 1.0, int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        typedef typename Calculus::Index  Index;
        typedef typename Calculus::SCell  SCell;
        typedef typename Calculus::Scalar Scalar;
        typedef typename Calculus::KSpace KSpace;
        typedef typename KSpace::Point    Point;
        typedef typename KSpace::Integer  Integer;
        double min_u = NumberTraits<Scalar>::castToDouble( u.myContainer[ 0 ] );
        double max_u = min_u;
        for ( Index index = 0; index < u.myContainer.rows(); index++)
          {
            double v = NumberTraits<Scalar>::castToDouble( u.myContainer[ index ] );
            min_u = std::min( min_u, v );
            max_u = std::max( max_u, v );
          }
        if ( min_u < cut_low ) min_u = cut_low;
        if ( max_u > cut_up  ) max_u = cut_up;
        for ( Index index = 0; index < u.myContainer.rows(); index++)
          {
            SCell cell = u.getSCell( index );
            double v = NumberTraits<Scalar>::castToDouble( u.myContainer[ index ] );
            double w = std::min( cut_up, std::max( cut_low, v ) );
            if ( min_u != max_u ) w = ( w - min_u ) / ( max_u - min_u );
            writePixel( image, calculus.myKSpace.sCoords( cell ), functor( w ), pixel_size );
          }
      }
 
      template <typename Calculus, typename Form1, typename Image, typename Function, typename Predicate>
      void form1ToImage
      ( const Calculus& calculus, 
        const Form1& v, bool primal,
        Image& image,
        const Function& functor,
        const Predicate& predicate,
        double cut_low = 0.0, double cut_up = 1.0, int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        typedef typename Calculus::Index  Index;
        typedef typename Calculus::SCell  SCell;
        typedef typename Calculus::Scalar Scalar;
        typedef typename Calculus::KSpace KSpace;
        typedef typename KSpace::Point    Point;
        typedef typename KSpace::Integer  Integer;
        double min_v = NumberTraits<Scalar>::castToDouble( v.myContainer[ 0 ] );
        double max_v = min_v;
        for ( Index index = 0; index < v.myContainer.rows(); index++)
          {
            double w = NumberTraits<Scalar>::castToDouble( v.myContainer[ index ] );
            min_v = std::min( min_v, w );
            max_v = std::max( max_v, w );
          }
        if ( min_v < cut_low ) min_v = cut_low;
        if ( max_v > cut_up  ) max_v = cut_up;
        for ( Index index = 0; index < v.myContainer.rows(); index++)
          {
            SCell cell = v.getSCell( index );
            double u = NumberTraits<Scalar>::castToDouble( v.myContainer[ index ] );
            if ( ! predicate( u ) ) continue; 
            double w = std::min( cut_up, std::max( cut_low, u ) );
            if ( min_v != max_v ) w = ( w - min_v ) / ( max_v - min_v );
            Point kpt = calculus.myKSpace.sKCoords( cell );
            if ( primal ) writePrimalLinel( image, kpt, functor( w ), pixel_size );
            else          writeDualLinel  ( image, kpt, functor( w ), pixel_size );
          }
      }
      
      template <typename Calculus, typename Image, typename Function, typename Predicate>
      void dualForm1ToImage
      ( const Calculus& calculus, 
        const typename Calculus::DualForm1& v, 
        Image& image,
        const Function& functor,
        const Predicate& predicate,
        double cut_low = 0.0, double cut_up = 1.0, int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        form1ToImage( calculus, v, false, image, functor, predicate,
                      cut_low, cut_up, pixel_size );
      }
      
      template <typename Calculus, typename Image, typename Function, typename Predicate>
      void primalForm1ToImage
      ( const Calculus& calculus, 
        const typename Calculus::PrimalForm1& v, 
        Image& image,
        const Function& functor,
        const Predicate& predicate,
        double cut_low = 0.0, double cut_up = 1.0, int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        form1ToImage( calculus, v, true, image, functor, predicate,
                      cut_low, cut_up, pixel_size );
      }
 
      template <typename Calculus, typename AnyForm2, typename Image, typename Function>
      void threeForms2ToImage
      ( const Calculus& calculus, 
        const AnyForm2& u0, 
        const AnyForm2& u1, 
        const AnyForm2& u2, 
        Image& image,
        const Function& functor,
        double cut_low = 0.0, double cut_up = 1.0, int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        typedef typename Calculus::Index  Index;
        typedef typename Calculus::SCell  SCell;
        typedef typename Calculus::Scalar Scalar;
        double min_u = NumberTraits<Scalar>::castToDouble( u0.myContainer[ 0 ] );
        double max_u = min_u;
        for ( Index index = 0; index < u0.myContainer.rows(); index++)
          {
            double v = NumberTraits<Scalar>::castToDouble( u0.myContainer[ index ] );
            min_u = std::min( min_u, v );
            max_u = std::max( max_u, v );
          }
        for ( Index index = 0; index < u1.myContainer.rows(); index++)
          {
            double v = NumberTraits<Scalar>::castToDouble( u1.myContainer[ index ] );
            min_u = std::min( min_u, v );
            max_u = std::max( max_u, v );
          }
        for ( Index index = 0; index < u2.myContainer.rows(); index++)
          {
            double v = NumberTraits<Scalar>::castToDouble( u2.myContainer[ index ] );
            min_u = std::min( min_u, v );
            max_u = std::max( max_u, v );
          }
        if ( min_u < cut_low ) min_u = cut_low;
        if ( max_u > cut_up  ) max_u = cut_up;
        for ( Index index = 0; index < u0.myContainer.rows(); index++)
          {
            SCell cell = u0.getSCell( index );
            double v0 = NumberTraits<Scalar>::castToDouble( u0.myContainer[ index ] );
            double w0 = std::min( cut_up, std::max( cut_low, v0 ) );
            if ( min_u != max_u ) w0 = ( w0 - min_u ) / ( max_u - min_u );
            double v1 = NumberTraits<Scalar>::castToDouble( u1.myContainer[ index ] );
            double w1 = std::min( cut_up, std::max( cut_low, v1 ) );
            if ( min_u != max_u ) w1 = ( w1 - min_u ) / ( max_u - min_u );
            double v2 = NumberTraits<Scalar>::castToDouble( u2.myContainer[ index ] );
            double w2 = std::min( cut_up, std::max( cut_low, v2 ) );
            if ( min_u != max_u ) w2 = ( w2 - min_u ) / ( max_u - min_u );
            writePixel( image, calculus.myKSpace.sCoords( cell ), functor( w0, w1, w2 ), pixel_size );
          }
      }
      
      template <typename Calculus, typename AnyForm2, typename Image>
      void form2ToGreyLevelImage
      ( const Calculus& calculus, 
        const AnyForm2& u, 
        Image& image,
        double cut_low = 0.0, double cut_up = 1.0,
        int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        form2ToImage( calculus, u, image,
                      [] ( double x ) { return (unsigned char) ( round( x * 255.0 ) ); },
                      cut_low, cut_up, pixel_size );
      }
 
      template <typename Calculus, typename Image>
      void primalForm1ToGreyLevelImage
      ( const Calculus& calculus, 
        const typename Calculus::PrimalForm1& v, 
        Image& image,
        double cut_low = 0.0, double cut_up = 1.0,
        int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        // Threshold is 0.25 instead of 0.5 because an edge connecting
        // two vertices with v=0 and v=1 should not belong to the
        // discontinuity set.
        primalForm1ToImage( calculus, v, image,
                            [] ( double x ) { return (unsigned char) ( round( x * 255.0 ) ); },
                            [] ( double x ) { return x < 0.25; },
                            cut_low, cut_up, pixel_size );
      }
      
      template <typename Calculus, typename Image>
      void dualForm1ToGreyLevelImage
      ( const Calculus& calculus, 
        const typename Calculus::DualForm1& v, 
        Image& image,
        double cut_low = 0.0, double cut_up = 1.0,
        int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        // Threshold is 0.25 instead of 0.5 because an edge connecting
        // two vertices with v=0 and v=1 should not belong to the
        // discontinuity set.
        dualForm1ToImage( calculus, v, image,
                          [] ( double x ) { return (unsigned char) ( round( x * 255.0 ) ); },
                          [] ( double x ) { return x < 0.25; },
                          cut_low, cut_up, pixel_size );
      }
 
      template <typename Calculus, typename Image>
      void primalForm1ToRGBColorImage
      ( const Calculus& calculus, 
        const typename Calculus::PrimalForm1& v, 
        Image& image, Color color,
        double cut_low = 0.0, double cut_up = 1.0,
        int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        // Threshold is 0.25 instead of 0.5 because an edge connecting
        // two vertices with v=0 and v=1 should not belong to the
        // discontinuity set.
        primalForm1ToImage( calculus, v, image,
                            [color] ( double x ) { return color; },
                            [] ( double x ) { return x < 0.25; },
                            cut_low, cut_up, pixel_size );
      }
 
      template <typename Calculus, typename Image>
      void dualForm1ToRGBColorImage
      ( const Calculus& calculus, 
        const typename Calculus::DualForm1& v, 
        Image& image, Color color,
        double cut_low = 0.0, double cut_up = 1.0,
        int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        // Threshold is 0.25 instead of 0.5 because an edge connecting
        // two vertices with v=0 and v=1 should not belong to the
        // discontinuity set.
        dualForm1ToImage( calculus, v, image,
                          [color] ( double x ) { return color; },
                          [] ( double x ) { return x < 0.25; },
                          cut_low, cut_up, pixel_size );
      }
 
      template <typename Calculus, typename AnyForm2, typename Image>
      void threeForms2ToRGBColorImage
      ( const Calculus& calculus, 
        const AnyForm2& u0, 
        const AnyForm2& u1, 
        const AnyForm2& u2, 
        Image& image,
        double cut_low = 0.0, double cut_up = 1.0,
        int pixel_size = 1 )
      {
        BOOST_CONCEPT_ASSERT(( concepts::CImage<Image> ));
        threeForms2ToImage
          ( calculus, u0, u1, u2, image,
            [] ( double r, double g, double b )
            { return Color( (unsigned char) ( round( r * 255.0 ) ),
                            (unsigned char) ( round( g * 255.0 ) ),
                            (unsigned char) ( round( b * 255.0 ) ) ); },
            cut_low, cut_up, pixel_size );
      }
 
    } // namespace dec
  } // namespace functions
 
 
  // template class DECImage2D
  template < typename TKSpace,
             typename TLinearAlgebra = EigenLinearAlgebraBackend >
  struct DECImage2D {
    typedef TKSpace                                        KSpace;
    typedef TLinearAlgebra                                 LinearAlgebra;
    BOOST_CONCEPT_ASSERT(( concepts::CCellularGridSpaceND< KSpace > )); 
    typedef typename KSpace::Space                         Space;
    typedef typename Space::Point                          Point;
    typedef typename Space::RealVector                     RealVector;
    typedef typename RealVector::Component                 Scalar;
    typedef typename KSpace::SCell                         SCell;
    typedef typename KSpace::Cell                          Cell;
    typedef typename KSpace::Surfel                        Surfel;
    typedef HyperRectDomain<Space>                         Domain;
    typedef DiscreteExteriorCalculus<2,2, LinearAlgebra>   Calculus;
    typedef DiscreteExteriorCalculusFactory<LinearAlgebra> CalculusFactory;
    typedef typename Calculus::Index                       Index;
    typedef typename Calculus::PrimalForm0                 PrimalForm0;
    typedef typename Calculus::PrimalForm1                 PrimalForm1;
    typedef typename Calculus::PrimalForm2                 PrimalForm2;
    typedef typename Calculus::PrimalIdentity0             PrimalIdentity0;
    typedef typename Calculus::PrimalIdentity1             PrimalIdentity1;
    typedef typename Calculus::PrimalIdentity2             PrimalIdentity2;
    typedef typename Calculus::PrimalDerivative0           PrimalDerivative0;
    typedef typename Calculus::PrimalDerivative1           PrimalDerivative1;
    typedef typename Calculus::DualDerivative0             DualDerivative0;
    typedef typename Calculus::DualDerivative1             DualDerivative1;
    typedef typename Calculus::PrimalAntiderivative1       PrimalAntiderivative1;
    typedef typename Calculus::PrimalAntiderivative2       PrimalAntiderivative2;
    typedef typename Calculus::DualAntiderivative1         DualAntiderivative1;
    typedef typename Calculus::DualAntiderivative2         DualAntiderivative2;
    typedef typename Calculus::PrimalHodge0                PrimalHodge0;
    typedef typename Calculus::PrimalHodge1                PrimalHodge1;
    typedef typename Calculus::PrimalHodge2                PrimalHodge2;
    typedef typename Calculus::DualHodge0                  DualHodge0;
    typedef typename Calculus::DualHodge1                  DualHodge1;
    typedef typename Calculus::DualHodge2                  DualHodge2;
    typedef typename LinearAlgebra::SolverSimplicialLLT    LinearAlgebraSolver;
    typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> 
                                                           SolverU;
    typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 1, PRIMAL, 1, PRIMAL> 
                                                           SolverV;
 
    BOOST_STATIC_ASSERT(( KSpace::dimension == 2 ));
 
    // ----------------------- Standard services ------------------------------
  public:
 
    ~DECImage2D() = default;
 
    DECImage2D( int _verbose = 1 )
      : verbose( _verbose ), 
        calculus(), 
        D0( calculus ), D1( calculus ), 
        dual_D0( calculus ), dual_D1( calculus ),
        primal_h0( calculus ), primal_h1( calculus ), primal_h2( calculus ),
        dual_h0( calculus ), dual_h1( calculus ), dual_h2( calculus )
    {}
    
    void init( Clone<KSpace> aKSpace )
    {
      calculus.myKSpace = aKSpace;
      const KSpace & K  = calculus.myKSpace;
      domain            = Domain( K.lowerBound(), K.upperBound() );
      Point  p0         = K.uKCoords( K.lowerCell() );
      Point  p1         = K.uKCoords( K.upperCell() );
      cell_domain       = Domain( p0, p1 );
 
      if ( verbose > 0 ) trace.beginBlock("building AT functionnals");
      // Adds all the cell
      for ( typename Domain::ConstIterator it = cell_domain.begin(), itE = cell_domain.end(); 
            it != itE; ++it )
        calculus.insertSCell( K.sCell( *it ) ); // ajoute toutes les cellules de Khalimsky.
      calculus.updateIndexes();
      if ( verbose > 1 ) trace.info() << calculus << std::endl;
      // Precomputes operators.
      if ( verbose > 1 ) trace.info() << "primal_D0" << std::endl;
      D0 = calculus.template derivative<0,PRIMAL>();
      if ( verbose > 1 ) trace.info() << "primal_D1" << std::endl;
      D1 = calculus.template derivative<1,PRIMAL>();
      if ( verbose > 1 ) trace.info() << "dual_D0" << std::endl;
      dual_D0   = calculus.template derivative<0,DUAL>();
      if ( verbose > 1 ) trace.info() << "dual_D1" << std::endl;
      dual_D1   = calculus.template derivative<1,DUAL>();
      if ( verbose > 1 ) trace.info() << "primal_h1" << std::endl;
      primal_h0 = calculus.template hodge<0,PRIMAL>();
      if ( verbose > 1 ) trace.info() << "primal_h1" << std::endl;
      primal_h1 = calculus.template hodge<1,PRIMAL>();
      if ( verbose > 1 ) trace.info() << "primal_h2" << std::endl;
      primal_h2 = calculus.template hodge<2,PRIMAL>();
      if ( verbose > 1 ) trace.info() << "dual_h1" << std::endl;
      dual_h0   = calculus.template hodge<0,DUAL>();
      if ( verbose > 1 ) trace.info() << "dual_h1" << std::endl;
      dual_h1   = calculus.template hodge<1,DUAL>();
      if ( verbose > 1 ) trace.info() << "dual_h2" << std::endl;
      dual_h2   = calculus.template hodge<2,DUAL>();
      if ( verbose > 0 ) trace.endBlock();
    }
 
    DECImage2D ( const DECImage2D & other ) = delete;
 
    DECImage2D ( DECImage2D && other ) = delete;
 
    DECImage2D & operator= ( const DECImage2D & other ) = delete;
 
    DECImage2D & operator= ( DECImage2D && other ) = delete;
 
    // ------------------------- Public Datas ------------------------------
  public:
 
    int verbose;
    Calculus calculus;
    Domain   domain;
    Domain   cell_domain;
    PrimalDerivative0 D0;
    PrimalDerivative1 D1;
    DualDerivative0   dual_D0;
    DualDerivative1   dual_D1;
    PrimalHodge0      primal_h0;
    PrimalHodge1      primal_h1;
    PrimalHodge2      primal_h2;
    DualHodge0        dual_h0;
    DualHodge1        dual_h1;
    DualHodge2        dual_h2;
 
  };
 
} // namespace DGtal
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined DECImageHelpers_h
 
#undef DECImageHelpers_RECURSES
#endif // else defined(DECImageHelpers_RECURSES)