ShroudsRegularization.h

 
#pragma once
 
#if defined(ShroudsRegularization_RECURSES)
#error Recursive header files inclusion detected in ShroudsRegularization.h
#else // defined(ShroudsRegularization_RECURSES)
#define ShroudsRegularization_RECURSES
 
#if !defined ShroudsRegularization_h
#define ShroudsRegularization_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/base/ConstAlias.h"
#include "DGtal/topology/CanonicSCellEmbedder.h"
#include "DGtal/topology/IndexedDigitalSurface.h"
#include "DGtal/topology/DigitalSurface2DSlice.h"
 
namespace DGtal {
  template <typename TDigitalSurfaceContainer>
  class ShroudsRegularization
  {
    BOOST_CONCEPT_ASSERT
    (( concepts::CDigitalSurfaceContainer< TDigitalSurfaceContainer > ) );
    
  public:
    typedef TDigitalSurfaceContainer           Container;
    typedef ShroudsRegularization< Container > Self;
    typedef typename Container::KSpace         KSpace;
    typedef typename KSpace::Space             Space;
    typedef typename Space::RealVector         RealVector;
    typedef typename Space::RealPoint          RealPoint;
    typedef typename RealVector::Component     Scalar;
    typedef IndexedDigitalSurface< Container > IdxDigitalSurface;
    typedef typename IdxDigitalSurface::Vertex IdxVertex;
    typedef typename IdxDigitalSurface::Surfel IdxSurfel;
    typedef IdxVertex                          Vertex;
    typedef std::vector< IdxSurfel >           IdxSurfelRange;
    typedef std::vector< Scalar >              Scalars;
    typedef std::vector< RealVector >          RealVectors;
    typedef std::vector< RealPoint >           RealPoints;
 
    enum class Regularization { AREA, SNAKE, SQUARED_CURVATURE };
    
    // ----------------------- Standard services ------------------------------
  public:
    
    ShroudsRegularization()
      : myPtrIdxSurface( nullptr ), myPtrK( nullptr )
    {}
    
    ShroudsRegularization( CountedPtr< IdxDigitalSurface > surface )
      : myPtrIdxSurface( surface ),
        myPtrK( &surface->container().space() ),
        myEpsilon( 0.0001 ), myAlpha( 1.0 ), myBeta( 1.0 )
    {
      precomputeTopology();
      init();
    }
 
    void init()
    {
      const auto embedder = CanonicSCellEmbedder<KSpace>( *myPtrK );
      const auto nbV      = myPtrIdxSurface->nbVertices();
      myT.resize( nbV );
      myInsV.resize( nbV );
      myOutV.resize( nbV );
      for ( Dimension i = 0; i < 3; ++i )
        {
          myPrevD[ i ].resize( nbV );
          myNextD[ i ].resize( nbV );
        }
      for ( Vertex v = 0; v < myT.size(); ++v )
        {
          const auto s = myPtrIdxSurface->surfel( v );
          const auto k = myPtrK->sOrthDir( s );
          myInsV[ v ]  = embedder( myPtrK->sDirectIncident( s, k ) );
          myOutV[ v ]  = embedder( myPtrK->sIndirectIncident( s, k ) );
        }
    }
 
    void setParams( double eps, double alpha = 1.0, double beta = 1.0 )
    {
      myEpsilon = eps;
      myAlpha   = alpha;
      myBeta    = beta;
    }
    std::tuple< double, double, double > getParams() const
    {
      return std::make_tuple( myEpsilon, myAlpha, myBeta );
    }
    
    
    // ----------------------- Accessor services ------------------------------
  public:
    
    RealPoint position( const Vertex v, const double t ) const
    {
      return (1-t) * myInsV[ v ] + t * myOutV[ v ];
    }
    
    RealPoint position( const Vertex v ) const
    {
      const auto t = myT[ v ];
      return (1-t) * myInsV[ v ] + t * myOutV[ v ];
    }
    
