do { if ( val >= 0.0 )                                               \
           fprintf (fout, "sigma%c = %5.2f   FWHM%c = %5.2f \n",             \
                    c, val, c, val * 2.0*sqrt(2.0*log(2.0)));                \
         else                                                                \
           fprintf (fout, "sigma%c = NO_VALUE   FWHM%c = NO_VALUE\n",c,c);   \
    } while (0)

do{ int pv ; (ds) = THD_open_dataset((name)) ;                                    \
       if( !ISVALID_3DIM_DATASET((ds)) ){                                            \
          fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; }        \
       if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || (ds)->daxes->nzz!=nz ){   \
          fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; }             \
       if( DSET_NUM_TIMES((ds)) > 1 ){                                               \
         fprintf(stderr,"*** Can't use time-dependent data: %s\n",(name));exit(1); } \
       THD_load_datablock( (ds)->dblk ) ;                                            \
       pv = DSET_PRINCIPAL_VALUE((ds)) ;                                             \
       if( DSET_ARRAY((ds),pv) == NULL ){                                            \
          fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }          \
       if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                                \
          fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1); }     \
       break ; } while (0)

do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
         default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
            case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
                            (ptr) = (void *)( fim + (ind) ) ;                 \
            } break ;                                                         \
            case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
                            (ptr) = (void *)( fim + (ind) ) ;                 \
            } break ;                                                         \
            case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
                             (ptr) = (void *)( fim + (ind) ) ;                \
            } break ; } break ; } while(0)

00308 {
00309 }

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 }

