#include "mrilib.h"

#undef  GOOD
#define GOOD(i) (mask==NULL || mask[i])

static int dontcheckplus = 0 ;
void FHWM_1dif_dontcheckplus( int i ){ dontcheckplus = i; }

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

THD_fvec3 mriarr_estimate_FWHM_1dif( MRI_IMARR *imar, byte *mask, int unif )
{
   int nar=IMARR_COUNT(imar) , ii ;
   THD_fvec3 sv ;
   float cx,cy,cz , fx,fy,fz ;
   int   nx,ny,nz , nvox , kk ;
   MRI_IMAGE *medim=NULL , *madim=NULL ;
   float     *medar=NULL , *madar=NULL , *sar=NULL ;

   unif = unif && (nar > 2) ;
   nvox = IMARR_SUBIM(imar,0)->nvox ;
   if( unif ){
     MRI_IMARR *qar = IMARR_medmad_bricks( imar ) ;
     medim = IMARR_SUBIM(qar,0); medar = MRI_FLOAT_PTR(medim);
     madim = IMARR_SUBIM(qar,1); madar = MRI_FLOAT_PTR(madim); FREE_IMARR(qar);
     for( kk=0 ; kk < nvox ; kk++ )
       if( madar[kk] != 0.0f ) madar[kk] = 1.0f / madar[kk] ;
   }
   cx = cy = cz = 0.0f ; nx = ny = nz = 0 ;
   for( ii=0 ; ii < nar ; ii++ ){
     if( unif ){
       sar = MRI_FLOAT_PTR( IMARR_SUBIM(imar,ii) ) ;
       for( kk=0 ; kk < nvox ; kk++ ) sar[kk] = (sar[kk]-medar[kk])*madar[kk] ;
     }
     sv = mri_estimate_FWHM_1dif( IMARR_SUBIM(imar,ii) , mask ) ;
     UNLOAD_FVEC3(sv,fx,fy,fz) ;
/*** INFO_message("  sub-brick[%d]: fx=%g fy=%g fz=%g",ii,fx,fy,fz) ; ***/
     if( fx > 0.0f ){ cx += fx ; nx++ ; }
     if( fy > 0.0f ){ cy += fy ; ny++ ; }
     if( fz > 0.0f ){ cz += fz ; nz++ ; }
   }
   cx = (nx==0) ? -1.0f : cx / nx ;
   cy = (ny==0) ? -1.0f : cy / ny ;
   cz = (nz==0) ? -1.0f : cz / nz ;
   LOAD_FVEC3(sv,cx,cy,cz) ;
   if( unif ){ mri_free(medim); mri_free(madim); }
   return sv ;
}

/*---------------------------------------------------------------------------*/
/*! Routine to estimate Gaussian FWHM of data brick, using first differences.
     - A negative return value indicates an error condition in that direction
       (e.g., FWHM(z) == -1.0 when nz == 1).
     - Adapted from original 3dFWHM.c by BD Ward.
     - 31 Oct 2006 -- BOO!
     - A sister function is in SUMA_GeomComp.c if changes are made to this
       one, please notify the half-gloved man.
-----------------------------------------------------------------------------*/

THD_fvec3 mri_estimate_FWHM_1dif( MRI_IMAGE *im , byte *mask )
{
  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 */
  double dx;                    /* voxel size along x-axis */
  double dy;                    /* voxel size along y-axis */
  double dz;                    /* voxel size along z-axis */
  int ix, jy, kz, ixyz2;
  double fsum, fsq, var , arg ;
  double dfdx, dfdxsum, dfdxsq, varxx;
  double dfdy, dfdysum, dfdysq, varyy;
  double dfdz, dfdzsum, dfdzsq, varzz;
  int count, countx, county, countz;

  float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
  THD_fvec3 fw_xyz ;
  MRI_IMAGE *lim ; float *fim ;

  /*----- initialize local variables -----*/

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;  /* load with error flags */

  if( im == NULL || mri_allzero(im) ) return fw_xyz ;
  lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
  fim = MRI_FLOAT_PTR(lim) ;
  nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;

  /*----- estimate the variance of the data -----*/

  fsum = 0.0; fsq = 0.0; count = 0;
  for (ixyz = 0;  ixyz < nxyz;  ixyz++){
    if( GOOD(ixyz) ){ count++; arg = fim[ixyz]; fsum += arg; fsq  += arg*arg; }
  }
  if( count < 9 || fsq <= 0.0 ){     /* no data? */
    if( lim != im ) mri_free(lim) ;
    return fw_xyz ;
  }
  var = (fsq - (fsum * fsum)/count) / (count-1.0);
  if( var <= 0.0 ){                  /* crappy data? */
    if( lim != im ) mri_free(lim) ;
    return fw_xyz ;
  }

  /*----- 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( GOOD(ixyz) ){

      arg = fim[ixyz] ;
      IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);

      if( ix+1 < nx ){
        ixyz2 = ixyz+1 ;
        if( dontcheckplus || GOOD(ixyz2) ){
          dfdx = (fim[ixyz2] - arg) ;
          dfdxsum += dfdx; dfdxsq  += dfdx * dfdx; countx++;
        }
     }

      if( jy+1 < ny ){
        ixyz2 = ixyz+nx ;
        if( dontcheckplus || GOOD(ixyz2) ){
          dfdy = (fim[ixyz2] - arg) ;
          dfdysum += dfdy; dfdysq  += dfdy * dfdy; county++;
        }
      }

      if( kz+1 < nz ){
        ixyz2 = ixyz+nxy ;
        if( dontcheckplus || GOOD(ixyz2) ){
          dfdz = (fim[ixyz2] - arg) ;
          dfdzsum += dfdz; dfdzsq  += dfdz * dfdz; countz++;
        }
      }
    }
  }

  /*----- estimate the variance of the partial derivatives -----*/

  varxx = (countx < 6) ? 0.0
                       : (dfdxsq - (dfdxsum * dfdxsum)/countx) / (countx-1.0);

  varyy = (county < 6) ? 0.0
                       : (dfdysq - (dfdysum * dfdysum)/county) / (county-1.0);

  varzz = (countz < 6) ? 0.0
                       : (dfdzsq - (dfdzsum * dfdzsum)/countz) / (countz-1.0);

  /*----- now estimate the FWHMs                                    -----*/
  /*---- the model is the correlation of the neighbors in x is
            exp[ -dx^2 / (4 * sx^2) ]
         which is the result of passing a Gaussian exp[-x^2/(2*sx^2)]
         convolution over a white noise field.  Then the covariance of
         the neighbors is
            2*V * { 1 - exp[ -dx^2 / (4*sx^2) ] }
         where V = variance of the noise (var).                     -----*/
  /*---- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor -----*/

  dx = lim->dx; dy = lim->dy; dz = lim->dz;

  arg = 1.0 - 0.5*(varxx/var);
/* ININFO_message("dx: varxx/var=%g arg=%g",varxx/var,arg) ; */
  if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) ) * dx;

  arg = 1.0 - 0.5*(varyy/var);
  if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) ) * dy;
