aimsalgo  5.1.2
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 
43 template <class T>
44 class AimsVFilter
45 {
46 public:
47 
48  enum VFilterType
49  {
50  NonOptimized,
51  Optimized
52  };
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 
59 private:
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 
69 template< class T > inline
70 carto::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!
89 template< class T > inline
90 carto::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!
274 template< class T > inline
275 carto::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
The template class to perform a V-filter.
Definition: vfilter.h:45
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
AimsVFilter::Optimized
@ Optimized
Definition: vfilter.h:51
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)