31 #if defined(ATSolver2D_RECURSES)
32 #error Recursive header files inclusion detected in ATSolver2D.h
33 #else // defined(ATSolver2D_RECURSES)
35 #define ATSolver2D_RECURSES
37 #if !defined ATSolver2D_h
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"
88 template <
typename TKSpace,
89 typename TLinearAlgebra = EigenLinearAlgebraBackend >
160 std::vector<PrimalForm2>
g2;
164 std::vector<PrimalForm2>
u2;
214 for (
Index index = 0; index < size2; ++
index) {
286 template <
typename VectorFieldInput,
typename SurfelRangeConstIterator>
289 SurfelRangeConstIterator itB, SurfelRangeConstIterator itE,
290 bool normalize =
false )
293 trace.
beginBlock(
"[ATSolver2D::initInputVectorFieldU2] Initializing input data" );
294 ASSERT( ! input.empty() );
296 if (
verbose >= 2 )
trace.
info() <<
"input g as " << N <<
" 2-forms." << std::endl;
300 for (
auto it = itB; it != itE; ++it, ++i )
304 g2[ k ].myContainer( idx ) = input[ i ][ k ];
308 trace.
info() <<
"Unknown u[:] = g[:] at beginning." << std::endl;
328 template <
typename ScalarFieldInput,
typename SurfelRangeConstIterator>
331 SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
334 trace.
beginBlock(
"[ATSolver2D::initInputVectorFieldU2] Initializing input data" );
335 ASSERT( ! input.empty() );
340 for (
auto it = itB; it != itE; ++it, ++i )
343 g2[ 0 ].myContainer( idx ) = input[ i ];
347 trace.
info() <<
"Unknown u[:] = g[:] at beginning." << std::endl;
369 template <
typename ScalarVector>
372 bool normalize =
false )
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;
380 const ScalarVector zero;
382 for (
Index index = 0; index <
size(2); index++)
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;
389 g2[ k ].myContainer( index ) = n[ k ];
393 trace.
info() <<
"Unknown u[:] = g[:] at beginning." << std::endl;
409 template <
typename Scalar>
414 if (
verbose >= 2 )
trace.
info() <<
"discontinuity 0-form v = 1." << std::endl;
421 for (
Index index = 0; index <
size(2); index++)
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;
430 if (
verbose >= 2 )
trace.
info() <<
"Unknown u[:] = g[:] at beginning." << std::endl;
463 void setUp(
double a,
double l,
const std::map<Surfel,Scalar>& weights )
470 trace.
info() <<
"Using variable weights for fitting (alpha term)" << std::endl;
473 for (
Index index = 0; index <
size( 2 ); index++)
475 const SCell& cell =
g2[ 0 ].getSCell( index );
476 const Scalar& w = weights[ cell ];
479 alpha_g2[ k ].myContainer( index ) *= w;
498 template <
typename AlphaWeights,
typename SurfelRangeConstIterator>
500 const AlphaWeights& weights,
501 SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
508 trace.
info() <<
"Using variable weights for fitting (alpha term)" << std::endl;
512 for (
auto it = itB; it != itE; ++it, ++i )
515 const Scalar w = weights[ i ];
518 alpha_g2[ k ].myContainer( idx ) *= w;
548 bool solve_ok =
true;
555 if (
verbose >= 2 )
trace.
info() <<
"Prefactoring matrix U associated to u" << std::endl;
558 for (
Dimension d = 0; d < u2.size(); ++d )
560 if (
verbose >= 2 )
trace.
info() <<
"Solving U u[" << d <<
"] = a g[" << d <<
"]" << std::endl;
564 solve_ok = solve_ok && solver_u2.
isValid();
571 for (
Dimension d = 0; d < u2.size(); ++d )
573 trace.
info() <<
"build metric u2" << std::endl;
578 if (
verbose >= 2 )
trace.
info() <<
"Prefactoring matrix V associated to v" << std::endl;
585 solve_ok = solve_ok && solver_v0.
isValid();
606 double n_oo_max = 1e-4,
607 unsigned int iter_max = 10 )
611 std::ostringstream sstr;
612 sstr <<
"******* Solving AT for epsilon = " << eps <<
" **********";
616 for (
unsigned int i = 0; i < iter_max; ++i )
620 << i <<
"/" << iter_max <<
" ---------------" << std::endl;
625 trace.
info() <<
"Variation |v^k+1 - v^k|_oo = " << std::get<0>( diffs_v )
628 trace.
info() <<
"Variation |v^k+1 - v^k|_2 = " << std::get<1>( diffs_v )
630 trace.
info() <<
"Variation |v^k+1 - v^k|_1 = " << std::get<2>( diffs_v )
634 if ( std::get<0>( diffs_v ) < 1e-4 )
break;
659 bool compute_smallest_epsilon_map =
false,
660 double n_oo_max = 1e-4,
661 unsigned int iter_max = 10 )
664 if ( epsr <= 1.0 ) epsr = 2.0;
668 for (
double eps = eps1; eps >= eps2; eps /= epsr )
671 if ( compute_smallest_epsilon_map )
685 for (
Index index = 0; index <
size( 2 ); index++)
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;
701 std::tuple<double,double,double>
diffV0()
const
708 return std::make_tuple( n_oo, n_2, n_1 );
722 std::tuple<double,double,double>
checkV0()
const
724 const double m1 = v0.myContainer.minCoeff();
725 const double m2 = v0.myContainer.maxCoeff();
726 const double ma = v0.myContainer.mean();
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 );
769 template <
typename VectorFieldOutput,
typename SurfelRangeConstIterator>
772 SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
776 for (
auto it = itB; it != itE; ++it, ++i )
779 ASSERT( output[ i ].
size() >= N );
781 output[ i ][ k ] = u2[ k ].myContainer( idx );
795 template <
typename ScalarFieldOutput,
typename SurfelRangeConstIterator>
798 SurfelRangeConstIterator itB, SurfelRangeConstIterator itE )
800 ASSERT( u2.size() == 1 &&
"[ATSolver2D::getOutputScalarFieldU2] "
801 "You try to output a scalar field from a vector field." );
803 for (
auto it = itB; it != itE; ++it, ++i )
806 output[ i ] = u2[ 0 ].myContainer( idx );
821 template <
typename ScalarFieldOutput,
typename CellRangeConstIterator>
824 CellRangeConstIterator itB, CellRangeConstIterator itE,
833 for (
auto it = itB; it != itE; ++it, ++i )
835 const Cell pointel = *it;
837 output[ i ] = v0.myContainer( idx );
842 for (
auto it = itB; it != itE; ++it, ++i )
844 const Cell linel = *it;
851 case CellOutputPolicy::Average: output[ i ] = 0.5 * ( v0.myContainer( idx0 ) + v0.myContainer( idx1 ) );
853 case CellOutputPolicy::Minimum: output[ i ] = std::min( v0.myContainer( idx0 ), v0.myContainer( idx1 ) );
855 case CellOutputPolicy::Maximum: output[ i ] =
std::max( v0.myContainer( idx0 ), v0.myContainer( idx1 ) );
862 for (
auto it = itB; it != itE; ++it, ++i )
864 const Cell face = *it;
878 case CellOutputPolicy::Average:
879 output[ i ] = 0.25 * ( v0.myContainer( idx00 ) + v0.myContainer( idx01 )
880 + v0.myContainer( idx10 ) + v0.myContainer( idx11 ) );
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 ) ) );
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 ) ) );
905 for (
const SCell surfel :
ptrCalculus->template getIndexedSCells<2, PRIMAL>() )
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() );
926 if ( features[ 1 ] <= threshold )
930 it->second = std::min(
epsilon, it->second );
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;
978 if (
verbose >= 2 )
trace.
info() <<
"derivative of primal 0-forms: primal_D0" << std::endl;
980 if (
verbose >= 2 )
trace.
info() <<
"derivative of primal 1-forms: primal_D1" << std::endl;
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;
987 if (
verbose >= 2 )
trace.
info() <<
"edge to face average operator: M12" << std::endl;
1007 template <
typename T>
1009 operator<< ( std::ostream & out,
const ATSolver2D<T> &
object );
1020 #endif // !defined ATSolver2D_h
1022 #undef ATSolver2D_RECURSES
1023 #endif // else defined(ATSolver2D_RECURSES)