/* ININFO_message("dy: varyy/var=%g arg=%g",varyy/var,arg) ; */

  arg = 1.0 - 0.5*(varzz/var);
  if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) ) * dz;
/* ININFO_message("dz: varzz/var=%g arg=%g",varzz/var,arg) ; */

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
  if( lim != im ) mri_free(lim) ;
  return fw_xyz ;
}

/*---------------------------------------------------------------------------*/
/*! Routine to estimate Gaussian FWHM of data brick, using differences
    between 1st and 2cd nearest neighbors.
     - A negative return value indicates an error condition in that direction
       (e.g., FWHM(z) == -1.0 when nz == 1).
-----------------------------------------------------------------------------*/

THD_fvec3 mri_estimate_FWHM_12dif( MRI_IMAGE *im , byte *mask )
{
  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 */
  double dx;                    /* voxel size along x-axis */
  double dy;                    /* voxel size along y-axis */
  double dz;                    /* voxel size along z-axis */
  int ix, jy, kz, qm,qp ;
  double dx1 , dx1sum , dx1sqq , vx1,vy1,vz1 , arg ;
  double dx2 , dx2sum , dx2sqq , vx2,vy2,vz2 ;
  double dy1 , dy1sum , dy1sqq ;
  double dy2 , dy2sum , dy2sqq ;
  double dz1 , dz1sum , dz1sqq ;
  double dz2 , dz2sum , dz2sqq ;
  int countx, county, countz;

  float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
  THD_fvec3 fw_xyz ;
  MRI_IMAGE *lim ; float *fim ;

  /*----- initialize local variables -----*/

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;  /* load with error flags */

  if( im == NULL ) return fw_xyz ;
  lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
  fim = MRI_FLOAT_PTR(lim) ;
  nx  = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;

  /*----- loop over voxels, compute differences, sum and sum squares -----*/

  countx = county = countz = 0 ;
  dx1sum = dx2sum = dy1sum = dy2sum = dz1sum = dz2sum = 0.0 ;
  dx1sqq = dx2sqq = dy1sqq = dy2sqq = dz1sqq = dz2sqq = 0.0 ;
  for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
    if( GOOD(ixyz) ){
      arg = fim[ixyz] ; IJK_TO_THREE (ixyz, ix,jy,kz, nx,nxy);

      if( ix-1 >= 0 && ix+1 < nx ){
        qp = ixyz+1 ; qm = ixyz-1 ;
        if( GOOD(qp) && GOOD(qm) ){
          dx1     = fim[qp]-arg ; dx2     = fim[qp]-fim[qm] ;
          dx1sum += dx1         ; dx2sum += dx2             ;
          dx1sqq += dx1*dx1     ; dx2sqq += dx2*dx2         ;
          countx++ ;
        }
      }

      if( jy-1 >= 0 && jy+1 < ny ){
        qp = ixyz+nx ; qm = ixyz-nx ;
        if( GOOD(qp) && GOOD(qm) ){
          dy1     = fim[qp]-arg ; dy2     = fim[qp]-fim[qm] ;
          dy1sum += dy1         ; dy2sum += dy2             ;
          dy1sqq += dy1*dy1     ; dy2sqq += dy2*dy2         ;
          county++ ;
        }
      }

      if( kz-1 >= 0 && kz+1 < nz ){
        qp = ixyz+nxy ; qm = ixyz-nxy ;
        if( GOOD(qp) && GOOD(qm) ){
          dz1     = fim[qp]-arg ; dz2     = fim[qp]-fim[qm] ;
          dz1sum += dz1         ; dz2sum += dz2             ;
          dz1sqq += dz1*dz1     ; dz2sqq += dz2*dz2         ;
          countz++ ;
        }
      }
    }
  }

  /*----- estimate variances of differences -----*/

  vx1 = (countx < 6) ? 0.0 : (dx1sqq - (dx1sum*dx1sum)/countx) / (countx-1.0) ;
  vy1 = (county < 6) ? 0.0 : (dy1sqq - (dy1sum*dy1sum)/county) / (county-1.0) ;
  vz1 = (countz < 6) ? 0.0 : (dz1sqq - (dz1sum*dz1sum)/countz) / (countz-1.0) ;

  vx2 = (countx < 6) ? 0.0 : (dx2sqq - (dx2sum*dx2sum)/countx) / (countx-1.0) ;
  vy2 = (county < 6) ? 0.0 : (dy2sqq - (dy2sum*dy2sum)/county) / (county-1.0) ;
  vz2 = (countz < 6) ? 0.0 : (dz2sqq - (dz2sum*dz2sum)/countz) / (countz-1.0) ;

#if 0
fprintf(stderr,"countx=%d dx1sum=%g dx1sqq=%g vx1=%g\n",countx,dx1sum,dx1sqq,vx1) ;
fprintf(stderr,"countx=%d dx2sum=%g dx2sqq=%g vx2=%g\n",countx,dx2sum,dx2sqq,vx2) ;
#endif

  /*----- now estimate the FWHMs -----*/

  dx = lim->dx; dy = lim->dy; dz = lim->dz;
  if( lim != im ) mri_free(lim) ;

  /*--- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor ---*/
  /*--- y = cbrt(12*sqrt(48-120*r+81*r*r)+108*r-80), and then
        x = y/6 - 4/(3*y) - 1/3
        is the real solution to the equation (1-x^4)/(1-x) = r > 1 ;
        here, r = vx2/vx1 = ratio of variances at Delta=2*dx vs. Delta=1*dx;
        x = exp[-dx**2/(4*sigma**2)] = correlation coefficient of neighbors;
        we solve for x, then use that to solve for sigma, and scale to FWHM --*/

  if( vx1 > 0.0 && vx2 > vx1 ){
    arg = vx2 / vx1 ;
    arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
    arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
    if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dx ;
  }

  if( vy1 > 0.0 && vy2 > vy1 ){
    arg = vy2 / vy1 ;
    arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
    arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
    if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dy ;
  }

  if( vz1 > 0.0 && vz2 > vz1 ){
    arg = vz2 / vz1 ;
    arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
    arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
    if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dz ;
  }

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
  return fw_xyz ;
}

/*---------------------------------------------------------------------------*/
/*! Routine to estimate Gaussian FWHM of data brick, using differences
    between 1st and 2cd nearest neighbors.
     - A negative return value indicates an error condition in that direction
       (e.g., FWHM(z) == -1.0 when nz == 1).
-----------------------------------------------------------------------------*/

