/*
afni/src/3dLFCD.c
*/

/* 3dLFCD was created from 3dAutoTCorrelate by
   R. Cameron Craddock */

// Look for OpenMP macro
#ifdef USE_OMP
#include <omp.h>
#endif

// Include libraries
#include "mrilib.h"
#include <sys/mman.h>
#include <sys/types.h>

// Define constants
#define SPEARMAN 1
#define QUADRANT 2
#define PEARSON  3
#define ETA2     4

#define MAX_NUM_TAGS 32
#define MAX_TAG_LEN 256

#define FACES 0
#define FACES_EDGES 1
#define FACES_EDGES_CORNERS 2

#define X 0
#define Y 1
#define Z 2

#define faces_neighborhood_num 6
int faces_neighborhood[18] = {
     0,  0,  1,
     0,  0, -1,
     0,  1,  0,
     0, -1,  0,
     1,  0,  0,
    -1,  0,  0
};

/*   x  y  z */
#define faces_edges_neighborhood_num 18
const int faces_edges_neighborhood[54] = {
     0, -1, -1,
    -1,  0, -1,
     0,  0, -1,
     1,  0, -1,
     0,  1, -1,
    -1, -1,  0,
    -1,  0,  0,
    -1,  1,  0,
     0, -1,  0,
     0,  1,  0,
     1, -1,  0,
     1,  0,  0,
     1,  1,  0,
     0, -1,  1,
    -1,  0,  1,
     0,  0,  1,
     1,  0,  1,
     0,  1,  1
};

/*   x  y  z */
#define faces_edges_corners_neighborhood_num 26
const int faces_edges_corners_neighborhood[78] = {
    -1, -1, -1,
    -1,  0, -1,
    -1,  1, -1,
     0, -1, -1,
     0,  0, -1,
     0,  1, -1,
     1, -1, -1,
     1,  0, -1,
     1,  1, -1,
    -1, -1,  0,
    -1,  0,  0,
    -1,  1,  0,
     0, -1,  0,
     0,  1,  0,
     1, -1,  0,
     1,  0,  0,
     1,  1,  0,
    -1, -1,  1,
    -1,  0,  1,
    -1,  1,  1,
     0, -1,  1,
     0,  0,  1,
     0,  1,  1,
     1, -1,  1,
     1,  0,  1,
     1,  1,  1,
};

/* CC macro for updating mem stats */
#define INC_MEM_STATS( INC, TAG ) \
    { \
        if( MEM_PROF == 1 ) \
        { \
            int ndx = 0; \
            while( ndx < mem_num_tags ) \
            { \
                if( strncmp( mem_tags[ndx], TAG, MAX_TAG_LEN ) == 0 ) \
                { \
                    break; \
                } \
                ndx++; \
            } \
            if(( ndx >= mem_num_tags ) && (ndx < MAX_NUM_TAGS)) \
            { \
                /* adding a new tag */ \
                strncpy( mem_tags[ ndx ], TAG, (MAX_TAG_LEN-1) ); \
                mem_allocated[ ndx ] = 0; \
                mem_freed[ ndx ] = 0; \
                mem_num_tags++; \
            } \
            if( ndx < MAX_NUM_TAGS ) \
            { \
                mem_allocated[ ndx ] += (long)(INC); \
                if ((long)(INC) > 1024 ) WARNING_message("Incrementing memory for %s by %ldB\n", TAG, (INC)); \
            } \
            else WARNING_message("No room in mem profiler for %s\n", TAG ); \
        } \
        total_mem += (long)(INC); \
        running_mem += (long)(INC); \
        if (running_mem > peak_mem) peak_mem = running_mem; \
    }

#define DEC_MEM_STATS( DEC, TAG ) \
    { \
        if( MEM_PROF == 1 ) \
        { \
            int ndx = 0; \
            while( ndx < mem_num_tags ) \
            { \
                if( strncmp( mem_tags[ndx], TAG, MAX_TAG_LEN ) == 0 ) \
                { \
                    break; \
                } \
                else ndx++ ; \
            } \
            if(( ndx >= mem_num_tags ) && (ndx < MAX_NUM_TAGS)) \
            { \
                WARNING_message("Could not find tag %s in mem profiler\n", TAG ); \
            } \
            else \
            { \
                mem_freed[ ndx ] += (long)(DEC); \
                if ((long)(DEC) > 1024 ) INFO_message("Free %ldB of memory for %s\n", (DEC), TAG); \
            } \
        } \
        running_mem -= (long)(DEC); \
    }

