DGtal  1.3.beta
DGtal::ATSolver2D< TKSpace, TLinearAlgebra > Class Template Reference

Aim: This class solves Ambrosio-Tortorelli functional on a two-dimensional digital space (a 2D grid or 2D digital surface) for a piecewise smooth scalar/vector function u represented as one/several 2-form(s) and a discontinuity function v represented as a 0-form. The 2-form(s) u is a regularized approximation of an input vector data g, while v represents the set of discontinuities of u. The norm chosen for u is the $$l_2$$-norm. More...

#include <DGtal/dec/ATSolver2D.h>

## Public Types

enum  CellOutputPolicy { Average, Minimum, Maximum }

typedef TKSpace KSpace

typedef TLinearAlgebra LinearAlgebra

typedef ATSolver2D< KSpace, LinearAlgebraSelf

typedef KSpace::Space Space

typedef Space::RealVector RealVector

typedef RealVector::Component Scalar

typedef KSpace::SCell SCell

typedef KSpace::Cell Cell

typedef KSpace::Surfel Surfel

typedef HyperRectDomain< SpaceDomain

typedef DiscreteExteriorCalculus< 2, dimension, LinearAlgebraCalculus

typedef KSpace::template SurfelMap< double >::Type SmallestEpsilonMap

typedef Calculus::Index Index

typedef Calculus::PrimalForm0 PrimalForm0

typedef Calculus::PrimalForm1 PrimalForm1

typedef Calculus::PrimalForm2 PrimalForm2

typedef Calculus::PrimalIdentity0 PrimalIdentity0

typedef Calculus::PrimalIdentity1 PrimalIdentity1

typedef Calculus::PrimalIdentity2 PrimalIdentity2

typedef Calculus::PrimalDerivative0 PrimalDerivative0

typedef Calculus::PrimalDerivative1 PrimalDerivative1

typedef Calculus::PrimalAntiderivative1 PrimalAntiderivative1

typedef Calculus::PrimalAntiderivative2 PrimalAntiderivative2

typedef Calculus::PrimalHodge0 PrimalHodge0

typedef Calculus::PrimalHodge1 PrimalHodge1

typedef Calculus::PrimalHodge2 PrimalHodge2

typedef KSpace::template SurfelMap< Index >::Type Surfel2IndexMap

typedef EigenLinearAlgebraBackend::SolverSimplicialLDLT LinearAlgebraSolver

typedef DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMALSolverU2

typedef DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMALSolverV0

## Public Member Functions

Standard services
ATSolver2D (ConstAlias< Calculus > aCalculus, int aVerbose=0)

ATSolver2D ()=delete

~ATSolver2D ()=default

ATSolver2D (const ATSolver2D &other)=default

ATSolver2D (ATSolver2D &&other)=default

ATSolver2Doperator= (const ATSolver2D &other)=default

ATSolver2Doperator= (ATSolver2D &&other)=default

Index size (const int order) const

Initialization services
template<typename VectorFieldInput , typename SurfelRangeConstIterator >
void initInputVectorFieldU2 (const VectorFieldInput &input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE, bool normalize=false)

template<typename ScalarFieldInput , typename SurfelRangeConstIterator >
void initInputScalarFieldU2 (const ScalarFieldInput &input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)

template<typename ScalarVector >
Index initInputVectorFieldU2 (const std::map< Surfel, ScalarVector > &input, bool normalize=false)

template<typename Scalar >
Index initInputScalarFieldU2 (const std::map< Surfel, Scalar > &input)

void setUp (double a, double l)

void setUp (double a, double l, const std::map< Surfel, Scalar > &weights)

template<typename AlphaWeights , typename SurfelRangeConstIterator >
void setUp (double a, double l, const AlphaWeights &weights, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)

void setEpsilon (double e)

Optimization services
bool solveOneAlternateStep ()

bool solveForEpsilon (double eps, double n_oo_max=1e-4, unsigned int iter_max=10)

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)

void normalizeU2 ()

std::tuple< double, double, double > diffV0 () const

Access services
std::tuple< double, double, double > checkV0 () const

PrimalForm0 getV0 () const

PrimalForm1 getV1 () const

PrimalForm2 getV2 () const

