/*****************************************************************************
   Major portions of this software are copyrighted by the Medical College
   of Wisconsin, 1994-2001, and are released under the Gnu General Public
   License, Version 2.  See the file README.Copyright for details.
******************************************************************************/

/*
  This program estimates the Filter Width Half Maximum (FWHM) of the Gaussian
  filter required to produce clustering in noise data equivalent to that
  of the input AFNI 3d dataset.

  File:    3dFWHM.c
  Author:  B. D. Ward
  Date:    20 February 1997

  Mod:     Added -mask option to restrict calculations to masked voxels only.
           Also, allow the '[]' sub-brick selector for input datasets.
  Date:    02 March 2000

  Mod:     Added call to AFNI_logger.
  Date:    15 August 2001

  Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
  Date:    02 December 2002

  Mod:     Output NO_VALUE when results cannot be computed.
  Date:    08 March 2004  [rickr]
*/

#undef ENABLE_THIS_PROGRAM      /* 03 Apr 2017 */

#ifdef ENABLE_THIS_PROGRAM

/*---------------------------------------------------------------------------*/

#define PROGRAM_NAME "3dFWHM"                        /* name of this program */
#define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
#define PROGRAM_INITIAL "20 February 1997"/* date of initial program release */
#define PROGRAM_LATEST  "08 March 2004"   /* date of latest program revision */

/*---------------------------------------------------------------------------*/

#include "mrilib.h"

#define MAX_NAME_LENGTH THD_MAX_NAME    /* max. string length for file names */

static int dontcheckplus = 0 ;  /* 09 Nov 2006 */

/*---------------------------------------------------------------------------*/

/** macro to open a dataset and make it ready for processing **/

#define DOPEN(ds,name)                                                               \
   do{ int pv ; (ds) = THD_open_dataset((name)) ;                                    \
       CHECK_OPEN_ERROR((ds),(name)) ;                                               \
       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); } \
       DSET_load((ds)) ;                                                             \
       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)


/*---------------------------------------------------------------------------*/

/** macro to return pointer to correct location in brick for current processing **/

#define SUB_POINTER(ds,vv,ind,ptr)                                            \
   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)


/*---------------------------------------------------------------------------*/
/*
  Data structure.
*/

typedef struct input_options
{
  char * infilename;            /* name of input file */
  char * maskfilename;          /* name of mask file */
  int nx;                       /* number of voxels along x-axis */
  int ny;                       /* number of voxels along y-axis */
  int nz;                       /* number of voxels along z-axis */
  int nxyz;                     /* total number of voxels */
  float dx;                     /* voxel size along x-axis */
  float dy;                     /* voxel size along y-axis */
  float dz;                     /* voxel size along z-axis */
  int quiet;                    /* set to 1 to suppress screen output */
  char * outfilename;           /* name of output file */
} input_options;


/*---------------------------------------------------------------------------*/
/*
  Routine to display 3dFWHM help menu.
*/

void display_help_menu()
{
  printf
    (
     "This program estimates the Filter Width Half Maximum (FWHM).  \n\n"
     "3dFWHM itself will no longer be upgraded.  Any future improvements\n"
     "will be made to 3dFWHMx.  **** PLEASE SWITCH TO THAT PROGRAM ****\n\n"
     "Usage: \n"
     "3dFWHM \n"
     "-dset file     file  = name of input AFNI 3d dataset \n"
     "[-mask mname]  mname = filename of 3d mask dataset   \n"
     "[-quiet]       suppress screen output                \n"
     "[-out file]    file  = name of output file           \n"
     "\n"
     "[-compat] = Be compatible with the older 3dFWHM, where if a\n"
     "            voxel is in the mask, then its neighbors are used\n"
     "            for differencing, even if they are not themselves in\n"
     "            the mask.  This was an error; now, neighbors must also\n"
     "            be in the mask to be used in the differencing.\n"
     "            Use '-compat' to use the older method [for comparison].\n"
     "         ** This change made 09 Nov 2006.\n"
     "\n"
     "ALSO SEE:\n"
     " - 3dFWHMx, which can deal with multi-brick datasets\n"
     " - 3dLocalstat -stat FWHM, which can estimate the FWHM at each voxel\n"
     "3dFWHM itself will no longer be upgraded.  Any future improvements\n"
     "will be made to 3dFWHMx.  **** PLEASE SWITCH TO THAT PROGRAM ****\n"
     "\n"
     "INPUT FILE RECOMMENDATIONS:\n"
     "For FMRI statistical purposes, you DO NOT want the FWHM to reflect\n"
     "the spatial structure of the underlying anatomy.  Rather, you want\n"
     "the FWHM to reflect the spatial structure of the noise.  This means\n"
     "that the input dataset should not have anatomical structure.  One\n"
     "good form of input is the output of '3dDeconvolve -errts', which is\n"
     "the residuals left over after the GLM fitted signal model is subtracted\n"
     "out from each voxel's time series.  If you don't want to go to that\n"
     "trouble, use the output of 3dDetrend for the same purpose.  But just\n"
     "giving a raw EPI dataset to this program will produce un-useful values.\n"
     "\n"
    );

   printf("\n" MASTER_SHORTHELP_STRING ) ;

  PRINT_COMPILE_DATE ; exit(0);
}