#define PRINT_MEM_STATS( TAG ) \
        if ( MEM_STAT == 1 ) \
        { \
            INFO_message("\n======\n== Mem Stats (%s): Running %fB, Total %3.3fMB, Peak %3.3fMB\n", \
            TAG, \
            (double)(running_mem), \
            (double)(total_mem/(1024.0*1024.0)), \
            (double)(peak_mem/(1024.0*1024.0))); \
            INFO_message("\n======\n== Mem Stats (%s): Running %3.3fMB, Total %3.3fMB, Peak %3.3fMB\n", \
            TAG, \
            (double)(running_mem/(1024.0*1024.0)), \
            (double)(total_mem/(1024.0*1024.0)), \
            (double)(peak_mem/(1024.0*1024.0))); \
            if( MEM_PROF ==  1 ) \
            { \
                int ndx = 0; \
                INFO_message("== Memory Profile\n"); \
                for( ndx=0; ndx < mem_num_tags; ndx++ ) \
                { \
                    INFO_message("%s: %ld allocated %ld freed\n", mem_tags[ndx], \
                        mem_allocated[ndx], mem_freed[ndx] ); \
                } \
            } \
        }


/* CC define nodes for list that will be used
   for region growing */
typedef struct _list_node list_node;

// Define list node structure
struct _list_node
{
    long vox_vol_ndx;
    long ix;
    long jy;
    long kz;

    list_node* next;
};


/* free the list of list_nodes */
#define FREE_LIST( list ) \
{ \
    list_node *pptr; \
    list_node *hptr = list; \
 \
    while(hptr != NULL ) \
    { \
        /* increment node pointers */ \
        pptr = hptr; \
        hptr = hptr->next; \
 \
        /* delete the node */ \
        if(pptr != NULL) \
        {  \
            /* -- update running memory estimate to reflect memory allocation */  \
            DEC_MEM_STATS(( sizeof(list_node)), "list nodes"); \
            free(pptr); \
            pptr=NULL; \
        } \
    } \
    list = NULL; \
}

/* freeing all of the allocated mem on an error can get a little messy. instead
   we can use this macro to check what has been allocated and kill it. this of 
   course requires strict discipline for initializing all pointers to NULL and 
   resetting them to NULL when free'd. i should be able to handle that */
#define CHECK_AND_FREE_ALL_ALLOCATED_MEM \
{ \
    /* eliminate DSETS */ \
        if ( mset != NULL ) \
        { \
            DSET_unload(mset) ; \
            DSET_delete(mset) ; \
            mset = NULL ; \
        } \
\
        if ( xset != NULL ) \
        { \
            DSET_unload(xset) ; \
            DSET_delete(xset) ; \
            xset = NULL ; \
        } \
\
        if ( cset != NULL ) \
        { \
            DSET_unload(cset) ; \
            DSET_delete(cset) ; \
            cset = NULL ; \
        } \
\
     /* free the xvectim */ \
        if ( xvectim != NULL ) \
        { \
            VECTIM_destroy(xvectim) ; \
            xvectim = NULL ; \
        } \
\
    /* free allocated mems */ \
        if( mask != NULL ) \
        { \
            free(mask) ; \
            mask = NULL ; \
        } \
        \
        if( vol_ndx_to_mask_ndx != NULL ) \
        { \
            free(vol_ndx_to_mask_ndx) ; \
            vol_ndx_to_mask_ndx = NULL ; \
        } \
}

/* being good and cleaning up before erroring out can be a pain, and messy, lets
   make a variadic macro that does it for us. */
#define ERROR_EXIT_CC( ... ) \
    { \
        CHECK_AND_FREE_ALL_ALLOCATED_MEM; \
        ERROR_exit( __VA_ARGS__ ); \
    }




/*----------------------------------------------------------------*/
/**** Include these hesre for potential optimization for OpenMP ****/
/*----------------------------------------------------------------*/
/*! Pearson correlation of x[] and y[] (x and y are NOT modified.
    And we know ahead of time that the time series have 0 mean
    and L2 norm 1.
*//*--------------------------------------------------------------*/

float zm_THD_pearson_corr( int n, float *x , float *y ) /* inputs are */
{                                                       /* zero mean  */
   register float xy ; register int ii ;                /* and norm=1 */
   if( n%2 == 0 ){
     xy = 0.0f ;
     for( ii=0 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
   } else {
     xy = x[0]*y[0] ;
     for( ii=1 ; ii < n ; ii+=2 ) xy += x[ii]*y[ii] + x[ii+1]*y[ii+1] ;
   }
   return xy ;
}


/* */

/*----------------------------------------------------------------*/
/* General correlation calculation. */

#if 0
float my_THD_pearson_corr( int n, float *x , float *y )
{
   float xv,yv,xy , vv,ww , xm,ym ;
   register int ii ;

   xm = ym = 0.0f ;
   for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; }
   xm /= n ; ym /= n ;
   xv = yv = xy = 0.0f ;
   for( ii=0 ; ii < n ; ii++ ){
     vv = x[ii]-xm ; ww = y[ii]-ym ; xv += vv*vv ; yv += ww*ww ; xy += vv*ww ;
   }

   if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
   return xy/sqrtf(xv*yv) ;
}
#endif

/*----------------------------------------------------------------*/
/*! eta^2 (Cohen, NeuroImage 2008)              25 Jun 2010 [rickr]
 *
 *  eta^2 = 1 -  SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ]
 *               ------------------------------------
 *               SUM[ (a_i - M  )^2 + (b_i - M  )^2 ]
 *
 *  where  o  a_i and b_i are the vector elements
 *         o  m_i = (a_i + b_i)/2
 *         o  M = mean across both vectors
 -----------------------------------------------------------------*/