PrimalForm2 getU2 (Dimension k) const

template<typename VectorFieldOutput , typename SurfelRangeConstIterator >
void getOutputVectorFieldU2 (VectorFieldOutput &output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)

template<typename ScalarFieldOutput , typename SurfelRangeConstIterator >
void getOutputScalarFieldU2 (ScalarFieldOutput &output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE)

template<typename ScalarFieldOutput , typename CellRangeConstIterator >
void getOutputScalarFieldV0 (ScalarFieldOutput &output, CellRangeConstIterator itB, CellRangeConstIterator itE, CellOutputPolicy policy=CellOutputPolicy::Average)

Interface services
void selfDisplay (std::ostream &out) const

bool isValid () const

## Data Fields

Surfel2IndexMap surfel2idx

SmallestEpsilonMap smallest_epsilon_map

double alpha

double lambda

double epsilon

bool normalize_u2
Indicates whether to normalize U (unit norm) at each iteration or not. More...

int verbose
Tells the verbose level. More...

## Static Public Attributes

static const Dimension dimension = KSpace::dimension

## Protected Member Functions

Hidden services
void initOperators ()
Initializes the operators. More...

## Protected Attributes

CountedConstPtrOrConstPtr< CalculusptrCalculus
A smart (or not) pointer to a calculus object. More...

PrimalDerivative0 primal_D0
the derivative operator for primal 0-forms More...

PrimalDerivative1 primal_D1
the derivative operator for primal 1-forms More...

PrimalDerivative0 M01
The primal vertex to edge average operator. More...

PrimalDerivative1 M12
The primal edge to face average operator. More...

The antiderivative of primal 2-forms. More...

PrimalIdentity2 alpha_Id2
The alpha-weighted identity operator for primal 2-forms (stored for performance) More...

PrimalIdentity0 l_1_over_4e_Id0
The 1/(4epsilon)-weighted identity operator for primal 0-forms (stored for performance) More...

std::vector< PrimalForm2g2
The N-array of input primal 2-forms g. More...

std::vector< PrimalForm2alpha_g2
The alpha-weighted N-array of input primal 2-forms g. More...

std::vector< PrimalForm2u2
The N-array of regularized primal 2-forms u. More...

PrimalForm0 v0

PrimalForm0 former_v0
The primal 0-form v at the previous iteration. More...

PrimalForm0 l_1_over_4e
The primal 0-form lambda/(4epsilon) (stored for performance) More...

## Detailed Description

### template<typename TKSpace, typename TLinearAlgebra = EigenLinearAlgebraBackend> class DGtal::ATSolver2D< TKSpace, TLinearAlgebra >

Aim: This class solves Ambrosio-Tortorelli functional on a two-dimensional digital space (a 2D grid or 2D digital surface) for a piecewise smooth scalar/vector function u represented as one/several 2-form(s) and a discontinuity function v represented as a 0-form. The 2-form(s) u is a regularized approximation of an input vector data g, while v represents the set of discontinuities of u. The norm chosen for u is the $$l_2$$-norm.

Description of template class 'ATSolver2D'

Template Parameters
 TKSpace any model of CCellularGridSpaceND, e.g KhalimskySpaceND TLinearAlgebra any back-end for performing linear algebra, default is EigenLinearAlgebraBackend.
// Typical use (with appropriate definitions for types and variables).
typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;
const auto calculus = CalculusFactory::createFromNSCells<2>( surfels.begin(), surfels.end() );
ATSolver2D< KSpace > at_solver(calculus, 1);
at_solver.initInputVectorFieldU2( normals, surfels.cbegin(), surfels.cend() );
at_solver.setUp( alpha_at, lambda_at );
at_solver.solveGammaConvergence( 2.0, 0.5, 2.0 );
at_solver.getOutputVectorFieldU2( normals, surfels.cbegin(), surfels.cend() );
exampleSurfaceATNormals.cpp

Definition at line 90 of file ATSolver2D.h.

## ◆ Calculus

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef DiscreteExteriorCalculus<2,dimension, LinearAlgebra> DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Calculus

Definition at line 114 of file ATSolver2D.h.

## ◆ Cell

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef KSpace::Cell DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Cell

