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