DGtal  1.5.beta
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 
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  (void)randomization; //parameter not used, avoiding warning
314  do {
315  std::tie( loo, l2 ) =
318  : reg == Regularization::SNAKE
321  if ( nb % 50 == 0 )
322  trace.info() << "[Shrouds iteration " << nb
323  << " E=" << energy( reg )
324  << "] dx <= " << loo << " l2=" << l2 << std::endl;
325  r *= 0.9;
326  nb += 1;
327  } while ( loo > max_loo && nb < maxNb );
328  return std::make_pair( loo, l2 );
329  }
330 
337  {
339  return energySquaredCurvature();
340  else if ( reg == Regularization::SNAKE )
341  return energySnake();
342  else
343  return energyArea();
344  }
345 
349 
352  double energySnake();
353 
356  double energyArea();
357 
367  std::pair<double,double> oneStepAreaMinimization( const double randomization = 0.0 );
368 
405  std::pair<double,double> oneStepSnakeMinimization
406  ( const double alpha = 1.0, const double beta = 1.0, const double randomization = 0.0 );
407 
462  std::pair<double,double> oneStepSquaredCurvatureMinimization
463  ( const double randomization = 0.0 );
464 
467 
469 
470  // -------------------------- internal methods ------------------------------
471  protected:
474 
478  {
479  typedef typename Container::Tracker Tracker;
480  typedef DigitalSurface2DSlice< Tracker > Slice;
481  myT = Scalars( myPtrIdxSurface->nbVertices(), 0.5 );
482  myInvalid = myT.size();
483  myOrthDir.resize( myT.size() );
484  for ( Dimension i = 0; i < 3; ++i )
485  {
486  myNext[ i ] = std::vector<Vertex>( myT.size(), myInvalid );
487  myPrev[ i ] = std::vector<Vertex>( myT.size(), myInvalid );
488  }
489  // for each vertex, extracts its two slices.
490  for ( Vertex v = 0; v < myT.size(); ++v )
491  {
492  auto surf = myPtrIdxSurface->surfel( v );
493  Dimension k = myPtrK->sOrthDir( surf );
494  myOrthDir[ v ] = k;
495  for ( Dimension i = 0; i < 3; ++i )
496  {
497  if ( k == i ) continue; // not a valid slice
498  if ( myNext[ i ][ v ] != myInvalid ) continue; // already computed
499  Tracker* tracker = myPtrIdxSurface->container().newTracker( surf );
500  Slice slice( tracker, i );
501  if ( ! slice.isClosed() ) {
502  trace.error() << "[ShroudsRegularization::precomputeTopology]"
503  << " Shrouds works solely on closed surfaces."
504  << std::endl;
505  return;
506  }
507  auto start = slice.cstart();
508  auto next = start;
509  auto prev = next++;
510  Vertex vp = v;
511  do
512  {
513  auto sp = *prev;
514  auto sn = *next;
515  Dimension in = myPtrK->sOrthDir( sn ) == k ? i : k;
516  Dimension ip = myPtrK->sOrthDir( sp ) == k ? i : k;
517  Vertex vn = myPtrIdxSurface->getVertex( sn );
518  myNext[ ip ][ vp ] = vn;
519  myPrev[ in ][ vn ] = vp;
520  prev = next++;
521  vp = vn;
522  }
523  while ( prev != start );
524  delete tracker;
525  }
526  }
527  }
528 
530 
531  // -------------------------- data ---------------------------------
532  private:
533 
537  const KSpace* myPtrK;
558  std::vector<Dimension> myOrthDir;
560  std::vector<Vertex> myNext[ 3 ];
562  std::vector<Vertex> myPrev[ 3 ];
569 
570  }; // end of class ShroudsRegularization
571 
584  template < typename TDigitalSurfaceContainer >
588  double eps = 0.00001 )
589  {
591  }
592 
593 } // namespace surfaces
594 
595 
597 // Includes inline functions/methods if necessary.
598 #include "DGtal/geometry/surfaces/ShroudsRegularization.ih"
599 // //
601 
602 #endif // !defined ShroudsRegularization_h
603 
604 #undef ShroudsRegularization_RECURSES
605 #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()
CountedPtr< SH3::DigitalSurface > surface
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
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...