/*---------------------------------------------------------------------------*/
/*
   Routine to print error message and stop.
*/

void FWHM_error (char * message)
{
   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
   exit(1);
}


/*---------------------------------------------------------------------------*/
/*
   Routine to get the dimensions of the 3d AFNI data sets.
*/

void get_dimensions (input_options * option_data)
{

   THD_3dim_dataset * dset=NULL;

   /*----- read first dataset to get dimensions, etc. -----*/

   dset = THD_open_dataset( option_data->infilename) ;
   CHECK_OPEN_ERROR(dset,option_data->infilename) ;

   /*----- voxel dimensions and data set dimensions -----*/
   option_data->dx = fabs(dset->daxes->xxdel) ;
   option_data->dy = fabs(dset->daxes->yydel) ;
   option_data->dz = fabs(dset->daxes->zzdel) ;
   option_data->nx = dset->daxes->nxx ;
   option_data->ny = dset->daxes->nyy ;
   option_data->nz = dset->daxes->nzz ;
   option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;

   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;

}


/*---------------------------------------------------------------------------*/
/*
  Routine to read one AFNI data set from the input file.
  The data is converted to floating point (in ffim).
*/

void read_afni_data (input_options * option_data,  char * filename,
		     float * ffim)
{
  int iv;                          /* index number of intensity sub-brick */
  THD_3dim_dataset * dset=NULL;    /* data set pointer */
  void * vfim = NULL;              /* image data pointer */
  int nx, ny, nz, nxyz;            /* data set dimensions in voxels */

  nx = option_data->nx;
  ny = option_data->ny;
  nz = option_data->nz;
  nxyz = option_data->nxyz;


  /*----- read in the data -----*/
  DOPEN(dset,filename) ;
  iv = DSET_PRINCIPAL_VALUE(dset) ;

  /*----- convert it to floats (in ffim) -----*/
  SUB_POINTER(dset,iv,0,vfim) ;
  EDIT_coerce_scale_type( nxyz , DSET_BRICK_FACTOR(dset,iv) ,
			  DSET_BRICK_TYPE(dset,iv),vfim ,      /* input  */
			  MRI_float               ,ffim  ) ;   /* output */

  THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
}


/*---------------------------------------------------------------------------*/
/*
  Routine to initialize the input options.
*/

void initialize_options (input_options * option_data)
{
  option_data->infilename = NULL;    /* name of input file */
  option_data->maskfilename = NULL;  /* name of mask file */
  option_data->quiet = 0;            /* generate screen output (default)  */
  option_data->outfilename = NULL;   /* name of output file */
}


/*---------------------------------------------------------------------------*/
/*
  Routine to get user specified input options.
*/

