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
1067 bComputed = true;
1068 }
1069 else if( currentInnerSpel == lastOuterSpel )
1070 {
1071 m = lastOuterSum;
1072 outerSum = m;
1073
1074 bComputed = true;
1075 }
1076 else
1077 {
1078 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
1079
1080 x = diffSpel[ 0 ];
1081 y = diffSpel[ 1 ];
1082 x2 = x * x;
1083 y2 = y * y;
1084 x2y2 = x2 + y2;
1085
1086 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 );
1087
1088 if( x2y2 != 4 && x2y2 != 8 ) /// Previous and current cells aren't adjacent. Compute on the full kernel
1089 {
1090 useLastResults = false;
1091#ifdef DGTAL_DEBUG_VERBOSE
1092 recount = true;
1093#endif
1094 }
1095 else if( offset == positionOfFullKernel ) /// Full kernel in 2D. Never reached because case already considered before. It's a honeypot.
1096 {
1097 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
1098 }
1099 else
1100 {
1101 useLastResults = true;
1102 }
1103 }
1104 }
1105
1106 if( !bComputed )
1107 {
1108 if( !useLastResults ) /// Computation on full kernel, we have no previous results
1109 {
1110 m = NumberTraits< Quantity >::ZERO;
1111
1112 if( isInitFullMasks )
1113 {
1114 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
1115 {
1116 auto preShiftedSpel = KPS::sSpel( *itm );
1117 preShiftedSpel.coordinates += shiftInnerSpel;
1118
1119 if( myKSpace.sIsInside( preShiftedSpel ) )
1120 {
1121 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1122
1123 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1124 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1125
1126 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1127 {
1128 m += 1.0;
1129 }
1130 }
1131 }
1132 }
1133 else if( isInitKernelAndMasks )
1134 {
1135 Domain domain = myKernel->getDomain();
1136 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
1137 {
1138 if( myKernel->operator()( *itm ))
1139 {
1140 auto preShiftedSpel = KPS::sSpel( *itm );
1141 preShiftedSpel.coordinates += shiftInnerSpel;
1142
1143 if( myKSpace.sIsInside( preShiftedSpel ) )
1144 {
1145 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1146
1147 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1148 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1149
1150 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1151 {
1152 m += 1.0;
1153 }
1154 }
1155 }
1156 }
1157 }
1158 else
1159 {
1160 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
1161 return false;
1162 }
1163
1164 }
1165 else /// Using lastInnerMoments
1166 {
1167 m = lastInnerSum;
1168
1169 /// Part to subtract from previous result.
1170 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
1171 {
1172 auto preShiftedSpel = KPS::sSpel( *itm );
1173 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
1174
1175 if( myKSpace.sIsInside( preShiftedSpel ) )
1176 {
1177 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1178
1179 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1180 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1181
1182 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1183 {
1184 m -= 1.0;
1185 }
1186 }
1187 }
1188
1189 /// Part to add from previous result.
1190 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
1191 {
1192 auto preShiftedSpel = KPS::sSpel( *itm );
1193 preShiftedSpel.coordinates += shiftInnerSpel;
1194
1195 if( myKSpace.sIsInside( preShiftedSpel ) )
1196 {
1197 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1198
1199 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1200 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1201
1202 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1203 {
1204 m += 1.0;
1205 }
1206 }
1207 }
1208 }
1209
1210 /// Computation of covariance Matrix
1211 innerSum = m;
1212 lastInnerSum = m;
1213 }
1214 }
1215
1216 /// Outer cell
1217 {
1218 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
1219 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
1220 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
1221
1222 ASSERT( currentInnerSpel != currentOuterSpel );
1223
1224 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
1225
1226 x = diffSpel[ 0 ];
1227 y = diffSpel[ 1 ];
1228 x2 = x * x;
1229 y2 = y * y;
1230 x2y2 = x2 + y2;
1231
1232 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 );
1233
1234 if( x2y2 != 4 && x2y2 != 8 ) /// Inside and outside cells aren't adjacent, but should be. It's a honeypot.
1235 {
1236 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
1237 }
1238 else if( offset == positionOfFullKernel ) /// Full kernel in 2D. Never reached because case already considered before. It's a honeypot.
1239 {
1240 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
1241 }
1242 else
1243 {
1244 /// Part to subtract from previous result.
1245 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
1246 {
1247 auto preShiftedSpel = KPS::sSpel( *itm );
1248 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
1249
1250 if( myKSpace.sIsInside( preShiftedSpel ) )
1251 {
1252 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1253
1254 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1255 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1256
1257 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1258 {
1259 m -= 1.0;
1260 }
1261 }
1262 }
1263
1264 /// Part to add from previous result.
1265 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
1266 {
1267 auto preShiftedSpel = KPS::sSpel( *itm );
1268 preShiftedSpel.coordinates += shiftOuterSpel;
1269
1270 if( myKSpace.sIsInside( preShiftedSpel ) )
1271 {
1272 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
1273
1274 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
1275 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
1276
1277 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
1278 {
1279 m += 1.0;
1280 }
1281 }
1282 }
1283 }
1284
1285 /// Computation of covariance Matrix
1286 outerSum = m;
1287 lastOuterSum = m;
1288 }
1289
1290 ASSERT (( lastInnerSum != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
1291
1292 lastInnerSpel = currentInnerSpel;
1293 lastOuterSpel = currentOuterSpel;
1294
1295#ifdef DGTAL_DEBUG_VERBOSE
1296 return recount;
1297#else
1298 return false;
1299#endif
1300}
1301
1302
1303template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
1304template< typename SurfelIterator >
1305bool
1306DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 2 >::core_evalCovarianceMatrix
1307( const SurfelIterator & it,
1308 CovarianceMatrix & innerMatrix,
1309 CovarianceMatrix & outerMatrix,
1310 bool useLastResults,
1311 Spel & lastInnerSpel,
1312 Spel & lastOuterSpel,
1313 Quantity * lastInnerMoments,
1314 Quantity * lastOuterMoments ) const
1315{
1316 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
1317
1318 using KPS = typename KSpace::PreCellularGridSpace;
1319
1320#ifdef DGTAL_DEBUG_VERBOSE
1321 Dimension recount = false;
1322#endif
1323
1324 typedef typename Functor::Quantity FQuantity;
1325 DGtal::Dimension nbMasks = myMasks->size() - 1;
1326 DGtal::Dimension positionOfFullKernel = 4;
1327
1328 Quantity m[ nbMoments ]; /// <! [ m00, m01, m10, m11, m02, m20 ]
1329 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
1330
1331 Spel currentInnerSpel, currentOuterSpel;
1332 Spel shiftedSpel;
1333 Point shiftInnerSpel, shiftOuterSpel;
1334 Point diffSpel;
1335
1336 bool bComputed = false; /// <! if the cell has already been computed, continue to the next
1337
1338 int x, y, x2, y2, x2y2;
1339 unsigned int offset;
1340
1341 /// Inner cell
1342 {
1343 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
1344 currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
1345 shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKernelKCoordsOrigin;
1346
1347 /// Check if we can use previous results
1348 if( useLastResults )
1349 {
1350 bComputed = false;
1351
1352 if( currentInnerSpel == lastInnerSpel )
1353 {
1354 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
1355 computeCovarianceMatrix( m, innerMatrix );
1356
1357 bComputed = true;
1358 }
1359 else if( currentInnerSpel == lastOuterSpel )
1360 {
1361 memcpy( m, lastOuterMoments, nbMoments * sizeof( Quantity ));
1362 computeCovarianceMatrix( m, outerMatrix );
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 outerSum = m;
2125
2126 bComputed = true;
2127 }
2128 else
2129 {
2130 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2131
2132 x = diffSpel[ 0 ];
2133 y = diffSpel[ 1 ];
2134 z = diffSpel[ 2 ];
2135 x2 = x * x;
2136 y2 = y * y;
2137 z2 = z * z;
2138 x2y2z2 = x2 + y2 + z2;
2139
2140 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2141
2142 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2143 {
2144 useLastResults = false;
2145#ifdef DGTAL_DEBUG_VERBOSE
2146 recount = true;
2147#endif
2148 }
2149 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2150 {
2151 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2152 }
2153 else
2154 {
2155 useLastResults = true;
2156 }
2157 }
2158 }
2159
2160 if( !bComputed )
2161 {
2162 if( !useLastResults ) /// Computation on full kernel, we have no previous results
2163 {
2164 m = NumberTraits< Quantity >::ZERO;
2165
2166 if( isInitFullMasks )
2167 {
2168 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2169 {
2170 auto preShiftedSpel = KPS::sSpel( *itm );
2171 preShiftedSpel.coordinates += shiftInnerSpel;
2172
2173 if( myKSpace.sIsInside( preShiftedSpel ) )
2174 {
2175 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2176
2177 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2178 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2179 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2180
2181 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2182 {
2183 m += 1.0;
2184 }
2185 }
2186 }
2187 }
2188 else if( isInitKernelAndMasks )
2189 {
2190 Domain domain = myKernel->getDomain();
2191 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2192 {
2193 if( myKernel->operator()( *itm ))
2194 {
2195 auto preShiftedSpel = KPS::sSpel( *itm );
2196 preShiftedSpel.coordinates += shiftInnerSpel;
2197
2198 if( myKSpace.sIsInside( preShiftedSpel ) )
2199 {
2200 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2201
2202 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2203 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2204 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2205
2206 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2207 {
2208 m += 1.0;
2209 }
2210 }
2211 }
2212 }
2213 }
2214 else
2215 {
2216 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2217 return false;
2218 }
2219 }
2220 else /// Using lastInnerMoments
2221 {
2222 m = lastInnerSum;
2223
2224 /// Part to subtract from previous result.
2225 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2226 {
2227 auto preShiftedSpel = KPS::sSpel( *itm );
2228 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2229
2230 if( myKSpace.sIsInside( preShiftedSpel ) )
2231 {
2232 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2233
2234 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2235 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2236 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2237
2238 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2239 {
2240 m -= 1.0;
2241 }
2242 }
2243 }
2244
2245 /// Part to add from previous result.
2246 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2247 {
2248 auto preShiftedSpel = KPS::sSpel( *itm );
2249 preShiftedSpel.coordinates += shiftInnerSpel;
2250
2251 if( myKSpace.sIsInside( preShiftedSpel ) )
2252 {
2253 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2254
2255 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2256 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2257 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2258
2259 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2260 {
2261 m += 1.0;
2262 }
2263 }
2264 }
2265 }
2266
2267 /// Computation of covariance Matrix
2268 innerSum = m;
2269 lastInnerSum = m;
2270 }
2271 }
2272
2273 /// Outer cell
2274 {
2275 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2276 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2277 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
2278
2279 ASSERT( currentInnerSpel != currentOuterSpel );
2280
2281 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2282
2283 x = diffSpel[ 0 ];
2284 y = diffSpel[ 1 ];
2285 z = diffSpel[ 2 ];
2286 x2 = x * x;
2287 y2 = y * y;
2288 z2 = z * z;
2289 x2y2z2 = x2 + y2 + z2;
2290
2291 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2292
2293 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.
2294 {
2295 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2296 }
2297 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2298 {
2299 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2300 }
2301 else
2302 {
2303 /// Part to subtract from previous result.
2304 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2305 {
2306 auto preShiftedSpel = KPS::sSpel( *itm );
2307 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2308
2309 if( myKSpace.sIsInside( preShiftedSpel ) )
2310 {
2311 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates
2312 );
2313
2314 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2315 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2316 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2317
2318 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2319 {
2320 m -= 1.0;
2321 }
2322 }
2323 }
2324
2325 /// Part to add from previous result.
2326 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2327 {
2328 auto preShiftedSpel = KPS::sSpel( *itm );
2329 preShiftedSpel.coordinates += shiftOuterSpel;
2330
2331 if( myKSpace.sIsInside( preShiftedSpel ) )
2332 {
2333 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2334
2335 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2336 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2337 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2338
2339 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2340 {
2341 m += 1.0;
2342 }
2343 }
2344 }
2345 }
2346
2347 /// Computation of covariance Matrix
2348 outerSum = m;
2349 lastOuterSum = m;
2350 }
2351
2352 ASSERT (( lastInnerSum != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2353
2354 lastInnerSpel = currentInnerSpel;
2355 lastOuterSpel = currentOuterSpel;
2356
2357#ifdef DGTAL_DEBUG_VERBOSE
2358 return recount;
2359#else
2360 return false;
2361#endif
2362}
2363
2364
2365template< typename Functor, typename KernelFunctor, typename KSpace, typename DigitalKernel >
2366template< typename SurfelIterator >
2367bool
2368DGtal::DigitalSurfaceConvolver< Functor, KernelFunctor, KSpace, DigitalKernel, 3 >::core_evalCovarianceMatrix
2369( const SurfelIterator & it,
2370 CovarianceMatrix & innerMatrix,
2371 CovarianceMatrix & outerMatrix,
2372 bool useLastResults,
2373 Spel & lastInnerSpel,
2374 Spel & lastOuterSpel,
2375 Quantity * lastInnerMoments,
2376 Quantity * lastOuterMoments ) const
2377{
2378 ASSERT ( isInitFullMasks == true || isInitKernelAndMasks == true );
2379
2380 using KPS = typename KSpace::PreCellularGridSpace;
2381
2382#ifdef DGTAL_DEBUG_VERBOSE
2383 Dimension recount = false;
2384#endif
2385
2386 typedef typename Functor::Quantity FQuantity;
2387 DGtal::Dimension nbMasks = static_cast<Dimension>(myMasks->size()) - 1;
2388 DGtal::Dimension positionOfFullKernel = 13;
2389
2390 Quantity m[ nbMoments ]; /// <! [ m000, m001, m010, m100, m011, m101, m110, m002, m020, m200 ]
2391 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2392
2393 Spel currentInnerSpel, currentOuterSpel;
2394 Spel shiftedSpel;
2395 Point shiftInnerSpel, shiftOuterSpel;
2396 Point diffSpel;
2397
2398 bool bComputed = false; /// <! if the cell has already been computed, continue to the next
2399
2400 int x, y, z, x2, y2, z2, x2y2z2;
2401 unsigned int offset;
2402
2403 /// Inner cell
2404 {
2405 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2406 currentInnerSpel = myKSpace.sDirectIncident( *it, kDim ); /// Spel on the border, but inside the shape
2407 shiftInnerSpel = myKSpace.sKCoords( currentInnerSpel ) - myKernelKCoordsOrigin;
2408
2409 /// Check if we can use previous results
2410 if( useLastResults )
2411 {
2412 bComputed = false;
2413
2414 if( currentInnerSpel == lastInnerSpel )
2415 {
2416 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2417 computeCovarianceMatrix( m, innerMatrix );
2418
2419 bComputed = true;
2420 }
2421 else if( currentInnerSpel == lastOuterSpel )
2422 {
2423 memcpy( m, lastOuterMoments, nbMoments * sizeof( Quantity ));
2424 computeCovarianceMatrix( m, outerMatrix );
2425
2426 bComputed = true;
2427 }
2428 else
2429 {
2430 diffSpel = myKSpace.sKCoords( currentInnerSpel ) - myKSpace.sKCoords( lastInnerSpel );
2431
2432 x = diffSpel[ 0 ];
2433 y = diffSpel[ 1 ];
2434 z = diffSpel[ 2 ];
2435 x2 = x * x;
2436 y2 = y * y;
2437 z2 = z * z;
2438 x2y2z2 = x2 + y2 + z2;
2439
2440 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1 ) * 9 );
2441
2442 if( x2y2z2 != 4 && x2y2z2 != 8 && !( x2 == 4 && y2 == 4 && z2 == 4 )) /// Previous and current cells aren't adjacent. Compute on the full kernel
2443 {
2444 useLastResults = false;
2445#ifdef DGTAL_DEBUG_VERBOSE
2446 recount = true;
2447#endif
2448 }
2449 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2450 {
2451 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2452 }
2453 else
2454 {
2455 useLastResults = true;
2456 }
2457 }
2458 }
2459
2460 if( !bComputed )
2461 {
2462 if( !useLastResults ) /// Computation on full kernel, we have no previous results
2463 {
2464 std::fill( m, m + nbMoments, NumberTraits< Quantity >::ZERO ); /// <! clear array
2465
2466 if( isInitFullMasks )
2467 {
2468 for( KernelConstIterator itm = myKernelMask->first; itm != myKernelMask->second; ++itm )
2469 {
2470 auto preShiftedSpel = KPS::sSpel( *itm );
2471 preShiftedSpel.coordinates += shiftInnerSpel;
2472
2473 if( myKSpace.sIsInside( preShiftedSpel ) )
2474 {
2475 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2476
2477 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2478 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2479 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2480
2481 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2482 {
2483 fillMoments( m, shiftedSpel, 1.0 );
2484 }
2485 }
2486 }
2487 }
2488 else if( isInitKernelAndMasks )
2489 {
2490 Domain domain = myKernel->getDomain();
2491 for( typename Domain::ConstIterator itm = domain.begin(), itend = domain.end(); itm != itend; ++itm )
2492 {
2493 if( myKernel->operator()( *itm ))
2494 {
2495 auto preShiftedSpel = KPS::sSpel( *itm );
2496 preShiftedSpel.coordinates += shiftInnerSpel;
2497
2498 if( myKSpace.sIsInside( preShiftedSpel ) )
2499 {
2500 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2501
2502 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2503 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2504 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2505
2506 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2507 {
2508 fillMoments( m, shiftedSpel, 1.0 );
2509 }
2510 }
2511 }
2512 }
2513 }
2514 else
2515 {
2516 trace.error() << "DigitalSurfaceConvolver: You need to init the convolver first." << std::endl;
2517 return false;
2518 }
2519 }
2520 else /// Using lastInnerMoments
2521 {
2522 memcpy( m, lastInnerMoments, nbMoments * sizeof( Quantity ));
2523
2524 /// Part to subtract from previous result.
2525 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2526 {
2527 auto preShiftedSpel = KPS::sSpel( *itm );
2528 preShiftedSpel.coordinates += shiftInnerSpel - diffSpel;
2529
2530 if( myKSpace.sIsInside( preShiftedSpel ) )
2531 {
2532 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2533
2534 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2535 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2536 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2537
2538 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2539 {
2540 fillMoments( m, shiftedSpel, -1.0 );
2541 }
2542 }
2543 }
2544
2545 /// Part to add from previous result.
2546 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2547 {
2548 auto preShiftedSpel = KPS::sSpel( *itm );
2549 preShiftedSpel.coordinates += shiftInnerSpel;
2550
2551 if( myKSpace.sIsInside( preShiftedSpel ) )
2552 {
2553 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2554
2555 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2556 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2557 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2558
2559 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2560 {
2561 fillMoments( m, shiftedSpel, 1.0 );
2562 }
2563 }
2564 }
2565 }
2566
2567 /// Computation of covariance Matrix
2568 computeCovarianceMatrix( m, innerMatrix );
2569 memcpy( lastInnerMoments, m, nbMoments * sizeof( Quantity ));
2570 }
2571 }
2572
2573 /// Outer cell
2574 {
2575 DGtal::Dimension kDim = myKSpace.sOrthDir( *it );
2576 currentOuterSpel = myKSpace.sIndirectIncident( *it, kDim );
2577 shiftOuterSpel = myKSpace.sKCoords( currentOuterSpel ) - myKernelKCoordsOrigin;
2578
2579 ASSERT( currentInnerSpel != currentOuterSpel );
2580
2581 diffSpel = myKSpace.sKCoords( currentOuterSpel ) - myKSpace.sKCoords( currentInnerSpel );
2582
2583 x = diffSpel[ 0 ];
2584 y = diffSpel[ 1 ];
2585 z = diffSpel[ 2 ];
2586 x2 = x * x;
2587 y2 = y * y;
2588 z2 = z * z;
2589 x2y2z2 = x2 + y2 + z2;
2590
2591 offset = (( x / 2 ) + 1 ) + ((( y / 2 ) + 1 ) * 3 ) + ((( z / 2 ) + 1) * 9 );
2592
2593 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.
2594 {
2595 trace.error() << "Error - Found that inside and outside cells aren't adjacent. That's not logic.\n";
2596 }
2597 else if( offset == positionOfFullKernel ) /// Full kernel in 3D. Never reached because case already considered before. It's a honeypot.
2598 {
2599 trace.error() << "Error - Ask to compute full kernel at the wrong moment. Maybe it's an offset computation bug ?\n";
2600 }
2601 else
2602 {
2603 /// Part to subtract from previous result.
2604 for( KernelConstIterator itm = (*myMasks)[ offset ].first, itend = (*myMasks)[ offset ].second; itm != itend; ++itm )
2605 {
2606 auto preShiftedSpel = KPS::sSpel( *itm );
2607 preShiftedSpel.coordinates += shiftOuterSpel - diffSpel;
2608
2609 if( myKSpace.sIsInside( preShiftedSpel ) )
2610 {
2611 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2612
2613 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2614 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2615 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2616
2617 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2618 {
2619 fillMoments( m, shiftedSpel, -1.0 );
2620 }
2621 }
2622 }
2623
2624 /// Part to add from previous result.
2625 for( KernelConstIterator itm = (*myMasks)[ nbMasks - offset ].first, itend = (*myMasks)[ nbMasks - offset ].second; itm != itend; ++itm )
2626 {
2627 auto preShiftedSpel = KPS::sSpel( *itm );
2628 preShiftedSpel.coordinates += shiftOuterSpel;
2629
2630 if( myKSpace.sIsInside( preShiftedSpel ) )
2631 {
2632 myKSpace.sSetKCoords( shiftedSpel, preShiftedSpel.coordinates );
2633
2634 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 0 ] & 0x1 );
2635 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 1 ] & 0x1 );
2636 ASSERT( myKSpace.sKCoords( shiftedSpel )[ 2 ] & 0x1 );
2637
2638 if( myFFunctor( shiftedSpel ) != NumberTraits< FQuantity >::ZERO )
2639 {
2640 fillMoments( m, shiftedSpel, 1.0 );
2641 }
2642 }
2643 }
2644 }
2645
2646 /// Computation of covariance Matrix
2647 computeCovarianceMatrix( m, outerMatrix );
2648 memcpy( lastOuterMoments, m, nbMoments * sizeof( Quantity ));
2649 }
2650
2651 ASSERT (( lastInnerMoments[ 0 ] != 0 )); // Maybe a problem here. Can be normal ... but needs to check twice.
2652
2653 lastInnerSpel = currentInnerSpel;
2654 lastOuterSpel = currentOuterSpel;
2655
2656#ifdef DGTAL_DEBUG_VERBOSE
2657 return recount;
2658#else
2659 return false;
2660#endif
2661}