00001 #include "mrilib.h"
00002 
00003 #define SPEARMAN 1
00004 #define QUADRANT 2
00005 #define PEARSON  3
00006 
00007 int main( int argc , char *argv[] )
00008 {
00009    THD_3dim_dataset *xset , *yset , *cset ;
00010    int nopt=1 , method=PEARSON , do_autoclip=0 ;
00011    int nvox , nvals , ii , polort=1 ;
00012    MRI_IMAGE *xsim , *ysim ;
00013    float     *xsar , *ysar , *car ;
00014    char * prefix = "Tcorr" ;
00015    byte * mmm=NULL ;
00016    MRI_IMAGE *im_ort=NULL ;            /* 13 Mar 2003 */
00017    int nort=0 ; float **fort=NULL ;
00018 
00019    /*----*/
00020 
00021    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00022       printf("Usage: 3dTcorrelate [options] xset yset\n"
00023              "Computes the correlation coefficient between corresponding voxel\n"
00024              "time series in two input 3D+time datasets 'xset' and 'yset', and\n"
00025              "stores the output in a new 1 sub-brick dataset.\n"
00026              "\n"
00027              "Options:\n"
00028              "  -pearson  = Correlation is the normal Pearson (product moment)\n"
00029              "                correlation coefficient [default].\n"
00030              "  -spearman = Correlation is the Spearman (rank) correlation\n"
00031              "                coefficient.\n"
00032              "  -quadrant = Correlation is the quadrant correlation coefficient.\n"
00033              "\n"
00034              "  -polort m = Remove polynomical trend of order 'm', for m=-1..3.\n"
00035              "                [default is m=1; removal is by least squares].\n"
00036              "                Using m=-1 means no detrending; this is only useful\n"
00037              "                for data/information that has been pre-processed.\n"
00038              "\n"
00039              "  -ort r.1D = Also detrend using the columns of the 1D file 'r.1D'.\n"
00040              "                Only one -ort option can be given.  If you want to use\n"
00041              "                more than one, create a temporary file using 1dcat.\n"
00042              "\n"
00043              "  -autoclip = Clip off low-intensity regions in the two datasets,\n"
00044              "  -automask =  so that the correlation is only computed between\n"
00045              "               high-intensity (presumably brain) voxels.  The\n"
00046              "               intensity level is determined the same way that\n"
00047              "               3dClipLevel works.\n"
00048              "\n"
00049              "  -prefix p = Save output into dataset with prefix 'p'\n"
00050              "               [default prefix is 'Tcorr'].\n"
00051              "\n"
00052              "Notes:\n"
00053              " * The output dataset is functional bucket type, with one\n"
00054              "    sub-brick, stored in floating point format.\n"
00055              " * Because both time series are detrended prior to correlation,\n"
00056              "    the results will not be identical to using FIM or FIM+ to\n"
00057              "    calculate correlations (whose ideal vector is not detrended).\n"
00058              " * This is a quick hack for Mike Beauchamp.  Thanks for you-know-what.\n"
00059              "\n"
00060              "-- RWCox - Aug 2001\n"
00061             ) ;
00062       exit(0) ;
00063    }
00064 
00065    mainENTRY("3dTCorrelate main"); machdep(); AFNI_logger("3dTcorrelate",argc,argv);
00066    PRINT_VERSION("3dTcorrelate") ;
00067 
00068    /*-- option processing --*/
00069 
00070    while( nopt < argc && argv[nopt][0] == '-' ){
00071 
00072       if( strcmp(argv[nopt],"-ort") == 0 ){           /* 13 Mar 2003 */
00073         if( im_ort != NULL ){
00074           fprintf(stderr,"** Can't have multiple -ort options!\n"); exit(1);
00075         }
00076         im_ort = mri_read_1D( argv[++nopt] ) ;
00077         if( im_ort == NULL ){
00078           fprintf(stderr,"** Can't read 1D file %s\n",argv[nopt]); exit(1);
00079         }
00080         nort = im_ort->ny ;
00081         fort = (float **) malloc( sizeof(float *)*nort ) ;
00082         for( ii=0 ; ii < nort ; ii++ )
00083           fort[ii] = MRI_FLOAT_PTR(im_ort) + (ii*im_ort->nx) ;
00084 
00085         nopt++ ; continue ;
00086       }
00087 
00088       if( strcmp(argv[nopt],"-autoclip") == 0 ||
00089           strcmp(argv[nopt],"-automask") == 0   ){
00090 
00091          do_autoclip = 1 ; nopt++ ; continue ;
00092       }
00093 
00094       if( strcmp(argv[nopt],"-pearson") == 0 ){
00095          method = PEARSON ; nopt++ ; continue ;
00096       }
00097 
00098       if( strcmp(argv[nopt],"-spearman") == 0 ){
00099          method = SPEARMAN ; nopt++ ; continue ;
00100       }
00101 
00102       if( strcmp(argv[nopt],"-quadrant") == 0 ){
00103          method = QUADRANT ; nopt++ ; continue ;
00104       }
00105 
00106       if( strcmp(argv[nopt],"-prefix") == 0 ){
00107          prefix = argv[++nopt] ;
00108          if( !THD_filename_ok(prefix) ){
00109             fprintf(stderr,"** Illegal value after -prefix!\n");exit(1);
00110          }
00111          nopt++ ; continue ;
00112       }
00113 
00114       if( strcmp(argv[nopt],"-polort") == 0 ){
00115          char *cpt ;
00116          int val = strtod(argv[++nopt],&cpt) ;
00117          if( *cpt != '\0' || val < -1 || val > 3 ){
00118             fprintf(stderr,"** Illegal value after -polort!\n");exit(1);
00119          }
00120          polort = val ; nopt++ ; continue ;
00121       }
00122 
00123       fprintf(stderr,"** Illegal option: %s\n",argv[nopt]) ; exit(1) ;
00124    }
00125 
00126    /*-- open datasets, check for legality --*/
00127 
00128    if( nopt+1 >= argc ){
00129       fprintf(stderr,"*** Need 2 datasets on command line!?\n"); exit(1);
00130    }
00131 
00132    xset = THD_open_dataset( argv[nopt] ) ;
00133    if( xset == NULL ){
00134       fprintf(stderr,"** Can't open dataset %s\n",argv[nopt]); exit(1);
00135    }
00136    if( DSET_NUM_TIMES(xset) < 2 ){
00137       fprintf(stderr,"** Input dataset %s is not 3D+time\n",argv[nopt]); exit(1);
00138    }
00139    yset = THD_open_dataset( argv[++nopt] ) ;
00140    if( yset == NULL ){
00141       fprintf(stderr,"** Can't open dataset %s\n",argv[nopt]); exit(1);
00142    }
00143    if( DSET_NUM_TIMES(yset) != DSET_NUM_TIMES(xset) ){
00144       fprintf(stderr,"** Input dataset %s is different length than %s\n",
00145               argv[nopt],argv[nopt-1]) ;
00146       exit(1) ;
00147    }
00148    if( DSET_NVOX(yset) != DSET_NVOX(xset) ){
00149       fprintf(stderr,"** Input dataset %s is different size than %s\n",
00150               argv[nopt],argv[nopt-1]) ;
00151       exit(1) ;
00152    }
00153    if( im_ort != NULL && im_ort->nx < DSET_NUM_TIMES(xset) ){
00154       fprintf(stderr,"** Input datsets are longer than -ort file!\n"); exit(1);
00155    }
00156    DSET_load(xset) ;
00157    if( !DSET_LOADED(xset) ){
00158       fprintf(stderr,"*** Can't read dataset bricks from %s\n",argv[nopt-1]);
00159       exit(1);
00160    }
00161    DSET_load(yset) ;
00162    if( !DSET_LOADED(yset) ){
00163       fprintf(stderr,"*** Can't read dataset bricks from %s\n",argv[nopt]);
00164       exit(1);
00165    }
00166 
00167    /*-- compute mask array, if desired --*/
00168 
00169    nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
00170 
00171    if( do_autoclip ){
00172       byte *xmm , *ymm ;
00173       xmm = THD_automask( xset ) ;
00174       ymm = THD_automask( yset ) ;
00175       mmm = (byte *) malloc(sizeof(byte)*nvox) ;
00176       for( ii=0 ; ii < nvox ; ii++ )
00177          mmm[ii] = ( xmm[ii] && ymm[ii] ) ;
00178       free(xmm) ; free(ymm) ;
00179       ii = THD_countmask( nvox , mmm ) ;
00180       fprintf(stderr,"++ %d voxels survive -autoclip\n",ii) ;
00181       if( ii == 0 ) exit(1) ;
00182    }
00183 
00184    /*-- create output dataset --*/
00185 
00186    cset = EDIT_empty_copy( xset ) ;
00187    EDIT_dset_items( cset ,
00188                       ADN_prefix    , prefix         ,
00189                       ADN_nvals     , 1              ,
00190                       ADN_ntt       , 0              , /* no time axis */
00191                       ADN_type      , HEAD_FUNC_TYPE ,
00192                       ADN_func_type , FUNC_BUCK_TYPE ,
00193                     ADN_none ) ;
00194 
00195    if( THD_is_file(DSET_HEADNAME(cset)) ){
00196       fprintf(stderr,"** Output dataset %s already exists!\n",
00197               DSET_HEADNAME(cset)) ;
00198       exit(1) ;
00199    }
00200 
00201    EDIT_BRICK_TO_FICO(cset,0,nvals,1,polort+1+nort) ;  /* stat params */
00202    EDIT_BRICK_FACTOR(cset,0,0.0) ;                     /* to be safe  */
00203 
00204    switch( method ){                                   /* looks nice  */
00205       default:
00206       case PEARSON:  EDIT_BRICK_LABEL(cset,0,"Pear.Corr.") ; break ;
00207       case SPEARMAN: EDIT_BRICK_LABEL(cset,0,"Spmn.Corr.") ; break ;
00208       case QUADRANT: EDIT_BRICK_LABEL(cset,0,"Quad.Corr.") ; break ;
00209    }
00210 
00211    EDIT_substitute_brick( cset , 0 , MRI_float , NULL ) ; /* make array  */
00212    car = DSET_ARRAY(cset,0) ;                             /* get array   */
00213 
00214    tross_Make_History( "3dTcorrelate" , argc,argv , cset ) ;
00215 
00216    /* loop over voxels, correlate */
00217    /* fprintf(stderr,"have %d voxels to work with, %d values/time series.\n", nvox, nvals);*/
00218    for( ii=0 ; ii < nvox ; ii++ ){
00219 
00220       if( mmm != NULL && mmm[ii] == 0 ){  /* the easy case */
00221          car[ii] = 0.0 ; continue ;
00222       }
00223 
00224       /* get time series */
00225 
00226       xsim = THD_extract_series(ii,xset,0) ; xsar = MRI_FLOAT_PTR(xsim) ;
00227       ysim = THD_extract_series(ii,yset,0) ; ysar = MRI_FLOAT_PTR(ysim) ;
00228 
00229       THD_generic_detrend( nvals,xsar, polort, nort,fort ) ;  /* 13 Mar 2003 */
00230       THD_generic_detrend( nvals,ysar, polort, nort,fort ) ;
00231 
00232       switch( method ){                    /* correlate */
00233          default:
00234          case PEARSON:  car[ii] = THD_pearson_corr ( nvals,xsar,ysar ); break;
00235          case SPEARMAN: car[ii] = THD_spearman_corr( nvals,xsar,ysar ); break;
00236          case QUADRANT: car[ii] = THD_quadrant_corr( nvals,xsar,ysar ); break;
00237       }
00238 
00239       mri_free(xsim) ; mri_free(ysim) ;    /* toss time series */
00240 
00241    } /* end of loop over voxels */
00242 
00243    /* toss the other trash */
00244 
00245    DSET_unload(xset); DSET_unload(yset); if( mmm != NULL ) free(mmm);
00246 
00247    /* finito */
00248 
00249    DSET_write(cset) ;
00250    fprintf(stderr,"++ Wrote dataset: %s\n",DSET_BRIKNAME(cset)) ;
00251    exit(0) ;
00252 }