Definition at line 111 of file ATSolver2D.h.

## ◆ Domain

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef HyperRectDomain DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Domain

Definition at line 113 of file ATSolver2D.h.

## ◆ Index

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::Index DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Index

Definition at line 116 of file ATSolver2D.h.

## ◆ KSpace

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef TKSpace DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::KSpace

Definition at line 95 of file ATSolver2D.h.

## ◆ LinearAlgebra

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef TLinearAlgebra DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::LinearAlgebra

Definition at line 96 of file ATSolver2D.h.

## ◆ LinearAlgebraSolver

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef EigenLinearAlgebraBackend::SolverSimplicialLDLT DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::LinearAlgebraSolver

Definition at line 138 of file ATSolver2D.h.

## ◆ PrimalAntiderivative1

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalAntiderivative1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalAntiderivative1

Definition at line 125 of file ATSolver2D.h.

## ◆ PrimalAntiderivative2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalAntiderivative2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalAntiderivative2

Definition at line 126 of file ATSolver2D.h.

## ◆ PrimalDerivative0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalDerivative0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalDerivative0

Definition at line 123 of file ATSolver2D.h.

## ◆ PrimalDerivative1

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalDerivative1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalDerivative1

Definition at line 124 of file ATSolver2D.h.

## ◆ PrimalForm0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalForm0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalForm0

Definition at line 117 of file ATSolver2D.h.

## ◆ PrimalForm1

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalForm1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalForm1

Definition at line 118 of file ATSolver2D.h.

## ◆ PrimalForm2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalForm2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalForm2

Definition at line 119 of file ATSolver2D.h.

## ◆ PrimalHodge0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalHodge0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalHodge0

Definition at line 127 of file ATSolver2D.h.

## ◆ PrimalHodge1

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalHodge1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalHodge1

Definition at line 128 of file ATSolver2D.h.

## ◆ PrimalHodge2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalHodge2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalHodge2

Definition at line 129 of file ATSolver2D.h.

## ◆ PrimalIdentity0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalIdentity0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalIdentity0

Definition at line 120 of file ATSolver2D.h.

## ◆ PrimalIdentity1

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalIdentity1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalIdentity1

Definition at line 121 of file ATSolver2D.h.

## ◆ PrimalIdentity2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Calculus::PrimalIdentity2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::PrimalIdentity2

Definition at line 122 of file ATSolver2D.h.

## ◆ RealVector

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef Space::RealVector DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::RealVector

Definition at line 108 of file ATSolver2D.h.

## ◆ Scalar

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef RealVector::Component DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Scalar

Definition at line 109 of file ATSolver2D.h.

## ◆ SCell

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef KSpace::SCell DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::SCell

Definition at line 110 of file ATSolver2D.h.

## ◆ Self

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef ATSolver2D< KSpace, LinearAlgebra > DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Self

Definition at line 97 of file ATSolver2D.h.

## ◆ SmallestEpsilonMap

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef KSpace::template SurfelMap::Type DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::SmallestEpsilonMap

Definition at line 115 of file ATSolver2D.h.

## ◆ SolverU2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef DiscreteExteriorCalculusSolver DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::SolverU2

Definition at line 139 of file ATSolver2D.h.

## ◆ SolverV0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef DiscreteExteriorCalculusSolver DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::SolverV0

Definition at line 140 of file ATSolver2D.h.

## ◆ Space

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef KSpace::Space DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Space

Definition at line 107 of file ATSolver2D.h.

## ◆ Surfel

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef KSpace::Surfel DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Surfel

Definition at line 112 of file ATSolver2D.h.

## ◆ Surfel2IndexMap

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 typedef KSpace::template SurfelMap::Type DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::Surfel2IndexMap

Definition at line 130 of file ATSolver2D.h.

## ◆ CellOutputPolicy

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>

Specifies how to merge the different values of 0-form v at cell vertices when outputing the 0-form v for a range of cells (either pointels, linels, surfels).

Enumerator
Average

compute average values at cell vertices

Minimum

compute minimum value at cell vertices,

Maximum

compute maximum value at cell vertices

Definition at line 103 of file ATSolver2D.h.

103  { Average,
104  Minimum,
105  Maximum,
106  };

