DGtal  1.5.beta
DECCommon.h
1 #if !defined(__DEC_TESTS_COMMON_H__)
2 #define __DEC_TESTS_COMMON_H__
3 
4 #include <list>
5 
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"
14 
15 template <typename Container>
16 bool
17 is_all_zero(const Container& container)
18 {
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)
22  return false;
23  return true;
24 }
25 
26 template <typename Container>
27 bool
28 equal(const Container& aa, const Container& bb)
29 {
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))
35  return false;
36  return true;
37 }
38 
39 template <typename Container, typename Value>
40 bool
41 is_identity(const Container& container, const Value& value)
42 {
43  for (typename Container::Index ii=0; ii<container.rows(); ii++)
44  for (typename Container::Index jj=0; jj<container.cols(); jj++)
45  {
46  const Value foo = static_cast<Value>(container.coeff(ii,jj));
47  if ((ii != jj && foo != 0) || (ii == jj && foo != value))
48  return false;
49  }
50  return true;
51 }
52 
53 //DC: changing order to int (to match with -1 specialization)
54 template <typename Calculus, int order>
56 {
57  BOOST_STATIC_ASSERT(( order <= Calculus::dimensionEmbedded ));
58 
59  static bool test(const Calculus& calculus)
60  {
61  DGtal::trace.info() << "testing identity operators " << order << std::endl;
62 
63  { // test identity operator
65  PrimalIdentity primal_identity = calculus.template identity<order, DGtal::PRIMAL>();
66  if (!is_identity(primal_identity.myContainer, 1)) return false;
67 
69  SolveForm input(calculus);
70  SolveForm output = primal_identity * input;
71  typedef typename Calculus::LinearAlgebraBackend LinearAlgebraBackend;
72  typedef typename LinearAlgebraBackend::SolverConjugateGradient LinearSolver;
74  Solver solver;
75  SolveForm input_solved = solver.compute(primal_identity).solve(output);
76  //if (input_solved != input) return false;
77 
79  DualIdentity dual_identity = calculus.template identity<order, DGtal::DUAL>();
80  if (!is_identity(dual_identity.myContainer, 1)) return false;
81  }
82 
83  typedef DGtal::LinearOperator<Calculus, order, DGtal::PRIMAL, Calculus::dimensionEmbedded-order, DGtal::DUAL> PrimalHodge;
84  typedef DGtal::LinearOperator<Calculus, Calculus::dimensionEmbedded-order, DGtal::DUAL, order, DGtal::PRIMAL> DualHodge;
85  const PrimalHodge primal_hodge = calculus.template hodge<order, DGtal::PRIMAL>();
86  const DualHodge dual_hodge= calculus.template hodge<Calculus::dimensionEmbedded-order, DGtal::DUAL>();
87 
88  DGtal::trace.info() << "testing primal to primal hodge composition order " << order << std::endl;
89 
90  { // test primal to primal composition
92  PrimalPrimal primal_primal = dual_hodge * primal_hodge;
93  if (!is_identity(primal_primal.myContainer, pow(-1, order*(Calculus::dimensionEmbedded-order)))) return false;
94  }
95 
96  DGtal::trace.info() << "testing dual to dual hodge composition order " << order << std::endl;
97 
98  { // test dual to dual composition
99  typedef DGtal::LinearOperator<Calculus, Calculus::dimensionEmbedded-order, DGtal::DUAL, Calculus::dimensionEmbedded-order, DGtal::DUAL> DualDual;
100  DualDual dual_dual = primal_hodge * dual_hodge;
101  if (!is_identity(dual_dual.myContainer, pow(-1, order*(Calculus::dimensionEmbedded-order)))) return false;
102  }
103 
105  }
106 };
107 
108 template <typename Calculus>
109 struct HodgeTester<Calculus, -1>
110 {
111  static bool test(const Calculus& calculus)
112  {
113  boost::ignore_unused_variable_warning( calculus );
114  return true;
115  }
116 };
117 
118 template <typename DigitalSet, typename LinearAlgebraBackend>
119 void
120 test_hodge(int domain_size)
121 {
122  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDigitalSet<DigitalSet> ));
123 
124  typedef typename DigitalSet::Domain Domain;
125  typedef typename DigitalSet::Point Point;
126  DGtal::trace.info() << "dimension=" << Point::dimension << std::endl;
127  Domain domain(Point(), Point::diagonal(domain_size-1));
128  DGtal::trace.info() << "domain=" << domain << std::endl;
129 
130  DigitalSet set(domain);
131  for (typename Domain::ConstIterator di=domain.begin(), die=domain.end(); di!=die; di++)
132  {
133  if (std::rand()%4!=0) continue;
134  const typename Domain::Point& point = *di;
135  set.insertNew(point);
136  }
137  DGtal::trace.info() << "domain.size()=" << domain.size() << std::endl;
138  DGtal::trace.info() << "set.size()=" << set.size() << std::endl;
139 
142  {
143  DGtal::trace.beginBlock("testing indexes");
144 
145  {
146  typename Calculus::Properties properties = calculus.getProperties();
147  DGtal::trace.info() << "properties.size()=" << properties.size() << std::endl;
148  }
149 
150  typedef typename Calculus::ConstIterator ConstIterator;
151  typedef typename Calculus::Cell Cell;
152  typedef typename Calculus::SCell SCell;
153  typedef typename Calculus::Index Index;
154  typedef typename Calculus::KSpace KSpace;
155 
156  bool test_result = true;
157  for (ConstIterator iter = calculus.begin(), iter_end = calculus.end(); test_result && iter!=iter_end; iter++)
158  {
159  const Cell& cell = iter->first;
160  const Index& index = calculus.getCellIndex(cell);
161  test_result &= (iter->second.index == index);
162  const SCell& signed_cell = calculus.myKSpace.signs(cell, iter->second.flipped ? KSpace::NEG : KSpace::POS);
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);
167  }
169 
170  FATAL_ERROR(test_result);
171  }
172 
173  {
174  DGtal::trace.beginBlock("testing laplace sign");
175 
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 );
179 
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 );
183 
185  }
186 
187  DGtal::trace.beginBlock("testing hodge");
190 
191  FATAL_ERROR(test_result);
192 }
193 
194 //DC Order->int (see above)
195 template <typename Calculus, int order>
197 {
198  BOOST_STATIC_ASSERT(( order < (int)Calculus::dimensionEmbedded - 1 ));
199 
200  static bool test(const Calculus& calculus)
201  {
202  DGtal::trace.info() << "testing primal derivative composition order " << order << std::endl;
203 
204  { // test primal composition
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;
212  }
213 
214  DGtal::trace.info() << "testing dual derivative composition order " << order << std::endl;
215 
216  { // test dual composition
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;
224  }
225 
226  /*
227  DGtal::trace.info() << "testing liebnitz rule order " << order << std::endl;
228 
229  {
230  typedef DGtal::LinearOperator<Calculus, order, DGtal::PRIMAL, order+1, DGtal::PRIMAL> Derivative;
231  Derivative derivative = calculus.template derivative<order, DGtal::PRIMAL>();
232 
233  typedef DGtal::KForm<Calculus, order, DGtal::PRIMAL> InputForm;
234  typedef DGtal::KForm<Calculus, order+1, DGtal::PRIMAL> OutputForm;
235  InputForm alpha(calculus), beta(calculus), gamma(calculus);
236 
237  for (int kk=0; kk<calculus.kFormLength(order, DGtal::PRIMAL); kk++)
238  {
239  const double ak = static_cast<double>(random())/RAND_MAX;
240  const double bk = static_cast<double>(random())/RAND_MAX;
241  alpha.myContainer(kk) = ak;
242  beta.myContainer(kk) = bk;
243  gamma.myContainer(kk) = ak*bk;
244  }
245 
246  }
247  */
248 
250  }
251 };
252 
253 template <typename Calculus>
254 struct DerivativeTester<Calculus, -1>
255 {
256  static bool test(const Calculus& )
257  {
258  return true;
259  }
260 };
261 
262 template <typename DigitalSet, typename LinearAlgebraBackend>
263 void
264 test_derivative(int domain_size)
265 {
266  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDigitalSet<DigitalSet> ));
267 
268  typedef typename DigitalSet::Domain Domain;
269  typedef typename DigitalSet::Point Point;
270  DGtal::trace.info() << "dimension=" << Point::dimension << std::endl;
271  Domain domain(Point(), Point::diagonal(domain_size-1));
272  DGtal::trace.info() << "domain=" << domain << std::endl;
273 
274  DigitalSet set(domain);
275  for (typename Domain::ConstIterator di=domain.begin(), die=domain.end(); di!=die; di++)
276  {
277  if (std::rand()%3!=0) continue;
278  const typename Domain::Point& point = *di;
279  set.insertNew(point);
280  }
281  DGtal::trace.info() << "domain.size()=" << domain.size() << std::endl;
282  DGtal::trace.info() << "set.size()=" << set.size() << std::endl;
283 
286 
287  {
288  DGtal::trace.beginBlock("testing derivative without border");
289  Calculus calculus = CalculusFactory::createFromDigitalSet(set, false);
290 
291  typename Calculus::Properties properties = calculus.getProperties();
292  DGtal::trace.info() << "properties.size()=" << properties.size() << std::endl;
293 
294  bool test_result = DerivativeTester<Calculus, (int)Calculus::dimensionEmbedded-2>::test(calculus);
295  FATAL_ERROR(test_result);
296 
298  }
299 
300  {
301  DGtal::trace.beginBlock("testing derivative with border");
302  Calculus calculus = CalculusFactory::createFromDigitalSet(set, true);
303 
304  typename Calculus::Properties properties = calculus.getProperties();
305  DGtal::trace.info() << "properties.size()=" << properties.size() << std::endl;
306 
307  bool test_result = DerivativeTester<Calculus, (int)Calculus::dimensionEmbedded-2>::test(calculus);
308  FATAL_ERROR(test_result);
309 
311  }
312 }
313 
314 template <typename LinearAlgebraBackend>
315 void
316 test_concepts()
317 {
318  DGtal::trace.beginBlock("concepts");
319 
320  //FIXME add 1d
321 
322  { // 2d
330 
333 
338 
345 
352  }
353 
354  { // 2d embedded in 3d
362 
365 
370 
377 
384  }
385 
386  { // 3d
396 
399 
406 
415 
424  }
425 
427 }
428 
429 template <typename LinearAlgebraBackend>
430 void
431 test_hodge_sign()
432 {
434 
435  DGtal::trace.beginBlock("testing hodge sign");
436 
437  {
440  const DGtal::Z2i::DigitalSet set(domain);
441  const Calculus calculus = CalculusFactory::createFromDigitalSet(set);
442  typedef DGtal::Z2i::Point Point;
443 
444  // primal point, dual cell
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 );
447  // primal horizontal edge, dual vertical edge
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 );
450  // primal vectical edge, dual horizontal edge
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 );
453  // primal cell, dual point
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 );
456  }
457 
458  {
461  const DGtal::Z3i::DigitalSet set(domain);
462  const Calculus calculus = CalculusFactory::createFromDigitalSet(set);
463  typedef DGtal::Z3i::Point Point;
464 
465  // primal point, dual cell
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 );
468  // primal edge, dual surfel
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 );
475  // primal surfel, dual edge
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 );
482  // primal cell, dual point
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 );
485  }
486 
488 }
489 
490 void
491 test_duality()
492 {
493  DGtal::trace.beginBlock("testing duality");
494 
497 
499 }
500 
501 template <typename LinearAlgebraBackend>
502 void
503 test_backend(const int& ntime, const int& maxdim)
504 {
505  test_duality();
506 
507  test_hodge_sign<LinearAlgebraBackend>();
508 
509  for (int kk=0; kk<ntime; kk++)
510  {
511  typedef DGtal::SpaceND<1, int> Space1;
512  typedef DGtal::HyperRectDomain<Space1> Domain1;
513  typedef DGtal::DigitalSetBySTLSet<Domain1> DigitalSet1;
514 
515  typedef DGtal::SpaceND<4, int> Space4;
516  typedef DGtal::HyperRectDomain<Space4> Domain4;
517  typedef DGtal::DigitalSetBySTLSet<Domain4> DigitalSet4;
518 
519  typedef DGtal::SpaceND<5, int> Space5;
520  typedef DGtal::HyperRectDomain<Space5> Domain5;
521  typedef DGtal::DigitalSetBySTLSet<Domain5> DigitalSet5;
522 
523  typedef DGtal::SpaceND<6, int> Space6;
524  typedef DGtal::HyperRectDomain<Space6> Domain6;
525  typedef DGtal::DigitalSetBySTLSet<Domain6> DigitalSet6;
526 
527  typedef DGtal::SpaceND<7, int> Space7;
528  typedef DGtal::HyperRectDomain<Space7> Domain7;
529  typedef DGtal::DigitalSetBySTLSet<Domain7> DigitalSet7;
530 
531  DGtal::trace.beginBlock("testing hodges");
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);
540 
541  DGtal::trace.beginBlock("testing derivatives");
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);
550  }
551 
552  test_concepts<LinearAlgebraBackend>();
553 }
554 
555 #endif
556 
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
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...
Aim: KForm represents discrete kforms in the dec package.
Definition: KForm.h:66
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
PolyCalculus * calculus
Z3i::SCell SCell
SMesh::Index Index
MyDigitalSurface::ConstIterator ConstIterator
Trace trace
Definition: Common.h:153
@ PRIMAL
Definition: Duality.h:61
@ DUAL
Definition: Duality.h:62
Aim: Represents a set of points within the given domain. This set of points is modifiable by the user...
Definition: CDigitalSet.h:141
Aim: Lift linear algebra container concept into the dec package.
static bool test(const Calculus &)
Definition: DECCommon.h:256
static bool test(const Calculus &calculus)
Definition: DECCommon.h:200
BOOST_STATIC_ASSERT((order<(int) Calculus::dimensionEmbedded - 1))
static bool test(const Calculus &calculus)
Definition: DECCommon.h:111
static bool test(const Calculus &calculus)
Definition: DECCommon.h:59
BOOST_STATIC_ASSERT((order<=Calculus::dimensionEmbedded))
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition: testBits.cpp:44
MyPointD Point
Definition: testClone2.cpp:383
KSpace::Cell Cell
bool test(const I &itb, const I &ite)
Domain domain
HyperRectDomain< Space > Domain
Z2i::DigitalSet DigitalSet