DGtal 1.4.0
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
1956
1957
1958template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1959inline
1960bool
1961DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::isValid() const
1962{
1963 return true;
1964}
1965
1966template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1967void
1968DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::fillMoments
1969( Quantity * aMomentMatrix,
1970 const Spel & aSpel,
1971 double direction ) const
1972{
1973 RealPoint current = myEmbedder( aSpel );
1974 double x = current[ 0 ];
1975 double y = current[ 1 ];
1976 double z = current[ 2 ];
1977
1978 aMomentMatrix[ 0 ] += direction * 1;
1979 aMomentMatrix[ 1 ] += direction * z;
1980 aMomentMatrix[ 2 ] += direction * y;
1981 aMomentMatrix[ 3 ] += direction * x;
1982 aMomentMatrix[ 4 ] += direction * y * z;
1983 aMomentMatrix[ 5 ] += direction * x * z;
1984 aMomentMatrix[ 6 ] += direction * x * y;
1985 aMomentMatrix[ 7 ] += direction * z * z;
1986 aMomentMatrix[ 8 ] += direction * y * y;
1987 aMomentMatrix[ 9 ] += direction * x * x;
1988
1989
1990}
1991
1992template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1993void
1994DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::computeCovarianceMatrix
1995( const Quantity * aMomentMatrix,
1996 CovarianceMatrix & aCovarianceMatrix ) const
1997{
1998 MatrixQuantity A;
1999 double B;
2000 MatrixQuantity C;
2001
2002 A.setComponent( 0, 0, aMomentMatrix[ 9 ] );
2003 A.setComponent( 0, 1, aMomentMatrix[ 6 ] );
2004 A.setComponent( 0, 2, aMomentMatrix[ 5 ] );
2005 A.setComponent( 1, 0, aMomentMatrix[ 6 ] );
2006 A.setComponent( 1, 1, aMomentMatrix[ 8 ] );
2007 A.setComponent( 1, 2, aMomentMatrix[ 4 ] );
2008 A.setComponent( 2, 0, aMomentMatrix[ 5 ] );
2009 A.setComponent( 2, 1, aMomentMatrix[ 4 ] );
2010 A.setComponent( 2, 2, aMomentMatrix[ 7 ] );
2011
2012 B = 1.0 / aMomentMatrix[ 0 ];
2013
2014 C.setComponent( 0, 0, aMomentMatrix[ 3 ] * aMomentMatrix[ 3 ] );
2015 C.setComponent( 0, 1, aMomentMatrix[ 3 ] * aMomentMatrix[ 2 ] );
2016 C.setComponent( 0, 2, aMomentMatrix[ 3 ] * aMomentMatrix[ 1 ] );
2017 C.setComponent( 1, 0, aMomentMatrix[ 2 ] * aMomentMatrix[ 3 ] );
2018 C.setComponent( 1, 1, aMomentMatrix[ 2 ] * aMomentMatrix[ 2 ] );
2019 C.setComponent( 1, 2, aMomentMatrix[ 2 ] * aMomentMatrix[ 1 ] );
2020 C.setComponent( 2, 0, aMomentMatrix[ 1 ] * aMomentMatrix[ 3 ] );
2021 C.setComponent( 2, 1, aMomentMatrix[ 1 ] * aMomentMatrix[ 2 ] );
2022 C.setComponent( 2, 2, aMomentMatrix[ 1 ] * aMomentMatrix[ 1 ] );
2023
2024 aCovarianceMatrix = A - C * B;
2025}
2026
2027#ifndef _MSC_VER
2028// For Visual Studio, to be defined as a static const, it has to be intialized into the header file
2029template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2030const int
2031DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::nbMoments = 10;
2032#endif //_MSC_VER
2033
2034template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2035typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Spel
2036DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerSpel = Spel();
2037
2038template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2039typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Spel
2040DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterSpel = Spel();
2041
2042template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2043typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2044DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerMoments[ 10 ] = {Quantity(0)};
2045
2046template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2047typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2048DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterMoments[ 10 ] = {Quantity(0)};
2049
2050template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2051typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2052DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultInnerSum = Quantity(0);
2053
2054template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2055typename DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::Quantity
2056DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::defaultOuterSum = Quantity(0);
2057
2058template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2059template< typename SurfelIterator >
2060bool
2061DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_eval
2062( const SurfelIterator & it,
2063 Quantity & innerSum,
2064 Quantity & outerSum,
2065 bool useLastResults,
2066 Spel & lastInnerSpel,
2067 Spel & lastOuterSpel,
2068 Quantity & lastInnerSum,
2069 Quantity & lastOuterSum ) const
2070{
2071 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2072
2073 using KPS = typename KSpace::PreCellularGridSpace;
2074
2075#ifdef DEBUG_VERBOSE
2076 Dimension recount = false;
2077#endif
2078
2079 typedef typename Functor::Quantity FQuantity;
2080 DGtal::Dimension nbMasks = static_cast<DGtal::Dimension>(myMasks->size()) - 1;
2081 DGtal::Dimension positionOfFullKernel = 13;
2082
2083 Quantity m = NumberTraits< Quantity >::ZERO;
2084
2085 Spel currentInnerSpel, currentOuterSpel;
2086 Spel shiftedSpel;
2087 Point shiftInnerSpel, shiftOuterSpel;
2088 Point diffSpel;
2089
2090 bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2091
2092 int x, y, z, x2, y2, z2, x2y2z2;
2093 unsigned int offset;
2094
2095 /// Inner cell
2096 {
2097 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2098 currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2099 shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2100
2101 /// Check if we can use previous results
2102 if( useLastResults )
2103 {
2104 bComputed = false;
2105
2106 if( currentInnerSpel == lastInnerSpel )
2107 {
2108 m = lastInnerSum;
2109 innerSum = m;
2110
2111 bComputed = true;
2112 }
2113 else if( currentInnerSpel == lastOuterSpel )
2114 {
2115 m = lastOuterSum;
2116 outerSum = m;
2117
2118 bComputed = true;
2119 }
2120 else
2121 {
2122 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2123
2124 x = diffSpel[ 0 ];
2125 y = diffSpel[ 1 ];
2126 z = diffSpel[ 2 ];
2127 x2 = x * x;
2128 y2 = y * y;
2129 z2 = z * z;
2130 x2y2z2 = x2 + y2 + z2;
2131
2132 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2133
2134 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2135 {
2136 useLastResults = false;
2137#ifdef DEBUG_VERBOSE
2138 recount = true;
2139#endif
2140 }
2141 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2142 {
2143 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2144 }
2145 else
2146 {
2147 useLastResults = true;
2148 }
2149 }
2150 }
2151
2152 if( !bComputed )
2153 {
2154 if( !useLastResults ) /// Computation on full kernel, we have no previous results
2155 {
2156 m = NumberTraits< Quantity >::ZERO;
2157
2158 if( isInitFullMasks )
2159 {
2160 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2161 {
2162 auto preShiftedSpel = KPS::sSpel( *itm );
2163 preShiftedSpel.coordinates += shiftInnerSpel;
2164
2165 if( myKSpace.sIsInside( preShiftedSpel ) )
2166 {
2167 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2168
2169 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2170 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2171 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2172
2173 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2174 {
2175 m += 1.0;
2176 }
2177 }
2178 }
2179 }
2180 else if( isInitKernelAndMasks )
2181 {
2182 Domain domain = myKernel->getDomain();
2183 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2184 {
2185 if( myKernel->operator()( *itm ))
2186 {
2187 auto preShiftedSpel = KPS::sSpel( *itm );
2188 preShiftedSpel.coordinates += shiftInnerSpel;
2189
2190 if( myKSpace.sIsInside( preShiftedSpel ) )
2191 {
2192 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2193
2194 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2195 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2196 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2197
2198 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2199 {
2200 m += 1.0;
2201 }
2202 }
2203 }
2204 }
2205 }
2206 else
2207 {
2208 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2209 return false;
2210 }
2211 }
2212 else /// Using lastInnerMoments
2213 {
2214 m = lastInnerSum;
2215
2216 /// Part to substract from previous result.
2217 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2218 {
2219 auto preShiftedSpel = KPS::sSpel( *itm );
2220 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2221
2222 if( myKSpace.sIsInside( preShiftedSpel ) )
2223 {
2224 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2225
2226 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2227 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2228 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2229
2230 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2231 {
2232 m -= 1.0;
2233 }
2234 }
2235 }
2236
2237 /// Part to add from previous result.
2238 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2239 {
2240 auto preShiftedSpel = KPS::sSpel( *itm );
2241 preShiftedSpel.coordinates += shiftInnerSpel;
2242
2243 if( myKSpace.sIsInside( preShiftedSpel ) )
2244 {
2245 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2246
2247 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2248 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2249 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2250
2251 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2252 {
2253 m += 1.0;
2254 }
2255 }
2256 }
2257 }
2258
2259 /// Computation of covariance Matrix
2260 innerSum = m;
2261 lastInnerSum = m;
2262 }
2263 }
2264
2265 /// Outter cell
2266 {
2267 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2268 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2269 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2270
2271 ASSERT( currentInnerSpel != currentOuterSpel );
2272
2273 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2274
2275 x = diffSpel[ 0 ];
2276 y = diffSpel[ 1 ];
2277 z = diffSpel[ 2 ];
2278 x2 = x * x;
2279 y2 = y * y;
2280 z2 = z * z;
2281 x2y2z2 = x2 + y2 + z2;
2282
2283 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2284
2285 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
2286 {
2287 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2288 }
2289 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2290 {
2291 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2292 }
2293 else
2294 {
2295 /// Part to substract from previous result.
2296 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2297 {
2298 auto preShiftedSpel = KPS::sSpel( *itm );
2299 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2300
2301 if( myKSpace.sIsInside( preShiftedSpel ) )
2302 {
2303 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates
2304 );
2305
2306 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2307 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2308 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2309
2310 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2311 {
2312 m -= 1.0;
2313 }
2314 }
2315 }
2316
2317 /// Part to add from previous result.
2318 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2319 {
2320 auto preShiftedSpel = KPS::sSpel( *itm );
2321 preShiftedSpel.coordinates += shiftOuterSpel;
2322
2323 if( myKSpace.sIsInside( preShiftedSpel ) )
2324 {
2325 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2326
2327 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2328 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2329 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2330
2331 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2332 {
2333 m += 1.0;
2334 }
2335 }
2336 }
2337 }
2338
2339 /// Computation of covariance Matrix
2340 outerSum = m;
2341 lastOuterSum = m;
2342 }
2343
2344 ASSERT (( lastInnerSum != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2345
2346 lastInnerSpel = currentInnerSpel;
2347 lastOuterSpel = currentOuterSpel;
2348
2349#ifdef DEBUG_VERBOSE
2350 return recount;
2351#else
2352 return false;
2353#endif
2354}
2355
2356
2357template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2358template< typename SurfelIterator >
2359bool
2360DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_evalCovarianceMatrix
2361( const SurfelIterator & it,
2362 CovarianceMatrix & innerMatrix,
2363 CovarianceMatrix & outerMatrix,
2364 bool useLastResults,
2365 Spel & lastInnerSpel,
2366 Spel & lastOuterSpel,
2367 Quantity * lastInnerMoments,
2368 Quantity * lastOuterMoments ) const
2369{
2370 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2371
2372 using KPS = typename KSpace::PreCellularGridSpace;
2373
2374#ifdef DEBUG_VERBOSE
2375 Dimension recount = false;
2376#endif
2377
2378 typedef typename Functor::Quantity FQuantity;
2379 DGtal::Dimension nbMasks = static_cast<Dimension>(myMasks->size()) - 1;
2380 DGtal::Dimension positionOfFullKernel = 13;
2381
2382 Quantity m[ nbMoments ]; /// <! [ m000, m001, m010, m100, m011, m101, m110, m002, m020, m200 ]
2383 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2384
2385 Spel currentInnerSpel, currentOuterSpel;
2386 Spel shiftedSpel;
2387 Point shiftInnerSpel, shiftOuterSpel;
2388 Point diffSpel;
2389
2390 bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2391
2392 int x, y, z, x2, y2, z2, x2y2z2;
2393 unsigned int offset;
2394
2395 /// Inner cell
2396 {
2397 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2398 currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2399 shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2400
2401 /// Check if we can use previous results
2402 if( useLastResults )
2403 {
2404 bComputed = false;
2405
2406 if( currentInnerSpel == lastInnerSpel )
2407 {
2408 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2409 computeCovarianceMatrix( m, innerMatrix );
2410
2411 bComputed = true;
2412 }
2413 else if( currentInnerSpel == lastOuterSpel )
2414 {
2415 memcpy( m, lastOuterMoments, nbMoments * sizeof( Quantity ));
2416 computeCovarianceMatrix( m, outerMatrix );
2417
2418 bComputed = true;
2419 }
2420 else
2421 {
2422 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2423
2424 x = diffSpel[ 0 ];
2425 y = diffSpel[ 1 ];
2426 z = diffSpel[ 2 ];
2427 x2 = x * x;
2428 y2 = y * y;
2429 z2 = z * z;
2430 x2y2z2 = x2 + y2 + z2;
2431
2432 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2433
2434 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2435 {
2436 useLastResults = false;
2437#ifdef DEBUG_VERBOSE
2438 recount = true;
2439#endif
2440 }
2441 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2442 {
2443 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2444 }
2445 else
2446 {
2447 useLastResults = true;
2448 }
2449 }
2450 }
2451
2452 if( !bComputed )
2453 {
2454 if( !useLastResults ) /// Computation on full kernel, we have no previous results
2455 {
2456 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2457
2458 if( isInitFullMasks )
2459 {
2460 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2461 {
2462 auto preShiftedSpel = KPS::sSpel( *itm );
2463 preShiftedSpel.coordinates += shiftInnerSpel;
2464
2465 if( myKSpace.sIsInside( preShiftedSpel ) )
2466 {
2467 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2468
2469 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2470 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2471 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2472
2473 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2474 {
2475 fillMoments( m, shiftedSpel, 1.0 );
2476 }
2477 }
2478 }
2479 }
2480 else if( isInitKernelAndMasks )
2481 {
2482 Domain domain = myKernel->getDomain();
2483 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2484 {
2485 if( myKernel->operator()( *itm ))
2486 {
2487 auto preShiftedSpel = KPS::sSpel( *itm );
2488 preShiftedSpel.coordinates += shiftInnerSpel;
2489
2490 if( myKSpace.sIsInside( preShiftedSpel ) )
2491 {
2492 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2493
2494 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2495 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2496 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2497
2498 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2499 {
2500 fillMoments( m, shiftedSpel, 1.0 );
2501 }
2502 }
2503 }
2504 }
2505 }
2506 else
2507 {
2508 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2509 return false;
2510 }
2511 }
2512 else /// Using lastInnerMoments
2513 {
2514 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2515
2516 /// Part to substract from previous result.
2517 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2518 {
2519 auto preShiftedSpel = KPS::sSpel( *itm );
2520 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2521
2522 if( myKSpace.sIsInside( preShiftedSpel ) )
2523 {
2524 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2525
2526 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2527 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2528 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2529
2530 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2531 {
2532 fillMoments( m, shiftedSpel, -1.0 );
2533 }
2534 }
2535 }
2536
2537 /// Part to add from previous result.
2538 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2539 {
2540 auto preShiftedSpel = KPS::sSpel( *itm );
2541 preShiftedSpel.coordinates += shiftInnerSpel;
2542
2543 if( myKSpace.sIsInside( preShiftedSpel ) )
2544 {
2545 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2546
2547 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2548 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2549 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2550
2551 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2552 {
2553 fillMoments( m, shiftedSpel, 1.0 );
2554 }
2555 }
2556 }
2557 }
2558
2559 /// Computation of covariance Matrix
2560 computeCovarianceMatrix( m, innerMatrix );
2561 memcpy( lastInnerMoments, m, nbMoments * sizeof( Quantity ));
2562 }
2563 }
2564
2565 /// Outter cell
2566 {
2567 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2568 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2569 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( myKernelSpelOrigin );
2570
2571 ASSERT( currentInnerSpel != currentOuterSpel );
2572
2573 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2574
2575 x = diffSpel[ 0 ];
2576 y = diffSpel[ 1 ];
2577 z = diffSpel[ 2 ];
2578 x2 = x * x;
2579 y2 = y * y;
2580 z2 = z * z;
2581 x2y2z2 = x2 + y2 + z2;
2582
2583 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2584
2585 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
2586 {
2587 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2588 }
2589 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2590 {
2591 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2592 }
2593 else
2594 {
2595 /// Part to substract from previous result.
2596 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2597 {
2598 auto preShiftedSpel = KPS::sSpel( *itm );
2599 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2600
2601 if( myKSpace.sIsInside( preShiftedSpel ) )
2602 {
2603 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2604
2605 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2606 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2607 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2608
2609 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2610 {
2611 fillMoments( m, shiftedSpel, -1.0 );
2612 }
2613 }
2614 }
2615
2616 /// Part to add from previous result.
2617 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2618 {
2619 auto preShiftedSpel = KPS::sSpel( *itm );
2620 preShiftedSpel.coordinates += shiftOuterSpel;
2621
2622 if( myKSpace.sIsInside( preShiftedSpel ) )
2623 {
2624 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2625
2626 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2627 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2628 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2629
2630 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2631 {
2632 fillMoments( m, shiftedSpel, 1.0 );
2633 }
2634 }
2635 }
2636 }
2637
2638 /// Computation of covariance Matrix
2639 computeCovarianceMatrix( m, outerMatrix );
2640 memcpy( lastOuterMoments, m, nbMoments * sizeof( Quantity ));
2641 }
2642
2643 ASSERT (( lastInnerMoments[ 0 ] != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2644
2645 lastInnerSpel = currentInnerSpel;
2646 lastOuterSpel = currentOuterSpel;
2647
2648#ifdef DEBUG_VERBOSE
2649 return recount;
2650#else
2651 return false;
2652#endif
2653}