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 , *cset ;
00010 int nopt=1 , method=PEARSON , do_autoclip=0 ;
00011 int nvox , nvals , ii,jj,kout , polort=1 , ix,jy,kz ;
00012 MRI_IMAGE *xsim , *ysim ;
00013 float *xsar , *ysar ;
00014 short *car ;
00015 char * prefix = "ATcorr" ;
00016 byte * mmm=NULL ;
00017 int nmask , abuc=1 ;
00018 char str[32] ;
00019
00020
00021
00022 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00023 printf("Usage: 3dAutoTcorrelate [options] dset\n"
00024 "Computes the correlation coefficient between each pair of\n"
00025 "voxels in the input dataset, and stores the output into\n"
00026 "a new anatomical bucket dataset.\n"
00027 "\n"
00028 "Options:\n"
00029 " -pearson = Correlation is the normal Pearson (product moment)\n"
00030 " correlation coefficient [default].\n"
00031 " -spearman = Correlation is the Spearman (rank) correlation\n"
00032 " coefficient.\n"
00033 " -quadrant = Correlation is the quadrant correlation coefficient.\n"
00034 "\n"
00035 " -polort m = Remove polynomical trend of order 'm', for m=-1..3.\n"
00036 " [default is m=1; removal is by least squares].\n"
00037 " Using m=-1 means no detrending; this is only useful\n"
00038 " for data/information that has been pre-processed.\n"
00039 "\n"
00040 " -autoclip = Clip off low-intensity regions in the dataset,\n"
00041 " -automask = so that the correlation is only computed between\n"
00042 " high-intensity (presumably brain) voxels. The\n"
00043 " intensity level is determined the same way that\n"
00044 " 3dClipLevel works.\n"
00045 "\n"
00046 " -prefix p = Save output into dataset with prefix 'p'\n"
00047 " [default prefix is 'ATcorr'].\n"
00048 "\n"
00049 " -time = Save output as a 3D+time dataset instead\n"
00050 " of a anat bucket.\n"
00051 "\n"
00052 "Notes:\n"
00053 " * The output dataset is anatomical bucket type of shorts.\n"
00054 " * The output file might be gigantic and you might run out\n"
00055 " of memory running this program. Use at your own risk!\n"
00056 " * The program prints out an estimate of its memory usage\n"
00057 " when it starts. It also prints out a progress 'meter'\n"
00058 " of 1 dot per 10 output sub-bricks.\n"
00059 " * This is a quick hack for Peter Bandettini. Now pay up.\n"
00060 "\n"
00061 "-- RWCox - Jan 31 2002\n"
00062 ) ;
00063 exit(0) ;
00064 }
00065
00066 mainENTRY("3dAutoTcorrelate main"); machdep(); PRINT_VERSION("3dAutoTcorrelate");
00067 AFNI_logger("3dAutoTcorrelate",argc,argv);
00068
00069
00070
00071 while( nopt < argc && argv[nopt][0] == '-' ){
00072
00073 if( strcmp(argv[nopt],"-time") == 0 ){
00074 abuc = 0 ; nopt++ ; continue ;
00075 }
00076
00077 if( strcmp(argv[nopt],"-autoclip") == 0 ||
00078 strcmp(argv[nopt],"-automask") == 0 ){
00079
00080 do_autoclip = 1 ; nopt++ ; continue ;
00081 }
00082
00083 if( strcmp(argv[nopt],"-pearson") == 0 ){
00084 method = PEARSON ; nopt++ ; continue ;
00085 }
00086
00087 if( strcmp(argv[nopt],"-spearman") == 0 ){
00088 method = SPEARMAN ; nopt++ ; continue ;
00089 }
00090
00091 if( strcmp(argv[nopt],"-quadrant") == 0 ){
00092 method = QUADRANT ; nopt++ ; continue ;
00093 }
00094
00095 if( strcmp(argv[nopt],"-prefix") == 0 ){
00096 prefix = argv[++nopt] ;
00097 if( !THD_filename_ok(prefix) ){
00098 fprintf(stderr,"** Illegal value after -prefix!\n");exit(1);
00099 }
00100 nopt++ ; continue ;
00101 }
00102
00103 if( strcmp(argv[nopt],"-polort") == 0 ){
00104 char *cpt ;
00105 int val = strtod(argv[++nopt],&cpt) ;
00106 if( *cpt != '\0' || val < -1 || val > 3 ){
00107 fprintf(stderr,"** Illegal value after -polort!\n");exit(1);
00108 }
00109 polort = val ; nopt++ ; continue ;
00110 }
00111
00112 fprintf(stderr,"** Illegal option: %s\n",argv[nopt]) ; exit(1) ;
00113 }
00114
00115
00116
00117 if( nopt >= argc ){
00118 fprintf(stderr,"*** Need a dataset on command line!?\n"); exit(1);
00119 }
00120
00121 xset = THD_open_dataset( argv[nopt] ) ;
00122 if( xset == NULL ){
00123 fprintf(stderr,"** Can't open dataset %s\n",argv[nopt]); exit(1);
00124 }
00125 if( DSET_NUM_TIMES(xset) < 2 ){
00126 fprintf(stderr,"** Input dataset %s is not 3D+time\n",argv[nopt]); exit(1);
00127 }
00128 DSET_load(xset) ;
00129 if( !DSET_LOADED(xset) ){
00130 fprintf(stderr,"*** Can't read dataset bricks from %s\n",argv[nopt-1]);
00131 exit(1);
00132 }
00133
00134
00135
00136 nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
00137
00138 if( do_autoclip ){
00139 mmm = THD_automask( xset ) ;
00140 nmask = THD_countmask( nvox , mmm ) ;
00141 fprintf(stderr,"++ %d voxels survive -autoclip\n",nmask) ;
00142 if( nmask < 2 ) exit(1) ;
00143 } else {
00144 nmask = nvox ;
00145 fprintf(stderr,"++ computing for all %d voxels\n",nmask) ;
00146 }
00147
00148
00149
00150 cset = EDIT_empty_copy( xset ) ;
00151
00152 if( abuc ){
00153 EDIT_dset_items( cset ,
00154 ADN_prefix , prefix ,
00155 ADN_nvals , nmask ,
00156 ADN_ntt , 0 ,
00157 ADN_type , HEAD_ANAT_TYPE ,
00158 ADN_func_type , ANAT_BUCK_TYPE ,
00159 ADN_none ) ;
00160 } else {
00161 EDIT_dset_items( cset ,
00162 ADN_prefix , prefix ,
00163 ADN_nvals , nmask ,
00164 ADN_ntt , nmask ,
00165 ADN_ttdel , 1.0 ,
00166 ADN_nsl , 0 ,
00167 ADN_type , HEAD_ANAT_TYPE ,
00168 ADN_func_type , ANAT_EPI_TYPE ,
00169 ADN_none ) ;
00170 }
00171
00172 if( THD_is_file(DSET_HEADNAME(cset)) ){
00173 fprintf(stderr,"** Output dataset %s already exists!\n",
00174 DSET_HEADNAME(cset)) ;
00175 exit(1) ;
00176 }
00177
00178 { double nb = (double)(xset->dblk->total_bytes) ;
00179 nb += (double)(nmask) * (double)(nvox) * sizeof(short) ;
00180 nb /= (1024.0*1024.0) ;
00181 fprintf(stderr,
00182 "++ Memory required = %.1f Mbytes for %d output sub-bricks\n",nb,nmask);
00183 if( nb > 2000.0 ){
00184 fprintf(stderr,"** ERROR: can't use more than 2000 Mbytes of memory!\n");
00185 exit(1) ;
00186 }
00187 }
00188
00189 tross_Make_History( "3dAutoTcorrelate" , argc,argv , cset ) ;
00190
00191
00192
00193 kout = 0 ;
00194 for( ii=0 ; ii < nvox ; ii++ ){
00195
00196 if( mmm != NULL && mmm[ii] == 0 ) continue ;
00197
00198 if( kout > 0 && kout%10 == 0 )
00199 fprintf(stderr, (kout%1000==0) ? "!"
00200 :(kout%100 ==0) ? ":" : "." ) ;
00201
00202
00203
00204 EDIT_substitute_brick( cset,kout, MRI_short, NULL ) ;
00205 car = DSET_ARRAY(cset,kout) ;
00206
00207 EDIT_BRICK_TO_FICO(cset,kout,nvals,1,polort+1) ;
00208 EDIT_BRICK_FACTOR(cset,kout,0.0001) ;
00209
00210 ix = DSET_index_to_ix(cset,ii) ;
00211 jy = DSET_index_to_jy(cset,ii) ;
00212 kz = DSET_index_to_kz(cset,ii) ;
00213 sprintf(str,"%03d_%03d_%03d",ix,jy,kz) ;
00214 EDIT_BRICK_LABEL(cset,kout,str) ;
00215
00216
00217
00218 xsim = THD_extract_series(ii,xset,0) ; xsar = MRI_FLOAT_PTR(xsim) ;
00219
00220 DETREND_polort(polort,nvals,xsar) ;
00221
00222 for( jj=0 ; jj < nvox ; jj++ ){
00223
00224 if( mmm != NULL && mmm[jj] == 0 ){
00225 car[jj] = 0 ; continue ;
00226 }
00227
00228 ysim = THD_extract_series(jj,xset,0) ; ysar = MRI_FLOAT_PTR(ysim) ;
00229 DETREND_polort(polort,nvals,ysar) ;
00230
00231 switch( method ){
00232 default:
00233 case PEARSON:
00234 car[jj] = rint(10000.0*THD_pearson_corr (nvals,xsar,ysar));
00235 break;
00236 case SPEARMAN:
00237 car[jj] = rint(10000.0*THD_spearman_corr(nvals,xsar,ysar));
00238 break;
00239 case QUADRANT:
00240 car[jj] = rint(10000.0*THD_quadrant_corr(nvals,xsar,ysar));
00241 break;
00242 }
00243
00244 mri_free(ysim) ;
00245
00246 }
00247
00248 mri_free(xsim) ;
00249 kout++ ;
00250
00251 }
00252
00253
00254
00255 DSET_unload(xset); if( mmm != NULL ) free(mmm);
00256
00257
00258
00259
00260 DSET_write(cset) ;
00261 fprintf(stderr,"++ Wrote output: %s\n",DSET_BRIKNAME(cset)) ;
00262 exit(0) ;
00263 }