DGtal 1.4.2
Loading...
Searching...
No Matches
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
52template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
53inline
54DGtal::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
67template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
68inline
69DGtal::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
84template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
85inline
86void
87DGtal::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
102template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
103inline
104void
105DGtal::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
120template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
121template< typename SurfelIterator >
122inline
123typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
124DGtal::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
137template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
138template< typename SurfelIterator, typename EvalFunctor >
139inline
140typename EvalFunctor::Value
141DGtal::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
157template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
158template< typename SurfelIterator, typename OutputIterator >
159inline
160void
161DGtal::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
209template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
210template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
211inline
212void
213DGtal::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
268template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
269template< typename SurfelIterator >
270inline
271typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::CovarianceMatrix
272DGtal::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
285template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
286template< typename SurfelIterator, typename EvalFunctor >
287inline
288typename EvalFunctor::Value
289DGtal::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
305template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
306template< typename SurfelIterator, typename OutputIterator >
307inline
308void
309DGtal::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
357template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
358template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
359inline
360void
361DGtal::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
415template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
416inline
417bool
418DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::isValid() const
419{
420 return true;
421}
422
423template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
424void
425DGtal::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
442template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
443void
444DGtal::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
476template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
477const int
478DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::nbMoments = 6;
479#endif //_MSC_VER
480
481template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
482typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Spel
483DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultInnerSpel = Spel();
484
485template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
486typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Spel
487DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultOuterSpel = Spel();
488
489template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
490typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
491DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultInnerMoments[ 6 ] = {Quantity(0)};
492
493template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
494typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
495DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultOuterMoments[ 6 ] = {Quantity(0)};
496
497template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
498typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
499DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultInnerSum = Quantity(0);
500
501template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
502typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::Quantity
503DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, dimension >::defaultOuterSum = Quantity(0);
504
505template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
506template< typename SurfelIterator >
507bool
508DGtal::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
531template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel, DGtal::Dimension dimension >
532template< typename SurfelIterator >
533bool
534DGtal::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
562template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
563inline
564DGtal::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
578template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel>
579inline
580DGtal::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
596template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
597inline
598void
599DGtal::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
614template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
615inline
616void
617DGtal::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
632template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
633template< typename SurfelIterator >
634inline
635typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
636DGtal::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
649template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
650template< typename SurfelIterator, typename EvalFunctor >
651inline
652typename EvalFunctor::Value
653DGtal::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
669template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
670template< typename SurfelIterator, typename OutputIterator >
671inline
672void
673DGtal::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
721template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
722template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
723inline
724void
725DGtal::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
780template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
781template< typename SurfelIterator >
782inline
783typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::CovarianceMatrix
784DGtal::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
797template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
798template< typename SurfelIterator, typename EvalFunctor >
799inline
800typename EvalFunctor::Value
801DGtal::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
817template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
818template< typename SurfelIterator, typename OutputIterator >
819inline
820void
821DGtal::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
869template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
870template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
871inline
872void
873DGtal::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
927template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
928inline
929bool
930DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::isValid() const
931{
932 return true;
933}
934
935template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
936void
937DGtal::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
954template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
955void
956DGtal::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
981template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
982const int
983DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::nbMoments = 6;
984#endif //_MSC_VER
985
986template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
987typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Spel
988DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultInnerSpel = Spel();
989
990template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
991typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Spel
992DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultOuterSpel = Spel();
993
994template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
995typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
996DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultInnerMoments[ 6 ] = {Quantity(0)};
997
998template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
999typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
1000DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultOuterMoments[ 6 ] = {Quantity(0)};
1001
1002template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1003typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
1004DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultInnerSum = Quantity(0);
1005
1006template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1007typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::Quantity
1008DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::defaultOuterSum = Quantity(0);
1009
1010template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1011template< typename SurfelIterator >
1012bool
1013DGtal::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
1299template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1300template< typename SurfelIterator >
1301bool
1302DGtal::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
1593template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1594inline
1595DGtal::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
1609template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel>
1610inline
1611DGtal::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
1627template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1628inline
1629void
1630DGtal::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
1645template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1646inline
1647void
1648DGtal::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
1666template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1667template< typename SurfelIterator >
1668inline
1669typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
1670DGtal::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
1683template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1684template< typename SurfelIterator, typename EvalFunctor >
1685inline
1686typename EvalFunctor::Value
1687DGtal::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
1703template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1704template< typename SurfelIterator, typename OutputIterator >
1705inline
1706void
1707DGtal::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
1755template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1756template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
1757inline
1758void
1759DGtal::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
1812template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1813template< typename SurfelIterator >
1814inline
1815typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::CovarianceMatrix
1816DGtal::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
1829template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1830template< typename SurfelIterator, typename EvalFunctor >
1831inline
1832typename EvalFunctor::Value
1833DGtal::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
1849template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1850template< typename SurfelIterator, typename OutputIterator >
1851inline
1852void
1853DGtal::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
1901template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1902template< typename SurfelIterator, typename OutputIterator, typename EvalFunctor >
1903inline
1904void
1905DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::evalCovarianceMatrix
1906( const SurfelIterator & itbegin,
1907 const SurfelIterator & itend,
1908 OutputIterator & result,
1909 EvalFunctor functor ) const
1910{
1911 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1912
1913 Dimension total = 0;
1914#ifdef DEBUG_VERBOSE
1915 Dimension recount = 0;
1916#endif
1917
1918 Quantity *lastInnerMoments = new Quantity[ nbMoments ];
1919 Quantity *lastOuterMoments = new Quantity[ nbMoments ];
1920
1921 CovarianceMatrix innerCovarianceMatrix, outerCovarianceMatrix;
1922 CovarianceMatrix resultCovarianceMatrix;
1923
1924 Spel lastInnerSpel, lastOuterSpel;
1925
1926 /// Iterate on all cells
1927 for( SurfelIterator it = itbegin; it != itend; ++it )
1928 {
1929 if( total != 0 )
1930 {
1931#ifdef DEBUG_VERBOSE
1932 bool hasJumped = core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1933 recount = ( hasJumped ) ? recount + 1 : recount;
1934#else
1935 core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, true, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1936#endif
1937 }
1938 else
1939 {
1940 core_evalCovarianceMatrix( it, innerCovarianceMatrix, outerCovarianceMatrix, false, lastInnerSpel, lastOuterSpel, lastInnerMoments, lastOuterMoments );
1941 }
1942
1943 double lambda = 0.5;
1944 resultCovarianceMatrix = innerCovarianceMatrix * lambda + outerCovarianceMatrix * ( 1.0 - lambda );
1945 *result++ = functor( resultCovarianceMatrix );
1946
1947 ++total;
1948 }
1949
1950#ifdef DEBUG_VERBOSE
1951 std::cout << "#total cells = " << total << std::endl;
1952 std::cout << "#recount = " << recount << std::endl;
1953#endif
1954
1955 delete[] lastOuterMoments;
1956 delete[] lastInnerMoments;
1957}
1958
1959
1960
1961template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1962inline
1963bool
1964DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::isValid() const
1965{
1966 return true;
1967}
1968
1969template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1970void
1971DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::fillMoments
1972( Quantity * aMomentMatrix,
1973 const Spel & aSpel,
1974 double direction ) const
1975{
1976 RealPoint current = myEmbedder( aSpel );
1977 double x = current[ 0 ];
1978 double y = current[ 1 ];
1979 double z = current[ 2 ];
1980
1981 aMomentMatrix[ 0 ] += direction * 1;
1982 aMomentMatrix[ 1 ] += direction * z;
1983 aMomentMatrix[ 2 ] += direction * y;
1984 aMomentMatrix[ 3 ] += direction * x;
1985 aMomentMatrix[ 4 ] += direction * y * z;
1986 aMomentMatrix[ 5 ] += direction * x * z;
1987 aMomentMatrix[ 6 ] += direction * x * y;
1988 aMomentMatrix[ 7 ] += direction * z * z;
1989 aMomentMatrix[ 8 ] += direction * y * y;
1990 aMomentMatrix[ 9 ] += direction * x * x;
1991
1992
1993}
1994
1995template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1996void
1997DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::computeCovarianceMatrix
1998( const Quantity * aMomentMatrix,
1999 CovarianceMatrix & aCovarianceMatrix ) const
2000{
2001 MatrixQuantity A;
2002 double B;
2003 MatrixQuantity C;
2004
2005 A.setComponent( 0, 0, aMomentMatrix[ 9 ] );
2006 A.setComponent( 0, 1, aMomentMatrix[ 6 ] );
2007 A.setComponent( 0, 2, aMomentMatrix[ 5 ] );
2008 A.setComponent( 1, 0, aMomentMatrix[ 6 ] );
2009 A.setComponent( 1, 1, aMomentMatrix[ 8 ] );
2010 A.setComponent( 1, 2, aMomentMatrix[ 4 ] );
2011 A.setComponent( 2, 0, aMomentMatrix[ 5 ] );
2012 A.setComponent( 2, 1, aMomentMatrix[ 4 ] );
2013 A.setComponent( 2, 2, aMomentMatrix[ 7 ] );
2014
2015 B = 1.0 / aMomentMatrix[ 0 ];
2016
2017 C.setComponent( 0, 0, aMomentMatrix[ 3 ] * aMomentMatrix[ 3 ] );
2018 C.setComponent( 0, 1, aMomentMatrix[ 3 ] * aMomentMatrix[ 2 ] );
2019 C.setComponent( 0, 2, aMomentMatrix[ 3 ] * aMomentMatrix[ 1 ] );
2020 C.setComponent( 1, 0, aMomentMatrix[ 2 ] * aMomentMatrix[ 3 ] );
2021 C.setComponent( 1, 1, aMomentMatrix[ 2 ] * aMomentMatrix[ 2 ] );
2022 C.setComponent( 1, 2, aMomentMatrix[ 2 ] * aMomentMatrix[ 1 ] );
2023 C.setComponent( 2, 0, aMomentMatrix[ 1 ] * aMomentMatrix[ 3 ] );
2024 C.setComponent( 2, 1, aMomentMatrix[ 1 ] * aMomentMatrix[ 2 ] );
2025 C.setComponent( 2, 2, aMomentMatrix[ 1 ] * aMomentMatrix[ 1 ] );
2026
2027 aCovarianceMatrix = A - C * B;
2028}
2029
2030#ifndef _MSC_VER
2031// For Visual Studio, to be defined as a static const, it has to be intialized into the header file
2032template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2033const int
2034DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::nbMoments = 10;
2035#endif //_MSC_VER
2036
2037template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2038typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Spel
2039DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerSpel = Spel();
2040
2041template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2042typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Spel
2043DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterSpel = Spel();
2044
2045template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2046typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2047DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerMoments[ 10 ] = {Quantity(0)};
2048
2049template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2050typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2051DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterMoments[ 10 ] = {Quantity(0)};
2052
2053template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2054typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2055DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerSum = Quantity(0);
2056
2057template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2058typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2059DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterSum = Quantity(0);
2060
2061template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2062template< typename SurfelIterator >
2063bool
2064DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_eval
2065( const SurfelIterator & it,
2066 Quantity & innerSum,
2067 Quantity & outerSum,
2068 bool useLastResults,
2069 Spel & lastInnerSpel,
2070 Spel & lastOuterSpel,
2071 Quantity & lastInnerSum,
2072 Quantity & lastOuterSum ) const
2073{
2074 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2075
2076 using KPS = typename KSpace::PreCellularGridSpace;
2077
2078#ifdef DEBUG_VERBOSE
2079 Dimension recount = false;
2080#endif
2081
2082 typedef typename Functor::Quantity FQuantity;
2083 DGtal::Dimension nbMasks = static_cast<DGtal::Dimension>(myMasks->size()) - 1;
2084 DGtal::Dimension positionOfFullKernel = 13;
2085
2086 Quantity m = NumberTraits< Quantity >::ZERO;
2087
2088 Spel currentInnerSpel, currentOuterSpel;
2089 Spel shiftedSpel;
2090 Point shiftInnerSpel, shiftOuterSpel;
2091 Point diffSpel;
2092
2093 bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2094
2095 int x, y, z, x2, y2, z2, x2y2z2;
2096 unsigned int offset;
2097
2098 /// Inner cell
2099 {
2100 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2101 currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2102 shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2103
2104 /// Check if we can use previous results
2105 if( useLastResults )
2106 {
2107 bComputed = false;
2108
2109 if( currentInnerSpel == lastInnerSpel )
2110 {
2111 m = lastInnerSum;
2112 innerSum = m;
2113
2114 bComputed = true;
2115 }
2116 else if( currentInnerSpel == lastOuterSpel )
2117 {
2118 m = lastOuterSum;
2119 outerSum = m;
2120
2121 bComputed = true;
2122 }
2123 else
2124 {
2125 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2126
2127 x = diffSpel[ 0 ];
2128 y = diffSpel[ 1 ];
2129 z = diffSpel[ 2 ];
2130 x2 = x * x;
2131 y2 = y * y;
2132 z2 = z * z;
2133 x2y2z2 = x2 + y2 + z2;
2134
2135 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2136
2137 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2138 {
2139 useLastResults = false;
2140#ifdef DEBUG_VERBOSE
2141 recount = true;
2142#endif
2143 }
2144 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2145 {
2146 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2147 }
2148 else
2149 {
2150 useLastResults = true;
2151 }
2152 }
2153 }
2154
2155 if( !bComputed )
2156 {
2157 if( !useLastResults ) /// Computation on full kernel, we have no previous results
2158 {
2159 m = NumberTraits< Quantity >::ZERO;
2160
2161 if( isInitFullMasks )
2162 {
2163 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2164 {
2165 auto preShiftedSpel = KPS::sSpel( *itm );
2166 preShiftedSpel.coordinates += shiftInnerSpel;
2167
2168 if( myKSpace.sIsInside( preShiftedSpel ) )
2169 {
2170 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2171
2172 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2173 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2174 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2175
2176 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2177 {
2178 m += 1.0;
2179 }
2180 }
2181 }
2182 }
2183 else if( isInitKernelAndMasks )
2184 {
2185 Domain domain = myKernel->getDomain();
2186 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2187 {
2188 if( myKernel->operator()( *itm ))
2189 {
2190 auto preShiftedSpel = KPS::sSpel( *itm );
2191 preShiftedSpel.coordinates += shiftInnerSpel;
2192
2193 if( myKSpace.sIsInside( preShiftedSpel ) )
2194 {
2195 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2196
2197 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2198 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2199 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2200
2201 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2202 {
2203 m += 1.0;
2204 }
2205 }
2206 }
2207 }
2208 }
2209 else
2210 {
2211 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2212 return false;
2213 }
2214 }
2215 else /// Using lastInnerMoments
2216 {
2217 m = lastInnerSum;
2218
2219 /// Part to substract from previous result.
2220 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2221 {
2222 auto preShiftedSpel = KPS::sSpel( *itm );
2223 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2224
2225 if( myKSpace.sIsInside( preShiftedSpel ) )
2226 {
2227 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2228
2229 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2230 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2231 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2232
2233 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2234 {
2235 m -= 1.0;
2236 }
2237 }
2238 }
2239
2240 /// Part to add from previous result.
2241 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2242 {
2243 auto preShiftedSpel = KPS::sSpel( *itm );
2244 preShiftedSpel.coordinates += shiftInnerSpel;
2245
2246 if( myKSpace.sIsInside( preShiftedSpel ) )
2247 {
2248 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2249
2250 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2251 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2252 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2253
2254 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2255 {
2256 m += 1.0;
2257 }
2258 }
2259 }
2260 }
2261
2262 /// Computation of covariance Matrix
2263 innerSum = m;
2264 lastInnerSum = m;
2265 }
2266 }
2267
2268 /// Outter cell
2269 {
2270 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2271 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2272 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2273
2274 ASSERT( currentInnerSpel != currentOuterSpel );
2275
2276 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2277
2278 x = diffSpel[ 0 ];
2279 y = diffSpel[ 1 ];
2280 z = diffSpel[ 2 ];
2281 x2 = x * x;
2282 y2 = y * y;
2283 z2 = z * z;
2284 x2y2z2 = x2 + y2 + z2;
2285
2286 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2287
2288 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
2289 {
2290 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2291 }
2292 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2293 {
2294 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2295 }
2296 else
2297 {
2298 /// Part to substract from previous result.
2299 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2300 {
2301 auto preShiftedSpel = KPS::sSpel( *itm );
2302 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2303
2304 if( myKSpace.sIsInside( preShiftedSpel ) )
2305 {
2306 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates
2307 );
2308
2309 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2310 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2311 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2312
2313 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2314 {
2315 m -= 1.0;
2316 }
2317 }
2318 }
2319
2320 /// Part to add from previous result.
2321 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2322 {
2323 auto preShiftedSpel = KPS::sSpel( *itm );
2324 preShiftedSpel.coordinates += shiftOuterSpel;
2325
2326 if( myKSpace.sIsInside( preShiftedSpel ) )
2327 {
2328 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2329
2330 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2331 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2332 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2333
2334 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2335 {
2336 m += 1.0;
2337 }
2338 }
2339 }
2340 }
2341
2342 /// Computation of covariance Matrix
2343 outerSum = m;
2344 lastOuterSum = m;
2345 }
2346
2347 ASSERT (( lastInnerSum != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2348
2349 lastInnerSpel = currentInnerSpel;
2350 lastOuterSpel = currentOuterSpel;
2351
2352#ifdef DEBUG_VERBOSE
2353 return recount;
2354#else
2355 return false;
2356#endif
2357}
2358
2359
2360template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2361template< typename SurfelIterator >
2362bool
2363DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_evalCovarianceMatrix
2364( const SurfelIterator & it,
2365 CovarianceMatrix & innerMatrix,
2366 CovarianceMatrix & outerMatrix,
2367 bool useLastResults,
2368 Spel & lastInnerSpel,
2369 Spel & lastOuterSpel,
2370 Quantity * lastInnerMoments,
2371 Quantity * lastOuterMoments ) const
2372{
2373 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2374
2375 using KPS = typename KSpace::PreCellularGridSpace;
2376
2377#ifdef DEBUG_VERBOSE
2378 Dimension recount = false;
2379#endif
2380
2381 typedef typename Functor::Quantity FQuantity;
2382 DGtal::Dimension nbMasks = static_cast<Dimension>(myMasks->size()) - 1;
2383 DGtal::Dimension positionOfFullKernel = 13;
2384
2385 Quantity m[ nbMoments ]; /// <! [ m000, m001, m010, m100, m011, m101, m110, m002, m020, m200 ]
2386 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2387
2388 Spel currentInnerSpel, currentOuterSpel;
2389 Spel shiftedSpel;
2390 Point shiftInnerSpel, shiftOuterSpel;
2391 Point diffSpel;
2392
2393 bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2394
2395 int x, y, z, x2, y2, z2, x2y2z2;
2396 unsigned int offset;
2397
2398 /// Inner cell
2399 {
2400 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2401 currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2402 shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2403
2404 /// Check if we can use previous results
2405 if( useLastResults )
2406 {
2407 bComputed = false;
2408
2409 if( currentInnerSpel == lastInnerSpel )
2410 {
2411 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2412 computeCovarianceMatrix( m, innerMatrix );
2413
2414 bComputed = true;
2415 }
2416 else if( currentInnerSpel == lastOuterSpel )
2417 {
2418 memcpy( m, lastOuterMoments, nbMoments * sizeof( Quantity ));
2419 computeCovarianceMatrix( m, outerMatrix );
2420
2421 bComputed = true;
2422 }
2423 else
2424 {
2425 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2426
2427 x = diffSpel[ 0 ];
2428 y = diffSpel[ 1 ];
2429 z = diffSpel[ 2 ];
2430 x2 = x * x;
2431 y2 = y * y;
2432 z2 = z * z;
2433 x2y2z2 = x2 + y2 + z2;
2434
2435 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2436
2437 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2438 {
2439 useLastResults = false;
2440#ifdef DEBUG_VERBOSE
2441 recount = true;
2442#endif
2443 }
2444 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2445 {
2446 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2447 }
2448 else
2449 {
2450 useLastResults = true;
2451 }
2452 }
2453 }
2454
2455 if( !bComputed )
2456 {
2457 if( !useLastResults ) /// Computation on full kernel, we have no previous results
2458 {
2459 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2460
2461 if( isInitFullMasks )
2462 {
2463 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2464 {
2465 auto preShiftedSpel = KPS::sSpel( *itm );
2466 preShiftedSpel.coordinates += shiftInnerSpel;
2467
2468 if( myKSpace.sIsInside( preShiftedSpel ) )
2469 {
2470 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2471
2472 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2473 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2474 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2475
2476 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2477 {
2478 fillMoments( m, shiftedSpel, 1.0 );
2479 }
2480 }
2481 }
2482 }
2483 else if( isInitKernelAndMasks )
2484 {
2485 Domain domain = myKernel->getDomain();
2486 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2487 {
2488 if( myKernel->operator()( *itm ))
2489 {
2490 auto preShiftedSpel = KPS::sSpel( *itm );
2491 preShiftedSpel.coordinates += shiftInnerSpel;
2492
2493 if( myKSpace.sIsInside( preShiftedSpel ) )
2494 {
2495 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2496
2497 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2498 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2499 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2500
2501 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2502 {
2503 fillMoments( m, shiftedSpel, 1.0 );
2504 }
2505 }
2506 }
2507 }
2508 }
2509 else
2510 {
2511 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2512 return false;
2513 }
2514 }
2515 else /// Using lastInnerMoments
2516 {
2517 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2518
2519 /// Part to substract from previous result.
2520 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2521 {
2522 auto preShiftedSpel = KPS::sSpel( *itm );
2523 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2524
2525 if( myKSpace.sIsInside( preShiftedSpel ) )
2526 {
2527 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2528
2529 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2530 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2531 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2532
2533 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2534 {
2535 fillMoments( m, shiftedSpel, -1.0 );
2536 }
2537 }
2538 }
2539
2540 /// Part to add from previous result.
2541 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2542 {
2543 auto preShiftedSpel = KPS::sSpel( *itm );
2544 preShiftedSpel.coordinates += shiftInnerSpel;
2545
2546 if( myKSpace.sIsInside( preShiftedSpel ) )
2547 {
2548 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2549
2550 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2551 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2552 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2553
2554 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2555 {
2556 fillMoments( m, shiftedSpel, 1.0 );
2557 }
2558 }
2559 }
2560 }
2561
2562 /// Computation of covariance Matrix
2563 computeCovarianceMatrix( m, innerMatrix );
2564 memcpy( lastInnerMoments, m, nbMoments * sizeof( Quantity ));
2565 }
2566 }
2567
2568 /// Outter cell
2569 {
2570 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2571 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2572 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2573
2574 ASSERT( currentInnerSpel != currentOuterSpel );
2575
2576 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2577
2578 x = diffSpel[ 0 ];
2579 y = diffSpel[ 1 ];
2580 z = diffSpel[ 2 ];
2581 x2 = x * x;
2582 y2 = y * y;
2583 z2 = z * z;
2584 x2y2z2 = x2 + y2 + z2;
2585
2586 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2587
2588 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
2589 {
2590 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2591 }
2592 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2593 {
2594 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2595 }
2596 else
2597 {
2598 /// Part to substract from previous result.
2599 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2600 {
2601 auto preShiftedSpel = KPS::sSpel( *itm );
2602 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2603
2604 if( myKSpace.sIsInside( preShiftedSpel ) )
2605 {
2606 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2607
2608 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2609 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2610 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2611
2612 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2613 {
2614 fillMoments( m, shiftedSpel, -1.0 );
2615 }
2616 }
2617 }
2618
2619 /// Part to add from previous result.
2620 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2621 {
2622 auto preShiftedSpel = KPS::sSpel( *itm );
2623 preShiftedSpel.coordinates += shiftOuterSpel;
2624
2625 if( myKSpace.sIsInside( preShiftedSpel ) )
2626 {
2627 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2628
2629 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2630 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2631 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2632
2633 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2634 {
2635 fillMoments( m, shiftedSpel, 1.0 );
2636 }
2637 }
2638 }
2639 }
2640
2641 /// Computation of covariance Matrix
2642 computeCovarianceMatrix( m, outerMatrix );
2643 memcpy( lastOuterMoments, m, nbMoments * sizeof( Quantity ));
2644 }
2645
2646 ASSERT (( lastInnerMoments[ 0 ] != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2647
2648 lastInnerSpel = currentInnerSpel;
2649 lastOuterSpel = currentOuterSpel;
2650
2651#ifdef DEBUG_VERBOSE
2652 return recount;
2653#else
2654 return false;
2655#endif
2656}