DGtal 1.4.2
Loading...
Searching...
No Matches
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
51namespace 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;
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;
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
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 >
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.
TEuclideanRing Component
Type for Vector elements.
Aim: Implements the Shrouds Regularization algorithm of Nielson et al nielson2003shrouds.
ShroudsRegularization(CountedPtr< IdxDigitalSurface > surface)
Regularization
The enum class specifying the possible shrouds regularization.
std::tuple< double, double, double > getParams() const
std::pair< Vertex, Dimension > prev(const std::pair< Vertex, Dimension > &v_i) const
std::pair< double, double > oneStepAreaMinimization(const double randomization=0.0)
IdxDigitalSurface::Vertex IdxVertex
RealPoint position(const Vertex v) const
RealPoint position(const Vertex v, const double t) const
Vertex myInvalid
the index of the invalid vertex.
Scalar myAlpha
The alpha parameter for Snake first order regularization (~ area)
ShroudsRegularization< Container > Self
std::vector< Vertex > myNext[3]
for each vertex, its successor on the slice of given axis direction.
std::tuple< Scalar, Scalar, Scalar > c2_all(const std::pair< Vertex, Dimension > &v_i) const
std::vector< Vertex > myPrev[3]
for each vertex, its predessor on the slice of given axis direction.
Dimension orthDir(const Vertex v) const
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 energy(const Regularization reg=Regularization::SQUARED_CURVATURE)
CountedPtr< IdxDigitalSurface > myPtrIdxSurface
the indexed digital surface (internal surface representation).
std::vector< RealVector > RealVectors
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.
std::vector< IdxSurfel > IdxSurfelRange
std::vector< RealPoint > RealPoints
TDigitalSurfaceContainer Container
Scalar myBeta
The beta parameter for Snake second order regularization (~ curvature)
void setParams(double eps, double alpha=1.0, double beta=1.0)
void enforceBounds()
Forces t to stay in ]0,1[.
IdxDigitalSurface::Surfel IdxSurfel
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.
IndexedDigitalSurface< Container > IdxDigitalSurface
std::pair< double, double > oneStepSquaredCurvatureMinimization(const double randomization=0.0)
std::pair< Vertex, Dimension > next(const std::pair< Vertex, Dimension > &v_i) const
std::pair< double, double > oneStepSnakeMinimization(const double alpha=1.0, const double beta=1.0, const double randomization=0.0)
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.
ShroudsRegularization< TDigitalSurfaceContainer > makeShroudsRegularization(CountedPtr< IndexedDigitalSurface< TDigitalSurfaceContainer > > surface, double eps=0.00001)
DGtal::uint32_t Dimension
Definition Common.h:136
Trace trace
Definition Common.h:153
Aim: A trivial embedder for signed cell, which corresponds to the canonic injection of cell centroids...
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...