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