00368 {
00369   int nx;                       /* number of voxels along x-axis */
00370   int ny;                       /* number of voxels along y-axis */
00371   int nz;                       /* number of voxels along z-axis */
00372   int nxy, nxyz;                /* total number of voxels */
00373   int ixyz;                     /* voxel index */
00374   float dx;                     /* voxel size along x-axis */
00375   float dy;                     /* voxel size along y-axis */
00376   float dz;                     /* voxel size along z-axis */
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   /*----- initialize local variables -----*/
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   /*----- estimate the variance of the data -----*/
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   /*----- estimate the partial derivatives -----*/
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   /*----- estimate the variance of the partial derivatives -----*/
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   /*----- now estimate the equivalent Gaussian filter width -----*/
00470   if ( var == 0.0 )
00471   {
00472     *sx = *sy = *sz = 0.0;    /* do not compute         08 Mar 2004 [rickr] */
00473   }
00474   else
00475   {
00476     arg = 1.0 - 0.5*(varxx/var);
00477     if ( (arg <= 0.0) || (varxx == 0.0) )
00478       *sx = -1.0;              /* flag value */
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;              /* flag value */
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;              /* flag value */
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       /* check flags instead...                      08 March 2004  [rickr]
00502        * printf ("   sx=%f    sy=%f    sz=%f \n", *sx, *sy, *sz);
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 }

00134 {
00135    fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00136    exit(1);
00137 }

00146 {
00147   
00148    THD_3dim_dataset * dset=NULL;
00149 
00150    /*----- read first dataset to get dimensions, etc. -----*/
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    /*----- voxel dimensions and data set dimensions -----*/
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 }

00227 {
00228   int nopt = 1;                  /* input option argument counter */
00229   int ival;                      /* integer input */
00230   float fval;                    /* float input */
00231   char message[MAX_NAME_LENGTH];            /* error message */
00232 
00233   
00234   /*----- does user request help menu? -----*/
00235   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00236   
00237   
00238   /*----- add to program log -----*/
00239   AFNI_logger (PROGRAM_NAME,argc,argv); 
00240 
00241 
00242   /*----- initialize the input options -----*/
00243   initialize_options (option_data);
00244     
00245 
00246   /*----- main loop over input options -----*/
00247   while (nopt < argc )
00248     {
00249       
00250       /*-----   -dset filename   -----*/
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       /*-----   -mask filename   -----*/
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       /*-----   -quiet q  -----*/
00275       if (strncmp(argv[nopt], "-quiet", 6) == 0)
00276         {
00277           option_data->quiet = 1;
00278           nopt++;
00279           continue;
00280         }
00281 
00282 
00283       /*-----   -out filename   -----*/
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       /*----- unknown command -----*/
00296       FWHM_error ("unrecognized command line option ");
00297     }
00298   
00299 }

00320 {
00321 
00322 
00323   /*----- allocate memory space for input options -----*/   
00324   *option_data = (input_options *) malloc(sizeof(input_options));
00325   if (*option_data == NULL)
00326     FWHM_error ("memory allocation error");
00327   
00328   /*----- get command line inputs -----*/
00329   get_options(argc, argv, *option_data);
00330 
00331   /*----- check for valid inputs -----*/
00332   check_for_valid_inputs (*option_data);
00333 
00334   /*-----  get data set dimensions -----*/
00335   get_dimensions (*option_data);
00336 
00337   /*----- allocate memory space for image data -----*/   
00338   *fim = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
00339   if (*fim == NULL)
00340     FWHM_error ("memory allocation error");
00341   
00342   /*----- read input data set -----*/
00343   read_afni_data (*option_data, (*option_data)->infilename, *fim);
00344 
00345 
00346   /*----- check for mask file -----*/
00347   if ((*option_data)->maskfilename != NULL)
00348     {
00349       /*----- allocate memory space for mask data -----*/   
00350       *fmask = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
00351       if (*fmask == NULL)  FWHM_error ("memory allocation error");
00352       
00353       /*----- read mask data set -----*/
00354       read_afni_data (*option_data, (*option_data)->maskfilename, *fmask);
00355       
00356     }
00357 
00358 }

00213 {
00214   option_data->infilename = NULL;    /* name of input file */
00215   option_data->maskfilename = NULL;  /* name of mask file */
00216   option_data->quiet = 0;            /* generate screen output (default)  */
00217   option_data->outfilename = NULL;   /* name of output file */
00218 }

00609 {
00610   input_options * option_data = NULL;
00611   float * fim = NULL;
00612   float * fmask = NULL;
00613   float sx, sy, sz;
00614 
00615   
00616   /*----- Identify software -----*/
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   /*----- program initialization -----*/
00626   initialize (argc, argv, &option_data, &fim, &fmask);
00627 
00628 
00629   /*----- estimate equivalent gaussian filter width -----*/
00630   estimate_gfw (option_data, fim, fmask, &sx, &sy, &sz );
00631   
00632   
00633   /*----- generate requested output -----*/
00634   output_results (option_data, sx, sy, sz);
00635   
00636   /*----- terminate program -----*/
00637   terminate (&option_data, &fim, &fmask);
00638   
00639   exit(0);
00640 }

00531 {
00532   char message[MAX_NAME_LENGTH];     /* error message */
00533   char filename[MAX_NAME_LENGTH];    /* output file name */ 
00534   FILE * fout;
00535 
00536 
00537   /*----- if output file has not been specified, use stdout -----*/
00538   if (option_data->outfilename == NULL)
00539     fout = stdout;
00540   else
00541     {
00542       /*----- see if output file already exists -----*/
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       /*----- open file for output -----*/
00552       fout = fopen (filename, "w");
00553       if (fout == NULL)
00554         { 
00555           FWHM_error ("unable to write file ");
00556         }
00557     }
00558   
00559   /*----- print out the results -----*/
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   /* check for flag value of -1.0                      08 Mar 2004 [rickr] */
00573   DISP_SIG_RESULTS('x', sx);
00574   DISP_SIG_RESULTS('y', sy);
00575   DISP_SIG_RESULTS('z', sz);
00576   
00577   /* report any failures */
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 }

00181 {
00182   int iv;                          /* index number of intensity sub-brick */
00183   THD_3dim_dataset * dset=NULL;    /* data set pointer */
00184   void * vfim = NULL;              /* image data pointer */
00185   int nx, ny, nz, nxyz;            /* data set dimensions in voxels */
00186   
00187   nx = option_data->nx;
00188   ny = option_data->ny;
00189   nz = option_data->nz;
00190   nxyz = option_data->nxyz;
00191   
00192   
00193   /*----- read in the data -----*/
00194   DOPEN(dset,filename) ;
00195   iv = DSET_PRINCIPAL_VALUE(dset) ;
00196   
00197   /*----- convert it to floats (in ffim) -----*/
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 ,      /* input  */
00201                           MRI_float               ,ffim  ) ;   /* output */
00202   
00203   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00204 }

00595 {
00596   free (*option_data);   *option_data = NULL;
00597   free (*fim);           *fim = NULL;
00598   if (*fmask != NULL)
00599     {  free (*fmask);       *fmask = NULL; }
00600 }