Doxygen Source Code Documentation
3dTcorrelate.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | SPEARMAN 1 |
#define | QUADRANT 2 |
#define | PEARSON 3 |
Functions | |
int | main (int argc, char *argv[]) |
Define Documentation
|
Definition at line 5 of file 3dTcorrelate.c. Referenced by main(). |
|
Definition at line 4 of file 3dTcorrelate.c. Referenced by main(). |
|
Definition at line 3 of file 3dTcorrelate.c. Referenced by main(). |
Function Documentation
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 7 of file 3dTcorrelate.c. References ADN_func_type, ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, ADN_type, AFNI_logger(), argc, DSET_ARRAY, DSET_BRIKNAME, DSET_HEADNAME, DSET_load, DSET_LOADED, DSET_NUM_TIMES, DSET_NVALS, DSET_NVOX, DSET_unload, DSET_write, EDIT_BRICK_FACTOR, EDIT_BRICK_LABEL, EDIT_BRICK_TO_FICO, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), free, FUNC_BUCK_TYPE, HEAD_FUNC_TYPE, machdep(), mainENTRY, malloc, mmm, MRI_FLOAT_PTR, mri_free(), mri_read_1D(), MRI_IMAGE::nx, MRI_IMAGE::ny, PEARSON, PRINT_VERSION, QUADRANT, SPEARMAN, strtod(), THD_automask(), THD_countmask(), THD_extract_series(), THD_filename_ok(), THD_generic_detrend(), THD_is_file(), THD_open_dataset(), THD_pearson_corr(), THD_quadrant_corr(), THD_spearman_corr(), and tross_Make_History().
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 } |