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