A.I.M.S algorithms


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