THD_fvec3 mri_estimate_FWHM_12dif_MAD( MRI_IMAGE *im , byte *mask )
{
  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, qm,qp ;
  float vx1,vy1,vz1 , arg ;
  float vx2,vy2,vz2 ;
  float wx1,wy1,wz1 , brg ;
  float wx2,wy2,wz2 ;
  int countx, county, countz;

  int ngood ; float *dx1ar,*dy1ar,*dz1ar , *dx2ar,*dy2ar,*dz2ar ;

  float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
  THD_fvec3 fw_xyz ;
  MRI_IMAGE *lim ; float *fim ;

  /*----- initialize local variables -----*/

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;  /* load with error flags */

  if( im == NULL ) return fw_xyz ;
  lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
  fim = MRI_FLOAT_PTR(lim) ;
  nx  = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;

  for( ngood=ixyz=0 ; ixyz < nxyz ; ixyz++ ) if( GOOD(ixyz) ) ngood++ ;
  if( ngood < 9 ) return fw_xyz ;

  dx1ar = (float *)malloc(sizeof(float)*ngood) ;
  dy1ar = (float *)malloc(sizeof(float)*ngood) ;
  dz1ar = (float *)malloc(sizeof(float)*ngood) ;
  dx2ar = (float *)malloc(sizeof(float)*ngood) ;
  dy2ar = (float *)malloc(sizeof(float)*ngood) ;
  dz2ar = (float *)malloc(sizeof(float)*ngood) ;

  /*----- loop over voxels, compute differences, sum and sum squares -----*/

  countx = county = countz = 0 ;
  for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
    if( GOOD(ixyz) ){
      arg = fim[ixyz] ; IJK_TO_THREE (ixyz, ix,jy,kz, nx,nxy);

      if( ix-1 >= 0 && ix+1 < nx ){
        qp = ixyz+1 ; qm = ixyz-1 ;
        if( GOOD(qp) && GOOD(qm) ){
          dx1ar[countx] = fim[qp]-arg ; dx2ar[countx] = fim[qp]-fim[qm] ; countx++ ;
        }
      }

      if( jy-1 >= 0 && jy+1 < ny ){
        qp = ixyz+nx ; qm = ixyz-nx ;
        if( GOOD(qp) && GOOD(qm) ){
          dy1ar[county] = fim[qp]-arg ; dy2ar[county] = fim[qp]-fim[qm] ; county++ ;
        }
      }

      if( kz-1 >= 0 && kz+1 < nz ){
        qp = ixyz+nxy ; qm = ixyz-nxy ;
        if( GOOD(qp) && GOOD(qm) ){
          dz1ar[countz] = fim[qp]-arg ; dz2ar[countz] = fim[qp]-fim[qm] ; countz++ ;
        }
      }
    }
  }

  /*----- estimate variances of differences -----*/

  qmedmadmeanad_float( countx, dx1ar, NULL, &vx1, &wx1 ); vx1 = vx1*vx1; wx1 = wx1*wx1 ;
  qmedmadmeanad_float( county, dy1ar, NULL, &vy1, &wy1 ); vy1 = vy1*vy1; wy1 = wy1*wy1 ;
  qmedmadmeanad_float( countz, dz1ar, NULL, &vz1, &wz1 ); vz1 = vz1*vz1; wz1 = wz1*wz1 ;
  qmedmadmeanad_float( countx, dx2ar, NULL, &vx2, &wx2 ); vx2 = vx2*vx2; wx2 = wx2*wx2 ;
  qmedmadmeanad_float( county, dy2ar, NULL, &vy2, &wy2 ); vy2 = vy2*vy2; wy2 = wy2*wy2 ;
  qmedmadmeanad_float( countz, dz2ar, NULL, &vz2, &wz2 ); vz2 = vz2*vz2; wz2 = wz2*wz2 ;

  free(dx1ar) ; free(dy1ar) ; free(dz1ar) ;
  free(dx2ar) ; free(dy2ar) ; free(dz2ar) ;

  /*----- now estimate the FWHMs -----*/

  dx = lim->dx; dy = lim->dy; dz = lim->dz;
  if( lim != im ) mri_free(lim) ;

  /*--- 2.35482 = sqrt(8*log(2)) = sigma-to-FWHM conversion factor ---*/
  /*--- y = cbrt(12*sqrt(48-120*r+81*r*r)+108*r-80), and then
        x = y/6 - 4/(3*y) - 1/3
        is the real solution to the equation (1-x^4)/(1-x) = r > 1 ;
        here, r = vx2/vx1 = ratio of variances at Delta=2*dx vs. Delta=1*dx;
        x = exp[-dx**2/(4*sigma**2)] = correlation coefficient of neighbors;
        we solve for x, then use that to solve for sigma, and scale to FWHM --*/

  /*--- Modified 11 Aug 2015, to use the mean absolute deviations (wx1 & wx2)
        instead of the median absolute deviations (vx1 & vx2) if the
        computation fails for the median absolute deviations (the backup) ---*/

  if( vx1 > 0.0 && vx2 > vx1 ){
    arg = vx2 / vx1 ; brg = wx2 / wx1 ;
    arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
    arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
    if( arg > 0.0 && arg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dx ;
    else if( brg > 1.0f ){
      brg = cbrt(12.0*sqrt(48.0 - 120.0*brg + 81.0*brg*brg) + 108.0*brg - 80.0) ;
      brg = brg/6.0 - 4.0/(3.0*brg) - 1.0/3.0 ;
      if( brg > 0.0 && brg < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(brg)) )*dx ;
    }
/* INFO_message("FWHM x: vx1=%f vx2=%f arg=%f sx=%f countx=%d",vx1,vx2,arg,sx,countx) ; */
  } else {
/* INFO_message("FWHM x: vx1=%f vx2=%f NOT GOOD",vx1,vx2) ; */
  }

  if( vy1 > 0.0 && vy2 > vy1 ){
    arg = vy2 / vy1 ; brg = wy2 / wy1 ;
    arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
    arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
    if( arg > 0.0 && arg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dy ;
    else if( brg > 1.0f ){
      brg = cbrt(12.0*sqrt(48.0 - 120.0*brg + 81.0*brg*brg) + 108.0*brg - 80.0) ;
      brg = brg/6.0 - 4.0/(3.0*brg) - 1.0/3.0 ;
      if( brg > 0.0 && brg < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(brg)) )*dy ;
    }
/* INFO_message("FWHM y: vy1=%f vy2=%f arg=%f sy=%f county=%d",vy1,vy2,arg,sy,county) ; */
  } else {
/* INFO_message("FWHM y: vy1=%f vy2=%f NOT GOOD",vy1,vy2) ; */
  }

  if( vz1 > 0.0 && vz2 > vz1 ){
    arg = vz2 / vz1 ; brg = wz2 / wz1 ;
    arg = cbrt(12.0*sqrt(48.0 - 120.0*arg + 81.0*arg*arg) + 108.0*arg - 80.0) ;
    arg = arg/6.0 - 4.0/(3.0*arg) - 1.0/3.0 ;
    if( arg > 0.0 && arg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(arg)) )*dz ;
    else if( brg > 1.0f ){
      brg = cbrt(12.0*sqrt(48.0 - 120.0*brg + 81.0*brg*brg) + 108.0*brg - 80.0) ;
      brg = brg/6.0 - 4.0/(3.0*brg) - 1.0/3.0 ;
      if( brg > 0.0 && brg < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(brg)) )*dz ;
    }