float my_THD_eta_squared( int n, float *x , float *y )
{
   float num , denom , gm , lm, vv, ww;
   register int ii ;

   gm = 0.0f ;
   for( ii=0 ; ii < n ; ii++ ){ gm += x[ii] + y[ii] ; }
   gm /= (2.0f*n) ;

   num = denom = 0.0f ;
   for( ii=0 ; ii < n ; ii++ ){
     lm = 0.5f * ( x[ii] + y[ii] ) ;
     vv = (x[ii]-lm); ww = (y[ii]-lm); num   += ( vv*vv + ww*ww );
     vv = (x[ii]-gm); ww = (y[ii]-gm); denom += ( vv*vv + ww*ww );
   }

   if( num < 0.0f || denom <= 0.0f || num >= denom ) return 0.0f ;
   return (1.0f - num/denom) ;
}

/*----------------------------------------------------------------------------*/
static void vstep_print(void)
{
   static int nn=0 ;
   static char xx[10] = "0123456789" ;
   fprintf(stderr , "%c" , xx[nn%10] ) ;
   if( nn%10 == 9) fprintf(stderr,",") ;
   nn++ ;
}

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

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

int main( int argc , char *argv[] )
{
    THD_3dim_dataset *xset = NULL;
    THD_3dim_dataset *cset = NULL;
    THD_3dim_dataset *mset = NULL ;
    int nopt=1 , method=PEARSON , do_autoclip=0 ;
    int nvox , nvals , ii, polort=1 ;
    char *prefix = "LFCD" ;
    byte *mask=NULL;
    int   nmask , abuc=1 ;
    char str[32] , *cpt = NULL;
    int *mask_ndx_to_vol_ndx = NULL;
    MRI_vectim *xvectim = NULL ;
    float (*corfun)(int,float *,float*) = NULL ;

    /* CC - we will have two subbricks: binary and weighted centrality */
    int nsubbriks = 2;
    int subbrik = 0;
    float * bodset;
    float * wodset;

    /* CC - added flags for thresholding correlations */
    double thresh = 0.0;

    double eps = 1e-9;

    /* CC - defines the neighborhood */
    int num_neighbors = faces_neighborhood_num;
    const int * neighborhood = faces_neighborhood;
  
    /* CC - variables to assist going back and forth between mask and volume indices */
    long * vol_ndx_to_mask_ndx = NULL;

    /* CC - variables for tracking memory usage stats */
    char mem_tags[MAX_NUM_TAGS][MAX_TAG_LEN];
    long mem_allocated[MAX_NUM_TAGS];
    long mem_freed[MAX_NUM_TAGS];
    long mem_num_tags = 0;
    long running_mem = 0;
    long peak_mem = 0;
    long total_mem = 0;
    int MEM_PROF = 0;
    int MEM_STAT = 0;
   /*----*/

   /* the number of zero variance voxels removed from the data */ 
   int rem_count = 0;

   AFNI_SETUP_OMP(0) ;  /* 24 Jun 2013 */

   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
      printf(
"Usage: 3dLFCD [options] dset\n"
"  Computes voxelwise local functional connectivity density as defined in:\n"
"      Tomasi, D and Volkow, PNAS, May 2010, 107 (21) 9885-9890; \n"
"        DOI: 10.1073/pnas.1001414107 \n\n"
"  The results are stored in a new 3D bucket dataset\n as floats to preserve\n"
"  their values. Local functional connectivity density (LFCD; as opposed to global\n"
"  functional connectivity density, see 3dDegreeCentrality), reflects\n"
"  the extent of the correlation of a voxel within its locally connected cluster.\n\n"
"  Conceptually the process involves: \n"
"      1. Calculating the correlation between voxel time series for\n"
"         every pair of voxels in the brain (as determined by masking)\n"
"      2. Applying a threshold to the resulting correlations to exclude\n"
"         those that might have arisen by chance\n"
"      3. Find the cluster of above-threshold voxels that are spatially\n"
"         connected to the target voxel.\n"
"      4. Count the number of voxels in the local cluster.\n"
"  Practically the algorithm is ordered differently to optimize for\n"
"  computational time and memory usage.\n\n"
"  The procedure described in the paper defines a voxels\n"
"  neighborhood to be the 6 voxels with which it shares a face.\n"
"  This definition can be changed to include edge and corner \n"
"  voxels using the -neighborhood flags below.\n\n"
"  LFCD is a localized variant of binarized degree centrality,\n"
"  the weighted alternative is calculated by changing step 4\n"
"  above to calculate the sum of the correlation coefficients\n"
"  between the seed region and the neighbors. 3dLFCD outputs\n"
"  both of these values (in separate briks), since they are\n"
"  so easy to calculate in tandem.\n\n"
"  You might prefer to calculate this on your data after\n"
"  spatial normalization, so that the range of values are\n"
"  consistent between datasets. Similarly the same brain mask\n"
"  should be used for all datasets that will be directly compared.\n\n"
"  The original paper used a correlation threshold = 0.6 and \n"
"  excluded all voxels with tSNR < 50. 3dLFCD does not discard\n"
"  voxels based on tSNR, this would need to be done beforehand.\n"
"\n"
"Options:\n"
"  -pearson  = Correlation is the normal Pearson (product moment)\n"
"               correlation coefficient [default].\n"
   #if 0
"  -spearman = Correlation is the Spearman (rank) correlation\n"
"               coefficient.\n"
"  -quadrant = Correlation is the quadrant correlation coefficient.\n"
   #else
"  -spearman AND -quadrant are disabled at this time :-(\n"
   #endif
"\n"
"  -thresh r = exclude correlations <= r from calculations\n"
"\n"
"  -faces    = define neighborhood to include face touching\n"
"              edges (default)\n"
"  -faces_edges = define neighborhood to include face and \n"
"                 edge touching voxels\n"
"  -faces_edges_corners = define neighborhood to include face, \n"
"                         edge, and corner touching voxels\n"
"\n"
"  -polort m = Remove polynomial trend of order 'm', for m=-1..3.\n"
"               [default is m=1; removal is by least squares].\n"
"               Using m=-1 means no detrending; this is only useful\n"
"               for data/information that has been pre-processed.\n"
"\n"
"  -autoclip = Clip off low-intensity regions in the dataset,\n"
"  -automask =  so that the correlation is only computed between\n"
"               high-intensity (presumably brain) voxels.  The\n"
"               mask is determined the same way that 3dAutomask works.\n"
"               This is done automatically if no mask is proveded.\n"
"\n"
"  -mask mmm = Mask to define 'in-brain' voxels. Reducing the number\n"
"               the number of voxels included in the calculation will\n"
"               significantly speedup the calculation. Consider using\n"
"               a mask to constrain the calculations to the grey matter\n"
"               rather than the whole brain. This is also preferable\n"
"               to using -autoclip or -automask.\n"
"\n"
"  -prefix p = Save output into dataset with prefix 'p', this file will\n"
"               contain bricks for both 'weighted' and 'binarized' lFCD\n"
"               [default prefix is 'LFCD'].\n"
"\n"
"Notes:\n"
" * The output dataset is a bucket type of floats.\n"
" * The program prints out an estimate of its memory used\n"
"    when it ends.  It also prints out a progress 'meter'\n"
"    to keep you pacified.\n"
"\n"
"-- RWCox - 31 Jan 2002 and 16 Jul 2010\n"
"-- Cameron Craddock - 13 Nov 2015 \n"
            ) ;
      PRINT_AFNI_OMP_USAGE("3dLFCD",NULL) ;
      PRINT_COMPILE_DATE ; exit(0) ;
   }

   mainENTRY("3dLFCD main"); machdep(); PRINT_VERSION("3dLFCD");
   AFNI_logger("3dLFCD",argc,argv);

   /*-- option processing --*/

   while( nopt < argc && argv[nopt][0] == '-' ){

        if( strcmp(argv[nopt],"-time") == 0 ){
            abuc = 0 ; nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-autoclip") == 0 ||
            strcmp(argv[nopt],"-automask") == 0   ){

            do_autoclip = 1 ; nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-mask") == 0 ){
            /* mset is opened here, but not loaded? */
            mset = THD_open_dataset(argv[++nopt]);
            CHECK_OPEN_ERROR(mset,argv[nopt]);
            nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-pearson") == 0 ){
            method = PEARSON ; nopt++ ; continue ;
        }

#if 0
        if( strcmp(argv[nopt],"-spearman") == 0 ){
            method = SPEARMAN ; nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-quadrant") == 0 ){
            method = QUADRANT ; nopt++ ; continue ;
        }
#endif

        if( strcmp(argv[nopt],"-eta2") == 0 ){
            method = ETA2 ; nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-prefix") == 0 ){
            prefix = strdup(argv[++nopt]) ;
            if( !THD_filename_ok(prefix) ){
                ERROR_EXIT_CC("Illegal value after -prefix!") ;
            }
            nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-thresh") == 0 ){
            double val = (double)strtod(argv[++nopt],&cpt) ;
            if( *cpt != '\0' || val >= 1.0 || val < 0.0 ){
                ERROR_EXIT_CC("Illegal value (%f) after -thresh!", val) ;
            }
            thresh = val ; nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-faces") == 0 ){
            num_neighbors = faces_neighborhood_num;
            neighborhood = faces_neighborhood;
            nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-faces_edges") == 0 ){
            num_neighbors = faces_edges_neighborhood_num;
            neighborhood = faces_edges_neighborhood;
            nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-faces_edges_neighborhood") == 0 ){
            num_neighbors = faces_edges_corners_neighborhood_num;
            neighborhood = faces_edges_corners_neighborhood;
            nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-polort") == 0 ){
            int val = (int)strtod(argv[++nopt],&cpt) ;
            if( *cpt != '\0' || val < -1 || val > 3 ){
                ERROR_EXIT_CC("Illegal value after -polort!") ;
            }
            polort = val ; nopt++ ; continue ;
        }

        if( strcmp(argv[nopt],"-mem_stat") == 0 ){
            MEM_STAT = 1 ; nopt++ ; continue ;
        }

        if( strncmp(argv[nopt],"-mem_profile",8) == 0 ){
            MEM_PROF = 1 ; nopt++ ; continue ;
        }

        ERROR_EXIT_CC("Illegal option: %s",argv[nopt]) ;
    }

    /*-- open dataset, check for legality --*/

    if( nopt >= argc ) ERROR_EXIT_CC("Need a dataset on command line!?") ;

    INFO_message("Calculating LFCD from %s using mask, r_threshold %lf and %d voxels in the neighborhood\n",
         argv[nopt], thresh, num_neighbors);

    xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);

    if( DSET_NVALS(xset) < 3 )
        ERROR_EXIT_CC("Input dataset %s does not have 3 or more sub-bricks!",
            argv[nopt]) ;
    DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;

    /*-- compute mask array, if desired --*/
    nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
    INC_MEM_STATS((nvox * nvals * sizeof(double)), "input dset");
    PRINT_MEM_STATS("inset");

    /* if a mask was specified make sure it is appropriate */
    if( mset ){

        if( DSET_NVOX(mset) != nvox )
            ERROR_EXIT_CC("Input and mask dataset differ in number of voxels!") ;
        mask  = THD_makemask(mset, 0, 1.0, 0.0) ;

        /* update running memory statistics to reflect loading the image */
        INC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );
        PRINT_MEM_STATS( "mset load" );

        /* iupdate statistics to reflect creating mask array */
        nmask = THD_countmask( nvox , mask ) ;
        INC_MEM_STATS( nmask * sizeof(byte), "mask array" );
        PRINT_MEM_STATS( "mask" );

        INFO_message("%d voxels in -mask dataset",nmask) ;
        if( nmask < 2 ) ERROR_EXIT_CC("Only %d voxels in -mask, exiting...",nmask);

        /* update running memory statistics to reflect loading the image */
        DEC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );
        PRINT_MEM_STATS( "mset unload" );

        /* free all memory associated with the mask datast */
        DSET_unload(mset) ;
        DSET_delete(mset) ;
        mset = NULL ;
    } 
    /* if automasking is requested, handle that now */
    else {
        mask  = THD_automask( xset ) ;
        nmask = THD_countmask( nvox , mask ) ;
        INC_MEM_STATS( nmask * sizeof(byte), "mask array" );
        PRINT_MEM_STATS( "mask" );
        INFO_message("No mask provided, created one using AFNI's automask procedure, %d voxels survive",nmask) ;
        if( nmask < 2 ) ERROR_EXIT_CC("Only %d voxels in -automask!",nmask);
    }
    
    /*-- create vectim from input dataset --*/
    INFO_message("vectim-izing input dataset") ;

    /*-- CC added in mask to reduce the size of xvectim -- */
    xvectim = THD_dset_to_vectim( xset , mask , 0 ) ;
    if( xvectim == NULL ) ERROR_EXIT_CC("Can't create vectim?!") ;

    /*-- CC update our memory stats to reflect vectim -- */
    INC_MEM_STATS((xvectim->nvec*sizeof(int)) +
                    ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
                    sizeof(MRI_vectim), "vectim");
    PRINT_MEM_STATS( "vectim" );

    /* iterate through and remove all voxels that have zero
       variance */
    for (ii=0; ii<xvectim->nvec; ii++)
    {
        double sum = 0.0;
        double sum_sq = 0.0;
        int    tt;

        float* xsar = VECTIM_PTR(xvectim,ii);

        for(tt=0; tt<xvectim->nvals; tt++)
        {
            sum += xsar[tt];
            sum_sq += xsar[tt] * xsar[tt];
        }

        if((sum_sq-sum*sum/(double)nvals)/(double)(nvals-1) < eps)
        {
            mask[xvectim->ivec[ii]] = 0;
            rem_count++;
        }
    }

    /* update running memory statistics to reflect freeing the vectim */
    DEC_MEM_STATS(((xvectim->nvec*sizeof(int)) +
                       ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
                       sizeof(MRI_vectim)), "vectim");

    /* toss some trash */
    VECTIM_destroy(xvectim) ;

    /* make another xvectim with the updated mask */
    xvectim = THD_dset_to_vectim( xset , mask , 0 ) ;
    if( xvectim == NULL ) ERROR_EXIT_CC("Can't create vectim?!") ;

    /*-- CC update our memory stats to reflect vectim -- */
    INC_MEM_STATS((xvectim->nvec*sizeof(int)) +
                    ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
                    sizeof(MRI_vectim), "vectim");
    PRINT_MEM_STATS( "vectim" ); 

    INFO_message("!! %d voxels remain after removing those with no variance (%d)\n", 
        xvectim->nvec, rem_count);


    /* unload DSET to reduce the amount of memory used */
    DSET_unload(xset);
    DEC_MEM_STATS((nvox * nvals * sizeof(double)), "input dset");

    /**  For the case of Pearson correlation, we make sure the  **/
    /**  data time series have their mean removed (polort >= 0) **/
    /**  and are normalized, so that correlation = dot product, **/
    /**  and we can use function zm_THD_pearson_corr for speed. **/

    switch( method ){
        default:
        case PEARSON: corfun = zm_THD_pearson_corr ; break ;
        case ETA2:    corfun = my_THD_eta_squared  ; break ;
    }

    if( method == ETA2 && polort >= 0 )
        WARNING_message("Polort for -eta2 should probably be -1...");
    

   /*--- CC the vectim contains a mapping between voxel index and mask index, 
         tap into that here to avoid duplicating memory usage, also create a
         mapping that goes the other way ---*/

    if( mask != NULL )
    {
        /* tap into the xvectim mapping */
        mask_ndx_to_vol_ndx = xvectim->ivec;

        /* create a mapping that goes the opposite way */
        if(( vol_ndx_to_mask_ndx = (long*)calloc( nvox, sizeof(long) )) == NULL)
        {
            ERROR_EXIT_CC( 
                "Could not allocate %d byte array for mask index mapping\n",
                nmask*sizeof(long));
        }

        /* --- CC account for the size of the mask */
        INC_MEM_STATS( nvox*sizeof(long), "ndx mask array" );

        /* --- create the reverse mapping */
        for ( ii = 0; ii < xvectim->nvec; ii++)
        {
            vol_ndx_to_mask_ndx[ mask_ndx_to_vol_ndx[ ii ] ] = ii;
        }

        /* --- CC free the mask */
        DEC_MEM_STATS( nmask*sizeof(byte), "mask array" );
        free(mask); mask=NULL;
        PRINT_MEM_STATS( "mask unload" );
    }

    /* -- CC configure detrending --*/
    if( polort < 0 && method == PEARSON )
    {
        polort = 0; 
        WARNING_message("Pearson correlation always uses polort >= 0");
    }
    if( polort >= 0 )
    {
        for( ii=0 ; ii < xvectim->nvec ; ii++ )
        {  
            /* remove polynomial trend */
            DETREND_polort(polort,nvals,VECTIM_PTR(xvectim,ii)) ;
        }
    }


    /* -- this procedure does not change time series that 
          have zero variance -- */
    if( method == PEARSON ) THD_vectim_normalize(xvectim) ;  /* L2 norm = 1 */

    /*-- create output dataset --*/
    cset = EDIT_empty_copy( xset ) ;

    /*-- configure the output dataset */
    if( abuc )
    {
        EDIT_dset_items( 
            cset ,
            ADN_prefix    , prefix         ,
            ADN_nvals     , nsubbriks      , /* 2 subbricks */
            ADN_ntt       , 0              , /* no time axis */
            ADN_type      , HEAD_ANAT_TYPE ,
            ADN_func_type , ANAT_BUCK_TYPE ,
            ADN_datum_all , MRI_float      ,
            ADN_none ) ;
    } 
    else 
    {
        EDIT_dset_items( cset ,
            ADN_prefix    , prefix         ,
            ADN_nvals     , nsubbriks      , /* 2 subbricks */
            ADN_ntt       , nsubbriks      ,  /* num times */
            ADN_ttdel     , 1.0            ,  /* fake TR */
            ADN_nsl       , 0              ,  /* no slice offsets */
            ADN_type      , HEAD_ANAT_TYPE ,
            ADN_func_type , ANAT_EPI_TYPE  ,
            ADN_datum_all , MRI_float      ,
            ADN_none ) ;
    }

    /* add history information to the hearder */
    tross_Make_History( "3dLFCD" , argc,argv , cset ) ;

    INFO_message("creating output dataset in memory") ;

    /* -- Configure the subbriks: Binary LFCD */
    subbrik = 0;
    EDIT_BRICK_TO_NOSTAT(cset,subbrik) ;                     /* stat params  */
    /* CC this sets the subbrik scaling factor, which we will probably want
      to do again after we calculate the voxel values */
    EDIT_BRICK_FACTOR(cset,subbrik,1.0) ;                 /* scale factor */

    sprintf(str,"Binary LFCD") ;

    EDIT_BRICK_LABEL(cset,subbrik,str) ;
    EDIT_substitute_brick(cset,subbrik,MRI_float,NULL) ;   /* make array   */

    /* copy measure data into the subbrik */
    bodset = DSET_ARRAY(cset,subbrik);
 
    /* -- Configure the subbriks: Weighted LFCD */
    subbrik = 1;
    EDIT_BRICK_TO_NOSTAT(cset,subbrik) ;                     /* stat params  */
    /* CC this sets the subbrik scaling factor, which we will probably want
       to do again after we calculate the voxel values */
    EDIT_BRICK_FACTOR(cset,subbrik,1.0) ;                 /* scale factor */

    sprintf(str,"Weighted LFCD") ;

    EDIT_BRICK_LABEL(cset,subbrik,str) ;
    EDIT_substitute_brick(cset,subbrik,MRI_float,NULL) ;   /* make array   */

    /* copy measure data into the subbrik */
    wodset = DSET_ARRAY(cset,subbrik);

    /* increment memory stats */
    INC_MEM_STATS( (DSET_NVOX(cset)*DSET_NVALS(cset)*sizeof(float)),
        "output dset");
    PRINT_MEM_STATS( "outset" );    

    /*-- tell the user what we are about to do --*/
    INFO_message( "Calculating LFCD with threshold = %f.\n", thresh);

    /*---------- loop over mask voxels, correlate ----------*/
    AFNI_OMP_START ;
