DGtal  1.4.2
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>
55 struct HodgeTester
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
64  typedef DGtal::LinearOperator<Calculus, order, DGtal::PRIMAL, order, DGtal::PRIMAL> PrimalIdentity;
65  PrimalIdentity primal_identity = calculus.template identity<order, DGtal::PRIMAL>();
66  if (!is_identity(primal_identity.myContainer, 1)) return false;
67 
68  typedef DGtal::KForm<Calculus, order, DGtal::PRIMAL> SolveForm;
69  SolveForm input(calculus);
70  SolveForm output = primal_identity * input;
71  typedef typename Calculus::LinearAlgebraBackend LinearAlgebraBackend;
72  typedef typename LinearAlgebraBackend::SolverConjugateGradient LinearSolver;
73  typedef DGtal::DiscreteExteriorCalculusSolver<Calculus, LinearSolver, order, DGtal::PRIMAL, order, DGtal::PRIMAL> Solver;
74  Solver solver;
75  SolveForm input_solved = solver.compute(primal_identity).solve(output);
76  //if (input_solved != input) return false;
77 
78  typedef DGtal::LinearOperator<Calculus, order, DGtal::DUAL, order, DGtal::DUAL> DualIdentity;
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
91  typedef DGtal::LinearOperator<Calculus, order, DGtal::PRIMAL, order, DGtal::PRIMAL> PrimalPrimal;
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 
104  return HodgeTester<Calculus, order-1>::test(calculus);
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 
140  typedef DGtal::DiscreteExteriorCalculus<Domain::Space::dimension, Domain::Space::dimension, LinearAlgebraBackend> Calculus;
141  Calculus calculus = DGtal::DiscreteExteriorCalculusFactory<LinearAlgebraBackend>::createFromDigitalSet(set);
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  }
168  DGtal::trace.endBlock();
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 
184  DGtal::trace.endBlock();
185  }
186 
187  DGtal::trace.beginBlock("testing hodge");
188  bool test_result = HodgeTester<Calculus, Calculus::dimensionEmbedded>::test(calculus);
189  DGtal::trace.endBlock();
190 
191  FATAL_ERROR(test_result);
192 }
193 
194 //DC Order->int (see above)
195 template <typename Calculus, int order>
196 struct DerivativeTester
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
205  typedef DGtal::LinearOperator<Calculus, order, DGtal::PRIMAL, order+1, DGtal::PRIMAL> FirstDerivative;
206  FirstDerivative first_derivative = calculus.template derivative<order, DGtal::PRIMAL>();
207  typedef DGtal::LinearOperator<Calculus, order+1, DGtal::PRIMAL, order+2, DGtal::PRIMAL> SecondDerivative;
208  SecondDerivative second_derivative = calculus.template derivative<order+1, DGtal::PRIMAL>();
209  typedef DGtal::LinearOperator<Calculus, order, DGtal::PRIMAL, order+2, DGtal::PRIMAL> DoubleDerivative;
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
217  typedef DGtal::LinearOperator<Calculus, order, DGtal::DUAL, order+1, DGtal::DUAL> FirstDerivative;
218  FirstDerivative first_derivative = calculus.template derivative<order, DGtal::DUAL>();
219  typedef DGtal::LinearOperator<Calculus, order+1, DGtal::DUAL, order+2, DGtal::DUAL> SecondDerivative;
220  SecondDerivative second_derivative = calculus.template derivative<order+1, DGtal::DUAL>();
221  typedef DGtal::LinearOperator<Calculus, order, DGtal::DUAL, order+2, DGtal::DUAL> DoubleDerivative;
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 
249  return DerivativeTester<Calculus, order-1>::test(calculus);
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 
284  typedef DGtal::DiscreteExteriorCalculusFactory<LinearAlgebraBackend> CalculusFactory;
285  typedef DGtal::DiscreteExteriorCalculus<Domain::Space::dimension, Domain::Space::dimension, LinearAlgebraBackend> Calculus;
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 
297  DGtal::trace.endBlock();
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 
310  DGtal::trace.endBlock();
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
323  typedef DGtal::DiscreteExteriorCalculus<2, 2, LinearAlgebraBackend> Calculus;
324  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm0> ));
325  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm1> ));
326  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm2> ));
327  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm0> ));
328  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm1> ));
329  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm2> ));
330 
331  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalVectorField> ));
332  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualVectorField> ));
333 
334  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative0> ));
335  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative1> ));
336  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative0> ));
337  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative1> ));
338 
339  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge0> ));
340  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge1> ));
341  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge2> ));
342  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge0> ));
343  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge1> ));
344  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge2> ));
345 
346  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity0> ));
347  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity1> ));
348  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity2> ));
349  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity0> ));
350  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity1> ));
351  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity2> ));
352  }
353 
354  { // 2d embedded in 3d
355  typedef DGtal::DiscreteExteriorCalculus<2, 3, LinearAlgebraBackend> Calculus;
356  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm0> ));
357  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm1> ));
358  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm2> ));
359  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm0> ));
360  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm1> ));
361  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm2> ));
362 
363  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalVectorField> ));
364  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualVectorField> ));
365 
366  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative0> ));
367  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative1> ));
368  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative0> ));
369  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative1> ));
370 
371  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge0> ));
372  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge1> ));
373  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge2> ));
374  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge0> ));
375  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge1> ));
376  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge2> ));
377 
378  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity0> ));
379  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity1> ));
380  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity2> ));
381  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity0> ));
382  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity1> ));
383  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity2> ));
384  }
385 
386  { // 3d
387  typedef DGtal::DiscreteExteriorCalculus<3, 3, LinearAlgebraBackend> Calculus;
388  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm0> ));
389  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm1> ));
390  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm2> ));
391  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalForm3> ));
392  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm0> ));
393  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm1> ));
394  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm2> ));
395  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualForm3> ));
396 
397  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalVectorField> ));
398  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualVectorField> ));
399 
400  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative0> ));
401  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative1> ));
402  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalDerivative2> ));
403  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative0> ));
404  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative1> ));
405  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualDerivative2> ));
406 
407  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge0> ));
408  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge1> ));
409  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge2> ));
410  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalHodge3> ));
411  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge0> ));
412  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge1> ));
413  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge2> ));
414  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualHodge3> ));
415 
416  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity0> ));
417  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity1> ));
418  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity2> ));
419  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::PrimalIdentity3> ));
420  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity0> ));
421  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity1> ));
422  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity2> ));
423  BOOST_CONCEPT_ASSERT(( DGtal::concepts::CDiscreteExteriorCalculusVectorSpace<typename Calculus::DualIdentity3> ));
424  }
425 
426  DGtal::trace.endBlock();
427 }
428 
429 template <typename LinearAlgebraBackend>
430 void
431 test_hodge_sign()
432 {
433  typedef DGtal::DiscreteExteriorCalculusFactory<LinearAlgebraBackend> CalculusFactory;
434 
435  DGtal::trace.beginBlock("testing hodge sign");
436 
437  {
438  typedef DGtal::DiscreteExteriorCalculus<2, 2, LinearAlgebraBackend> Calculus;
439  const DGtal::Z2i::Domain domain;
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  {
459  typedef DGtal::DiscreteExteriorCalculus<3, 3, LinearAlgebraBackend> Calculus;
460  const DGtal::Z3i::Domain domain;
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 
487  DGtal::trace.endBlock();
488 }
489 
490 void
491 test_duality()
492 {
493  DGtal::trace.beginBlock("testing duality");
494 
495  BOOST_STATIC_ASSERT(( DGtal::OppositeDuality<DGtal::PRIMAL>::duality == DGtal::DUAL ));
496  BOOST_STATIC_ASSERT(( DGtal::OppositeDuality<DGtal::DUAL>::duality == DGtal::PRIMAL ));
497 
498  DGtal::trace.endBlock();
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);
539  DGtal::trace.endBlock();
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);
549  DGtal::trace.endBlock();
550  }
551 
552  test_concepts<LinearAlgebraBackend>();
553 }
554 
555 #endif
556 
DGtal::DigitalSetByAssociativeContainer
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Definition: DigitalSetByAssociativeContainer.h:90
DGtal::DigitalSetBySTLSet
Aim: A container class for storing sets of digital points within some given domain.
Definition: DigitalSetBySTLSet.h:85
DGtal::DiscreteExteriorCalculusFactory
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Definition: DiscreteExteriorCalculusFactory.h:62
DGtal::DiscreteExteriorCalculusFactory::createFromDigitalSet
static DiscreteExteriorCalculus< TDigitalSet::Point::dimension, TDigitalSet::Point::dimension, TLinearAlgebraBackend, TInteger > createFromDigitalSet(const TDigitalSet &set, const bool add_border=true)
DGtal::DiscreteExteriorCalculusSolver
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
Definition: DiscreteExteriorCalculusSolver.h:70
DGtal::DiscreteExteriorCalculusSolver::solve
SolutionKForm solve(const InputKForm &input_kform) const
DGtal::DiscreteExteriorCalculusSolver::compute
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition: DiscreteExteriorCalculus.h:98
DGtal::KForm
Aim: KForm represents discrete kforms in the dec package.
Definition: KForm.h:66
DGtal::LinearOperator
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
Definition: LinearOperator.h:69
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
calculus
PolyCalculus * calculus
Definition: dgtalCalculus-geodesic.cpp:70
SCell
Z3i::SCell SCell
Definition: digitalPolyhedronBuilder3D.cpp:79
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
ConstIterator
MyDigitalSurface::ConstIterator ConstIterator
Definition: greedy-plane-segmentation-ex2.cpp:93
DGtal::trace
Trace trace
Definition: Common.h:153
DGtal::PRIMAL
@ PRIMAL
Definition: Duality.h:61
DGtal::DUAL
@ DUAL
Definition: Duality.h:62
DGtal::concepts::CDigitalSet
Aim: Represents a set of points within the given domain. This set of points is modifiable by the user...
Definition: CDigitalSet.h:141
DGtal::concepts::CDiscreteExteriorCalculusVectorSpace
Aim: Lift linear algebra container concept into the dec package.
Definition: CDiscreteExteriorCalculusVectorSpace.h:90
DerivativeTester< Calculus, -1 >::test
static bool test(const Calculus &)
Definition: DECCommon.h:256
DerivativeTester::test
static bool test(const Calculus &calculus)
Definition: DECCommon.h:200
DerivativeTester::BOOST_STATIC_ASSERT
BOOST_STATIC_ASSERT((order<(int) Calculus::dimensionEmbedded - 1))
HodgeTester< Calculus, -1 >::test
static bool test(const Calculus &calculus)
Definition: DECCommon.h:111
HodgeTester::test
static bool test(const Calculus &calculus)
Definition: DECCommon.h:59
HodgeTester::BOOST_STATIC_ASSERT
BOOST_STATIC_ASSERT((order<=Calculus::dimensionEmbedded))
index
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition: testBits.cpp:44
Point
MyPointD Point
Definition: testClone2.cpp:383
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
test
bool test(const I &itb, const I &ite)
Definition: testMostCenteredMSEstimator.cpp:71
domain
Domain domain
Definition: testProjection.cpp:88
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
DigitalSet
Z2i::DigitalSet DigitalSet
Definition: testVoronoiMapComplete.cpp:40