A.I.M.S algorithms


diffusionSmoother_d.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 
36 #ifndef AIMS_DIFFUSION_SMOOTHER_D_H
37 #define AIMS_DIFFUSION_SMOOTHER_D_H
38 
39 #include <cstdlib>
43 #include <aims/data/volumemanip.h>
44 #include <stdlib.h>
45 #include <math.h>
46 #include <float.h>
47 
48 namespace aims {
49 
50 
51 template<typename T>
52 AimsData<float> BaseDiffusionSmoother<T>::laplacian =
54 
55 template<typename T>
57 {
58  AimsData<float> laplacian(3, 3, 3, 1);
59  laplacian = 0.;
60  /* WARNING: this kernel is only suitable for isotropic, 1mm voxels.
61  Moreover I don't explain the 1/2 factor on the whole kernel
62  (I would expect coefs 1 -6 1 instead of 0.5 -3 0.5)
63  Anyway, this kernel is now completely rebuilt in doSmoothing() to
64  take into account voxel size, so it is not used as is.
65  (Denis, 2014/02/14)
66  */
67  laplacian(1,1,0) = laplacian(1,0,1) = laplacian(0,1,1) =
68  laplacian(2,1,1) = laplacian(1,2,1) = laplacian(1,1,2) = 0.5;
69  laplacian(1,1,1) = -3;
70 
71  return laplacian;
72 }
73 
74 template<class T>
76 {
77  if (maxiter < 0)
78  {
79  std::cerr << "diffusionConvolution Smoother: "
80  << "maxiter negative value is forbidden" << std::endl;
81  exit(EXIT_FAILURE);
82  }
83 }
84 
85 template<class T>
87  int maxiter, bool verbose)
88 {
89  this->check(maxiter);
91  // rebuild the kernel accounting for possibly anisotropic voxel size
92  kernel(0,1,1) = kernel(2,1,1) = this->_dt / ( ima.sizeX() * ima.sizeX() );
93  kernel(1,0,1) = kernel(1,2,1) = this->_dt / ( ima.sizeY() * ima.sizeY() );
94  kernel(1,1,0) = kernel(1,1,2) = this->_dt / ( ima.sizeZ() * ima.sizeZ() );
95  kernel(1,1,1) = -this->_dt * (2. / (ima.sizeX() * ima.sizeX())
96  + 2. / (ima.sizeY() * ima.sizeY())
97  + 2. / (ima.sizeZ() * ima.sizeZ()) );
98  //
99  // 1 (center) + kernel * dt
100  kernel(1,1,1) = 1. + kernel(1,1,1);
101  /*
102  kernel(1,1,1) = 1. + kernel(1,1,1) * this->_dt;
103  kernel(1,1,0) = kernel(1,0,1) = kernel(0,1,1) =
104  kernel(2,1,1) = kernel(1,2,1) = kernel(1,1,2) = 0.5 * this->_dt;
105  */
106 
107  AimsData<float> *tmp1, *tmp2, *swap;
108  AimsConvolution<float> convolution;
109 
111  AimsData< float > imaF(ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT()),
112  ima2(ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT());
113  conv.convert(ima, imaF);
114  tmp1 = &imaF; tmp2 = &ima2;
115 
116  for (int i = 0; i < maxiter; i++)
117  {
118  if (verbose) std::cout << "+" << std::flush;
119  if( _hasConstantSources )
120  {
121  int x, y, z;
122  ForEach3d( (*tmp1), x, y, z )
123  {
124  T & v = _constantSources( x, y, z );
125  if( v != _constantSourcesBackground )
126  (*tmp1)( x, y, z ) = v;
127  }
128  }
129  (*tmp2) = convolution.doit((*tmp1), kernel);
130  swap = tmp1;
131  tmp1 = tmp2;
132  tmp2 = swap;
133  }
134 
135  if( _hasConstantSources )
136  {
137  int x, y, z;
138  ForEach3d( (*tmp1), x, y, z )
139  {
140  T & v = _constantSources( x, y, z );
141  if( v != _constantSourcesBackground )
142  (*tmp1)( x, y, z ) = v;
143  }
144  }
145 
146  if (verbose) std::cout << "Finished" << std::endl;
148  AimsData<T> ima3( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
149  conv2.convert( (*tmp1), ima3);
150 
151  return ima3;
152 }
153 
154 
155 template<typename T>
157  const T & bgr )
158 {
159  _constantSources = sources;
160  _constantSourcesBackground = bgr;
161  _hasConstantSources = true;
162 }
163 
164 
165 template <typename T>
167 {
168  _hasConstantSources = false;
169  _constantSources = AimsData<T>();
170 }
171 
172 
173 template<typename T> AimsData<T>
175 doSmoothing(const AimsData<T> & ima, int maxiter, bool verbose)
176 {
177  this->check(maxiter);
179  // rebuild the kernel accounting for possibly anisotropic voxel size
180  kernel(0,1,1) = kernel(2,1,1) = this->_dt / ( ima.sizeX() * ima.sizeX() );
181  kernel(1,0,1) = kernel(1,2,1) = this->_dt / ( ima.sizeY() * ima.sizeY() );
182  kernel(1,1,0) = kernel(1,1,2) = this->_dt / ( ima.sizeZ() * ima.sizeZ() );
183  kernel(1,1,1) = -this->_dt * (2. / (ima.sizeX() * ima.sizeX())
184  + 2. / (ima.sizeY() * ima.sizeY())
185  + 2. / (ima.sizeZ() * ima.sizeZ()) );
186  //
187  // 1 (center) + kernel * dt
188  kernel(1,1,1) = 1. + kernel(1,1,1);
189  /*
190  kernel(1,1,1) = 1. + kernel(1,1,1) * this->_dt;
191  kernel(1,1,0) = kernel(1,0,1) = kernel(0,1,1) =
192  kernel(2,1,1) = kernel(1,2,1) = kernel(1,1,2) = 0.5 * this->_dt;
193  */
194 
195  AimsData<float> *tmp1, *tmp2, *swap;
196  AimsMaskedConvolution<float> convolution(*(this->_mask),
197  this->_background, this->_safe);
198 
200  AimsData< float > imaF(ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT()),
201  ima2(ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT());
202  conv.convert(ima, imaF);
203  tmp1 = &imaF; tmp2 = &ima2;
204 
205  for (int i = 0; i < maxiter; i++)
206  {
207  if (verbose) std::cout << "+" << std::flush;
208  if (_has_neumann_condition)
209  update_neumann_conditions(*tmp1);
210  (*tmp2) = convolution.doit((*tmp1), kernel);
211  swap = tmp1;
212  tmp1 = tmp2;
213  tmp2 = swap;
214  }
215 
216  if (verbose) std::cout << "Finished" << std::endl;
218  AimsData<T> ima3( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
219  conv2.convert( (*tmp1), ima3);
220 
221  return ima3;
222 }
223 
224 
225 template<typename T> inline
227 add_neumann_condition(const Point3df &p)
228 {
229  _has_neumann_condition = true;
230  _neumann_conditions.push_back(p);
231 }
232 
233 
234 
235 template<typename T>
237 update_neumann_conditions(AimsData<float> &ima)
238 {
239  int i, d;
240  Point3df dim(ima.dimX(), ima.dimY(), ima.dimZ());
241 
242  for (unsigned int i = 0; i < _neumann_conditions.size(); ++i)
243  {
244  //copy first encountred unmasked value among neighbours
245  const Point3df &p = _neumann_conditions[i];
246  for (i = 0; i < 3; ++i)
247  for (d = -1; d <= 1; d += 2)
248  {
249  Point3df p2(p);
250  p2[i] += d;
251  if (p2[i] < 0 or p2[i] >= dim[i]) continue;
252  if ((*(this->_mask))(p2[0], p2[1], p2[2]) != 0)
253  {
254  ima(p[0], p[1], p[2]) =ima(p2[0], p2[1], p2[2]);
255  i = d = 3;
256  break;
257  }
258  }
259  }
260 }
261 
262 
263 template<typename T> AimsData<T>
264 MaskedDiffusionSmoother<T, AimsData<short> >::
265 doSmoothing(const AimsData<T> & ima, int maxiter, bool verbose)
266 {
267  this->check(maxiter);
268 
269  AimsData<float> *tmp1, *tmp2, *swap;
270 
272  AimsData< float > imaF(ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT()),
273  ima2(ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT());
274  conv.convert(ima, imaF);
275  tmp1 = &imaF; tmp2 = &ima2;
276 
277  for (int i = 0; i < maxiter; i++)
278  {
279  if (verbose) std::cout << "+" << std::flush;
280  convolution((*tmp1), (*tmp2));
281  swap = tmp1;
282  tmp1 = tmp2;
283  tmp2 = swap;
284  }
285 
286  if (verbose) std::cout << "Finished" << std::endl;
288  AimsData<T> ima3( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
289  conv2.convert( (*tmp1), ima3);
290 
291  return ima3;
292 }
293 
294 template<typename T> void MaskedDiffusionSmoother<T, AimsData<short> >::
295 convolution(const AimsData<float> &ima1, AimsData<float> &ima2) const
296 {
297  int i, d, x, y, z, t, n;
298  int dx = ima1.dimX();
299  int dy = ima1.dimY();
300  int dz = ima1.dimZ();
301  int dt = ima1.dimT();
302  Point3df dim(ima1.dimX(), ima1.dimY(), ima1.dimZ());
303 
304  for (t = 0; t < dt; t++)
305  for (z = 0; z < dz; z++)
306  for (y = 0; y < dy; y++)
307  for (x = 0; x < dx; x++)
308  {
309  std::cout << "x,y,z = " << x <<", " << y << ", " << z << std::endl;
310  std::cout << "mask = " << ((*this->_mask)(x, y, z, t)) << std::endl;
311  if ((*this->_mask)(x, y, z, t) == this->_background or
312  (*this->_mask)(x, y, z, t) == this->_neumann_value)
313  {
314  ima2(x, y, z, t) = ima1(x, y, z, t);
315  continue;
316  }
317  n = 0;
318  T val = (T) 0.;
319 
320  for (i = 0; i < 3; ++i)
321  for (d = -1; d <= 1; d += 2)
322  {
323  Point3df p2(x, y, z);
324  p2[i] += d;
325  if (p2[i] < 0 or p2[i] >= dim[i]) continue;
326  if ((*this->_mask)(p2[0], p2[1], p2[2], t)
327  == this->_neumann_value) continue;
328  n++;
329  val += ima1(p2[0], p2[1], p2[2], t);
330  }
331  val -= n * ima1(x, y, z, t);
332  ima2(x, y, z, t) = ima1(x, y, z, t) + this->_dt * val;
333  }
334 }
335 //FIXME : manque une version en connexit√© 26
336 
337 }
338 
339 #endif
Make convolution only on a specified mask.
Definition: convol.h:98
int dimX() const
static AimsData< float > init_laplacian(void)
float sizeY() const
Heat diffusion with a volume of given datasource (Dirichlet conditions)
#define ForEach3d(thing, x, y, z)
int dimT() const
void setConstantSources(const AimsData< T > &, const T &background)
The template class to make convolutions.
Definition: convol.h:52
virtual void convert(const INP &in, OUTP &out) const
float sizeZ() const
virtual AimsData< T > doSmoothing(const AimsData< T > &ima, int maxiter, bool verbose=false)
float sizeX() const
int dimZ() const
void swap(scoped_ptr< T > &x, scoped_ptr< T > &y)
int dimY() const
AimsData< T > doit(AimsData< T > &, AimsData< T > &)
Definition: convol.h:120