#pragma omp parallel if( nmask > 999 )
    {
        int dx = 0;
        int dy = 0;
        int dz = 0;
        int ix = 0;
        int jy = 0;
        int kz = 0;
        int target_mask_ndx = 0;
        int lout,ithr,nthr,vstep,vii ;
        int neighbor_index;
        float *xsar , *ysar ;
        list_node* current_node = NULL ;
        list_node* new_node = NULL ;
        list_node* recycled_nodes = NULL ;
        list_node* boundary_list = NULL ;
        long* seen_voxels;
        double car = 0.0 ; 

        /*-- get information about who we are --*/
#ifdef USE_OMP
        ithr = omp_get_thread_num() ;
        nthr = omp_get_num_threads() ;
        if( ithr == 0 ) INFO_message("%d OpenMP threads started",nthr) ;
#else
        ithr = 0 ; nthr = 1 ;
#endif

        /*-- For the progress tracker, we want to print out 50 numbers,
            figure out a number of loop iterations that will make this easy */
        vstep = (int)( nmask / (nthr*50.0f) + 0.901f ) ; vii = 0 ;
        if((MEM_STAT==0) && (ithr == 0 )) fprintf(stderr,"Looping:") ;

        /* allocate a vector to track previously seen voxels */
        if((seen_voxels = (long*)calloc(xvectim->nvec,sizeof(long)))==NULL)
        {
            ERROR_EXIT_CC( "Could not allocate memory for seen_voxels!\n");
        }

#pragma omp for schedule(static, 1)
        for( lout=0 ; lout < xvectim->nvec ; lout++ )
        {  /*----- outer voxel loop -----*/

            if( ithr == 0 && vstep > 2 )
            { 
                vii++ ; 
                if( vii%vstep == vstep/2 && MEM_STAT == 0 ) vstep_print(); 
            }

            /* get ref time series from this voxel */
            xsar = VECTIM_PTR(xvectim,lout) ;

            /* initialize the boundary with the target voxel */
            if( recycled_nodes == NULL)
            {
                /* this looks like it is redundant, but I want the first if 
                   statement to run in the critical section and the second
                   if statement to run after it */
#pragma omp critical(mem_alloc)
                if(( new_node = (list_node*)malloc(sizeof(list_node)) ) != NULL )
                {
                    INC_MEM_STATS((sizeof(list_node)), "list nodes");
                }
                if( new_node == NULL )
                {
                    WARNING_message( "Could not allocate list node\n");
                    continue;
                }
            }
            else
            {
                new_node = recycled_nodes;
                recycled_nodes = recycled_nodes->next;
            }

            /* determine the full ndx for the seed target */
            new_node->vox_vol_ndx = mask_ndx_to_vol_ndx[ lout ];

            /* add source, dest, correlation to 1D file */
            new_node->ix = DSET_index_to_ix(cset, new_node->vox_vol_ndx) ;
            new_node->jy = DSET_index_to_jy(cset, new_node->vox_vol_ndx) ;
            new_node->kz = DSET_index_to_kz(cset, new_node->vox_vol_ndx) ;

            /* push the new node onto the boundary list,
               if the boundary list is empty, there is a problem */
            if (boundary_list != NULL)
            {
                WARNING_message("Boundary list not empty!\n");
            }
            new_node->next = boundary_list;
            boundary_list = new_node;
            
            /* reset the seen_voxels map */
            memset(seen_voxels, 0, xvectim->nvec*sizeof(long));

            /* iterate over the boundary until it is empty */
            while( boundary_list != NULL )
            {

                /* pop a node off of the list */
                current_node = boundary_list;
                boundary_list = boundary_list->next;

                /* iterate through a box around the current voxel */
                for ( neighbor_index = 0; neighbor_index < num_neighbors; neighbor_index++ )
                {
                    dx = neighborhood[3*neighbor_index+X];
                    dy = neighborhood[3*neighbor_index+Y];
                    dz = neighborhood[3*neighbor_index+Z];

                    ix = ( current_node->ix + (dx-1) );
                    /* make sure that we are in bounds */
                    if (( ix < 0 ) || ( ix > DSET_NX(xset) )) 
                    {
                        continue;
                    }

                    jy = ( current_node->jy + (dy-1) );
                    /* make sure that we are in bounds */                        
                    if (( jy < 0 ) || ( jy > DSET_NY(xset) ))
                    {
                        continue;
                    }

                    kz = ( current_node->kz + (dz-1) );
                    /* make sure that we are in bounds */
                    if (( kz < 0 ) || (kz > DSET_NZ(xset))) 
                    {
                        continue;
                    }

                    /* get the index of this voxel */
                    target_mask_ndx = 
                        vol_ndx_to_mask_ndx[ DSET_ixyz_to_index(xset,ix,jy,kz) ];

                    /* if the voxel is in the mask, and hasn't been
                        seen before, evaluate it for inclusion in the
                        boundary */
                    if(( target_mask_ndx != 0 ) && ( seen_voxels[ target_mask_ndx ] != 1 ))
                    {

                        /* indicate that we have seen the voxel */
                        seen_voxels[ target_mask_ndx ] = 1;

                        /* extract the time series */
                        ysar = VECTIM_PTR(xvectim,target_mask_ndx) ;

                        /* calculate the correlation */
                        car = (double)(corfun(nvals,xsar,ysar)) ;

                        /* if correlation is above threshold, add
                            it to the LFCD stats and to the boundary */
                        if ( car > thresh )
                        {

                            if( recycled_nodes == NULL)
                            {
        /* this looks like it is redundant, but I want the first if 
            statement to run in the critical section and the second
            if statement to run after it */
#pragma omp critical(mem_alloc)
                                if(( new_node = (list_node*)malloc(sizeof(list_node)) ) != NULL )
                                {
                                    INC_MEM_STATS((sizeof(list_node)), "list nodes");
                                }
                                if( new_node == NULL )
                                {
                                    WARNING_message( "Could not allocate list node\n");
                                    continue;
                                }
                            }
                            else
                            {
                                new_node = recycled_nodes;
                                recycled_nodes = recycled_nodes->next;
                            }

                            /* determine the full ndx for the seed target */
                            new_node->vox_vol_ndx = mask_ndx_to_vol_ndx[ target_mask_ndx ]; 

                            /* add source, dest, correlation to 1D file */
                            new_node->ix = ix ;
                            new_node->jy = jy ;
                            new_node->kz = kz ;

                            /* add the node to the boundary */
                            new_node->next = boundary_list;
                            boundary_list = new_node;

                            /* now increment the LFCD measures, this is
                                done in a critical section */
#pragma omp critical(dataupdate)
                            {
                                wodset[ mask_ndx_to_vol_ndx[ lout ]] += car;
                                bodset[ mask_ndx_to_vol_ndx[ lout ]] += 1;
                            }                                  
                        } /* if car > thresh */
                    } /* if vox is in mask and hasn't been seen */
                } /* for neighbor_index */

                /* we have finished processing this node, so recycle it */
                current_node->next = recycled_nodes;
                recycled_nodes = current_node;

            } /* while( boundary_list != NULL ) */
        } /* for lout */
 
        /* clean up the memory that is local to this thread */
        if (seen_voxels != NULL) free( seen_voxels );

        /* clean out the recycled list */
#pragma omp critical(mem_alloc)
        FREE_LIST( recycled_nodes )
        PRINT_MEM_STATS( "Free recycled nodes" );

    } /* end OpenMP */
    AFNI_OMP_END ;

    /* update the user so that they know what we are up to */
    INFO_message ("AFNI_OMP finished\n");
 
    INFO_message("Done..\n") ;

    if( vol_ndx_to_mask_ndx != NULL )
    {
        DEC_MEM_STATS(nvox*sizeof(long), "ndx mask array");
        free(vol_ndx_to_mask_ndx);
        vol_ndx_to_mask_ndx = NULL;
        PRINT_MEM_STATS( "Free ndx mask array" );
    }

    /* update running memory statistics to reflect freeing the vectim */
    DEC_MEM_STATS(((xvectim->nvec*sizeof(int)) +
                       ((xvectim->nvec)*(xvectim->nvals))*sizeof(float) +
                       sizeof(MRI_vectim)), "vectim");

    /* toss some trash */
    VECTIM_destroy(xvectim) ;


    DSET_delete(xset) ;

    PRINT_MEM_STATS( "Vectim and Input dset Unload" );

    /* finito */
    INFO_message("Writing output dataset to disk [%s bytes]",
                commaized_integer_string(cset->dblk->total_bytes)) ;

    /* write the dataset */
    DSET_write(cset) ;
    WROTE_DSET(cset) ;

    /* increment our memory stats, since we are relying on the header for this
       information, we update the stats before actually freeing the memory */
    DEC_MEM_STATS( (DSET_NVOX(cset)*DSET_NVALS(cset)*sizeof(float)), "output dset");

    /* free up the output dataset memory */
    DSET_unload(cset) ;
    DSET_delete(cset) ;

    PRINT_MEM_STATS( "Unload Output Dset" );

    /* force a print */
    MEM_STAT = 1;
    PRINT_MEM_STATS( "Fin" );

    exit(0) ;
}
