1#if !defined(__DEC_TESTS_COMMON_H__)
2#define __DEC_TESTS_COMMON_H__
6#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
7#include "DGtal/base/Common.h"
8#include "DGtal/helpers/StdDefs.h"
9#include "DGtal/math/linalg/EigenSupport.h"
10#include "DGtal/dec/DiscreteExteriorCalculus.h"
11#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
12#include "DGtal/dec/CDiscreteExteriorCalculusVectorSpace.h"
13#include "DGtal/kernel/sets/DigitalSetBySTLSet.h"
15template <
typename Container>
17is_all_zero(
const Container& container)
19 for (
typename Container::Index ii=0; ii<container.rows(); ii++)
20 for (
typename Container::Index jj=0; jj<container.cols(); jj++)
21 if (container.coeff(ii,jj) != 0)
26template <
typename Container>
28equal(
const Container& aa,
const Container& bb)
30 if (aa.rows() != bb.rows())
return false;
31 if (aa.cols() != bb.cols())
return false;
32 for (
typename Container::Index ii=0; ii<aa.rows(); ii++)
33 for (
typename Container::Index jj=0; jj<aa.cols(); jj++)
34 if (aa.coeff(ii,jj) != bb.coeff(ii,jj))
39template <
typename Container,
typename Value>
41is_identity(
const Container& container,
const Value& value)
43 for (
typename Container::Index ii=0; ii<container.rows(); ii++)
44 for (
typename Container::Index jj=0; jj<container.cols(); jj++)
46 const Value foo =
static_cast<Value>(container.coeff(ii,jj));
47 if ((ii != jj && foo != 0) || (ii == jj && foo != value))
54template <
typename Calculus,
int order>
59 static bool test(
const Calculus& calculus)
65 PrimalIdentity primal_identity = calculus.template identity<order, DGtal::PRIMAL>();
66 if (!is_identity(primal_identity.myContainer, 1))
return false;
69 SolveForm input(calculus);
70 SolveForm output = primal_identity * input;
71 typedef typename Calculus::LinearAlgebraBackend LinearAlgebraBackend;
72 typedef typename LinearAlgebraBackend::SolverConjugateGradient LinearSolver;
75 SolveForm input_solved = solver.
compute(primal_identity).
solve(output);
79 DualIdentity dual_identity = calculus.template identity<order, DGtal::DUAL>();
80 if (!is_identity(dual_identity.myContainer, 1))
return false;
85 const PrimalHodge primal_hodge = calculus.template hodge<order, DGtal::PRIMAL>();
86 const DualHodge dual_hodge= calculus.template hodge<Calculus::dimensionEmbedded-order,
DGtal::DUAL>();
88 DGtal::trace.
info() <<
"testing primal to primal hodge composition order " << order << std::endl;
92 PrimalPrimal primal_primal = dual_hodge * primal_hodge;
93 if (!is_identity(primal_primal.myContainer, pow(-1, order*(Calculus::dimensionEmbedded-order))))
return false;
96 DGtal::trace.
info() <<
"testing dual to dual hodge composition order " << order << std::endl;
100 DualDual dual_dual = primal_hodge * dual_hodge;
101 if (!is_identity(dual_dual.myContainer, pow(-1, order*(Calculus::dimensionEmbedded-order))))
return false;
108template <
typename Calculus>
111 static bool test(
const Calculus& calculus)
113 boost::ignore_unused_variable_warning( calculus );
118template <
typename DigitalSet,
typename LinearAlgebraBackend>
120test_hodge(
int domain_size)
133 if (std::rand()%4!=0)
continue;
135 set.insertNew(point);
146 typename Calculus::Properties properties = calculus.getProperties();
147 DGtal::trace.
info() <<
"properties.size()=" << properties.size() << std::endl;
151 typedef typename Calculus::Cell
Cell;
152 typedef typename Calculus::SCell
SCell;
153 typedef typename Calculus::Index
Index;
154 typedef typename Calculus::KSpace
KSpace;
156 bool test_result =
true;
157 for (
ConstIterator iter = calculus.begin(), iter_end = calculus.end(); test_result && iter!=iter_end; iter++)
159 const Cell& cell = iter->first;
160 const Index& index = calculus.getCellIndex(cell);
161 test_result &= (iter->second.index == index);
163 const SCell& primal_signed_cell = calculus.getSCell(calculus.myKSpace.uDim(cell),
DGtal::PRIMAL, index);
164 test_result &= (signed_cell == primal_signed_cell);
165 const SCell& dual_signed_cell = calculus.getSCell(Calculus::dimensionEmbedded-calculus.myKSpace.uDim(cell),
DGtal::DUAL, index);
166 test_result &= (signed_cell == dual_signed_cell);
170 FATAL_ERROR(test_result);
176 const typename Calculus::PrimalIdentity0 primal_laplace = calculus.template laplace<DGtal::PRIMAL>();
177 DGtal::trace.
info() <<
"primal_laplace_trace=" << primal_laplace.myContainer.diagonal().sum() << std::endl;
178 FATAL_ERROR( ( primal_laplace.myContainer.diagonal().array() >= 0 ).prod() ==
true );
180 const typename Calculus::DualIdentity0 dual_laplace = calculus.template laplace<DGtal::DUAL>();
181 DGtal::trace.
info() <<
"dual_laplace_trace=" << dual_laplace.myContainer.diagonal().sum() << std::endl;
182 FATAL_ERROR( ( dual_laplace.myContainer.diagonal().array() >= 0 ).prod() ==
true );
191 FATAL_ERROR(test_result);
195template <
typename Calculus,
int order>
200 static bool test(
const Calculus& calculus)
202 DGtal::trace.
info() <<
"testing primal derivative composition order " << order << std::endl;
206 FirstDerivative first_derivative = calculus.template derivative<order, DGtal::PRIMAL>();
208 SecondDerivative second_derivative = calculus.template derivative<order+1, DGtal::PRIMAL>();
210 DoubleDerivative double_derivative = second_derivative * first_derivative;
211 if (!is_all_zero(double_derivative.myContainer))
return false;
214 DGtal::trace.
info() <<
"testing dual derivative composition order " << order << std::endl;
218 FirstDerivative first_derivative = calculus.template derivative<order, DGtal::DUAL>();
220 SecondDerivative second_derivative = calculus.template derivative<order+1, DGtal::DUAL>();
222 DoubleDerivative double_derivative = second_derivative * first_derivative;
223 if (!is_all_zero(double_derivative.myContainer))
return false;
253template <
typename Calculus>
256 static bool test(
const Calculus& )
262template <
typename DigitalSet,
typename LinearAlgebraBackend>
264test_derivative(
int domain_size)
277 if (std::rand()%3!=0)
continue;
279 set.insertNew(point);
289 Calculus calculus = CalculusFactory::createFromDigitalSet(set,
false);
291 typename Calculus::Properties properties = calculus.getProperties();
292 DGtal::trace.
info() <<
"properties.size()=" << properties.size() << std::endl;
295 FATAL_ERROR(test_result);
302 Calculus calculus = CalculusFactory::createFromDigitalSet(set,
true);
304 typename Calculus::Properties properties = calculus.getProperties();
305 DGtal::trace.
info() <<
"properties.size()=" << properties.size() << std::endl;
308 FATAL_ERROR(test_result);
314template <
typename LinearAlgebraBackend>
429template <
typename LinearAlgebraBackend>
441 const Calculus calculus = CalculusFactory::createFromDigitalSet(set);
445 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,0)),
DGtal::PRIMAL ) == 1 );
446 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,0)),
DGtal::DUAL ) == 1 );
448 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,0)),
DGtal::PRIMAL ) == 1 );
449 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,0)),
DGtal::DUAL ) == -1 );
451 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,1)),
DGtal::PRIMAL ) == 1 );
452 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,1)),
DGtal::DUAL ) == -1 );
454 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,1)),
DGtal::PRIMAL ) == 1 );
455 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,1)),
DGtal::DUAL ) == 1 );
462 const Calculus calculus = CalculusFactory::createFromDigitalSet(set);
466 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,0,0)),
DGtal::PRIMAL ) == 1 );
467 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,0,0)),
DGtal::DUAL ) == 1 );
469 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,0,0)),
DGtal::PRIMAL ) == 1 );
470 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,0,0)),
DGtal::DUAL ) == 1 );
471 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,1,0)),
DGtal::PRIMAL ) == 1 );
472 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,1,0)),
DGtal::DUAL ) == 1 );
473 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,0,1)),
DGtal::PRIMAL ) == 1 );
474 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,0,1)),
DGtal::DUAL ) == 1 );
476 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,1,0)),
DGtal::PRIMAL ) == 1 );
477 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,1,0)),
DGtal::DUAL ) == 1 );
478 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,1,1)),
DGtal::PRIMAL ) == 1 );
479 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(0,1,1)),
DGtal::DUAL ) == 1 );
480 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,0,1)),
DGtal::PRIMAL ) == 1 );
481 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,0,1)),
DGtal::DUAL ) == 1 );
483 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,1,1)),
DGtal::PRIMAL ) == 1 );
484 FATAL_ERROR( calculus.hodgeSign( calculus.myKSpace.uFirst(
Point(1,1,1)),
DGtal::DUAL ) == 1 );
501template <
typename LinearAlgebraBackend>
503test_backend(
const int& ntime,
const int& maxdim)
507 test_hodge_sign<LinearAlgebraBackend>();
509 for (
int kk=0; kk<ntime; kk++)
532 if (maxdim>=1) test_hodge<DigitalSet1, LinearAlgebraBackend>(10);
533 if (maxdim>=2) test_hodge<DGtal::Z2i::DigitalSet, LinearAlgebraBackend>(5);
534 if (maxdim>=3) test_hodge<DGtal::Z3i::DigitalSet, LinearAlgebraBackend>(5);
535 if (maxdim>=4) test_hodge<DigitalSet4, LinearAlgebraBackend>(4);
536 if (maxdim>=5) test_hodge<DigitalSet5, LinearAlgebraBackend>(3);
537 if (maxdim>=6) test_hodge<DigitalSet6, LinearAlgebraBackend>(3);
538 if (maxdim>=7) test_hodge<DigitalSet7, LinearAlgebraBackend>(2);
542 if (maxdim>=1) test_derivative<DigitalSet1, LinearAlgebraBackend>(10);
543 if (maxdim>=2) test_derivative<DGtal::Z2i::DigitalSet, LinearAlgebraBackend>(5);
544 if (maxdim>=3) test_derivative<DGtal::Z3i::DigitalSet, LinearAlgebraBackend>(5);
545 if (maxdim>=4) test_derivative<DigitalSet4, LinearAlgebraBackend>(4);
546 if (maxdim>=5) test_derivative<DigitalSet5, LinearAlgebraBackend>(3);
547 if (maxdim>=6) test_derivative<DigitalSet6, LinearAlgebraBackend>(3);
548 if (maxdim>=7) test_derivative<DigitalSet7, LinearAlgebraBackend>(2);
552 test_concepts<LinearAlgebraBackend>();
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
TDomain Domain
Domain type.
Aim: A container class for storing sets of digital points within some given domain.
Aim: This class provides static members to create DEC structures from various other DGtal structures.
static DiscreteExteriorCalculus< TDigitalSet::Point::dimension, TDigitalSet::Point::dimension, TLinearAlgebraBackend, TInteger > createFromDigitalSet(const TDigitalSet &set, const bool add_border=true)
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
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...
static const constexpr Sign NEG
static const constexpr Sign POS
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
void beginBlock(const std::string &keyword="")
MyDigitalSurface::ConstIterator ConstIterator
Aim: Represents a set of points within the given domain. This set of points is modifiable by the user...
Aim: Lift linear algebra container concept into the dec package.
static bool test(const Calculus &)
static bool test(const Calculus &calculus)
BOOST_STATIC_ASSERT((order<(int) Calculus::dimensionEmbedded - 1))
static bool test(const Calculus &calculus)
static bool test(const Calculus &calculus)
BOOST_STATIC_ASSERT((order<=Calculus::dimensionEmbedded))
bool test(const I &itb, const I &ite)
HyperRectDomain< Space > Domain
Z2i::DigitalSet DigitalSet