## ◆ ATSolver2D() [1/4]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::ATSolver2D ( ConstAlias< Calculus > aCalculus, int aVerbose = 0 )
inline

Prepare an AT-solver from a valid calculus.

Parameters
 aCalculus any valid calculus aVerbose tells how the solver displays computing information: 0 none, 1 more, 2 even more...
DiscreteExteriorCalculusFactory for creating calculus objects.

Definition at line 202 of file ATSolver2D.h.

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  }

## ◆ ATSolver2D() [2/4]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::ATSolver2D ( )
delete

Default constructor.

## ◆ ~ATSolver2D()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::~ATSolver2D ( )
default

Destructor.

## ◆ ATSolver2D() [3/4]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::ATSolver2D ( const ATSolver2D< TKSpace, TLinearAlgebra > & other )
default

Copy constructor.

Parameters
 other the object to clone.

## ◆ ATSolver2D() [4/4]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::ATSolver2D ( ATSolver2D< TKSpace, TLinearAlgebra > && other )
default

Move constructor.

Parameters
 other the object to move.

## ◆ checkV0()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 std::tuple DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::checkV0 ( ) const
inline

Debug method for checking if v is a scalar field between 0 and 1.

Returns
the tuple (min(v), average(v), max(v))

Definition at line 722 of file ATSolver2D.h.

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  }

## ◆ diffV0()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 std::tuple DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::diffV0 ( ) const
inline

Computes the norms loo, l2, l1 of (v - former_v), i.e. the evolution of discontinuity function v.

Returns
a tuple (n_infty,n_2,n_1) giving the loo/l2/l1-norm of (v - former_v)

Definition at line 701 of file ATSolver2D.h.

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  }

## ◆ getOutputScalarFieldU2()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename ScalarFieldOutput , typename SurfelRangeConstIterator >
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getOutputScalarFieldU2 ( ScalarFieldOutput & output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
inline

Given a range of surfels [itB,itE), returns in output the regularized scalar field u.

Template Parameters
 ScalarFieldOutput the type of scalar field for output values (RandomAccess container) SurfelRangeConstIterator the type of iterator for traversing a range of surfels
Parameters
 [out] output the vector of output values (a scalar field), which should be of size length(itB,itE) itB the start of the range of surfels. itE past the end of the range of surfels.

Definition at line 797 of file ATSolver2D.h.

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  }

## ◆ getOutputScalarFieldV0()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename ScalarFieldOutput , typename CellRangeConstIterator >
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getOutputScalarFieldV0 ( ScalarFieldOutput & output, CellRangeConstIterator itB, CellRangeConstIterator itE, CellOutputPolicy policy = CellOutputPolicy::Average )
inline

Given a range of pointels, linels or 2-cells [itB,itE), returns in output the feature vector v (the average of v for linels/surfels).

Template Parameters
 ScalarFieldOutput the type of scalar field for output values (RandomAccess container) CellRangeConstIterator the type of iterator for traversing a range of cells
Parameters
 [out] output the vector of output scalar values (a scalar field), which should be of size length(itB,itE) [in] itB the start of the range of cells. [in] itE past the end of the range of cells. [in] policy the chosen policy for outputing v values for a given cell.

Definition at line 823 of file ATSolver2D.h.

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  }

## ◆ getOutputVectorFieldU2()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename VectorFieldOutput , typename SurfelRangeConstIterator >
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getOutputVectorFieldU2 ( VectorFieldOutput & output, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
inline

Given a range of surfels [itB,itE), returns in output the regularized vector field u.

Template Parameters
 VectorFieldOutput the type of vector field for output values (RandomAccess container) SurfelRangeConstIterator the type of iterator for traversing a range of surfels
Parameters
 [out] output the vector of output values (a scalar or vector field depending on input). itB the start of the range of surfels. itE past the end of the range of surfels.

Definition at line 771 of file ATSolver2D.h.

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  }

Referenced by DGtal::ShortcutsGeometry< TKSpace >::getATVectorFieldApproximation(), and main().

## ◆ getU2()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getU2 ( Dimension k ) const
inline
Parameters
 k an integer such that 0 <= k < u2.size()
