#ifdef USE_OMP #include #endif #include "mrilib.h" #include #include #define SPEARMAN 1 #define QUADRANT 2 #define PEARSON 3 #define ETA2 4 #define ALLOW_MMAP #ifndef ALLOW_MMAP # define do_mmap 0 #else static int do_mmap = 0 ; #endif /*----------------------------------------------------------------*/ /**** Include these here 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++ ; } /*-----------------------------------------------------------------------------*/ #ifdef ALLOW_MMAP static char *cbrik = NULL ; static int64_t ncbrik = 0 ; static char *cbrik_name = NULL ; #undef MFLAG #if defined(MAP_NOCACHE) # define MFLAG (MAP_SHARED | MAP_NOCACHE) #elif defined(MAP_NORESERVE) # define MFLAG (MAP_SHARED | MAP_NORESERVE) #else # define MFLAG MAP_SHARED #endif #include void AC_sigfunc(int sig) /** signal handler for fatal errors **/ { char *sname ; static volatile int fff=0 ; if( fff ) _exit(1) ; else fff = 1 ; switch(sig){ default: sname = "unknown" ; break ; case SIGINT: sname = "SIGINT" ; break ; case SIGPIPE: sname = "SIGPIPE" ; break ; case SIGSEGV: sname = "SIGSEGV" ; break ; case SIGBUS: sname = "SIGBUS" ; break ; case SIGTERM: sname = "SIGTERM" ; break ; } fprintf(stderr,"\n** 3dAutoTcorrelate: Fatal Signal %d (%s) received\n",sig,sname) ; if( cbrik != NULL ){ fprintf(stderr,"** Un-mmap-ing %s file (but not deleting it)\n",cbrik_name) ; munmap(cbrik,ncbrik) ; cbrik = NULL ; } exit(1) ; } #endif /* ALLOW_MMAP */ /*----------------------------------------------------------------------------*/ int main( int argc , char *argv[] ) { THD_3dim_dataset *xset , *cset, *mset=NULL, *msetinner=NULL ; int nopt=1 , method=PEARSON , do_autoclip=0 ; int nvox , nvals , ii,kout , polort=1 , ix,jy,kz ; char *prefix = "ATcorr" ; byte *mask=NULL, *maskinner=NULL; int nmask , nmaskinner , abuc=1 ; int all_source=0; /* output all source voxels 25 Jun 2010 [rickr] */ char str[32] , *cpt ; int *imap ; MRI_vectim *xvectim ; float (*corfun)(int,float *,float*) = NULL ; FILE *fout1D=NULL; /*----*/ AFNI_SETUP_OMP(0) ; /* 24 Jun 2013 */ if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ printf( "Usage: 3dAutoTcorrelate [options] dset\n" "Computes the correlation coefficient between the time series of each\n" "pair of voxels in the input dataset, and stores the output into a\n" "new anatomical bucket dataset [scaled to shorts to save memory space].\n" "\n" "*** Also see program 3dTcorrMap ***\n" "\n" "Options:\n" " -pearson = Correlation is the normal Pearson (product moment)\n" " correlation coefficient [default].\n" " -eta2 = Output is eta^2 measure from Cohen et al., NeuroImage, 2008:\n" " http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2705206/\n" " http://dx.doi.org/10.1016/j.neuroimage.2008.01.066\n" " ** '-eta2' is intended to be used to measure the similarity\n" " between 2 correlation maps; therefore, this option is\n" " to be used in a second stage analysis, where the input\n" " dataset is the output of running 3dAutoTcorrelate with\n" " the '-pearson' option -- the voxel 'time series' from\n" " that first stage run is the correlation map of that\n" " voxel with all other voxels.\n" " ** '-polort -1' is recommended with this option!\n" " ** Odds are you do not want use this option if the dataset\n" " on which eta^2 is to be computed was generated with\n" " options -mask_only_targets or -mask_source.\n" " In this program, the eta^2 is computed between pseudo-\n" " timeseries (the 4th dimension of the dataset).\n" " If you want to compute eta^2 between sub-bricks then use\n" " 3ddot -eta2 instead.\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" " -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" "\n" " -mask mmm = Mask of both 'source' and 'target' voxels.\n" " ** Restricts computations to those in the mask. Output\n" " volumes are restricted to masked voxels. Also, only\n" " masked voxels will have non-zero output.\n" " ** A dataset with 1000 voxels would lead to output of\n" " 1000 sub-bricks. With a '-mask' of 50 voxels, the\n" " output dataset have 50 sub-bricks, where the 950\n" " unmasked voxels would be all zero in all 50 sub-bricks\n" " (unless option '-mask_only_targets' is also used).\n" " ** The mask is encoded in the output dataset header in the\n" " attribute named 'AFNI_AUTOTCORR_MASK' (cf. 3dMaskToASCII).\n" "\n" " -mask_only_targets = Provide output for all voxels.\n" " ** Used with '-mask': every voxel is correlated with each\n" " of the mask voxels. In the example above, there would\n" " be 50 output sub-bricks; the n-th output sub-brick\n" " would contain the correlations of the n-th voxel in\n" " the mask with ALL 1000 voxels in the dataset (rather\n" " than with just the 50 voxels in the mask).\n" "\n" " -mask_source sss = Provide output for voxels only in mask sss\n" " ** For each seed in mask mm, compute correlations only with \n" " non-zero voxels in sss. If you have 250 non-zero voxels \n" " in sss, then the output will still have 50 sub-bricks, but\n" " each n-th sub-brick will have non-zero values at the 250\n" " non-zero voxels in sss\n" " Do not use this option along with -mask_only_targets\n" "\n" " -prefix p = Save output into dataset with prefix 'p'\n" " [default prefix is 'ATcorr'].\n" " -out1D FILE.1D = Save output in a text file formatted thusly:\n" " Row 1 contains the 1D indices of non zero voxels in the \n" " mask from option -mask.\n" " Column 1 contains the 1D indices of non zero voxels in the\n" " mask from option -mask_source\n" " The rest of the matrix contains the correlation/eta2 \n" " values. Each column k corresponds to sub-brick k in \n" " the output volume p.\n" " To see 1D indices in AFNI, right click on the top left\n" " corner of the AFNI controller - where coordinates are\n" " shown - and chose voxel indices.\n" " A 1D index (ijk) is computed from the 3D (i,j,k) indices:\n" " ijk = i + j*Ni + k*Ni*Nj , with Ni and Nj being the\n" " number of voxels in the slice orientation and given by:\n" " 3dinfo -ni -nj YOUR_VOLUME_HERE\n" " This option can only be used in conjunction with \n" " options -mask and -mask_source. Otherwise it makes little\n" " sense to write a potentially enormous text file.\n" "\n" " -time = Mark output as a 3D+time dataset instead of an anat bucket.\n" #ifdef ALLOW_MMAP "\n" " -mmap = Write .BRIK results to disk directly using Unix mmap().\n" " This trick can speed the program up when the amount\n" " of memory required to hold the output is very large.\n" " ** In many case, the amount of time needed to write\n" " the results to disk is longer than the CPU time.\n" " This option can shorten the disk write time.\n" " ** If the program crashes, you'll have to manually\n" " remove the .BRIK file, which will have been created\n" " before the loop over voxels and written into during\n" " that loop, rather than being written all at once\n" " at the end of the analysis, as is usually the case.\n" " ** If the amount of memory needed is bigger than the\n" " RAM on your system, this program will be very slow\n" " with or without '-mmap'.\n" " ** This option won't work with NIfTI-1 (.nii) output!\n" #else "\n" " -mmap is disabled at this time :-(\n" #endif "\n" "Example: correlate every voxel in mask_in+tlrc with only those voxels in\n" " mask_out+tlrc (the rest of each volume is zero, for speed).\n" " Assume detrending was already done along with other pre-processing.\n" " The output will have one volume per masked voxel in mask_in+tlrc.\n" " Volumes will be labeled by the ijk index triples of mask_in+tlrc.\n" "\n" " 3dAutoTcorrelate -mask_source mask_out+tlrc -mask mask_in+tlrc \\\n" " -polort -1 -prefix test_corr clean_epi+tlrc\n" "\n" "Notes:\n" " * The output dataset is anatomical bucket type of shorts\n" " (unless '-time' is used).\n" " * Values are scaled so that a correlation (or eta-squared)\n" " of 1 corresponds to a value of 10000.\n" " * The output file might be gigantic and you might run out\n" " of memory running this program. Use at your own risk!\n" " ++ If you get an error message like\n" " *** malloc error for dataset sub-brick\n" " this means that the program ran out of memory when making\n" " the output dataset.\n" #ifdef ALLOW_MMAP " ++ If this happens, you can try to use the '-mmap' option,\n" " and if you are lucky, the program may actually run.\n" #endif " * The program prints out an estimate of its memory usage\n" " when it starts. It also prints out a progress 'meter'\n" " to keep you pacified.\n" " * This is a quick hack for Peter Bandettini. Now pay up.\n" " * OpenMP-ized for Hang Joon Jo. Where's my baem-sul?\n" "\n" "-- RWCox - 31 Jan 2002 and 16 Jul 2010\n" ) ; PRINT_AFNI_OMP_USAGE("3dAutoTcorrelate",NULL) ; PRINT_COMPILE_DATE ; exit(0) ; } mainENTRY("3dAutoTcorrelate main"); machdep(); PRINT_VERSION("3dAutoTcorrelate"); AFNI_logger("3dAutoTcorrelate",argc,argv); /*-- option processing --*/ while( nopt < argc && argv[nopt][0] == '-' ){ #ifdef ALLOW_MMAP if( strcmp(argv[nopt],"-mmap") == 0 ){ /* Jul 2010 */ do_mmap = 1 ; nopt++ ; continue ; } #endif 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 = THD_open_dataset(argv[++nopt]); CHECK_OPEN_ERROR(mset,argv[nopt]); nopt++ ; continue ; } if( strcmp(argv[nopt],"-mask_source") == 0 ){ msetinner = THD_open_dataset(argv[++nopt]); CHECK_OPEN_ERROR(msetinner,argv[nopt]); nopt++ ; continue ; } if( strcmp(argv[nopt],"-out1D") == 0 ){ if (!(fout1D = fopen(argv[++nopt], "w"))) { ERROR_message("Failed to open %s for writing", argv[nopt]); exit(1); } nopt++ ; continue ; } if( strcmp(argv[nopt],"-mask_only_targets") == 0 ){ all_source = 1 ; 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("Illegal value after -prefix!") ; } nopt++ ; continue ; } if( strcmp(argv[nopt],"-polort") == 0 ){ int val = (int)strtod(argv[++nopt],&cpt) ; if( *cpt != '\0' || val < -1 || val > 3 ){ ERROR_exit("Illegal value after -polort!") ; } polort = val ; nopt++ ; continue ; } ERROR_exit("Illegal option: %s",argv[nopt]) ; } #ifdef ALLOW_MMAP cpt = strstr(prefix,".nii") ; /* Jul 2010 */ if( do_mmap && cpt != NULL ){ *cpt = '\0' ; if( *prefix == '\0' ) prefix = "ATcorr" ; WARNING_message("-mmap ==> can't use NIfTI-1 ==> output prefix = %s",prefix) ; } #endif /*-- open dataset, check for legality --*/ if( nopt >= argc ) ERROR_exit("Need a dataset on command line!?") ; if (fout1D && (!mset || !msetinner)) { fclose(fout1D); ERROR_message("Must use -mask and -mask_source with -out1D"); exit(1); } xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]); if( DSET_NVALS(xset) < 3 ) ERROR_exit("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) ; if( mset ){ if( DSET_NVOX(mset) != nvox ) ERROR_exit("Input and mask dataset differ in number of voxels!") ; mask = THD_makemask(mset, 0, 1.0, 0.0) ; nmask = THD_countmask( nvox , mask ) ; INFO_message("%d voxels in -mask dataset",nmask) ; if( nmask < 2 ) ERROR_exit("Only %d voxels in -mask, exiting...",nmask); DSET_unload(mset) ; } else if( do_autoclip ){ mask = THD_automask( xset ) ; nmask = THD_countmask( nvox , mask ) ; INFO_message("%d voxels survive -autoclip",nmask) ; if( nmask < 2 ) ERROR_exit("Only %d voxels in -automask!",nmask); } else { nmask = nvox ; INFO_message("computing for all %d voxels",nmask) ; } if (msetinner) { if( DSET_NVOX(msetinner) != nvox ) ERROR_exit("Input and mask_source dataset differ in number of voxels!"); maskinner = THD_makemask(msetinner, 0, 1.0, 0.0) ; nmaskinner = THD_countmask( nvox , maskinner ) ; INFO_message("%d voxels in -mask dataset",nmaskinner) ; if( nmaskinner < 2 ) ERROR_exit("Only %d voxels in -mask, exiting...",nmaskinner); DSET_unload(msetinner) ; DSET_delete(msetinner); msetinner = NULL; } else { if (mask) { maskinner = mask; nmaskinner = -nvox; /* a flag to be sure we remember this deed */ } } if( method == ETA2 && polort >= 0 ) WARNING_message("Polort for -eta2 should probably be -1..."); /*-- create vectim from input dataset --*/ /** 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 ; } INFO_message("vectim-izing input dataset") ; xvectim = THD_dset_to_vectim( xset , NULL , 0 ) ; if( xvectim == NULL ) ERROR_exit("Can't create vectim?!") ; DSET_unload(xset) ; if( polort < 0 && method == PEARSON ){ polort = 0; WARNING_message("Pearson correlation always uses polort >= 0"); } if( polort >= 0 ){ for( ii=0 ; ii < nvox ; ii++ ){ /* remove polynomial trend */ DETREND_polort(polort,nvals,VECTIM_PTR(xvectim,ii)) ; } } if( method == PEARSON ) THD_vectim_normalize(xvectim) ; /* L2 norm = 1 */ /*-- create output dataset --*/ cset = EDIT_empty_copy( xset ) ; if( abuc ){ EDIT_dset_items( cset , ADN_prefix , prefix , ADN_nvals , nmask , ADN_ntt , 0 , /* no time axis */ ADN_type , HEAD_ANAT_TYPE , ADN_func_type , ANAT_BUCK_TYPE , ADN_datum_all , MRI_short , ADN_none ) ; } else { EDIT_dset_items( cset , ADN_prefix , prefix , ADN_nvals , nmask , ADN_ntt , nmask , /* 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_short , ADN_none ) ; } if( THD_deathcon() && THD_is_file(DSET_HEADNAME(cset)) ) ERROR_exit("Output dataset %s already exists!",DSET_HEADNAME(cset)) ; { double nb = (double)(xset->dblk->total_bytes) +(double)(cset->dblk->total_bytes) ; nb /= 1000000 ; INFO_message( "Memory required = %.1f Mbytes for %d output sub-bricks",nb,nmask); if( nb > 2000.0 && (sizeof(void *) < 8 || sizeof(size_t) < 8) ) ERROR_exit("Can't run on a 32-bit system!") ; #ifdef ALLOW_MMAP if( nb > 1000.0 && !do_mmap ) INFO_message("If you run out of memory, try using the -mmap option.") ; #endif } tross_Make_History( "3dAutoTcorrelate" , argc,argv , cset ) ; /*--- make output dataset first ---*/ imap = (int *)calloc(sizeof(int),nmask) ; if( !do_mmap ) ININFO_message("creating output dataset in memory") ; for( kout=ii=0 ; ii < nvox ; ii++ ){ if( mask != NULL && mask[ii] == 0 ) continue ; /* skip it */ EDIT_BRICK_TO_NOSTAT(cset,kout) ; /* stat params */ EDIT_BRICK_FACTOR(cset,kout,0.0001) ; /* scale factor */ ix = DSET_index_to_ix(cset,ii) ; /* brick label */ jy = DSET_index_to_jy(cset,ii) ; kz = DSET_index_to_kz(cset,ii) ; sprintf(str,"%03d_%03d_%03d",ix,jy,kz) ; EDIT_BRICK_LABEL(cset,kout,str) ; if( !do_mmap ) EDIT_substitute_brick(cset,kout,MRI_short,NULL) ; /* make array */ imap[kout] = ii ; kout++ ; } #ifdef ALLOW_MMAP if( do_mmap ){ FILE *fp ; int64_t qq,dq ; char buf[2]={0,0} ; signal(SIGINT ,AC_sigfunc) ; /* setup signal handler */ signal(SIGBUS ,AC_sigfunc) ; /* for fatal errors */ signal(SIGSEGV,AC_sigfunc) ; signal(SIGTERM,AC_sigfunc) ; cbrik_name = strdup(DSET_BRIKNAME(cset)) ; /* save name for later use */ ncbrik = cset->dblk->total_bytes ; ININFO_message("creating output file %s via mmap() [%s bytes -- about %s]", cbrik_name , commaized_integer_string(ncbrik) , approximate_number_string((double)ncbrik) ) ; fp = fopen(cbrik_name,"w+") ; if( fp == NULL ) ERROR_exit("Can't create output BRIK file :-(") ; fseek( fp , ncbrik-1 , SEEK_SET ) ; fwrite( &buf , 1 , 2 , fp ) ; fflush(fp) ; cbrik = mmap( 0 , (size_t)ncbrik , PROT_READ | PROT_WRITE , MFLAG , fileno(fp) , 0 ) ; fclose(fp) ; if( cbrik == (char *)(-1) ){ remove(cbrik_name) ; ERROR_exit("Can't mmap output BRIK file :-(") ; } ININFO_message("Testing output BRIK") ; dq = ncbrik / 65536 ; if( dq < 1024 ) dq = 1024 ; for( qq=0 ; qq < ncbrik ; qq+=dq ) cbrik[qq] = 0 ; cbrik[ncbrik-1] = 0 ; ININFO_message("... done with test") ; } #endif /*---------- loop over mask voxels, correlate ----------*/ AFNI_OMP_START ; #pragma omp parallel if( nmask > 999 ) { int ii,jj,kout , ithr,nthr , vstep,vii ; float *xsar , *ysar ; short *car=NULL ; #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 #ifdef ALLOW_MMAP if( do_mmap ) car = (short *)malloc(sizeof(short)*nvox) ; #endif vstep = (int)( nmask / (nthr*50.0f) + 0.901f ) ; vii = 0 ; if( ithr == 0 ) fprintf(stderr,"Looping:") ; #pragma omp for for( kout=0 ; kout < nmask ; kout++ ){ /*----- outer voxel loop -----*/ ii = imap[kout] ; /* ii= source voxel (we know that ii is in the mask) */ if( ithr == 0 && vstep > 2 ) /* allow small dsets 16 Jun 2011 [rickr] */ { vii++ ; if( vii%vstep == vstep/2 ) vstep_print(); } /* get ref time series from this voxel */ xsar = VECTIM_PTR(xvectim,ii) ; if( !do_mmap ) car = DSET_ARRAY(cset,kout) ; for( jj=0 ; jj < nvox ; jj++ ){ /*----- inner loop over voxels -----*/ car[jj] = 0 ; /* default value */ /* skip unmasked voxels, unless want correlations from all source voxels */ if( !all_source && maskinner != NULL && maskinner[jj] == 0 ) continue ; #if 0 if( jj == kout ){ car[jj] = 10000 ; continue ; } /* correlation = 1.0 */ #endif ysar = VECTIM_PTR(xvectim,jj) ; car[jj] = (short)(10000.49f*corfun(nvals,xsar,ysar)) ; } /* end of inner loop over voxels */ #ifdef ALLOW_MMAP if( do_mmap ){ /* copy results to disk mmap now */ short *cout = ((short *)(cbrik)) + (int64_t)(nvox)*(int64_t)(kout) ; AAmemcpy( cout , car , sizeof(short)*nvox ) ; } #endif } /* end of outer loop over ref voxels */ if( do_mmap ) free(car) ; /* workspace for each thread */ if( ithr == 0 ) fprintf(stderr,".") ; } /* end OpenMP */ AFNI_OMP_END ; /*---------- Finish up ---------*/ /* write mask info (if any) to output dataset header */ if( mask != NULL ){ char *maskstring = mask_to_b64string(nvox,mask) ; THD_set_string_atr( cset->dblk , "AFNI_AUTOTCORR_MASK" , maskstring ) ; free(maskstring) ; free(mask) ; if( ISVALID_DSET(mset) ){ THD_set_string_atr( cset->dblk , "AFNI_AUTOCORR_MASK_IDCODE" , DSET_IDCODE_STR(mset) ) ; THD_set_string_atr( cset->dblk , "AFNI_AUTOCORR_MASK_NAME" , DSET_HEADNAME(mset) ) ; DSET_delete(mset) ; } } /* toss some trash */ VECTIM_destroy(xvectim) ; /* if did mmap(), finish that off as well */ if( !do_mmap ){ fprintf(stderr,"Done..\n") ; } else { #ifdef ALLOW_MMAP char *ptr = cbrik ; fprintf(stderr,"Done.[un-mmap-ing") ; for( ii=0 ; ii < nmask ; ii++ ){ /* map sub-bricks */ mri_fix_data_pointer( ptr , DSET_BRICK(cset,ii) ) ; ptr += DSET_BRICK_BYTES(cset,ii) ; } THD_load_statistics(cset) ; /* so can do stats */ fprintf(stderr,".") ; munmap(cbrik,ncbrik) ; cbrik = NULL ; fprintf(stderr,"].\n") ; #endif } if (fout1D) { /* write results to ASCII file, hopefully when masked */ if (fout1D && (!mask || !maskinner)) { ERROR_message("Option -1Dout restricted to commands using" "-mask and -mask_source options."); } else { float *far = (float *)calloc(nmask, sizeof(float)); int jj; fprintf(fout1D,"#Text output of:\n#"); for (ii=0; iidblk->total_bytes)) ; THD_set_write_compression(COMPRESS_NONE) ; AFNI_setenv("AFNI_AUTOGZIP=NO") ; DSET_write(cset) ; } else { INFO_message("Writing output dataset header") ; DSET_write_header(cset) ; } WROTE_DSET(cset) ; exit(0) ; }