DGtal  1.2.0
ShroudsRegularization.h
1 
17 #pragma once
18 
31 #if defined(ShroudsRegularization_RECURSES)
32 #error Recursive header files inclusion detected in ShroudsRegularization.h
33 #else // defined(ShroudsRegularization_RECURSES)
35 #define ShroudsRegularization_RECURSES
36 
37 #if !defined ShroudsRegularization_h
39 #define ShroudsRegularization_h
40 
42 // Inclusions
43 #include <iostream>
44 #include "DGtal/base/Common.h"
45 #include "DGtal/base/ConstAlias.h"
46 #include "DGtal/topology/CanonicSCellEmbedder.h"
47 #include "DGtal/topology/IndexedDigitalSurface.h"
48 #include "DGtal/topology/DigitalSurface2DSlice.h"
50 
51 namespace DGtal {
72  template <typename TDigitalSurfaceContainer>
74  {
77 
78  public:
79  typedef TDigitalSurfaceContainer Container;
81  typedef typename Container::KSpace KSpace;
82  typedef typename KSpace::Space Space;
83  typedef typename Space::RealVector RealVector;
84  typedef typename Space::RealPoint RealPoint;
85  typedef typename RealVector::Component Scalar;
89  typedef IdxVertex Vertex;
90  typedef std::vector< IdxSurfel > IdxSurfelRange;
91  typedef std::vector< Scalar > Scalars;
92  typedef std::vector< RealVector > RealVectors;
93  typedef std::vector< RealPoint > RealPoints;
94 
97 
98  // ----------------------- Standard services ------------------------------
99  public:
102 
105  : myPtrIdxSurface( nullptr ), myPtrK( nullptr )
106  {}
107 
114  : myPtrIdxSurface( surface ),
115  myPtrK( &surface->container().space() ),
116  myEpsilon( 0.0001 ), myAlpha( 1.0 ), myBeta( 1.0 )
117  {
119  init();
120  }
121 
128  void init()
129  {
130  const auto embedder = CanonicSCellEmbedder<KSpace>( *myPtrK );
131  const auto nbV = myPtrIdxSurface->nbVertices();
132  myT.resize( nbV );
133  myInsV.resize( nbV );
134  myOutV.resize( nbV );
135  for ( Dimension i = 0; i < 3; ++i )
136  {
137  myPrevD[ i ].resize( nbV );
138  myNextD[ i ].resize( nbV );
139  }
140  for ( Vertex v = 0; v < myT.size(); ++v )
141  {
142  const auto s = myPtrIdxSurface->surfel( v );
143  const auto k = myPtrK->sOrthDir( s );
144  myInsV[ v ] = embedder( myPtrK->sDirectIncident( s, k ) );
145  myOutV[ v ] = embedder( myPtrK->sIndirectIncident( s, k ) );
146  }
147  }
148 
154  void setParams( double eps, double alpha = 1.0, double beta = 1.0 )
155  {
156  myEpsilon = eps;
157  myAlpha = alpha;
158  myBeta = beta;
159  }
167  std::tuple< double, double, double > getParams() const
168  {
169  return std::make_tuple( myEpsilon, myAlpha, myBeta );
170  }
171 
173 
174  // ----------------------- Accessor services ------------------------------
175  public:
178 
182  RealPoint position( const Vertex v, const double t ) const
183  {
184  return (1-t) * myInsV[ v ] + t * myOutV[ v ];
185  }
186 
189  RealPoint position( const Vertex v ) const
190  {
191  const auto t = myT[ v ];
192  return (1-t) * myInsV[ v ] + t * myOutV[ v ];
193  }
194 
197  {
198  RealPoints result( myT.size() );
199  for ( Vertex v = 0; v < myT.size(); ++v )
200  result[ v ] = position( v );
201  return result;
202  }
203 
208  Dimension orthDir( const Vertex v ) const
209  {
210  return myOrthDir[ v ];
211  }
212 
216  std::pair<Vertex,Dimension> next( const std::pair<Vertex,Dimension> & v_i ) const
217  {
218  const Vertex vn = myNext[ v_i.second ][ v_i.first ];
219  const Dimension in = myOrthDir[ vn ] == v_i.second
220  ? myOrthDir[ v_i.first ] : v_i.second;
221  return std::make_pair( vn, in );
222  }
223 
227  std::pair<Vertex,Dimension> prev( const std::pair<Vertex,Dimension> &v_i ) const
228  {
229  const Vertex vp = myPrev[ v_i.second ][ v_i.first ];
230  const Dimension ip = myOrthDir[ vp ] == v_i.second
231  ? myOrthDir[ v_i.first ] : v_i.second;
232  return std::make_pair( vp, ip );
233  }
234 
236 
237  // ----------------------- Geometric services ------------------------------
238  public:
241 
244  {
245  Scalar d = 0.0;
246  // Vertex n = 0;
247  for ( Dimension i = 0; i < 3; ++i )
248  for ( Vertex v = 0; v < myT.size(); ++v )
249  {
250  if ( myNext[ i ][ v ] == myInvalid ) continue; // not a valid slice
251  myNextD[ i ][ v ] = ( position( myNext[ i ][ v ] ) - position( v ) ).norm();
252  myPrevD[ i ][ v ] = ( position( myPrev[ i ][ v ] ) - position( v ) ).norm();
253  d += myNextD[ i ][ v ];
254  // n += 1;
255  }
256  }
257 
261  Scalar c1( const std::pair<Vertex,Dimension> &v_i ) const
262  {
263  const Scalar din = myNextD[ v_i.second ][ v_i.first ];
264  const Scalar dip = myPrevD[ v_i.second ][ v_i.first ];
265  return 1.0 / (din + dip);
266  }
267 
271  std::tuple<Scalar,Scalar,Scalar> c2_all( const std::pair<Vertex,Dimension> &v_i ) const
272  {
273  const Scalar din = myNextD[ v_i.second ][ v_i.first ];
274  const Scalar dip = myPrevD[ v_i.second ][ v_i.first ];
275  return std::make_tuple( 2.0 / ( din * ( din + dip ) ),
276  2.0 / ( din * dip ),
277  2.0 / ( dip * ( din + dip ) ) );
278  }
279 
281 
282  // -------------------------- regularization services ----------------------------
283  public:
286 
303  std::pair<double,double>
305  const double randomization = 0.0,
306  const double max_loo = 0.0001,
307  const int maxNb = 100 )
308  {
309  double loo = 0.0;
310  double l2 = 0.0;
311  int nb = 0;
312  double r = 0.5;
313  do {
314  std::tie( loo, l2 ) =
317  : reg == Regularization::SNAKE
320  if ( nb % 50 == 0 )
321  trace.info() << "[Shrouds iteration " << nb
322  << " E=" << energy( reg )
323  << "] dx <= " << loo << " l2=" << l2 << std::endl;
324  r *= 0.9;
325  nb += 1;
326  } while ( loo > max_loo && nb < maxNb );
327  return std::make_pair( loo, l2 );
328  }
329 
336  {
338  return energySquaredCurvature();
339  else if ( reg == Regularization::SNAKE )
340  return energySnake();
341  else
342  return energyArea();
343  }
344 
348 
351  double energySnake();
352 
355  double energyArea();
356 
366  std::pair<double,double> oneStepAreaMinimization( const double randomization = 0.0 );
367 
404  std::pair<double,double> oneStepSnakeMinimization
405  ( const double alpha = 1.0, const double beta = 1.0, const double randomization = 0.0 );
406 
461  std::pair<double,double> oneStepSquaredCurvatureMinimization
462  ( const double randomization = 0.0 );
463 
466 
468 
469  // -------------------------- internal methods ------------------------------
470  protected:
473 
477  {
478  typedef typename Container::Tracker Tracker;
479  typedef DigitalSurface2DSlice< Tracker > Slice;
480  myT = Scalars( myPtrIdxSurface->nbVertices(), 0.5 );
481  myInvalid = myT.size();
482  myOrthDir.resize( myT.size() );
483  for ( Dimension i = 0; i < 3; ++i )
484  {
485  myNext[ i ] = std::vector<Vertex>( myT.size(), myInvalid );
486  myPrev[ i ] = std::vector<Vertex>( myT.size(), myInvalid );
487  }
488  // for each vertex, extracts its two slices.
489  for ( Vertex v = 0; v < myT.size(); ++v )
490  {
491  auto surf = myPtrIdxSurface->surfel( v );
492  Dimension k = myPtrK->sOrthDir( surf );
493  myOrthDir[ v ] = k;
494  for ( Dimension i = 0; i < 3; ++i )
495  {
496  if ( k == i ) continue; // not a valid slice
497  if ( myNext[ i ][ v ] != myInvalid ) continue; // already computed
498  Tracker* tracker = myPtrIdxSurface->container().newTracker( surf );
499  Slice slice( tracker, i );
500  if ( ! slice.isClosed() ) {
501  trace.error() << "[ShroudsRegularization::precomputeTopology]"
502  << " Shrouds works solely on closed surfaces."
503  << std::endl;
504  return;
505  }
506  auto start = slice.cstart();
507  auto next = start;
508  auto prev = next++;
509  Vertex vp = v;
510  do
511  {
512  auto sp = *prev;
513  auto sn = *next;
514  Dimension in = myPtrK->sOrthDir( sn ) == k ? i : k;
515  Dimension ip = myPtrK->sOrthDir( sp ) == k ? i : k;
516  Vertex vn = myPtrIdxSurface->getVertex( sn );
517  myNext[ ip ][ vp ] = vn;
518  myPrev[ in ][ vn ] = vp;
519  prev = next++;
520  vp = vn;
521  }
522  while ( prev != start );
523  delete tracker;
524  }
525  }
526  }
527 
529 
530  // -------------------------- data ---------------------------------
531  private:
532 
536  const KSpace* myPtrK;
557  std::vector<Dimension> myOrthDir;
559  std::vector<Vertex> myNext[ 3 ];
561  std::vector<Vertex> myPrev[ 3 ];
568 
569  }; // end of class ShroudsRegularization
570 
583  template < typename TDigitalSurfaceContainer >
587  double eps = 0.00001 )
588  {
590  }
591 
592 } // namespace surfaces
593 
594 
596 // Includes inline functions/methods if necessary.
597 #include "DGtal/geometry/surfaces/ShroudsRegularization.ih"
598 // //
600 
601 #endif // !defined ShroudsRegularization_h
602 
603 #undef ShroudsRegularization_RECURSES
604 #endif // else defined(ShroudsRegularization_RECURSES)
Aim: Smart pointer based on reference counts.
Definition: CountedPtr.h:80
Aim: Represents a 2-dimensional slice in a DigitalSurface. In a sense, it is a 4-connected contour,...
Aim: Represents a digital surface with the topology of its dual surface. Its aim is to mimick the sta...
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:614
Aim: Implements the Shrouds Regularization algorithm of Nielson et al .
ShroudsRegularization(CountedPtr< IdxDigitalSurface > surface)
std::vector< RealVector > RealVectors
Regularization
The enum class specifying the possible shrouds regularization.
std::tuple< double, double, double > getParams() const
RealPoint position(const Vertex v) const
RealPoint position(const Vertex v, const double t) const
Vertex myInvalid
the index of the invalid vertex.
std::vector< IdxSurfel > IdxSurfelRange
Scalar myAlpha
The alpha parameter for Snake first order regularization (~ area)
std::tuple< Scalar, Scalar, Scalar > c2_all(const std::pair< Vertex, Dimension > &v_i) const
std::vector< Vertex > myNext[3]
for each vertex, its successor on the slice of given axis direction.
std::pair< double, double > oneStepSquaredCurvatureMinimization(const double randomization=0.0)
std::vector< Vertex > myPrev[3]
for each vertex, its predessor on the slice of given axis direction.
Dimension orthDir(const Vertex v) const
double energy(const Regularization reg=Regularization::SQUARED_CURVATURE)
std::pair< double, double > oneStepAreaMinimization(const double randomization=0.0)
CountedPtr< IdxDigitalSurface > myPtrIdxSurface
the indexed digital surface (internal surface representation).
std::pair< Vertex, Dimension > next(const std::pair< Vertex, Dimension > &v_i) const
std::vector< RealPoint > RealPoints
std::pair< double, double > oneStepSnakeMinimization(const double alpha=1.0, const double beta=1.0, const double randomization=0.0)
const KSpace * myPtrK
A const pointer to the cellular space in which lives the digital surface.
void parameterize()
Computes the distances between the vertices along slices.
ShroudsRegularization< Container > Self
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)
Scalar myBeta
The beta parameter for Snake second order regularization (~ curvature)
std::pair< Vertex, Dimension > prev(const std::pair< Vertex, Dimension > &v_i) const
void setParams(double eps, double alpha=1.0, double beta=1.0)
void enforceBounds()
Forces t to stay in ]0,1[.
Scalar c1(const std::pair< Vertex, Dimension > &v_i) const
BOOST_CONCEPT_ASSERT((concepts::CDigitalSurfaceContainer< TDigitalSurfaceContainer >))
std::vector< Dimension > myOrthDir
the direction axis of each dual edge.
IdxDigitalSurface::Surfel IdxSurfel
IndexedDigitalSurface< Container > IdxDigitalSurface
TDigitalSurfaceContainer Container
IdxDigitalSurface::Vertex IdxVertex
ShroudsRegularization()
Default constructor. The object is not valid.
std::ostream & error()
std::ostream & info()
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:137
Trace trace
Definition: Common.h:154
ShroudsRegularization< TDigitalSurfaceContainer > makeShroudsRegularization(CountedPtr< IndexedDigitalSurface< TDigitalSurfaceContainer > > surface, double eps=0.00001)
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Aim: The digital surface container concept describes a minimal set of inner types and methods so as t...