void get_options (int argc, char ** argv, input_options * option_data)
{
  int nopt = 1;                  /* input option argument counter */
  int ival;                      /* integer input */
  float fval;                    /* float input */
  char message[MAX_NAME_LENGTH];            /* error message */


  /*----- does user request help menu? -----*/
  if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();


  /*----- add to program log -----*/
  AFNI_logger (PROGRAM_NAME,argc,argv);


  /*----- initialize the input options -----*/
  initialize_options (option_data);


  /*----- main loop over input options -----*/
  while (nopt < argc )
    {

     if( strncmp(argv[nopt],"-compat",6) == 0 ){        /* 09 Nov 2006 */
       dontcheckplus = 1 ; nopt++ ; continue ;
     }


      /*-----   -dset filename   -----*/
      if (strncmp(argv[nopt], "-dset", 5) == 0 ||
          strncmp(argv[nopt], "-input",5) == 0   )
	{
	  nopt++;
	  if (nopt >= argc)  FWHM_error ("need argument after -dset ");
	  option_data->infilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
	  strcpy (option_data->infilename, argv[nopt]);
	  nopt++;
	  continue;
	}


      /*-----   -mask filename   -----*/
      if (strncmp(argv[nopt], "-mask", 5) == 0)
	{
	  nopt++;
	  if (nopt >= argc)  FWHM_error ("need argument after -mask ");
	  option_data->maskfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
	  strcpy (option_data->maskfilename, argv[nopt]);
	  nopt++;
	  continue;
	}


      /*-----   -quiet q  -----*/
      if (strncmp(argv[nopt], "-quiet", 6) == 0)
	{
	  option_data->quiet = 1;
	  nopt++;
	  continue;
	}


      /*-----   -out filename   -----*/
      if (strncmp(argv[nopt], "-out", 4) == 0)
	{
	  nopt++;
	  if (nopt >= argc)  FWHM_error ("need argument after -out ");
	  option_data->outfilename = malloc (sizeof(char) * MAX_NAME_LENGTH);
	  strcpy (option_data->outfilename, argv[nopt]);
	  nopt++;
	  continue;
	}


      /*----- unknown command -----*/
      FWHM_error ("unrecognized command line option ");
    }

}


/*---------------------------------------------------------------------------*/
/*
  Routine to check for valid inputs.
*/

void check_for_valid_inputs (input_options *option_data)
{
}


/*---------------------------------------------------------------------------*/
/*
  Routine to perform program initialization.
*/


void initialize (int argc, char ** argv,
		 input_options ** option_data, float ** fim, float ** fmask)
{


  /*----- allocate memory space for input options -----*/
  *option_data = (input_options *) malloc(sizeof(input_options));
  if (*option_data == NULL)
    FWHM_error ("memory allocation error");

  /*----- get command line inputs -----*/
  get_options(argc, argv, *option_data);

  /*----- check for valid inputs -----*/
  check_for_valid_inputs (*option_data);

  /*-----  get data set dimensions -----*/
  get_dimensions (*option_data);

  /*----- allocate memory space for image data -----*/
  *fim = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
  if (*fim == NULL)
    FWHM_error ("memory allocation error");

  /*----- read input data set -----*/
  read_afni_data (*option_data, (*option_data)->infilename, *fim);


  /*----- check for mask file -----*/
  if ((*option_data)->maskfilename != NULL)
    {
      /*----- allocate memory space for mask data -----*/
      *fmask = (float *) malloc( (*option_data)->nxyz * sizeof(float) );
      if (*fmask == NULL)  FWHM_error ("memory allocation error");

      /*----- read mask data set -----*/
      read_afni_data (*option_data, (*option_data)->maskfilename, *fmask);

    }

}


/*---------------------------------------------------------------------------*/
/*
  Routine to estimate the Gaussian filter width required to generate the data.
*/