Returns
the k-th piecewise smooth function u as a primal 2-form.

Definition at line 753 of file ATSolver2D.h.

754  {
755  return u2[ k ];
756  }

## ◆ getV0()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getV0 ( ) const
inline
Returns
the discontinuity function v as a primal 0-form.

Definition at line 734 of file ATSolver2D.h.

735  {
736  return v0;
737  }

## ◆ getV1()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getV1 ( ) const
inline
Returns
the discontinuity function v as a primal 1-form.

Definition at line 740 of file ATSolver2D.h.

741  {
742  return M01*v0;
743  }

## ◆ getV2()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::getV2 ( ) const
inline
Returns
the discontinuity function u as a primal 2-form.

Definition at line 746 of file ATSolver2D.h.

747  {
748  return M12*M01*v0;
749  }

## ◆ initInputScalarFieldU2() [1/2]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename ScalarFieldInput , typename SurfelRangeConstIterator >
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::initInputScalarFieldU2 ( const ScalarFieldInput & input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
inline

Given a range of surfels [itB,itE) and an input scalar field, initializes the AT 2-forms u and g. The 0-form v is itself initialized to 1 everywhere.

Template Parameters
 ScalarFieldInput the type of scalar field for input values (RandomAccess container) SurfelRangeConstIterator the type of iterator for traversing a range of surfels
Parameters
 [in] input the input scalar field (a vector of scalar values) itB the start of the range of surfels. itE past the end of the range of surfels.

Definition at line 330 of file ATSolver2D.h.

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  }

## ◆ initInputScalarFieldU2() [2/2]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename Scalar >
 Index DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::initInputScalarFieldU2 ( const std::map< Surfel, Scalar > & input )
inline

Given a map Surfel -> Scalar, initializes forms g, u and v of the AT solver. Note that there is only one 2-form u/g.

Parameters
 input any map Surfel -> Scalar
Returns
the number of cells that were initialized.
Template Parameters
 Scalar any type representing a scalar (float, double)

Definition at line 411 of file ATSolver2D.h.

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  }

## ◆ initInputVectorFieldU2() [1/2]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename ScalarVector >
 Index DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::initInputVectorFieldU2 ( const std::map< Surfel, ScalarVector > & input, bool normalize = false )
inline

Given a map Surfel -> ScalarVector, initializes forms g, u and v of the AT solver. Note that there are as many 2-forms u/g as the number of dimensions of ScalarVector.

Parameters
 input any map Surfel -> ScalarVector normalize when 'true', the input is supposed to be a unit vector field and the solver will output a unit regularized vector field at the end of each minimization step.
Returns
the number of cells that were initialized.
Template Parameters
 ScalarVector any type representing a vector/array of scalars (float, double)

Definition at line 371 of file ATSolver2D.h.

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  }

## ◆ initInputVectorFieldU2() [2/2]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename VectorFieldInput , typename SurfelRangeConstIterator >
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::initInputVectorFieldU2 ( const VectorFieldInput & input, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE, bool normalize = false )
inline

Given a range of surfels [itB,itE) and an input vector field, initializes the AT 2-forms u and g. The 0-form v is itself initialized to 1 everywhere.

Template Parameters
 VectorFieldInput the type of vector field for input values (RandomAccess container) SurfelRangeConstIterator the type of iterator for traversing a range of surfels
Parameters
 [in] input the input vector field (a vector of vector values) [in] itB the start of the range of surfels. [in] itE past the end of the range of surfels. [in] normalize when 'true', the input is supposed to be a unit vector field and the solver will output a unit regularized vector field at the end of each minimization step.

Definition at line 288 of file ATSolver2D.h.

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  }

Referenced by DGtal::ShortcutsGeometry< TKSpace >::getATVectorFieldApproximation(), and main().

## ◆ initOperators()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::initOperators ( )
inlineprotected

Initializes the operators.

Definition at line 975 of file ATSolver2D.h.

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;
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  }

## ◆ isValid()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 bool DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::isValid ( ) const
inline

Checks the validity/consistency of the object.

Returns
'true' if the object is valid, 'false' otherwise.

Definition at line 961 of file ATSolver2D.h.

962  {
963  return true;
964  }