/* INFO_message("FWHM z: vz1=%f vz2=%f arg=%f sz=%f countz=%d",vz1,vz2,arg,sz,countz) ; */
  } else {
/* INFO_message("FWHM z: vz1=%f vz2=%f NOT GOOD",vz1,vz2) ; */
  }

/* INFO_message("12dif_MAD: sxyz = %g %g %g",sx,sy,sz) ; */
  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
/* INFO_message("12dif_MAD: fw = %g %g %g",fw_xyz.xyz[0],fw_xyz.xyz[1],fw_xyz.xyz[2]) ; */
  return fw_xyz ;
}

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

static THD_fvec3 (*fester)(MRI_IMAGE *, byte *) = mri_estimate_FWHM_1dif ;

void mri_fwhm_setfester( THD_fvec3 (*func)(MRI_IMAGE *, byte *) )
{
  if( func == NULL ) func = mri_estimate_FWHM_1dif ;
  fester = func ;
}

/*---------------------------------------------------------------------------*/
/*! Get FWHM estimates for each sub-brick in a dataset.
    Output image is 3xN where N=# of sub-bricks.
-----------------------------------------------------------------------------*/

MRI_IMAGE * THD_estimate_FWHM_all( THD_3dim_dataset *dset,
                                   byte *mask, int demed , int unif )
{
   int nvals , ii,nvox ;
   MRI_IMAGE *outim=NULL , *medim=NULL , *madim=NULL ;
   float     *outar,fac  , *medar=NULL , *madar=NULL ;

ENTRY("THD_estimate_FWHM_all") ;

   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
   DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;

   nvals = DSET_NVALS(dset) ;
   outim = mri_new( 3 , nvals , MRI_float ) ;
   outar = MRI_FLOAT_PTR(outim) ;
   nvox  = DSET_NVOX(dset) ;

   if( unif ){
     MRI_IMARR *imar ;
     demed = 1 ;
     imar  = THD_medmad_bricks(dset) ;
     medim = IMARR_SUBIM(imar,0) ; medar = MRI_FLOAT_PTR(medim) ;
     madim = IMARR_SUBIM(imar,1) ; madar = MRI_FLOAT_PTR(madim) ;
     FREE_IMARR(imar) ;
     for( ii=0 ; ii < nvox ; ii++ )
       if( madar[ii] > 0.0f ) madar[ii] = 1.0f / madar[ii] ;
   } else if( demed ){
     medim = THD_median_brick(dset) ;
     medar = MRI_FLOAT_PTR(medim) ;
   }

AFNI_OMP_START;
#pragma omp parallel if( nvals > 4 )
 { int iv,ii ; MRI_IMAGE *bim ; float *bar ; THD_fvec3 fw ;
#pragma omp for
   for( iv=0 ; iv < nvals ; iv++ ){
     if( mri_allzero(DSET_BRICK(dset,iv)) ){
       outar[0+3*iv] = outar[1+3*iv] = outar[2+3*iv] = 0.0f ; continue ;
     }
     bim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,iv), DSET_BRICK(dset,iv) );
     if( demed ){
       bar = MRI_FLOAT_PTR(bim) ;
       for( ii=0 ; ii < nvox ; ii++ ) bar[ii] -= medar[ii] ;
       if( unif )
       for( ii=0 ; ii < nvox ; ii++ ) bar[ii] *= madar[ii] ;
     }
     fw = fester( bim , mask ) ; mri_free(bim) ;
     UNLOAD_FVEC3( fw , outar[0+3*iv] , outar[1+3*iv] , outar[2+3*iv] ) ;
   }
 }
AFNI_OMP_END;

   if( demed ) mri_free(medim) ;
   if( unif  ) mri_free(madim) ;

   RETURN(outim) ;
}

/*========================================================================*/
static double mom12_stdev_fac = 1.0 ;

void mri_fwhm_mom12_set_stdev_fac( double fff ){ mom12_stdev_fac = fff ; }

/*------------------------------------------------------------------------*/
/*! New method using first and second moments of absolute values of
    differences, instead of only second moments.  [05 Aug 2015]
*//*----------------------------------------------------------------------*/

THD_fvec3 mri_FWHM_1dif_mom12( MRI_IMAGE *im , byte *mask )
{
  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 */
  double dx;                    /* voxel size along x-axis */
  double dy;                    /* voxel size along y-axis */
  double dz;                    /* voxel size along z-axis */
  int ix, jy, kz, ixyz2;
  double fsum, fsq, var , vfac , arg ;
  double dfdx, dfdy, dfdz ;
  double dfdxs1, dfdxs2, dfdxs3, dfdxs4, dfdxna ;
  double dfdys1, dfdys2, dfdys3, dfdys4, dfdyna ;
  double dfdzs1, dfdzs2, dfdzs3, dfdzs4, dfdzna ;
  int count=0, countx=0, county=0, countz=0 ;

  float sx=-1.0f,sy=-1.0f,sz=-1.0f ;
  MRI_IMAGE *lim ; float *fim ;
  THD_fvec3 fw_xyz ;

  /*----- initialize local variables -----*/

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;  /* load with error flags */

  if( im == NULL || mri_allzero(im) ) return fw_xyz ;
  lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
  fim = MRI_FLOAT_PTR(lim) ;
  nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;

  /*----- estimate the variance of the data itself -----*/

  fsum = 0.0; fsq = 0.0; count = 0;
  for (ixyz = 0;  ixyz < nxyz;  ixyz++){
    if( GOOD(ixyz) ){ count++; arg = fim[ixyz]; fsum += arg; fsq  += arg*arg; }
  }
  if( count < 9 || fsq <= 0.0 ){     /* no data? */
    if( lim != im ) mri_free(lim) ;
    return fw_xyz ;
  }
  var = (fsq - (fsum * fsum)/count) / (count-1.0);
  if( var <= 0.0 ){                  /* crappy data? */
    if( lim != im ) mri_free(lim) ;
    return fw_xyz ;
  }
  vfac = 1.0 / sqrt(2.0*var) ;

  /*----- estimate moments of the partial derivatives -----*/
  /*((((( 3rd and 4th moments are not used at present )))))*/

  dfdxs1 = dfdxs2 = dfdxs3 = dfdxs4 = 0.0 ; dfdxna = 0.0 ;
  dfdys1 = dfdys2 = dfdys3 = dfdys4 = 0.0 ; dfdyna = 0.0 ;
  dfdzs1 = dfdzs2 = dfdzs3 = dfdzs4 = 0.0 ; dfdzna = 0.0 ;
  countx = county = countz = 0 ;

  for (ixyz = 0;  ixyz < nxyz;  ixyz++){
    if( GOOD(ixyz) ){

      arg = fim[ixyz] ;
      IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);

      if( ix+1 < nx ){
        ixyz2 = ixyz+1 ;
        if( dontcheckplus || GOOD(ixyz2) ){
          dfdx = vfac*(fim[ixyz2] - arg) ; dfdxna += dfdx ; dfdx = fabs(dfdx) ;
          dfdxs1 += dfdx; dfdxs2 += dfdx*dfdx;
          dfdxs3 += dfdx*dfdx*dfdx; dfdxs4 += dfdx*dfdx*dfdx*dfdx;
          countx++;
        }
     }

      if( jy+1 < ny ){
        ixyz2 = ixyz+nx ;
        if( dontcheckplus || GOOD(ixyz2) ){
          dfdy = vfac*(fim[ixyz2] - arg) ; dfdyna += dfdy ; dfdy = fabs(dfdy) ;
          dfdys1 += dfdy; dfdys2 += dfdy*dfdy;
          dfdys3 += dfdy*dfdy*dfdy; dfdys4 += dfdy*dfdy*dfdy*dfdy;
          county++;
        }
      }

      if( kz+1 < nz ){
        ixyz2 = ixyz+nxy ;
        if( dontcheckplus || GOOD(ixyz2) ){
          dfdz = vfac*(fim[ixyz2] - arg) ; dfdzna += dfdz ; dfdz = fabs(dfdz) ;
          dfdzs1 += dfdz; dfdzs2 += dfdz*dfdz;
          dfdzs3 += dfdz*dfdz*dfdz; dfdzs4 += dfdz*dfdz*dfdz*dfdz;
          countz++;
        }
      }
    }
  }
  if( countx == 0 ) countx = 1 ;  /* patch for possible stupidity in data */
  if( county == 0 ) county = 1 ;
  if( county == 0 ) county = 1 ;

  dfdxs1 /= countx; dfdxs2 /= countx; dfdxs3 /= countx; dfdxs4 /= countx;
  dfdys1 /= county; dfdys2 /= county; dfdys3 /= county; dfdys4 /= county;
  dfdzs1 /= countz; dfdzs2 /= countz; dfdzs3 /= countz; dfdzs4 /= countz;

  dfdxna /= countx ; dfdyna /= county ; dfdzna /= countz ;

  dx = lim->dx; dy = lim->dy; dz = lim->dz;
  if( lim != im ) mri_free(lim) ;