    RealPoints positions() const
    {
      RealPoints result( myT.size() );
      for ( Vertex v = 0; v < myT.size(); ++v )
        result[ v ] = position( v );
      return result;
    }
 
    Dimension orthDir( const Vertex v ) const
    {
      return myOrthDir[ v ];
    }
    
    std::pair<Vertex,Dimension> next( const std::pair<Vertex,Dimension> & v_i ) const
    {
      const Vertex    vn = myNext[ v_i.second ][ v_i.first ];
      const Dimension in = myOrthDir[ vn ] == v_i.second
        ? myOrthDir[ v_i.first ] : v_i.second;
      return std::make_pair( vn, in );
    }
 
    std::pair<Vertex,Dimension> prev( const std::pair<Vertex,Dimension> &v_i ) const
    {
      const Vertex    vp = myPrev[ v_i.second ][ v_i.first ];
      const Dimension ip = myOrthDir[ vp ] == v_i.second
        ? myOrthDir[ v_i.first ] : v_i.second;
      return std::make_pair( vp, ip );
    }
    
 
    // ----------------------- Geometric services ------------------------------
  public:
 
    void parameterize()
    {
      Scalar d = 0.0;
      // Vertex n = 0;
      for ( Dimension i = 0; i < 3; ++i )
        for ( Vertex v = 0; v < myT.size(); ++v )
          {
            if ( myNext[ i ][ v ] == myInvalid )  continue; // not a valid slice
            myNextD[ i ][ v ] = ( position( myNext[ i ][ v ] ) - position( v ) ).norm();
            myPrevD[ i ][ v ] = ( position( myPrev[ i ][ v ] ) - position( v ) ).norm();
            d += myNextD[ i ][ v ];
            //      n += 1;
          }
    }
 
    Scalar c1( const std::pair<Vertex,Dimension> &v_i ) const
    {
      const Scalar din = myNextD[ v_i.second ][ v_i.first ];
      const Scalar dip = myPrevD[ v_i.second ][ v_i.first ];
      return 1.0 / (din + dip);
    }
    
    std::tuple<Scalar,Scalar,Scalar> c2_all( const std::pair<Vertex,Dimension> &v_i ) const
    {
      const Scalar din = myNextD[ v_i.second ][ v_i.first ];
      const Scalar dip = myPrevD[ v_i.second ][ v_i.first ];
      return std::make_tuple( 2.0 / ( din * ( din + dip ) ),
                              2.0 / ( din * dip ),
                              2.0 / ( dip * ( din + dip ) ) );
    }
    
 
    // -------------------------- regularization services ----------------------------
  public:
 
    std::pair<double,double>
    regularize( const Regularization   reg = Regularization::SQUARED_CURVATURE,
                const double randomization = 0.0,
                const double       max_loo = 0.0001,
                const int            maxNb = 100 )
    {
      double loo      = 0.0;
      double  l2      = 0.0;
      int     nb      = 0;
      double   r      = 0.5;
      do {
        std::tie( loo, l2 ) =
          reg == Regularization::SQUARED_CURVATURE
          ? oneStepSquaredCurvatureMinimization( r )
          : reg == Regularization::SNAKE
          ? oneStepSnakeMinimization( myAlpha, myBeta, r )
          : oneStepAreaMinimization( r );
        if ( nb % 50 == 0 )
          trace.info() << "[Shrouds iteration " << nb
                       << " E=" << energy( reg )
                       << "] dx <= " << loo << " l2=" << l2 << std::endl;
        r  *= 0.9;
        nb += 1;
      } while ( loo > max_loo && nb < maxNb );
      return std::make_pair( loo, l2 );
    }
 
    double energy( const Regularization reg = Regularization::SQUARED_CURVATURE )
    {
      if ( reg == Regularization::SQUARED_CURVATURE )
        return energySquaredCurvature();
      else if ( reg == Regularization::SNAKE )
        return energySnake();
      else
        return energyArea();
    }
 
    double energySquaredCurvature();
 
    double energySnake();
 
    double energyArea();
    