## ◆ normalizeU2()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::normalizeU2 ( )
inline

Forces the normalization of the vector u, meaning for all index i, $$\sum_{k=0}^{K-1} u[k][i]^2 = 1$$. Can be useful in some applications where you are looking for unitary vector field.

Definition at line 683 of file ATSolver2D.h.

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  }

## ◆ operator=() [1/2]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 ATSolver2D& DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::operator= ( ATSolver2D< TKSpace, TLinearAlgebra > && other )
default

Move assignment operator.

Parameters
 other the object to move.
Returns
a reference on 'this'.

## ◆ operator=() [2/2]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 ATSolver2D& DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::operator= ( const ATSolver2D< TKSpace, TLinearAlgebra > & other )
default

Copy assignment operator.

Parameters
 other the object to copy.
Returns
a reference on 'this'.

## ◆ selfDisplay()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::selfDisplay ( std::ostream & out ) const
inline

Writes/Displays the object on an output stream.

Parameters
 out the output stream where the object is written.

Definition at line 948 of file ATSolver2D.h.

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  }

## ◆ setEpsilon()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::setEpsilon ( double e )
inline

Initializes the epsilon parameter of AT and precomputes the assaociated forms and operators.

Parameters
 e the epsilon parameter in AT

Definition at line 525 of file ATSolver2D.h.

526  {
527  epsilon = e;
529  l_1_over_4e_Id0 = (lambda/4./epsilon)*ptrCalculus->template identity<0, PRIMAL>();
530  }

## ◆ setUp() [1/3]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::setUp ( double a, double l )
inline

Initializes the alpha and lambda parameters of AT.

Parameters
 a the global alpha parameter l the global lambda parameter
Note
all 2-cells have the same weight for the data term.

Definition at line 444 of file ATSolver2D.h.

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  }

## ◆ setUp() [2/3]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
template<typename AlphaWeights , typename SurfelRangeConstIterator >
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::setUp ( double a, double l, const AlphaWeights & weights, SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
inline

Initializes the alpha and lambda parameters of AT, with weights on the 2-cells for the data terms.

Template Parameters
 AlphaWeights the type of RandomAccess container for alpha weight values SurfelRangeConstIterator the type of iterator for traversing a range of surfels
Parameters
 [in] a the global alpha parameter [in] l the global lambda parameter [in] weights the vector of alpha weights for each surfel of the range [itB,itE) [in] itB the start of the range of surfels. [in] itE past the end of the range of surfels.
Note
Useful for inpainting applications or for adaptive piecewise smooth reconstruction.

Definition at line 499 of file ATSolver2D.h.

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  }

## ◆ setUp() [3/3]

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::setUp ( double a, double l, const std::map< Surfel, Scalar > & weights )
inline

Initializes the alpha and lambda parameters of AT, with weights on the 2-cells for the data terms.

Parameters
 a the global alpha parameter l the global lambda parameter weights the map Surfel -> Scalar that gives the weight of each 2-cell in the data terms.
Note
Useful for inpainting applications or for adaptive piecewise smooth reconstruction.

Definition at line 463 of file ATSolver2D.h.

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  }

## ◆ size()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 Index DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::size ( const int order ) const
inline
Parameters
 order the dimension of cells (0,1,2)
Returns
the number of cells with dimension order

Definition at line 258 of file ATSolver2D.h.

259  {
260  return ptrCalculus->kFormLength(order, PRIMAL);
261  }

## ◆ solveForEpsilon()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 bool DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::solveForEpsilon ( double eps, double n_oo_max = 1e-4, unsigned int iter_max = 10 )
inline

Solves the alternate minimization of AT for a given eps. Solves for u then for v till convergence.

Parameters
 eps the epsilon parameter at which AT is solved. n_oo_max the alternate minimization will stop when the loo-norm of $$v^{k+1} - v^k$$ is below this bound. iter_max the alternate minimization will stop when the number of minimization steps exceeds iter_max.
Returns
true if everything went fine, false if there was a problem in the optimization.
Note
Use diffV0 to check if you are close to a critical point of AT.

Definition at line 605 of file ATSolver2D.h.

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;
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  }

## ◆ solveGammaConvergence()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 bool DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::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 )
inline