#if 0
  ININFO_message(" var= %13.6g",var) ;
  ININFO_message(" dx:  %13.6g  %13.6g  %13.6g  %13.6g  na = %13.6g  #x = %d",
                 dfdxs1 , dfdxs2 , dfdxs3 , dfdxs4 , dfdxna , countx ) ;
  ININFO_message(" dy:  %13.6g  %13.6g  %13.6g  %13.6g  na = %13.6g  #y = %d",
                 dfdys1 , dfdys2 , dfdys3 , dfdys4 , dfdyna , county ) ;
  ININFO_message(" dz:  %13.6g  %13.6g  %13.6g  %13.6g  na = %13.6g  #z = %d",
                 dfdzs1 , dfdzs2 , dfdzs3 , dfdzs4 , dfdzna , countz ) ;
#endif

#undef  SQHPI
#define SQHPI 1.253314  /* sqrt(PI/2) */
#undef  SQEPI
#define SQEPI 0.626657  /* sqrt(PI/8) -- needed for 3rd moment */

  /*--------------------------------------------------------------------------*/
  /*   Model: first differences are distributed like
         ef = Normal( 0 , 2*V*tau^2 )
       where tau = sqrt{ 1 - exp[-del^2/(4*s^2)] }
             V    = variance (var) of data
             del  = spacing of data in mm
             s    = sigma of smoothing kernel
                  = del * sqrt( -1/(4*log(1-tau^2))
                    [s(tau) is a monotonic decreasing function of tau]
             FWHM = sqrt(8*log(2))*s = 2.35482*s
       and we assume that tau in turn is (independently) random, distributed
       with density w(tau) over the interval [0,1].  That is, the density
       function of df = ef/sqrt(2*V) is
             |1
          int|  w(tau) 1/[tau*sqrt(2*pi)] * exp[-df^2/(2*tau^2)] d(tau)
             |0
       Unless w(tau) is a delta function, df will be distributed non-normally.
       Then the moments of the absolute value of df are related to the moments
       of tau by        |1
         E[|df|^p] = int|  w(tau) E[ |Normal(0,tau^2)|^p ] d(tau)
                        |0
                   = 2^(p/2) * Gamma((p+1)/2) / sqrt(pi) * E[tau^p]
         E[|df|]   = sqrt(2/pi) * E[tau]
         E[|df|^2] =              E[tau^2]
         E[|df|^3] = sqrt(8/pi) * E[tau^3]
         E[|df|^4] =          3 * E[tau^4] etc.
       (At this time, we don't use the 3rd or 4th moments, but someday ...?)
       The "standard" (Forman) way to estimate FWHN (mri_estimate_FWHM_1dif)
       assumes s is constant; e.g., that w(tau) = delta(tau-q), in which
       case E[|df|^2] = q^2, and so only the second moment is needed to
       estimate tau=q (and thus FWHM).  In this routine, we use E[|df|] to
       estimate E[tau], and then subtract off one standard deviation
       of tau sqrt(E[tau^2]-E[tau]^2) to weight the estimate towards
       smaller tau and thus larger FWHM -- since larger FWHM will produce
       more "false positive" clusters. Completely ad hoc, but it's something. */
  /*--------------------------------------------------------------------------*/

  { double mn, sg , ff ;
    mn = SQHPI*dfdxs1 ;
    sg = dfdxs2 - mn*mn ;
    if( mn > 0.0 && sg > 0.0 ){
      sg = sqrt(sg) ;
      for( ff=mom12_stdev_fac ; sg*ff >= mn ; ff *= 0.888 ) ; /*nada*/
      mn -= ff*sg ;
    }
    mn *= mn ;
    if( mn > 0.0 && mn < 1.0 ) sx = 2.35482*sqrt( -1.0 / (4.0*log(1.0-mn)) ) * dx;

    mn = SQHPI*dfdys1 ;
    sg = dfdys2 - mn*mn ;
    if( mn > 0.0 && sg > 0.0 ){
      sg = sqrt(sg) ;
      for( ff=mom12_stdev_fac ; sg*ff >= mn ; ff *= 0.888 ) ; /*nada*/
      mn -= ff*sg ;
    }
    mn *= mn ;
    if( mn > 0.0 && mn < 1.0 ) sy = 2.35482*sqrt( -1.0 / (4.0*log(1.0-mn)) ) * dy;

    mn = SQHPI*dfdzs1 ;
    sg = dfdzs2 - mn*mn ;
    if( mn > 0.0 && sg > 0.0 ){
      sg = sqrt(sg) ;
      for( ff=mom12_stdev_fac ; sg*ff >= mn ; ff *= 0.888 ) ; /*nada*/
      mn -= ff*sg ;
    }
    mn *= mn ;
    if( mn > 0.0 && mn < 1.0 ) sz = 2.35482*sqrt( -1.0 / (4.0*log(1.0-mn)) ) * dz;
  }

  LOAD_FVEC3(fw_xyz,sx,sy,sz) ;
  return fw_xyz ;
}

