DGtal  1.3.beta
ATSolver2D.h
1 
17 #pragma once
18 
31 #if defined(ATSolver2D_RECURSES)
32 #error Recursive header files inclusion detected in ATSolver2D.h
33 #else // defined(ATSolver2D_RECURSES)
34 
35 #define ATSolver2D_RECURSES
36 
37 #if !defined ATSolver2D_h
38 
39 #define ATSolver2D_h
40 
42 // Inclusions
43 #include <iostream>
44 #include <sstream>
45 #include <tuple>
46 #include "DGtal/base/Common.h"
47 #include "DGtal/base/ConstAlias.h"
48 #include "DGtal/math/linalg/EigenSupport.h"
49 #include "DGtal/dec/DiscreteExteriorCalculus.h"
50 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
51 #include "DGtal/dec/DECHelpers.h"
53 
54 namespace DGtal
55 {
56 
58  // template class ATSolver2D
88  template < typename TKSpace,
89  typename TLinearAlgebra = EigenLinearAlgebraBackend >
90  class ATSolver2D
91  {
92  // ----------------------- Standard services ------------------------------
93  public:
94 
95  typedef TKSpace KSpace;
96  typedef TLinearAlgebra LinearAlgebra;
97  typedef ATSolver2D< KSpace, LinearAlgebra > Self;
98 
99  static const Dimension dimension = KSpace::dimension;
100 
103  enum CellOutputPolicy { Average,
104  Minimum,
105  Maximum,
106  };
107  typedef typename KSpace::Space Space;
108  typedef typename Space::RealVector RealVector;
109  typedef typename RealVector::Component Scalar;
110  typedef typename KSpace::SCell SCell;
111  typedef typename KSpace::Cell Cell;
112  typedef typename KSpace::Surfel Surfel;
113  typedef HyperRectDomain<Space> Domain;
114  typedef DiscreteExteriorCalculus<2,dimension, LinearAlgebra> Calculus;
115  typedef typename KSpace::template SurfelMap<double>::Type SmallestEpsilonMap;
116  typedef typename Calculus::Index Index;
117  typedef typename Calculus::PrimalForm0 PrimalForm0;
118  typedef typename Calculus::PrimalForm1 PrimalForm1;
119  typedef typename Calculus::PrimalForm2 PrimalForm2;
120  typedef typename Calculus::PrimalIdentity0 PrimalIdentity0;
121  typedef typename Calculus::PrimalIdentity1 PrimalIdentity1;
122  typedef typename Calculus::PrimalIdentity2 PrimalIdentity2;
123  typedef typename Calculus::PrimalDerivative0 PrimalDerivative0;
124  typedef typename Calculus::PrimalDerivative1 PrimalDerivative1;
125  typedef typename Calculus::PrimalAntiderivative1 PrimalAntiderivative1;
126  typedef typename Calculus::PrimalAntiderivative2 PrimalAntiderivative2;
127  typedef typename Calculus::PrimalHodge0 PrimalHodge0;
128  typedef typename Calculus::PrimalHodge1 PrimalHodge1;
129  typedef typename Calculus::PrimalHodge2 PrimalHodge2;
130  typedef typename KSpace::template SurfelMap<Index>::Type Surfel2IndexMap;
131 
132  // SparseLU is so much faster than SparseQR
133  // SimplicialLLT is much faster than SparseLU
134  // SimplicialLDLT is as fast as SimplicialLLT but more robust
135  // typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
136  // typedef EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;
137  // typedef EigenLinearAlgebraBackend::SolverSimplicialLLT LinearAlgebraSolver;
138  typedef EigenLinearAlgebraBackend::SolverSimplicialLDLT LinearAlgebraSolver;
139  typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL> SolverU2;
140  typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> SolverV0;
141 
142  protected:
144  CountedConstPtrOrConstPtr< Calculus > ptrCalculus;
146  PrimalDerivative0 primal_D0;
148  PrimalDerivative1 primal_D1;
150  PrimalDerivative0 M01;
152  PrimalDerivative1 M12;
154  PrimalAntiderivative2 primal_AD2;
156  PrimalIdentity2 alpha_Id2;
158  PrimalIdentity0 l_1_over_4e_Id0;
160  std::vector<PrimalForm2> g2;
162  std::vector<PrimalForm2> alpha_g2;
164  std::vector<PrimalForm2> u2;
167  PrimalForm0 v0;
169  PrimalForm0 former_v0;
171  PrimalForm0 l_1_over_4e;
172 
173  public:
174  // The map Surfel -> Index that gives the index of the surfel in 2-forms.
175  Surfel2IndexMap surfel2idx;
178  SmallestEpsilonMap smallest_epsilon_map;
181  double alpha;
184  double lambda;
187  double epsilon;
189  bool normalize_u2;
191  int verbose;
192 
193  // ----------------------- Standard services ------------------------------
196 
202  ATSolver2D( ConstAlias< Calculus > aCalculus, int aVerbose = 0 )
203  : ptrCalculus( aCalculus ),
204  primal_D0( *ptrCalculus ), primal_D1( *ptrCalculus ),
205  M01( *ptrCalculus ), M12( *ptrCalculus ), primal_AD2( *ptrCalculus ),
206  alpha_Id2( *ptrCalculus ), l_1_over_4e_Id0( *ptrCalculus ),
207  g2(), alpha_g2(), u2(), v0( *ptrCalculus ), former_v0( *ptrCalculus ),
208  l_1_over_4e( *ptrCalculus ), verbose( aVerbose )
209  {
210  if ( verbose >= 2 )
211  trace.info() << "[ATSolver::ATSolver] " << *ptrCalculus << std::endl;
212  initOperators();
213  const auto size2 = ptrCalculus->kFormLength( 2, PRIMAL );
214  for ( Index index = 0; index < size2; ++index) {
215  const auto& calculus_cell = ptrCalculus->getSCell( 2, PRIMAL, index );
216  surfel2idx[ calculus_cell ] = index;
217  }
218  }
219 
223  ATSolver2D() = delete;
224 
228  ~ATSolver2D() = default;
229 
234  ATSolver2D ( const ATSolver2D & other ) = default;
235 
240  ATSolver2D ( ATSolver2D && other ) = default;
241 
247  ATSolver2D & operator= ( const ATSolver2D & other ) = default;
248 
254  ATSolver2D & operator= ( ATSolver2D && other ) = default;
255 
258  Index size( const int order ) const
259  {
260  return ptrCalculus->kFormLength(order, PRIMAL);
261  }
262 
264 
265 
266  // ----------------------- Initialization services ------------------------------
267  public:
270 
286  template <typename VectorFieldInput, typename SurfelRangeConstIterator>
287  void
288  initInputVectorFieldU2( const VectorFieldInput& input,
289  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE,
290  bool normalize = false )
291  {
292  if ( verbose >= 1 )
293  trace.beginBlock( "[ATSolver2D::initInputVectorFieldU2] Initializing input data" );
294  ASSERT( ! input.empty() );
295  const Dimension N = input[ 0 ].size();
296  if ( verbose >= 2 ) trace.info() << "input g as " << N << " 2-forms." << std::endl;
297  g2 = std::vector<PrimalForm2>( N, PrimalForm2( *ptrCalculus ) );
298  alpha_g2 = std::vector<PrimalForm2>( N, PrimalForm2( *ptrCalculus ) );
299  Index i = 0;
300  for ( auto it = itB; it != itE; ++it, ++i )
301  {
302  Index idx = surfel2idx[ *it ];
303  for ( Dimension k = 0; k < N; ++k )
304  g2[ k ].myContainer( idx ) = input[ i ][ k ];
305  }
306  // u = g at the beginning
307  if ( verbose >= 2 )
308  trace.info() << "Unknown u[:] = g[:] at beginning." << std::endl;
309  u2 = g2;
310  // v = 1 at the beginning
311  if ( verbose >= 2 ) trace.info() << "Unknown v = 1" << std::endl;
312  v0 = PrimalForm0::ones(*ptrCalculus);
313  if ( verbose >= 1 ) trace.endBlock();
314  normalize_u2 = normalize;
315  }
316 
328  template <typename ScalarFieldInput, typename SurfelRangeConstIterator>
329  void
330  initInputScalarFieldU2( const ScalarFieldInput& input,
331  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
332  {
333  if ( verbose >= 1 )
334  trace.beginBlock( "[ATSolver2D::initInputVectorFieldU2] Initializing input data" );
335  ASSERT( ! input.empty() );
336  if ( verbose >= 2 ) trace.info() << "input g as one 2-form." << std::endl;
337  g2 = std::vector<PrimalForm2>( 1, PrimalForm2( *ptrCalculus ) );
338  alpha_g2 = std::vector<PrimalForm2>( 1, PrimalForm2( *ptrCalculus ) );
339  Index i = 0;
340  for ( auto it = itB; it != itE; ++it, ++i )
341  {
342  Index idx = surfel2idx[ *it ];
343  g2[ 0 ].myContainer( idx ) = input[ i ];
344  }
345  // u = g at the beginning
346  if ( verbose >= 2 )
347  trace.info() << "Unknown u[:] = g[:] at beginning." << std::endl;
348  u2 = g2;
349  // v = 1 at the beginning
350  if ( verbose >= 2 ) trace.info() << "Unknown v = 1" << std::endl;
351  v0 = PrimalForm0::ones(*ptrCalculus);
352  if ( verbose >= 1 ) trace.endBlock();
353  normalize_u2 = false;
354  }
355 
356 
369  template <typename ScalarVector>
370  Index
371  initInputVectorFieldU2( const std::map<Surfel,ScalarVector>& input,
372  bool normalize = false )
373  {
374  if ( verbose >= 1 ) trace.beginBlock( "[ATSolver2D::initVectorInput] Initializing input data" );
375  if ( verbose >= 2 ) trace.info() << "discontinuity 0-form v = 1." << std::endl;
376  const Dimension N = ScalarVector().size();
377  if ( verbose >= 2 ) trace.info() << "input g as " << N << " 2-forms." << std::endl;
378  g2 = std::vector<PrimalForm2>( N, PrimalForm2( *ptrCalculus ) );
379  alpha_g2 = std::vector<PrimalForm2>( N, PrimalForm2( *ptrCalculus ) );
380  const ScalarVector zero;
381  Index nbok = 0;
382  for ( Index index = 0; index < size(2); index++)
383  {
384  const SCell& cell = g2[ 0 ].getSCell( index );
385  const auto it = input.find( cell );
386  const auto n = ( it != input.end() ) ? it->second : zero;
387  nbok += ( it != input.end() ) ? 1 : 0;
388  for ( Dimension k = 0; k < N; ++k )
389  g2[ k ].myContainer( index ) = n[ k ];
390  }
391  // u = g at the beginning
392  if ( verbose >= 2 )
393  trace.info() << "Unknown u[:] = g[:] at beginning." << std::endl;
394  u2 = g2;
395  // v = 1 at the beginning
396  if ( verbose >= 2 ) trace.info() << "Unknown v = 1" << std::endl;
397  v0 = PrimalForm0::ones(*ptrCalculus);
398  if ( verbose >= 1 ) trace.endBlock();
399  normalize_u2 = normalize;
400  return nbok;
401  }
402 
409  template <typename Scalar>
410  Index
411  initInputScalarFieldU2( const std::map<Surfel,Scalar>& input )
412  {
413  if ( verbose >= 1 ) trace.beginBlock( "[ATSolver2D::initScalarInput] Initializing input data" );
414  if ( verbose >= 2 ) trace.info() << "discontinuity 0-form v = 1." << std::endl;
415  v0 = PrimalForm0::ones(*ptrCalculus);
416  if ( verbose >= 2 ) trace.info() << "input g as one 2-form." << std::endl;
417  g2 = std::vector<PrimalForm2>( 1, PrimalForm2( *ptrCalculus ) );
418  alpha_g2 = std::vector<PrimalForm2>( 1, PrimalForm2( *ptrCalculus ) );
419  const Scalar zero;
420  Index nbok = 0;
421  for ( Index index = 0; index < size(2); index++)
422  {
423  const SCell& cell = g2[ 0 ].getSCell( index );
424  const auto it = input.find( cell );
425  const auto n = ( it != input.end() ) ? *it : zero;
426  nbok += ( it != input.end() ) ? 1 : 0;
427  g2[ 0 ].myContainer( index ) = n;
428  }
429  // u = g at the beginning
430  if ( verbose >= 2 ) trace.info() << "Unknown u[:] = g[:] at beginning." << std::endl;
431  u2 = g2;
432  // v = 1 at the beginning
433  if ( verbose >= 2 ) trace.info() << "Unknown v = 1" << std::endl;
434  v0 = PrimalForm0::ones(*ptrCalculus);
435  if ( verbose >= 1 ) trace.endBlock();
436  normalize_u2 = false;
437  return nbok;
438  }
439 
444  void setUp( double a, double l )
445  {
446  const Dimension N = g2.size();
447  alpha = a;
448  lambda = l;
449  alpha_Id2 = alpha * ptrCalculus->template identity<2, PRIMAL>();
450  for ( Dimension k = 0; k < N; ++k )
451  alpha_g2[ k ] = alpha * g2[ k ];
452  }
453 
463  void setUp( double a, double l, const std::map<Surfel,Scalar>& weights )
464  {
465  const Dimension N = g2.size();
466  alpha = a;
467  lambda = l;
468  PrimalForm2 w_form( *ptrCalculus );
469  if ( verbose >= 2 )
470  trace.info() << "Using variable weights for fitting (alpha term)" << std::endl;
471  for ( Dimension k = 0; k < N; ++k )
472  alpha_g2[ k ] = alpha * g2[ k ];
473  for ( Index index = 0; index < size( 2 ); index++)
474  {
475  const SCell& cell = g2[ 0 ].getSCell( index );
476  const Scalar& w = weights[ cell ];
477  w_form.myContainer( index ) = w;
478  for ( Dimension k = 0; k < N; ++k )
479  alpha_g2[ k ].myContainer( index ) *= w;
480  }
481  alpha_Id2 = alpha * diagonal( w_form );
482  }
483 
498  template <typename AlphaWeights, typename SurfelRangeConstIterator>
499  void setUp( double a, double l,
500  const AlphaWeights& weights,
501  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
502  {
503  const Dimension N = g2.size();
504  alpha = a;
505  lambda = l;
506  PrimalForm2 w_form( *ptrCalculus );
507  if ( verbose >= 2 )
508  trace.info() << "Using variable weights for fitting (alpha term)" << std::endl;
509  for ( Dimension k = 0; k < N; ++k )
510  alpha_g2[ k ] = alpha * g2[ k ];
511  Index i = 0;
512  for ( auto it = itB; it != itE; ++it, ++i )
513  {
514  const Index idx = surfel2idx[ *it ];
515  const Scalar w = weights[ i ];
516  w_form.myContainer( idx ) = w;
517  for ( Dimension k = 0; k < N; ++k )
518  alpha_g2[ k ].myContainer( idx ) *= w;
519  }
520  alpha_Id2 = alpha * diagonal( w_form );
521  }
522 
525  void setEpsilon( double e )
526  {
527  epsilon = e;
528  l_1_over_4e = (lambda/4./epsilon)*PrimalForm0::ones(*ptrCalculus);
529  l_1_over_4e_Id0 = (lambda/4./epsilon)*ptrCalculus->template identity<0, PRIMAL>();
530  }
531 
533 
534  // ----------------------- Optimization services ------------------------------
535  public:
538 
546  bool solveOneAlternateStep()
547  {
548  bool solve_ok = true;
549  if ( verbose >= 1 ) trace.beginBlock("Solving for u as a 2-form");
550  PrimalForm1 v1_squared = M01*v0;
551  v1_squared.myContainer.array() = v1_squared.myContainer.array().square();
552  const PrimalIdentity2 ope_u2 = alpha_Id2
553  + primal_AD2.transpose() * dec_helper::diagonal( v1_squared ) * primal_AD2;
554 
555  if ( verbose >= 2 ) trace.info() << "Prefactoring matrix U associated to u" << std::endl;
556  SolverU2 solver_u2;
557  solver_u2.compute( ope_u2 );
558  for ( Dimension d = 0; d < u2.size(); ++d )
559  {
560  if ( verbose >= 2 ) trace.info() << "Solving U u[" << d << "] = a g[" << d << "]" << std::endl;
561  u2[ d ] = solver_u2.solve( alpha_g2[ d ] );
562  if ( verbose >= 2 ) trace.info() << " => " << ( solver_u2.isValid() ? "OK" : "ERROR" )
563  << " " << solver_u2.myLinearAlgebraSolver.info() << std::endl;
564  solve_ok = solve_ok && solver_u2.isValid();
565  }
566  if ( normalize_u2 ) normalizeU2();
567  if ( verbose >= 1 ) trace.endBlock();
568  if ( verbose >= 1 ) trace.beginBlock("Solving for v");
569  former_v0 = v0;
570  PrimalForm1 squared_norm_d_u2 = PrimalForm1::zeros(*ptrCalculus);
571  for ( Dimension d = 0; d < u2.size(); ++d )
572  squared_norm_d_u2.myContainer.array() += (primal_AD2 * u2[ d ] ).myContainer.array().square();
573  trace.info() << "build metric u2" << std::endl;
574  const PrimalIdentity0 ope_v0 = l_1_over_4e_Id0
575  + (lambda * epsilon) * primal_D0.transpose() * primal_D0
576  + M01.transpose() * dec_helper::diagonal( squared_norm_d_u2 ) * M01;
577 
578  if ( verbose >= 2 ) trace.info() << "Prefactoring matrix V associated to v" << std::endl;
579  SolverV0 solver_v0;
580  solver_v0.compute( ope_v0 );
581  if ( verbose >= 2 ) trace.info() << "Solving V v = l/4e * 1" << std::endl;
582  v0 = solver_v0.solve( l_1_over_4e );
583  if ( verbose >= 2 ) trace.info() << " => " << ( solver_v0.isValid() ? "OK" : "ERROR" )
584  << " " << solver_v0.myLinearAlgebraSolver.info() << std::endl;
585  solve_ok = solve_ok && solver_v0.isValid();
586  if ( verbose >= 1 ) trace.endBlock();
587  return solve_ok;
588  }
589 
605  bool solveForEpsilon( double eps,
606  double n_oo_max = 1e-4,
607  unsigned int iter_max = 10 )
608  {
609  bool ok = true;
610  if ( verbose >= 1 ) {
611  std::ostringstream sstr;
612  sstr << "******* Solving AT for epsilon = " << eps << " **********";
613  trace.beginBlock( sstr.str() );
614  }
615  setEpsilon( eps );
616  for ( unsigned int i = 0; i < iter_max; ++i )
617  {
618  if ( verbose >= 1 )
619  trace.info() << "---------- Iteration "
620  << i << "/" << iter_max << " ---------------" << std::endl;
621  solveOneAlternateStep( );
622  auto norms_v = checkV0();
623  auto diffs_v = diffV0();
624  if ( verbose >= 1 ) {
625  trace.info() << "Variation |v^k+1 - v^k|_oo = " << std::get<0>( diffs_v )
626  << std::endl;
627  if ( verbose >= 2 ) {
628  trace.info() << "Variation |v^k+1 - v^k|_2 = " << std::get<1>( diffs_v )
629  << std::endl;
630  trace.info() << "Variation |v^k+1 - v^k|_1 = " << std::get<2>( diffs_v )
631  << std::endl;
632  }
633  }
634  if ( std::get<0>( diffs_v ) < 1e-4 ) break;
635  }
636  if ( verbose >= 1 ) trace.endBlock();
637  return ok;
638  }
639 
656  bool solveGammaConvergence( double eps1 = 2.0,
657  double eps2 = 0.25,
658  double epsr = 2.0,
659  bool compute_smallest_epsilon_map = false,
660  double n_oo_max = 1e-4,
661  unsigned int iter_max = 10 )
662  {
663  bool ok = true;
664  if ( epsr <= 1.0 ) epsr = 2.0;
665  if ( verbose >= 1 )
666  trace.beginBlock( "#### Solve AT by Gamma-convergence ##########" );
667  if ( compute_smallest_epsilon_map ) smallest_epsilon_map.clear();
668  for ( double eps = eps1; eps >= eps2; eps /= epsr )
669  {
670  solveForEpsilon( eps, n_oo_max, iter_max );
671  if ( compute_smallest_epsilon_map )
672  updateSmallestEpsilonMap( 0.5 );
673  }
674  if ( verbose >= 1 )
675  trace.endBlock();
676  return ok;
677  }
678 
683  void normalizeU2()
684  {
685  for ( Index index = 0; index < size( 2 ); index++)
686  {
687  double n2 = 0.0;
688  for ( unsigned int d = 0; d < u2.size(); ++d )
689  n2 += u2[ d ].myContainer( index ) * u2[ d ].myContainer( index );
690  double norm = sqrt( n2 );
691  if (norm == 0.0) continue;
692  for ( unsigned int d = 0; d < u2.size(); ++d )
693  u2[ d ].myContainer( index ) /= norm;
694  }
695  }
696 
701  std::tuple<double,double,double> diffV0() const
702  {
703  PrimalForm0 delta = v0 - former_v0;
704  delta.myContainer = delta.myContainer.cwiseAbs();
705  const double n_oo = delta.myContainer.maxCoeff();
706  const double n_2 = std::sqrt(delta.myContainer.squaredNorm()/delta.myContainer.size());
707  const double n_1 = delta.myContainer.mean();
708  return std::make_tuple( n_oo, n_2, n_1 );
709  }
710 
712 
713 
714  // ----------------------- Access services ------------------------------
715  public:
718 
722  std::tuple<double,double,double> checkV0() const
723  {
724  const double m1 = v0.myContainer.minCoeff();
725  const double m2 = v0.myContainer.maxCoeff();
726  const double ma = v0.myContainer.mean();
727  if ( verbose >= 1 )
728  trace.info() << "0-form v (should be in [0,1]): min=" << m1 << " avg=" << ma << " max=" << m2 << std::endl;
729  return std::make_tuple( m1, m2, ma );
730  }
731 
732 
734  PrimalForm0 getV0() const
735  {
736  return v0;
737  }
738 
740  PrimalForm1 getV1() const
741  {
742  return M01*v0;
743  }
744 
746  PrimalForm2 getV2() const
747  {
748  return M12*M01*v0;
749  }
750 
753  PrimalForm2 getU2( Dimension k ) const
754  {
755  return u2[ k ];
756  }
757 
769  template <typename VectorFieldOutput, typename SurfelRangeConstIterator>
770  void
771  getOutputVectorFieldU2( VectorFieldOutput& output,
772  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
773  {
774  const Dimension N = u2.size();
775  Index i = 0;
776  for ( auto it = itB; it != itE; ++it, ++i )
777  {
778  Index idx = surfel2idx[ *it ];
779  ASSERT( output[ i ].size() >= N );
780  for ( Dimension k = 0; k < N; ++k )
781  output[ i ][ k ] = u2[ k ].myContainer( idx );
782  }
783  }
784 
795  template <typename ScalarFieldOutput, typename SurfelRangeConstIterator>
796  void
797  getOutputScalarFieldU2( ScalarFieldOutput& output,
798  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
799  {
800  ASSERT( u2.size() == 1 && "[ATSolver2D::getOutputScalarFieldU2] "
801  "You try to output a scalar field from a vector field." );
802  Index i = 0;
803  for ( auto it = itB; it != itE; ++it, ++i )
804  {
805  Index idx = surfel2idx[ *it ];
806  output[ i ] = u2[ 0 ].myContainer( idx );
807  }
808  }
809 
821  template <typename ScalarFieldOutput, typename CellRangeConstIterator>
822  void
823  getOutputScalarFieldV0( ScalarFieldOutput& output,
824  CellRangeConstIterator itB, CellRangeConstIterator itE,
825  CellOutputPolicy policy = CellOutputPolicy::Average )
826  {
827  const KSpace& K = ptrCalculus->myKSpace;
828  const Dimension k = K.uDim( *itB );
829  ASSERT( k <= 2 );
830  Index i = 0;
831  if ( k == 0 )
832  {
833  for ( auto it = itB; it != itE; ++it, ++i )
834  {
835  const Cell pointel = *it;
836  const Index idx = ptrCalculus->getCellIndex( pointel );
837  output[ i ] = v0.myContainer( idx );
838  }
839  }
840  else if ( k == 1 )
841  {
842  for ( auto it = itB; it != itE; ++it, ++i )
843  {
844  const Cell linel = *it;
845  const Dimension d = * K.uDirs( linel );
846  const Cell p0 = K.uIncident( linel, d, false );
847  const Cell p1 = K.uIncident( linel, d, true );
848  const Index idx0 = ptrCalculus->getCellIndex( p0 );
849  const Index idx1 = ptrCalculus->getCellIndex( p1 );
850  switch (policy) {
851  case CellOutputPolicy::Average: output[ i ] = 0.5 * ( v0.myContainer( idx0 ) + v0.myContainer( idx1 ) );
852  break;
853  case CellOutputPolicy::Minimum: output[ i ] = std::min( v0.myContainer( idx0 ), v0.myContainer( idx1 ) );
854  break;
855  case CellOutputPolicy::Maximum: output[ i ] = std::max( v0.myContainer( idx0 ), v0.myContainer( idx1 ) );
856  break;
857  }
858  }
859  }
860  else if ( k == 2 )
861  {
862  for ( auto it = itB; it != itE; ++it, ++i )
863  {
864  const Cell face = *it;
865  const Dimension d = * K.uDirs( face );
866  const Cell l0 = K.uIncident( face, d, false );
867  const Cell l1 = K.uIncident( face, d, true );
868  const Dimension j = * K.uDirs( l0 );
869  const Cell p00 = K.uIncident( l0, j, false );
870  const Cell p01 = K.uIncident( l0, j, true );
871  const Cell p10 = K.uIncident( l1, j, false );
872  const Cell p11 = K.uIncident( l1, j, true );
873  const Index idx00 = ptrCalculus->getCellIndex( p00 );
874  const Index idx01 = ptrCalculus->getCellIndex( p01 );
875  const Index idx10 = ptrCalculus->getCellIndex( p10 );
876  const Index idx11 = ptrCalculus->getCellIndex( p11 );
877  switch (policy) {
878  case CellOutputPolicy::Average:
879  output[ i ] = 0.25 * ( v0.myContainer( idx00 ) + v0.myContainer( idx01 )
880  + v0.myContainer( idx10 ) + v0.myContainer( idx11 ) );
881  break;
882  case CellOutputPolicy::Minimum:
883  output[ i ] = std::min( std::min( v0.myContainer( idx00 ), v0.myContainer( idx01 ) ),
884  std::min( v0.myContainer( idx10 ), v0.myContainer( idx11 ) ) );
885  break;
886  case CellOutputPolicy::Maximum:
887  output[ i ] = std::max( std::max( v0.myContainer( idx00 ), v0.myContainer( idx01 ) ),
888  std::max( v0.myContainer( idx10 ), v0.myContainer( idx11 ) ) );
889  break;
890  }
891  }
892  }
893  }
894 
902  void updateSmallestEpsilonMap( const double threshold = .5 )
903  {
904  const KSpace& K = ptrCalculus->myKSpace;
905  for ( const SCell surfel : ptrCalculus->template getIndexedSCells<2, PRIMAL>() )
906  {
907  const Cell face = K.unsigns( surfel );
908  const Dimension k1 = * K.uDirs( face );
909  const Cell l0 = K.uIncident( face, k1, false );
910  const Cell l1 = K.uIncident( face, k1, true );
911  const Dimension k2 = * K.uDirs( l0 );
912  const Cell ll0 = K.uIncident( face, k2, false );
913  const Cell ll1 = K.uIncident( face, k2, true );
914  const Cell p00 = K.uIncident( l0, k2, false );
915  const Cell p01 = K.uIncident( l0, k2, true );
916  const Cell p10 = K.uIncident( l1, k2, false );
917  const Cell p11 = K.uIncident( l1, k2, true );
918 
919  std::vector<double> features( 4 );
920  features[ 0 ] = v0.myContainer( ptrCalculus->getCellIndex( p00 ) );
921  features[ 1 ] = v0.myContainer( ptrCalculus->getCellIndex( p01 ) );
922  features[ 2 ] = v0.myContainer( ptrCalculus->getCellIndex( p10 ) );
923  features[ 3 ] = v0.myContainer( ptrCalculus->getCellIndex( p11 ) );
924  std::sort( features.begin(), features.end() );
925 
926  if ( features[ 1 ] <= threshold )
927  {
928  auto it = smallest_epsilon_map.find( surfel );
929  if ( it != smallest_epsilon_map.end() )
930  it->second = std::min( epsilon, it->second );
931  else smallest_epsilon_map[ surfel ] = epsilon;
932  }
933  }
934  }
935 
937 
938 
939  // ----------------------- Interface --------------------------------------
940  public:
943 
948  void selfDisplay ( std::ostream & out ) const
949  {
950  auto cv = checkV0();
951  out << "[ATSolver2D] v is between min/avg/max:"
952  << std::get<0>(cv) << "/"
953  << std::get<1>(cv) << "/"
954  << std::get<2>(cv) << std::endl;
955  }
956 
961  bool isValid() const
962  {
963  return true;
964  }
965 
967 
968  // ------------------------- Hidden services ------------------------------
969  protected:
970 
973 
975  void initOperators()
976  {
977  if ( verbose >= 1 ) trace.beginBlock( "[ATSolver2D::initOperators] Solver initialization" );
978  if ( verbose >= 2 ) trace.info() << "derivative of primal 0-forms: primal_D0" << std::endl;
979  primal_D0 = ptrCalculus->template derivative<0,PRIMAL>();
980  if ( verbose >= 2 ) trace.info() << "derivative of primal 1-forms: primal_D1" << std::endl;
981  primal_D1 = ptrCalculus->template derivative<1,PRIMAL>();
982  if ( verbose >= 2 ) trace.info() << "antiderivative of primal 2-forms: primal_AD2" << std::endl;
983  primal_AD2 = ptrCalculus->template antiderivative<2,PRIMAL>();
984  if ( verbose >= 2 ) trace.info() << "vertex to edge average operator: M01" << std::endl;
985  M01 = primal_D0;
986  M01.myContainer = .5 * M01.myContainer.cwiseAbs();
987  if ( verbose >= 2 ) trace.info() << "edge to face average operator: M12" << std::endl;
988  M12 = primal_D1;
989  M12.myContainer = .25 * M12.myContainer.cwiseAbs();
990  if ( verbose >= 1 ) trace.endBlock();
991  }
992 
994 
995  // ------------------------- Internals ------------------------------------
996  private:
997 
998  }; // end of class ATSolver2D
999 
1000 
1007  template <typename T>
1008  std::ostream&
1009  operator<< ( std::ostream & out, const ATSolver2D<T> & object );
1010 
1011 } // namespace DGtal
1012 
1013 
1015 // Includes inline functions.
1016 
1017 // //
1019 
1020 #endif // !defined ATSolver2D_h
1021 
1022 #undef ATSolver2D_RECURSES
1023 #endif // else defined(ATSolver2D_RECURSES)
DGtal::ATSolver2D::checkV0
std::tuple< double, double, double > checkV0() const
Definition: ATSolver2D.h:722
DGtal::ATSolver2D::Minimum
@ Minimum
compute minimum value at cell vertices,
Definition: ATSolver2D.h:104
DGtal::ATSolver2D::PrimalForm0
Calculus::PrimalForm0 PrimalForm0
Definition: ATSolver2D.h:117
DGtal::ATSolver2D::getOutputVectorFieldU2
void getOutputVectorFieldU2(VectorFieldOutput &output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:771
DGtal::ATSolver2D::alpha_Id2
PrimalIdentity2 alpha_Id2
The alpha-weighted identity operator for primal 2-forms (stored for performance)
Definition: ATSolver2D.h:156
DGtal::ATSolver2D::LinearAlgebraSolver
EigenLinearAlgebraBackend::SolverSimplicialLDLT LinearAlgebraSolver
Definition: ATSolver2D.h:138
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition: DiscreteExteriorCalculus.h:97
DGtal::ATSolver2D::setUp
void setUp(double a, double l)
Definition: ATSolver2D.h:444
DGtal::ConstAlias
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:186
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::ATSolver2D::getV2
PrimalForm2 getV2() const
Definition: ATSolver2D.h:746
DGtal::KForm::ones
static KForm< TCalculus, order, duality > ones(ConstAlias< Calculus > calculus)
DGtal::KhalimskySpaceND::uDim
Dimension uDim(const Cell &p) const
Return the dimension of the cell [p].
DGtal::ATSolver2D::PrimalForm1
Calculus::PrimalForm1 PrimalForm1
Definition: ATSolver2D.h:118
DGtal::ATSolver2D::M01
PrimalDerivative0 M01
The primal vertex to edge average operator.
Definition: ATSolver2D.h:150
DGtal::ATSolver2D::l_1_over_4e_Id0
PrimalIdentity0 l_1_over_4e_Id0
The 1/(4epsilon)-weighted identity operator for primal 0-forms (stored for performance)
Definition: ATSolver2D.h:158
DGtal::DiscreteExteriorCalculus::Index
LinearAlgebraBackend::DenseVector::Index Index
Definition: DiscreteExteriorCalculus.h:107
DGtal::ATSolver2D::g2
std::vector< PrimalForm2 > g2
The N-array of input primal 2-forms g.
Definition: ATSolver2D.h:160
max
int max(int a, int b)
Definition: testArithmeticalDSS.cpp:1108
DGtal::ATSolver2D::PrimalAntiderivative2
Calculus::PrimalAntiderivative2 PrimalAntiderivative2
Definition: ATSolver2D.h:126
DGtal::ATSolver2D::PrimalIdentity1
Calculus::PrimalIdentity1 PrimalIdentity1
Definition: ATSolver2D.h:121
DGtal::ATSolver2D::PrimalDerivative0
Calculus::PrimalDerivative0 PrimalDerivative0
Definition: ATSolver2D.h:123
DGtal::ATSolver2D::alpha_g2
std::vector< PrimalForm2 > alpha_g2
The alpha-weighted N-array of input primal 2-forms g.
Definition: ATSolver2D.h:162
DGtal::ATSolver2D::initInputScalarFieldU2
void initInputScalarFieldU2(const ScalarFieldInput &input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:330
DGtal::ATSolver2D::getU2
PrimalForm2 getU2(Dimension k) const
Definition: ATSolver2D.h:753
DGtal::ATSolver2D::M12
PrimalDerivative1 M12
The primal edge to face average operator.
Definition: ATSolver2D.h:152
DGtal::ATSolver2D::normalize_u2
bool normalize_u2
Indicates whether to normalize U (unit norm) at each iteration or not.
Definition: ATSolver2D.h:189
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::ATSolver2D::SCell
KSpace::SCell SCell
Definition: ATSolver2D.h:110
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::ATSolver2D::getOutputScalarFieldU2
void getOutputScalarFieldU2(ScalarFieldOutput &output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:797
DGtal::DiscreteExteriorCalculusSolver::compute
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::ATSolver2D::Surfel
KSpace::Surfel Surfel
Definition: ATSolver2D.h:112
DGtal::ATSolver2D::ptrCalculus
CountedConstPtrOrConstPtr< Calculus > ptrCalculus
A smart (or not) pointer to a calculus object.
Definition: ATSolver2D.h:144
DGtal::DiscreteExteriorCalculusSolver::myLinearAlgebraSolver
LinearAlgebraSolver myLinearAlgebraSolver
Definition: DiscreteExteriorCalculusSolver.h:100
DGtal::ATSolver2D::PrimalHodge0
Calculus::PrimalHodge0 PrimalHodge0
Definition: ATSolver2D.h:127
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::ATSolver2D::isValid
bool isValid() const
Definition: ATSolver2D.h:961
DGtal::ATSolver2D::alpha
double alpha
Definition: ATSolver2D.h:181
DGtal::LinearOperator
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
Definition: LinearOperator.h:68
DGtal::ATSolver2D::operator=
ATSolver2D & operator=(const ATSolver2D &other)=default
DGtal::ATSolver2D::PrimalIdentity2
Calculus::PrimalIdentity2 PrimalIdentity2
Definition: ATSolver2D.h:122
DGtal::ATSolver2D::CellOutputPolicy
CellOutputPolicy
Definition: ATSolver2D.h:103
DGtal::ATSolver2D::PrimalIdentity0
Calculus::PrimalIdentity0 PrimalIdentity0
Definition: ATSolver2D.h:120
DGtal::ATSolver2D::PrimalHodge2
Calculus::PrimalHodge2 PrimalHodge2
Definition: ATSolver2D.h:129
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:208
DGtal::ATSolver2D::u2
std::vector< PrimalForm2 > u2
The N-array of regularized primal 2-forms u.
Definition: ATSolver2D.h:164
DGtal::ATSolver2D::setEpsilon
void setEpsilon(double e)
Definition: ATSolver2D.h:525
DGtal::ATSolver2D::epsilon
double epsilon
Definition: ATSolver2D.h:187
DGtal::ATSolver2D::solveGammaConvergence
bool solveGammaConvergence(double eps1=2.0, double eps2=0.25, double epsr=2.0, bool compute_smallest_epsilon_map=false, double n_oo_max=1e-4, unsigned int iter_max=10)
Definition: ATSolver2D.h:656
DGtal::PointVector::Component
TEuclideanRing Component
Type for Vector elements.
Definition: PointVector.h:614
DGtal::operator<<
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::ATSolver2D::primal_D0
PrimalDerivative0 primal_D0
the derivative operator for primal 0-forms
Definition: ATSolver2D.h:146
DGtal::SpaceND
Definition: SpaceND.h:95
DGtal::PRIMAL
@ PRIMAL
Definition: Duality.h:61
DGtal::ATSolver2D::initOperators
void initOperators()
Initializes the operators.
Definition: ATSolver2D.h:975
DGtal::DiscreteExteriorCalculusSolver::isValid
bool isValid() const
DGtal::Trace::info
std::ostream & info()
DGtal::KhalimskySpaceND::dimension
static const constexpr Dimension dimension
Definition: KhalimskySpaceND.h:430
DGtal::ATSolver2D::Domain
HyperRectDomain< Space > Domain
Definition: ATSolver2D.h:113
DGtal::LinearOperator::transpose
TransposedLinearOperator transpose() const
DGtal::ATSolver2D::solveOneAlternateStep
bool solveOneAlternateStep()
Definition: ATSolver2D.h:546
DGtal::ATSolver2D::Surfel2IndexMap
KSpace::template SurfelMap< Index >::Type Surfel2IndexMap
Definition: ATSolver2D.h:130
DGtal::CountedConstPtrOrConstPtr
Aim: Smart or simple const pointer on T. It can be a smart pointer based on reference counts or a sim...
Definition: CountedConstPtrOrConstPtr.h:94
DGtal::ATSolver2D::size
Index size(const int order) const
Definition: ATSolver2D.h:258
DGtal::ATSolver2D::setUp
void setUp(double a, double l, const AlphaWeights &weights, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:499
DGtal::ATSolver2D::Average
@ Average
compute average values at cell vertices
Definition: ATSolver2D.h:103
DGtal::ATSolver2D::former_v0
PrimalForm0 former_v0
The primal 0-form v at the previous iteration.
Definition: ATSolver2D.h:169
DGtal::ATSolver2D
Aim: This class solves Ambrosio-Tortorelli functional on a two-dimensional digital space (a 2D grid o...
Definition: ATSolver2D.h:90
DGtal::ATSolver2D::Self
ATSolver2D< KSpace, LinearAlgebra > Self
Definition: ATSolver2D.h:97
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::KForm::myContainer
Container myContainer
Definition: KForm.h:131
index
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition: testBits.cpp:44
DGtal::ATSolver2D::Maximum
@ Maximum
compute maximum value at cell vertices
Definition: ATSolver2D.h:105
DGtal::ATSolver2D::updateSmallestEpsilonMap
void updateSmallestEpsilonMap(const double threshold=.5)
Definition: ATSolver2D.h:902
DGtal::ATSolver2D::Index
Calculus::Index Index
Definition: ATSolver2D.h:116
DGtal::ATSolver2D::Calculus
DiscreteExteriorCalculus< 2, dimension, LinearAlgebra > Calculus
Definition: ATSolver2D.h:114
DGtal::ATSolver2D::Scalar
RealVector::Component Scalar
Definition: ATSolver2D.h:109
DGtal::ATSolver2D::normalizeU2
void normalizeU2()
Definition: ATSolver2D.h:683
DGtal::ATSolver2D::smallest_epsilon_map
SmallestEpsilonMap smallest_epsilon_map
Definition: ATSolver2D.h:178
DGtal::ATSolver2D::surfel2idx
Surfel2IndexMap surfel2idx
Definition: ATSolver2D.h:175
DGtal::ATSolver2D::PrimalHodge1
Calculus::PrimalHodge1 PrimalHodge1
Definition: ATSolver2D.h:128
DGtal::ATSolver2D::SolverU2
DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL > SolverU2
Definition: ATSolver2D.h:139
DGtal::DiscreteExteriorCalculusSolver::solve
SolutionKForm solve(const InputKForm &input_kform) const
DGtal::ATSolver2D::primal_AD2
PrimalAntiderivative2 primal_AD2
The antiderivative of primal 2-forms.
Definition: ATSolver2D.h:154
DGtal::KhalimskySpaceND::uDirs
DirIterator uDirs(const Cell &p) const
Given an unsigned cell [p], returns an iterator to iterate over each coordinate the cell spans.
DGtal::ATSolver2D::initInputVectorFieldU2
Index initInputVectorFieldU2(const std::map< Surfel, ScalarVector > &input, bool normalize=false)
Definition: ATSolver2D.h:371
DGtal::ATSolver2D::Space
KSpace::Space Space
Definition: ATSolver2D.h:107
DGtal::ATSolver2D::diffV0
std::tuple< double, double, double > diffV0() const
Definition: ATSolver2D.h:701
DGtal::ATSolver2D::ATSolver2D
ATSolver2D(ConstAlias< Calculus > aCalculus, int aVerbose=0)
Definition: ATSolver2D.h:202
DGtal::EigenLinearAlgebraBackend::SolverSimplicialLDLT
Eigen::SimplicialLDLT< SparseMatrix > SolverSimplicialLDLT
Definition: EigenSupport.h:107
DGtal::ATSolver2D::SolverV0
DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL > SolverV0
Definition: ATSolver2D.h:140
DGtal::ATSolver2D::~ATSolver2D
~ATSolver2D()=default
DGtal::ATSolver2D::selfDisplay
void selfDisplay(std::ostream &out) const
Definition: ATSolver2D.h:948
DGtal::KForm::zeros
static KForm< TCalculus, order, duality > zeros(ConstAlias< Calculus > calculus)
DGtal::ATSolver2D::getV1
PrimalForm1 getV1() const
Definition: ATSolver2D.h:740
DGtal::ATSolver2D::ATSolver2D
ATSolver2D()=delete
DGtal::ATSolver2D::SmallestEpsilonMap
KSpace::template SurfelMap< double >::Type SmallestEpsilonMap
Definition: ATSolver2D.h:115
DGtal::ATSolver2D::l_1_over_4e
PrimalForm0 l_1_over_4e
The primal 0-form lambda/(4epsilon) (stored for performance)
Definition: ATSolver2D.h:171
DGtal::KhalimskySpaceND::uIncident
Cell uIncident(const Cell &c, Dimension k, bool up) const
Return the forward or backward unsigned cell incident to [c] along axis [k], depending on [up].
DGtal::ATSolver2D::PrimalDerivative1
Calculus::PrimalDerivative1 PrimalDerivative1
Definition: ATSolver2D.h:124
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
DGtal::ATSolver2D::solveForEpsilon
bool solveForEpsilon(double eps, double n_oo_max=1e-4, unsigned int iter_max=10)
Definition: ATSolver2D.h:605
DGtal::ATSolver2D::LinearAlgebra
TLinearAlgebra LinearAlgebra
Definition: ATSolver2D.h:96
DGtal::DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL >
DGtal::ATSolver2D::PrimalAntiderivative1
Calculus::PrimalAntiderivative1 PrimalAntiderivative1
Definition: ATSolver2D.h:125
DGtal::KForm
Aim: KForm represents discrete kforms in the dec package.
Definition: KForm.h:65
DGtal::ATSolver2D::initInputScalarFieldU2
Index initInputScalarFieldU2(const std::map< Surfel, Scalar > &input)
Definition: ATSolver2D.h:411
DGtal::ATSolver2D::KSpace
TKSpace KSpace
Definition: ATSolver2D.h:95
DGtal::KhalimskySpaceND::unsigns
Cell unsigns(const SCell &p) const
Creates an unsigned cell from a signed one.
DGtal::ATSolver2D::v0
PrimalForm0 v0
Definition: ATSolver2D.h:167
DGtal::ATSolver2D::PrimalForm2
Calculus::PrimalForm2 PrimalForm2
Definition: ATSolver2D.h:119
DGtal::ATSolver2D::verbose
int verbose
Tells the verbose level.
Definition: ATSolver2D.h:191
DGtal::ATSolver2D::dimension
static const Dimension dimension
Definition: ATSolver2D.h:99
DGtal::ATSolver2D::lambda
double lambda
Definition: ATSolver2D.h:184
DGtal::ATSolver2D::setUp
void setUp(double a, double l, const std::map< Surfel, Scalar > &weights)
Definition: ATSolver2D.h:463
DGtal::ATSolver2D::Cell
KSpace::Cell Cell
Definition: ATSolver2D.h:111
DGtal::ATSolver2D::initInputVectorFieldU2
void initInputVectorFieldU2(const VectorFieldInput &input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE, bool normalize=false)
Definition: ATSolver2D.h:288
DGtal::ATSolver2D::primal_D1
PrimalDerivative1 primal_D1
the derivative operator for primal 1-forms
Definition: ATSolver2D.h:148
DGtal::ATSolver2D::getOutputScalarFieldV0
void getOutputScalarFieldV0(ScalarFieldOutput &output, CellRangeConstIterator itB, CellRangeConstIterator itE, CellOutputPolicy policy=CellOutputPolicy::Average)
Definition: ATSolver2D.h:823
DGtal::ATSolver2D::RealVector
Space::RealVector RealVector
Definition: ATSolver2D.h:108
DGtal::ATSolver2D::getV0
PrimalForm0 getV0() const
Definition: ATSolver2D.h:734
DGtal::KhalimskyCell< dim, Integer >
DGtal::dec_helper::diagonal
DGtal::LinearOperator< Calculus, dim, duality, dim, duality > diagonal(const DGtal::KForm< Calculus, dim, duality > &kform)
Definition: DECHelpers.h:60
DGtal::LinearOperator::myContainer
Container myContainer
Definition: LinearOperator.h:117
Generated on Mon Jun 20 2022 18:23:18 for DGtal by   doxygen 1.8.17