DGtal  1.5.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)
35 #define ATSolver2D_RECURSES
36 
37 #if !defined ATSolver2D_h
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;
98 
99  static const Dimension dimension = KSpace::dimension;
100 
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;
115  typedef typename KSpace::template SurfelMap<double>::Type SmallestEpsilonMap;
116  typedef typename Calculus::Index Index;
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;
141 
142  protected:
160  std::vector<PrimalForm2> g2;
162  std::vector<PrimalForm2> alpha_g2;
164  std::vector<PrimalForm2> u2;
172 
173  public:
174  // The map Surfel -> Index that gives the index of the surfel in 2-forms.
181  double alpha;
184  double lambda;
187  double epsilon;
191  int verbose;
192 
193  // ----------------------- Standard services ------------------------------
196 
202  ATSolver2D( ConstAlias< Calculus > aCalculus, int aVerbose = 0 )
203  : ptrCalculus( aCalculus ),
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;
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;
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;
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;
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;
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 = (Dimension)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;
529  l_1_over_4e_Id0 = (lambda/4./epsilon)*ptrCalculus->template identity<0, PRIMAL>();
530  }
531 
533 
534  // ----------------------- Optimization services ------------------------------
535  public:
538 
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
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  (void)n_oo_max;//paramerter not used
610 
611  bool ok = true;
612  if ( verbose >= 1 ) {
613  std::ostringstream sstr;
614  sstr << "******* Solving AT for epsilon = " << eps << " **********";
615  trace.beginBlock( sstr.str() );
616  }
617  setEpsilon( eps );
618  for ( unsigned int i = 0; i < iter_max; ++i )
619  {
620  if ( verbose >= 1 )
621  trace.info() << "---------- Iteration "
622  << i << "/" << iter_max << " ---------------" << std::endl;
624  auto diffs_v = diffV0();
625  if ( verbose >= 1 ) {
626  trace.info() << "Variation |v^k+1 - v^k|_oo = " << std::get<0>( diffs_v )
627  << std::endl;
628  if ( verbose >= 2 ) {
629  trace.info() << "Variation |v^k+1 - v^k|_2 = " << std::get<1>( diffs_v )
630  << std::endl;
631  trace.info() << "Variation |v^k+1 - v^k|_1 = " << std::get<2>( diffs_v )
632  << std::endl;
633  }
634  }
635  if ( std::get<0>( diffs_v ) < 1e-4 ) break;
636  }
637  if ( verbose >= 1 ) trace.endBlock();
638  return ok;
639  }
640 
657  bool solveGammaConvergence( double eps1 = 2.0,
658  double eps2 = 0.25,
659  double epsr = 2.0,
660  bool compute_smallest_epsilon_map = false,
661  double n_oo_max = 1e-4,
662  unsigned int iter_max = 10 )
663  {
664  bool ok = true;
665  if ( epsr <= 1.0 ) epsr = 2.0;
666  if ( verbose >= 1 )
667  trace.beginBlock( "#### Solve AT by Gamma-convergence ##########" );
668  if ( compute_smallest_epsilon_map ) smallest_epsilon_map.clear();
669  for ( double eps = eps1; eps >= eps2; eps /= epsr )
670  {
671  solveForEpsilon( eps, n_oo_max, iter_max );
672  if ( compute_smallest_epsilon_map )
674  }
675  if ( verbose >= 1 )
676  trace.endBlock();
677  return ok;
678  }
679 
684  void normalizeU2()
685  {
686  for ( Index index = 0; index < size( 2 ); index++)
687  {
688  double n2 = 0.0;
689  for ( unsigned int d = 0; d < u2.size(); ++d )
690  n2 += u2[ d ].myContainer( index ) * u2[ d ].myContainer( index );
691  double norm = sqrt( n2 );
692  if (norm == 0.0) continue;
693  for ( unsigned int d = 0; d < u2.size(); ++d )
694  u2[ d ].myContainer( index ) /= norm;
695  }
696  }
697 
702  std::tuple<double,double,double> diffV0() const
703  {
704  PrimalForm0 delta = v0 - former_v0;
705  delta.myContainer = delta.myContainer.cwiseAbs();
706  const double n_oo = delta.myContainer.maxCoeff();
707  const double n_2 = std::sqrt(delta.myContainer.squaredNorm()/delta.myContainer.size());
708  const double n_1 = delta.myContainer.mean();
709  return std::make_tuple( n_oo, n_2, n_1 );
710  }
711 
713 
714 
715  // ----------------------- Access services ------------------------------
716  public:
719 
723  std::tuple<double,double,double> checkV0() const
724  {
725  const double m1 = v0.myContainer.minCoeff();
726  const double m2 = v0.myContainer.maxCoeff();
727  const double ma = v0.myContainer.mean();
728  if ( verbose >= 1 )
729  trace.info() << "0-form v (should be in [0,1]): min=" << m1 << " avg=" << ma << " max=" << m2 << std::endl;
730  return std::make_tuple( m1, m2, ma );
731  }
732 
733 
736  {
737  return v0;
738  }
739 
742  {
743  return M01*v0;
744  }
745 
748  {
749  return M12*M01*v0;
750  }
751 
755  {
756  return u2[ k ];
757  }
758 
770  template <typename VectorFieldOutput, typename SurfelRangeConstIterator>
771  void
772  getOutputVectorFieldU2( VectorFieldOutput& output,
773  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
774  {
775  const Dimension N = u2.size();
776  Index i = 0;
777  for ( auto it = itB; it != itE; ++it, ++i )
778  {
779  Index idx = surfel2idx[ *it ];
780  ASSERT( output[ i ].size() >= N );
781  for ( Dimension k = 0; k < N; ++k )
782  output[ i ][ k ] = u2[ k ].myContainer( idx );
783  }
784  }
785 
796  template <typename ScalarFieldOutput, typename SurfelRangeConstIterator>
797  void
798  getOutputScalarFieldU2( ScalarFieldOutput& output,
799  SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
800  {
801  ASSERT( u2.size() == 1 && "[ATSolver2D::getOutputScalarFieldU2] "
802  "You try to output a scalar field from a vector field." );
803  Index i = 0;
804  for ( auto it = itB; it != itE; ++it, ++i )
805  {
806  Index idx = surfel2idx[ *it ];
807  output[ i ] = u2[ 0 ].myContainer( idx );
808  }
809  }
810 
822  template <typename ScalarFieldOutput, typename CellRangeConstIterator>
823  void
824  getOutputScalarFieldV0( ScalarFieldOutput& output,
825  CellRangeConstIterator itB, CellRangeConstIterator itE,
826  CellOutputPolicy policy = CellOutputPolicy::Average )
827  {
828  const KSpace& K = ptrCalculus->myKSpace;
829  const Dimension k = K.uDim( *itB );
830  ASSERT( k <= 2 );
831  Index i = 0;
832  if ( k == 0 )
833  {
834  for ( auto it = itB; it != itE; ++it, ++i )
835  {
836  const Cell pointel = *it;
837  const Index idx = ptrCalculus->getCellIndex( pointel );
838  output[ i ] = v0.myContainer( idx );
839  }
840  }
841  else if ( k == 1 )
842  {
843  for ( auto it = itB; it != itE; ++it, ++i )
844  {
845  const Cell linel = *it;
846  const Dimension d = * K.uDirs( linel );
847  const Cell p0 = K.uIncident( linel, d, false );
848  const Cell p1 = K.uIncident( linel, d, true );
849  const Index idx0 = ptrCalculus->getCellIndex( p0 );
850  const Index idx1 = ptrCalculus->getCellIndex( p1 );
851  switch (policy) {
852  case CellOutputPolicy::Average: output[ i ] = 0.5 * ( v0.myContainer( idx0 ) + v0.myContainer( idx1 ) );
853  break;
854  case CellOutputPolicy::Minimum: output[ i ] = std::min( v0.myContainer( idx0 ), v0.myContainer( idx1 ) );
855  break;
856  case CellOutputPolicy::Maximum: output[ i ] = std::max( v0.myContainer( idx0 ), v0.myContainer( idx1 ) );
857  break;
858  }
859  }
860  }
861  else if ( k == 2 )
862  {
863  for ( auto it = itB; it != itE; ++it, ++i )
864  {
865  const Cell face = *it;
866  const Dimension d = * K.uDirs( face );
867  const Cell l0 = K.uIncident( face, d, false );
868  const Cell l1 = K.uIncident( face, d, true );
869  const Dimension j = * K.uDirs( l0 );
870  const Cell p00 = K.uIncident( l0, j, false );
871  const Cell p01 = K.uIncident( l0, j, true );
872  const Cell p10 = K.uIncident( l1, j, false );
873  const Cell p11 = K.uIncident( l1, j, true );
874  const Index idx00 = ptrCalculus->getCellIndex( p00 );
875  const Index idx01 = ptrCalculus->getCellIndex( p01 );
876  const Index idx10 = ptrCalculus->getCellIndex( p10 );
877  const Index idx11 = ptrCalculus->getCellIndex( p11 );
878  switch (policy) {
879  case CellOutputPolicy::Average:
880  output[ i ] = 0.25 * ( v0.myContainer( idx00 ) + v0.myContainer( idx01 )
881  + v0.myContainer( idx10 ) + v0.myContainer( idx11 ) );
882  break;
883  case CellOutputPolicy::Minimum:
884  output[ i ] = std::min( std::min( v0.myContainer( idx00 ), v0.myContainer( idx01 ) ),
885  std::min( v0.myContainer( idx10 ), v0.myContainer( idx11 ) ) );
886  break;
887  case CellOutputPolicy::Maximum:
888  output[ i ] = std::max( std::max( v0.myContainer( idx00 ), v0.myContainer( idx01 ) ),
889  std::max( v0.myContainer( idx10 ), v0.myContainer( idx11 ) ) );
890  break;
891  }
892  }
893  }
894  }
895 
903  void updateSmallestEpsilonMap( const double threshold = .5 )
904  {
905  const KSpace& K = ptrCalculus->myKSpace;
906  for ( const SCell surfel : ptrCalculus->template getIndexedSCells<2, PRIMAL>() )
907  {
908  const Cell face = K.unsigns( surfel );
909  const Dimension k1 = * K.uDirs( face );
910  const Cell l0 = K.uIncident( face, k1, false );
911  const Cell l1 = K.uIncident( face, k1, true );
912  const Dimension k2 = * K.uDirs( l0 );
913  const Cell ll0 = K.uIncident( face, k2, false );
914  const Cell ll1 = K.uIncident( face, k2, true );
915  const Cell p00 = K.uIncident( l0, k2, false );
916  const Cell p01 = K.uIncident( l0, k2, true );
917  const Cell p10 = K.uIncident( l1, k2, false );
918  const Cell p11 = K.uIncident( l1, k2, true );
919 
920  std::vector<double> features( 4 );
921  features[ 0 ] = v0.myContainer( ptrCalculus->getCellIndex( p00 ) );
922  features[ 1 ] = v0.myContainer( ptrCalculus->getCellIndex( p01 ) );
923  features[ 2 ] = v0.myContainer( ptrCalculus->getCellIndex( p10 ) );
924  features[ 3 ] = v0.myContainer( ptrCalculus->getCellIndex( p11 ) );
925  std::sort( features.begin(), features.end() );
926 
927  if ( features[ 1 ] <= threshold )
928  {
929  auto it = smallest_epsilon_map.find( surfel );
930  if ( it != smallest_epsilon_map.end() )
931  it->second = std::min( epsilon, it->second );
932  else smallest_epsilon_map[ surfel ] = epsilon;
933  }
934  }
935  }
936 
938 
939 
940  // ----------------------- Interface --------------------------------------
941  public:
944 
949  void selfDisplay ( std::ostream & out ) const
950  {
951  auto cv = checkV0();
952  out << "[ATSolver2D] v is between min/avg/max:"
953  << std::get<0>(cv) << "/"
954  << std::get<1>(cv) << "/"
955  << std::get<2>(cv) << std::endl;
956  }
957 
962  bool isValid() const
963  {
964  return true;
965  }
966 
968 
969  // ------------------------- Hidden services ------------------------------
970  protected:
971 
974 
977  {
978  if ( verbose >= 1 ) trace.beginBlock( "[ATSolver2D::initOperators] Solver initialization" );
979  if ( verbose >= 2 ) trace.info() << "derivative of primal 0-forms: primal_D0" << std::endl;
980  primal_D0 = ptrCalculus->template derivative<0,PRIMAL>();
981  if ( verbose >= 2 ) trace.info() << "derivative of primal 1-forms: primal_D1" << std::endl;
982  primal_D1 = ptrCalculus->template derivative<1,PRIMAL>();
983  if ( verbose >= 2 ) trace.info() << "antiderivative of primal 2-forms: primal_AD2" << std::endl;
984  primal_AD2 = ptrCalculus->template antiderivative<2,PRIMAL>();
985  if ( verbose >= 2 ) trace.info() << "vertex to edge average operator: M01" << std::endl;
986  M01 = primal_D0;
987  M01.myContainer = .5 * M01.myContainer.cwiseAbs();
988  if ( verbose >= 2 ) trace.info() << "edge to face average operator: M12" << std::endl;
989  M12 = primal_D1;
990  M12.myContainer = .25 * M12.myContainer.cwiseAbs();
991  if ( verbose >= 1 ) trace.endBlock();
992  }
993 
995 
996  // ------------------------- Internals ------------------------------------
997  private:
998 
999  }; // end of class ATSolver2D
1000 
1001 
1008  template <typename T>
1009  std::ostream&
1010  operator<< ( std::ostream & out, const ATSolver2D<T> & object );
1011 
1012 } // namespace DGtal
1013 
1014 
1016 // Includes inline functions.
1017 
1018 // //
1020 
1021 #endif // !defined ATSolver2D_h
1022 
1023 #undef ATSolver2D_RECURSES
1024 #endif // else defined(ATSolver2D_RECURSES)
Aim: This class solves Ambrosio-Tortorelli functional on a two-dimensional digital space (a 2D grid o...
Definition: ATSolver2D.h:91
bool solveForEpsilon(double eps, double n_oo_max=1e-4, unsigned int iter_max=10)
Definition: ATSolver2D.h:605
Index initInputVectorFieldU2(const std::map< Surfel, ScalarVector > &input, bool normalize=false)
Definition: ATSolver2D.h:371
KSpace::Surfel Surfel
Definition: ATSolver2D.h:112
Calculus::PrimalHodge2 PrimalHodge2
Definition: ATSolver2D.h:129
Index initInputScalarFieldU2(const std::map< Surfel, Scalar > &input)
Definition: ATSolver2D.h:411
void getOutputScalarFieldU2(ScalarFieldOutput &output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:798
KSpace::Space Space
Definition: ATSolver2D.h:107
void setUp(double a, double l, const AlphaWeights &weights, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:499
void getOutputScalarFieldV0(ScalarFieldOutput &output, CellRangeConstIterator itB, CellRangeConstIterator itE, CellOutputPolicy policy=CellOutputPolicy::Average)
Definition: ATSolver2D.h:824
Calculus::PrimalDerivative0 PrimalDerivative0
Definition: ATSolver2D.h:123
static const Dimension dimension
Definition: ATSolver2D.h:99
Calculus::PrimalHodge0 PrimalHodge0
Definition: ATSolver2D.h:127
bool solveOneAlternateStep()
Definition: ATSolver2D.h:546
PrimalForm2 getU2(Dimension k) const
Definition: ATSolver2D.h:754
void initInputScalarFieldU2(const ScalarFieldInput &input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:330
std::vector< PrimalForm2 > alpha_g2
The alpha-weighted N-array of input primal 2-forms g.
Definition: ATSolver2D.h:162
std::vector< PrimalForm2 > u2
The N-array of regularized primal 2-forms u.
Definition: ATSolver2D.h:164
void setUp(double a, double l)
Definition: ATSolver2D.h:444
PrimalIdentity2 alpha_Id2
The alpha-weighted identity operator for primal 2-forms (stored for performance)
Definition: ATSolver2D.h:156
Calculus::PrimalIdentity1 PrimalIdentity1
Definition: ATSolver2D.h:121
ATSolver2D(ConstAlias< Calculus > aCalculus, int aVerbose=0)
Definition: ATSolver2D.h:202
HyperRectDomain< Space > Domain
Definition: ATSolver2D.h:113
Calculus::PrimalForm0 PrimalForm0
Definition: ATSolver2D.h:117
int verbose
Tells the verbose level.
Definition: ATSolver2D.h:191
void updateSmallestEpsilonMap(const double threshold=.5)
Definition: ATSolver2D.h:903
RealVector::Component Scalar
Definition: ATSolver2D.h:109
KSpace::SCell SCell
Definition: ATSolver2D.h:110
Calculus::PrimalAntiderivative1 PrimalAntiderivative1
Definition: ATSolver2D.h:125
EigenLinearAlgebraBackend::SolverSimplicialLDLT LinearAlgebraSolver
Definition: ATSolver2D.h:138
TLinearAlgebra LinearAlgebra
Definition: ATSolver2D.h:96
void setEpsilon(double e)
Definition: ATSolver2D.h:525
void initInputVectorFieldU2(const VectorFieldInput &input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE, bool normalize=false)
Definition: ATSolver2D.h:288
PrimalDerivative1 M12
The primal edge to face average operator.
Definition: ATSolver2D.h:152
SmallestEpsilonMap smallest_epsilon_map
Definition: ATSolver2D.h:178
Calculus::PrimalForm1 PrimalForm1
Definition: ATSolver2D.h:118
PrimalForm0 l_1_over_4e
The primal 0-form lambda/(4epsilon) (stored for performance)
Definition: ATSolver2D.h:171
PrimalForm0 getV0() const
Definition: ATSolver2D.h:735
PrimalAntiderivative2 primal_AD2
The antiderivative of primal 2-forms.
Definition: ATSolver2D.h:154
Index size(const int order) const
Definition: ATSolver2D.h:258
Calculus::PrimalIdentity2 PrimalIdentity2
Definition: ATSolver2D.h:122
PrimalForm1 getV1() const
Definition: ATSolver2D.h:741
Calculus::PrimalDerivative1 PrimalDerivative1
Definition: ATSolver2D.h:124
std::tuple< double, double, double > checkV0() const
Definition: ATSolver2D.h:723
PrimalForm0 v0
Definition: ATSolver2D.h:167
DiscreteExteriorCalculus< 2, dimension, LinearAlgebra > Calculus
Definition: ATSolver2D.h:114
void initOperators()
Initializes the operators.
Definition: ATSolver2D.h:976
void setUp(double a, double l, const std::map< Surfel, Scalar > &weights)
Definition: ATSolver2D.h:463
Calculus::PrimalAntiderivative2 PrimalAntiderivative2
Definition: ATSolver2D.h:126
ATSolver2D(ATSolver2D &&other)=default
KSpace::template SurfelMap< double >::Type SmallestEpsilonMap
Definition: ATSolver2D.h:115
Surfel2IndexMap surfel2idx
Definition: ATSolver2D.h:175
bool isValid() const
Definition: ATSolver2D.h:962
Calculus::PrimalHodge1 PrimalHodge1
Definition: ATSolver2D.h:128
ATSolver2D & operator=(const ATSolver2D &other)=default
Space::RealVector RealVector
Definition: ATSolver2D.h:108
KSpace::template SurfelMap< Index >::Type Surfel2IndexMap
Definition: ATSolver2D.h:130
Calculus::PrimalIdentity0 PrimalIdentity0
Definition: ATSolver2D.h:120
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:657
Calculus::PrimalForm2 PrimalForm2
Definition: ATSolver2D.h:119
PrimalForm2 getV2() const
Definition: ATSolver2D.h:747
PrimalDerivative1 primal_D1
the derivative operator for primal 1-forms
Definition: ATSolver2D.h:148
DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL > SolverV0
Definition: ATSolver2D.h:140
PrimalIdentity0 l_1_over_4e_Id0
The 1/(4epsilon)-weighted identity operator for primal 0-forms (stored for performance)
Definition: ATSolver2D.h:158
ATSolver2D< KSpace, LinearAlgebra > Self
Definition: ATSolver2D.h:97
void selfDisplay(std::ostream &out) const
Definition: ATSolver2D.h:949
bool normalize_u2
Indicates whether to normalize U (unit norm) at each iteration or not.
Definition: ATSolver2D.h:189
Calculus::Index Index
Definition: ATSolver2D.h:116
PrimalForm0 former_v0
The primal 0-form v at the previous iteration.
Definition: ATSolver2D.h:169
KSpace::Cell Cell
Definition: ATSolver2D.h:111
DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL > SolverU2
Definition: ATSolver2D.h:139
std::vector< PrimalForm2 > g2
The N-array of input primal 2-forms g.
Definition: ATSolver2D.h:160
~ATSolver2D()=default
std::tuple< double, double, double > diffV0() const
Definition: ATSolver2D.h:702
void getOutputVectorFieldU2(VectorFieldOutput &output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)
Definition: ATSolver2D.h:772
PrimalDerivative0 primal_D0
the derivative operator for primal 0-forms
Definition: ATSolver2D.h:146
@ Maximum
compute maximum value at cell vertices
Definition: ATSolver2D.h:105
@ Average
compute average values at cell vertices
Definition: ATSolver2D.h:103
@ Minimum
compute minimum value at cell vertices,
Definition: ATSolver2D.h:104
CountedConstPtrOrConstPtr< Calculus > ptrCalculus
A smart (or not) pointer to a calculus object.
Definition: ATSolver2D.h:144
PrimalDerivative0 M01
The primal vertex to edge average operator.
Definition: ATSolver2D.h:150
ATSolver2D(const ATSolver2D &other)=default
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
Aim: Smart or simple const pointer on T. It can be a smart pointer based on reference counts or a sim...
SolutionKForm solve(const InputKForm &input_kform) const
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
LinearAlgebraBackend::DenseVector::Index Index
Aim: KForm represents discrete kforms in the dec package.
Definition: KForm.h:66
static KForm< TCalculus, order, duality > ones(ConstAlias< Calculus > calculus)
Container myContainer
Definition: KForm.h:131
static KForm< TCalculus, order, duality > zeros(ConstAlias< Calculus > calculus)
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
TransposedLinearOperator transpose() const
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
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
SMesh::Index Index
DGtal::LinearOperator< Calculus, dim, duality, dim, duality > diagonal(const DGtal::KForm< Calculus, dim, duality > &kform)
Definition: DECHelpers.h:60
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
@ PRIMAL
Definition: Duality.h:61
Eigen::SimplicialLDLT< SparseMatrix > SolverSimplicialLDLT
Definition: EigenSupport.h:106
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
int max(int a, int b)
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition: testBits.cpp:44
KSpace K