1 #include "mrilib.h"
2
3 /** if using OpenMP, this file should be #include-d into the main program! **/
4
5 #ifdef USE_OMP
6 #include <omp.h>
7 #include "cs_qmed.c"
8 #endif
9
10 /* default fill and unfill values */
11 static float fillvalue = 1.0;
12 static float unfillvalue = 1.0;
13 static float maskvalue = -1.0;
14 static float maskvalue2 = -2.0;
15
16 /*--------------------------------------------------------------------------*/
17 /*! Input = 1D float array, and an NSTAT_ code to compute some statistic.
18 Output = statistic's value.
19 *//*------------------------------------------------------------------------*/
20
mri_nstat(int code,int npt,float * far,float voxval,MCW_cluster * nbhd)21 float mri_nstat( int code , int npt , float *far , float voxval, MCW_cluster *nbhd )
22 {
23 register float outval ; float val ;
24
25 outval = 0.0f ;
26
27 if( npt <= 0 || far == NULL ) return outval ;
28
29 switch( code ){
30
31 case NSTAT_NUM: outval = (float)npt ; break ; /* quite easy */
32
33 /* another easy case - filled */
34 case NSTAT_FILLED:
35 if(npt==nbhd->num_pt)
36 outval = fillvalue;
37 break ;
38
39 /* and one more simple case - not filled */
40 case NSTAT_UNFILLED:
41 if(npt<nbhd->num_pt)
42 outval = unfillvalue;
43 break ;
44
45 /* check if neighborhood contains a specified mask value */
46 case NSTAT_MASKED:{
47 register int ii;
48 outval = 0.0f;
49 for( ii=0; ii < npt; ii++) {
50 if(far[ii] == maskvalue) {
51 outval = unfillvalue;
52 break;
53 }
54 }
55 }
56 break ;
57
58 /* check if neighborhood contains another specified mask value */
59 case NSTAT_MASKED2: {
60 register int ii;
61 outval = 0.0f;
62 for( ii=0; ii < npt; ii++) {
63 if(far[ii] == maskvalue2) {
64 outval = unfillvalue;
65 break;
66 }
67 }
68 }
69 break ;
70
71 default:
72 case NSTAT_SUM:
73 case NSTAT_MEAN:{
74 register int ii ;
75 for( ii=0 ; ii < npt ; ii++ ) outval += far[ii] ;
76 if( code != NSTAT_SUM ) outval /= npt ;
77 }
78 break ;
79
80 case NSTAT_FNZNUM:
81 case NSTAT_NZNUM:{
82 register int ii ;
83 for( ii=0 ; ii < npt ; ii++ ) if (far[ii] != 0.0f) outval += 1 ;
84 if( code != NSTAT_NZNUM) outval /= npt ; /* fractional number of points that are non-zero */
85 }
86 break ;
87
88
89 case NSTAT_SIGMA: /* these 3 need the mean and variance sums */
90 case NSTAT_CVAR:
91 case NSTAT_VAR:{
92 register float mm,vv ; register int ii ;
93 if( npt == 1 ) break ; /* will return 0.0 */
94 for( mm=0.0,ii=0 ; ii < npt ; ii++ ) mm += far[ii] ;
95 mm /= npt ;
96 for( vv=0.0,ii=0 ; ii < npt ; ii++ ) vv += (far[ii]-mm)*(far[ii]-mm) ;
97 vv /= (npt-1) ;
98 if( code == NSTAT_SIGMA ) outval = sqrt(vv) ;
99 else if( code == NSTAT_VAR ) outval = vv ;
100 else if( mm != 0.0f ) outval = sqrt(vv) / fabsf(mm) ;
101 }
102 break ;
103
104 case NSTAT_MEDIAN:
105 qmedmad_float( npt , far , &val , NULL ) ; outval = val ;
106 break ;
107
108 case NSTAT_MAD:
109 qmedmad_float( npt , far , NULL , &val ) ; outval = val ;
110 break ;
111
112 case NSTAT_MODE:
113 outval = qmode_float( npt , far);
114 break ;
115
116 case NSTAT_NZMODE:
117 outval = qnzmode_float( npt , far);
118 break ;
119
120 case NSTAT_P2SKEW:
121 /* Pearson's second skewness coefficient */
122 {
123 register float mm,vv, sig, mean; register int ii ;
124 if( npt == 1 ) break ; /* will return 0.0 */
125 for( mm=0.0,ii=0 ; ii < npt ; ii++ ) mm += far[ii] ;
126 mm /= npt ; mean = mm;
127 for( vv=0.0,ii=0 ; ii < npt ; ii++ )
128 vv += (far[ii]-mm)*(far[ii]-mm) ;
129 vv /= (npt-1) ;
130 sig = sqrt(vv) ;
131 if (sig) {
132 qmedmad_float( npt , far , &val , NULL ) ;
133 outval = 3.0 * (mean - val) / sig;
134 } else outval = 0.0;
135 }
136 break ;
137
138 case NSTAT_KURT:
139 /* Kurtosis estimate, unbiased under normality condition */
140 {
141 register double mm,vv,vv2,pp, sig; register int ii ;
142 if( npt < 4 ) break ; /* will return 0.0 */
143 for( mm=0.0,ii=0 ; ii < npt ; ii++ ) mm += far[ii] ;
144 mm /= npt ;
145 for( vv=0.0,vv2=0.0,ii=0 ; ii < npt ; ii++ ) {
146 pp = (far[ii]-mm)*(far[ii]-mm) ;
147 vv += pp; vv2 += pp*pp;
148 }
149 if (vv != 0.0f) {
150 ii = npt-1;
151 vv2 = (vv2/(vv*vv))*(npt+1.0)*npt*ii - 3.0*ii*ii;
152 outval = (float)(vv2/((npt-2)*(npt-3)));
153 } else outval = 0.0;
154 }
155 break ;
156
157 case NSTAT_MAX:{
158 register int ii ;
159 outval = far[0] ;
160 for( ii=1 ; ii < npt ; ii++ ) if( far[ii] > outval ) outval = far[ii] ;
161 }
162 break ;
163
164 case NSTAT_MIN:{
165 register int ii ;
166 outval = far[0] ;
167 for( ii=1 ; ii < npt ; ii++ ) if( far[ii] < outval ) outval = far[ii] ;
168 }
169 break ;
170
171 case NSTAT_ABSMAX:{
172 register int ii ; register float vv ;
173 outval = fabsf(far[0]) ;
174 for( ii=1 ; ii < npt ; ii++ ){
175 vv = fabsf(far[ii]) ; if( vv > outval ) outval = vv ;
176 }
177 }
178 break ;
179
180 case NSTAT_RANK:{
181 register int ii ;
182 qsort_float(npt, far);
183 outval = 1.0 ;
184 for( ii=1 ; ii < npt ; ii++ ){
185 if (voxval > far[ii]) outval = ii;
186 else break ;
187 }
188 }
189 break ;
190
191 case NSTAT_FRANK:{
192 register int ii ;
193 outval = 1.0 ;
194 if (npt) {
195 qsort_float(npt, far);
196 for( ii=1 ; ii < npt ; ii++ ){
197 if (voxval > far[ii]) outval = ii;
198 else break ;
199 }
200 outval /= npt;
201 }
202 }
203 break ;
204
205 }
206
207 return outval ;
208 }
209
210 /*--------------------------------------------------------------------------*/
211 /*!
212 A specialized function for speeding up computations for segmentation
213 computes 5 statistic values: mean, median, sigma, mad, and skew and
214 stores them in fv5
215 */
216 /*------------------------------------------------------------------------*/
217
mri_nstat_mMP2S(int npt,float * far,float voxval,float * fv5)218 int mri_nstat_mMP2S( int npt , float *far , float voxval, float *fv5)
219 {
220 /* fv5[5]={mean, median, sigma, MAD, skew} */
221 register float mm,vv;
222 register int ii ;
223
224 fv5[0] = fv5[1] = fv5[2] = fv5[3] = fv5[4] = 0.0;
225 if( npt <= 0 || far == NULL ) return 0 ;
226 if ( npt == 1 ) {
227 fv5[0] = fv5[1] = voxval ;
228 return 1;
229 }
230
231 for( mm=0.0,ii=0 ; ii < npt ; ii++ ) mm += far[ii] ;
232 mm /= npt ; fv5[0] = mm;
233 for( vv=0.0,ii=0 ; ii < npt ; ii++ )
234 vv += (far[ii]-mm)*(far[ii]-mm) ;
235 vv /= (npt-1) ;
236 fv5[2] = sqrt(vv) ;
237 if (fv5[2]) {
238 qmedmad_float( npt , far , fv5+1 , fv5+3 ) ;
239 fv5[4] = 3.0 * (fv5[0] - fv5[1]) / fv5[2];
240 } else fv5[4] = 0.0;
241
242
243 return 1 ;
244 }
245
246 /* Compute differences fom the central value far[0] */
mri_nstat_diffs(int npt,float * far,float * fv6,int doabs)247 int mri_nstat_diffs( int npt , float *far , float *fv6, int doabs)
248 {
249 /* fv6[6]={average_diff, min_diff, max_diff,
250 average_adiff, min_adiff, max_adiff} */
251 register float mm,vv, vvmin, vvmax;
252 register int ii ;
253
254 fv6[0] = fv6[1] = fv6[2] = fv6[3] = fv6[4] = fv6[5] = 0.0;
255 if( npt <= 0 || far == NULL ) return 0 ;
256 if ( npt == 1 ) { /* Nothing to do, return quietly though */
257 return 1;
258 }
259
260 if (doabs==0) {
261 vv = (far[1]-far[0]);
262 mm = vvmin = vvmax = vv;
263 for( ii=2 ; ii < npt ; ii++ ) {
264 vv = (far[ii]-far[0]);
265 if (vv < vvmin) {
266 vvmin = vv;
267 } else if (vv > vvmax) {
268 vvmax = vv;
269 }
270 mm += vv ;
271 }
272 mm /= (npt-1) ;
273 fv6[0] = mm;
274 fv6[1] = vvmin;
275 fv6[2] = vvmax;
276 } else if (doabs==1) {
277 vv = ABS(far[1]-far[0]);
278 mm = vvmin = vvmax = vv;
279 for( ii=2 ; ii < npt ; ii++ ) {
280 vv = ABS(far[ii]-far[0]);
281 if (vv < vvmin) {
282 vvmin = vv;
283 } else if (vv > vvmax) {
284 vvmax = vv;
285 }
286 mm += vv ;
287 }
288 mm /= (npt-1) ;
289 fv6[0] = mm;
290 fv6[1] = vvmin;
291 fv6[2] = vvmax;
292 } else {
293 register float mma,vva, vvmina, vvmaxa;
294 vv = (far[1]-far[0]); vva = ABS(vv);
295 mm = vvmin = vvmax = vv;
296 mma = vvmina = vvmaxa = vva;
297 for( ii=2 ; ii < npt ; ii++ ) {
298 vv = (far[ii]-far[0]); vva = ABS(vv);
299 if (vv < vvmin) {
300 vvmin = vv;
301 } else if (vv > vvmax) {
302 vvmax = vv;
303 }
304 mm += vv ;
305 if (vva < vvmina) {
306 vvmina = vva;
307 } else if (vva > vvmaxa) {
308 vvmaxa = vva;
309 }
310 mma += vva ;
311 }
312 mm /= (npt-1) ;
313 mma /= (npt-1) ;
314 fv6[0] = mm;
315 fv6[1] = vvmin;
316 fv6[2] = vvmax;
317 fv6[3] = mma;
318 fv6[4] = vvmina;
319 fv6[5] = vvmaxa;
320 }
321
322 return 1 ;
323 }
324
325 /*--------------------------------------------------------------------------*/
326
327 #if 0
328 static int fwhm_use_variance = 1 ;
329 void mri_nstat_fwhmxyz_usevar( int i ){ fwhm_use_variance = i; }
330 #endif
331
332 #undef INMASK
333 #define INMASK(i) (mask==NULL || mask[i])
334
335 /*--------------------------------------------------------------------------*/
336 /*! FWHM parameters in a neigbhorhood of a point. */
337
mri_nstat_fwhmxyz(int xx,int yy,int zz,MRI_IMAGE * im,byte * mask,MCW_cluster * nbhd)338 THD_fvec3 mri_nstat_fwhmxyz( int xx, int yy, int zz,
339 MRI_IMAGE *im, byte *mask, MCW_cluster *nbhd )
340 {
341 float *far ;
342 int npt , nx,ny,nz,nxy , aa,bb,cc, kk,ii,pp ;
343 THD_fvec3 fw_xyz ;
344 double fsum, fsq, var , arg ;
345 double dfdx, dfdxsum, dfdxsq, varxx;
346 double dfdy, dfdysum, dfdysq, varyy;
347 double dfdz, dfdzsum, dfdzsq, varzz;
348 double dx,dy,dz ;
349 float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
350 int count, countx, county, countz;
351
352 LOAD_FVEC3(fw_xyz,-1,-1,-1) ; /* load with bad values */
353
354 if( im == NULL || im->kind != MRI_float || nbhd == NULL
355 || nbhd->num_pt < 19) return fw_xyz;
356
357 far = MRI_FLOAT_PTR(im) ;
358 nx = im->nx; ny = im->ny; nz = im->nz; nxy = nx*ny; npt = nbhd->num_pt;
359 kk = xx + yy*nx + zz*nxy ;
360 if( npt < 6 || kk < 0 || kk >= nxy*nz || !INMASK(kk) ) return fw_xyz ;
361
362 /*----- estimate the variance of the local data -----*/
363
364 fsum = 0.0; fsq = 0.0; count = 0 ;
365 for( ii=0 ; ii < npt ; ii++ ){
366 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
367 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
368 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
369 kk = aa + bb*nx + cc*nxy ;
370 if( INMASK(kk) ){
371 count++; arg = far[kk]; fsum += arg; fsq += arg*arg;
372 }
373 }
374 if( count < 6 || fsq <= 0.0 ) return fw_xyz ;
375 var = (fsq - (fsum * fsum)/count) / (count-1.0);
376 if( var <= 0.0 ) return fw_xyz ;
377
378 /*----- estimate the partial derivatives -----*/
379
380 dfdxsum = 0.0; dfdysum = 0.0; dfdzsum = 0.0;
381 dfdxsq = 0.0; dfdysq = 0.0; dfdzsq = 0.0;
382 countx = 0; county = 0; countz = 0;
383 for( ii=0 ; ii < npt ; ii++ ){
384 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
385 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
386 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
387 kk = aa + bb*nx + cc*nxy ; if( !INMASK(kk) ) continue ;
388 arg = far[kk] ;
389 if( aa+1 < nx ){
390 pp = kk+1 ;
391 if( INMASK(pp) ){
392 dfdx = (far[pp] - arg) ;
393 dfdxsum += dfdx; dfdxsq += dfdx * dfdx; countx++ ;
394 }
395 }
396 if( bb+1 < ny ){
397 pp = kk+nx ;
398 if( INMASK(pp) ){
399 dfdy = (far[pp] - arg) ;
400 dfdysum += dfdy; dfdysq += dfdy * dfdy; county++ ;
401 }
402 }
403 if( cc+1 < nz ){
404 pp = kk+nxy ;
405 if( INMASK(pp) ){
406 dfdz = (far[pp] - arg) ;
407 dfdzsum += dfdz; dfdzsq += dfdz * dfdz; countz++ ;
408 }
409 }
410 }
411
412 /*----- estimate the variance of the partial derivatives -----*/
413
414 varxx = (countx < 6) ? 0.0
415 : (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1.0);
416
417 varyy = (county < 6) ? 0.0
418 : (dfdysq - (dfdysum * dfdysum)/county) / (county-1.0);
419
420 varzz = (countz < 6) ? 0.0
421 : (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1.0);
422
423 /*---- now estimate the FWHMs ----*/
424 /*---- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor ----*/
425
426 dx = im->dx; dy = im->dy; dz = im->dz;
427
428 arg = 1.0 - 0.5*(varxx/var);
429 if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0/(4.0*log(arg)) )*dx;
430
431 arg = 1.0 - 0.5*(varyy/var);
432 if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0/(4.0*log(arg)) )*dy;
433
434 arg = 1.0 - 0.5*(varzz/var);
435 if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0/(4.0*log(arg)) )*dz;
436
437 LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
438 return fw_xyz ;
439 }
440
441 /*--------------------------------------------------------------------------*/
442
mri_nstat_fwhmbar(int xx,int yy,int zz,MRI_IMAGE * im,byte * mask,MCW_cluster * nbhd)443 float mri_nstat_fwhmbar( int xx, int yy, int zz,
444 MRI_IMAGE *im, byte *mask, MCW_cluster *nbhd )
445 {
446 THD_fvec3 fw ;
447 float fx,fy,fz , sum ; int nsum ;
448
449 fw = mri_nstat_fwhmxyz( xx,yy,zz , im,mask,nbhd ) ;
450 UNLOAD_FVEC3(fw,fx,fy,fz) ;
451
452 sum = 0.0f ; nsum = 0 ;
453 if( fx > 0.0f ){ sum += fx ; nsum++ ; }
454 if( fy > 0.0f ){ sum += fy ; nsum++ ; }
455 if( fz > 0.0f ){ sum += fz ; nsum++ ; }
456 if( nsum > 0 ) sum /= nsum ;
457 return sum ;
458 }
459
460 /*--------------------------------------------------------------------------*/
461 /*! FWHM parameters in a neigbhorhood of a point -- another way. */
462
mri_nstat_fwhmxyz_12dif(int xx,int yy,int zz,MRI_IMAGE * im,byte * mask,MCW_cluster * nbhd,float * ws)463 THD_fvec3 mri_nstat_fwhmxyz_12dif( int xx, int yy, int zz,
464 MRI_IMAGE *im, byte *mask, MCW_cluster *nbhd,
465 float *ws )
466 {
467 int nx; /* number of voxels along x-axis */
468 int ny; /* number of voxels along y-axis */
469 int nz; /* number of voxels along z-axis */
470 int nxy, nxyz; /* total number of voxels */
471 float dx; /* voxel size along x-axis */
472 float dy; /* voxel size along y-axis */
473 float dz; /* voxel size along z-axis */
474 int aa,bb,cc, ii,kk, qm,qp ;
475 float vx1,vy1,vz1 , arg ;
476 float vx2,vy2,vz2 ;
477 int countx, county, countz , npt ;
478 float *dx1ar,*dy1ar,*dz1ar , *dx2ar,*dy2ar,*dz2ar , *fim ;
479 float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
480 THD_fvec3 fw_xyz ;
481
482 /*----- initialize local variables -----*/
483
484 LOAD_FVEC3(fw_xyz,sx,sy,sz) ; /* load with error flags */
485
486 if( im == NULL || im->kind != MRI_float || nbhd == NULL ) return fw_xyz ;
487 nx = im->nx; ny = im->ny; nz = im->nz; nxy = nx*ny; nxyz = nxy*nz;
488 kk = xx + yy*nx + zz*nxy ; npt = nbhd->num_pt ;
489 if( npt < 9 || kk < 0 || kk >= nxyz || !INMASK(kk) ) return fw_xyz ;
490 fim = MRI_FLOAT_PTR(im) ;
491
492 if( ws == NULL ){
493 #pragma omp critical (MALLOC)
494 { dx1ar = (float *)malloc(sizeof(float)*npt) ;
495 dy1ar = (float *)malloc(sizeof(float)*npt) ;
496 dz1ar = (float *)malloc(sizeof(float)*npt) ;
497 dx2ar = (float *)malloc(sizeof(float)*npt) ;
498 dy2ar = (float *)malloc(sizeof(float)*npt) ;
499 dz2ar = (float *)malloc(sizeof(float)*npt) ;
500 }
501 } else {
502 dx1ar = ws + 0*npt ; dy1ar = ws + 1*npt ;
503 dz1ar = ws + 2*npt ; dx2ar = ws + 3*npt ;
504 dy2ar = ws + 4*npt ; dz2ar = ws + 5*npt ;
505 }
506
507 /*----- loop over voxels, compute differences, sum and sum squares -----*/
508
509 countx = county = countz = 0 ;
510 for( ii=0 ; ii < npt ; ii++ ){
511 aa = xx + nbhd->i[ii] ; if( aa < 0 || aa >= nx ) continue ;
512 bb = yy + nbhd->j[ii] ; if( bb < 0 || bb >= ny ) continue ;
513 cc = zz + nbhd->k[ii] ; if( cc < 0 || cc >= nz ) continue ;
514 kk = aa + bb*nx + cc*nxy ; if( !INMASK(kk) ) continue ;
515
516 arg = fim[kk] ;
517
518 if( aa-1 >= 0 && aa+1 < nx ){
519 qp = kk+1 ; qm = kk-1 ;
520 if( INMASK(qp) && INMASK(qm) ){
521 dx1ar[countx] = fim[qp]-arg ; dx2ar[countx] = fim[qp]-fim[qm] ; countx++ ;
522 }
523 }
524
525 if( bb-1 >= 0 && bb+1 < ny ){
526 qp = kk+nx ; qm = kk-nx ;
527 if( INMASK(qp) && INMASK(qm) ){
528 dy1ar[county] = fim[qp]-arg ; dy2ar[county] = fim[qp]-fim[qm] ; county++ ;
529 }
530 }
531
532 if( cc-1 >= 0 && cc+1 < nz ){
533 qp = kk+nxy ; qm = kk-nxy ;
534 if( INMASK(qp) && INMASK(qm) ){
535 dz1ar[countz] = fim[qp]-arg ; dz2ar[countz] = fim[qp]-fim[qm] ; countz++ ;
536 }
537 }
538 }
539
540 /*----- estimate variances of differences -----*/
541
542 qmedmad_float( countx , dx1ar , NULL , &vx1 ) ; vx1 = vx1*vx1 ;
543 qmedmad_float( county , dy1ar , NULL , &vy1 ) ; vy1 = vy1*vy1 ;
544 qmedmad_float( countz , dz1ar , NULL , &vz1 ) ; vz1 = vz1*vz1 ;
545 qmedmad_float( countx , dx2ar , NULL , &vx2 ) ; vx2 = vx2*vx2 ;
546 qmedmad_float( county , dy2ar , NULL , &vy2 ) ; vy2 = vy2*vy2 ;
547 qmedmad_float( countz , dz2ar , NULL , &vz2 ) ; vz2 = vz2*vz2 ;
548
549 if( ws == NULL ){
550 #pragma omp critical (MALLOC)
551 { free(dx1ar); free(dy1ar); free(dz1ar); free(dx2ar); free(dy2ar); free(dz2ar); }
552 }
553
554 /*----- now estimate the FWHMs -----*/
555
556 dx = im->dx; dy = im->dy; dz = im->dz;
557
558 /*--- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor ---*/
559 /*--- y = cbrt(12*sqrt(48-120*r+81*r*r)+108*r-80), and then
560 x = y/6 - 4/(3*y) - 1/3
561 is the real solution to the equation (1-x^4)/(1-x) = r > 1 ;
562 here, r = vx2/vx1 = ratio of variances at Delta=2*dx vs. Delta=1*dx;
563 x = exp[-dx**2/(4*sigma**2)] = correlation coefficient of neighbors;
564 we solve for x, then use that to solve for sigma, and scale to FWHM --*/
565
566 if( vx1 > 0.0 && vx2 > vx1 ){
567 arg = vx2 / vx1 ;
568 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
569 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
570 if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dx ;
571 }
572
573 if( vy1 > 0.0 && vy2 > vy1 ){
574 arg = vy2 / vy1 ;
575 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
576 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
577 if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dy ;
578 }
579
580 if( vz1 > 0.0 && vz2 > vz1 ){
581 arg = vz2 / vz1 ;
582 arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
583 arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
584 if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dz ;
585 }
586
587 LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
588 return fw_xyz ;
589 }
590
591 /*--------------------------------------------------------------------------*/
592
mri_nstat_fwhmbar12(int xx,int yy,int zz,MRI_IMAGE * im,byte * mask,MCW_cluster * nbhd,float * ws)593 float mri_nstat_fwhmbar12( int xx, int yy, int zz,
594 MRI_IMAGE *im, byte *mask, MCW_cluster *nbhd,float *ws )
595 {
596 THD_fvec3 fw ;
597 float fx,fy,fz , sum ; int nsum ;
598
599 fw = mri_nstat_fwhmxyz_12dif( xx,yy,zz , im,mask,nbhd,ws ) ;
600 UNLOAD_FVEC3(fw,fx,fy,fz) ;
601
602 sum = 0.0f ; nsum = 0 ;
603 if( fx > 0.0f ){ sum += fx ; nsum++ ; }
604 if( fy > 0.0f ){ sum += fy ; nsum++ ; }
605 if( fz > 0.0f ){ sum += fz ; nsum++ ; }
606 if( nsum > 0 ) sum /= nsum ;
607 return sum ;
608 }
609
610 /*--------------------------------------------------------------------------*/
611 #if 0
612 /*--------------------------------------------------------------------------*/
613 /*! Compute a local statistic at each voxel of an image, possibly with
614 a mask; 'local' is defined with a neighborhood; 'statistic' is defined
615 by an NSTAT_ code.
616 ----------------------------------------------------------------------------*/
617
618 MRI_IMAGE * mri_localstat( MRI_IMAGE *im, byte *mask, MCW_cluster *nbhd, int code )
619 {
620 MRI_IMAGE *outim , *nbim ;
621 float *outar ;
622 int ii,jj,kk , nx,ny,nz , ijk ;
623
624 ENTRY("mri_localstat") ;
625
626 if( im == NULL || nbhd == NULL ) RETURN(NULL) ;
627
628 outim = mri_new_conforming( im , MRI_float ) ;
629 outar = MRI_FLOAT_PTR(outim) ;
630 nx = outim->nx ; ny = outim->ny ; nz = outim->nz ;
631
632 for( ijk=kk=0 ; kk < nz ; kk++ ){
633 for( jj=0 ; jj < ny ; jj++ ){
634 for( ii=0 ; ii < nx ; ii++ ){
635 nbim = mri_get_nbhd( im , mask , ii,jj,kk , nbhd ) ;
636 outar[ijk++] = mri_nstat( code , nbim ) ;
637 mri_free(nbim) ;
638 }}}
639
640 RETURN(outim) ;
641 }
642 #endif
643
644 /*--------------------------------------------------------------------------*/
645
646 static int verb=0 , vn=0 ;
647 static int localstat_datum = MRI_float;
648
THD_localstat_datum(int i)649 void THD_localstat_datum(int i) {
650 localstat_datum=i;
651 if ( localstat_datum != MRI_byte &&
652 localstat_datum != MRI_short &&
653 localstat_datum != MRI_float) {
654 fprintf(stderr , "Warning: Datum can only be one of "
655 "MRI_byte, MRI_short or MRI_float\n"
656 "Setting datum to float default.\n");
657 localstat_datum = MRI_float;
658 }
659 }
660
THD_localstat_verb(int i)661 void THD_localstat_verb(int i){ verb=i; vn=0; }
662
663 #ifndef USE_OMP
vstep_print(void)664 static void vstep_print(void)
665 {
666 static char xx[10] = "0123456789" ;
667 fprintf(stderr , "%c" , xx[vn%10] ) ;
668 if( vn%10 == 9) fprintf(stderr,".") ;
669 vn++ ;
670 }
671 #endif
672
673 /*--------------------------------------------------------------------------*/
674 /*
675 Function to turn ijk 1D index on iset into its equivalent i,j,k on a
676 different gridset gset
677
678 ZSS Gov. Shutdown Imminent 2011
679 */
DSET_1Dindex_to_regrid_ijk(THD_3dim_dataset * iset,int ijk,THD_3dim_dataset * gset,int * ii,int * jj,int * kk)680 int DSET_1Dindex_to_regrid_ijk( THD_3dim_dataset *iset, int ijk,
681 THD_3dim_dataset *gset,
682 int *ii, int *jj, int *kk)
683 {
684 THD_fvec3 m_ncoord, m_ndicom;
685 THD_ivec3 m_nind3;
686
687 /* turn ijk on oset to RAI */
688 m_nind3.ijk[0] = DSET_index_to_ix(iset,ijk) ;
689 m_nind3.ijk[1] = DSET_index_to_jy(iset,ijk) ;
690 m_nind3.ijk[2] = DSET_index_to_kz(iset,ijk) ;
691 m_ncoord = THD_3dind_to_3dmm(iset,m_nind3);
692 /* setenv OMP_NUM_THREADS 1 when you uncomment debugging lines */
693 /*
694 if (ijk < 10) fprintf(stderr,"LR ijk %d [%d %d %d] [%f %f %f]\n",
695 ijk, m_nind3.ijk[0], m_nind3.ijk[1], m_nind3.ijk[2],
696 m_ncoord.xyz[0], m_ncoord.xyz[1], m_ncoord.xyz[2]);
697 */
698 m_ndicom = THD_3dmm_to_dicomm(iset,m_ncoord);
699 /* now go back to new grid ijk */
700 m_ncoord = THD_dicomm_to_3dmm(gset, m_ndicom);
701 m_nind3 = THD_3dmm_to_3dind(gset,m_ncoord);
702 *ii = m_nind3.ijk[0];
703 *jj = m_nind3.ijk[1];
704 *kk = m_nind3.ijk[2];
705 /*
706 if (ijk < 10) fprintf(stderr,"HR [%d %d %d] [%f %f %f]\n",
707 m_nind3.ijk[0], m_nind3.ijk[1], m_nind3.ijk[2],
708 m_ncoord.xyz[0], m_ncoord.xyz[1], m_ncoord.xyz[2]);
709 */
710 return(1);
711 }
712
713 /*
714 Return a downsampled grid, or just an empty copy
715 */
THD_reduced_grid_copy(THD_3dim_dataset * dset,float * redx)716 THD_3dim_dataset * THD_reduced_grid_copy(THD_3dim_dataset *dset, float *redx)
717 {
718 float dx_o,dy_o,dz_o ;
719 THD_3dim_dataset *oset=NULL;
720
721 if( dset == NULL) return(NULL);
722
723 if (!redx) {
724 oset = EDIT_empty_copy( dset ) ;
725 } else { /* create a new grid */
726 if (verb) {
727 INFO_message("Reducing output grid by %f %f %f",
728 redx[0], redx[1], redx[2]);
729 }
730 dx_o = fabs(DSET_DX(dset)*redx[0]);
731 dy_o = fabs(DSET_DY(dset)*redx[1]);
732 dz_o = fabs(DSET_DZ(dset)*redx[2]);
733 if (!(oset = r_new_resam_dset( dset, dset, dx_o, dy_o, dz_o,
734 NULL, 0, NULL, 0, 1))) {
735 ERROR_message("Failed to reduce output grid");
736 return(NULL);
737 }
738 }
739 return(oset);
740 }
741
742 /*--------------------------------------------------------------------------*/
743
THD_localstat(THD_3dim_dataset * dset,byte * mask,MCW_cluster * nbhd,int ncode,int * code,float codeparam[][MAX_CODE_PARAMS+1],float * redx,int resam_mode)744 THD_3dim_dataset * THD_localstat( THD_3dim_dataset *dset , byte *mask ,
745 MCW_cluster *nbhd , int ncode, int *code,
746 float codeparam[][MAX_CODE_PARAMS+1],
747 float *redx, int resam_mode)
748 {
749 THD_3dim_dataset *oset=NULL;
750 int iv,cc , nvin,nvout , nxyz_o , need_nbar=0,need_ws12=0 , npt ;
751 float **aar ;
752 MRI_IMAGE *dsim ;
753 float dx,dy,dz , fac ;
754 float *brick=NULL, voxval=0.0;
755 #ifndef USE_OMP
756 int vstep ;
757 #endif
758
759 ENTRY("THD_localstat") ;
760 /*fprintf(stderr,"%p, %p, %d, %p", dset, nbhd, ncode, code);*/
761 if( dset == NULL || nbhd == NULL || ncode < 1 || code == NULL ) RETURN(NULL);
762 /*fprintf(stderr,"%d\n", nbhd->num_pt);*/
763 npt = nbhd->num_pt ; if( npt == 0 ) RETURN(NULL);
764
765 /* check for stupid reduction parameters allowing = 1.0 for testing purposes*/
766 if (redx && redx[0] < 1.0 && redx[1]<1.0 && redx[2] <1.0) redx = NULL;
767
768 if (!(oset = THD_reduced_grid_copy(dset, redx))) {
769 ERROR_message("Failed to create output dset");
770 return(NULL);
771 }
772
773 nvin = DSET_NVALS( dset ) ;
774 nvout = nvin * ncode ;
775 EDIT_dset_items( oset ,
776 ADN_nvals , nvout ,
777 ADN_datum_all , MRI_float ,
778 ADN_nsl , 0 ,
779 ADN_brick_fac , NULL ,
780 ADN_prefix , "localstat" ,
781 ADN_none ) ;
782
783 nxyz_o = DSET_NX(oset)*DSET_NY(oset)*DSET_NZ(oset) ;
784
785 dx = fabsf(DSET_DX(dset)) ; if( dx <= 0.0f ) dx = 1.0f ;
786 dy = fabsf(DSET_DY(dset)) ; if( dy <= 0.0f ) dy = 1.0f ;
787 dz = fabsf(DSET_DZ(dset)) ; if( dz <= 0.0f ) dz = 1.0f ;
788
789 #ifndef USE_OMP
790 vstep = (verb && nxyz_o > 9999) ? nxyz_o/50 : 0 ;
791 #endif
792
793 aar = (float **)malloc(sizeof(float *)*ncode) ; /* output array of arrays */
794
795 for( cc=0 ; cc < ncode ; cc++ ){
796 if( code[cc] < NSTAT_FWHMx ){ need_nbar = 1; } /* need nbhd image */
797 if( code[cc] == NSTAT_FWHMbar12 ){ need_ws12 = 1; }
798 }
799
800 /** loop over input sub-bricks **/
801
802 for( iv=0 ; iv < nvin ; iv++ ){
803
804 #ifdef USE_OMP
805 if( verb && nxyz_o > 999 ) INFO_message("Start sub-brick [%d]",iv) ;
806 #endif
807
808 for( cc=0 ; cc < ncode ; cc++ ){ /* create output sub-bricks */
809 aar[cc] = (float *)malloc(sizeof(float)*nxyz_o) ;
810 if( aar[cc] == NULL )
811 ERROR_exit("THD_localstat: out of memory at iv=%d cc=%d",iv,cc);
812 }
813
814 /* extract copy of float-ized brick */
815
816 dsim = THD_extract_float_brick( iv , dset ) ;
817 dsim->dx = dx ; dsim->dy = dy ; dsim->dz = dz ;
818 brick = MRI_FLOAT_PTR(dsim);
819
820 /** loop over voxels **/
821
822 #ifndef USE_OMP
823 if( vstep ) fprintf(stderr,"++ voxel loop [%d]:",iv) ;
824 #endif
825
826 AFNI_OMP_START ;
827 #pragma omp parallel if( nxyz_o > 1111 ) /* parallelization: 13 Jul 2009 */
828 {
829 int ijk,kk,jj,ii,cc ;
830 THD_fvec3 fwv ;
831 double perc[MAX_CODE_PARAMS], mpv[MAX_CODE_PARAMS] ; /* no longer static */
832 float *nbar , fv6[6]={0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
833 int nbar_num=0, Hinit = 0;
834 float *ws12 ;
835
836 /* 17 Jul 2009: create workspace for neighborhood data */
837
838 nbar = (need_nbar) ? (float *)malloc(sizeof(float)*npt ) : NULL ;
839 ws12 = (need_ws12) ? (float *)malloc(sizeof(float)*npt*6) : NULL ;
840 #pragma omp for
841 for( ijk=0 ; ijk < nxyz_o ; ijk++ ){ /* parallelized loop */
842 if (!redx) { /* no grid change */
843 ii = DSET_index_to_ix(dset,ijk) ; /* convert ijk to voxel indexes */
844 jj = DSET_index_to_jy(dset,ijk) ;
845 kk = DSET_index_to_kz(dset,ijk) ;
846 } else {
847 /* get ii, jj, kk on original resolution */
848 DSET_1Dindex_to_regrid_ijk(oset, ijk, dset, &ii, &jj, &kk);
849 }
850
851 #ifndef USE_OMP
852 if( vstep && ijk%vstep==vstep-1 ) vstep_print() ;
853 #endif
854
855 if( need_nbar ) /* extract vector of data from voxel neighborhood */
856 nbar_num = mri_get_nbhd_array( dsim , mask,ii,jj,kk , nbhd , nbar ) ;
857
858 for( cc=0 ; cc < ncode ; cc++ ){ /* loop over desired statistics */
859
860 if( code[cc] == NSTAT_FWHMbar ){ /* 1 FWHM measurement */
861
862 aar[cc][ijk] = mri_nstat_fwhmbar( ii,jj,kk , dsim,mask,nbhd ) ;
863
864 } else if( code[cc] == NSTAT_FWHMbar12 ){
865
866 aar[cc][ijk] = mri_nstat_fwhmbar12( ii,jj,kk , dsim,mask,nbhd,ws12 ) ;
867
868 } else if( code[cc] == NSTAT_FWHMx ){ /* 3 FWHM measurements */
869
870 fwv = mri_nstat_fwhmxyz( ii,jj,kk , dsim,mask,nbhd ) ;
871 UNLOAD_FVEC3( fwv, aar[cc][ijk],aar[cc+1][ijk],aar[cc+2][ijk] ) ;
872 cc += 2 ; /* skip FWHMy and FWHMz codes that follow */
873
874 } else if( code[cc] == NSTAT_PERCENTILE ){ /* percentiles */
875
876 int N_mp, pp;
877 float *sfar=NULL;
878
879 if( codeparam[cc][0] < 1 )
880 ERROR_exit("THD_localstat: No percentile parameters set!");
881 N_mp = (int)codeparam[cc][0];
882 if( N_mp > MAX_CODE_PARAMS )
883 ERROR_exit( "THD_localstat: Cannot exceed %d params but have %d!" ,
884 MAX_CODE_PARAMS, N_mp);
885
886 for( pp=0 ; pp < N_mp ; ++pp )
887 mpv[pp] = (double)codeparam[cc][1+pp]/100.0;
888
889 if( nbar != NULL && nbar_num > 0 ){
890
891 if( !(sfar = (float *)Percentate( nbar ,
892 NULL , nbar_num ,
893 MRI_float , mpv , N_mp ,
894 0 , perc , 1, 1, 1 ) ) ) {
895
896 ERROR_exit("Failed to compute percentiles.");
897 }
898
899 for (pp=0; pp<N_mp; ++pp) aar[cc+pp][ijk] = (float)perc[pp];
900 } else {
901 for( pp=0 ; pp < N_mp; ++pp ) aar[cc+pp][ijk] = 0.0;
902 }
903
904 cc += (N_mp-1) ; /* number of sub-bricks added, minus 1 */
905
906 } else if( code[cc] == NSTAT_HIST ){ /* histograms */
907 static int N, iout;
908 static double W, min, max;
909 int pp, ib, cnt;
910
911 if (!Hinit) { /* init */
912 if( codeparam[cc][0] < 1 )
913 ERROR_exit("THD_localstat: No histogram parameters set!");
914 if( (int)codeparam[cc][0] != 4 )
915 ERROR_exit("THD_localstat: Expecting only 4 params, have %d!" ,
916 (int)codeparam[cc][0]);
917
918 min = codeparam[cc][1];
919 max = codeparam[cc][2];
920 N = (int)codeparam[cc][3];
921 iout= (int)codeparam[cc][4];
922 W = (max - min)/(double)N;
923 Hinit = 1;
924 }
925 for (ib=0; ib<N;++ib) aar[cc+ib][ijk] = 0;
926 if (iout) { /* ignore outliers */
927 for (pp=0, cnt=0; pp<nbar_num;++pp) {
928 if (nbar[pp]>=min && nbar[pp]<=max) {
929 ib = (int)((nbar[pp]-min)/W); if (ib==N) ib = N-1;
930 ++aar[cc+ib][ijk];
931 ++cnt;
932 }
933 }
934 } else {
935 cnt = nbar_num;
936 for (pp=0; pp<nbar_num;++pp) {
937 ib = (int)((nbar[pp]-min)/W);
938 if (ib>=N) ib = N-1;
939 else if (ib < 0) ib = 0;
940 ++aar[cc+ib][ijk];
941 }
942 }
943 for (ib=0; ib<N;++ib) aar[cc+ib][ijk] /= (float)cnt;
944
945 cc += (N-1) ; /* number of sub-bricks added, minus 1 */
946
947 } else if( code[cc] == NSTAT_LIST ){ /* Just all the neighbors mam*/
948 int pp;
949
950 if (nbar) {
951 for (pp=0; pp<nbar_num; ++pp) aar[cc+pp][ijk] = (float)nbar[pp];
952 }
953 cc += (nbhd->num_pt-1) ; /* number of sub-bricks added, minus 1
954 Do not use nbar_num-1 because at
955 volume edges you will not have
956 as many neighbors as elsewhere*/
957
958 } else if( code[cc] == NSTAT_diffs0 ){ /*3 values */
959 mri_nstat_diffs( nbar_num , nbar, fv6, 0 ) ;
960 aar[cc ][ijk] = fv6[0]; /* average difference */
961 aar[cc+1][ijk] = fv6[1]; /* minimum difference */
962 aar[cc+2][ijk] = fv6[2]; /* maximum difference */
963 cc += 2 ; /* skip redundant codes that follow */
964 } else if( code[cc] == NSTAT_adiffs0 ){ /*3 values */
965 mri_nstat_diffs( nbar_num , nbar, fv6, 1) ;
966 aar[cc ][ijk] = fv6[0]; /* average absolute difference */
967 aar[cc+1][ijk] = fv6[1]; /* minimum absolute difference */
968 aar[cc+2][ijk] = fv6[2]; /* maximum absolute difference */
969 cc += 2 ; /* skip redundant codes that follow */
970 } else if( code[cc] == NSTAT_mMP2s0 ){ /*3 values, median, MAD, P2Skew*/
971 mri_nstat_mMP2S( nbar_num , nbar, brick[ijk], fv6 ) ;
972 aar[cc ][ijk] = fv6[1]; /* median */
973 aar[cc+1][ijk] = fv6[3]; /* MAD */
974 aar[cc+2][ijk] = fv6[4]; /* Skew */
975 cc += 2 ; /* skip redundant codes that follow */
976 } else if( code[cc] == NSTAT_mmMP2s0 ){
977 /*4 values, mean, median, MAD, P2Skew*/
978 mri_nstat_mMP2S( nbar_num , nbar, brick[ijk], fv6 ) ;
979 aar[cc ][ijk] = fv6[0]; /* mean */
980 aar[cc+1][ijk] = fv6[1]; /* median */
981 aar[cc+2][ijk] = fv6[3]; /* MAD */
982 aar[cc+3][ijk] = fv6[4]; /* Skew */
983 cc += 3 ; /* skip redundant codes that follow */
984 } else { /* the "usual" (catchall) case */
985
986 aar[cc][ijk] = mri_nstat( code[cc] , nbar_num , nbar, brick[ijk], nbhd ) ;
987
988 }
989
990 } /* end of loop over cc */
991
992 } /** end of voxel loop **/
993
994 if( nbar != NULL ) free(nbar) ;
995 if( ws12 != NULL ) free(ws12) ;
996
997 } /* end OpenMP */
998 AFNI_OMP_END ;
999
1000 #ifndef USE_OMP
1001 if( vstep ) fprintf(stderr,"\n") ;
1002 #endif
1003
1004 if( dsim != NULL ){ mri_free(dsim); dsim = NULL; }
1005
1006 /* put data arrays from aar into the dataset */
1007
1008 for( cc=0 ; cc < ncode ; cc++ ) {
1009 /* EDIT_substitute_brick( oset , iv*ncode+cc , MRI_float , aar[cc] ) ; */
1010 EDIT_substscale_brick( oset , iv*ncode+cc , MRI_float , aar[cc],
1011 localstat_datum, -1.0);
1012 if( localstat_datum != MRI_float ) free(aar[cc]) ; /* 13 Jul 2009 */
1013 }
1014
1015 } /** end of sub-brick loop **/
1016
1017 free((void *)aar) ;
1018
1019 if ( resam_mode >= FIRST_RESAM_TYPE) {
1020 THD_3dim_dataset *tout=NULL;
1021 INFO_message("Restoring grid with %d resampling mode", resam_mode);
1022 /* resample back to original grid */
1023 if (!(tout = r_new_resam_dset( oset, dset, 0.0, 0.0, 0.0,
1024 NULL, resam_mode, NULL, 1, 1))) {
1025 ERROR_exit("Failed to reduce output grid");
1026 }
1027 DSET_delete(oset) ; oset = tout; tout = NULL;
1028 /* apply mask to resampled output */
1029 if (mask) THD_applydsetmask(oset,mask);
1030 }
1031 RETURN(oset) ;
1032 }
1033
set_mri_nstat_fillvalue(float tf)1034 void set_mri_nstat_fillvalue(float tf)
1035 {
1036 fillvalue = tf;
1037 }
1038
set_mri_nstat_unfillvalue(float tf)1039 void set_mri_nstat_unfillvalue(float tf)
1040 {
1041 unfillvalue = tf;
1042 }
1043
set_mri_nstat_maskvalue(float tf)1044 void set_mri_nstat_maskvalue(float tf)
1045 {
1046 maskvalue = tf;
1047 }
1048
set_mri_nstat_maskvalue2(float tf)1049 void set_mri_nstat_maskvalue2(float tf)
1050 {
1051 maskvalue2 = tf;
1052 }
1053