Doxygen Source Code Documentation
3ddot.c File Reference
#include "mrilib.h"
#include <string.h>
Go to the source code of this file.
Functions | |
double | DSET_cor (THD_3dim_dataset *, THD_3dim_dataset *, byte *, int) |
int | main (int argc, char *argv[]) |
Function Documentation
|
Definition at line 135 of file 3ddot.c. References THD_3dim_dataset::dblk, DSET_ARRAY, DSET_BRICK_TYPE, DSET_NVOX, EDIT_coerce_type(), free, malloc, mmm, PURGE_DSET, and THD_load_datablock().
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 /* load bricks */ 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 /* 29 Feb 2000: remove mean? */ 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 ; /* ERROR */ 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 /* compute sums */ 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 /* toss trash */ 00203 00204 if( fxar_new ) free( fxar ) ; 00205 if( fyar_new ) free( fyar ) ; 00206 00207 /* compute result */ 00208 00209 dxy = sumxx * sumyy ; if( dxy <= 0.0 ) return 0.0 ; 00210 dxy = sumxy / sqrt(dxy) ; 00211 return dxy ; 00212 } |
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 12 of file 3ddot.c. References argc, DSET_BRICK_TYPE, DSET_cor(), DSET_delete, DSET_NUM_TIMES, DSET_NVOX, MASTER_SHORTHELP_STRING, mmm, strtod(), THD_countmask(), THD_makemask(), and THD_open_dataset().
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 /*-- read command line arguments --*/ 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 /* should have at least 2 more arguments */ 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 /* make a byte mask from mask dataset */ 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 } |