void estimate_gfw (input_options * option_data, float * fim, float * fmask,
		   float * sx, float * sy, float * sz)
{
  int nx;                       /* number of voxels along x-axis */
  int ny;                       /* number of voxels along y-axis */
  int nz;                       /* number of voxels along z-axis */
  int nxy, nxyz;                /* total number of voxels */
  int ixyz;                     /* voxel index */
  float dx;                     /* voxel size along x-axis */
  float dy;                     /* voxel size along y-axis */
  float dz;                     /* voxel size along z-axis */
  int ix, jy, kz, ixyz2;
  float fsum, fsq, var;
  float dfdx, dfdxsum, dfdxsq, varxx;
  float dfdy, dfdysum, dfdysq, varyy;
  float dfdz, dfdzsum, dfdzsq, varzz;
  int count, countx, county, countz;
  float arg;


  /*----- initialize local variables -----*/
  nx = option_data->nx;
  ny = option_data->ny;
  nz = option_data->nz;
  dx = option_data->dx;
  dy = option_data->dy;
  dz = option_data->dz;
  nxyz = option_data->nxyz;
  nxy = nx * ny;

  if( fmask == NULL ) dontcheckplus = 1 ;

  /*----- estimate the variance of the data -----*/
  fsum = 0.0;
  fsq = 0.0;
  count = 0;
  for (ixyz = 0;  ixyz < nxyz;  ixyz++)
    {
      if (fmask != NULL)
	if (fmask[ixyz] == 0.0)  continue;

      count++;
      fsum += fim[ixyz];
      fsq  += fim[ixyz] * fim[ixyz];
    }
  var = (fsq - (fsum * fsum)/count) / (count-1);


  /*----- estimate the partial derivatives -----*/
  dfdxsum = 0.0;   dfdysum = 0.0;   dfdzsum = 0.0;
  dfdxsq = 0.0;    dfdysq  = 0.0;   dfdzsq = 0.0;
  countx = 0;      county = 0;      countz = 0;
  for (ixyz = 0;  ixyz < nxyz;  ixyz++)
    {
      if (fmask != NULL)
	if (fmask[ixyz] == 0.0)  continue;

      IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);

      if (ix+1 < nx)
	{
	  ixyz2 = THREE_TO_IJK (ix+1, jy, kz, nx, nxy);
     if( dontcheckplus || fmask[ixyz2] != 0.0 ){
	    dfdx = (fim[ixyz2] - fim[ixyz]) / 1.0;
	    dfdxsum += dfdx;
	    dfdxsq  += dfdx * dfdx;
	    countx += 1;
     }
	}

      if (jy+1 < ny)
	{
	  ixyz2 = THREE_TO_IJK (ix, jy+1, kz, nx, nxy);
     if( dontcheckplus || fmask[ixyz2] != 0.0 ){
	    dfdy = (fim[ixyz2] - fim[ixyz]) / 1.0;
	    dfdysum += dfdy;
	    dfdysq  += dfdy * dfdy;
	    county += 1;
     }
	}

      if (kz+1 < nz)
	{
	  ixyz2 = THREE_TO_IJK (ix, jy, kz+1, nx, nxy);
     if( dontcheckplus || fmask[ixyz2] != 0.0 ){
	    dfdz = (fim[ixyz2] - fim[ixyz]) / 1.0;
	    dfdzsum += dfdz;
	    dfdzsq  += dfdz * dfdz;
	    countz += 1;
     }
	}

     }

  /*----- estimate the variance of the partial derivatives -----*/
  if (countx < 2)
    varxx = 0.0;
  else
    varxx = (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1);

  if (county < 2)
    varyy = 0.0;
  else
    varyy = (dfdysq - (dfdysum * dfdysum)/county) / (county-1);

  if (countz < 2)
    varzz = 0.0;
  else
    varzz = (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1);

  /*----- now estimate the equivalent Gaussian filter width -----*/
  if ( var == 0.0 )
  {
    *sx = *sy = *sz = 0.0;    /* do not compute         08 Mar 2004 [rickr] */
  }
  else
  {
    arg = 1.0 - 0.5*(varxx/var);
    if ( (arg <= 0.0) || (varxx == 0.0) )
      *sx = -1.0;              /* flag value */
    else
      *sx = sqrt( -1.0 / (4.0*log(arg)) ) * dx;

    arg = 1.0 - 0.5*(varyy/var);
    if ( (arg <= 0.0) || (varyy == 0.0) )
      *sy = -1.0;              /* flag value */
    else
      *sy = sqrt( -1.0 / (4.0*log(arg)) ) * dy;

    arg = 1.0 - 0.5*(varzz/var);
    if ( (arg <= 0.0) || (varzz == 0.0) )
      *sz = -1.0;              /* flag value */
    else
      *sz = sqrt( -1.0 / (4.0*log(arg)) ) * dz;
  }

  if (!(option_data->quiet))
    {
      printf ("count=%d \n", count);
      printf ("var  =%f \n", var);
      printf ("varxx=%f varyy=%f varzz=%f \n", varxx, varyy, varzz);

      /* check flags instead...                      08 March 2004  [rickr]
       * printf ("   sx=%f    sy=%f    sz=%f \n", *sx, *sy, *sz);
       */

      if ( *sx >= 0 ) printf("   sx=%f ", *sx);
      else            printf("   sx=NO_VALUE ");
      if ( *sy >= 0 ) printf("   sy=%f ", *sy);
      else            printf("   sy=NO_VALUE ");
      if ( *sz >= 0 ) printf("   sz=%f ", *sz);
      else            printf("   sz=NO_VALUE ");
      putchar('\n');
    }
}

