Doxygen Source Code Documentation
3dToutcount.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Functions | |
int | main (int argc, char *argv[]) |
Function Documentation
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 3 of file 3dToutcount.c. References ADN_brick_fac, ADN_datum_all, ADN_func_type, ADN_none, ADN_prefix, AFNI_logger(), argc, calloc, cl1_solve(), DSET_ARRAY, DSET_BRIKNAME, DSET_delete, DSET_HEADNAME, DSET_load, DSET_LOADED, DSET_NUM_TIMES, DSET_NVOX, DSET_unload, DSET_write, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), far, fit, free, FUNC_FIM_TYPE, ISVALID_DSET, log10qg(), machdep(), mainENTRY, malloc, MASTER_SHORTHELP_STRING, mmm, mmvox, MRI_FLOAT_PTR, mri_free(), PRINT_VERSION, qginv(), qmed_float(), qmedmad_float(), ref, strtod(), THD_automask(), THD_countmask(), THD_extract_series(), THD_filename_ok(), THD_insert_series(), THD_is_file(), THD_makemask(), THD_open_dataset(), tross_Copy_History(), tross_Make_History(), and var.
00004 { 00005 THD_3dim_dataset * dset , * oset=NULL ; 00006 int nvals , iv , nxyz , ii,jj , iarg , saveit=0 , oot , ic,cc ; 00007 int * count ; 00008 float qthr=0.001 , alph,fmed,fmad , fbot,ftop,fsig , sq2p,cls ; 00009 MRI_IMAGE * flim ; 00010 float * far , * var ; 00011 byte * mmm=NULL ; 00012 int mmvox=0 ; 00013 char * prefix=NULL ; 00014 int do_autoclip=0 , npass=0 , do_range=0 ; /* 12 Aug 2001 */ 00015 00016 int polort=0 , nref ; /* 07 Aug 2002 */ 00017 float **ref ; 00018 float *fit ; 00019 00020 /*----- Read command line -----*/ 00021 00022 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ 00023 printf("Usage: 3dToutcount [options] dataset\n" 00024 "Calculates number of 'outliers' a 3D+time dataset, at each\n" 00025 "time point, and writes the results to stdout.\n" 00026 "\n" 00027 "Options:\n" 00028 " -mask mset = Only count voxels in the mask dataset.\n" 00029 " -qthr q = Use 'q' instead of 0.001 in the calculation\n" 00030 " of alpha (below): 0 < q < 1.\n" 00031 "\n" 00032 " -autoclip }= Clip off 'small' voxels (as in 3dClipLevel);\n" 00033 " -automask }= you can't use this with -mask!\n" 00034 "\n" 00035 " -range = Print out median+3.5*MAD of outlier count with\n" 00036 " each time point; use with 1dplot as in\n" 00037 " 3dToutcount -range fred+orig | 1dplot -stdin -one\n" 00038 " -save ppp = Make a new dataset, and save the outlier Q in each\n" 00039 " voxel, where Q is calculated from voxel value v by\n" 00040 " Q = -log10(qg(abs((v-median)/(sqrt(PI/2)*MAD))))\n" 00041 " or Q = 0 if v is 'close' to the median (not an outlier).\n" 00042 " That is, 10**(-Q) is roughly the p-value of value v\n" 00043 " under the hypothesis that the v's are iid normal.\n" 00044 " The prefix of the new dataset (float format) is 'ppp'.\n" 00045 "\n" 00046 " -polort nn = Detrend each voxel time series with polynomials of\n" 00047 " order 'nn' prior to outlier estimation. Default\n" 00048 " value of nn=0, which means just remove the median.\n" 00049 " Detrending is done with L1 regression, not L2.\n" 00050 "\n" 00051 "OUTLIERS are defined as follows:\n" 00052 " * The trend and MAD of each time series are calculated.\n" 00053 " - MAD = median absolute deviation\n" 00054 " = median absolute value of time series minus trend.\n" 00055 " * In each time series, points that are 'far away' from the\n" 00056 " trend are called outliers, where 'far' is defined by\n" 00057 " alpha * sqrt(PI/2) * MAD\n" 00058 " alpha = qginv(0.001/N) (inverse of reversed Gaussian CDF)\n" 00059 " N = length of time series\n" 00060 " * Some outliers are to be expected, but if a large fraction of the\n" 00061 " voxels in a volume are called outliers, you should investigate\n" 00062 " the dataset more fully.\n" 00063 "\n" 00064 "Since the results are written to stdout, you probably want to redirect\n" 00065 "them to a file or another program, as in this example:\n" 00066 " 3dToutcount -automask v1+orig | 1dplot -stdin\n" 00067 "\n" 00068 "NOTE: also see program 3dTqual for a similar quality check.\n" 00069 ) ; 00070 printf("\n" MASTER_SHORTHELP_STRING ) ; 00071 exit(0) ; 00072 } 00073 00074 mainENTRY("3dToutcount main"); machdep(); AFNI_logger("3dToutcount",argc,argv); 00075 PRINT_VERSION("3dToutcount"); 00076 00077 iarg = 1 ; 00078 while( iarg < argc && argv[iarg][0] == '-' ){ 00079 00080 if( strcmp(argv[iarg],"-autoclip") == 0 || 00081 strcmp(argv[iarg],"-automask") == 0 ){ 00082 00083 if( mmm != NULL ){ 00084 fprintf(stderr,"** ERROR: can't use -autoclip/mask with -mask!\n"); 00085 exit(1) ; 00086 } 00087 do_autoclip = 1 ; iarg++ ; continue ; 00088 } 00089 00090 if( strcmp(argv[iarg],"-range") == 0 ){ /* 12 Aug 2001 */ 00091 do_range = 1 ; iarg++ ; continue ; 00092 } 00093 00094 if( strcmp(argv[iarg],"-save") == 0 ){ 00095 prefix = argv[++iarg] ; 00096 if( !THD_filename_ok(prefix) ){ 00097 fprintf(stderr,"** ERROR: -save prefix is not good!\n"); exit(1); 00098 } 00099 saveit = 1 ; iarg++ ; continue ; 00100 } 00101 00102 if( strcmp(argv[iarg],"-qthr") == 0 ){ 00103 qthr = strtod( argv[++iarg] , NULL ) ; 00104 if( qthr <= 0.0 || qthr >= 0.999 ){ 00105 fprintf(stderr,"** ERROR: -qthr value is illegal!\n"); exit(1); 00106 } 00107 iarg++ ; continue ; 00108 } 00109 00110 if( strcmp(argv[iarg],"-mask") == 0 ){ 00111 int mcount ; 00112 THD_3dim_dataset * mask_dset ; 00113 if( do_autoclip ){ 00114 fprintf(stderr,"** ERROR: can't use -mask with -autoclip/mask!\n"); 00115 exit(1) ; 00116 } 00117 mask_dset = THD_open_dataset(argv[++iarg]) ; 00118 if( mask_dset == NULL ){ 00119 fprintf(stderr,"** ERROR: can't open -mask dataset!\n"); exit(1); 00120 } 00121 if( mmm != NULL ){ 00122 fprintf(stderr,"** ERROR: can't have 2 -mask options!\n"); exit(1); 00123 } 00124 mmm = THD_makemask( mask_dset , 0 , 1.0,-1.0 ) ; 00125 mmvox = DSET_NVOX( mask_dset ) ; 00126 mcount = THD_countmask( mmvox , mmm ) ; 00127 fprintf(stderr,"++ %d voxels in the mask\n",mcount) ; 00128 if( mcount <= 5 ){ 00129 fprintf(stderr,"** Mask is too small!\n");exit(1); 00130 } 00131 DSET_delete(mask_dset) ; iarg++ ; continue ; 00132 } 00133 00134 if( strcmp(argv[iarg],"-polort") == 0 ){ 00135 polort = strtol( argv[++iarg] , NULL , 10 ) ; 00136 if( polort < 0 || polort > 3){ 00137 fprintf(stderr,"** Illegal value of polort!\n"); exit(1); 00138 } 00139 iarg++ ; continue ; 00140 } 00141 00142 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ; 00143 } 00144 00145 /*----- read input dataset -----*/ 00146 00147 if( iarg >= argc ){ 00148 fprintf(stderr,"** No input dataset!?\n"); exit(1); 00149 } 00150 00151 dset = THD_open_dataset( argv[iarg] ) ; 00152 if( !ISVALID_DSET(dset) ){ 00153 fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]); exit(1); 00154 } 00155 nvals = DSET_NUM_TIMES(dset) ; 00156 if( nvals < 5 ){ 00157 fprintf(stderr,"** Can't use dataset with < 5 time points per voxel!\n") ; 00158 exit(1) ; 00159 } 00160 nxyz = DSET_NVOX(dset) ; 00161 if( mmm != NULL && mmvox != nxyz ){ 00162 fprintf(stderr,"** Mask and input datasets not the same size!\n") ; 00163 exit(1) ; 00164 } 00165 DSET_load(dset) ; 00166 if( !DSET_LOADED(dset) ){ 00167 fprintf(stderr,"** Can't read dataset from disk!\n") ; 00168 exit(1) ; 00169 } 00170 00171 /*-- 12 Aug 2001: compute clip, if desired --*/ 00172 00173 if( do_autoclip && mmm == NULL ){ 00174 mmm = THD_automask( dset ) ; 00175 } 00176 00177 /*-- setup to save a new dataset, if desired --*/ 00178 00179 if( saveit ){ 00180 float * bar ; 00181 oset = EDIT_empty_copy( dset ) ; 00182 EDIT_dset_items( oset , 00183 ADN_prefix , prefix , 00184 ADN_brick_fac , NULL , 00185 ADN_datum_all , MRI_float , 00186 ADN_func_type , FUNC_FIM_TYPE , 00187 ADN_none ) ; 00188 00189 if( THD_is_file(DSET_HEADNAME(oset)) ){ 00190 fprintf(stderr,"** ERROR: -save dataset already exists!\n"); exit(1); 00191 } 00192 00193 tross_Copy_History( oset , dset ) ; 00194 tross_Make_History( "3dToutcount" , argc , argv , oset ) ; 00195 00196 for( iv=0 ; iv < nvals ; iv++ ){ 00197 EDIT_substitute_brick( oset , iv , MRI_float , NULL ) ; 00198 bar = DSET_ARRAY(oset,iv) ; 00199 for( ii=0 ; ii < nxyz ; ii++ ) bar[ii] = 0.0 ; 00200 } 00201 } 00202 00203 /*-- setup to count --*/ 00204 00205 sq2p = sqrt(0.5*PI) ; 00206 alph = qginv(qthr/nvals) * sq2p ; 00207 count = (int *) calloc( sizeof(int) , nvals ) ; 00208 var = (float *) malloc( sizeof(float) * nvals ) ; 00209 00210 /* 07 Aug 2002: make polort refs */ 00211 00212 nref = polort+1 ; 00213 ref = (float **) malloc( sizeof(float *) * nref ) ; 00214 for( jj=0 ; jj < nref ; jj++ ) 00215 ref[jj] = (float *) malloc( sizeof(float) * nvals ) ; 00216 00217 fit = (float *) malloc( sizeof(float) * nref ) ; 00218 00219 /* r(t) = 1 */ 00220 00221 for( iv=0 ; iv < nvals ; iv++ ) ref[0][iv] = 1.0 ; 00222 00223 jj = 1 ; 00224 if( polort > 0 ){ 00225 00226 /* r(t) = t - tmid */ 00227 00228 float tm = 0.5 * (nvals-1.0) ; float fac = 2.0 / nvals ; 00229 for( iv=0 ; iv < nvals ; iv++ ) ref[1][iv] = (iv-tm)*fac ; 00230 jj = 2 ; 00231 00232 /* r(t) = (t-tmid)**jj */ 00233 00234 for( ; jj <= polort ; jj++ ) 00235 for( iv=0 ; iv < nvals ; iv++ ) 00236 ref[jj][iv] = pow( (iv-tm)*fac , (double)jj ) ; 00237 } 00238 00239 /*--- loop over voxels and count ---*/ 00240 00241 for( cc=ii=0 ; ii < nxyz ; ii++ ){ 00242 if( mmm != NULL && mmm[ii] == 0 ) continue ; /* masked out */ 00243 00244 npass++ ; 00245 00246 flim = THD_extract_series( ii , dset , 0 ) ; /* get data */ 00247 far = MRI_FLOAT_PTR(flim) ; 00248 memcpy(var,far,sizeof(float)*nvals ) ; /* copy data */ 00249 00250 if( polort == 0 ){ /* the old way */ 00251 00252 fmed = qmed_float( nvals , far ) ; 00253 for( iv=0 ; iv < nvals ; iv++ ){ 00254 var[iv] = var[iv] - fmed ; /* remove median = resid */ 00255 far[iv] = fabs(var[iv]) ; /* abs value of resid */ 00256 } 00257 00258 } else { /* 07 Aug 2002: detrend */ 00259 00260 float val ; 00261 cls = cl1_solve( nvals , nref , far , ref , fit,0 ); /* get fit */ 00262 if( cls < 0.0 ) continue ; /* bad! should not happen */ 00263 for( iv=0 ; iv < nvals ; iv++ ){ /* detrend */ 00264 val = 0.0 ; 00265 for( jj=0 ; jj < nref ; jj++ ) /* fitted value */ 00266 val += fit[jj] * ref[jj][iv] ; 00267 00268 var[iv] = var[iv]-val ; /* remove fitted value = resid */ 00269 far[iv] = fabs(var[iv]) ; /* abs value of resid */ 00270 } 00271 } 00272 00273 /* find median of detrended data */ 00274 00275 fmad = qmed_float( nvals , far ) ; 00276 ftop = alph*fmad ; fbot = -ftop ; 00277 00278 if( fmad > 0.0 ){ 00279 if( saveit ) fsig = 1.0/(sq2p*fmad) ; 00280 for( ic=iv=0 ; iv < nvals ; iv++ ){ 00281 oot = (var[iv] < fbot || var[iv] > ftop ) ; 00282 if( oot ){ count[iv]++ ; cc++ ; } 00283 if( saveit ){ 00284 if( oot ){ far[iv] = -log10qg(fabs(var[iv]*fsig)); ic++; } 00285 else { far[iv] = 0.0 ; } 00286 } 00287 } 00288 00289 if( ic > 0 ) 00290 THD_insert_series( ii,oset, nvals,MRI_float,far , 0 ) ; 00291 } 00292 00293 mri_free(flim) ; 00294 } 00295 00296 if( saveit && cc > 0 ){ 00297 DSET_write( oset ) ; 00298 fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(oset)) ; 00299 } 00300 00301 if( do_range ){ 00302 float *ff = (float *)malloc(sizeof(float)*nvals) , cmed,cmad ; 00303 int ctop ; 00304 00305 for( iv=0 ; iv < nvals ; iv++ ) ff[iv] = count[iv] ; 00306 qmedmad_float( nvals,ff , &cmed,&cmad ) ; free(ff) ; 00307 ctop = (int)(cmed+3.5*cmad+0.499) ; 00308 for( iv=0 ; iv < nvals ; iv++ ) printf("%6d %d\n",count[iv],ctop) ; 00309 } else { 00310 for( iv=0 ; iv < nvals ; iv++ ) printf("%6d\n",count[iv]) ; 00311 } 00312 00313 #if 0 00314 DSET_unload(dset) ; free(count) ; free(var) ; if(mmm!=NULL)free(mmm) ; 00315 #endif 00316 00317 if( npass < nxyz ) 00318 fprintf(stderr,"++ %d voxels passed mask/clip\n",npass) ; 00319 00320 exit(0) ; 00321 } |