/*========================================================================*/
/*========================================================================*/
/* Return the cluster of points over which we estimate the ACF.
   Use the central NCLU_BASE set of points, and then fill out
   up to NCLU_GOAL from a selection of the outer points,
   so that we don't end up with a vast cluster and slow estimation.
*//*----------------------------------------------------------------------*/

#undef  NCLU_GOAL
#define NCLU_GOAL 666
#undef  NCLU_BASE
#define NCLU_BASE 111

static int do_2D_ACF = 0 ;

void set_ACF_2D( int nn ){ do_2D_ACF = nn ; }

static MCW_cluster * get_ACF_cluster( float dx, float dy, float dz, float radius )
{
   MCW_cluster *clout=NULL , *cltemp=NULL ;
   int pp , dp ;

   /* a ball of desired radius */

   cltemp = MCW_spheremask( dx,dy,dz , radius ) ;

   /* if not too big, just use this */

   if( cltemp->num_pt < NCLU_GOAL ) return cltemp ;

   /* keep the central NCLU_BASE points */

   INIT_CLUSTER(clout) ;
   for( pp=0 ; clout->num_pt < NCLU_BASE && pp < cltemp->num_pt ; pp++ ){
     if( do_2D_ACF && cltemp->k[pp] != 0 ) continue ;
     ADDTO_CLUSTER(clout,cltemp->i[pp],cltemp->j[pp],cltemp->k[pp],0.0f) ;
   }

   /* keep a selection of the outer points, so we don't get too many */

   dp = (int)rintf( (cltemp->num_pt-NCLU_BASE) / (float)(NCLU_GOAL-NCLU_BASE) ) ;
   for( ; pp < cltemp->num_pt ; ){
     if( do_2D_ACF && cltemp->k[pp] != 0 ){
       pp++ ;
     } else {
       ADDTO_CLUSTER(clout,cltemp->i[pp],cltemp->j[pp],cltemp->k[pp],0.0f) ;
       pp += dp ;
     }
   }

   KILL_CLUSTER(cltemp) ;
   MCW_radsort_cluster(clout,dx,dy,dz) ;
   return clout ;
}
/*----------------------------------------------------------------------------*/
/* Estimate the shape of the spatial autocorrelation function   [09 Nov 2015] */

int mri_estimate_ACF( MRI_IMAGE *im , byte *mask , MCW_cluster *clout )
{
   float dx,dy,dz ;
   MRI_IMAGE *lim ; float *fim ;
   int nx,ny,nz,nxy,nxyz , ixyz,ip,jp,kp,ppp,qqq , ix,jy,kz ;
   double fsum, fsq, fbar, fvar , arg ;
   int count , nclu ;
   double *acfsqq ;
   int    *nacf ;

ENTRY("mri_estimate_ACF") ;

   if( im    == NULL || mri_allzero(im)   ) RETURN(-1) ;
   if( clout == NULL || clout->num_pt < 5 ) RETURN(-1) ;

   for( ppp=0 ; ppp < clout->num_pt ; ppp++ ) clout->mag[ppp] = 0.0f ;

   /*-- setup to process the input image --*/

   lim = (im->kind == MRI_float) ? im : mri_to_float(im) ;
   fim = MRI_FLOAT_PTR(lim) ;
   nx = lim->nx; ny = lim->ny; nz = lim->nz; nxy = nx*ny; nxyz = nx*ny*nz;

   /*----- estimate the variance of the data -----*/

   fsum = 0.0; fsq = 0.0; count = 0;
   for (ixyz = 0;  ixyz < nxyz;  ixyz++){
     if( GOOD(ixyz) ){ count++; arg = fim[ixyz]; fsum += arg; fsq += arg*arg; }
   }
   if( count < 9 || fsq <= 0.0 ){     /* no data? */
     if( lim != im ) mri_free(lim) ;
     RETURN(-1) ;
   }
   fbar = fsum / count ;
   fvar = (fsq - (fsum * fsum)/count) / (count-1.0);
   if( fvar <= 0.0 ){
     if( lim != im ) mri_free(lim) ;
     RETURN(-1) ;
   }

   /*--- space to accumulate results ---*/

   nclu   = clout->num_pt ;
   acfsqq = (double *)calloc(sizeof(double),nclu) ;
   nacf   = (int    *)calloc(sizeof(int)   ,nclu) ;

   /*--- loop over voxels, sum up cross products ---*/

   for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){

     if( !GOOD(ixyz) ) continue ;  /* not in mask */

     arg = fim[ixyz] - fbar ;
     IJK_TO_THREE(ixyz,ix,jy,kz,nx,nxy) ;

     /* loop over offsets */

     { int ppp , ip,jp,kp,qqq ;
       for( ppp=1 ; ppp < nclu ; ppp++ ){
         ip = ix + clout->i[ppp] ; if( ip >= nx || ip < 0 ) continue ;
         jp = jy + clout->j[ppp] ; if( jp >= ny || jp < 0 ) continue ;
         kp = kz + clout->k[ppp] ; if( kp >= nz || kp < 0 ) continue ;
         qqq = THREE_TO_IJK(ip,jp,kp,nx,nxy) ; if( !GOOD(qqq) ) continue ;
         acfsqq[ppp] += (fim[qqq]-fbar)*arg ; nacf[ppp]++ ;
       }
     }

   }

   /*--- load results into clout ---*/

   for( ppp=1 ; ppp < nclu ; ppp++ ){
     if( nacf[ppp] > 5 )
       clout->mag[ppp] = (float)( acfsqq[ppp] / ( fvar * (nacf[ppp]-1.0) ) ) ;
   }
   clout->mag[0] = 1.0f ;

   free(acfsqq) ; free(nacf) ;
   if( lim != im ) mri_free(lim) ;
   RETURN(0) ;
}

/*----------------------------------------------------------------------------*/
/*! Get average ACF estimate for all sub-bricks                  [09 Nov 2015]
*//*--------------------------------------------------------------------------*/

MCW_cluster * THD_estimate_ACF( THD_3dim_dataset *dset,
                                byte *mask, int demed, int unif, float radius )
{
   int iv , nvals , ii,nvox , nout ;
   MRI_IMAGE *medim=NULL , *madim=NULL ;
   float     *medar=NULL , *madar=NULL ;
   MCW_cluster *clout=NULL , **cltemp ;

ENTRY("THD_estimate_ACF") ;

   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
   DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;

   nvals = DSET_NVALS(dset) ;
   nvox  = DSET_NVOX(dset) ;

   if( unif ){
     MRI_IMARR *imar ;
     demed = 1 ;
     imar  = THD_medmad_bricks(dset) ;
     medim = IMARR_SUBIM(imar,0) ; medar = MRI_FLOAT_PTR(medim) ;
     madim = IMARR_SUBIM(imar,1) ; madar = MRI_FLOAT_PTR(madim) ;
     FREE_IMARR(imar) ;
     for( ii=0 ; ii < nvox ; ii++ )
       if( madar[ii] > 0.0f ) madar[ii] = 1.0f / madar[ii] ;
   } else if( demed ){
     medim = THD_median_brick(dset) ;
     medar = MRI_FLOAT_PTR(medim) ;
   }

   set_ACF_2D( DSET_NZ(dset)==1 ) ;

   clout = get_ACF_cluster( fabsf(DSET_DX(dset)) , fabsf(DSET_DY(dset)) ,
                            fabsf(DSET_DZ(dset)) , radius ) ;
   if( clout == NULL ){
     if( medim != NULL ) mri_free(medim) ;
     if( madim != NULL ) mri_free(madim) ;
     RETURN(NULL) ;
   }

   for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] = 0.0f ;
   cltemp = (MCW_cluster **)malloc(sizeof(MCW_cluster *)*nvals) ;

