aimsalgo 6.0.0
Neuroimaging image processing
vfilter.h
Go to the documentation of this file.
1/* This software and supporting documentation are distributed by
2 * Institut Federatif de Recherche 49
3 * CEA/NeuroSpin, Batiment 145,
4 * 91191 Gif-sur-Yvette cedex
5 * France
6 *
7 * This software is governed by the CeCILL-B license under
8 * French law and abiding by the rules of distribution of free software.
9 * You can use, modify and/or redistribute the software under the
10 * terms of the CeCILL-B license as circulated by CEA, CNRS
11 * and INRIA at the following URL "http://www.cecill.info".
12 *
13 * As a counterpart to the access to the source code and rights to copy,
14 * modify and redistribute granted by the license, users are provided only
15 * with a limited warranty and the software's author, the holder of the
16 * economic rights, and the successive licensors have only limited
17 * liability.
18 *
19 * In this respect, the user's attention is drawn to the risks associated
20 * with loading, using, modifying and/or developing or reproducing the
21 * software by the user in light of its specific status of free software,
22 * that may mean that it is complicated to manipulate, and that also
23 * therefore means that it is reserved for developers and experienced
24 * professionals having in-depth computer knowledge. Users are therefore
25 * encouraged to load and test the software's suitability as regards their
26 * requirements in conditions enabling the security of their systems and/or
27 * data to be ensured and, more generally, to use and operate it in the
28 * same conditions as regards security.
29 *
30 * The fact that you are presently reading this means that you have had
31 * knowledge of the CeCILL-B license and that you accept its terms.
32 */
33
34
35#ifndef AIMS_SIGNALFILTER_VFILTER_H
36#define AIMS_SIGNALFILTER_VFILTER_H
37
38#include <cartodata/volume/volume.h>
39
40
43template <class T>
44class AimsVFilter
45{
46public:
47
53
54 AimsVFilter( int s=2, VFilterType t=Optimized ) : _size( s ), _type( t ) { }
55 virtual ~AimsVFilter() { }
56
57 carto::VolumeRef< T > doit( carto::VolumeRef< T >& );
58
59private:
60
61 carto::VolumeRef< T > optimized( carto::VolumeRef< T >& );
62 carto::VolumeRef< T > nonOptimized( carto::VolumeRef< T >& );
63
64 int _size;
65 VFilterType _type;
66};
67
68
69template< class T > inline
70carto::VolumeRef< T > AimsVFilter< T >::doit( carto::VolumeRef< T >& d )
71{
72 carto::VolumeRef< T > res;
73
74 switch( _type )
75 {
76 case NonOptimized:
77 res = nonOptimized( d );
78 break;
79 case Optimized:
80 res = optimized( d );
81 break;
82 }
83
84 return res;
85}
86
87
88// This is the non-optimized version!
89template< class T > inline
90carto::VolumeRef< T > AimsVFilter< T >::nonOptimized( carto::VolumeRef< T >& d )
91{
92 int x, y, z, t, i, j, k, index;
93 double sum, sum2, varmin, var[8], moy[8];
94 T temp;
95
96 int tm1 = _size - 1;
97 double t3 = (double)( _size * _size * _size );
98
99 int dx = d->getSizeX();
100 int dy = d->getSizeY();
101 int dz = d->getSizeZ();
102 int dt = d->getSizeT();
103
104 carto::VolumeRef< T > res( dx, dy, dz, dt );
105 res.setVoxelSize( d->getVoxelSize() );
106
107 for ( t=0; t<dt; t++ )
108 for ( z=0; z<dz; z++ )
109 for ( y=0; y<dy; y++ )
110 for ( x=0; x<dx; x++ )
111 {
112 // Initialization
113 for ( i=0; i<8; i++ )
114 {
115 var[i] = 1.0e20;
116 moy[i] = 0.0;
117 }
118
119 // First cube (-1,-1,-1)
120 if ( ( x >= tm1 ) && ( y >= tm1 ) && ( z >= tm1 ) )
121 {
122 sum = sum2 = 0.0;
123 for ( k=z-tm1; k<z+1; k++ )
124 for ( j=y-tm1; j<y+1; j++ )
125 for ( i=x-tm1; i<x+1; i++ )
126 {
127 temp = d( i, j, k, t );
128 sum += (double)temp;
129 sum2 += (double)( temp * temp );
130 }
131
132 var[ 0 ] = sum2 - ( sum * sum ) / t3;
133 moy[ 0 ] = sum / t3;
134 }
135
136 // Second cube (-1,1,-1)
137 if ( ( x >= tm1 ) && ( y < dy-tm1 ) && ( z >= tm1 ) )
138 {
139 sum = sum2 = 0.0;
140 for ( k=z-tm1; k<z+1; k++ )
141 for ( j=y; j<y+_size; j++ )
142 for ( i=x-tm1; i<x+1; i++ )
143 {
144 temp = d( i, j, k, t );
145 sum += (double)temp;
146 sum2 += (double)( temp * temp );
147 }
148
149 var[ 1 ] = sum2 - ( sum * sum ) / t3;
150 moy[ 1 ] = sum / t3;
151 }
152
153 // Third cube (1,1,-1)
154 if ( ( x < dx-tm1 ) && ( y < dy-tm1 ) && ( z >= tm1 ) )
155 {
156 sum = sum2 = 0.0;
157 for ( k=z-tm1; k<z+1; k++ )
158 for ( j=y; j<y+_size; j++ )
159 for ( i=x; i<x+_size; i++ )
160 {
161 temp = d( i, j, k, t );
162 sum += (double)temp;
163 sum2 += (double)( temp * temp );
164 }
165
166 var[ 2 ] = sum2 - ( sum * sum ) / t3;
167 moy[ 2 ] = sum / t3;
168 }
169
170 // Fourth cube (1,-1,-1)
171 if ( ( x < dx-tm1 ) && ( y >= tm1 ) && ( z >= tm1 ) )
172 {
173 sum = sum2 = 0.0;
174 for ( k=z-tm1; k<z+1; k++ )
175 for ( j=y-tm1; j<y+1; j++ )
176 for ( i=x; i<x+_size; i++ )
177 {
178 temp = d( i, j, k, t );
179 sum += (double)temp;
180 sum2 += (double)( temp * temp );
181 }
182
183 var[ 3 ] = sum2 - ( sum * sum ) / t3;
184 moy[ 3 ] = sum / t3;
185 }
186
187 // Fifth cube (-1,-1,1)
188 if ( ( x >= tm1 ) && ( y >= tm1 ) && ( z < dz-tm1 ) )
189 {
190 sum = sum2 = 0.0;
191 for ( k=z; k<z+_size; k++ )
192 for ( j=y-tm1; j<y+1; j++ )
193 for ( i=x-tm1; i<x+1; i++ )
194 {
195 temp = d( i, j, k, t );
196 sum += (double)temp;
197 sum2 += (double)( temp * temp );
198 }
199
200 var[ 4 ] = sum2 - ( sum * sum ) / t3;
201 moy[ 4 ] = sum / t3;
202 }
203
204 // Sixth cube (-1,1,1)
205 if ( ( x >= tm1 ) && ( y < dy-tm1 ) && ( z < dz-tm1 ) )
206 {
207 sum = sum2 = 0.0;
208 for ( k=z; k<z+_size; k++ )
209 for ( j=y; j<y+_size; j++ )
210 for ( i=x-tm1; i<x+1; i++ )
211 {
212 temp = d( i, j, k, t );
213 sum += (double)temp;
214 sum2 += (double)( temp * temp );
215 }
216
217 var[ 5 ] = sum2 - ( sum * sum ) / t3;
218 moy[ 5 ] = sum / t3;
219 }
220
221 // Seventh cube (1,1,1)
222 if ( ( x < dx-tm1 ) && ( y < dy-tm1 ) && ( z < dz-tm1 ) )
223 {
224 sum = sum2 = 0.0;
225 for ( k=z; k<z+_size; k++ )
226 for ( j=y; j<y+_size; j++ )
227 for ( i=x; i<x+_size; i++ )
228 {
229 temp = d( i, j, k, t );
230 sum += (double)temp;
231 sum2 += (double)( temp * temp );
232 }
233
234 var[ 6 ] = sum2 - ( sum * sum ) / t3;
235 moy[ 6 ] = sum / t3;
236 }
237
238 // Eighth cube (1,-1,1)
239 if ( ( x < dx-tm1 ) && ( y >= tm1 ) && ( z < dz-tm1 ) )
240 {
241 sum = sum2 = 0.0;
242 for ( k=z; k<z+_size; k++ )
243 for ( j=y-tm1; j<y+1; j++ )
244 for ( i=x; i<x+_size; i++ )
245 {
246 temp = d( i, j, k, t );
247 sum += (double)temp;
248 sum2 += (double)( temp * temp );
249 }
250
251 var[ 7 ] = sum2 - ( sum * sum ) / t3;
252 moy[ 7 ] = sum / t3;
253 }
254
255 /* Recherche du meilleur cube */
256 index = 0;
257 varmin = var[ 0 ];
258 for ( i=1; i<8; i++ )
259 if ( var[ i ] < varmin )
260 {
261 varmin = var[ i ];
262 index = i;
263 }
264
265 /* Affectation du point */
266 res( x, y, z, t ) = (T)moy[ index ];
267 }
268
269 return res;
270}
271
272
273// This is the optimized version!
274template< class T > inline
275carto::VolumeRef< T > AimsVFilter< T >::optimized( carto::VolumeRef< T >& d )
276{
277 int i, j, k, x, y, z, loc;
278 double somme, carre, somme1, carre1, vars, *pvar, moyenne, temp;
279 int dx = d->getSizeX();
280 int dy = d->getSizeY();
281 int dz = d->getSizeZ();
282 int tm1 = _size - 1;
283 double size2 = (double)( _size * _size * _size );
284
285 carto::VolumeRef< T > vf( dx, dy, dz );
286 vf.setVoxelSize( d->getVoxelSize() );
287
288 carto::VolumeRef< double > table( dx, dy, _size );
289
290 for ( z=_size; z--; )
291 for ( y=dy; y--; )
292 for ( x=dx; x--; ) table( x, y, z ) = -1.0;
293
294 for (z=0; z+tm1<dz; z++ )
295 {
296 i = ( z + tm1 ) % _size;
297
298 // Current table slice initialization
299 for ( y=0; y<dy; y++ )
300 for ( x=0; x<dx; x++ ) table( x, y, i ) = -1.0;
301
302 // Compute the first cube for each slice
303 somme = carre = 0.0;
304 for ( k=z; k<z+_size; k++ )
305 for ( j=0; j<_size; j++ )
306 for ( i=0; i<_size; i++ )
307 {
308 temp = (double)d( i, j, k );
309 somme += temp;
310 carre += temp * temp;
311 }
312
313 vars = carre - somme * somme / size2;
314 moyenne = somme / size2;
315 somme1 = somme;
316 carre1 = carre;
317
318 pvar = &table( 0, 0, z % _size );
319 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
320 {
321 vf( 0, 0, z ) = (T)moyenne;
322 *pvar = vars;
323 }
324 pvar = &table( 0, tm1, z % _size );
325 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
326 {
327 vf( 0, tm1, z ) = (T)moyenne;
328 *pvar = vars;
329 }
330 pvar = &table( tm1, 0, z % _size );
331 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
332 {
333 vf( tm1, 0, z ) = (T)moyenne;
334 *pvar = vars;
335 }
336 pvar = &table( tm1, tm1, z % _size );
337 if ( (vars < *pvar ) || ( *pvar == -1.0 ) )
338 {
339 vf( tm1, tm1, z ) = (T)moyenne;
340 *pvar = vars;
341 }
342 pvar = &table( 0, 0, ( z + tm1 ) % _size );
343 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
344 {
345 vf( 0, 0, z + tm1 ) = (T)moyenne;
346 *pvar = vars;
347 }
348 pvar = &table( 0, tm1, ( z + tm1 ) % _size );
349 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
350 {
351 vf( 0, tm1, z + tm1 ) = (T)moyenne;
352 *pvar = vars;
353 }
354 pvar = &table( tm1, 0, ( z + tm1 ) % _size );
355 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
356 {
357 vf( tm1, 0, z + tm1 ) = (T)moyenne;
358 *pvar = vars;
359 }
360 vf( tm1, tm1, z + tm1 ) = (T)moyenne;
361 table( tm1, tm1, ( z + tm1 ) % _size ) = vars;
362
363 // First X scan
364 for ( x=1; x+tm1<dx; x++ )
365 {
366 for ( k=z; k<z+_size; k++ )
367 for ( j=0; j<_size; j++ )
368 {
369 loc = x-1;
370 temp = (double)d( loc, j, k );
371 somme -= temp;
372 carre -= temp * temp;
373 loc = x + tm1;
374 temp = (double)d( loc, j, k );
375 somme += temp;
376 carre += temp * temp;
377 }
378 vars = carre - somme * somme / size2;
379 moyenne = somme / size2;
380
381 pvar = &table( x, 0, z % _size );
382 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
383 {
384 vf( x, 0, z ) = (T)moyenne;
385 *pvar = vars;
386 }
387 pvar = &table( x, tm1, z % _size );
388 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
389 {
390 vf( x, tm1, z ) = (T)moyenne;
391 *pvar = vars;
392 }
393 pvar = &table( x + tm1, 0, z % _size );
394 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
395 {
396 vf( x + tm1, 0, z ) = (T)moyenne;
397 *pvar = vars;
398 }
399 pvar = &table( x + tm1, tm1, z % _size );
400 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
401 {
402 vf( x + tm1, tm1, z ) = (T)moyenne;
403 *pvar = vars;
404 }
405 pvar = &table( x, 0, ( z + tm1 ) % _size );
406 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
407 {
408 vf( x, 0, z + tm1 ) = (T)moyenne;
409 *pvar = vars;
410 }
411 pvar = &table( x, tm1, ( z + tm1 ) % _size );
412 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
413 {
414 vf( x, tm1, z + tm1 ) = (T)moyenne;
415 *pvar = vars;
416 }
417 pvar = &table( x + tm1, 0, ( z + tm1 ) % _size );
418 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
419 {
420 vf( x + tm1, 0, z + tm1 ) = (T)moyenne;
421 *pvar = vars;
422 }
423 vf( x + tm1, tm1, z + tm1 ) = (T)moyenne;
424 table( x + tm1, tm1, ( z + tm1 ) % _size ) = vars;
425 }
426
427 // Computation loop
428 for ( y=1; y+tm1<dy; y++ )
429 {
430 somme = somme1;
431 carre = carre1;
432
433 /*** Compute the new cube for each Y ***/
434 for ( k=z; k<z+_size; k++)
435 for ( i=0; i<_size; i++)
436 {
437 loc = y - 1;
438 temp = (double)d( i, loc, k );
439 somme -= temp;
440 carre -= temp * temp;
441 loc = y + tm1;
442 temp = (double)d( i, loc, k );
443 somme += temp;
444 carre += temp * temp;
445 }
446
447 vars = carre - somme * somme / size2;
448 moyenne = somme / size2;
449 somme1 = somme;
450 carre1 = carre;
451
452 pvar = &table( 0, y, z % _size );
453 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
454 {
455 vf( 0, y, z ) = (T)moyenne;
456 *pvar = vars;
457 }
458 pvar = &table( 0, y + tm1, z % _size );
459 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
460 {
461 vf( 0, y + tm1, z ) = (T)moyenne;
462 *pvar = vars;
463 }
464 pvar = &table( tm1, y, z % _size );
465 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
466 {
467 vf( tm1, y, z ) = (T)moyenne;
468 *pvar = vars;
469 }
470 pvar = &table( tm1, y + tm1, z % _size );
471 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
472 {
473 vf( tm1, y + tm1, z ) = (T)moyenne;
474 *pvar = vars;
475 }
476 pvar = &table( 0, y, ( z + tm1 ) % _size );
477 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
478 {
479 vf( 0, y, z + tm1 ) = (T)moyenne;
480 *pvar = vars;
481 }
482 pvar = &table( 0, y + tm1, ( z + tm1 ) % _size );
483 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
484 {
485 vf( 0, y + tm1, z + tm1 ) = (T)moyenne;
486 *pvar = vars;
487 }
488 pvar = &table( tm1, y, ( z + tm1 ) % _size );
489 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
490 {
491 vf( tm1, y, z + tm1 ) = (T)moyenne;
492 *pvar = vars;
493 }
494 vf( tm1, y + tm1, z + tm1 ) = (T)moyenne;
495 table( tm1, y + tm1, ( z + tm1 ) % _size ) = vars;
496
497 // Main X scan
498 for ( x=1; x+tm1<dx; x++ )
499 {
500 for ( k=z; k<z+_size; k++ )
501 for ( j=y; j<y+_size; j++ )
502 {
503 loc = x - 1;
504 temp = (double)d( loc, j, k);
505 somme -= temp;
506 carre -= temp * temp;
507 loc = x + tm1;
508 temp = (double)d( loc, j, k);
509 somme += temp;
510 carre += temp * temp;
511 }
512
513 vars = carre - somme * somme / size2;
514 moyenne = somme / size2;
515
516 pvar = &table( x, y, z % _size );
517 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
518 {
519 vf( x, y, z ) = (T)moyenne;
520 *pvar = vars;
521 }
522 pvar = &table( x, y + tm1, z % _size );
523 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
524 {
525 vf( x, y + tm1, z ) = (T)moyenne;
526 *pvar = vars;
527 }
528 pvar = &table( x + tm1, y, z % _size );
529 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
530 {
531 vf( x + tm1, y, z ) = (T)moyenne;
532 *pvar = vars;
533 }
534 pvar = &table( x + tm1, y + tm1, z % _size );
535 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
536 {
537 vf( x + tm1, y + tm1, z ) = (T)moyenne;
538 *pvar = vars;
539 }
540 pvar = &table( x, y, ( z + tm1 ) % _size );
541 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
542 {
543 vf( x, y, z + tm1 ) = (T)moyenne;
544 *pvar = vars;
545 }
546 pvar = &table( x, y + tm1, ( z + tm1 ) % _size );
547 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
548 {
549 vf( x, y + tm1, z + tm1 ) = (T)moyenne;
550 *pvar = vars;
551 }
552 pvar = &table( x + tm1, y, ( z + tm1 ) % _size );
553 if ( ( vars < *pvar ) || ( *pvar == -1.0 ) )
554 {
555 vf( x + tm1, y, z + tm1 ) = (T)moyenne;
556 *pvar = vars;
557 }
558 vf( x + tm1, y + tm1, z + tm1 ) = (T)moyenne;
559 table( x + tm1, y + tm1, ( z + tm1 ) % _size ) = vars;
560 }
561 }
562 }
563
564 return vf;
565}
566
567#endif
AimsVFilter::~AimsVFilter
virtual ~AimsVFilter()
Definition vfilter.h:55
AimsVFilter::doit
carto::VolumeRef< T > doit(carto::VolumeRef< T > &)
Definition vfilter.h:70
AimsVFilter::AimsVFilter
AimsVFilter(int s=2, VFilterType t=Optimized)
Definition vfilter.h:54
AimsVFilter::NonOptimized
@ NonOptimized
Definition vfilter.h:50
carto::VolumeRef::getSizeZ
int getSizeZ() const
carto::VolumeRef::getVoxelSize
std::vector< float > getVoxelSize() const
carto::VolumeRef::getSizeT
int getSizeT() const
carto::VolumeRef::getSizeY
int getSizeY() const
carto::VolumeRef::getSizeX
int getSizeX() const
sum
DataTypeTraits< T >::LongType sum(const Volume< T > &vol)