aimsalgo  5.0.5
Neuroimaging image processing
coupledDiffusion2DSmoother_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_COUPLEDDIFFUSION2DSMOOTHER_D_H
36 #define AIMS_PRIMALSKETCH_COUPLEDDIFFUSION2DSMOOTHER_D_H
37 
38 #include <cstdlib>
39 #include <aims/primalsketch/coupledDiffusion2DSmoother.h>
40 #include <aims/signalfilter/convol.h>
41 #include <cmath>
42 #include <cfloat>
43 
44 //
45 // WARNING : so far, this class implements a very particular type
46 // of coupled diffusion with soft constraints, i.e. :
47 // du/dt=Lapl(u) + div( (grad(u).grad(v)).grad(v) ) - (u-Cu)
48 // dv/dt=Lapl(v) + div( (grad(u).grad(v)).grad(u) ) - (v-Cv);
49 //
50 // which is a heat equation plus a term that keeps
51 // isocontours of both images as orthogonal as possible
52 //
53 // this is diffusion with constraints, that is with constant
54 // heat sources. So, not much use for anybody but me...
55 //
56 
57 namespace aims
58 {
59 
60 template<class T> std::pair<AimsData<T>, AimsData<T> >
61 CoupledDiffusion2DSmoother<T>::doSmoothing(const std::pair<AimsData<T>,
62  AimsData<T> > & ima,
63  const std::pair<AimsData<T>,
64  AimsData<T> > & constraints,
65  int maxiter,
66  bool /*verbose*/)
67 {
68  if ( maxiter >= 0)
69  {
70  AimsData<T> ima1, ima2;
71  ima1=ima.first;
72  ima2=ima.second;
73 
74  int PAS=20; //10
75  float epsilon=0.01;
76 
77  if ((ima.first.dimZ()>1) || (ima.second.dimZ()>1))
78  {
79  std::cerr << "coupledDiffusion2DSmoother: only for 2D images !!"
80  << std::endl;
81  exit(1);
82  }
83  if ((!hasSameDim(ima.first, ima.second))
84  || (!hasSameDim(ima.first, constraints.first))
85  || (!hasSameDim(ima.first, constraints.second)))
86  {
87  std::cerr << "coupledDiffusion2DSmoother: images do not all have the same size..." << std::endl;
88  exit(1);
89  }
90 
91  int sx=ima.first.dimX(), sy=ima.first.dimY(), x, y;
92  AimsData<float> *tmp1_1, *tmp1_2, *swap1;
93  AimsData<float> *tmp2_1, *tmp2_2, *swap2;
94  int i;
95  float lapl1, lapl2, div1, div2, lapx1, lapy1, lapx2, lapy2;
96 
97  std::vector<std::pair<int, int> > vcont1, vcont2;
98  std::vector<std::pair<int, int> >::iterator pt1, pt2;
99 
100  carto::Converter< AimsData<T>, AimsData<float> > conv;
101  AimsData< float > imaF1( sx, sy ),imaF2( sx, sy ),
102  imaB1( sx, sy ),imaB2( sx, sy ),
103  grad1_x( sx, sy ), grad1_y( sx, sy ),
104  grad2_x( sx, sy ), grad2_y( sx, sy );
105  conv.convert( ima.first, imaF1 );
106  conv.convert( ima.second, imaF2 );
107  float cdiff1, cdiff2, diff1, diff2, diffMax1, diffMax2;
108  tmp1_1=&imaF1; tmp1_2=&imaB1;
109  tmp2_1=&imaF2; tmp2_2=&imaB2;
110 
111  std::cout << "Initalizing images with constraints" << std::endl;
112 
113  for (y=0; y<sy; y++)
114  for (x=0; x<sx; x++)
115  {
116  if ((fabs(constraints.first(x,y)-1) < epsilon) || (fabs(constraints.first(x,y)-80)<epsilon))
117  {
118  imaF1(x,y)=(float) constraints.first(x,y);
119  vcont1.push_back(std::pair<int,int>(x,y));
120  }
121  if ((fabs(constraints.second(x,y)-1) < epsilon) || (fabs(constraints.second(x,y)-80)<epsilon))
122  {
123  imaF2(x,y)=(float) constraints.second(x,y);
124  vcont2.push_back(std::pair<int,int>(x,y));
125  }
126  }
127 
128 // cout << "Contraintes 1 :" << endl;
129 // for (pt1=vcont1.begin(); pt1!=vcont1.end(); ++pt1)
130 // cout << "(" << (*pt1).first << "," << (*pt1).second << ") "; fflush(stdout);
131 // cout << endl;
132 // cout << "Contraintes 2 :" << endl;
133 // for (pt2=vcont2.begin(); pt2!=vcont2.end(); ++pt2)
134 // cout << "(" << (*pt2).first << "," << (*pt2).second << ") "; fflush(stdout);
135 // cout << endl;
136  std::cout << "Starting " << maxiter
137  << " iterations of diffusion process" << std::endl;
138 
139  int sz=(maxiter/PAS) + 1;
140  AimsData<T> debug1(sx, sy, sz), debug2(sx, sy, sz);
141  AimsData<T> grad1x(sx, sy, sz), grad1y(sx, sy, sz), grad2x(sx, sy, sz), grad2y(sx, sy, sz);
142  AimsData<T> gugv(sx, sy), scal(sx, sy, sz), gu(sx, sy, sz), gv(sx, sy, sz);
143 
144  for (y=0; y<sy; y++)
145  for (x=0; x<sx; x++)
146  {
147  debug1(x, y, 0)=(*tmp1_1)(x, y);
148  debug2(x, y, 0)=(*tmp2_1)(x, y);
149  }
150  int z=1;
151  diffMax1=diffMax2=0.0;
152 
153  for (i=0; i<maxiter; i++)
154  {
155  if ((i%1)==0) //100
156  {
157  std::cout << "(t=" << i*_dt << ") -> diff=(" << diffMax1 << ","
158  << diffMax2 << ") - "<< std::endl;;
159  }
160 // std::cout << "(G"; fflush(stdout);
161  diffMax1=diffMax2=0.0;
162  // Calcul des gradients de chaque image, et du produit scalaire
163  for (y=0; y<sy; y++)
164  for (x=0; x<sx; x++)
165  {
166  if (x==0)
167  {
168  grad1_x(x,y)=((*tmp1_1)(x+1,y) - (*tmp1_1)(x,y))/2.0;
169  grad2_x(x,y)=((*tmp2_1)(x+1,y) - (*tmp2_1)(x,y))/2.0;
170  }
171  else if (x==sx-1)
172  {
173  grad1_x(x,y)=((*tmp1_1)(x,y) - (*tmp1_1)(x-1,y))/2.0;
174  grad2_x(x,y)=((*tmp2_1)(x,y) - (*tmp2_1)(x-1,y))/2.0;
175  }
176  else
177  {
178  grad1_x(x,y)=((*tmp1_1)(x+1,y) - (*tmp1_1)(x-1,y))/2.0;
179  grad2_x(x,y)=((*tmp2_1)(x+1,y) - (*tmp2_1)(x-1,y))/2.0;
180  }
181  if (y==0)
182  {
183  grad1_y(x,y)=((*tmp1_1)(x,y+1) - (*tmp1_1)(x,y))/2.0;
184  grad2_y(x,y)=((*tmp2_1)(x,y+1) - (*tmp2_1)(x,y))/2.0;
185  }
186  else if (y==sy-1)
187  {
188  grad1_y(x,y)=((*tmp1_1)(x,y) - (*tmp1_1)(x,y-1))/2.0;
189  grad2_y(x,y)=((*tmp2_1)(x,y) - (*tmp2_1)(x,y-1))/2.0;
190  }
191  else
192  {
193  grad1_y(x,y)=((*tmp1_1)(x,y+1) - (*tmp1_1)(x,y-1))/2.0;
194  grad2_y(x,y)=((*tmp2_1)(x,y+1) - (*tmp2_1)(x,y-1))/2.0;
195  }
196 
197  gugv(x,y)=grad1_x(x,y)*grad2_x(x,y) + grad1_y(x,y)*grad2_y(x,y);
198 
199  }
200 
201  // lissage de gugv pour stabilit� du sch�ma num�rique...
202 
203  Gaussian2DSmoothing<float> g2d(2.0,2.0);
204  AimsData<float> g1x(sx,sy), g2x(sx,sy), g1y(sx,sy), g2y(sx,sy);
205  g1x=g2d.doit(grad1_x);
206  g2x=g2d.doit(grad2_x);
207  g1y=g2d.doit(grad1_y);
208  g2y=g2d.doit(grad2_y);
209 
210  for (y=0; y<sy; y++)
211  for (x=0; x<sx; x++)
212  {
213  double prod=g1x(x,y)*g2x(x,y) + g1y(x,y)*g2y(x,y);
214  gugv(x,y)=float(2*prod*exp(-prod*prod)); // Fonction de conductance !!!!!!!
215  }
216 
217  // passe de calcul du terme total
218 
219  float div1x, div1y, div2x, div2y;
220  for (y=0; y<(sy); y++)
221  for (x=0; x<(sx); x++)
222  {
223  if (x==0)
224  {
225  lapx1=(*tmp1_1)(x+1,y) + (*tmp1_1)(x,y) -2*(*tmp1_1)(x,y);
226  lapx2=(*tmp2_1)(x+1,y) + (*tmp2_1)(x,y) -2*(*tmp2_1)(x,y);
227  div1x=( (g1x(x+1,y)*g2x(x+1,y) + g1y(x+1,y)*g2x(x+1,y) )*g2x(x+1,y)
228  - (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g2x(x,y) )/2.0;
229  div2x=( (g1x(x+1,y)*g2x(x+1,y) + g1y(x+1,y)*g2x(x+1,y) )*g1x(x+1,y)
230  - (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g1x(x,y) )/2.0;
231  }
232  else if (x==(sx-1))
233  {
234  lapx1=(*tmp1_1)(x,y) + (*tmp1_1)(x-1,y) -2*(*tmp1_1)(x,y);
235  lapx2=(*tmp2_1)(x,y) + (*tmp2_1)(x-1,y) -2*(*tmp2_1)(x,y);
236  div1x=( (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g2x(x,y)
237  - (g1x(x-1,y)*g2x(x-1,y) + g1y(x-1,y)*g2x(x-1,y) )*g2x(x-1,y) )/2.0;
238  div2x=( (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g1x(x,y)
239  - (g1x(x-1,y)*g2x(x-1,y) + g1y(x-1,y)*g2x(x-1,y) )*g1x(x-1,y) )/2.0;
240  }
241  else
242  {
243  lapx1=(*tmp1_1)(x+1,y) + (*tmp1_1)(x-1,y) -2*(*tmp1_1)(x,y);
244  lapx2=(*tmp2_1)(x+1,y) + (*tmp2_1)(x-1,y) -2*(*tmp2_1)(x,y);
245  div1x=( (g1x(x+1,y)*g2x(x+1,y) + g1y(x+1,y)*g2x(x+1,y) )*g2x(x+1,y)
246  - (g1x(x-1,y)*g2x(x-1,y) + g1y(x-1,y)*g2x(x-1,y) )*g2x(x-1,y) )/2.0;
247  div2x=( (g1x(x+1,y)*g2x(x+1,y) + g1y(x+1,y)*g2x(x+1,y) )*g1x(x+1,y)
248  - (g1x(x-1,y)*g2x(x-1,y) + g1y(x-1,y)*g2x(x-1,y) )*g1x(x-1,y) )/2.0;
249  }
250  if (y==0)
251  {
252  lapy1=(*tmp1_1)(x,y+1) + (*tmp1_1)(x,y) -2*(*tmp1_1)(x,y);
253  lapy2=(*tmp2_1)(x,y+1) + (*tmp2_1)(x,y) -2*(*tmp2_1)(x,y);
254  div1y=( (g1x(x,y+1)*g2x(x,y+1) + g1y(x,y+1)*g2x(x,y+1) )*g2y(x,y+1)
255  - (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g2y(x,y) )/2.0;
256  div2y=( (g1x(x,y+1)*g2x(x,y+1) + g1y(x,y+1)*g2x(x,y+1) )*g1y(x,y+1)
257  - (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g1y(x,y) )/2.0;
258  }
259  else if (y==(sy-1))
260  {
261  lapy1=(*tmp1_1)(x,y) + (*tmp1_1)(x,y-1) -2*(*tmp1_1)(x,y);
262  lapy2=(*tmp2_1)(x,y) + (*tmp2_1)(x,y-1) -2*(*tmp2_1)(x,y);
263  div1y=( (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g2y(x,y)
264  - (g1x(x,y-1)*g2x(x,y-1) + g1y(x,y-1)*g2x(x,y-1) )*g2y(x,y-1) )/2.0;
265  div2y=( (g1x(x,y)*g2x(x,y) + g1y(x,y)*g2x(x,y) )*g1y(x,y)
266  - (g1x(x,y-1)*g2x(x,y-1) + g1y(x,y-1)*g2x(x,y-1) )*g1y(x,y-1) )/2.0;
267  }
268  else
269  {
270  lapy1=(*tmp1_1)(x,y+1) + (*tmp1_1)(x,y-1) -2*(*tmp1_1)(x,y);
271  lapy2=(*tmp2_1)(x,y+1) + (*tmp2_1)(x,y-1) -2*(*tmp2_1)(x,y);
272  div1y=( (g1x(x,y+1)*g2x(x,y+1) + g1y(x,y+1)*g2x(x,y+1) )*g2y(x,y+1)
273  - (g1x(x,y-1)*g2x(x,y-1) + g1y(x,y-1)*g2x(x,y-1) )*g2y(x,y-1) )/2.0;
274  div2y=( (g1x(x,y+1)*g2x(x,y+1) + g1y(x,y+1)*g2x(x,y+1) )*g1y(x,y+1)
275  - (g1x(x,y-1)*g2x(x,y-1) + g1y(x,y-1)*g2x(x,y-1) )*g1y(x,y-1) )/2.0;
276  }
277 
278  lapl1=lapx1+lapy1;
279  lapl2=lapx2+lapy2;
280  div1 = div1x+div1y;
281  div2 = div2x+div2y;
282 
283  if ((fabs(constraints.first(x,y))>epsilon) && (fabs(constraints.first(x,y)-0) > epsilon) && (fabs(constraints.first(x,y)-80) > epsilon))
284  cdiff1=(*tmp1_1)(x,y) - constraints.first(x,y);
285  else
286  cdiff1=0.0;
287  if ((fabs(constraints.second(x,y)) > epsilon) && (fabs(constraints.second(x,y)-0) > epsilon) && (fabs(constraints.second(x,y)-80)> epsilon))
288  cdiff2=(*tmp2_1)(x,y) - constraints.second(x,y);
289  else
290  cdiff2=0.0;
291  (*tmp1_2)(x,y) = (*tmp1_1)(x,y) + _dt*(_alpha*lapl1 + _beta*div1 - _gamma*cdiff1);
292  (*tmp2_2)(x,y) = (*tmp2_1)(x,y) + _dt*(_alpha*lapl2 + _beta*div2 - _gamma*cdiff2);
293  diff1 = fabs((*tmp1_2)(x,y) - (*tmp1_1)(x,y));
294  diff2 = fabs((*tmp2_2)(x,y) - (*tmp2_1)(x,y));
295 
296  diffMax1 += diff1; diffMax2 += diff2;
297  }
298 // std::cout << " OK) "; fflush(stdout);
299  diffMax1=diffMax1/float(sx*sy);
300  diffMax2=diffMax2/float(sx*sy);
301  // on fixe les contraintes � leur valeur initiale
302  for (pt1=vcont1.begin(); pt1!=vcont1.end(); ++pt1)
303  (*tmp1_2)((*pt1).first, (*pt1).second)= (float) constraints.first((*pt1).first, (*pt1).second);
304  for (pt2=vcont2.begin(); pt2!=vcont2.end(); ++pt2)
305  (*tmp2_2)((*pt2).first, (*pt2).second)= (float) constraints.second((*pt2).first, (*pt2).second);
306  swap1=tmp1_1;
307  tmp1_1=tmp1_2;
308  tmp1_2=swap1;
309  swap2=tmp2_1;
310  tmp2_1=tmp2_2;
311  tmp2_2=swap2;
312 
313  if ((i%PAS)==0)
314  {
315  for (y=0; y<sy; y++)
316  for (x=0; x<sx; x++)
317  {
318  debug1(x, y, z)=(*tmp1_1)(x, y);
319  debug2(x, y, z)=(*tmp2_1)(x, y);
320  grad1x(x, y, z)=grad1_x(x,y);
321  grad1y(x, y, z)=grad1_y(x,y);
322  grad2x(x, y, z)=grad2_x(x,y);
323  grad2y(x, y, z)=grad2_y(x,y);
324  scal(x, y, z)=gugv(x,y);
325  }
326  z++;
327  }
328  }
329 
330  std::cout << "Finished" << std::endl;
331  carto::Converter< AimsData<float>, AimsData<T> > conv2;
332  AimsData<T> out1( sx, sy), out2(sx,sy);
333  conv2.convert( (*tmp1_1),out1);
334  conv2.convert( (*tmp2_1),out2);
335 
336  //Pour debug : evolution des iso-contours
337 
338  std::cout << "Computing and writing iso-contours" << std::endl;
339  AimsData<uint8_t> iso( sx, sy, sz);
340  for (z=0; z<sz; z++)
341  for (y=0; y<sy-1; y++)
342  for (x=0; x<sx-1; x++)
343  {
344  iso(x,y,z)=0;
345  }
346  for (z=0; z<sz; z++)
347  for (y=0; y<sy-1; y++)
348  for (x=0; x<sx-1; x++)
349  {
350  int value;
351  for (value=0; value < 80; value=value+5)
352  {
353 // if ((constraints.first(x,y)>0) || (constraints.second(x,y)>0))
354 // iso(x,y,z)=2;
355  if ((((debug1(x,y,z)-float(value))*(debug1(x,y+1,z)-float(value)))<=0)
356  || (((debug1(x,y,z)-float(value))*(debug1(x+1,y,z)-float(value)))<=0))
357  iso(x,y,z)=1;
358  if ((((debug2(x,y,z)-float(value))*(debug2(x,y+1,z)-float(value)))<=0)
359  || (((debug2(x,y,z)-float(value))*(debug2(x+1,y,z)-float(value)))<=0))
360  iso(x,y,z)=1;
361  }
362  }
363 
364  Writer<AimsData<uint8_t> > writerIso( "grille.ima" );
365  writerIso.write(iso);
366 
367  Writer<AimsData<float> > writerD1( "evolution1.ima" );
368  Writer<AimsData<float> > writerD2( "evolution2.ima" );
369  writerD1.write(debug1);
370  writerD2.write(debug2);
371  Writer<AimsData<float> > writerG1x( "grad1x.ima" );
372  Writer<AimsData<float> > writerG2x( "grad2x.ima" );
373  writerG1x.write(grad1x);
374  writerG2x.write(grad2x);
375  Writer<AimsData<float> > writerG1y( "grad1y.ima" );
376  Writer<AimsData<float> > writerG2y( "grad2y.ima" );
377  writerG1y.write(grad1y);
378  writerG2y.write(grad2y);
379  Writer<AimsData<float> > writerGuGv( "gugv.ima" );
380  writerGuGv.write(scal);
381 
382  return (std::pair<AimsData<T>, AimsData<T> >(out1, out2) );
383  }
384  else
385  {
386  std::cerr << "coupledDiffusion2DSmoother: must have tIn < tOut"
387  << std::endl;
388  exit(1);
389  }
390 
391 }
392 
393 }
394 
395 #endif
aims::CoupledDiffusion2DSmoother::doSmoothing
std::pair< AimsData< T >, AimsData< T > > doSmoothing(const std::pair< AimsData< T >, AimsData< T > > &ima, const std::pair< AimsData< T >, AimsData< T > > &constraint, int maxiter, bool verbose=false)
Definition: coupledDiffusion2DSmoother_d.h:61
carto::Converter::convert
virtual void convert(const INP &in, OUTP &out) const
aims::Writer::write
virtual bool write(const T &obj, bool ascii=false, const std::string *format=0)
aims::hasSameDim
bool hasSameDim(const AimsData< T > &v1, const AimsData< T > &v2)
Gaussian2DSmoothing::doit
AimsData< T > doit(const AimsData< T > &)
Definition: g2dsmooth.h:70
aims::util::epsilon
T epsilon()
Definition: dynamic_d.h:50
iterator
T * iterator