AFNI_OMP_START;
#pragma omp parallel if( nvals > 4 )
 { int iv,gg ; MRI_IMAGE *bim ;
#pragma omp for
   for( iv=0 ; iv < nvals ; iv++ ){
     bim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,iv), DSET_BRICK(dset,iv) );
     bim->dx = DSET_DX(dset) ; bim->dy = DSET_DY(dset) ; bim->dz = DSET_DZ(dset) ;
     if( demed ){
       float *bar = MRI_FLOAT_PTR(bim) ; int ii ;
       for( ii=0 ; ii < nvox ; ii++ ) bar[ii] -= medar[ii] ;
       if( unif )
        for( ii=0 ; ii < nvox ; ii++ ) bar[ii] *= madar[ii] ;
     }
     COPY_CLUSTER(cltemp[iv],clout) ;
     gg = mri_estimate_ACF( bim , mask , cltemp[iv] ) ;
     if( gg < 0 ){
       KILL_CLUSTER(cltemp[iv]) ; cltemp[iv] = NULL ;
     }
   }
 }
AFNI_OMP_END;

   for( nout=iv=0 ; iv < nvals ; iv++ ){
     if( cltemp[iv] != NULL ){
       for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] += cltemp[iv]->mag[ii] ;
       KILL_CLUSTER(cltemp[iv]) ; nout++ ;
     }
   }
   free(cltemp) ;

   if( medim != NULL ) mri_free(medim) ;
   if( madim != NULL ) mri_free(madim) ;

   if( nout > 1 ){
     for( ii=0 ; ii < clout->num_pt ; ii++ )
       clout->mag[ii] /= nout ;
   } else if( nout == 0 ){
     ERROR_message("ACF estimation fails :-(") ;
     KILL_CLUSTER(clout) ; clout = NULL ;
   }

   RETURN(clout) ;
}

/*----------------------------------------------------------------------------*/
/*! Get average ACF FWHM estimate for all images in the array    [30 Dec 2015]
*//*--------------------------------------------------------------------------*/

float mriarr_estimate_FWHM_acf( MRI_IMARR *imar, byte *mask, int unif, float radius )
{
   int iv , nvals , ii,nvox , nout , demed=0 ;
   MRI_IMAGE *medim=NULL , *madim=NULL ;
   float     *medar=NULL , *madar=NULL , dx,dy,dz , fwhm_out=0.0f ;
   MCW_cluster *clout=NULL , **cltemp ;

ENTRY("mriarr_estimate_FWHM_acf") ;

   if( imar == NULL ) RETURN(0.0f) ;

   nvals = IMARR_COUNT(imar) ; if( nvals < 2 ) unif = 0 ;
   nvox  = IMARR_SUBIM(imar,0)->nvox ;
   dx    = IMARR_SUBIM(imar,0)->dx ;
   dy    = IMARR_SUBIM(imar,0)->dy ;
   dz    = IMARR_SUBIM(imar,0)->dz ;

   if( unif ){
     MRI_IMARR *qmar ;
     demed = 1 ;
     qmar  = IMARR_medmad_bricks(imar) ;
     medim = IMARR_SUBIM(qmar,0) ; medar = MRI_FLOAT_PTR(medim) ;
     madim = IMARR_SUBIM(qmar,1) ; madar = MRI_FLOAT_PTR(madim) ;
     FREE_IMARR(qmar) ;
     for( ii=0 ; ii < nvox ; ii++ )
       if( madar[ii] > 0.0f ) madar[ii] = 1.0f / madar[ii] ;
   }

   set_ACF_2D( IMARR_SUBIM(imar,0)->nz == 1 ) ;

   clout = get_ACF_cluster( dx,dy,dz, radius ) ;
   if( clout == NULL ){
     if( medim != NULL ) mri_free(medim) ;
     if( madim != NULL ) mri_free(madim) ;
     RETURN(0.0f) ;
   }

   for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] = 0.0f ;
   cltemp = (MCW_cluster **)malloc(sizeof(MCW_cluster *)*nvals) ;

AFNI_OMP_START;
#pragma omp parallel if( nvals > 4 )
 { int iv,gg ; MRI_IMAGE *bim ;
#pragma omp for
   for( iv=0 ; iv < nvals ; iv++ ){
     bim = mri_to_float( IMARR_SUBIM(imar,iv) ) ;
     bim->dx = dx ; bim->dy = dy ; bim->dz = dz ;
     if( demed ){
       float *bar = MRI_FLOAT_PTR(bim) ; int ii ;
       for( ii=0 ; ii < nvox ; ii++ ) bar[ii] -= medar[ii] ;
       if( unif )
        for( ii=0 ; ii < nvox ; ii++ ) bar[ii] *= madar[ii] ;
     }
     COPY_CLUSTER(cltemp[iv],clout) ;
     gg = mri_estimate_ACF( bim , mask , cltemp[iv] ) ;
     if( gg < 0 ){
       KILL_CLUSTER(cltemp[iv]) ; cltemp[iv] = NULL ;
     }
   }
 }
AFNI_OMP_END;

   for( nout=iv=0 ; iv < nvals ; iv++ ){
     if( cltemp[iv] != NULL ){
       for( ii=0 ; ii < clout->num_pt ; ii++ ) clout->mag[ii] += cltemp[iv]->mag[ii] ;
       KILL_CLUSTER(cltemp[iv]) ; nout++ ;
     }
   }
   free(cltemp) ;

   if( medim != NULL ) mri_free(medim) ;
   if( madim != NULL ) mri_free(madim) ;

   if( nout > 1 ){
     for( ii=0 ; ii < clout->num_pt ; ii++ )
       clout->mag[ii] /= nout ;
   } else if( nout == 0 ){
     ERROR_message("ACF estimation fails :-(") ;
     KILL_CLUSTER(clout) ; clout = NULL ;
   }

   if( clout != NULL ){
     float_quad acf_Epar ;
     acf_Epar = ACF_cluster_to_modelE( clout , dx,dy,dz ) ;
     fwhm_out = acf_Epar.d ;
     KILL_CLUSTER(clout) ; clout = NULL ;
   }

   RETURN(fwhm_out) ;
}

/*----------------------------------------------------------------------------*/
/* Function for Powell minimization for ACF modelE fit. */

