1 #include "mrilib.h"
2
3 #undef GOOD
4 #define GOOD(i) (mask==NULL || mask[i])
5
6 static int dontcheckplus = 0 ;
FHWM_1dif_dontcheckplus(int i)7 void FHWM_1dif_dontcheckplus( int i ){ dontcheckplus = i; }
8
9 /*---------------------------------------------------------------------------*/
10
mriarr_estimate_FWHM_1dif(MRI_IMARR * imar,byte * mask,int unif)11 THD_fvec3 mriarr_estimate_FWHM_1dif( MRI_IMARR *imar, byte *mask, int unif )
12 {
13 int nar=IMARR_COUNT(imar) , ii ;
14 THD_fvec3 sv ;
15 float cx,cy,cz , fx,fy,fz ;
16 int nx,ny,nz , nvox , kk ;
17 MRI_IMAGE *medim=NULL , *madim=NULL ;
18 float *medar=NULL , *madar=NULL , *sar=NULL ;
19
20 unif = unif && (nar > 2) ;
21 nvox = IMARR_SUBIM(imar,0)->nvox ;
22 if( unif ){
23 MRI_IMARR *qar = IMARR_medmad_bricks( imar ) ;
24 medim = IMARR_SUBIM(qar,0); medar = MRI_FLOAT_PTR(medim);
25 madim = IMARR_SUBIM(qar,1); madar = MRI_FLOAT_PTR(madim); FREE_IMARR(qar);
26 for( kk=0 ; kk < nvox ; kk++ )
27 if( madar[kk] != 0.0f ) madar[kk] = 1.0f / madar[kk] ;
28 }
29 cx = cy = cz = 0.0f ; nx = ny = nz = 0 ;
30 for( ii=0 ; ii < nar ; ii++ ){
31 if( unif ){
32 sar = MRI_FLOAT_PTR( IMARR_SUBIM(imar,ii) ) ;
33 for( kk=0 ; kk < nvox ; kk++ ) sar[kk] = (sar[kk]-medar[kk])*madar[kk] ;
34 }
35 sv = mri_estimate_FWHM_1dif( IMARR_SUBIM(imar,ii) , mask ) ;
36 UNLOAD_FVEC3(sv,fx,fy,fz) ;
37 /*** INFO_message(" sub-brick[%d]: fx=%g fy=%g fz=%g",ii,fx,fy,fz) ; ***/
38 if( fx > 0.0f ){ cx += fx ; nx++ ; }
39 if( fy > 0.0f ){ cy += fy ; ny++ ; }
40 if( fz > 0.0f ){ cz += fz ; nz++ ; }
41 }
42 cx = (nx==0) ? -1.0f : cx / nx ;
43 cy = (ny==0) ? -1.0f : cy / ny ;
44 cz = (nz==0) ? -1.0f : cz / nz ;
45 LOAD_FVEC3(sv,cx,cy,cz) ;
46 if( unif ){ mri_free(medim); mri_free(madim); }
47 return sv ;
48 }
49
50 /*---------------------------------------------------------------------------*/
51 /*! Routine to estimate Gaussian FWHM of data brick, using first differences.
52 - A negative return value indicates an error condition in that direction
53 (e.g., FWHM(z) == -1.0 when nz == 1).
54 - Adapted from original 3dFWHM.c by BD Ward.
55 - 31 Oct 2006 -- BOO!
56 - A sister function is in SUMA_GeomComp.c if changes are made to this
57 one, please notify the half-gloved man.
58 -----------------------------------------------------------------------------*/
59
mri_estimate_FWHM_1dif(MRI_IMAGE * im,byte * mask)60 THD_fvec3 mri_estimate_FWHM_1dif( MRI_IMAGE *im , byte *mask )
61 {
62 int nx; /* number of voxels along x-axis */
63 int ny; /* number of voxels along y-axis */
64 int nz; /* number of voxels along z-axis */
65 int nxy, nxyz; /* total number of voxels */
66 int ixyz; /* voxel index */
67 double dx; /* voxel size along x-axis */
68 double dy; /* voxel size along y-axis */
69 double dz; /* voxel size along z-axis */
70 int ix, jy, kz, ixyz2;
71 double fsum, fsq, var , arg ;
72 double dfdx, dfdxsum, dfdxsq, varxx;
73 double dfdy, dfdysum, dfdysq, varyy;
74 double dfdz, dfdzsum, dfdzsq, varzz;
75 int count, countx, county, countz;
76
77 float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
78 THD_fvec3 fw_xyz ;
79 MRI_IMAGE *lim ; float *fim ;
80
81 /*----- initialize local variables -----*/
82
83 LOAD_FVEC3(fw_xyz,sx,sy,sz) ; /* load with error flags */
84
85 if( im == NULL || mri_allzero(im) ) return fw_xyz ;
86 lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
87 fim = MRI_FLOAT_PTR(lim) ;
88 nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;
89
90 /*----- estimate the variance of the data -----*/
91
92 fsum = 0.0; fsq = 0.0; count = 0;
93 for (ixyz = 0; ixyz < nxyz; ixyz++){
94 if( GOOD(ixyz) ){ count++; arg = fim[ixyz]; fsum += arg; fsq += arg*arg; }
95 }
96 if( count < 9 || fsq <= 0.0 ){ /* no data? */
97 if( lim != im ) mri_free(lim) ;
98 return fw_xyz ;
99 }
100 var = (fsq - (fsum * fsum)/count) / (count-1.0);
101 if( var <= 0.0 ){ /* crappy data? */
102 if( lim != im ) mri_free(lim) ;
103 return fw_xyz ;
104 }
105
106 /*----- estimate the partial derivatives -----*/
107
108 dfdxsum = 0.0; dfdysum = 0.0; dfdzsum = 0.0;
109 dfdxsq = 0.0; dfdysq = 0.0; dfdzsq = 0.0;
110 countx = 0; county = 0; countz = 0;
111 for (ixyz = 0; ixyz < nxyz; ixyz++){
112 if( GOOD(ixyz) ){
113
114 arg = fim[ixyz] ;
115 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
116
117 if( ix+1 < nx ){
118 ixyz2 = ixyz+1 ;
119 if( dontcheckplus || GOOD(ixyz2) ){
120 dfdx = (fim[ixyz2] - arg) ;
121 dfdxsum += dfdx; dfdxsq += dfdx * dfdx; countx++;
122 }
123 }
124
125 if( jy+1 < ny ){
126 ixyz2 = ixyz+nx ;
127 if( dontcheckplus || GOOD(ixyz2) ){
128 dfdy = (fim[ixyz2] - arg) ;
129 dfdysum += dfdy; dfdysq += dfdy * dfdy; county++;
130 }
131 }
132
133 if( kz+1 < nz ){
134 ixyz2 = ixyz+nxy ;
135 if( dontcheckplus || GOOD(ixyz2) ){
136 dfdz = (fim[ixyz2] - arg) ;
137 dfdzsum += dfdz; dfdzsq += dfdz * dfdz; countz++;
138 }
139 }
140 }
141 }
142
143 /*----- estimate the variance of the partial derivatives -----*/
144
145 varxx = (countx < 6) ? 0.0
146 : (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1.0);
147
148 varyy = (county < 6) ? 0.0
149 : (dfdysq - (dfdysum * dfdysum)/county) / (county-1.0);
150
151 varzz = (countz < 6) ? 0.0
152 : (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1.0);
153
154 /*----- now estimate the FWHMs -----*/
155 /*---- the model is the correlation of the neighbors in x is
156 exp[ -dx^2 / (4 * sx^2) ]
157 which is the result of passing a Gaussian exp[-x^2/(2*sx^2)]
158 convolution over a white noise field. Then the covariance of
159 the neigbhors is
160 2*V * { 1 - exp[ -dx^2 / (4*sx^2) ] }
161 where V = variance of the noise (var). -----*/
162 /*---- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor -----*/
163
164 dx = lim->dx; dy = lim->dy; dz = lim->dz;
165
166 arg = 1.0 - 0.5*(varxx/var);
167 /* ININFO_message("dx: varxx/var=%g arg=%g",varxx/var,arg) ; */
168 if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) ) * dx;
169
170 arg = 1.0 - 0.5*(varyy/var);
171 if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) ) * dy;
172 /* ININFO_message("dy: varyy/var=%g arg=%g",varyy/var,arg) ; */
173
174 arg = 1.0 - 0.5*(varzz/var);
175 if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) ) * dz;
176 /* ININFO_message("dz: varzz/var=%g arg=%g",varzz/var,arg) ; */
177
178 LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
179 if( lim != im ) mri_free(lim) ;
180 return fw_xyz ;
181 }
182
183 /*---------------------------------------------------------------------------*/
184 /*! Routine to estimate Gaussian FWHM of data brick, using differences
185 between 1st and 2cd nearest neighbors.
186 - A negative return value indicates an error condition in that direction
187 (e.g., FWHM(z) == -1.0 when nz == 1).
188 -----------------------------------------------------------------------------*/
189
mri_estimate_FWHM_12dif(MRI_IMAGE * im,byte * mask)190 THD_fvec3 mri_estimate_FWHM_12dif( MRI_IMAGE *im , byte *mask )
191 {
192 int nx; /* number of voxels along x-axis */
193 int ny; /* number of voxels along y-axis */
194 int nz; /* number of voxels along z-axis */
195 int nxy, nxyz; /* total number of voxels */
196 int ixyz; /* voxel index */
197 double dx; /* voxel size along x-axis */
198 double dy; /* voxel size along y-axis */
199 double dz; /* voxel size along z-axis */
200 int ix, jy, kz, qm,qp ;
201 double dx1 , dx1sum , dx1sqq , vx1,vy1,vz1 , arg ;
202 double dx2 , dx2sum , dx2sqq , vx2,vy2,vz2 ;
203 double dy1 , dy1sum , dy1sqq ;
204 double dy2 , dy2sum , dy2sqq ;
205 double dz1 , dz1sum , dz1sqq ;
206 double dz2 , dz2sum , dz2sqq ;
207 int countx, county, countz;
208
209 float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
210 THD_fvec3 fw_xyz ;
211 MRI_IMAGE *lim ; float *fim ;
212
213 /*----- initialize local variables -----*/
214
215 LOAD_FVEC3(fw_xyz,sx,sy,sz) ; /* load with error flags */
216
217 if( im == NULL ) return fw_xyz ;
218 lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
219 fim = MRI_FLOAT_PTR(lim) ;
220 nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;
221
222 /*----- loop over voxels, compute differences, sum and sum squares -----*/
223
224 countx = county = countz = 0 ;
225 dx1sum = dx2sum = dy1sum = dy2sum = dz1sum = dz2sum = 0.0 ;
226 dx1sqq = dx2sqq = dy1sqq = dy2sqq = dz1sqq = dz2sqq = 0.0 ;
227 for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
228 if( GOOD(ixyz) ){
229 arg = fim[ixyz] ; IJK_TO_THREE (ixyz, ix,jy,kz, nx,nxy);
230
231 if( ix-1 >= 0 && ix+1 < nx ){
232 qp = ixyz+1 ; qm = ixyz-1 ;
233 if( GOOD(qp) && GOOD(qm) ){
234 dx1 = fim[qp]-arg ; dx2 = fim[qp]-fim[qm] ;
235 dx1sum += dx1 ; dx2sum += dx2 ;
236 dx1sqq += dx1*dx1 ; dx2sqq += dx2*dx2 ;
237 countx++ ;
238 }
239 }
240
241 if( jy-1 >= 0 && jy+1 < ny ){
242 qp = ixyz+nx ; qm = ixyz-nx ;
243 if( GOOD(qp) && GOOD(qm) ){
244 dy1 = fim[qp]-arg ; dy2 = fim[qp]-fim[qm] ;
245 dy1sum += dy1 ; dy2sum += dy2 ;
246 dy1sqq += dy1*dy1 ; dy2sqq += dy2*dy2 ;
247 county++ ;
248 }
249 }
250
251 if( kz-1 >= 0 && kz+1 < nz ){
252 qp = ixyz+nxy ; qm = ixyz-nxy ;
253 if( GOOD(qp) && GOOD(qm) ){
254 dz1 = fim[qp]-arg ; dz2 = fim[qp]-fim[qm] ;
255 dz1sum += dz1 ; dz2sum += dz2 ;
256 dz1sqq += dz1*dz1 ; dz2sqq += dz2*dz2 ;
257 countz++ ;
258 }
259 }
260 }
261 }
262
263 /*----- estimate variances of differences -----*/
264
265 vx1 = (countx < 6) ? 0.0 : (dx1sqq - (dx1sum*dx1sum)/countx) / (countx-1.0) ;
266 vy1 = (county < 6) ? 0.0 : (dy1sqq - (dy1sum*dy1sum)/county) / (county-1.0) ;
267 vz1 = (countz < 6) ? 0.0 : (dz1sqq - (dz1sum*dz1sum)/countz) / (countz-1.0) ;
268
269 vx2 = (countx < 6) ? 0.0 : (dx2sqq - (dx2sum*dx2sum)/countx) / (countx-1.0) ;
270 vy2 = (county < 6) ? 0.0 : (dy2sqq - (dy2sum*dy2sum)/county) / (county-1.0) ;
271 vz2 = (countz < 6) ? 0.0 : (dz2sqq - (dz2sum*dz2sum)/countz) / (countz-1.0) ;
272
273 #if 0
274 fprintf(stderr,"countx=%d dx1sum=%g dx1sqq=%g vx1=%g\n",countx,dx1sum,dx1sqq,vx1) ;
275 fprintf(stderr,"countx=%d dx2sum=%g dx2sqq=%g vx2=%g\n",countx,dx2sum,dx2sqq,vx2) ;
276 #endif
277
278 /*----- now estimate the FWHMs -----*/
279
280 dx = lim->dx; dy = lim->dy; dz = lim->dz;
281 if( lim != im ) mri_free(lim) ;
282
283 /*--- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor ---*/
284 /*--- y = cbrt(12*sqrt(48-120*r+81*r*r)+108*r-80), and then
285 x = y/6 - 4/(3*y) - 1/3
286 is the real solution to the equation (1-x^4)/(1-x) = r > 1 ;
287 here, r = vx2/vx1 = ratio of variances at Delta=2*dx vs. Delta=1*dx;
288 x = exp[-dx**2/(4*sigma**2)] = correlation coefficient of neighbors;
289 we solve for x, then use that to solve for sigma, and scale to FWHM --*/
290
291 if( vx1 > 0.0 && vx2 > vx1 ){
292 arg = vx2 / vx1 ;
293 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
294 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
295 if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dx ;
296 }
297
298 if( vy1 > 0.0 && vy2 > vy1 ){
299 arg = vy2 / vy1 ;
300 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
301 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
302 if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dy ;
303 }
304
305 if( vz1 > 0.0 && vz2 > vz1 ){
306 arg = vz2 / vz1 ;
307 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
308 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
309 if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dz ;
310 }
311
312 LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
313 return fw_xyz ;
314 }
315
316 /*---------------------------------------------------------------------------*/
317 /*! Routine to estimate Gaussian FWHM of data brick, using differences
318 between 1st and 2cd nearest neighbors.
319 - A negative return value indicates an error condition in that direction
320 (e.g., FWHM(z) == -1.0 when nz == 1).
321 -----------------------------------------------------------------------------*/
322
mri_estimate_FWHM_12dif_MAD(MRI_IMAGE * im,byte * mask)323 THD_fvec3 mri_estimate_FWHM_12dif_MAD( MRI_IMAGE *im , byte *mask )
324 {
325 int nx; /* number of voxels along x-axis */
326 int ny; /* number of voxels along y-axis */
327 int nz; /* number of voxels along z-axis */
328 int nxy, nxyz; /* total number of voxels */
329 int ixyz; /* voxel index */
330 float dx; /* voxel size along x-axis */
331 float dy; /* voxel size along y-axis */
332 float dz; /* voxel size along z-axis */
333 int ix, jy, kz, qm,qp ;
334 float vx1,vy1,vz1 , arg ;
335 float vx2,vy2,vz2 ;
336 float wx1,wy1,wz1 , brg ;
337 float wx2,wy2,wz2 ;
338 int countx, county, countz;
339
340 int ngood ; float *dx1ar,*dy1ar,*dz1ar , *dx2ar,*dy2ar,*dz2ar ;
341
342 float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
343 THD_fvec3 fw_xyz ;
344 MRI_IMAGE *lim ; float *fim ;
345
346 /*----- initialize local variables -----*/
347
348 LOAD_FVEC3(fw_xyz,sx,sy,sz) ; /* load with error flags */
349
350 if( im == NULL ) return fw_xyz ;
351 lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
352 fim = MRI_FLOAT_PTR(lim) ;
353 nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;
354
355 for( ngood=ixyz=0 ; ixyz < nxyz ; ixyz++ ) if( GOOD(ixyz) ) ngood++ ;
356 if( ngood < 9 ) return fw_xyz ;
357
358 dx1ar = (float *)malloc(sizeof(float)*ngood) ;
359 dy1ar = (float *)malloc(sizeof(float)*ngood) ;
360 dz1ar = (float *)malloc(sizeof(float)*ngood) ;
361 dx2ar = (float *)malloc(sizeof(float)*ngood) ;
362 dy2ar = (float *)malloc(sizeof(float)*ngood) ;
363 dz2ar = (float *)malloc(sizeof(float)*ngood) ;
364
365 /*----- loop over voxels, compute differences, sum and sum squares -----*/
366
367 countx = county = countz = 0 ;
368 for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
369 if( GOOD(ixyz) ){
370 arg = fim[ixyz] ; IJK_TO_THREE (ixyz, ix,jy,kz, nx,nxy);
371
372 if( ix-1 >= 0 && ix+1 < nx ){
373 qp = ixyz+1 ; qm = ixyz-1 ;
374 if( GOOD(qp) && GOOD(qm) ){
375 dx1ar[countx] = fim[qp]-arg ; dx2ar[countx] = fim[qp]-fim[qm] ; countx++ ;
376 }
377 }
378
379 if( jy-1 >= 0 && jy+1 < ny ){
380 qp = ixyz+nx ; qm = ixyz-nx ;
381 if( GOOD(qp) && GOOD(qm) ){
382 dy1ar[county] = fim[qp]-arg ; dy2ar[county] = fim[qp]-fim[qm] ; county++ ;
383 }
384 }
385
386 if( kz-1 >= 0 && kz+1 < nz ){
387 qp = ixyz+nxy ; qm = ixyz-nxy ;
388 if( GOOD(qp) && GOOD(qm) ){
389 dz1ar[countz] = fim[qp]-arg ; dz2ar[countz] = fim[qp]-fim[qm] ; countz++ ;
390 }
391 }
392 }
393 }
394
395 /*----- estimate variances of differences -----*/
396
397 qmedmadmeanad_float( countx, dx1ar, NULL, &vx1, &wx1 ); vx1 = vx1*vx1; wx1 = wx1*wx1 ;
398 qmedmadmeanad_float( county, dy1ar, NULL, &vy1, &wy1 ); vy1 = vy1*vy1; wy1 = wy1*wy1 ;
399 qmedmadmeanad_float( countz, dz1ar, NULL, &vz1, &wz1 ); vz1 = vz1*vz1; wz1 = wz1*wz1 ;
400 qmedmadmeanad_float( countx, dx2ar, NULL, &vx2, &wx2 ); vx2 = vx2*vx2; wx2 = wx2*wx2 ;
401 qmedmadmeanad_float( county, dy2ar, NULL, &vy2, &wy2 ); vy2 = vy2*vy2; wy2 = wy2*wy2 ;
402 qmedmadmeanad_float( countz, dz2ar, NULL, &vz2, &wz2 ); vz2 = vz2*vz2; wz2 = wz2*wz2 ;
403
404 free(dx1ar) ; free(dy1ar) ; free(dz1ar) ;
405 free(dx2ar) ; free(dy2ar) ; free(dz2ar) ;
406
407 /*----- now estimate the FWHMs -----*/
408
409 dx = lim->dx; dy = lim->dy; dz = lim->dz;
410 if( lim != im ) mri_free(lim) ;
411
412 /*--- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor ---*/
413 /*--- y = cbrt(12*sqrt(48-120*r+81*r*r)+108*r-80), and then
414 x = y/6 - 4/(3*y) - 1/3
415 is the real solution to the equation (1-x^4)/(1-x) = r > 1 ;
416 here, r = vx2/vx1 = ratio of variances at Delta=2*dx vs. Delta=1*dx;
417 x = exp[-dx**2/(4*sigma**2)] = correlation coefficient of neighbors;
418 we solve for x, then use that to solve for sigma, and scale to FWHM --*/
419
420 /*--- Modified 11 Aug 2015, to use the mean absolute deviations (wx1 & wx2)
421 instead of the median absolute deviations (vx1 & vx2) if the
422 computation fails for the median absolute deviations (the backup) ---*/
423
424 if( vx1 > 0.0 && vx2 > vx1 ){
425 arg = vx2 / vx1 ; brg = wx2 / wx1 ;
426 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
427 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
428 if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dx ;
429 else if( brg > 1.0f ){
430 brg = cbrt(12.0*sqrt(48.0 - 120.0*brg + 81.0*brg*brg) + 108.0*brg - 80.0) ;
431 brg = brg/6.0 - 4.0/(3.0*brg) - 1.0/3.0 ;
432 if( brg > 0.0 && brg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(brg)) )*dx ;
433 }
434 /* INFO_message("FWHM x: vx1=%f vx2=%f arg=%f sx=%f countx=%d",vx1,vx2,arg,sx,countx) ; */
435 } else {
436 /* INFO_message("FWHM x: vx1=%f vx2=%f NOT GOOD",vx1,vx2) ; */
437 }
438
439 if( vy1 > 0.0 && vy2 > vy1 ){
440 arg = vy2 / vy1 ; brg = wy2 / wy1 ;
441 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
442 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
443 if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dy ;
444 else if( brg > 1.0f ){
445 brg = cbrt(12.0*sqrt(48.0 - 120.0*brg + 81.0*brg*brg) + 108.0*brg - 80.0) ;
446 brg = brg/6.0 - 4.0/(3.0*brg) - 1.0/3.0 ;
447 if( brg > 0.0 && brg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(brg)) )*dy ;
448 }
449 /* INFO_message("FWHM y: vy1=%f vy2=%f arg=%f sy=%f county=%d",vy1,vy2,arg,sy,county) ; */
450 } else {
451 /* INFO_message("FWHM y: vy1=%f vy2=%f NOT GOOD",vy1,vy2) ; */
452 }
453
454 if( vz1 > 0.0 && vz2 > vz1 ){
455 arg = vz2 / vz1 ; brg = wz2 / wz1 ;
456 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
457 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
458 if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dz ;
459 else if( brg > 1.0f ){
460 brg = cbrt(12.0*sqrt(48.0 - 120.0*brg + 81.0*brg*brg) + 108.0*brg - 80.0) ;
461 brg = brg/6.0 - 4.0/(3.0*brg) - 1.0/3.0 ;
462 if( brg > 0.0 && brg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(brg)) )*dz ;
463 }
464 /* INFO_message("FWHM z: vz1=%f vz2=%f arg=%f sz=%f countz=%d",vz1,vz2,arg,sz,countz) ; */
465 } else {
466 /* INFO_message("FWHM z: vz1=%f vz2=%f NOT GOOD",vz1,vz2) ; */
467 }
468
469 /* INFO_message("12dif_MAD: sxyz = %g %g %g",sx,sy,sz) ; */
470 LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
471 /* INFO_message("12dif_MAD: fw = %g %g %g",fw_xyz.xyz[0],fw_xyz.xyz[1],fw_xyz.xyz[2]) ; */
472 return fw_xyz ;
473 }
474
475 /*---------------------------------------------------------------------------*/
476
477 static THD_fvec3 (*fester)(MRI_IMAGE *, byte *) = mri_estimate_FWHM_1dif ;
478
mri_fwhm_setfester(THD_fvec3 (* func)(MRI_IMAGE *,byte *))479 void mri_fwhm_setfester( THD_fvec3 (*func)(MRI_IMAGE *, byte *) )
480 {
481 if( func == NULL ) func = mri_estimate_FWHM_1dif ;
482 fester = func ;
483 }
484
485 /*---------------------------------------------------------------------------*/
486 /*! Get FWHM estimates for each sub-brick in a dataset.
487 Output image is 3xN where N=# of sub-bricks.
488 -----------------------------------------------------------------------------*/
489
THD_estimate_FWHM_all(THD_3dim_dataset * dset,byte * mask,int demed,int unif)490 MRI_IMAGE * THD_estimate_FWHM_all( THD_3dim_dataset *dset,
491 byte *mask, int demed , int unif )
492 {
493 int nvals , ii,nvox ;
494 MRI_IMAGE *outim=NULL , *medim=NULL , *madim=NULL ;
495 float *outar,fac , *medar=NULL , *madar=NULL ;
496
497 ENTRY("THD_estimate_FWHM_all") ;
498
499 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
500 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
501
502 nvals = DSET_NVALS(dset) ;
503 outim = mri_new( 3 , nvals , MRI_float ) ;
504 outar = MRI_FLOAT_PTR(outim) ;
505 nvox = DSET_NVOX(dset) ;
506
507 if( unif ){
508 MRI_IMARR *imar ;
509 demed = 1 ;
510 imar = THD_medmad_bricks(dset) ;
511 medim = IMARR_SUBIM(imar,0) ; medar = MRI_FLOAT_PTR(medim) ;
512 madim = IMARR_SUBIM(imar,1) ; madar = MRI_FLOAT_PTR(madim) ;
513 FREE_IMARR(imar) ;
514 for( ii=0 ; ii < nvox ; ii++ )
515 if( madar[ii] > 0.0f ) madar[ii] = 1.0f / madar[ii] ;
516 } else if( demed ){
517 medim = THD_median_brick(dset) ;
518 medar = MRI_FLOAT_PTR(medim) ;
519 }
520
521 AFNI_OMP_START;
522 #pragma omp parallel if( nvals > 4 )
523 { int iv,ii ; MRI_IMAGE *bim ; float *bar ; THD_fvec3 fw ;
524 #pragma omp for
525 for( iv=0 ; iv < nvals ; iv++ ){
526 if( mri_allzero(DSET_BRICK(dset,iv)) ){
527 outar[0+3*iv] = outar[1+3*iv] = outar[2+3*iv] = 0.0f ; continue ;
528 }
529 bim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,iv), DSET_BRICK(dset,iv) );
530 if( demed ){
531 bar = MRI_FLOAT_PTR(bim) ;
532 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] -= medar[ii] ;
533 if( unif )
534 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] *= madar[ii] ;
535 }
536 fw = fester( bim , mask ) ; mri_free(bim) ;
537 UNLOAD_FVEC3( fw , outar[0+3*iv] , outar[1+3*iv] , outar[2+3*iv] ) ;
538 }
539 }
540 AFNI_OMP_END;
541
542 if( demed ) mri_free(medim) ;
543 if( unif ) mri_free(madim) ;
544
545 RETURN(outim) ;
546 }
547
548 /*========================================================================*/
549 static double mom12_stdev_fac = 1.0 ;
550
mri_fwhm_mom12_set_stdev_fac(double fff)551 void mri_fwhm_mom12_set_stdev_fac( double fff ){ mom12_stdev_fac = fff ; }
552
553 /*------------------------------------------------------------------------*/
554 /*! New method using first and second moments of absolute values of
555 differences, instead of only second moments. [05 Aug 2015]
556 *//*----------------------------------------------------------------------*/
557
mri_FWHM_1dif_mom12(MRI_IMAGE * im,byte * mask)558 THD_fvec3 mri_FWHM_1dif_mom12( MRI_IMAGE *im , byte *mask )
559 {
560 int nx; /* number of voxels along x-axis */
561 int ny; /* number of voxels along y-axis */
562 int nz; /* number of voxels along z-axis */
563 int nxy, nxyz; /* total number of voxels */
564 int ixyz; /* voxel index */
565 double dx; /* voxel size along x-axis */
566 double dy; /* voxel size along y-axis */
567 double dz; /* voxel size along z-axis */
568 int ix, jy, kz, ixyz2;
569 double fsum, fsq, var , vfac , arg ;
570 double dfdx, dfdy, dfdz ;
571 double dfdxs1, dfdxs2, dfdxs3, dfdxs4, dfdxna ;
572 double dfdys1, dfdys2, dfdys3, dfdys4, dfdyna ;
573 double dfdzs1, dfdzs2, dfdzs3, dfdzs4, dfdzna ;
574 int count=0, countx=0, county=0, countz=0 ;
575
576 float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
577 MRI_IMAGE *lim ; float *fim ;
578 THD_fvec3 fw_xyz ;
579
580 /*----- initialize local variables -----*/
581
582 LOAD_FVEC3(fw_xyz,sx,sy,sz) ; /* load with error flags */
583
584 if( im == NULL || mri_allzero(im) ) return fw_xyz ;
585 lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
586 fim = MRI_FLOAT_PTR(lim) ;
587 nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;
588
589 /*----- estimate the variance of the data itself -----*/
590
591 fsum = 0.0; fsq = 0.0; count = 0;
592 for (ixyz = 0; ixyz < nxyz; ixyz++){
593 if( GOOD(ixyz) ){ count++; arg = fim[ixyz]; fsum += arg; fsq += arg*arg; }
594 }
595 if( count < 9 || fsq <= 0.0 ){ /* no data? */
596 if( lim != im ) mri_free(lim) ;
597 return fw_xyz ;
598 }
599 var = (fsq - (fsum * fsum)/count) / (count-1.0);
600 if( var <= 0.0 ){ /* crappy data? */
601 if( lim != im ) mri_free(lim) ;
602 return fw_xyz ;
603 }
604 vfac = 1.0 / sqrt(2.0*var) ;
605
606 /*----- estimate moments of the partial derivatives -----*/
607 /*((((( 3rd and 4th moments are not used at present )))))*/
608
609 dfdxs1 = dfdxs2 = dfdxs3 = dfdxs4 = 0.0 ; dfdxna = 0.0 ;
610 dfdys1 = dfdys2 = dfdys3 = dfdys4 = 0.0 ; dfdyna = 0.0 ;
611 dfdzs1 = dfdzs2 = dfdzs3 = dfdzs4 = 0.0 ; dfdzna = 0.0 ;
612 countx = county = countz = 0 ;
613
614 for (ixyz = 0; ixyz < nxyz; ixyz++){
615 if( GOOD(ixyz) ){
616
617 arg = fim[ixyz] ;
618 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
619
620 if( ix+1 < nx ){
621 ixyz2 = ixyz+1 ;
622 if( dontcheckplus || GOOD(ixyz2) ){
623 dfdx = vfac*(fim[ixyz2] - arg) ; dfdxna += dfdx ; dfdx = fabs(dfdx) ;
624 dfdxs1 += dfdx; dfdxs2 += dfdx*dfdx;
625 dfdxs3 += dfdx*dfdx*dfdx; dfdxs4 += dfdx*dfdx*dfdx*dfdx;
626 countx++;
627 }
628 }
629
630 if( jy+1 < ny ){
631 ixyz2 = ixyz+nx ;
632 if( dontcheckplus || GOOD(ixyz2) ){
633 dfdy = vfac*(fim[ixyz2] - arg) ; dfdyna += dfdy ; dfdy = fabs(dfdy) ;
634 dfdys1 += dfdy; dfdys2 += dfdy*dfdy;
635 dfdys3 += dfdy*dfdy*dfdy; dfdys4 += dfdy*dfdy*dfdy*dfdy;
636 county++;
637 }
638 }
639
640 if( kz+1 < nz ){
641 ixyz2 = ixyz+nxy ;
642 if( dontcheckplus || GOOD(ixyz2) ){
643 dfdz = vfac*(fim[ixyz2] - arg) ; dfdzna += dfdz ; dfdz = fabs(dfdz) ;
644 dfdzs1 += dfdz; dfdzs2 += dfdz*dfdz;
645 dfdzs3 += dfdz*dfdz*dfdz; dfdzs4 += dfdz*dfdz*dfdz*dfdz;
646 countz++;
647 }
648 }
649 }
650 }
651 if( countx == 0 ) countx = 1 ; /* patch for possible stupidity in data */
652 if( county == 0 ) county = 1 ;
653 if( county == 0 ) county = 1 ;
654
655 dfdxs1 /= countx; dfdxs2 /= countx; dfdxs3 /= countx; dfdxs4 /= countx;
656 dfdys1 /= county; dfdys2 /= county; dfdys3 /= county; dfdys4 /= county;
657 dfdzs1 /= countz; dfdzs2 /= countz; dfdzs3 /= countz; dfdzs4 /= countz;
658
659 dfdxna /= countx ; dfdyna /= county ; dfdzna /= countz ;
660
661 dx = lim->dx; dy = lim->dy; dz = lim->dz;
662 if( lim != im ) mri_free(lim) ;
663
664 #if 0
665 ININFO_message(" var= %13.6g",var) ;
666 ININFO_message(" dx: %13.6g %13.6g %13.6g %13.6g na = %13.6g #x = %d",
667 dfdxs1 , dfdxs2 , dfdxs3 , dfdxs4 , dfdxna , countx ) ;
668 ININFO_message(" dy: %13.6g %13.6g %13.6g %13.6g na = %13.6g #y = %d",
669 dfdys1 , dfdys2 , dfdys3 , dfdys4 , dfdyna , county ) ;
670 ININFO_message(" dz: %13.6g %13.6g %13.6g %13.6g na = %13.6g #z = %d",
671 dfdzs1 , dfdzs2 , dfdzs3 , dfdzs4 , dfdzna , countz ) ;
672 #endif
673
674 #undef SQHPI
675 #define SQHPI 1.253314 /* sqrt(PI/2) */
676 #undef SQEPI
677 #define SQEPI 0.626657 /* sqrt(PI/8) -- needed for 3rd moment */
678
679 /*--------------------------------------------------------------------------*/
680 /* Model: first differences are distributed like
681 ef = Normal( 0 , 2*V*tau^2 )
682 where tau = sqrt{ 1 - exp[-del^2/(4*s^2)] }
683 V = variance (var) of data
684 del = spacing of data in mm
685 s = sigma of smoothing kernel
686 = del * sqrt( -1/(4*log(1-tau^2))
687 [s(tau) is a monotonic decreasing function of tau]
688 FWHM = sqrt(8*log(2))*s = 2.35482*s
689 and we assume that tau in turn is (independently) random, distributed
690 with density w(tau) over the interval [0,1]. That is, the density
691 function of df = ef/sqrt(2*V) is
692 |1
693 int| w(tau) 1/[tau*sqrt(2*pi)] * exp[-df^2/(2*tau^2)] d(tau)
694 |0
695 Unless w(tau) is a delta function, df will be distributed non-normally.
696 Then the moments of the absolute value of df are related to the moments
697 of tau by |1
698 E[|df|^p] = int| w(tau) E[ |Normal(0,tau^2)|^p ] d(tau)
699 |0
700 = 2^(p/2) * Gamma((p+1)/2) / sqrt(pi) * E[tau^p]
701 E[|df|] = sqrt(2/pi) * E[tau]
702 E[|df|^2] = E[tau^2]
703 E[|df|^3] = sqrt(8/pi) * E[tau^3]
704 E[|df|^4] = 3 * E[tau^4] etc.
705 (At this time, we don't use the 3rd or 4th moments, but someday ...?)
706 The "standard" (Forman) way to estimate FWHN (mri_estimate_FWHM_1dif)
707 assumes s is constant; e.g., that w(tau) = delta(tau-q), in which
708 case E[|df|^2] = q^2, and so only the second moment is needed to
709 estimate tau=q (and thus FWHM). In this routine, we use E[|df|] to
710 estimate E[tau], and then subtract off one standard deviation
711 of tau sqrt(E[tau^2]-E[tau]^2) to weight the estimate towards
712 smaller tau and thus larger FWHM -- since larger FWHM will produce
713 more "false positive" clusters. Completely ad hoc, but it's something. */
714 /*--------------------------------------------------------------------------*/
715
716 { double mn, sg , ff ;
717 mn = SQHPI*dfdxs1 ;
718 sg = dfdxs2 - mn*mn ;
719 if( mn > 0.0 && sg > 0.0 ){
720 sg = sqrt(sg) ;
721 for( ff=mom12_stdev_fac ; sg*ff >= mn ; ff *= 0.888 ) ; /*nada*/
722 mn -= ff*sg ;
723 }
724 mn *= mn ;
725 if( mn > 0.0 && mn < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(1.0-mn)) ) * dx;
726
727 mn = SQHPI*dfdys1 ;
728 sg = dfdys2 - mn*mn ;
729 if( mn > 0.0 && sg > 0.0 ){
730 sg = sqrt(sg) ;
731 for( ff=mom12_stdev_fac ; sg*ff >= mn ; ff *= 0.888 ) ; /*nada*/
732 mn -= ff*sg ;
733 }
734 mn *= mn ;
735 if( mn > 0.0 && mn < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(1.0-mn)) ) * dy;
736
737 mn = SQHPI*dfdzs1 ;
738 sg = dfdzs2 - mn*mn ;
739 if( mn > 0.0 && sg > 0.0 ){
740 sg = sqrt(sg) ;
741 for( ff=mom12_stdev_fac ; sg*ff >= mn ; ff *= 0.888 ) ; /*nada*/
742 mn -= ff*sg ;
743 }
744 mn *= mn ;
745 if( mn > 0.0 && mn < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(1.0-mn)) ) * dz;
746 }
747
748 LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
749 return fw_xyz ;
750 }
751
752 /*========================================================================*/
753 /*========================================================================*/
754 /* Return the cluster of points over which we estimage the ACF.
755 Use the central NCLU_BASE set of points, and then fill out
756 up to NCLU_GOAL from a selection of the outer points,
757 so that we don't end up with a vast cluster and slow estimation.
758 *//*----------------------------------------------------------------------*/
759
760 #undef NCLU_GOAL
761 #define NCLU_GOAL 666
762 #undef NCLU_BASE
763 #define NCLU_BASE 111
764
765 static int do_2D_ACF = 0 ;
766
set_ACF_2D(int nn)767 void set_ACF_2D( int nn ){ do_2D_ACF = nn ; }
768
get_ACF_cluster(float dx,float dy,float dz,float radius)769 static MCW_cluster * get_ACF_cluster( float dx, float dy, float dz, float radius )
770 {
771 MCW_cluster *clout=NULL , *cltemp=NULL ;
772 int pp , dp ;
773
774 /* a ball of desired radius */
775
776 cltemp = MCW_spheremask( dx,dy,dz , radius ) ;
777
778 /* if not too big, just use this */
779
780 if( cltemp->num_pt < NCLU_GOAL ) return cltemp ;
781
782 /* keep the central NCLU_BASE points */
783
784 INIT_CLUSTER(clout) ;
785 for( pp=0 ; clout->num_pt < NCLU_BASE && pp < cltemp->num_pt ; pp++ ){
786 if( do_2D_ACF && cltemp->k[pp] != 0 ) continue ;
787 ADDTO_CLUSTER(clout,cltemp->i[pp],cltemp->j[pp],cltemp->k[pp],0.0f) ;
788 }
789
790 /* keep a selection of the outer points, so we don't get too many */
791
792 dp = (int)rintf( (cltemp->num_pt-NCLU_BASE) / (float)(NCLU_GOAL-NCLU_BASE) ) ;
793 for( ; pp < cltemp->num_pt ; ){
794 if( do_2D_ACF && cltemp->k[pp] != 0 ){
795 pp++ ;
796 } else {
797 ADDTO_CLUSTER(clout,cltemp->i[pp],cltemp->j[pp],cltemp->k[pp],0.0f) ;
798 pp += dp ;
799 }
800 }
801
802 KILL_CLUSTER(cltemp) ;
803 MCW_radsort_cluster(clout,dx,dy,dz) ;
804 return clout ;
805 }
806 /*----------------------------------------------------------------------------*/
807 /* Estimate the shape of the spatial autocorrelation function [09 Nov 2015] */
808
mri_estimate_ACF(MRI_IMAGE * im,byte * mask,MCW_cluster * clout)809 int mri_estimate_ACF( MRI_IMAGE *im , byte *mask , MCW_cluster *clout )
810 {
811 float dx,dy,dz ;
812 MRI_IMAGE *lim ; float *fim ;
813 int nx,ny,nz,nxy,nxyz , ixyz,ip,jp,kp,ppp,qqq , ix,jy,kz ;
814 double fsum, fsq, fbar, fvar , arg ;
815 int count , nclu ;
816 double *acfsqq ;
817 int *nacf ;
818
819 ENTRY("mri_estimate_ACF") ;
820
821 if( im == NULL || mri_allzero(im) ) RETURN(-1) ;
822 if( clout == NULL || clout->num_pt < 5 ) RETURN(-1) ;
823
824 for( ppp=0 ; ppp < clout->num_pt ; ppp++ ) clout->mag[ppp] = 0.0f ;
825
826 /*-- setup to process the input image --*/
827
828 lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
829 fim = MRI_FLOAT_PTR(lim) ;
830 nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;
831
832 /*----- estimate the variance of the data -----*/
833
834 fsum = 0.0; fsq = 0.0; count = 0;
835 for (ixyz = 0; ixyz < nxyz; ixyz++){
836 if( GOOD(ixyz) ){ count++; arg = fim[ixyz]; fsum += arg; fsq += arg*arg; }
837 }
838 if( count < 9 || fsq <= 0.0 ){ /* no data? */
839 if( lim != im ) mri_free(lim) ;
840 RETURN(-1) ;
841 }
842 fbar = fsum / count ;
843 fvar = (fsq - (fsum * fsum)/count) / (count-1.0);
844 if( fvar <= 0.0 ){
845 if( lim != im ) mri_free(lim) ;
846 RETURN(-1) ;
847 }
848
849 /*--- space to accumulate results ---*/
850
851 nclu = clout->num_pt ;
852 acfsqq = (double *)calloc(sizeof(double),nclu) ;
853 nacf = (int *)calloc(sizeof(int) ,nclu) ;
854
855 /*--- loop over voxels, sum up cross products ---*/
856
857 for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
858
859 if( !GOOD(ixyz) ) continue ; /* not in mask */
860
861 arg = fim[ixyz] - fbar ;
862 IJK_TO_THREE(ixyz,ix,jy,kz,nx,nxy) ;
863
864 /* loop over offsets */
865
866 { int ppp , ip,jp,kp,qqq ;
867 for( ppp=1 ; ppp < nclu ; ppp++ ){
868 ip = ix + clout->i[ppp] ; if( ip >= nx || ip < 0 ) continue ;
869 jp = jy + clout->j[ppp] ; if( jp >= ny || jp < 0 ) continue ;
870 kp = kz + clout->k[ppp] ; if( kp >= nz || kp < 0 ) continue ;
871 qqq = THREE_TO_IJK(ip,jp,kp,nx,nxy) ; if( !GOOD(qqq) ) continue ;
872 acfsqq[ppp] += (fim[qqq]-fbar)*arg ; nacf[ppp]++ ;
873 }
874 }
875
876 }
877
878 /*--- load results into clout ---*/
879
880 for( ppp=1 ; ppp < nclu ; ppp++ ){
881 if( nacf[ppp] > 5 )
882 clout->mag[ppp] = (float)( acfsqq[ppp] / ( fvar * (nacf[ppp]-1.0) ) ) ;
883 }
884 clout->mag[0] = 1.0f ;
885
886 free(acfsqq) ; free(nacf) ;
887 if( lim != im ) mri_free(lim) ;
888 RETURN(0) ;
889 }
890
891 /*----------------------------------------------------------------------------*/
892 /*! Get average ACF estimate for all sub-bricks [09 Nov 2015]
893 *//*--------------------------------------------------------------------------*/
894
THD_estimate_ACF(THD_3dim_dataset * dset,byte * mask,int demed,int unif,float radius)895 MCW_cluster * THD_estimate_ACF( THD_3dim_dataset *dset,
896 byte *mask, int demed, int unif, float radius )
897 {
898 int iv , nvals , ii,nvox , nout ;
899 MRI_IMAGE *medim=NULL , *madim=NULL ;
900 float *medar=NULL , *madar=NULL ;
901 MCW_cluster *clout=NULL , **cltemp ;
902
903 ENTRY("THD_estimate_ACF") ;
904
905 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
906 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
907
908 nvals = DSET_NVALS(dset) ;
909 nvox = DSET_NVOX(dset) ;
910
911 if( unif ){
912 MRI_IMARR *imar ;
913 demed = 1 ;
914 imar = THD_medmad_bricks(dset) ;
915 medim = IMARR_SUBIM(imar,0) ; medar = MRI_FLOAT_PTR(medim) ;
916 madim = IMARR_SUBIM(imar,1) ; madar = MRI_FLOAT_PTR(madim) ;
917 FREE_IMARR(imar) ;
918 for( ii=0 ; ii < nvox ; ii++ )
919 if( madar[ii] > 0.0f ) madar[ii] = 1.0f / madar[ii] ;
920 } else if( demed ){
921 medim = THD_median_brick(dset) ;
922 medar = MRI_FLOAT_PTR(medim) ;
923 }
924
925 set_ACF_2D( DSET_NZ(dset)==1 ) ;
926
927 clout = get_ACF_cluster( fabsf(DSET_DX(dset)) , fabsf(DSET_DY(dset)) ,
928 fabsf(DSET_DZ(dset)) , radius ) ;
929 if( clout == NULL ){
930 if( medim != NULL ) mri_free(medim) ;
931 if( madim != NULL ) mri_free(madim) ;
932 RETURN(NULL) ;
933 }
934
935 for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] = 0.0f ;
936 cltemp = (MCW_cluster **)malloc(sizeof(MCW_cluster *)*nvals) ;
937
938 AFNI_OMP_START;
939 #pragma omp parallel if( nvals > 4 )
940 { int iv,gg ; MRI_IMAGE *bim ;
941 #pragma omp for
942 for( iv=0 ; iv < nvals ; iv++ ){
943 bim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,iv), DSET_BRICK(dset,iv) );
944 bim->dx = DSET_DX(dset) ; bim->dy = DSET_DY(dset) ; bim->dz = DSET_DZ(dset) ;
945 if( demed ){
946 float *bar = MRI_FLOAT_PTR(bim) ; int ii ;
947 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] -= medar[ii] ;
948 if( unif )
949 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] *= madar[ii] ;
950 }
951 COPY_CLUSTER(cltemp[iv],clout) ;
952 gg = mri_estimate_ACF( bim , mask , cltemp[iv] ) ;
953 if( gg < 0 ){
954 KILL_CLUSTER(cltemp[iv]) ; cltemp[iv] = NULL ;
955 }
956 }
957 }
958 AFNI_OMP_END;
959
960 for( nout=iv=0 ; iv < nvals ; iv++ ){
961 if( cltemp[iv] != NULL ){
962 for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] += cltemp[iv]->mag[ii] ;
963 KILL_CLUSTER(cltemp[iv]) ; nout++ ;
964 }
965 }
966 free(cltemp) ;
967
968 if( medim != NULL ) mri_free(medim) ;
969 if( madim != NULL ) mri_free(madim) ;
970
971 if( nout > 1 ){
972 for( ii=0 ; ii < clout->num_pt ; ii++ )
973 clout->mag[ii] /= nout ;
974 } else if( nout == 0 ){
975 ERROR_message("ACF estimation fails :-(") ;
976 KILL_CLUSTER(clout) ; clout = NULL ;
977 }
978
979 RETURN(clout) ;
980 }
981
982 /*----------------------------------------------------------------------------*/
983 /*! Get average ACF FWHM estimate for all images in the array [30 Dec 2015]
984 *//*--------------------------------------------------------------------------*/
985
mriarr_estimate_FWHM_acf(MRI_IMARR * imar,byte * mask,int unif,float radius)986 float mriarr_estimate_FWHM_acf( MRI_IMARR *imar, byte *mask, int unif, float radius )
987 {
988 int iv , nvals , ii,nvox , nout , demed=0 ;
989 MRI_IMAGE *medim=NULL , *madim=NULL ;
990 float *medar=NULL , *madar=NULL , dx,dy,dz , fwhm_out=0.0f ;
991 MCW_cluster *clout=NULL , **cltemp ;
992
993 ENTRY("mriarr_estimate_FWHM_acf") ;
994
995 if( imar == NULL ) RETURN(0.0f) ;
996
997 nvals = IMARR_COUNT(imar) ; if( nvals < 2 ) unif = 0 ;
998 nvox = IMARR_SUBIM(imar,0)->nvox ;
999 dx = IMARR_SUBIM(imar,0)->dx ;
1000 dy = IMARR_SUBIM(imar,0)->dy ;
1001 dz = IMARR_SUBIM(imar,0)->dz ;
1002
1003 if( unif ){
1004 MRI_IMARR *qmar ;
1005 demed = 1 ;
1006 qmar = IMARR_medmad_bricks(imar) ;
1007 medim = IMARR_SUBIM(qmar,0) ; medar = MRI_FLOAT_PTR(medim) ;
1008 madim = IMARR_SUBIM(qmar,1) ; madar = MRI_FLOAT_PTR(madim) ;
1009 FREE_IMARR(qmar) ;
1010 for( ii=0 ; ii < nvox ; ii++ )
1011 if( madar[ii] > 0.0f ) madar[ii] = 1.0f / madar[ii] ;
1012 }
1013
1014 set_ACF_2D( IMARR_SUBIM(imar,0)->nz == 1 ) ;
1015
1016 clout = get_ACF_cluster( dx,dy,dz, radius ) ;
1017 if( clout == NULL ){
1018 if( medim != NULL ) mri_free(medim) ;
1019 if( madim != NULL ) mri_free(madim) ;
1020 RETURN(0.0f) ;
1021 }
1022
1023 for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] = 0.0f ;
1024 cltemp = (MCW_cluster **)malloc(sizeof(MCW_cluster *)*nvals) ;
1025
1026 AFNI_OMP_START;
1027 #pragma omp parallel if( nvals > 4 )
1028 { int iv,gg ; MRI_IMAGE *bim ;
1029 #pragma omp for
1030 for( iv=0 ; iv < nvals ; iv++ ){
1031 bim = mri_to_float( IMARR_SUBIM(imar,iv) ) ;
1032 bim->dx = dx ; bim->dy = dy ; bim->dz = dz ;
1033 if( demed ){
1034 float *bar = MRI_FLOAT_PTR(bim) ; int ii ;
1035 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] -= medar[ii] ;
1036 if( unif )
1037 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] *= madar[ii] ;
1038 }
1039 COPY_CLUSTER(cltemp[iv],clout) ;
1040 gg = mri_estimate_ACF( bim , mask , cltemp[iv] ) ;
1041 if( gg < 0 ){
1042 KILL_CLUSTER(cltemp[iv]) ; cltemp[iv] = NULL ;
1043 }
1044 }
1045 }
1046 AFNI_OMP_END;
1047
1048 for( nout=iv=0 ; iv < nvals ; iv++ ){
1049 if( cltemp[iv] != NULL ){
1050 for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] += cltemp[iv]->mag[ii] ;
1051 KILL_CLUSTER(cltemp[iv]) ; nout++ ;
1052 }
1053 }
1054 free(cltemp) ;
1055
1056 if( medim != NULL ) mri_free(medim) ;
1057 if( madim != NULL ) mri_free(madim) ;
1058
1059 if( nout > 1 ){
1060 for( ii=0 ; ii < clout->num_pt ; ii++ )
1061 clout->mag[ii] /= nout ;
1062 } else if( nout == 0 ){
1063 ERROR_message("ACF estimation fails :-(") ;
1064 KILL_CLUSTER(clout) ; clout = NULL ;
1065 }
1066
1067 if( clout != NULL ){
1068 float_quad acf_Epar ;
1069 acf_Epar = ACF_cluster_to_modelE( clout , dx,dy,dz ) ;
1070 fwhm_out = acf_Epar.d ;
1071 KILL_CLUSTER(clout) ; clout = NULL ;
1072 }
1073
1074 RETURN(fwhm_out) ;
1075 }
1076
1077 /*----------------------------------------------------------------------------*/
1078 /* Function for Powell minimization for ACF modelE fit. */
1079
1080 static int nar=0 ;
1081 static float *aar=NULL , *rar=NULL ;
1082
ACF_modelE_costfunc(int npar,double * par)1083 double ACF_modelE_costfunc( int npar , double *par )
1084 {
1085 double aa=par[0] , bb=par[1] , cc=par[2] ;
1086 double a1=1.0-aa , bq=0.5/(bb*bb) , ci=1.0/cc ;
1087 double sum=0.0 , fval ;
1088 int ii ;
1089
1090 for( ii=0 ; ii < nar ; ii++ ){
1091 fval = aa * exp(-bq*rar[ii]*rar[ii]) + a1 * exp(-ci*rar[ii]) - aar[ii] ;
1092 sum += fval * fval ;
1093 }
1094
1095 return sum ;
1096 }
1097
1098 /*----------------------------------------------------------------------------*/
1099
1100 static MRI_IMAGE *ACF_im=NULL ;
ACF_get_1D(void)1101 MRI_IMAGE * ACF_get_1D(void){ return ACF_im ; }
1102
1103 /*----------------------------------------------------------------------------*/
1104 /*! Convert ACF cluster output to function vs radius [09 Nov 2015]
1105 Exponential model function vs. radius r is
1106 ACF(r) = a*exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)
1107 where a,b,c are the estimated parameters returned in the float_quad.
1108 With restrictions that 0 <= a <= 1 ; b > 0 ; c > 0
1109 The last parameter return (d) is the FWHM of the curve fit.
1110 *//*--------------------------------------------------------------------------*/
1111
ACF_cluster_to_modelE(MCW_cluster * acf,float dx,float dy,float dz)1112 float_quad ACF_cluster_to_modelE( MCW_cluster *acf, float dx, float dy, float dz )
1113 {
1114 int nclu, nr, pp, qq, ss, dq ;
1115 float xx , yy , zz , rmax , dr , vp,rp , vs ;
1116 float apar , bpar , cpar , dpar , sig , *acar ;
1117 float_quad fqout = { -1.0f , -1.0f , -1.0f , -1.0f } ;
1118 double xpar[3] , xbot[3] , xtop[3] ;
1119
1120 ENTRY("ACF_cluster_to_modelE") ;
1121
1122 if( ACF_im != NULL ){ mri_free(ACF_im) ; ACF_im = NULL ; }
1123
1124 if( acf == NULL ) RETURN(fqout) ;
1125 nclu = acf->num_pt ; if( nclu < 5 ) RETURN(fqout) ;
1126
1127 if( dx <= 0.0f ) dx = 1.0f ;
1128 if( dy <= 0.0f ) dy = 1.0f ;
1129 if( dz <= 0.0f ) dz = 1.0f ;
1130
1131 /* load rar and aar arrays with radius and ACF values */
1132
1133 #undef CRAD
1134 #define CRAD(a) ( xx = dx * abs(acf->i[a]) , \
1135 yy = dy * abs(acf->j[a]) , \
1136 zz = dz * abs(acf->k[a]) , sqrt(xx*xx+yy*yy+zz*zz) )
1137
1138 aar = (float *)malloc(sizeof(float)*nclu) ;
1139 rar = (float *)malloc(sizeof(float)*nclu) ;
1140
1141 for( pp=0 ; pp < nclu ; pp++ ){
1142 rar[pp] = CRAD(pp) ;
1143 aar[pp] = acf->mag[pp] ;
1144 }
1145
1146 qsort_floatfloat( nclu , rar , aar ) ; /* to ensure increasing rar */
1147
1148 /* collapse data with very similar rar values */
1149
1150 for( pp=0 ; pp < nclu-1 ; pp++ ){
1151 vp = aar[pp] ; rp = rar[pp]*1.01f ;
1152 for( qq=pp+1 ; qq < nclu && rar[qq] <= rp ; qq++ ) ; /*nada*/
1153 if( qq > pp+1 ){
1154 for( vs=0.0f,ss=pp ; ss < qq ; ss++ ) vs += aar[ss] ;
1155 aar[pp] = vs/(qq-pp) ;
1156 dq = qq-(pp+1) ;
1157 for( ss=qq ; ss < nclu ; ss++ ){
1158 aar[ss-dq] = aar[ss]; rar[ss-dq] = rar[ss];
1159 }
1160 nclu -= dq ;
1161 }
1162 }
1163
1164 /* crude estimates for model parameters */
1165
1166 nar = nclu ;
1167
1168 for( pp=0 ; pp < nclu && aar[pp] > 0.5f ; pp++ ) ; /*nada*/
1169 if( pp == nclu ){
1170 ERROR_message("largest ACF found is %g -- too big for model fit!",aar[nclu-1]) ;
1171 free(aar) ; free(rar) ; aar = rar = NULL ; nar = 0 ;
1172 RETURN(fqout) ;
1173 }
1174 apar = 1.0 / logf(aar[pp]) ;
1175 bpar = sqrtf(-0.5*apar) * rar[pp] ;
1176 cpar = -rar[pp] * apar;
1177 apar = 0.5f ;
1178
1179 /* nonlinear optimization of parameters */
1180
1181 xpar[0] = apar ; xpar[1] = bpar ; xpar[2] = cpar ;
1182 xbot[0] = 0.006 ; xbot[1] = 0.05*bpar ; xbot[2] = 0.05*cpar ;
1183 xtop[0] = 0.994 ; xtop[1] = 5.55*bpar ; xtop[2] = 5.55*cpar ;
1184
1185 #if 0
1186 pp = powell_newuoa_con( 3 , xpar , xbot , xtop ,
1187 99 , 0.05 , 0.0005 , 999 , ACF_modelE_costfunc ) ;
1188 #else
1189 pp = powell_newuoa_constrained( 3 , xpar , NULL , xbot , xtop ,
1190 666 , 44 , 9 ,
1191 0.05 , 0.0005 , 999 , ACF_modelE_costfunc ) ;
1192 #endif
1193
1194 if( pp < 0 ){
1195 ERROR_message("optimization of ACF model fit fails :-(") ;
1196 free(aar) ; free(rar) ; aar = rar = NULL ; nar = 0 ;
1197 RETURN(fqout) ;
1198 }
1199
1200 ACF_im = mri_new( nar , 4 , MRI_float ) ;
1201 acar = MRI_FLOAT_PTR(ACF_im) ;
1202 apar = xpar[0] ; bpar = xpar[1] ; cpar = xpar[2] ;
1203 for( pp=0 ; pp < nar ; pp++ ){
1204 vs = apar * expf( -0.5f*rar[pp]*rar[pp]/(bpar*bpar) )
1205 + (1.0f-apar) * expf( -rar[pp]/cpar ) ;
1206 acar[pp] = rar[pp] ; /* radius */
1207 acar[pp+nar] = aar[pp] ; /* empirical ACF */
1208 acar[pp+2*nar] = vs ; /* model ACF */
1209 }
1210
1211 { floatvec *fv ;
1212 MAKE_floatvec(fv,nar) ;
1213 for( pp=0 ; pp < nar ; pp++ ) fv->ar[pp] = acar[pp+2*nar] ;
1214 vs = interp_inverse_floatvec( fv , 0.5f ) ;
1215 for( pp=0 ; pp < nar ; pp++ ) fv->ar[pp] = rar[pp] ;
1216 dpar = 2.0*interp_floatvec( fv , vs ) ; sig = FWHM_TO_SIGMA(dpar) ;
1217 KILL_floatvec(fv) ;
1218 }
1219
1220 for( pp=0 ; pp < nar ; pp++ ){ /* add in Gaussian ACF for fun */
1221 acar[pp+3*nar] = expf(-0.5f*rar[pp]*rar[pp]/(sig*sig)) ;
1222 }
1223
1224 free(aar) ; free(rar) ; aar = rar = NULL ; nar = 0 ;
1225
1226 fqout.a = xpar[0] ; fqout.b = xpar[1] ; fqout.c = xpar[2] ; fqout.d = dpar ;
1227 RETURN(fqout) ;
1228 }
1229