Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
3ddot.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include <string.h>
00009
00010 double DSET_cor( THD_3dim_dataset *, THD_3dim_dataset *, byte *, int ) ;
00011
00012 int main( int argc , char * argv[] )
00013 {
00014 double dxy ;
00015 int narg , ndset , nvox , demean=0 ;
00016 THD_3dim_dataset * xset , * yset , * mask_dset=NULL ;
00017 float mask_bot=666.0 , mask_top=-666.0 ;
00018 byte * mmm=NULL ;
00019
00020
00021
00022 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){
00023 printf("Usage: 3ddot [options] dset1 dset2\n"
00024 "Output = correlation coefficient between 2 dataset bricks\n"
00025 " - you can use sub-brick selectors on the dsets\n"
00026 " - the result is a number printed to stdout"
00027 "Options:\n"
00028 " -mask mset Means to use the dataset 'mset' as a mask:\n"
00029 " Only voxels with nonzero values in 'mset'\n"
00030 " will be averaged from 'dataset'. Note\n"
00031 " that the mask dataset and the input dataset\n"
00032 " must have the same number of voxels.\n"
00033 " -mrange a b Means to further restrict the voxels from\n"
00034 " 'mset' so that only those mask values\n"
00035 " between 'a' and 'b' (inclusive) will\n"
00036 " be used. If this option is not given,\n"
00037 " all nonzero values from 'mset' are used.\n"
00038 " Note that if a voxel is zero in 'mset', then\n"
00039 " it won't be included, even if a < 0 < b.\n"
00040 " -demean Means to remove the mean from each volume\n"
00041 " prior to computing the correlation.\n"
00042 ) ;
00043
00044 printf("\n" MASTER_SHORTHELP_STRING ) ;
00045
00046 exit(0) ;
00047 }
00048
00049 narg = 1 ;
00050 while( narg < argc && argv[narg][0] == '-' ){
00051
00052 if( strncmp(argv[narg],"-demean",5) == 0 ){
00053 demean++ ;
00054 narg++ ; continue ;
00055 }
00056
00057 if( strncmp(argv[narg],"-mask",5) == 0 ){
00058 if( mask_dset != NULL ){
00059 fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
00060 }
00061 if( narg+1 >= argc ){
00062 fprintf(stderr,"*** -mask option requires a following argument!\n");
00063 exit(1) ;
00064 }
00065 mask_dset = THD_open_dataset( argv[++narg] ) ;
00066 if( mask_dset == NULL ){
00067 fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
00068 }
00069 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
00070 fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
00071 exit(1) ;
00072 }
00073 narg++ ; continue ;
00074 }
00075
00076 if( strncmp(argv[narg],"-mrange",5) == 0 ){
00077 if( narg+2 >= argc ){
00078 fprintf(stderr,"*** -mrange option requires 2 following arguments!\n")
00079 ;
00080 exit(1) ;
00081 }
00082 mask_bot = strtod( argv[++narg] , NULL ) ;
00083 mask_top = strtod( argv[++narg] , NULL ) ;
00084 if( mask_top < mask_top ){
00085 fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
00086 }
00087 narg++ ; continue ;
00088 }
00089
00090 fprintf(stderr,"*** Unknown option: %s\n",argv[narg]) ; exit(1) ;
00091 }
00092
00093
00094
00095 ndset = argc - narg ;
00096 if( ndset <= 1 ){
00097 fprintf(stderr,"*** No input datasets!?\n") ; exit(1) ;
00098 }
00099
00100 xset = THD_open_dataset( argv[narg++] ) ;
00101 yset = THD_open_dataset( argv[narg++] ) ;
00102 if( xset == NULL || yset == NULL ){
00103 fprintf(stderr,"*** cannot open both input datasets!\n") ; exit(1) ;
00104 }
00105 if( DSET_NUM_TIMES(xset) > 1 || DSET_NUM_TIMES(yset) > 1 ){
00106 fprintf(stderr,"*** cannot use time-dependent datasets!\n") ; exit(1) ;
00107 }
00108 nvox = DSET_NVOX(xset) ;
00109 if( nvox != DSET_NVOX(yset) ){
00110 fprintf(stderr,"*** input datasets dimensions don't match!\n");exit(1);
00111 }
00112
00113
00114
00115 if( mask_dset != NULL ){
00116 int mcount ;
00117 if( DSET_NVOX(mask_dset) != nvox ){
00118 fprintf(stderr,"*** Input and mask datasets are not same dimensions!\n");
00119 exit(1) ;
00120 }
00121 mmm = THD_makemask( mask_dset , 0 , mask_bot,mask_top ) ;
00122 mcount = THD_countmask( nvox , mmm ) ;
00123 fprintf(stderr,"+++ %d voxels in the mask\n",mcount) ;
00124 if( mcount <= 5 ){
00125 fprintf(stderr,"*** Mask is too small!\n");exit(1);
00126 }
00127 DSET_delete(mask_dset) ;
00128 }
00129
00130 dxy = DSET_cor( xset , yset , mmm , demean ) ;
00131 printf("%g\n",dxy) ;
00132 exit(0) ;
00133 }
00134
00135 double DSET_cor( THD_3dim_dataset *xset,
00136 THD_3dim_dataset *yset, byte *mmm , int dm )
00137 {
00138 double sumxx , sumyy , sumxy , tx,ty , dxy ;
00139 void * xar , * yar ;
00140 float * fxar , * fyar ;
00141 int ii , nxyz , ivx,ivy , itypx,itypy , fxar_new,fyar_new , nnn ;
00142
00143 nxyz = DSET_NVOX(xset) ;
00144
00145
00146
00147 THD_load_datablock( xset->dblk ) ;
00148 ivx = 0 ;
00149 itypx = DSET_BRICK_TYPE(xset,ivx) ;
00150 xar = DSET_ARRAY(xset,ivx) ; if( xar == NULL ) return 0.0 ;
00151 if( itypx == MRI_float ){
00152 fxar = (float *) xar ; fxar_new = 0 ;
00153 } else {
00154 fxar = (float *) malloc( sizeof(float) * nxyz ) ; fxar_new = 1 ;
00155 EDIT_coerce_type( nxyz , itypx,xar , MRI_float,fxar ) ;
00156 PURGE_DSET( xset ) ;
00157 }
00158
00159 THD_load_datablock( yset->dblk ) ;
00160 ivy = 0 ;
00161 itypy = DSET_BRICK_TYPE(yset,ivy) ;
00162 yar = DSET_ARRAY(yset,ivy) ; if( yar == NULL ) return 0.0 ;
00163 if( itypy == MRI_float ){
00164 fyar = (float *) yar ; fyar_new = 0 ;
00165 } else {
00166 fyar = (float *) malloc( sizeof(float) * nxyz ) ; fyar_new = 1 ;
00167 EDIT_coerce_type( nxyz , itypy,yar , MRI_float,fyar ) ;
00168 PURGE_DSET( yset ) ;
00169 }
00170
00171
00172
00173 if( dm ){
00174 sumxx = sumyy = 0.0 ;
00175 for( nnn=ii=0 ; ii < nxyz ; ii++ ){
00176 if( mmm == NULL || mmm[ii] ){
00177 sumxx += fxar[ii] ; sumyy += fyar[ii] ; nnn++ ;
00178 }
00179 }
00180 if( nnn < 5 ) return 0.0 ;
00181 sumxx /= nnn ; sumyy /= nnn ;
00182 for( ii=0 ; ii < nxyz ; ii++ ){
00183 if( mmm == NULL || mmm[ii] ){
00184 fxar[ii] -= sumxx ; fyar[ii] -= sumyy ;
00185 }
00186 }
00187 }
00188
00189
00190
00191 sumxx = sumyy = sumxy = 0.0 ;
00192
00193 for( ii=0 ; ii < nxyz ; ii++ ){
00194 if( mmm == NULL || mmm[ii] ){
00195 tx = fxar[ii] ; ty = fyar[ii] ;
00196 sumxx += tx * tx ;
00197 sumyy += ty * ty ;
00198 sumxy += tx * ty ;
00199 }
00200 }
00201
00202
00203
00204 if( fxar_new ) free( fxar ) ;
00205 if( fyar_new ) free( fyar ) ;
00206
00207
00208
00209 dxy = sumxx * sumyy ; if( dxy <= 0.0 ) return 0.0 ;
00210 dxy = sumxy / sqrt(dxy) ;
00211 return dxy ;
00212 }