aimsalgo 6.0.0
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
57namespace aims
58{
59
60template<class T> std::pair<carto::VolumeRef<T>, carto::VolumeRef<T> >
61CoupledDiffusion2DSmoother<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
carto::Converter::convert
virtual void convert(const INP &in, OUTP &out) const