Solves AT by progressively decreasing epsilon from eps1 to eps2. AT is solved with solveForEpsilon at each epsilon.

Parameters
 eps1 the first epsilon parameter at which AT is solved. eps2 the last epsilon parameter at which AT is solved. epsr the ratio (>1) used to decrease progressively epsilon. compute_smallest_epsilon_map when 'true' determines for each surfel the smallest epsilon for which it is a discontinuity. n_oo_max the alternate minimization will stop when the loo-norm of $$v^{k+1} - v^k$$ is below this bound. iter_max the alternate minimization will stop when the number of minimization steps exceeds iter_max.
Returns
true if everything went fine, false if there was a problem in the optimization.

Definition at line 656 of file ATSolver2D.h.

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 )
673  }
674  if ( verbose >= 1 )
675  trace.endBlock();
676  return ok;
677  }

## ◆ solveOneAlternateStep()

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 bool DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::solveOneAlternateStep ( )
inline

Solves one step of the alternate minimization of AT. Solves for u then for v.

Returns
true if everything went fine, false if there was a problem in the optimization.
Note
Use diffV0 to check if you are close to a critical point of AT.

Definition at line 546 of file ATSolver2D.h.

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
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  }

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 void DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::updateSmallestEpsilonMap ( const double threshold = .5 )
inline

Computes the map that stores for each surfel the smallest epsilon for which the surfel was in the discontinuity zone (more precisely, the surfel has at least two vertices that belongs to the set of discontinuity).

Parameters
 [in] threshold the threshold for discontinuity function v (below u is discontinuous, above u is continuous)

Definition at line 902 of file ATSolver2D.h.

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  }

## ◆ alpha

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 double DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::alpha

The global coefficient alpha giving the smoothness of the reconstruction (the smaller, the smoother)

Definition at line 181 of file ATSolver2D.h.

## ◆ alpha_g2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 std::vector DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::alpha_g2
protected

The alpha-weighted N-array of input primal 2-forms g.

Definition at line 162 of file ATSolver2D.h.

## ◆ alpha_Id2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalIdentity2 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::alpha_Id2
protected

The alpha-weighted identity operator for primal 2-forms (stored for performance)

Definition at line 156 of file ATSolver2D.h.

## ◆ dimension

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 const Dimension DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::dimension = KSpace::dimension
static

Definition at line 99 of file ATSolver2D.h.

## ◆ epsilon

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 double DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::epsilon

The global coefficient epsilon giving the width of the discontinuities (the smaller, the thinner)

Definition at line 187 of file ATSolver2D.h.

## ◆ former_v0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::former_v0
protected

The primal 0-form v at the previous iteration.

Definition at line 169 of file ATSolver2D.h.

## ◆ g2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 std::vector DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::g2
protected

The N-array of input primal 2-forms g.

Definition at line 160 of file ATSolver2D.h.

## ◆ l_1_over_4e

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::l_1_over_4e
protected

The primal 0-form lambda/(4epsilon) (stored for performance)

Definition at line 171 of file ATSolver2D.h.

## ◆ l_1_over_4e_Id0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalIdentity0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::l_1_over_4e_Id0
protected

The 1/(4epsilon)-weighted identity operator for primal 0-forms (stored for performance)

Definition at line 158 of file ATSolver2D.h.

## ◆ lambda

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 double DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::lambda

The global coefficient lambda giving the length of discontinuities (the smaller, the more discontinuities)

Definition at line 184 of file ATSolver2D.h.

## ◆ M01

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalDerivative0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::M01
protected

The primal vertex to edge average operator.

Definition at line 150 of file ATSolver2D.h.

## ◆ M12

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalDerivative1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::M12
protected

The primal edge to face average operator.

Definition at line 152 of file ATSolver2D.h.

## ◆ normalize_u2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 bool DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::normalize_u2

Indicates whether to normalize U (unit norm) at each iteration or not.

Definition at line 189 of file ATSolver2D.h.

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
protected

The antiderivative of primal 2-forms.

Definition at line 154 of file ATSolver2D.h.

## ◆ primal_D0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalDerivative0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::primal_D0
protected

the derivative operator for primal 0-forms

Definition at line 146 of file ATSolver2D.h.