static int    nar=0 ;
static float *aar=NULL , *rar=NULL ;

double ACF_modelE_costfunc( int npar , double *par )
{
   double aa=par[0] , bb=par[1]      , cc=par[2] ;
   double a1=1.0-aa , bq=0.5/(bb*bb) , ci=1.0/cc ;
   double sum=0.0 , fval ;
   int ii ;

   for( ii=0 ; ii < nar ; ii++ ){
     fval = aa * exp(-bq*rar[ii]*rar[ii]) + a1 * exp(-ci*rar[ii]) - aar[ii] ;
     sum += fval * fval ;
   }

   return sum ;
}

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

static MRI_IMAGE *ACF_im=NULL ;
MRI_IMAGE * ACF_get_1D(void){ return ACF_im ; }

/*----------------------------------------------------------------------------*/
/*! Convert ACF cluster output to function vs radius             [09 Nov 2015]
    Exponential model function vs. radius r is
      ACF(r) = a*exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)
    where a,b,c are the estimated parameters returned in the float_quad.
    With restrictions that   0 <= a <= 1 ; b > 0 ; c > 0
    The last parameter return (d) is the FWHM of the curve fit.
*//*--------------------------------------------------------------------------*/

float_quad ACF_cluster_to_modelE( MCW_cluster *acf, float dx, float dy, float dz )
{
   int nclu, nr, pp, qq, ss, dq ;
   float xx , yy , zz , rmax , dr , vp,rp , vs ;
   float apar , bpar , cpar , dpar , sig , *acar ;
   float_quad fqout = { -1.0f , -1.0f , -1.0f , -1.0f } ;
   double xpar[3] , xbot[3] , xtop[3] ;

ENTRY("ACF_cluster_to_modelE") ;

   if( ACF_im != NULL ){ mri_free(ACF_im) ; ACF_im = NULL ; }

                        if( acf == NULL ) RETURN(fqout) ;
   nclu = acf->num_pt ; if( nclu < 5    ) RETURN(fqout) ;

   if( dx <= 0.0f ) dx = 1.0f ;
   if( dy <= 0.0f ) dy = 1.0f ;
   if( dz <= 0.0f ) dz = 1.0f ;

   /* load rar and aar arrays with radius and ACF values */

#undef  CRAD
#define CRAD(a) ( xx = dx * abs(acf->i[a]) ,   \
                  yy = dy * abs(acf->j[a]) ,   \
                  zz = dz * abs(acf->k[a]) , sqrt(xx*xx+yy*yy+zz*zz) )

   aar = (float *)malloc(sizeof(float)*nclu) ;
   rar = (float *)malloc(sizeof(float)*nclu) ;

   for( pp=0 ; pp < nclu ; pp++ ){
     rar[pp] = CRAD(pp) ;
     aar[pp] = acf->mag[pp] ;
   }

   qsort_floatfloat( nclu , rar , aar ) ;  /* to ensure increasing rar */

   /* collapse data with very similar rar values */

   for( pp=0 ; pp < nclu-1 ; pp++ ){
     vp = aar[pp] ; rp = rar[pp]*1.01f ;
     for( qq=pp+1 ; qq < nclu && rar[qq] <= rp ; qq++ ) ; /*nada*/
     if( qq > pp+1 ){
       for( vs=0.0f,ss=pp ; ss < qq ; ss++ ) vs += aar[ss] ;
       aar[pp] = vs/(qq-pp) ;
       dq = qq-(pp+1) ;
       for( ss=qq ; ss < nclu ; ss++ ){
         aar[ss-dq] = aar[ss]; rar[ss-dq] = rar[ss];
       }
       nclu -= dq ;
     }
   }

   /* crude estimates for model parameters */

   nar = nclu ;

   for( pp=0 ; pp < nclu && aar[pp] > 0.5f ; pp++ ) ; /*nada*/
   if( pp == nclu ){
     ERROR_message("largest ACF found is %g -- too big for model fit!",aar[nclu-1]) ;
     free(aar) ; free(rar) ; aar = rar = NULL ; nar = 0 ;
     RETURN(fqout) ;
   }
   apar = 1.0 / logf(aar[pp]) ;
   bpar = sqrtf(-0.5*apar) * rar[pp] ;
   cpar = -rar[pp] * apar;
   apar = 0.5f ;

   /* nonlinear optimization of parameters */

   xpar[0] = apar  ; xpar[1] = bpar      ; xpar[2] = cpar      ;
   xbot[0] = 0.006 ; xbot[1] = 0.05*bpar ; xbot[2] = 0.05*cpar ;
   xtop[0] = 0.994 ; xtop[1] = 5.55*bpar ; xtop[2] = 5.55*cpar ;

#if 0
   pp = powell_newuoa_con( 3 , xpar , xbot , xtop ,
                           99 , 0.05 , 0.0005 , 999 , ACF_modelE_costfunc ) ;
#else
   pp = powell_newuoa_constrained( 3 , xpar , NULL , xbot , xtop ,
                                   666 , 44 , 9 ,
                                   0.05 , 0.0005 , 999 , ACF_modelE_costfunc ) ;
#endif

   if( pp < 0 ){
     ERROR_message("optimization of ACF model fit fails :-(") ;
     free(aar) ; free(rar) ; aar = rar = NULL ; nar = 0 ;
     RETURN(fqout) ;
   }

   ACF_im = mri_new( nar , 4 , MRI_float ) ;
   acar = MRI_FLOAT_PTR(ACF_im) ;
   apar = xpar[0] ; bpar = xpar[1] ; cpar = xpar[2] ;
   for( pp=0 ; pp < nar ; pp++ ){
     vs =         apar * expf( -0.5f*rar[pp]*rar[pp]/(bpar*bpar) )
         + (1.0f-apar) * expf( -rar[pp]/cpar ) ;
     acar[pp]       = rar[pp] ;  /* radius */
     acar[pp+nar]   = aar[pp] ;  /* empirical ACF */
     acar[pp+2*nar] = vs ;       /* model ACF */
   }

   { floatvec *fv ;
     MAKE_floatvec(fv,nar) ;
     for( pp=0 ; pp < nar ; pp++ ) fv->ar[pp] = acar[pp+2*nar] ;
     vs = interp_inverse_floatvec( fv , 0.5f ) ;
     for( pp=0 ; pp < nar ; pp++ ) fv->ar[pp] = rar[pp] ;
     dpar = 2.0*interp_floatvec( fv , vs ) ; sig = FWHM_TO_SIGMA(dpar) ;
     KILL_floatvec(fv) ;
   }

   for( pp=0 ; pp < nar ; pp++ ){  /* add in Gaussian ACF for fun */
     acar[pp+3*nar] = expf(-0.5f*rar[pp]*rar[pp]/(sig*sig)) ;
   }

   free(aar) ; free(rar) ; aar = rar = NULL ; nar = 0 ;

   fqout.a = xpar[0] ; fqout.b = xpar[1] ; fqout.c = xpar[2] ; fqout.d = dpar ;
   RETURN(fqout) ;
}
