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