/* 						     08 March 2004  [rickr] */
#define DISP_SIG_RESULTS(c,val) 					     \
    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)


/*---------------------------------------------------------------------------*/
/*
  Routine to generate requested output.
*/

void output_results (input_options * option_data, float sx, float sy, float sz)
{
  char message[MAX_NAME_LENGTH];     /* error message */
  char filename[MAX_NAME_LENGTH];    /* output file name */
  FILE * fout;


  /*----- if output file has not been specified, use stdout -----*/
  if (option_data->outfilename == NULL)
    fout = stdout;
  else
    {
      /*----- see if output file already exists -----*/
      strcpy (filename, option_data->outfilename);
      fout = fopen (filename, "r");
      if (fout != NULL)
	{
	  sprintf (message, "file %s already exists. ", filename);
	  FWHM_error (message);
	}

      /*----- open file for output -----*/
      fout = fopen (filename, "w");
      if (fout == NULL)
	{
	  FWHM_error ("unable to write file ");
	}
    }

  /*----- print out the results -----*/
  fprintf (fout, "\n\n");
  fprintf (fout, "Gaussian filter widths: \n");

#if 0
  fprintf (fout, "sigmax = %5.2f   FWHMx = %5.2f \n",
  	   sx, sx * 2.0*sqrt(2.0*log(2.0)));
  fprintf (fout, "sigmay = %5.2f   FWHMy = %5.2f \n",
	   sy, sy * 2.0*sqrt(2.0*log(2.0)));
  fprintf (fout, "sigmaz = %5.2f   FWHMz = %5.2f \n\n",
	   sz, sz * 2.0*sqrt(2.0*log(2.0)));
#endif

  /* check for flag value of -1.0                      08 Mar 2004 [rickr] */
  DISP_SIG_RESULTS('x', sx);
  DISP_SIG_RESULTS('y', sy);
  DISP_SIG_RESULTS('z', sz);

  /* report any failures */
  if ( sx < 0.0 || sy < 0.0 || sz < 0.0 )
    fprintf(stderr,
	"\n"
	"** failure: some filter widths were not able to be estimated with\n"
        "            this tool, such results were reported as 'NO_VALUE'\n");

  fclose(fout);

}


/*---------------------------------------------------------------------------*/
/*
  Routine to terminate program.
*/

void terminate (input_options ** option_data, float ** fim, float ** fmask)
{
  free (*option_data);   *option_data = NULL;
  free (*fim);           *fim = NULL;
  if (*fmask != NULL)
    {  free (*fmask);       *fmask = NULL; }
}

#else   /* ENABLE_THIS_PROGRAM */

#include <stdio.h>
#include <stdlib.h>

#endif  /* ENABLE_THIS_PROGRAM */

/*---------------------------------------------------------------------------*/
/*
  Calculation of FWHM.
*/

int main (int argc, char ** argv)
{

#ifdef ENABLE_THIS_PROGRAM
  input_options * option_data = NULL;
  float * fim = NULL;
  float * fmask = NULL;
  float sx, sy, sz;
#endif


  /*----- Identify software -----*/

#ifndef ENABLE_THIS_PROGRAM
  fprintf(stderr,"\n\n") ;
  fprintf(stderr,"*** Program 3dFWHM is obsolete!  Use 3dFWHMx instead! ***\n\n\n") ;
  exit(0) ;
#else

  PRINT_VERSION("3dFWHM") ;
  AUTHOR(PROGRAM_AUTHOR);

  mainENTRY("3dFWHM main") ;
  machdep() ;

  /*----- program initialization -----*/
  initialize (argc, argv, &option_data, &fim, &fmask);

  /*----- estimate equivalent gaussian filter width -----*/
  estimate_gfw (option_data, fim, fmask, &sx, &sy, &sz );

  /*----- generate requested output -----*/
  output_results (option_data, sx, sy, sz);

  /*----- terminate program -----*/
  terminate (&option_data, &fim, &fmask);

#endif /* ENABLE_THIS_PROGRAM */

  exit(0);
}