    std::pair<double,double> oneStepAreaMinimization( const double randomization = 0.0 );
 
    std::pair<double,double> oneStepSnakeMinimization
    ( const double alpha = 1.0, const double beta = 1.0, const double randomization = 0.0 );
 
    std::pair<double,double> oneStepSquaredCurvatureMinimization
    ( const double randomization = 0.0 );
 
    void enforceBounds();
 
    
    // -------------------------- internal methods ------------------------------
  protected:
    
    void precomputeTopology()
    {
      typedef typename Container::Tracker      Tracker;
      typedef DigitalSurface2DSlice< Tracker > Slice;
      myT       = Scalars( myPtrIdxSurface->nbVertices(), 0.5 );
      myInvalid = myT.size();
      myOrthDir.resize( myT.size() );
      for ( Dimension i = 0; i < 3; ++i )
      {
        myNext[ i ] = std::vector<Vertex>( myT.size(), myInvalid );
        myPrev[ i ] = std::vector<Vertex>( myT.size(), myInvalid );
      }
      // for each vertex, extracts its two slices.
      for ( Vertex v = 0; v < myT.size(); ++v )
        {
          auto surf      = myPtrIdxSurface->surfel( v );
          Dimension k    = myPtrK->sOrthDir( surf );
          myOrthDir[ v ] = k;
          for ( Dimension i = 0; i < 3; ++i )
            {
              if ( k == i )   continue; // not a valid slice
              if ( myNext[ i ][ v ] != myInvalid ) continue; // already computed
              Tracker* tracker = myPtrIdxSurface->container().newTracker( surf );
              Slice    slice( tracker, i );
              if ( ! slice.isClosed() ) {
                trace.error() << "[ShroudsRegularization::precomputeTopology]"
                              << " Shrouds works solely on closed surfaces."
                              << std::endl;
                return;
              }
              auto start = slice.cstart();
              auto  next = start;
              auto  prev = next++;
              Vertex  vp = v;
              do
                {
                  auto sp = *prev;
                  auto sn = *next;
                  Dimension in = myPtrK->sOrthDir( sn ) == k ? i : k;
                  Dimension ip = myPtrK->sOrthDir( sp ) == k ? i : k;
                  Vertex vn = myPtrIdxSurface->getVertex( sn );
                  myNext[ ip ][ vp ] = vn;
                  myPrev[ in ][ vn ] = vp;
                  prev = next++;
                  vp   = vn;
                }
              while ( prev != start );      
              delete tracker;
            }
        }
    }
 
    
    // -------------------------- data ---------------------------------
  private:
 
    CountedPtr<IdxDigitalSurface>      myPtrIdxSurface;
    const KSpace*                      myPtrK;
    Scalar                             myEpsilon;
    Scalar                             myAlpha;
    Scalar                             myBeta;
    Vertex                             myInvalid;
    Scalars                            myT;
    RealPoints                         myInsV;
    RealPoints                         myOutV;
    std::vector<Dimension>             myOrthDir;
    std::vector<Vertex>                myNext[ 3 ];
    std::vector<Vertex>                myPrev[ 3 ];
    Scalars                            myNextD[ 3 ];
    Scalars                            myPrevD[ 3 ];
    
  }; // end of class ShroudsRegularization
  
  template < typename TDigitalSurfaceContainer > 
  ShroudsRegularization<TDigitalSurfaceContainer>
  makeShroudsRegularization
  ( CountedPtr< IndexedDigitalSurface< TDigitalSurfaceContainer > > surface,
    double eps = 0.00001 )
  {
    return ShroudsRegularization<TDigitalSurfaceContainer>( surface, eps );
  }
  
} // namespace surfaces
 
 
// Includes inline functions/methods if necessary.
#include "DGtal/geometry/surfaces/ShroudsRegularization.ih"
//                                                                           //
 
#endif // !defined ShroudsRegularization_h
 
#undef ShroudsRegularization_RECURSES
#endif // else defined(ShroudsRegularization_RECURSES)