## ◆ primal_D1

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalDerivative1 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::primal_D1
protected

the derivative operator for primal 1-forms

Definition at line 148 of file ATSolver2D.h.

## ◆ ptrCalculus

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 CountedConstPtrOrConstPtr< Calculus > DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::ptrCalculus
protected

## ◆ smallest_epsilon_map

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 SmallestEpsilonMap DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::smallest_epsilon_map

The map Surfel -> double telling the smallest epsilon for which the surfel was a discontinuity.

Definition at line 178 of file ATSolver2D.h.

## ◆ surfel2idx

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 Surfel2IndexMap DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::surfel2idx

## ◆ u2

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 std::vector DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::u2
protected

The N-array of regularized primal 2-forms u.

Definition at line 164 of file ATSolver2D.h.

## ◆ v0

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 PrimalForm0 DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::v0
protected

The primal 0-form v representing the set of discontinuities S (more precisely $$1 - \chi_S$$)

Definition at line 167 of file ATSolver2D.h.

## ◆ verbose

template<typename TKSpace , typename TLinearAlgebra = EigenLinearAlgebraBackend>
 int DGtal::ATSolver2D< TKSpace, TLinearAlgebra >::verbose

The documentation for this class was generated from the following file:
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::alpha_Id2
PrimalIdentity2 alpha_Id2
The alpha-weighted identity operator for primal 2-forms (stored for performance)
Definition: ATSolver2D.h:156
DGtal::Trace::endBlock
double endBlock()
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::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::alpha_g2
std::vector< PrimalForm2 > alpha_g2
The alpha-weighted N-array of input primal 2-forms g.
Definition: ATSolver2D.h:162
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
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::DiscreteExteriorCalculusSolver::compute
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::ATSolver2D::ptrCalculus
CountedConstPtrOrConstPtr< Calculus > ptrCalculus
A smart (or not) pointer to a calculus object.
Definition: ATSolver2D.h:144
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::ATSolver2D::PrimalIdentity2
Calculus::PrimalIdentity2 PrimalIdentity2
Definition: ATSolver2D.h:122
DGtal::ATSolver2D::PrimalIdentity0
Calculus::PrimalIdentity0 PrimalIdentity0
Definition: ATSolver2D.h:120
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::primal_D0
PrimalDerivative0 primal_D0
the derivative operator for primal 0-forms
Definition: ATSolver2D.h:146
DGtal::PRIMAL
@ PRIMAL
Definition: Duality.h:61
DGtal::ATSolver2D::initOperators
void initOperators()
Initializes the operators.
Definition: ATSolver2D.h:975
DGtal::Trace::info
std::ostream & info()
DGtal::LinearOperator::transpose
TransposedLinearOperator transpose() const
DGtal::ATSolver2D::solveOneAlternateStep
bool solveOneAlternateStep()
Definition: ATSolver2D.h:546
DGtal::ATSolver2D::size
Index size(const int order) const
Definition: ATSolver2D.h:258
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::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
Definition: ATSolver2D.h:902
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::SolverU2
DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL > SolverU2
Definition: ATSolver2D.h:139
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::diffV0
std::tuple< double, double, double > diffV0() const
Definition: ATSolver2D.h:701
DGtal::ATSolver2D::SolverV0
DiscreteExteriorCalculusSolver< Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL > SolverV0
Definition: ATSolver2D.h:140
DGtal::KForm::zeros
static KForm< TCalculus, order, duality > zeros(ConstAlias< Calculus > calculus)
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::solveForEpsilon
bool solveForEpsilon(double eps, double n_oo_max=1e-4, unsigned int iter_max=10)
Definition: ATSolver2D.h:605
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::lambda
double lambda
Definition: ATSolver2D.h:184
DGtal::ATSolver2D::primal_D1
PrimalDerivative1 primal_D1
the derivative operator for primal 1-forms
Definition: ATSolver2D.h:148
DGtal::functions::dec::diagonal
DGtal::LinearOperator< Calculus, dim, duality, dim, duality > diagonal(const DGtal::KForm< Calculus, dim, duality > &kform)
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::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64
DGtal::LinearOperator::myContainer
Container myContainer
Definition: LinearOperator.h:117