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