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