00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #define PROGRAM_NAME "3dFWHM"
00035 #define PROGRAM_AUTHOR "B. Douglas Ward"
00036 #define PROGRAM_INITIAL "20 February 1997"
00037 #define PROGRAM_LATEST "08 March 2004"
00038
00039
00040
00041 #include "mrilib.h"
00042
00043 #define MAX_NAME_LENGTH THD_MAX_NAME
00044
00045
00046
00047
00048
00049 #define DOPEN(ds,name) \
00050 do{ int pv ; (ds) = THD_open_dataset((name)) ; \
00051 if( !ISVALID_3DIM_DATASET((ds)) ){ \
00052 fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00053 if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || (ds)->daxes->nzz!=nz ){ \
00054 fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; } \
00055 if( DSET_NUM_TIMES((ds)) > 1 ){ \
00056 fprintf(stderr,"*** Can't use time-dependent data: %s\n",(name));exit(1); } \
00057 THD_load_datablock( (ds)->dblk ) ; \
00058 pv = DSET_PRINCIPAL_VALUE((ds)) ; \
00059 if( DSET_ARRAY((ds),pv) == NULL ){ \
00060 fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \
00061 if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \
00062 fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1); } \
00063 break ; } while (0)
00064
00065
00066
00067
00068
00069
00070 #define SUB_POINTER(ds,vv,ind,ptr) \
00071 do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \
00072 default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \
00073 case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \
00074 (ptr) = (void *)( fim + (ind) ) ; \
00075 } break ; \
00076 case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \
00077 (ptr) = (void *)( fim + (ind) ) ; \
00078 } break ; \
00079 case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \
00080 (ptr) = (void *)( fim + (ind) ) ; \
00081 } break ; } break ; } while(0)
00082
00083
00084
00085
00086
00087
00088
00089 typedef struct input_options
00090 {
00091 char * infilename;
00092 char * maskfilename;
00093 int nx;
00094 int ny;
00095 int nz;
00096 int nxyz;
00097 float dx;
00098 float dy;
00099 float dz;
00100 int quiet;
00101 char * outfilename;
00102 } input_options;
00103
00104
00105
00106
00107
00108
00109
00110 void display_help_menu()
00111 {
00112 printf
00113 (
00114 "This program estimates the Filter Width Half Maximum (FWHM). \n\n"
00115 "Usage: \n"
00116 "3dFWHM \n"
00117 "-dset file file = name of input AFNI 3d dataset \n"
00118 "[-mask mname] mname = filename of 3d mask dataset \n"
00119 "[-quiet] suppress screen output \n"
00120 "[-out file] file = name of output file \n"
00121 );
00122
00123 printf("\n" MASTER_SHORTHELP_STRING ) ;
00124
00125 exit(0);
00126 }
00127
00128
00129
00130
00131
00132
00133 void FWHM_error (char * message)
00134 {
00135 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00136 exit(1);
00137 }
00138
00139
00140
00141
00142
00143
00144
00145 void get_dimensions (input_options * option_data)
00146 {
00147
00148 THD_3dim_dataset * dset=NULL;
00149
00150
00151
00152 dset = THD_open_dataset( option_data->infilename) ;
00153 if( ! ISVALID_3DIM_DATASET(dset) ){
00154 fprintf(stderr,"*** Unable to open dataset file %s\n",
00155 option_data->infilename);
00156 exit(1) ;
00157 }
00158
00159
00160 option_data->dx = fabs(dset->daxes->xxdel) ;
00161 option_data->dy = fabs(dset->daxes->yydel) ;
00162 option_data->dz = fabs(dset->daxes->zzdel) ;
00163 option_data->nx = dset->daxes->nxx ;
00164 option_data->ny = dset->daxes->nyy ;
00165 option_data->nz = dset->daxes->nzz ;
00166 option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;
00167
00168 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00169
00170 }
00171
00172
00173
00174
00175
00176
00177
00178
00179 void read_afni_data (input_options * option_data, char * filename,
00180 float * ffim)
00181 {
00182 int iv;
00183 THD_3dim_dataset * dset=NULL;
00184 void * vfim = NULL;
00185 int nx, ny, nz, nxyz;
00186
00187 nx = option_data->nx;
00188 ny = option_data->ny;
00189 nz = option_data->nz;
00190 nxyz = option_data->nxyz;
00191
00192
00193
00194 DOPEN(dset,filename) ;
00195 iv = DSET_PRINCIPAL_VALUE(dset) ;
00196
00197
00198 SUB_POINTER(dset,iv,0,vfim) ;
00199 EDIT_coerce_scale_type( nxyz , DSET_BRICK_FACTOR(dset,iv) ,
00200 DSET_BRICK_TYPE(dset,iv),vfim ,
00201 MRI_float ,ffim ) ;
00202
00203 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00204 }
00205
00206
00207
00208
00209
00210
00211
00212 void initialize_options (input_options * option_data)
00213 {
00214 option_data->infilename = NULL;
00215 option_data->maskfilename = NULL;
00216 option_data->quiet = 0;
00217 option_data->outfilename = NULL;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226 void get_options (int argc, char ** argv, input_options * option_data)
00227 {
00228 int nopt = 1;
00229 int ival;
00230 float fval;
00231 char message[MAX_NAME_LENGTH];
00232
00233
00234
00235 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00236
00237
00238
00239 AFNI_logger (PROGRAM_NAME,argc,argv);
00240
00241
00242
00243 initialize_options (option_data);
00244
00245
00246
00247 while (nopt < argc )
00248 {
00249
00250
00251 if (strncmp(argv[nopt], "-dset", 5) == 0)
00252 {
00253 nopt++;
00254 if (nopt >= argc) FWHM_error ("need argument after -dset ");
00255 option_data->infilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00256 strcpy (option_data->infilename, argv[nopt]);
00257 nopt++;
00258 continue;
00259 }
00260
00261
00262
00263 if (strncmp(argv[nopt], "-mask", 5) == 0)
00264 {
00265 nopt++;
00266 if (nopt >= argc) FWHM_error ("need argument after -mask ");
00267 option_data->maskfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00268 strcpy (option_data->maskfilename, argv[nopt]);
00269 nopt++;
00270 continue;
00271 }
00272
00273
00274
00275 if (strncmp(argv[nopt], "-quiet", 6) == 0)
00276 {
00277 option_data->quiet = 1;
00278 nopt++;
00279 continue;
00280 }
00281
00282
00283
00284 if (strncmp(argv[nopt], "-out", 4) == 0)
00285 {
00286 nopt++;
00287 if (nopt >= argc) FWHM_error ("need argument after -out ");
00288 option_data->outfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
00289 strcpy (option_data->outfilename, argv[nopt]);
00290 nopt++;
00291 continue;
00292 }
00293
00294
00295
00296 FWHM_error ("unrecognized command line option ");
00297 }
00298
00299 }
00300
00301
00302
00303
00304
00305
00306
00307 void check_for_valid_inputs (input_options *option_data)
00308 {
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318 void initialize (int argc, char ** argv,
00319 input_options ** option_data, float ** fim, float ** fmask)
00320 {
00321
00322
00323
00324 *option_data = (input_options *) malloc(sizeof(input_options));
00325 if (*option_data == NULL)
00326 FWHM_error ("memory allocation error");
00327
00328
00329 get_options(argc, argv, *option_data);
00330
00331
00332 check_for_valid_inputs (*option_data);
00333
00334
00335 get_dimensions (*option_data);
00336
00337
00338 *fim = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
00339 if (*fim == NULL)
00340 FWHM_error ("memory allocation error");
00341
00342
00343 read_afni_data (*option_data, (*option_data)->infilename, *fim);
00344
00345
00346
00347 if ((*option_data)->maskfilename != NULL)
00348 {
00349
00350 *fmask = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
00351 if (*fmask == NULL) FWHM_error ("memory allocation error");
00352
00353
00354 read_afni_data (*option_data, (*option_data)->maskfilename, *fmask);
00355
00356 }
00357
00358 }
00359
00360
00361
00362
00363
00364
00365
00366 void estimate_gfw (input_options * option_data, float * fim, float * fmask,
00367 float * sx, float * sy, float * sz)
00368 {
00369 int nx;
00370 int ny;
00371 int nz;
00372 int nxy, nxyz;
00373 int ixyz;
00374 float dx;
00375 float dy;
00376 float dz;
00377 int ix, jy, kz, ixyz2;
00378 float fsum, fsq, var;
00379 float dfdx, dfdxsum, dfdxsq, varxx;
00380 float dfdy, dfdysum, dfdysq, varyy;
00381 float dfdz, dfdzsum, dfdzsq, varzz;
00382 int count, countx, county, countz;
00383 float arg;
00384
00385
00386
00387 nx = option_data->nx;
00388 ny = option_data->ny;
00389 nz = option_data->nz;
00390 dx = option_data->dx;
00391 dy = option_data->dy;
00392 dz = option_data->dz;
00393 nxyz = option_data->nxyz;
00394 nxy = nx * ny;
00395
00396
00397
00398 fsum = 0.0;
00399 fsq = 0.0;
00400 count = 0;
00401 for (ixyz = 0; ixyz < nxyz; ixyz++)
00402 {
00403 if (fmask != NULL)
00404 if (fmask[ixyz] == 0.0) continue;
00405
00406 count++;
00407 fsum += fim[ixyz];
00408 fsq += fim[ixyz] * fim[ixyz];
00409 }
00410 var = (fsq - (fsum * fsum)/count) / (count-1);
00411
00412
00413
00414 dfdxsum = 0.0; dfdysum = 0.0; dfdzsum = 0.0;
00415 dfdxsq = 0.0; dfdysq = 0.0; dfdzsq = 0.0;
00416 countx = 0; county = 0; countz = 0;
00417 for (ixyz = 0; ixyz < nxyz; ixyz++)
00418 {
00419 if (fmask != NULL)
00420 if (fmask[ixyz] == 0.0) continue;
00421
00422 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
00423
00424 if (ix+1 < nx)
00425 {
00426 ixyz2 = THREE_TO_IJK (ix+1, jy, kz, nx, nxy);
00427 dfdx = (fim[ixyz2] - fim[ixyz]) / 1.0;
00428 dfdxsum += dfdx;
00429 dfdxsq += dfdx * dfdx;
00430 countx += 1;
00431 }
00432
00433 if (jy+1 < ny)
00434 {
00435 ixyz2 = THREE_TO_IJK (ix, jy+1, kz, nx, nxy);
00436 dfdy = (fim[ixyz2] - fim[ixyz]) / 1.0;
00437 dfdysum += dfdy;
00438 dfdysq += dfdy * dfdy;
00439 county += 1;
00440 }
00441
00442 if (kz+1 < nz)
00443 {
00444 ixyz2 = THREE_TO_IJK (ix, jy, kz+1, nx, nxy);
00445 dfdz = (fim[ixyz2] - fim[ixyz]) / 1.0;
00446 dfdzsum += dfdz;
00447 dfdzsq += dfdz * dfdz;
00448 countz += 1;
00449 }
00450
00451 }
00452
00453
00454 if (countx < 2)
00455 varxx = 0.0;
00456 else
00457 varxx = (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1);
00458
00459 if (county < 2)
00460 varyy = 0.0;
00461 else
00462 varyy = (dfdysq - (dfdysum * dfdysum)/county) / (county-1);
00463
00464 if (countz < 2)
00465 varzz = 0.0;
00466 else
00467 varzz = (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1);
00468
00469
00470 if ( var == 0.0 )
00471 {
00472 *sx = *sy = *sz = 0.0;
00473 }
00474 else
00475 {
00476 arg = 1.0 - 0.5*(varxx/var);
00477 if ( (arg <= 0.0) || (varxx == 0.0) )
00478 *sx = -1.0;
00479 else
00480 *sx = sqrt( -1.0 / (4.0*log(arg)) ) * dx;
00481
00482 arg = 1.0 - 0.5*(varyy/var);
00483 if ( (arg <= 0.0) || (varyy == 0.0) )
00484 *sy = -1.0;
00485 else
00486 *sy = sqrt( -1.0 / (4.0*log(arg)) ) * dy;
00487
00488 arg = 1.0 - 0.5*(varzz/var);
00489 if ( (arg <= 0.0) || (varzz == 0.0) )
00490 *sz = -1.0;
00491 else
00492 *sz = sqrt( -1.0 / (4.0*log(arg)) ) * dz;
00493 }
00494
00495 if (!(option_data->quiet))
00496 {
00497 printf ("count=%d \n", count);
00498 printf ("var =%f \n", var);
00499 printf ("varxx=%f varyy=%f varzz=%f \n", varxx, varyy, varzz);
00500
00501
00502
00503
00504
00505 if ( *sx >= 0 ) printf(" sx=%f ", *sx);
00506 else printf(" sx=NO_VALUE ");
00507 if ( *sy >= 0 ) printf(" sy=%f ", *sy);
00508 else printf(" sy=NO_VALUE ");
00509 if ( *sz >= 0 ) printf(" sz=%f ", *sz);
00510 else printf(" sz=NO_VALUE ");
00511 putchar('\n');
00512 }
00513 }
00514
00515
00516 #define DISP_SIG_RESULTS(c,val) \
00517 do { if ( val >= 0.0 ) \
00518 fprintf (fout, "sigma%c = %5.2f FWHM%c = %5.2f \n", \
00519 c, val, c, val * 2.0*sqrt(2.0*log(2.0))); \
00520 else \
00521 fprintf (fout, "sigma%c = NO_VALUE FWHM%c = NO_VALUE\n",c,c); \
00522 } while (0)
00523
00524
00525
00526
00527
00528
00529
00530 void output_results (input_options * option_data, float sx, float sy, float sz)
00531 {
00532 char message[MAX_NAME_LENGTH];
00533 char filename[MAX_NAME_LENGTH];
00534 FILE * fout;
00535
00536
00537
00538 if (option_data->outfilename == NULL)
00539 fout = stdout;
00540 else
00541 {
00542
00543 strcpy (filename, option_data->outfilename);
00544 fout = fopen (filename, "r");
00545 if (fout != NULL)
00546 {
00547 sprintf (message, "file %s already exists. ", filename);
00548 FWHM_error (message);
00549 }
00550
00551
00552 fout = fopen (filename, "w");
00553 if (fout == NULL)
00554 {
00555 FWHM_error ("unable to write file ");
00556 }
00557 }
00558
00559
00560 fprintf (fout, "\n\n");
00561 fprintf (fout, "Gaussian filter widths: \n");
00562
00563 #if 0
00564 fprintf (fout, "sigmax = %5.2f FWHMx = %5.2f \n",
00565 sx, sx * 2.0*sqrt(2.0*log(2.0)));
00566 fprintf (fout, "sigmay = %5.2f FWHMy = %5.2f \n",
00567 sy, sy * 2.0*sqrt(2.0*log(2.0)));
00568 fprintf (fout, "sigmaz = %5.2f FWHMz = %5.2f \n\n",
00569 sz, sz * 2.0*sqrt(2.0*log(2.0)));
00570 #endif
00571
00572
00573 DISP_SIG_RESULTS('x', sx);
00574 DISP_SIG_RESULTS('y', sy);
00575 DISP_SIG_RESULTS('z', sz);
00576
00577
00578 if ( sx < 0.0 || sy < 0.0 || sz < 0.0 )
00579 fprintf(stderr,
00580 "\n"
00581 "** failure: some filter widths were not able to be estimated with\n"
00582 " this tool, such results were reported as 'NO_VALUE'\n");
00583
00584 fclose(fout);
00585
00586 }
00587
00588
00589
00590
00591
00592
00593
00594 void terminate (input_options ** option_data, float ** fim, float ** fmask)
00595 {
00596 free (*option_data); *option_data = NULL;
00597 free (*fim); *fim = NULL;
00598 if (*fmask != NULL)
00599 { free (*fmask); *fmask = NULL; }
00600 }
00601
00602
00603
00604
00605
00606
00607
00608 int main (int argc, char ** argv)
00609 {
00610 input_options * option_data = NULL;
00611 float * fim = NULL;
00612 float * fmask = NULL;
00613 float sx, sy, sz;
00614
00615
00616
00617 printf ("\n\n");
00618 printf ("Program: %s \n", PROGRAM_NAME);
00619 printf ("Author: %s \n", PROGRAM_AUTHOR);
00620 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00621 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00622 printf ("\n");
00623
00624
00625
00626 initialize (argc, argv, &option_data, &fim, &fmask);
00627
00628
00629
00630 estimate_gfw (option_data, fim, fmask, &sx, &sy, &sz );
00631
00632
00633
00634 output_results (option_data, sx, sy, sz);
00635
00636
00637 terminate (&option_data, &fim, &fmask);
00638
00639 exit(0);
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654