1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2001, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 /*
8   This program estimates the Filter Width Half Maximum (FWHM) of the Gaussian
9   filter required to produce clustering in noise data equivalent to that
10   of the input AFNI 3d dataset.
11 
12   File:    3dFWHM.c
13   Author:  B. D. Ward
14   Date:    20 February 1997
15 
16   Mod:     Added -mask option to restrict calculations to masked voxels only.
17            Also, allow the '[]' sub-brick selector for input datasets.
18   Date:    02 March 2000
19 
20   Mod:     Added call to AFNI_logger.
21   Date:    15 August 2001
22 
23   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
24   Date:    02 December 2002
25 
26   Mod:     Output NO_VALUE when results cannot be computed.
27   Date:    08 March 2004  [rickr]
28 */
29 
30 #undef ENABLE_THIS_PROGRAM      /* 03 Apr 2017 */
31 
32 #ifdef ENABLE_THIS_PROGRAM
33 
34 /*---------------------------------------------------------------------------*/
35 
36 #define PROGRAM_NAME "3dFWHM"                        /* name of this program */
37 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
38 #define PROGRAM_INITIAL "20 February 1997"/* date of initial program release */
39 #define PROGRAM_LATEST  "08 March 2004"   /* date of latest program revision */
40 
41 /*---------------------------------------------------------------------------*/
42 
43 #include "mrilib.h"
44 
45 #define MAX_NAME_LENGTH THD_MAX_NAME    /* max. string length for file names */
46 
47 static int dontcheckplus = 0 ;  /* 09 Nov 2006 */
48 
49 /*---------------------------------------------------------------------------*/
50 
51 /** macro to open a dataset and make it ready for processing **/
52 
53 #define DOPEN(ds,name)                                                               \
54    do{ int pv ; (ds) = THD_open_dataset((name)) ;                                    \
55        CHECK_OPEN_ERROR((ds),(name)) ;                                               \
56        if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || (ds)->daxes->nzz!=nz ){   \
57           fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; }             \
58        if( DSET_NUM_TIMES((ds)) > 1 ){                                               \
59          fprintf(stderr,"*** Can't use time-dependent data: %s\n",(name));exit(1); } \
60        DSET_load((ds)) ;                                                             \
61        pv = DSET_PRINCIPAL_VALUE((ds)) ;                                             \
62        if( DSET_ARRAY((ds),pv) == NULL ){                                            \
63           fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }          \
64        if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                                \
65           fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1); }     \
66        break ; } while (0)
67 
68 
69 /*---------------------------------------------------------------------------*/
70 
71 /** macro to return pointer to correct location in brick for current processing **/
72 
73 #define SUB_POINTER(ds,vv,ind,ptr)                                            \
74    do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
75          default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
76             case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
77                             (ptr) = (void *)( fim + (ind) ) ;                 \
78             } break ;                                                         \
79             case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
80                             (ptr) = (void *)( fim + (ind) ) ;                 \
81             } break ;                                                         \
82             case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
83                              (ptr) = (void *)( fim + (ind) ) ;                \
84             } break ; } break ; } while(0)
85 
86 
87 /*---------------------------------------------------------------------------*/
88 /*
89   Data structure.
90 */
91 
92 typedef struct input_options
93 {
94   char * infilename;            /* name of input file */
95   char * maskfilename;          /* name of mask file */
96   int nx;                       /* number of voxels along x-axis */
97   int ny;                       /* number of voxels along y-axis */
98   int nz;                       /* number of voxels along z-axis */
99   int nxyz;                     /* total number of voxels */
100   float dx;                     /* voxel size along x-axis */
101   float dy;                     /* voxel size along y-axis */
102   float dz;                     /* voxel size along z-axis */
103   int quiet;                    /* set to 1 to suppress screen output */
104   char * outfilename;           /* name of output file */
105 } input_options;
106 
107 
108 /*---------------------------------------------------------------------------*/
109 /*
110   Routine to display 3dFWHM help menu.
111 */
112 
display_help_menu()113 void display_help_menu()
114 {
115   printf
116     (
117      "This program estimates the Filter Width Half Maximum (FWHM).  \n\n"
118      "3dFWHM itself will no longer be upgraded.  Any future improvements\n"
119      "will be made to 3dFWHMx.  **** PLEASE SWITCH TO THAT PROGRAM ****\n\n"
120      "Usage: \n"
121      "3dFWHM \n"
122      "-dset file     file  = name of input AFNI 3d dataset \n"
123      "[-mask mname]  mname = filename of 3d mask dataset   \n"
124      "[-quiet]       suppress screen output                \n"
125      "[-out file]    file  = name of output file           \n"
126      "\n"
127      "[-compat] = Be compatible with the older 3dFWHM, where if a\n"
128      "            voxel is in the mask, then its neighbors are used\n"
129      "            for differencing, even if they are not themselves in\n"
130      "            the mask.  This was an error; now, neighbors must also\n"
131      "            be in the mask to be used in the differencing.\n"
132      "            Use '-compat' to use the older method [for comparison].\n"
133      "         ** This change made 09 Nov 2006.\n"
134      "\n"
135      "ALSO SEE:\n"
136      " - 3dFWHMx, which can deal with multi-brick datasets\n"
137      " - 3dLocalstat -stat FWHM, which can estimate the FWHM at each voxel\n"
138      "3dFWHM itself will no longer be upgraded.  Any future improvements\n"
139      "will be made to 3dFWHMx.  **** PLEASE SWITCH TO THAT PROGRAM ****\n"
140      "\n"
141      "INPUT FILE RECOMMENDATIONS:\n"
142      "For FMRI statistical purposes, you DO NOT want the FWHM to reflect\n"
143      "the spatial structure of the underlying anatomy.  Rather, you want\n"
144      "the FWHM to reflect the spatial structure of the noise.  This means\n"
145      "that the input dataset should not have anatomical structure.  One\n"
146      "good form of input is the output of '3dDeconvolve -errts', which is\n"
147      "the residuals left over after the GLM fitted signal model is subtracted\n"
148      "out from each voxel's time series.  If you don't want to go to that\n"
149      "trouble, use the output of 3dDetrend for the same purpose.  But just\n"
150      "giving a raw EPI dataset to this program will produce un-useful values.\n"
151      "\n"
152     );
153 
154    printf("\n" MASTER_SHORTHELP_STRING ) ;
155 
156   PRINT_COMPILE_DATE ; exit(0);
157 }
158 
159 /*---------------------------------------------------------------------------*/
160 /*
161    Routine to print error message and stop.
162 */
163 
FWHM_error(char * message)164 void FWHM_error (char * message)
165 {
166    fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
167    exit(1);
168 }
169 
170 
171 /*---------------------------------------------------------------------------*/
172 /*
173    Routine to get the dimensions of the 3d AFNI data sets.
174 */
175 
get_dimensions(input_options * option_data)176 void get_dimensions (input_options * option_data)
177 {
178 
179    THD_3dim_dataset * dset=NULL;
180 
181    /*----- read first dataset to get dimensions, etc. -----*/
182 
183    dset = THD_open_dataset( option_data->infilename) ;
184    CHECK_OPEN_ERROR(dset,option_data->infilename) ;
185 
186    /*----- voxel dimensions and data set dimensions -----*/
187    option_data->dx = fabs(dset->daxes->xxdel) ;
188    option_data->dy = fabs(dset->daxes->yydel) ;
189    option_data->dz = fabs(dset->daxes->zzdel) ;
190    option_data->nx = dset->daxes->nxx ;
191    option_data->ny = dset->daxes->nyy ;
192    option_data->nz = dset->daxes->nzz ;
193    option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;
194 
195    THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
196 
197 }
198 
199 
200 /*---------------------------------------------------------------------------*/
201 /*
202   Routine to read one AFNI data set from the input file.
203   The data is converted to floating point (in ffim).
204 */
205 
read_afni_data(input_options * option_data,char * filename,float * ffim)206 void read_afni_data (input_options * option_data,  char * filename,
207 		     float * ffim)
208 {
209   int iv;                          /* index number of intensity sub-brick */
210   THD_3dim_dataset * dset=NULL;    /* data set pointer */
211   void * vfim = NULL;              /* image data pointer */
212   int nx, ny, nz, nxyz;            /* data set dimensions in voxels */
213 
214   nx = option_data->nx;
215   ny = option_data->ny;
216   nz = option_data->nz;
217   nxyz = option_data->nxyz;
218 
219 
220   /*----- read in the data -----*/
221   DOPEN(dset,filename) ;
222   iv = DSET_PRINCIPAL_VALUE(dset) ;
223 
224   /*----- convert it to floats (in ffim) -----*/
225   SUB_POINTER(dset,iv,0,vfim) ;
226   EDIT_coerce_scale_type( nxyz , DSET_BRICK_FACTOR(dset,iv) ,
227 			  DSET_BRICK_TYPE(dset,iv),vfim ,      /* input  */
228 			  MRI_float               ,ffim  ) ;   /* output */
229 
230   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
231 }
232 
233 
234 /*---------------------------------------------------------------------------*/
235 /*
236   Routine to initialize the input options.
237 */
238 
initialize_options(input_options * option_data)239 void initialize_options (input_options * option_data)
240 {
241   option_data->infilename = NULL;    /* name of input file */
242   option_data->maskfilename = NULL;  /* name of mask file */
243   option_data->quiet = 0;            /* generate screen output (default)  */
244   option_data->outfilename = NULL;   /* name of output file */
245 }
246 
247 
248 /*---------------------------------------------------------------------------*/
249 /*
250   Routine to get user specified input options.
251 */
252 
get_options(int argc,char ** argv,input_options * option_data)253 void get_options (int argc, char ** argv, input_options * option_data)
254 {
255   int nopt = 1;                  /* input option argument counter */
256   int ival;                      /* integer input */
257   float fval;                    /* float input */
258   char message[MAX_NAME_LENGTH];            /* error message */
259 
260 
261   /*----- does user request help menu? -----*/
262   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
263 
264 
265   /*----- add to program log -----*/
266   AFNI_logger (PROGRAM_NAME,argc,argv);
267 
268 
269   /*----- initialize the input options -----*/
270   initialize_options (option_data);
271 
272 
273   /*----- main loop over input options -----*/
274   while (nopt < argc )
275     {
276 
277      if( strncmp(argv[nopt],"-compat",6) == 0 ){        /* 09 Nov 2006 */
278        dontcheckplus = 1 ; nopt++ ; continue ;
279      }
280 
281 
282       /*-----   -dset filename   -----*/
283       if (strncmp(argv[nopt], "-dset", 5) == 0 ||
284           strncmp(argv[nopt], "-input",5) == 0   )
285 	{
286 	  nopt++;
287 	  if (nopt >= argc)  FWHM_error ("need argument after -dset ");
288 	  option_data->infilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
289 	  strcpy (option_data->infilename, argv[nopt]);
290 	  nopt++;
291 	  continue;
292 	}
293 
294 
295       /*-----   -mask filename   -----*/
296       if (strncmp(argv[nopt], "-mask", 5) == 0)
297 	{
298 	  nopt++;
299 	  if (nopt >= argc)  FWHM_error ("need argument after -mask ");
300 	  option_data->maskfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
301 	  strcpy (option_data->maskfilename, argv[nopt]);
302 	  nopt++;
303 	  continue;
304 	}
305 
306 
307       /*-----   -quiet q  -----*/
308       if (strncmp(argv[nopt], "-quiet", 6) == 0)
309 	{
310 	  option_data->quiet = 1;
311 	  nopt++;
312 	  continue;
313 	}
314 
315 
316       /*-----   -out filename   -----*/
317       if (strncmp(argv[nopt], "-out", 4) == 0)
318 	{
319 	  nopt++;
320 	  if (nopt >= argc)  FWHM_error ("need argument after -out ");
321 	  option_data->outfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
322 	  strcpy (option_data->outfilename, argv[nopt]);
323 	  nopt++;
324 	  continue;
325 	}
326 
327 
328       /*----- unknown command -----*/
329       FWHM_error ("unrecognized command line option ");
330     }
331 
332 }
333 
334 
335 /*---------------------------------------------------------------------------*/
336 /*
337   Routine to check for valid inputs.
338 */
339 
check_for_valid_inputs(input_options * option_data)340 void check_for_valid_inputs (input_options *option_data)
341 {
342 }
343 
344 
345 /*---------------------------------------------------------------------------*/
346 /*
347   Routine to perform program initialization.
348 */
349 
350 
initialize(int argc,char ** argv,input_options ** option_data,float ** fim,float ** fmask)351 void initialize (int argc, char ** argv,
352 		 input_options ** option_data, float ** fim, float ** fmask)
353 {
354 
355 
356   /*----- allocate memory space for input options -----*/
357   *option_data = (input_options *) malloc(sizeof(input_options));
358   if (*option_data == NULL)
359     FWHM_error ("memory allocation error");
360 
361   /*----- get command line inputs -----*/
362   get_options(argc, argv, *option_data);
363 
364   /*----- check for valid inputs -----*/
365   check_for_valid_inputs (*option_data);
366 
367   /*-----  get data set dimensions -----*/
368   get_dimensions (*option_data);
369 
370   /*----- allocate memory space for image data -----*/
371   *fim = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
372   if (*fim == NULL)
373     FWHM_error ("memory allocation error");
374 
375   /*----- read input data set -----*/
376   read_afni_data (*option_data, (*option_data)->infilename, *fim);
377 
378 
379   /*----- check for mask file -----*/
380   if ((*option_data)->maskfilename != NULL)
381     {
382       /*----- allocate memory space for mask data -----*/
383       *fmask = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
384       if (*fmask == NULL)  FWHM_error ("memory allocation error");
385 
386       /*----- read mask data set -----*/
387       read_afni_data (*option_data, (*option_data)->maskfilename, *fmask);
388 
389     }
390 
391 }
392 
393 
394 /*---------------------------------------------------------------------------*/
395 /*
396   Routine to estimate the Gaussian filter width required to generate the data.
397 */
398 
estimate_gfw(input_options * option_data,float * fim,float * fmask,float * sx,float * sy,float * sz)399 void estimate_gfw (input_options * option_data, float * fim, float * fmask,
400 		   float * sx, float * sy, float * sz)
401 {
402   int nx;                       /* number of voxels along x-axis */
403   int ny;                       /* number of voxels along y-axis */
404   int nz;                       /* number of voxels along z-axis */
405   int nxy, nxyz;                /* total number of voxels */
406   int ixyz;                     /* voxel index */
407   float dx;                     /* voxel size along x-axis */
408   float dy;                     /* voxel size along y-axis */
409   float dz;                     /* voxel size along z-axis */
410   int ix, jy, kz, ixyz2;
411   float fsum, fsq, var;
412   float dfdx, dfdxsum, dfdxsq, varxx;
413   float dfdy, dfdysum, dfdysq, varyy;
414   float dfdz, dfdzsum, dfdzsq, varzz;
415   int count, countx, county, countz;
416   float arg;
417 
418 
419   /*----- initialize local variables -----*/
420   nx = option_data->nx;
421   ny = option_data->ny;
422   nz = option_data->nz;
423   dx = option_data->dx;
424   dy = option_data->dy;
425   dz = option_data->dz;
426   nxyz = option_data->nxyz;
427   nxy = nx * ny;
428 
429   if( fmask == NULL ) dontcheckplus = 1 ;
430 
431   /*----- estimate the variance of the data -----*/
432   fsum = 0.0;
433   fsq = 0.0;
434   count = 0;
435   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
436     {
437       if (fmask != NULL)
438 	if (fmask[ixyz] == 0.0)  continue;
439 
440       count++;
441       fsum += fim[ixyz];
442       fsq  += fim[ixyz] * fim[ixyz];
443     }
444   var = (fsq - (fsum * fsum)/count) / (count-1);
445 
446 
447   /*----- estimate the partial derivatives -----*/
448   dfdxsum = 0.0;   dfdysum = 0.0;   dfdzsum = 0.0;
449   dfdxsq = 0.0;    dfdysq  = 0.0;   dfdzsq = 0.0;
450   countx = 0;      county = 0;      countz = 0;
451   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
452     {
453       if (fmask != NULL)
454 	if (fmask[ixyz] == 0.0)  continue;
455 
456       IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
457 
458       if (ix+1 < nx)
459 	{
460 	  ixyz2 = THREE_TO_IJK (ix+1, jy, kz, nx, nxy);
461      if( dontcheckplus || fmask[ixyz2] != 0.0 ){
462 	    dfdx = (fim[ixyz2] - fim[ixyz]) / 1.0;
463 	    dfdxsum += dfdx;
464 	    dfdxsq  += dfdx * dfdx;
465 	    countx += 1;
466      }
467 	}
468 
469       if (jy+1 < ny)
470 	{
471 	  ixyz2 = THREE_TO_IJK (ix, jy+1, kz, nx, nxy);
472      if( dontcheckplus || fmask[ixyz2] != 0.0 ){
473 	    dfdy = (fim[ixyz2] - fim[ixyz]) / 1.0;
474 	    dfdysum += dfdy;
475 	    dfdysq  += dfdy * dfdy;
476 	    county += 1;
477      }
478 	}
479 
480       if (kz+1 < nz)
481 	{
482 	  ixyz2 = THREE_TO_IJK (ix, jy, kz+1, nx, nxy);
483      if( dontcheckplus || fmask[ixyz2] != 0.0 ){
484 	    dfdz = (fim[ixyz2] - fim[ixyz]) / 1.0;
485 	    dfdzsum += dfdz;
486 	    dfdzsq  += dfdz * dfdz;
487 	    countz += 1;
488      }
489 	}
490 
491      }
492 
493   /*----- estimate the variance of the partial derivatives -----*/
494   if (countx < 2)
495     varxx = 0.0;
496   else
497     varxx = (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1);
498 
499   if (county < 2)
500     varyy = 0.0;
501   else
502     varyy = (dfdysq - (dfdysum * dfdysum)/county) / (county-1);
503 
504   if (countz < 2)
505     varzz = 0.0;
506   else
507     varzz = (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1);
508 
509   /*----- now estimate the equivalent Gaussian filter width -----*/
510   if ( var == 0.0 )
511   {
512     *sx = *sy = *sz = 0.0;    /* do not compute         08 Mar 2004 [rickr] */
513   }
514   else
515   {
516     arg = 1.0 - 0.5*(varxx/var);
517     if ( (arg <= 0.0) || (varxx == 0.0) )
518       *sx = -1.0;              /* flag value */
519     else
520       *sx = sqrt( -1.0 / (4.0*log(arg)) ) * dx;
521 
522     arg = 1.0 - 0.5*(varyy/var);
523     if ( (arg <= 0.0) || (varyy == 0.0) )
524       *sy = -1.0;              /* flag value */
525     else
526       *sy = sqrt( -1.0 / (4.0*log(arg)) ) * dy;
527 
528     arg = 1.0 - 0.5*(varzz/var);
529     if ( (arg <= 0.0) || (varzz == 0.0) )
530       *sz = -1.0;              /* flag value */
531     else
532       *sz = sqrt( -1.0 / (4.0*log(arg)) ) * dz;
533   }
534 
535   if (!(option_data->quiet))
536     {
537       printf ("count=%d \n", count);
538       printf ("var  =%f \n", var);
539       printf ("varxx=%f varyy=%f varzz=%f \n", varxx, varyy, varzz);
540 
541       /* check flags instead...                      08 March 2004  [rickr]
542        * printf ("   sx=%f    sy=%f    sz=%f \n", *sx, *sy, *sz);
543        */
544 
545       if ( *sx >= 0 ) printf("   sx=%f ", *sx);
546       else            printf("   sx=NO_VALUE ");
547       if ( *sy >= 0 ) printf("   sy=%f ", *sy);
548       else            printf("   sy=NO_VALUE ");
549       if ( *sz >= 0 ) printf("   sz=%f ", *sz);
550       else            printf("   sz=NO_VALUE ");
551       putchar('\n');
552     }
553 }
554 
555 /* 						     08 March 2004  [rickr] */
556 #define DISP_SIG_RESULTS(c,val) 					     \
557     do { if ( val >= 0.0 ) 						     \
558 	   fprintf (fout, "sigma%c = %5.2f   FWHM%c = %5.2f \n", 	     \
559   	            c, val, c, val * 2.0*sqrt(2.0*log(2.0)));		     \
560 	 else								     \
561 	   fprintf (fout, "sigma%c = NO_VALUE   FWHM%c = NO_VALUE\n",c,c);   \
562     } while (0)
563 
564 
565 /*---------------------------------------------------------------------------*/
566 /*
567   Routine to generate requested output.
568 */
569 
output_results(input_options * option_data,float sx,float sy,float sz)570 void output_results (input_options * option_data, float sx, float sy, float sz)
571 {
572   char message[MAX_NAME_LENGTH];     /* error message */
573   char filename[MAX_NAME_LENGTH];    /* output file name */
574   FILE * fout;
575 
576 
577   /*----- if output file has not been specified, use stdout -----*/
578   if (option_data->outfilename == NULL)
579     fout = stdout;
580   else
581     {
582       /*----- see if output file already exists -----*/
583       strcpy (filename, option_data->outfilename);
584       fout = fopen (filename, "r");
585       if (fout != NULL)
586 	{
587 	  sprintf (message, "file %s already exists. ", filename);
588 	  FWHM_error (message);
589 	}
590 
591       /*----- open file for output -----*/
592       fout = fopen (filename, "w");
593       if (fout == NULL)
594 	{
595 	  FWHM_error ("unable to write file ");
596 	}
597     }
598 
599   /*----- print out the results -----*/
600   fprintf (fout, "\n\n");
601   fprintf (fout, "Gaussian filter widths: \n");
602 
603 #if 0
604   fprintf (fout, "sigmax = %5.2f   FWHMx = %5.2f \n",
605   	   sx, sx * 2.0*sqrt(2.0*log(2.0)));
606   fprintf (fout, "sigmay = %5.2f   FWHMy = %5.2f \n",
607 	   sy, sy * 2.0*sqrt(2.0*log(2.0)));
608   fprintf (fout, "sigmaz = %5.2f   FWHMz = %5.2f \n\n",
609 	   sz, sz * 2.0*sqrt(2.0*log(2.0)));
610 #endif
611 
612   /* check for flag value of -1.0                      08 Mar 2004 [rickr] */
613   DISP_SIG_RESULTS('x', sx);
614   DISP_SIG_RESULTS('y', sy);
615   DISP_SIG_RESULTS('z', sz);
616 
617   /* report any failures */
618   if ( sx < 0.0 || sy < 0.0 || sz < 0.0 )
619     fprintf(stderr,
620 	"\n"
621 	"** failure: some filter widths were not able to be estimated with\n"
622         "            this tool, such results were reported as 'NO_VALUE'\n");
623 
624   fclose(fout);
625 
626 }
627 
628 
629 /*---------------------------------------------------------------------------*/
630 /*
631   Routine to terminate program.
632 */
633 
terminate(input_options ** option_data,float ** fim,float ** fmask)634 void terminate (input_options ** option_data, float ** fim, float ** fmask)
635 {
636   free (*option_data);   *option_data = NULL;
637   free (*fim);           *fim = NULL;
638   if (*fmask != NULL)
639     {  free (*fmask);       *fmask = NULL; }
640 }
641 
642 #else   /* ENABLE_THIS_PROGRAM */
643 
644 #include <stdio.h>
645 #include <stdlib.h>
646 
647 #endif  /* ENABLE_THIS_PROGRAM */
648 
649 /*---------------------------------------------------------------------------*/
650 /*
651   Calculation of FWHM.
652 */
653 
main(int argc,char ** argv)654 int main (int argc, char ** argv)
655 {
656 
657 #ifdef ENABLE_THIS_PROGRAM
658   input_options * option_data = NULL;
659   float * fim = NULL;
660   float * fmask = NULL;
661   float sx, sy, sz;
662 #endif
663 
664 
665   /*----- Identify software -----*/
666 
667 #ifndef ENABLE_THIS_PROGRAM
668   fprintf(stderr,"\n\n") ;
669   fprintf(stderr,"*** Program 3dFWHM is obsolete!  Use 3dFWHMx instead! ***\n\n\n") ;
670   exit(0) ;
671 #else
672 
673   PRINT_VERSION("3dFWHM") ;
674   AUTHOR(PROGRAM_AUTHOR);
675 
676   mainENTRY("3dFWHM main") ;
677   machdep() ;
678 
679   /*----- program initialization -----*/
680   initialize (argc, argv, &option_data, &fim, &fmask);
681 
682   /*----- estimate equivalent gaussian filter width -----*/
683   estimate_gfw (option_data, fim, fmask, &sx, &sy, &sz );
684 
685   /*----- generate requested output -----*/
686   output_results (option_data, sx, sy, sz);
687 
688   /*----- terminate program -----*/
689   terminate (&option_data, &fim, &fmask);
690 
691 #endif /* ENABLE_THIS_PROGRAM */
692 
693   exit(0);
694 }
695