DGtal  1.1.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)
34 
35 #define ShroudsRegularization_RECURSES
36 
37 #if !defined ShroudsRegularization_h
38 
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)
DGtal::ShroudsRegularization::oneStepSquaredCurvatureMinimization
std::pair< double, double > oneStepSquaredCurvatureMinimization(const double randomization=0.0)
DGtal::ShroudsRegularization::position
RealPoint position(const Vertex v) const
Definition: ShroudsRegularization.h:189
DGtal::ShroudsRegularization::KSpace
Container::KSpace KSpace
Definition: ShroudsRegularization.h:81
KSpace
Z3i::KSpace KSpace
Definition: sphereCotangentLaplaceOperator.cpp:70
DGtal::ShroudsRegularization::IdxSurfel
IdxDigitalSurface::Surfel IdxSurfel
Definition: ShroudsRegularization.h:88
DGtal::ShroudsRegularization::myEpsilon
Scalar myEpsilon
Definition: ShroudsRegularization.h:539
DGtal::IndexedDigitalSurface::Vertex
VertexIndex Vertex
Definition: IndexedDigitalSurface.h:121
DGtal::ShroudsRegularization::myInvalid
Vertex myInvalid
the index of the invalid vertex.
Definition: ShroudsRegularization.h:545
DGtal::Trace::error
std::ostream & error()
DGtal::ShroudsRegularization::BOOST_CONCEPT_ASSERT
BOOST_CONCEPT_ASSERT((concepts::CDigitalSurfaceContainer< TDigitalSurfaceContainer >))
DGtal::ShroudsRegularization::position
RealPoint position(const Vertex v, const double t) const
Definition: ShroudsRegularization.h:182
DGtal::IndexedDigitalSurface
Aim: Represents a digital surface with the topology of its dual surface. Its aim is to mimick the sta...
Definition: IndexedDigitalSurface.h:97
DGtal::ShroudsRegularization::energySnake
double energySnake()
DGtal::ShroudsRegularization::myInsV
RealPoints myInsV
Definition: ShroudsRegularization.h:552
DGtal::ShroudsRegularization::energyArea
double energyArea()
DGtal::trace
Trace trace
Definition: Common.h:150
DGtal::makeShroudsRegularization
ShroudsRegularization< TDigitalSurfaceContainer > makeShroudsRegularization(CountedPtr< IndexedDigitalSurface< TDigitalSurfaceContainer > > surface, double eps=0.00001)
Definition: ShroudsRegularization.h:586
DGtal::ShroudsRegularization::myPtrK
const KSpace * myPtrK
A const pointer to the cellular space in which lives the digital surface.
Definition: ShroudsRegularization.h:536
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:133
DGtal::ShroudsRegularization::myBeta
Scalar myBeta
The beta parameter for Snake second order regularization (~ curvature)
Definition: ShroudsRegularization.h:543
DGtal::ShroudsRegularization::IdxDigitalSurface
IndexedDigitalSurface< Container > IdxDigitalSurface
Definition: ShroudsRegularization.h:86
DGtal::ShroudsRegularization::prev
std::pair< Vertex, Dimension > prev(const std::pair< Vertex, Dimension > &v_i) const
Definition: ShroudsRegularization.h:227
DGtal::ShroudsRegularization::oneStepAreaMinimization
std::pair< double, double > oneStepAreaMinimization(const double randomization=0.0)
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:209
DGtal::ShroudsRegularization::regularize
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)
Definition: ShroudsRegularization.h:304
DGtal::ShroudsRegularization::oneStepSnakeMinimization
std::pair< double, double > oneStepSnakeMinimization(const double alpha=1.0, const double beta=1.0, const double randomization=0.0)
DGtal::ShroudsRegularization::energy
double energy(const Regularization reg=Regularization::SQUARED_CURVATURE)
Definition: ShroudsRegularization.h:335
DGtal::concepts::CDigitalSurfaceContainer
Aim: The digital surface container concept describes a minimal set of inner types and methods so as t...
Definition: CDigitalSurfaceContainer.h:128
DGtal::ShroudsRegularization::parameterize
void parameterize()
Computes the distances between the vertices along slices.
Definition: ShroudsRegularization.h:243
DGtal::PointVector::Component
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:614
DGtal::SpaceND
Definition: SpaceND.h:96
DGtal::ShroudsRegularization::IdxSurfelRange
std::vector< IdxSurfel > IdxSurfelRange
Definition: ShroudsRegularization.h:90
DGtal::Trace::info
std::ostream & info()
DGtal::ShroudsRegularization::RealPoints
std::vector< RealPoint > RealPoints
Definition: ShroudsRegularization.h:93
DGtal::ShroudsRegularization::RealPoint
Space::RealPoint RealPoint
Definition: ShroudsRegularization.h:84
DGtal::ShroudsRegularization::getParams
std::tuple< double, double, double > getParams() const
Definition: ShroudsRegularization.h:167
DGtal::ShroudsRegularization::energySquaredCurvature
double energySquaredCurvature()
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition: ClosedIntegerHalfPlane.h:49
DGtal::ShroudsRegularization::init
void init()
Definition: ShroudsRegularization.h:128
DGtal::ShroudsRegularization::orthDir
Dimension orthDir(const Vertex v) const
Definition: ShroudsRegularization.h:208
DGtal::ShroudsRegularization::myPtrIdxSurface
CountedPtr< IdxDigitalSurface > myPtrIdxSurface
the indexed digital surface (internal surface representation).
Definition: ShroudsRegularization.h:534
DGtal::ShroudsRegularization::IdxVertex
IdxDigitalSurface::Vertex IdxVertex
Definition: ShroudsRegularization.h:87
DGtal::ShroudsRegularization::c1
Scalar c1(const std::pair< Vertex, Dimension > &v_i) const
Definition: ShroudsRegularization.h:261
DGtal::ShroudsRegularization
Aim: Implements the Shrouds Regularization algorithm of Nielson et al .
Definition: ShroudsRegularization.h:74
DGtal::ShroudsRegularization::Space
KSpace::Space Space
Definition: ShroudsRegularization.h:82
DGtal::ShroudsRegularization::c2_all
std::tuple< Scalar, Scalar, Scalar > c2_all(const std::pair< Vertex, Dimension > &v_i) const
Definition: ShroudsRegularization.h:271
DGtal::ShroudsRegularization::Regularization::SNAKE
@ SNAKE
DGtal::ShroudsRegularization::Scalar
RealVector::Component Scalar
Definition: ShroudsRegularization.h:85
DGtal::CountedPtr
Aim: Smart pointer based on reference counts.
Definition: CountedPtr.h:80
DGtal::ShroudsRegularization::RealVectors
std::vector< RealVector > RealVectors
Definition: ShroudsRegularization.h:92
DGtal::ShroudsRegularization::myOutV
RealPoints myOutV
Definition: ShroudsRegularization.h:555
DGtal::ShroudsRegularization::myPrevD
Scalars myPrevD[3]
Definition: ShroudsRegularization.h:567
DGtal::ShroudsRegularization::RealVector
Space::RealVector RealVector
Definition: ShroudsRegularization.h:83
DGtal::ShroudsRegularization::positions
RealPoints positions() const
Definition: ShroudsRegularization.h:196
DGtal::ShroudsRegularization::Scalars
std::vector< Scalar > Scalars
Definition: ShroudsRegularization.h:91
DGtal::ShroudsRegularization::setParams
void setParams(double eps, double alpha=1.0, double beta=1.0)
Definition: ShroudsRegularization.h:154
DGtal::ShroudsRegularization::myNextD
Scalars myNextD[3]
Definition: ShroudsRegularization.h:564
DGtal::ShroudsRegularization::myT
Scalars myT
Definition: ShroudsRegularization.h:549
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
DGtal::ShroudsRegularization::Self
ShroudsRegularization< Container > Self
Definition: ShroudsRegularization.h:80
DGtal::DigitalSurface2DSlice
Aim: Represents a 2-dimensional slice in a DigitalSurface. In a sense, it is a 4-connected contour,...
Definition: DigitalSurface2DSlice.h:93
DGtal::ShroudsRegularization::Regularization
Regularization
The enum class specifying the possible shrouds regularization.
Definition: ShroudsRegularization.h:96
DGtal::ShroudsRegularization::ShroudsRegularization
ShroudsRegularization()
Default constructor. The object is not valid.
Definition: ShroudsRegularization.h:104
DGtal::ShroudsRegularization::ShroudsRegularization
ShroudsRegularization(CountedPtr< IdxDigitalSurface > surface)
Definition: ShroudsRegularization.h:113
DGtal::ShroudsRegularization::precomputeTopology
void precomputeTopology()
Definition: ShroudsRegularization.h:476
DGtal::ShroudsRegularization::myAlpha
Scalar myAlpha
The alpha parameter for Snake first order regularization (~ area)
Definition: ShroudsRegularization.h:541
DGtal::ShroudsRegularization::myNext
std::vector< Vertex > myNext[3]
for each vertex, its successor on the slice of given axis direction.
Definition: ShroudsRegularization.h:559
DGtal::ShroudsRegularization::Container
TDigitalSurfaceContainer Container
Definition: ShroudsRegularization.h:79
DGtal::ShroudsRegularization::enforceBounds
void enforceBounds()
Forces t to stay in ]0,1[.
DGtal::CanonicSCellEmbedder< KSpace >
DGtal::ShroudsRegularization::Vertex
IdxVertex Vertex
Definition: ShroudsRegularization.h:89
DGtal::ShroudsRegularization::Regularization::AREA
@ AREA
DGtal::ShroudsRegularization::next
std::pair< Vertex, Dimension > next(const std::pair< Vertex, Dimension > &v_i) const
Definition: ShroudsRegularization.h:216
DGtal::ShroudsRegularization::myPrev
std::vector< Vertex > myPrev[3]
for each vertex, its predessor on the slice of given axis direction.
Definition: ShroudsRegularization.h:561
DGtal::ShroudsRegularization::Regularization::SQUARED_CURVATURE
@ SQUARED_CURVATURE
DGtal::ShroudsRegularization::myOrthDir
std::vector< Dimension > myOrthDir
the direction axis of each dual edge.
Definition: ShroudsRegularization.h:557