A.I.M.S algorithms


peronaMalikSmoother_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 #ifndef AIMS_PRIMALSKETCH_PERONAMALIKSMOOTHER_D_H
36 #define AIMS_PRIMALSKETCH_PERONAMALIKSMOOTHER_D_H
37 
38 #include <cstdlib>
42 
43 namespace aims
44 {
45 
46  template<class T> AimsData<T> PeronaMalikSmoother<T>::doSmoothing(const AimsData<T> & ima, int maxiter, bool verbose)
47  {
48  if (maxiter >= 0)
49  {
50  AimsData<float> kernel(3, 3, 3, 1);
51  AimsData<float> *tmp1, *tmp2, *swap, grad;
52  int nbIter, i;
53  double nbFloat;
54  int x,y,z, count=0;
55  int sx=ima.dimX(), sy=ima.dimY(), sz=ima.dimZ();
56  std::map<float, int> valueSet;
57  std::map<float,int>::iterator itVal;
58  float nbThresh;
59  float gxm, gxp, gym, gyp, gzm, gzp, cxm, cxp, cym, cyp, czm, czp;
60 
61  GaussianGradient<float> gaussGrad(_sigma, _sigma, _sigma);
62  GaussianJacobian<float> gaussJac(_sigma, _sigma, _sigma);
64 
66  AimsData< float > imaF( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() ),
67  ima2( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
68  conv.convert( ima, imaF );
69 
70 
71  std::cout << "Evaluating gradient bound : " << std::endl;
72 
73  grad=gaussGrad.doit(imaF);
74  Writer<AimsData<float> > writerG( "gradient.ima" );
75  writerG.write( grad );
76 
77  // Evaluation of the gradient parameter from K
78 
79  std::cout << "Ordering values" << std::endl;
80 
81  nbThresh=sx*sz*sy*_K;
82 
83  for (z=0; z<sz; z++)
84  for (y=0; y<sy; y++)
85  for (x=0; x<sx; x++)
86  {
87  if (valueSet.find(grad(x,y,z)) != valueSet.end())
88  valueSet[grad(x,y,z)]++;
89  else
90  valueSet.insert(std::pair<float,int>(grad(x,y,z), 1));
91  }
92  itVal=valueSet.begin();
93 
94  std::cout << "Looking for bound with K= " << _K << " at value number " << nbThresh << std::endl;
95  for (; itVal!=valueSet.end(); ++itVal)
96  {
97  count+=(*itVal).second;
98  if ((count>=nbThresh) && ((count-1)<nbThresh))
99  _gradK=(*itVal).first;
100  }
101  std::cout << "Gradient bound has been set to " << _gradK << std::endl;
102 
103  tmp1=&imaF; tmp2=&ima2;
104  std::cout << "Starting " << maxiter << " iterations of diffusion process" << std::endl;
105 
106  for (i=0; i<maxiter; i++)
107  {
108  std::cout << "+" << std::flush;
109  grad=gaussGrad.doit((*tmp1));
110 
111  for (z=0; z<sz; z++)
112  for (y=0; y<sy; y++)
113  for (x=0; x<sx; x++)
114  {
115  if (x==0)
116  {
117  gxm=0; gxp=(*tmp1)(x+1,y,z) - (*tmp1)(x,y,z);
118  cxm=0; cxp=conductance( (grad(x+1,y,z)+grad(x,y,z))/2.0);
119  }
120  else if (x==sx-1)
121  {
122  gxm=(*tmp1)(x-1,y,z) - (*tmp1)(x,y,z); gxp=0;
123  cxm=conductance( (grad(x-1,y,z)+grad(x,y,z))/2.0); cxp=0;
124  }
125  else
126  {
127  gxm=(*tmp1)(x-1,y,z) - (*tmp1)(x,y,z); gxp=(*tmp1)(x+1,y,z) - (*tmp1)(x,y,z);
128  cxm=conductance( (grad(x-1,y,z)+grad(x,y,z))/2.0); cxp=conductance( (grad(x+1,y,z)+grad(x,y,z))/2.0);
129  }
130 
131  if (y==0)
132  {
133  gym=0; gyp=(*tmp1)(x,y+1,z) - (*tmp1)(x,y,z);
134  cym=0; cyp=conductance( (grad(x,y+1,z)+grad(x,y,z))/2.0);
135  }
136  else if (y==sy-1)
137  {
138  gym=(*tmp1)(x,y-1,z) - (*tmp1)(x,y,z); gyp=0;
139  cym=conductance( (grad(x,y-1,z)+grad(x,y,z))/2.0); cyp=0;
140  }
141  else
142  {
143  gym=(*tmp1)(x,y-1,z) - (*tmp1)(x,y,z); gyp=(*tmp1)(x,y+1,z) - (*tmp1)(x,y,z);
144  cym=conductance( (grad(x,y-1,z)+grad(x,y,z))/2.0); cyp=conductance( (grad(x,y+1,z)+grad(x,y,z))/2.0);
145  }
146 
147  if (z==0)
148  {
149  gzm=0; gzp=(*tmp1)(x,y,z+1) - (*tmp1)(x,y,z);
150  czm=0; czp=conductance( (grad(x,y,z+1)+grad(x,y,z))/2.0);
151  }
152  else if (z==sz-1)
153  {
154  gzm=(*tmp1)(x,y,z-1) - (*tmp1)(x,y,z); gzp=0;
155  czm=conductance( (grad(x,y,z-1)+grad(x,y,z))/2.0); czp=0;
156  }
157  else
158  {
159  gzm=(*tmp1)(x,y,z-1) - (*tmp1)(x,y,z); gzp=(*tmp1)(x,y,z+1) - (*tmp1)(x,y,z);
160  czm=conductance( (grad(x,y,z-1)+grad(x,y,z))/2.0); czp=conductance( (grad(x,y,z+1)+grad(x,y,z))/2.0);
161  }
162 
163  (*tmp2)(x,y,z)=(*tmp1)(x,y,z) + _dt *
164  ( cxm*gxm + cxp*gxp
165  +cym*gym + cyp*gyp
166  +czm*gzm + czp*gzp );
167  }
168  swap=tmp1;
169  tmp1=tmp2;
170  tmp2=swap;
171  }
172  std::cout << "Finished" << std::endl;
174  AimsData<T> ima3( ima.dimX(), ima.dimY(), ima.dimZ(), ima.dimT() );
175  conv2.convert( (*tmp1), ima3);
176 
177  return ima3;
178  }
179  else
180  {
181  std::cerr << "PeronaMalik smoother must have tIn < tOut" << std::endl;
182  exit(EXIT_FAILURE);
183  }
184  }
185 
186 }
187 
188 #endif
AimsData< T > doSmoothing(const AimsData< T > &ima, int maxiter, bool verbose=false)
int dimX() const
int dimT() const
virtual bool write(const T &obj, bool ascii=false, const std::string *format=0)
virtual void convert(const INP &in, OUTP &out) const
int dimZ() const
void swap(scoped_ptr< T > &x, scoped_ptr< T > &y)
int dimY() const
AimsData< float > doit(const AimsData< T > &)
Definition: ggradient.h:74