DGtal  1.5.beta
DigitalSurfaceConvolver.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file DigitalSurfaceConvolver.ih
19  * @author Jeremy Levallois (\c jeremy.levallois@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), INSA-Lyon, France
21  * LAboratoire de MAthématiques - LAMA (CNRS, UMR 5127), Université de Savoie, France
22  *
23  * @date 2012/03/27
24  *
25  * Implementation of inline methods defined in DigitalSurfaceConvolver.h
26  *
27  * This file is part of the DGtal library.
28  */
29 
30 ///////////////////////////////////////////////////////////////////////////////
31 // IMPLEMENTATION of inline methods.
32 ///////////////////////////////////////////////////////////////////////////////
33 
34 //////////////////////////////////////////////////////////////////////////////
35 #include <cstdlib>
36 //////////////////////////////////////////////////////////////////////////////
37 
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 // IMPLEMENTATION of inline methods.
41 ///////////////////////////////////////////////////////////////////////////////
42 #include "DGtal/geometry/surfaces/DigitalSurfaceConvolver.h"
43 #include "DGtal/kernel/NumberTraits.h"
44 ///////////////////////////////////////////////////////////////////////////////
45 // ----------------------- Standard services ------------------------------
46 
47 
48 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
49 ///////////////////////////////////////////////////// nD /////////////////////////////////////////////////////
50 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
51 
52 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
53 inline
54 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >
55 ::DigitalSurfaceConvolver( ConstAlias< Functor > f,
56  ConstAlias< KernelFunctor > g,
57  ConstAlias< KSpace > space)
58  : myFFunctor( f ),
59  myGFunctor( g ),
60  myKSpace( space ),
61  isInitFullMasks( false ),
62  isInitKernelAndMasks( false )
63 {
64  myEmbedder = Embedder( myKSpace );
65 }
66 
67 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
68 inline
69 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >
70 ::DigitalSurfaceConvolver( const DigitalSurfaceConvolver& other )
71  : myFFunctor( other.myFFunctor ),
72  myGFunctor( other.myGFunctor ),
73  myKSpace( other.myKSpace ),
74  myEmbedder( other.myEmbedder ),
75  isInitFullMasks( other.isInitFullMasks ),
76  isInitKernelAndMasks( other.isInitKernelAndMasks ),
77  myMasks( other.myMasks ),
78  myKernel( other.myKernel ),
79  myKernelMask( other.myKernelMask ),
80  myKernelSpelOrigin( other.myKernelSpelOrigin )
81 {
82 }
83 
84 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
85 inline
86 void
87 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::init
88 ( const Point & pOrigin,
89  ConstAlias< PairIterators > fullKernel,
90  ConstAlias< std::vector< PairIterators > > masks )
91 {
92  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
93  myKernelMask = &fullKernel;
94  myMasks = &masks;
95 
96  // ASSERT ( myMasks->size () == 9 );
97 
98  isInitFullMasks = true;
99  isInitKernelAndMasks = false;
100 }
101 
102 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
103 inline
104 void
105 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::init
106 ( const Point & pOrigin,
107  ConstAlias< DigitalKernel > fullKernel,
108  ConstAlias< std::vector< PairIterators > > masks )
109 {
110  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
111  myKernel = &fullKernel;
112  myMasks = &masks;
113 
114  // ASSERT ( myMasks->size () == 9 );
115 
116  isInitFullMasks = false;
117  isInitKernelAndMasks = true;
118 }
119 
120 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
121 template< typename SurfelIterator >
122 inline
123 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
124 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::eval
125 ( const SurfelIterator & it ) const
126 {
127  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
128 
129  Quantity innerSum, outerSum;
130 
131  core_eval( it, innerSum, outerSum, false );
132 
133  double lambda = 0.5;
134  return ( innerSum * lambda + outerSum * ( 1.0 - lambda ));
135 }
136 
137 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
138 template< typename SurfelIterator, typename EvalFunctor >
139 inline
140 typename EvalFunctor::Value
141 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::eval
142 ( const SurfelIterator & it,
143  EvalFunctor functor ) const
144 {
145  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
146 
147  Quantity innerSum, outerSum;
148  Quantity resultQuantity;
149 
150  core_eval( it, innerSum, outerSum, false );
151 
152  double lambda = 0.5;
153  resultQuantity = innerSum * lambda + outerSum * ( 1.0 - lambda );
154  return functor( resultQuantity );
155 }
156 
157 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
158 template< typename SurfelIterator, typename OutputIterator >
159 inline
160 void
161 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::eval
162 ( const SurfelIterator & itbegin,
163  const SurfelIterator & itend,
164  OutputIterator & result ) const
165 {
166  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
167 
168  Dimension total = 0;
169 #ifdef DEBUG_VERBOSE
170  Dimension recount = 0;
171 #endif
172 
173  Quantity lastInnerSum;
174  Quantity lastOuterSum;
175 
176  Quantity innerSum, outerSum;
177 
178  Spel lastInnerSpel, lastOuterSpel;
179 
180  /// Iterate on all cells
181  for( SurfelIterator it = itbegin; it != itend; ++it )
182  {
183  if( total != 0 )
184  {
185 #ifdef DEBUG_VERBOSE
186  bool hasJumped = core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
187  recount = ( hasJumped ) ? recount + 1 : recount;
188 #else
189  core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
190 #endif
191  }
192  else
193  {
194  core_eval( it, innerSum, outerSum, false, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
195  }
196 
197  double lambda = 0.5;
198  *result++ = ( innerSum * lambda + outerSum * ( 1.0 - lambda ));
199 
200  ++total;
201  }
202 
203 #ifdef DEBUG_VERBOSE
204  std::cout << "#total cells = " << total << std::endl;
205  std::cout << "#recount = " << recount << std::endl;
206 #endif
207 }
208 
209 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
210 template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
211 inline
212 void
213 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::eval
214 ( const SurfelIterator & itbegin,
215  const SurfelIterator & itend,
216  OutputIterator & result,
217  EvalFunctor functor ) const
218 {
219  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
220 
221  Dimension total = 0;
222 #ifdef DEBUG_VERBOSE
223  Dimension recount = 0;
224 #endif
225 
226  Quantity lastInnerSum;
227  Quantity lastOuterSum;
228 
229  Quantity innerSum, outerSum;
230  Quantity resultQuantity;
231 
232  Spel lastInnerSpel, lastOuterSpel;
233 
234  /// Iterate on all cells
235  for( SurfelIterator it = itbegin; it != itend; ++it )
236  {
237  if( total != 0 )
238  {
239 #ifdef DEBUG_VERBOSE
240  bool hasJumped = core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
241  recount = ( hasJumped ) ? recount + 1 : recount;
242 #else
243  core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
244 #endif
245  }
246  else
247  {
248  core_eval( it, innerSum, outerSum, false, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
249  }
250 
251  double lambda = 0.5;
252  resultQuantity = innerSum * lambda + outerSum * ( 1.0 - lambda );
253  *result++ = functor( resultQuantity );
254 
255  ++total;
256  }
257 
258 #ifdef DEBUG_VERBOSE
259  std::cout << "#total cells = " << total << std::endl;
260  std::cout << "#recount = " << recount << std::endl;
261 #endif
262 }
263 
264 
265 
266 
267 
268 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
269 template< typename SurfelIterator >
270 inline
271 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::CovarianceMatrix
272 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::evalCovarianceMatrix
273 ( const SurfelIterator & it ) const
274 {
275  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
276 
277  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
278 
279  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false );
280 
281  double lambda = 0.5;
282  return ( innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda ));
283 }
284 
285 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
286 template< typename SurfelIterator, typename EvalFunctor >
287 inline
288 typename EvalFunctor::Value
289 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::evalCovarianceMatrix
290 ( const SurfelIterator & it,
291  EvalFunctor functor ) const
292 {
293  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
294 
295  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
296  CovarianceMatrix resultCovarianceMatrix;
297 
298  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false );
299 
300  double lambda = 0.5;
301  resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
302  return functor( resultCovarianceMatrix );
303 }
304 
305 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
306 template< typename SurfelIterator, typename OutputIterator >
307 inline
308 void
309 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::evalCovarianceMatrix
310 ( const SurfelIterator & itbegin,
311  const SurfelIterator & itend,
312  OutputIterator & result ) const
313 {
314  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
315 
316  Dimension total = 0;
317 #ifdef DEBUG_VERBOSE
318  Dimension recount = 0;
319 #endif
320 
321  std::vector<Quantity> lastInnerMoments(nbMoments );
322  std::vector<Quantity> lastOuterMoments(nbMoments );
323 
324  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
325 
326  Spel lastInnerSpel, lastOuterSpel;
327 
328  /// Iterate on all cells
329  for( SurfelIterator it = itbegin; it != itend; ++it )
330  {
331  if( total != 0 )
332  {
333 #ifdef DEBUG_VERBOSE
334  bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
335  recount = ( hasJumped ) ? recount + 1 : recount;
336 #else
337  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
338 #endif
339  }
340  else
341  {
342  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
343  }
344 
345  double lambda = 0.5;
346  *result++ = ( innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda ));
347 
348  ++total;
349  }
350 
351 #ifdef DEBUG_VERBOSE
352  std::cout << "#total cells = " << total << std::endl;
353  std::cout << "#recount = " << recount << std::endl;
354 #endif
355 }
356 
357 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
358 template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
359 inline
360 void
361 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::evalCovarianceMatrix
362 ( const SurfelIterator & itbegin,
363  const SurfelIterator & itend,
364  OutputIterator & result,
365  EvalFunctor functor ) const
366 {
367  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
368 
369  Dimension total = 0;
370 #ifdef DEBUG_VERBOSE
371  Dimension recount = 0;
372 #endif
373 
374  std::vector<Quantity> lastInnerMoments( nbMoments );
375  std::vector<Quantity> lastOuterMoments( nbMoments );
376 
377  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
378  CovarianceMatrix resultCovarianceMatrix;
379 
380  Spel lastInnerSpel, lastOuterSpel;
381 
382  /// Iterate on all cells
383  for( SurfelIterator it = itbegin; it != itend; ++it )
384  {
385  if( total != 0 )
386  {
387 #ifdef DEBUG_VERBOSE
388  bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
389  recount = ( hasJumped ) ? recount + 1 : recount;
390 #else
391  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
392 #endif
393  }
394  else
395  {
396  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
397  }
398 
399  double lambda = 0.5;
400  resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
401  *result++ = functor( resultCovarianceMatrix );
402 
403  ++total;
404  }
405 
406 #ifdef DEBUG_VERBOSE
407  std::cout << "#total cells = " << total << std::endl;
408  std::cout << "#recount = " << recount << std::endl;
409 #endif
410 }
411 
412 
413 
414 
415 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
416 inline
417 bool
418 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::isValid() const
419 {
420  return true;
421 }
422 
423 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
424 void
425 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::fillMoments
426 ( Quantity * aMomentMatrix,
427  const Spel & aSpel,
428  double direction ) const
429 {
430  RealPoint current = myEmbedder( aSpel );
431  double x = current[ 0 ];
432  double y = current[ 1 ];
433 
434  aMomentMatrix[ 0 ] += direction * 1;
435  aMomentMatrix[ 1 ] += direction * y;
436  aMomentMatrix[ 2 ] += direction * x;
437  aMomentMatrix[ 3 ] += direction * x * y;
438  aMomentMatrix[ 4 ] += direction * y * y;
439  aMomentMatrix[ 5 ] += direction * x * x;
440 }
441 
442 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
443 void
444 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::computeCovarianceMatrix
445 ( const Quantity * aMomentMatrix,
446  CovarianceMatrix & aCovarianceMatrix ) const
447 {
448  /* Moment matrix:
449  *
450  * [ sum(1)
451  * sum(y) sum (x)
452  * sum(x*y) sum(y*y) sum(x*x)
453  * ]
454  */
455  MatrixQuantity A;
456  double B;
457  MatrixQuantity C;
458 
459  A.setComponent( 0, 0, aMomentMatrix[ 5 ] );
460  A.setComponent( 0, 1, aMomentMatrix[ 3 ] );
461  A.setComponent( 1, 0, aMomentMatrix[ 3 ] );
462  A.setComponent( 1, 1, aMomentMatrix[ 4 ] );
463 
464  B = 1.0 / aMomentMatrix[ 0 ];
465 
466  C.setComponent( 0, 0, aMomentMatrix[ 2 ] * aMomentMatrix[ 2 ] );
467  C.setComponent( 0, 1, aMomentMatrix[ 2 ] * aMomentMatrix[ 1 ] );
468  C.setComponent( 1, 0, aMomentMatrix[ 1 ] * aMomentMatrix[ 2 ] );
469  C.setComponent( 1, 1, aMomentMatrix[ 1 ] * aMomentMatrix[ 1 ] );
470 
471  aCovarianceMatrix = A - C * B;
472 }
473 
474 #ifndef _MSC_VER
475 // For Visual Studio, to be defined as a static const, it has to be intialized into the header file
476 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
477 const int
478 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::nbMoments = 6;
479 #endif //_MSC_VER
480 
481 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
482 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Spel
483 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultInnerSpel = Spel();
484 
485 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
486 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Spel
487 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultOuterSpel = Spel();
488 
489 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
490 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
491 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultInnerMoments[ 6 ] = {Quantity(0)};
492 
493 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
494 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
495 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultOuterMoments[ 6 ] = {Quantity(0)};
496 
497 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
498 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
499 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultInnerSum = Quantity(0);
500 
501 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
502 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
503 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultOuterSum = Quantity(0);
504 
505 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
506 template< typename SurfelIterator >
507 bool
508 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::core_eval
509 ( const SurfelIterator & it,
510  Quantity & innerSum,
511  Quantity & outerSum,
512  bool useLastResults,
513  Spel & lastInnerSpel,
514  Spel & lastOuterSpel,
515  Quantity & lastInnerSum,
516  Quantity & lastOuterSum ) const
517 {
518  boost::ignore_unused_variable_warning( it );
519  boost::ignore_unused_variable_warning( innerSum );
520  boost::ignore_unused_variable_warning( outerSum);
521  boost::ignore_unused_variable_warning(useLastResults);
522  boost::ignore_unused_variable_warning(lastInnerSum);
523  boost::ignore_unused_variable_warning(lastOuterSum);
524  boost::ignore_unused_variable_warning(lastInnerSpel);
525  boost::ignore_unused_variable_warning(lastOuterSpel);
526  trace.error() << "Unavailable yet." << std::endl;
527  return false;
528 }
529 
530 
531 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
532 template< typename SurfelIterator >
533 bool
534 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::core_evalCovarianceMatrix
535 ( const SurfelIterator & it,
536  CovarianceMatrix & innerMatrix,
537  CovarianceMatrix & outerMatrix,
538  bool useLastResults,
539  Spel & lastInnerSpel,
540  Spel & lastOuterSpel,
541  Quantity * lastInnerMoments,
542  Quantity * lastOuterMoments ) const
543 {
544  boost::ignore_unused_variable_warning(it);
545  boost::ignore_unused_variable_warning(innerMatrix);
546  boost::ignore_unused_variable_warning(outerMatrix);
547  boost::ignore_unused_variable_warning(useLastResults);
548  boost::ignore_unused_variable_warning(lastOuterMoments);
549  boost::ignore_unused_variable_warning(lastInnerMoments);
550  boost::ignore_unused_variable_warning(lastInnerSpel);
551  boost::ignore_unused_variable_warning(lastOuterSpel);
552  trace.error() << "Unavailable yet." << std::endl;
553  return false;
554 }
555 
556 
557 
558 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
559 ///////////////////////////////////////////////////// 2D /////////////////////////////////////////////////////
560 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
561 
562 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
563 inline
564 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >
565 ::DigitalSurfaceConvolver( ConstAlias< Functor > f,
566  ConstAlias< KernelFunctor > g,
567  ConstAlias< KSpace > space)
568  : dimension( 2 ),
569  myFFunctor( f ),
570  myGFunctor( g ),
571  myKSpace( space ),
572  isInitFullMasks( false ),
573  isInitKernelAndMasks( false )
574 {
575  myEmbedder = Embedder( myKSpace );
576 }
577 
578 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel>
579 inline
580 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >
581 ::DigitalSurfaceConvolver( const DigitalSurfaceConvolver& other )
582  : dimension( 2 ),
583  myFFunctor( other.myFFunctor ),
584  myGFunctor( other.myGFunctor ),
585  myKSpace( other.myKSpace ),
586  myEmbedder( other.myEmbedder ),
587  isInitFullMasks( other.isInitFullMasks ),
588  isInitKernelAndMasks( other.isInitKernelAndMasks ),
589  myMasks( other.myMasks ),
590  myKernel( other.myKernel ),
591  myKernelMask( other.myKernelMask ),
592  myKernelSpelOrigin( other.myKernelSpelOrigin )
593 {
594 }
595 
596 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
597 inline
598 void
599 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::init
600 ( const Point & pOrigin,
601  ConstAlias< PairIterators > fullKernel,
602  ConstAlias< std::vector< PairIterators > > masks )
603 {
604  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
605  myKernelMask = &fullKernel;
606  myMasks = &masks;
607 
608  ASSERT ( myMasks->size () == 9 );
609 
610  isInitFullMasks = true;
611  isInitKernelAndMasks = false;
612 }
613 
614 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
615 inline
616 void
617 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::init
618 ( const Point & pOrigin,
619  ConstAlias< DigitalKernel > fullKernel,
620  ConstAlias< std::vector< PairIterators > > masks )
621 {
622  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
623  myKernel = &fullKernel;
624  myMasks = &masks;
625 
626  ASSERT ( myMasks->size () == 9 );
627 
628  isInitFullMasks = false;
629  isInitKernelAndMasks = true;
630 }
631 
632 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
633 template< typename SurfelIterator >
634 inline
635 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
636 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::eval
637 ( const SurfelIterator & it ) const
638 {
639  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
640 
641  Quantity innerSum, outerSum;
642 
643  core_eval( it, innerSum, outerSum, false );
644 
645  double lambda = 0.5;
646  return ( innerSum * lambda + outerSum * ( 1.0 - lambda ));
647 }
648 
649 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
650 template< typename SurfelIterator, typename EvalFunctor >
651 inline
652 typename EvalFunctor::Value
653 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::eval
654 ( const SurfelIterator & it,
655  EvalFunctor functor ) const
656 {
657  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
658 
659  Quantity innerSum, outerSum;
660  Quantity resultQuantity;
661 
662  core_eval( it, innerSum, outerSum, false );
663 
664  double lambda = 0.5;
665  resultQuantity = innerSum * lambda + outerSum * ( 1.0 - lambda );
666  return functor( resultQuantity );
667 }
668 
669 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
670 template< typename SurfelIterator, typename OutputIterator >
671 inline
672 void
673 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::eval
674 ( const SurfelIterator & itbegin,
675  const SurfelIterator & itend,
676  OutputIterator & result ) const
677 {
678  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
679 
680  Dimension total = 0;
681 #ifdef DEBUG_VERBOSE
682  Dimension recount = 0;
683 #endif
684 
685  Quantity lastInnerSum;
686  Quantity lastOuterSum;
687 
688  Quantity innerSum, outerSum;
689 
690  Spel lastInnerSpel, lastOuterSpel;
691 
692  /// Iterate on all cells
693  for( SurfelIterator it = itbegin; it != itend; ++it )
694  {
695  if( total != 0 )
696  {
697 #ifdef DEBUG_VERBOSE
698  bool hasJumped = core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
699  recount = ( hasJumped ) ? recount + 1 : recount;
700 #else
701  core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
702 #endif
703  }
704  else
705  {
706  core_eval( it, innerSum, outerSum, false, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
707  }
708 
709  double lambda = 0.5;
710  *result++ = ( innerSum * lambda + outerSum * ( 1.0 - lambda ));
711 
712  ++total;
713  }
714 
715 #ifdef DEBUG_VERBOSE
716  std::cout << "#total cells = " << total << std::endl;
717  std::cout << "#recount = " << recount << std::endl;
718 #endif
719 }
720 
721 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
722 template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
723 inline
724 void
725 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::eval
726 ( const SurfelIterator & itbegin,
727  const SurfelIterator & itend,
728  OutputIterator & result,
729  EvalFunctor functor ) const
730 {
731  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
732 
733  Dimension total = 0;
734 #ifdef DEBUG_VERBOSE
735  Dimension recount = 0;
736 #endif
737 
738  Quantity lastInnerSum;
739  Quantity lastOuterSum;
740 
741  Quantity innerSum, outerSum;
742  Quantity resultQuantity;
743 
744  Spel lastInnerSpel, lastOuterSpel;
745 
746  /// Iterate on all cells
747  for( SurfelIterator it = itbegin; it != itend; ++it )
748  {
749  if( total != 0 )
750  {
751 #ifdef DEBUG_VERBOSE
752  bool hasJumped = core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
753  recount = ( hasJumped ) ? recount + 1 : recount;
754 #else
755  core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
756 #endif
757  }
758  else
759  {
760  core_eval( it, innerSum, outerSum, false, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
761  }
762 
763  double lambda = 0.5;
764  resultQuantity = innerSum * lambda + outerSum * ( 1.0 - lambda );
765  *result++ = functor( resultQuantity );
766 
767  ++total;
768  }
769 
770 #ifdef DEBUG_VERBOSE
771  std::cout << "#total cells = " << total << std::endl;
772  std::cout << "#recount = " << recount << std::endl;
773 #endif
774 }
775 
776 
777 
778 
779 
780 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
781 template< typename SurfelIterator >
782 inline
783 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::CovarianceMatrix
784 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::evalCovarianceMatrix
785 ( const SurfelIterator & it ) const
786 {
787  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
788 
789  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
790 
791  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false );
792 
793  double lambda = 0.5;
794  return ( innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda ));
795 }
796 
797 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
798 template< typename SurfelIterator, typename EvalFunctor >
799 inline
800 typename EvalFunctor::Value
801 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::evalCovarianceMatrix
802 ( const SurfelIterator & it,
803  EvalFunctor functor ) const
804 {
805  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
806 
807  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
808  CovarianceMatrix resultCovarianceMatrix;
809 
810  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false );
811 
812  double lambda = 0.5;
813  resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
814  return functor( resultCovarianceMatrix );
815 }
816 
817 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
818 template< typename SurfelIterator, typename OutputIterator >
819 inline
820 void
821 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::evalCovarianceMatrix
822 ( const SurfelIterator & itbegin,
823  const SurfelIterator & itend,
824  OutputIterator & result ) const
825 {
826  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
827 
828  Dimension total = 0;
829 #ifdef DEBUG_VERBOSE
830  Dimension recount = 0;
831 #endif
832 
833  std::vector<Quantity> lastInnerMoments( nbMoments );
834  std::vector<Quantity> lastOuterMoments( nbMoments );
835 
836  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
837 
838  Spel lastInnerSpel, lastOuterSpel;
839 
840  /// Iterate on all cells
841  for( SurfelIterator it = itbegin; it != itend; ++it )
842  {
843  if( total != 0 )
844  {
845 #ifdef DEBUG_VERBOSE
846  bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
847  recount = ( hasJumped ) ? recount + 1 : recount;
848 #else
849  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
850 #endif
851  }
852  else
853  {
854  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
855  }
856 
857  double lambda = 0.5;
858  *result++ = ( innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda ));
859 
860  ++total;
861  }
862 
863 #ifdef DEBUG_VERBOSE
864  std::cout << "#total cells = " << total << std::endl;
865  std::cout << "#recount = " << recount << std::endl;
866 #endif
867 }
868 
869 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
870 template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
871 inline
872 void
873 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::evalCovarianceMatrix
874 ( const SurfelIterator & itbegin,
875  const SurfelIterator & itend,
876  OutputIterator & result,
877  EvalFunctor functor ) const
878 {
879  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
880 
881  Dimension total = 0;
882 #ifdef DEBUG_VERBOSE
883  Dimension recount = 0;
884 #endif
885 
886  std::vector<Quantity> lastInnerMoments( nbMoments );
887  std::vector<Quantity> lastOuterMoments( nbMoments );
888 
889  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
890  CovarianceMatrix resultCovarianceMatrix;
891 
892  Spel lastInnerSpel, lastOuterSpel;
893 
894  /// Iterate on all cells
895  for( SurfelIterator it = itbegin; it != itend; ++it )
896  {
897  if( total != 0 )
898  {
899 #ifdef DEBUG_VERBOSE
900  bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
901  recount = ( hasJumped ) ? recount + 1 : recount;
902 #else
903  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
904 #endif
905  }
906  else
907  {
908  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
909  }
910 
911  double lambda = 0.5;
912  resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
913  *result++ = functor( resultCovarianceMatrix );
914 
915  ++total;
916  }
917 
918 #ifdef DEBUG_VERBOSE
919  std::cout << "#total cells = " << total << std::endl;
920  std::cout << "#recount = " << recount << std::endl;
921 #endif
922 }
923 
924 
925 
926 
927 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
928 inline
929 bool
930 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::isValid() const
931 {
932  return true;
933 }
934 
935 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
936 void
937 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::fillMoments
938 ( Quantity * aMomentMatrix,
939  const Spel & aSpel,
940  double direction ) const
941 {
942  RealPoint current = myEmbedder( aSpel );
943  double x = current[ 0 ];
944  double y = current[ 1 ];
945 
946  aMomentMatrix[ 0 ] += direction * 1;
947  aMomentMatrix[ 1 ] += direction * y;
948  aMomentMatrix[ 2 ] += direction * x;
949  aMomentMatrix[ 3 ] += direction * x * y;
950  aMomentMatrix[ 4 ] += direction * y * y;
951  aMomentMatrix[ 5 ] += direction * x * x;
952 }
953 
954 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
955 void
956 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::computeCovarianceMatrix
957 ( const Quantity * aMomentMatrix,
958  CovarianceMatrix & aCovarianceMatrix ) const
959 {
960  MatrixQuantity A;
961  double B;
962  MatrixQuantity C;
963 
964  A.setComponent( 0, 0, aMomentMatrix[ 5 ] );
965  A.setComponent( 0, 1, aMomentMatrix[ 3 ] );
966  A.setComponent( 1, 0, aMomentMatrix[ 3 ] );
967  A.setComponent( 1, 1, aMomentMatrix[ 4 ] );
968 
969  B = 1.0 / aMomentMatrix[ 0 ];
970 
971  C.setComponent( 0, 0, aMomentMatrix[ 2 ] * aMomentMatrix[ 2 ] );
972  C.setComponent( 0, 1, aMomentMatrix[ 2 ] * aMomentMatrix[ 1 ] );
973  C.setComponent( 1, 0, aMomentMatrix[ 1 ] * aMomentMatrix[ 2 ] );
974  C.setComponent( 1, 1, aMomentMatrix[ 1 ] * aMomentMatrix[ 1 ] );
975 
976  aCovarianceMatrix = A - C * B;
977 }
978 
979 #ifndef _MSC_VER
980 // For Visual Studio, to be defined as a static const, it has to be intialized into the header file
981 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
982 const int
983 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::nbMoments = 6;
984 #endif //_MSC_VER
985 
986 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
987 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Spel
988 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultInnerSpel = Spel();
989 
990 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
991 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Spel
992 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultOuterSpel = Spel();
993 
994 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
995 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
996 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultInnerMoments[ 6 ] = {Quantity(0)};
997 
998 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
999 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
1000 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultOuterMoments[ 6 ] = {Quantity(0)};
1001 
1002 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1003 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
1004 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultInnerSum = Quantity(0);
1005 
1006 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1007 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
1008 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultOuterSum = Quantity(0);
1009 
1010 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1011 template< typename SurfelIterator >
1012 bool
1013 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::core_eval
1014 ( const SurfelIterator & it,
1015  Quantity & innerSum,
1016  Quantity & outerSum,
1017  bool useLastResults,
1018  Spel & lastInnerSpel,
1019  Spel & lastOuterSpel,
1020  Quantity & lastInnerSum,
1021  Quantity & lastOuterSum ) const
1022 {
1023  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1024 
1025  using KPS = typename KSpace::PreCellularGridSpace;
1026 
1027 #ifdef DEBUG_VERBOSE
1028  Dimension recount = false;
1029 #endif
1030 
1031  typedef typename Functor::Quantity FQuantity;
1032  DGtal::Dimension nbMasks =static_cast<DGtal::Dimension>( (long)myMasks->size() - 1);
1033  DGtal::Dimension positionOfFullKernel = 4;
1034 
1035  Quantity m = NumberTraits< Quantity >::ZERO;
1036 
1037  Spel currentInnerSpel, currentOuterSpel;
1038  Spel shiftedSpel;
1039  Point shiftInnerSpel, shiftOuterSpel;
1040  Point diffSpel;
1041 
1042  bool bComputed = false; /// <! if the cell has already been computed, continue to the next
1043 
1044  int x, y, x2, y2, x2y2;
1045  unsigned int offset;
1046 
1047  /// Inner cell
1048  {
1049  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
1050  currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
1051  shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
1052 
1053  /// Check if we can use previous results
1054  if( useLastResults )
1055  {
1056  bComputed = false;
1057 
1058  if( currentInnerSpel == lastInnerSpel )
1059  {
1060  m = lastInnerSum;
1061  innerSum = m;
1062 
1063  bComputed = true;
1064  }
1065  else if( currentInnerSpel == lastOuterSpel )
1066  {
1067  m = lastOuterSum;
1068  outerSum = m;
1069 
1070  bComputed = true;
1071  }
1072  else
1073  {
1074  diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
1075 
1076  x = diffSpel[ 0 ];
1077  y = diffSpel[ 1 ];
1078  x2 = x * x;
1079  y2 = y * y;
1080  x2y2 = x2 + y2;
1081 
1082  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 );
1083 
1084  if( x2y2 != 4 && x2y2 != 8 ) /// Previous and current cells aren't adjacent. Compute on the full kernel
1085  {
1086  useLastResults = false;
1087 #ifdef DEBUG_VERBOSE
1088  recount = true;
1089 #endif
1090  }
1091  else if( offset == positionOfFullKernel ) /// Full kernel in 2D. Never reached because case already considered before. It's a honeypot.
1092  {
1093  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
1094  }
1095  else
1096  {
1097  useLastResults = true;
1098  }
1099  }
1100  }
1101 
1102  if( !bComputed )
1103  {
1104  if( !useLastResults ) /// Computation on full kernel, we have no previous results
1105  {
1106  m = NumberTraits< Quantity >::ZERO;
1107 
1108  if( isInitFullMasks )
1109  {
1110  for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
1111  {
1112  auto preShiftedSpel = KPS::sSpel( *itm );
1113  preShiftedSpel.coordinates += shiftInnerSpel;
1114 
1115  if( myKSpace.sIsInside( preShiftedSpel ) )
1116  {
1117  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1118 
1119  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1120  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1121 
1122  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1123  {
1124  m += 1.0;
1125  }
1126  }
1127  }
1128  }
1129  else if( isInitKernelAndMasks )
1130  {
1131  Domain domain = myKernel->getDomain();
1132  for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
1133  {
1134  if( myKernel->operator()( *itm ))
1135  {
1136  auto preShiftedSpel = KPS::sSpel( *itm );
1137  preShiftedSpel.coordinates += shiftInnerSpel;
1138 
1139  if( myKSpace.sIsInside( preShiftedSpel ) )
1140  {
1141  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1142 
1143  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1144  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1145 
1146  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1147  {
1148  m += 1.0;
1149  }
1150  }
1151  }
1152  }
1153  }
1154  else
1155  {
1156  trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
1157  return false;
1158  }
1159 
1160  }
1161  else /// Using lastInnerMoments
1162  {
1163  m = lastInnerSum;
1164 
1165  /// Part to substract from previous result.
1166  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
1167  {
1168  auto preShiftedSpel = KPS::sSpel( *itm );
1169  preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
1170 
1171  if( myKSpace.sIsInside( preShiftedSpel ) )
1172  {
1173  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1174 
1175  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1176  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1177 
1178  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1179  {
1180  m -= 1.0;
1181  }
1182  }
1183  }
1184 
1185  /// Part to add from previous result.
1186  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
1187  {
1188  auto preShiftedSpel = KPS::sSpel( *itm );
1189  preShiftedSpel.coordinates += shiftInnerSpel;
1190 
1191  if( myKSpace.sIsInside( preShiftedSpel ) )
1192  {
1193  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1194 
1195  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1196  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1197 
1198  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1199  {
1200  m += 1.0;
1201  }
1202  }
1203  }
1204  }
1205 
1206  /// Computation of covariance Matrix
1207  innerSum = m;
1208  lastInnerSum = m;
1209  }
1210  }
1211 
1212  /// Outter cell
1213  {
1214  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
1215  currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
1216  shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
1217 
1218  ASSERT( currentInnerSpel != currentOuterSpel );
1219 
1220  diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
1221 
1222  x = diffSpel[ 0 ];
1223  y = diffSpel[ 1 ];
1224  x2 = x * x;
1225  y2 = y * y;
1226  x2y2 = x2 + y2;
1227 
1228  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 );
1229 
1230  if( x2y2 != 4 && x2y2 != 8 ) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
1231  {
1232  trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
1233  }
1234  else if( offset == positionOfFullKernel ) /// Full kernel in 2D. Never reached because case already considered before. It's a honeypot.
1235  {
1236  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
1237  }
1238  else
1239  {
1240  /// Part to substract from previous result.
1241  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
1242  {
1243  auto preShiftedSpel = KPS::sSpel( *itm );
1244  preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
1245 
1246  if( myKSpace.sIsInside( preShiftedSpel ) )
1247  {
1248  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1249 
1250  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1251  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1252 
1253  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1254  {
1255  m -= 1.0;
1256  }
1257  }
1258  }
1259 
1260  /// Part to add from previous result.
1261  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
1262  {
1263  auto preShiftedSpel = KPS::sSpel( *itm );
1264  preShiftedSpel.coordinates += shiftOuterSpel;
1265 
1266  if( myKSpace.sIsInside( preShiftedSpel ) )
1267  {
1268  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1269 
1270  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1271  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1272 
1273  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1274  {
1275  m += 1.0;
1276  }
1277  }
1278  }
1279  }
1280 
1281  /// Computation of covariance Matrix
1282  outerSum = m;
1283  lastOuterSum = m;
1284  }
1285 
1286  ASSERT (( lastInnerSum != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
1287 
1288  lastInnerSpel = currentInnerSpel;
1289  lastOuterSpel = currentOuterSpel;
1290 
1291 #ifdef DEBUG_VERBOSE
1292  return recount;
1293 #else
1294  return false;
1295 #endif
1296 }
1297 
1298 
1299 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1300 template< typename SurfelIterator >
1301 bool
1302 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::core_evalCovarianceMatrix
1303 ( const SurfelIterator & it,
1304  CovarianceMatrix & innerMatrix,
1305  CovarianceMatrix & outerMatrix,
1306  bool useLastResults,
1307  Spel & lastInnerSpel,
1308  Spel & lastOuterSpel,
1309  Quantity * lastInnerMoments,
1310  Quantity * lastOuterMoments ) const
1311 {
1312  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1313 
1314  using KPS = typename KSpace::PreCellularGridSpace;
1315 
1316 #ifdef DEBUG_VERBOSE
1317  Dimension recount = false;
1318 #endif
1319 
1320  typedef typename Functor::Quantity FQuantity;
1321  DGtal::Dimension nbMasks = myMasks->size() - 1;
1322  DGtal::Dimension positionOfFullKernel = 4;
1323 
1324  Quantity m[ nbMoments ]; /// <! [ m00, m01, m10, m11, m02, m20 ]
1325  std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
1326 
1327  Spel currentInnerSpel, currentOuterSpel;
1328  Spel shiftedSpel;
1329  Point shiftInnerSpel, shiftOuterSpel;
1330  Point diffSpel;
1331 
1332  bool bComputed = false; /// <! if the cell has already been computed, continue to the next
1333 
1334  int x, y, x2, y2, x2y2;
1335  unsigned int offset;
1336 
1337  /// Inner cell
1338  {
1339  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
1340  currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
1341  shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
1342 
1343  /// Check if we can use previous results
1344  if( useLastResults )
1345  {
1346  bComputed = false;
1347 
1348  if( currentInnerSpel == lastInnerSpel )
1349  {
1350  memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
1351  computeCovarianceMatrix( m, innerMatrix );
1352 
1353  bComputed = true;
1354  }
1355  else if( currentInnerSpel == lastOuterSpel )
1356  {
1357  memcpy( m, lastOuterMoments, nbMoments * sizeof( Quantity ));
1358  computeCovarianceMatrix( m, outerMatrix );
1359 
1360  bComputed = true;
1361  }
1362  else
1363  {
1364  diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
1365 
1366  x = diffSpel[ 0 ];
1367  y = diffSpel[ 1 ];
1368  x2 = x * x;
1369  y2 = y * y;
1370  x2y2 = x2 + y2;
1371 
1372  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 );
1373 
1374  if( x2y2 != 4 && x2y2 != 8 ) /// Previous and current cells aren't adjacent. Compute on the full kernel
1375  {
1376  useLastResults = false;
1377 #ifdef DEBUG_VERBOSE
1378  recount = true;
1379 #endif
1380  }
1381  else if( offset == positionOfFullKernel ) /// Full kernel in 2D. Never reached because case already considered before. It's a honeypot.
1382  {
1383  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
1384  }
1385  else
1386  {
1387  useLastResults = true;
1388  }
1389  }
1390  }
1391 
1392  if( !bComputed )
1393  {
1394  if( !useLastResults ) /// Computation on full kernel, we have no previous results
1395  {
1396  std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
1397 
1398  if( isInitFullMasks )
1399  {
1400  for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
1401  {
1402  auto preShiftedSpel = KPS::sSpel( *itm );
1403  preShiftedSpel.coordinates += shiftInnerSpel;
1404 
1405  if( myKSpace.sIsInside( preShiftedSpel ) )
1406  {
1407  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1408 
1409  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1410  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1411 
1412  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1413  {
1414  fillMoments( m, shiftedSpel, 1.0 );
1415  }
1416  }
1417  }
1418  }
1419  else if( isInitKernelAndMasks )
1420  {
1421  Domain domain = myKernel->getDomain();
1422  for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
1423  {
1424  if( myKernel->operator()( *itm ))
1425  {
1426  auto preShiftedSpel = KPS::sSpel( *itm );
1427  preShiftedSpel.coordinates += shiftInnerSpel;
1428 
1429  if( myKSpace.sIsInside( preShiftedSpel ) )
1430  {
1431  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1432 
1433  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1434  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1435 
1436  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1437  {
1438  fillMoments( m, shiftedSpel, 1.0 );
1439  }
1440  }
1441  }
1442  }
1443  }
1444  else
1445  {
1446  trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
1447  return false;
1448  }
1449  }
1450  else /// Using lastInnerMoments
1451  {
1452  memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
1453 
1454  /// Part to substract from previous result.
1455  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
1456  {
1457  auto preShiftedSpel = KPS::sSpel( *itm );
1458  preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
1459 
1460  if( myKSpace.sIsInside( preShiftedSpel ) )
1461  {
1462  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1463 
1464  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1465  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1466 
1467  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1468  {
1469  fillMoments( m, shiftedSpel, -1.0 );
1470  }
1471  }
1472  }
1473 
1474  /// Part to add from previous result.
1475  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
1476  {
1477  auto preShiftedSpel = KPS::sSpel( *itm );
1478  preShiftedSpel.coordinates += shiftInnerSpel;
1479 
1480  if( myKSpace.sIsInside( preShiftedSpel ) )
1481  {
1482  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1483 
1484  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1485  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1486 
1487  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1488  {
1489  fillMoments( m, shiftedSpel, 1.0 );
1490  }
1491  }
1492  }
1493  }
1494 
1495  /// Computation of covariance Matrix
1496  computeCovarianceMatrix( m, innerMatrix );
1497  memcpy( lastInnerMoments, m, nbMoments * sizeof( Quantity ));
1498  }
1499  }
1500 
1501  /// Outter cell
1502  {
1503  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
1504  currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
1505  shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
1506 
1507  ASSERT( currentInnerSpel != currentOuterSpel );
1508 
1509  diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
1510 
1511  x = diffSpel[ 0 ];
1512  y = diffSpel[ 1 ];
1513  x2 = x * x;
1514  y2 = y * y;
1515  x2y2 = x2 + y2;
1516 
1517  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 );
1518 
1519  if( x2y2 != 4 && x2y2 != 8 ) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
1520  {
1521  trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
1522  }
1523  else if( offset == positionOfFullKernel ) /// Full kernel in 2D. Never reached because case already considered before. It's a honeypot.
1524  {
1525  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
1526  }
1527  else
1528  {
1529  /// Part to substract from previous result.
1530  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
1531  {
1532  auto preShiftedSpel = KPS::sSpel( *itm );
1533  preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
1534  if( myKSpace.sIsInside( preShiftedSpel ) )
1535  {
1536  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1537 
1538  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1539  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1540 
1541  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1542  {
1543  fillMoments( m, shiftedSpel, -1.0 );
1544  }
1545  }
1546  }
1547 
1548  /// Part to add from previous result.
1549  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
1550  {
1551  auto preShiftedSpel = KPS::sSpel( *itm );
1552  preShiftedSpel.coordinates += shiftOuterSpel;
1553 
1554  if( myKSpace.sIsInside( preShiftedSpel ) )
1555  {
1556  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1557 
1558  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1559  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1560 
1561  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1562  {
1563  fillMoments( m, shiftedSpel, 1.0 );
1564  }
1565  }
1566  }
1567  }
1568 
1569  /// Computation of covariance Matrix
1570  computeCovarianceMatrix( m, outerMatrix );
1571  memcpy( lastOuterMoments, m, nbMoments * sizeof( Quantity ));
1572  }
1573 
1574  ASSERT (( lastInnerMoments[ 0 ] != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
1575 
1576  lastInnerSpel = currentInnerSpel;
1577  lastOuterSpel = currentOuterSpel;
1578 
1579 #ifdef DEBUG_VERBOSE
1580  return recount;
1581 #else
1582  return false;
1583 #endif
1584 }
1585 
1586 
1587 
1588 
1589 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
1590 ///////////////////////////////////////////////////// 3D /////////////////////////////////////////////////////
1591 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
1592 
1593 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1594 inline
1595 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >
1596 ::DigitalSurfaceConvolver( ConstAlias< Functor > f,
1597  ConstAlias< KernelFunctor > g,
1598  ConstAlias< KSpace > space)
1599  : dimension( 3 ),
1600  myFFunctor( f ),
1601  myGFunctor( g ),
1602  myKSpace( space ),
1603  isInitFullMasks( false ),
1604  isInitKernelAndMasks( false )
1605 {
1606  myEmbedder = Embedder( myKSpace );
1607 }
1608 
1609 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel>
1610 inline
1611 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >
1612 ::DigitalSurfaceConvolver( const DigitalSurfaceConvolver& other )
1613  : dimension( 3 ),
1614  myFFunctor( other.myFFunctor ),
1615  myGFunctor( other.myGFunctor ),
1616  myKSpace( other.myKSpace ),
1617  myEmbedder( other.myEmbedder ),
1618  isInitFullMasks( other.isInitFullMasks ),
1619  isInitKernelAndMasks( other.isInitKernelAndMasks ),
1620  myMasks( other.myMasks ),
1621  myKernel( other.myKernel ),
1622  myKernelMask( other.myKernelMask ),
1623  myKernelSpelOrigin( other.myKernelSpelOrigin )
1624 {
1625 }
1626 
1627 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1628 inline
1629 void
1630 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::init
1631 ( const Point & pOrigin,
1632  ConstAlias< PairIterators > fullKernel,
1633  ConstAlias< std::vector< PairIterators > > masks )
1634 {
1635  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
1636  myKernelMask = &fullKernel;
1637  myMasks = &masks;
1638 
1639  ASSERT ( myMasks->size () == 27 );
1640 
1641  isInitFullMasks = true;
1642  isInitKernelAndMasks = false;
1643 }
1644 
1645 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1646 inline
1647 void
1648 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::init
1649 ( const Point & pOrigin,
1650  ConstAlias< DigitalKernel > fullKernel,
1651  ConstAlias< std::vector< PairIterators > > masks )
1652 {
1653  myKernelSpelOrigin = myKSpace.sSpel( pOrigin );
1654  myKernel = &fullKernel;
1655  myMasks = &masks;
1656 
1657  ASSERT ( myMasks->size () == 27 );
1658 
1659  isInitFullMasks = false;
1660  isInitKernelAndMasks = true;
1661 }
1662 
1663 
1664 
1665 
1666 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1667 template< typename SurfelIterator >
1668 inline
1669 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
1670 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::eval
1671 ( const SurfelIterator & it ) const
1672 {
1673  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1674 
1675  Quantity innerSum, outerSum;
1676 
1677  core_eval( it, innerSum, outerSum, false );
1678 
1679  double lambda = 0.5;
1680  return ( innerSum * lambda + outerSum * ( 1.0 - lambda ));
1681 }
1682 
1683 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1684 template< typename SurfelIterator, typename EvalFunctor >
1685 inline
1686 typename EvalFunctor::Value
1687 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::eval
1688 ( const SurfelIterator & it,
1689  EvalFunctor functor ) const
1690 {
1691  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1692 
1693  Quantity innerSum, outerSum;
1694  Quantity resultQuantity;
1695 
1696  core_eval( it, innerSum, outerSum, false );
1697 
1698  double lambda = 0.5;
1699  resultQuantity = innerSum * lambda + outerSum * ( 1.0 - lambda );
1700  return functor( resultQuantity );
1701 }
1702 
1703 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1704 template< typename SurfelIterator, typename OutputIterator >
1705 inline
1706 void
1707 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::eval
1708 ( const SurfelIterator & itbegin,
1709  const SurfelIterator & itend,
1710  OutputIterator & result ) const
1711 {
1712  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1713 
1714  Dimension total = 0;
1715 #ifdef DEBUG_VERBOSE
1716  Dimension recount = 0;
1717 #endif
1718 
1719  Quantity lastInnerSum;
1720  Quantity lastOuterSum;
1721 
1722  Quantity innerSum, outerSum;
1723 
1724  Spel lastInnerSpel, lastOuterSpel;
1725 
1726  /// Iterate on all cells
1727  for( SurfelIterator it = itbegin; it != itend; ++it )
1728  {
1729  if( total != 0 )
1730  {
1731 #ifdef DEBUG_VERBOSE
1732  bool hasJumped = core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
1733  recount = ( hasJumped ) ? recount + 1 : recount;
1734 #else
1735  core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
1736 #endif
1737  }
1738  else
1739  {
1740  core_eval( it, innerSum, outerSum, false, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
1741  }
1742 
1743  double lambda = 0.5;
1744  *result++ = ( innerSum * lambda + outerSum * ( 1.0 - lambda ));
1745 
1746  ++total;
1747  }
1748 
1749 #ifdef DEBUG_VERBOSE
1750  std::cout << "#total cells = " << total << std::endl;
1751  std::cout << "#recount = " << recount << std::endl;
1752 #endif
1753 }
1754 
1755 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1756 template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
1757 inline
1758 void
1759 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::eval
1760 ( const SurfelIterator & itbegin,
1761  const SurfelIterator & itend,
1762  OutputIterator & result,
1763  EvalFunctor functor ) const
1764 {
1765  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1766 
1767  Dimension total = 0;
1768 #ifdef DEBUG_VERBOSE
1769  Dimension recount = 0;
1770 #endif
1771 
1772  Quantity lastInnerSum;
1773  Quantity lastOuterSum;
1774 
1775  Quantity innerSum, outerSum;
1776  Quantity resultQuantity;
1777 
1778  Spel lastInnerSpel, lastOuterSpel;
1779 
1780  /// Iterate on all cells
1781  for( SurfelIterator it = itbegin; it != itend; ++it )
1782  {
1783  if( total != 0 )
1784  {
1785 #ifdef DEBUG_VERBOSE
1786  bool hasJumped = core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
1787  recount = ( hasJumped ) ? recount + 1 : recount;
1788 #else
1789  core_eval( it, innerSum, outerSum, true, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
1790 #endif
1791  }
1792  else
1793  {
1794  core_eval( it, innerSum, outerSum, false, lastInnerSpel, lastOuterSpel, lastInnerSum, lastOuterSum );
1795  }
1796 
1797  double lambda = 0.5;
1798  resultQuantity = innerSum * lambda + outerSum * ( 1.0 - lambda );
1799  *result++ = functor( resultQuantity );
1800 
1801  ++total;
1802  }
1803 
1804 #ifdef DEBUG_VERBOSE
1805  std::cout << "#total cells = " << total << std::endl;
1806  std::cout << "#recount = " << recount << std::endl;
1807 #endif
1808 }
1809 
1810 
1811 
1812 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1813 template< typename SurfelIterator >
1814 inline
1815 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::CovarianceMatrix
1816 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::evalCovarianceMatrix
1817 ( const SurfelIterator & it ) const
1818 {
1819  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1820 
1821  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
1822 
1823  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false );
1824 
1825  double lambda = 0.5;
1826  return ( innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda ));
1827 }
1828 
1829 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1830 template< typename SurfelIterator, typename EvalFunctor >
1831 inline
1832 typename EvalFunctor::Value
1833 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::evalCovarianceMatrix
1834 ( const SurfelIterator & it,
1835  EvalFunctor functor ) const
1836 {
1837  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1838 
1839  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
1840  CovarianceMatrix resultCovarianceMatrix;
1841 
1842  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false );
1843 
1844  double lambda = 0.5;
1845  resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
1846  return functor( resultCovarianceMatrix );
1847 }
1848 
1849 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1850 template< typename SurfelIterator, typename OutputIterator >
1851 inline
1852 void
1853 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::evalCovarianceMatrix
1854 ( const SurfelIterator & itbegin,
1855  const SurfelIterator & itend,
1856  OutputIterator & result ) const
1857 {
1858  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1859 
1860  Dimension total = 0;
1861 #ifdef DEBUG_VERBOSE
1862  Dimension recount = 0;
1863 #endif
1864 
1865  std::vector<Quantity> lastInnerMoments( nbMoments );
1866  std::vector<Quantity> lastOuterMoments( nbMoments );
1867 
1868  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
1869 
1870  Spel lastInnerSpel, lastOuterSpel;
1871 
1872  /// Iterate on all cells
1873  for( SurfelIterator it = itbegin; it != itend; ++it )
1874  {
1875  if( total != 0 )
1876  {
1877 #ifdef DEBUG_VERBOSE
1878  bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1879  recount = ( hasJumped ) ? recount + 1 : recount;
1880 #else
1881  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1882 #endif
1883  }
1884  else
1885  {
1886  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1887  }
1888 
1889  double lambda = 0.5;
1890  *result++ = ( innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda ));
1891 
1892  ++total;
1893  }
1894 
1895 #ifdef DEBUG_VERBOSE
1896  std::cout << "#total cells = " << total << std::endl;
1897  std::cout << "#recount = " << recount << std::endl;
1898 #endif
1899 }
1900 
1901 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1902 template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
1903 inline
1904 void
1905 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::evalCovarianceMatrix
1906 ( const SurfelIterator & itbegin,
1907  const SurfelIterator & itend,
1908  OutputIterator & result,
1909  EvalFunctor functor ) const
1910 {
1911  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1912 
1913  Dimension total = 0;
1914 #ifdef DEBUG_VERBOSE
1915  Dimension recount = 0;
1916 #endif
1917 
1918  Quantity *lastInnerMoments = new Quantity[ nbMoments ];
1919  Quantity *lastOuterMoments = new Quantity[ nbMoments ];
1920 
1921  CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
1922  CovarianceMatrix resultCovarianceMatrix;
1923 
1924  Spel lastInnerSpel, lastOuterSpel;
1925 
1926  /// Iterate on all cells
1927  for( SurfelIterator it = itbegin; it != itend; ++it )
1928  {
1929  if( total != 0 )
1930  {
1931 #ifdef DEBUG_VERBOSE
1932  bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1933  recount = ( hasJumped ) ? recount + 1 : recount;
1934 #else
1935  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1936 #endif
1937  }
1938  else
1939  {
1940  core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1941  }
1942 
1943  double lambda = 0.5;
1944  resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
1945  *result++ = functor( resultCovarianceMatrix );
1946 
1947  ++total;
1948  }
1949 
1950 #ifdef DEBUG_VERBOSE
1951  std::cout << "#total cells = " << total << std::endl;
1952  std::cout << "#recount = " << recount << std::endl;
1953 #endif
1954 
1955  delete[] lastOuterMoments;
1956  delete[] lastInnerMoments;
1957 }
1958 
1959 
1960 
1961 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1962 inline
1963 bool
1964 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::isValid() const
1965 {
1966  return true;
1967 }
1968 
1969 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1970 void
1971 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::fillMoments
1972 ( Quantity * aMomentMatrix,
1973  const Spel & aSpel,
1974  double direction ) const
1975 {
1976  RealPoint current = myEmbedder( aSpel );
1977  double x = current[ 0 ];
1978  double y = current[ 1 ];
1979  double z = current[ 2 ];
1980 
1981  aMomentMatrix[ 0 ] += direction * 1;
1982  aMomentMatrix[ 1 ] += direction * z;
1983  aMomentMatrix[ 2 ] += direction * y;
1984  aMomentMatrix[ 3 ] += direction * x;
1985  aMomentMatrix[ 4 ] += direction * y * z;
1986  aMomentMatrix[ 5 ] += direction * x * z;
1987  aMomentMatrix[ 6 ] += direction * x * y;
1988  aMomentMatrix[ 7 ] += direction * z * z;
1989  aMomentMatrix[ 8 ] += direction * y * y;
1990  aMomentMatrix[ 9 ] += direction * x * x;
1991 
1992 
1993 }
1994 
1995 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1996 void
1997 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::computeCovarianceMatrix
1998 ( const Quantity * aMomentMatrix,
1999  CovarianceMatrix & aCovarianceMatrix ) const
2000 {
2001  MatrixQuantity A;
2002  double B;
2003  MatrixQuantity C;
2004 
2005  A.setComponent( 0, 0, aMomentMatrix[ 9 ] );
2006  A.setComponent( 0, 1, aMomentMatrix[ 6 ] );
2007  A.setComponent( 0, 2, aMomentMatrix[ 5 ] );
2008  A.setComponent( 1, 0, aMomentMatrix[ 6 ] );
2009  A.setComponent( 1, 1, aMomentMatrix[ 8 ] );
2010  A.setComponent( 1, 2, aMomentMatrix[ 4 ] );
2011  A.setComponent( 2, 0, aMomentMatrix[ 5 ] );
2012  A.setComponent( 2, 1, aMomentMatrix[ 4 ] );
2013  A.setComponent( 2, 2, aMomentMatrix[ 7 ] );
2014 
2015  B = 1.0 / aMomentMatrix[ 0 ];
2016 
2017  C.setComponent( 0, 0, aMomentMatrix[ 3 ] * aMomentMatrix[ 3 ] );
2018  C.setComponent( 0, 1, aMomentMatrix[ 3 ] * aMomentMatrix[ 2 ] );
2019  C.setComponent( 0, 2, aMomentMatrix[ 3 ] * aMomentMatrix[ 1 ] );
2020  C.setComponent( 1, 0, aMomentMatrix[ 2 ] * aMomentMatrix[ 3 ] );
2021  C.setComponent( 1, 1, aMomentMatrix[ 2 ] * aMomentMatrix[ 2 ] );
2022  C.setComponent( 1, 2, aMomentMatrix[ 2 ] * aMomentMatrix[ 1 ] );
2023  C.setComponent( 2, 0, aMomentMatrix[ 1 ] * aMomentMatrix[ 3 ] );
2024  C.setComponent( 2, 1, aMomentMatrix[ 1 ] * aMomentMatrix[ 2 ] );
2025  C.setComponent( 2, 2, aMomentMatrix[ 1 ] * aMomentMatrix[ 1 ] );
2026 
2027  aCovarianceMatrix = A - C * B;
2028 }
2029 
2030 #ifndef _MSC_VER
2031 // For Visual Studio, to be defined as a static const, it has to be intialized into the header file
2032 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2033 const int
2034 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::nbMoments = 10;
2035 #endif //_MSC_VER
2036 
2037 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2038 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Spel
2039 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerSpel = Spel();
2040 
2041 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2042 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Spel
2043 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterSpel = Spel();
2044 
2045 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2046 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2047 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerMoments[ 10 ] = {Quantity(0)};
2048 
2049 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2050 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2051 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterMoments[ 10 ] = {Quantity(0)};
2052 
2053 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2054 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2055 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerSum = Quantity(0);
2056 
2057 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2058 typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2059 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterSum = Quantity(0);
2060 
2061 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2062 template< typename SurfelIterator >
2063 bool
2064 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_eval
2065 ( const SurfelIterator & it,
2066  Quantity & innerSum,
2067  Quantity & outerSum,
2068  bool useLastResults,
2069  Spel & lastInnerSpel,
2070  Spel & lastOuterSpel,
2071  Quantity & lastInnerSum,
2072  Quantity & lastOuterSum ) const
2073 {
2074  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2075 
2076  using KPS = typename KSpace::PreCellularGridSpace;
2077 
2078 #ifdef DEBUG_VERBOSE
2079  Dimension recount = false;
2080 #endif
2081 
2082  typedef typename Functor::Quantity FQuantity;
2083  DGtal::Dimension nbMasks = static_cast<DGtal::Dimension>(myMasks->size()) - 1;
2084  DGtal::Dimension positionOfFullKernel = 13;
2085 
2086  Quantity m = NumberTraits< Quantity >::ZERO;
2087 
2088  Spel currentInnerSpel, currentOuterSpel;
2089  Spel shiftedSpel;
2090  Point shiftInnerSpel, shiftOuterSpel;
2091  Point diffSpel;
2092 
2093  bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2094 
2095  int x, y, z, x2, y2, z2, x2y2z2;
2096  unsigned int offset;
2097 
2098  /// Inner cell
2099  {
2100  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2101  currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2102  shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2103 
2104  /// Check if we can use previous results
2105  if( useLastResults )
2106  {
2107  bComputed = false;
2108 
2109  if( currentInnerSpel == lastInnerSpel )
2110  {
2111  m = lastInnerSum;
2112  innerSum = m;
2113 
2114  bComputed = true;
2115  }
2116  else if( currentInnerSpel == lastOuterSpel )
2117  {
2118  m = lastOuterSum;
2119  outerSum = m;
2120 
2121  bComputed = true;
2122  }
2123  else
2124  {
2125  diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2126 
2127  x = diffSpel[ 0 ];
2128  y = diffSpel[ 1 ];
2129  z = diffSpel[ 2 ];
2130  x2 = x * x;
2131  y2 = y * y;
2132  z2 = z * z;
2133  x2y2z2 = x2 + y2 + z2;
2134 
2135  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2136 
2137  if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2138  {
2139  useLastResults = false;
2140 #ifdef DEBUG_VERBOSE
2141  recount = true;
2142 #endif
2143  }
2144  else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2145  {
2146  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2147  }
2148  else
2149  {
2150  useLastResults = true;
2151  }
2152  }
2153  }
2154 
2155  if( !bComputed )
2156  {
2157  if( !useLastResults ) /// Computation on full kernel, we have no previous results
2158  {
2159  m = NumberTraits< Quantity >::ZERO;
2160 
2161  if( isInitFullMasks )
2162  {
2163  for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2164  {
2165  auto preShiftedSpel = KPS::sSpel( *itm );
2166  preShiftedSpel.coordinates += shiftInnerSpel;
2167 
2168  if( myKSpace.sIsInside( preShiftedSpel ) )
2169  {
2170  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2171 
2172  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2173  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2174  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2175 
2176  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2177  {
2178  m += 1.0;
2179  }
2180  }
2181  }
2182  }
2183  else if( isInitKernelAndMasks )
2184  {
2185  Domain domain = myKernel->getDomain();
2186  for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2187  {
2188  if( myKernel->operator()( *itm ))
2189  {
2190  auto preShiftedSpel = KPS::sSpel( *itm );
2191  preShiftedSpel.coordinates += shiftInnerSpel;
2192 
2193  if( myKSpace.sIsInside( preShiftedSpel ) )
2194  {
2195  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2196 
2197  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2198  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2199  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2200 
2201  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2202  {
2203  m += 1.0;
2204  }
2205  }
2206  }
2207  }
2208  }
2209  else
2210  {
2211  trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2212  return false;
2213  }
2214  }
2215  else /// Using lastInnerMoments
2216  {
2217  m = lastInnerSum;
2218 
2219  /// Part to substract from previous result.
2220  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2221  {
2222  auto preShiftedSpel = KPS::sSpel( *itm );
2223  preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2224 
2225  if( myKSpace.sIsInside( preShiftedSpel ) )
2226  {
2227  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2228 
2229  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2230  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2231  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2232 
2233  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2234  {
2235  m -= 1.0;
2236  }
2237  }
2238  }
2239 
2240  /// Part to add from previous result.
2241  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2242  {
2243  auto preShiftedSpel = KPS::sSpel( *itm );
2244  preShiftedSpel.coordinates += shiftInnerSpel;
2245 
2246  if( myKSpace.sIsInside( preShiftedSpel ) )
2247  {
2248  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2249 
2250  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2251  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2252  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2253 
2254  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2255  {
2256  m += 1.0;
2257  }
2258  }
2259  }
2260  }
2261 
2262  /// Computation of covariance Matrix
2263  innerSum = m;
2264  lastInnerSum = m;
2265  }
2266  }
2267 
2268  /// Outter cell
2269  {
2270  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2271  currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2272  shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2273 
2274  ASSERT( currentInnerSpel != currentOuterSpel );
2275 
2276  diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2277 
2278  x = diffSpel[ 0 ];
2279  y = diffSpel[ 1 ];
2280  z = diffSpel[ 2 ];
2281  x2 = x * x;
2282  y2 = y * y;
2283  z2 = z * z;
2284  x2y2z2 = x2 + y2 + z2;
2285 
2286  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2287 
2288  if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
2289  {
2290  trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2291  }
2292  else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2293  {
2294  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2295  }
2296  else
2297  {
2298  /// Part to substract from previous result.
2299  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2300  {
2301  auto preShiftedSpel = KPS::sSpel( *itm );
2302  preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2303 
2304  if( myKSpace.sIsInside( preShiftedSpel ) )
2305  {
2306  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates
2307  );
2308 
2309  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2310  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2311  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2312 
2313  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2314  {
2315  m -= 1.0;
2316  }
2317  }
2318  }
2319 
2320  /// Part to add from previous result.
2321  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2322  {
2323  auto preShiftedSpel = KPS::sSpel( *itm );
2324  preShiftedSpel.coordinates += shiftOuterSpel;
2325 
2326  if( myKSpace.sIsInside( preShiftedSpel ) )
2327  {
2328  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2329 
2330  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2331  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2332  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2333 
2334  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2335  {
2336  m += 1.0;
2337  }
2338  }
2339  }
2340  }
2341 
2342  /// Computation of covariance Matrix
2343  outerSum = m;
2344  lastOuterSum = m;
2345  }
2346 
2347  ASSERT (( lastInnerSum != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2348 
2349  lastInnerSpel = currentInnerSpel;
2350  lastOuterSpel = currentOuterSpel;
2351 
2352 #ifdef DEBUG_VERBOSE
2353  return recount;
2354 #else
2355  return false;
2356 #endif
2357 }
2358 
2359 
2360 template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2361 template< typename SurfelIterator >
2362 bool
2363 DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_evalCovarianceMatrix
2364 ( const SurfelIterator & it,
2365  CovarianceMatrix & innerMatrix,
2366  CovarianceMatrix & outerMatrix,
2367  bool useLastResults,
2368  Spel & lastInnerSpel,
2369  Spel & lastOuterSpel,
2370  Quantity * lastInnerMoments,
2371  Quantity * lastOuterMoments ) const
2372 {
2373  ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2374 
2375  using KPS = typename KSpace::PreCellularGridSpace;
2376 
2377 #ifdef DEBUG_VERBOSE
2378  Dimension recount = false;
2379 #endif
2380 
2381  typedef typename Functor::Quantity FQuantity;
2382  DGtal::Dimension nbMasks = static_cast<Dimension>(myMasks->size()) - 1;
2383  DGtal::Dimension positionOfFullKernel = 13;
2384 
2385  Quantity m[ nbMoments ]; /// <! [ m000, m001, m010, m100, m011, m101, m110, m002, m020, m200 ]
2386  std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2387 
2388  Spel currentInnerSpel, currentOuterSpel;
2389  Spel shiftedSpel;
2390  Point shiftInnerSpel, shiftOuterSpel;
2391  Point diffSpel;
2392 
2393  bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2394 
2395  int x, y, z, x2, y2, z2, x2y2z2;
2396  unsigned int offset;
2397 
2398  /// Inner cell
2399  {
2400  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2401  currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2402  shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2403 
2404  /// Check if we can use previous results
2405  if( useLastResults )
2406  {
2407  bComputed = false;
2408 
2409  if( currentInnerSpel == lastInnerSpel )
2410  {
2411  memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2412  computeCovarianceMatrix( m, innerMatrix );
2413 
2414  bComputed = true;
2415  }
2416  else if( currentInnerSpel == lastOuterSpel )
2417  {
2418  memcpy( m, lastOuterMoments, nbMoments * sizeof( Quantity ));
2419  computeCovarianceMatrix( m, outerMatrix );
2420 
2421  bComputed = true;
2422  }
2423  else
2424  {
2425  diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2426 
2427  x = diffSpel[ 0 ];
2428  y = diffSpel[ 1 ];
2429  z = diffSpel[ 2 ];
2430  x2 = x * x;
2431  y2 = y * y;
2432  z2 = z * z;
2433  x2y2z2 = x2 + y2 + z2;
2434 
2435  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2436 
2437  if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2438  {
2439  useLastResults = false;
2440 #ifdef DEBUG_VERBOSE
2441  recount = true;
2442 #endif
2443  }
2444  else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2445  {
2446  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2447  }
2448  else
2449  {
2450  useLastResults = true;
2451  }
2452  }
2453  }
2454 
2455  if( !bComputed )
2456  {
2457  if( !useLastResults ) /// Computation on full kernel, we have no previous results
2458  {
2459  std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2460 
2461  if( isInitFullMasks )
2462  {
2463  for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2464  {
2465  auto preShiftedSpel = KPS::sSpel( *itm );
2466  preShiftedSpel.coordinates += shiftInnerSpel;
2467 
2468  if( myKSpace.sIsInside( preShiftedSpel ) )
2469  {
2470  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2471 
2472  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2473  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2474  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2475 
2476  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2477  {
2478  fillMoments( m, shiftedSpel, 1.0 );
2479  }
2480  }
2481  }
2482  }
2483  else if( isInitKernelAndMasks )
2484  {
2485  Domain domain = myKernel->getDomain();
2486  for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2487  {
2488  if( myKernel->operator()( *itm ))
2489  {
2490  auto preShiftedSpel = KPS::sSpel( *itm );
2491  preShiftedSpel.coordinates += shiftInnerSpel;
2492 
2493  if( myKSpace.sIsInside( preShiftedSpel ) )
2494  {
2495  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2496 
2497  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2498  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2499  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2500 
2501  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2502  {
2503  fillMoments( m, shiftedSpel, 1.0 );
2504  }
2505  }
2506  }
2507  }
2508  }
2509  else
2510  {
2511  trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2512  return false;
2513  }
2514  }
2515  else /// Using lastInnerMoments
2516  {
2517  memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2518 
2519  /// Part to substract from previous result.
2520  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2521  {
2522  auto preShiftedSpel = KPS::sSpel( *itm );
2523  preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2524 
2525  if( myKSpace.sIsInside( preShiftedSpel ) )
2526  {
2527  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2528 
2529  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2530  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2531  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2532 
2533  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2534  {
2535  fillMoments( m, shiftedSpel, -1.0 );
2536  }
2537  }
2538  }
2539 
2540  /// Part to add from previous result.
2541  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2542  {
2543  auto preShiftedSpel = KPS::sSpel( *itm );
2544  preShiftedSpel.coordinates += shiftInnerSpel;
2545 
2546  if( myKSpace.sIsInside( preShiftedSpel ) )
2547  {
2548  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2549 
2550  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2551  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2552  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2553 
2554  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2555  {
2556  fillMoments( m, shiftedSpel, 1.0 );
2557  }
2558  }
2559  }
2560  }
2561 
2562  /// Computation of covariance Matrix
2563  computeCovarianceMatrix( m, innerMatrix );
2564  memcpy( lastInnerMoments, m, nbMoments * sizeof( Quantity ));
2565  }
2566  }
2567 
2568  /// Outter cell
2569  {
2570  DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2571  currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2572  shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2573 
2574  ASSERT( currentInnerSpel != currentOuterSpel );
2575 
2576  diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2577 
2578  x = diffSpel[ 0 ];
2579  y = diffSpel[ 1 ];
2580  z = diffSpel[ 2 ];
2581  x2 = x * x;
2582  y2 = y * y;
2583  z2 = z * z;
2584  x2y2z2 = x2 + y2 + z2;
2585 
2586  offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2587 
2588  if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
2589  {
2590  trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2591  }
2592  else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2593  {
2594  trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2595  }
2596  else
2597  {
2598  /// Part to substract from previous result.
2599  for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2600  {
2601  auto preShiftedSpel = KPS::sSpel( *itm );
2602  preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2603 
2604  if( myKSpace.sIsInside( preShiftedSpel ) )
2605  {
2606  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2607 
2608  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2609  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2610  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2611 
2612  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2613  {
2614  fillMoments( m, shiftedSpel, -1.0 );
2615  }
2616  }
2617  }
2618 
2619  /// Part to add from previous result.
2620  for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2621  {
2622  auto preShiftedSpel = KPS::sSpel( *itm );
2623  preShiftedSpel.coordinates += shiftOuterSpel;
2624 
2625  if( myKSpace.sIsInside( preShiftedSpel ) )
2626  {
2627  myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2628 
2629  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2630  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2631  ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2632 
2633  if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2634  {
2635  fillMoments( m, shiftedSpel, 1.0 );
2636  }
2637  }
2638  }
2639  }
2640 
2641  /// Computation of covariance Matrix
2642  computeCovarianceMatrix( m, outerMatrix );
2643  memcpy( lastOuterMoments, m, nbMoments * sizeof( Quantity ));
2644  }
2645 
2646  ASSERT (( lastInnerMoments[ 0 ] != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2647 
2648  lastInnerSpel = currentInnerSpel;
2649  lastOuterSpel = currentOuterSpel;
2650 
2651 #ifdef DEBUG_VERBOSE
2652  return recount;
2653 #else
2654  return false;
2655 #endif
2656 }