Doxygen Source Code Documentation
3dTqual.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | SPEARMAN 1 |
#define | QUADRANT 2 |
Functions | |
void | rank_order_float (int n, float *a) |
float | spearman_rank_prepare (int n, float *a) |
float | quadrant_corr_prepare (int n, float *a) |
float | spearman_rank_corr (int n, float *x, float rv, float *r) |
float | quadrant_corr (int n, float *x, float rv, float *r) |
int | main (int argc, char *argv[]) |
Define Documentation
|
Definition at line 135 of file 3dTqual.c. Referenced by main(). |
|
Definition at line 134 of file 3dTqual.c. Referenced by main(). |
Function Documentation
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 137 of file 3dTqual.c. References AFNI_logger(), argc, DSET_load, DSET_LOADED, DSET_NUM_TIMES, DSET_NVALS, DSET_NVOX, DSET_unload_one, machdep(), mainENTRY, malloc, MRI_FLOAT_PTR, mri_free(), PRINT_VERSION, qmedmad_float(), QUADRANT, quadrant_corr(), quadrant_corr_prepare(), SPEARMAN, spearman_rank_corr(), spearman_rank_prepare(), strtod(), THD_cliplevel(), THD_extract_float_brick(), THD_median_brick(), THD_open_dataset(), and var.
00138 { 00139 THD_3dim_dataset *dset ; 00140 int nopt=1 , method=SPEARMAN , do_range=0 , do_autoclip=0 ; 00141 int nvox , nvals , ii,iv,jj ; 00142 float *medar, *var , rv , *corr , cmed,cmad,cbot,ctop , clip_val=0.0 ; 00143 MRI_IMAGE *tsim , *medim ; 00144 int nkeep , *keep ; 00145 float *mkeep , *tkeep ; 00146 00147 /*----*/ 00148 00149 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ 00150 printf("Usage: 3dTqual [options] dataset\n" 00151 "Computes a `quality index' for each sub-brick in a 3D+time dataset.\n" 00152 "The output is a 1D time series with the index for each sub-brick.\n" 00153 "The results are written to stdout.\n" 00154 "\n" 00155 "Note that small values of the index are 'good', indicating that\n" 00156 "the sub-brick is not very different from the norm. The purpose\n" 00157 "of this program is to provide a crude way of screening FMRI\n" 00158 "time series for sporadic abnormal images, such as might be\n" 00159 "caused by large subject head motion or scanner glitches.\n" 00160 "\n" 00161 "Do not take the results of this program too literally. It\n" 00162 "is intended as a GUIDE to help you find data problems, and no\n" 00163 "more. It is not an assurance that the dataset is good, and\n" 00164 "it may indicate problems where nothing is wrong.\n" 00165 "\n" 00166 "Sub-bricks with index values much higher than others should be\n" 00167 "examined for problems. How you determine what 'much higher' means\n" 00168 "is mostly up to you. I suggest graphical inspection of the indexes\n" 00169 "(cf. EXAMPLE, infra). As a guide, the program will print (stderr)\n" 00170 "the median quality index and the range median-3.5*MAD .. median+3.5*MAD\n" 00171 "(MAD=Median Absolute Deviation). Values well outside this range might\n" 00172 "be considered suspect; if the quality index were normally distributed,\n" 00173 "then values outside this range would occur only about 1%% of the time.\n" 00174 "\n" 00175 "OPTIONS:\n" 00176 " -spearman = Quality index is 1 minus the Spearman (rank)\n" 00177 " correlation coefficient of each sub-brick\n" 00178 " with the median sub-brick.\n" 00179 " [This is the default method.]\n" 00180 " -quadrant = Similar to -spearman, but using 1 minus the\n" 00181 " quadrant correlation coefficient as the\n" 00182 " quality index.\n" 00183 "\n" 00184 " -autoclip = Clip off low-intensity regions in the median sub-brick,\n" 00185 " -automask = so that the correlation is only computed between\n" 00186 " high-intensity (presumably brain) voxels. The\n" 00187 " intensity level is determined the same way that\n" 00188 " 3dClipLevel works. This prevents the vast number\n" 00189 " of nearly 0 voxels outside the brain from biasing\n" 00190 " the correlation coefficient calculations.\n" 00191 #if 1 00192 "\n" 00193 " -clip val = Clip off values below 'val' in the median sub-brick.\n" 00194 #endif 00195 "\n" 00196 " -range = Print the median-3.5*MAD and median+3.5*MAD values\n" 00197 " out with EACH quality index, so that they\n" 00198 " can be plotted (cf. Example, infra).\n" 00199 " Notes: * These values are printed to stderr in any case.\n" 00200 " * This is only useful for plotting with 1dplot.\n" 00201 " * The lower value median-3.5*MAD is never allowed\n" 00202 " to go below 0.\n" 00203 "\n" 00204 "EXAMPLE:\n" 00205 " 3dTqual -range -automask fred+orig | 1dplot -one -stdin\n" 00206 "will calculate the time series of quality indexes and plot them\n" 00207 "to an X11 window, along with the median+/-3.5*MAD bands.\n" 00208 "\n" 00209 "NOTE: cf. program 3dToutcount for a somewhat different quality check.\n" 00210 "\n" 00211 "-- RWCox - Aug 2001\n" 00212 ) ; 00213 exit(0) ; 00214 } 00215 00216 mainENTRY("3dTqual main"); machdep(); AFNI_logger("3dTqual",argc,argv); 00217 PRINT_VERSION("3dTqual"); 00218 00219 /*-- option processing --*/ 00220 00221 while( nopt < argc && argv[nopt][0] == '-' ){ 00222 00223 if( strcmp(argv[nopt],"-autoclip") == 0 || 00224 strcmp(argv[nopt],"-automask") == 0 ){ 00225 00226 do_autoclip = 1 ; nopt++ ; continue ; 00227 } 00228 00229 #if 1 00230 if( strcmp(argv[nopt],"-clip") == 0 ){ 00231 do_autoclip = 0 ; 00232 clip_val = strtod(argv[++nopt],NULL) ; 00233 if( clip_val <= 0.0 ){ 00234 fprintf(stderr,"** value after -clip is illegal!\n"); exit(1); 00235 } 00236 nopt++ ;continue ; 00237 } 00238 #endif 00239 00240 if( strcmp(argv[nopt],"-range") == 0 ){ 00241 do_range = 1 ; nopt++ ; continue ; 00242 } 00243 00244 if( strcmp(argv[nopt],"-spearman") == 0 ){ 00245 method = SPEARMAN ; nopt++ ; continue ; 00246 } 00247 00248 if( strcmp(argv[nopt],"-quadrant") == 0 ){ 00249 method = QUADRANT ; nopt++ ; continue ; 00250 } 00251 00252 fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]); exit(1); 00253 } 00254 00255 /*-- open dataset --*/ 00256 00257 if( nopt >= argc ){ 00258 fprintf(stderr,"*** No dataset on command line!?\n"); exit(1); 00259 } 00260 00261 dset = THD_open_dataset( argv[nopt] ) ; 00262 if( dset == NULL ){ 00263 fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1); 00264 } 00265 if( DSET_NUM_TIMES(dset) < 2 ){ 00266 fprintf(stderr,"*** Input dataset is not 3D+time\n"); exit(1); 00267 } 00268 DSET_load(dset) ; 00269 if( !DSET_LOADED(dset) ){ 00270 fprintf(stderr,"*** Can't read dataset bricks\n"); exit(1); 00271 } 00272 00273 /*-- compute median brick, and order it for correlation --*/ 00274 00275 nvox = DSET_NVOX(dset) ; 00276 nvals = DSET_NVALS(dset) ; 00277 00278 medim = THD_median_brick( dset ) ; if( medim == NULL ) exit(1) ; 00279 medar = MRI_FLOAT_PTR(medim) ; 00280 00281 if( do_autoclip ){ 00282 clip_val = THD_cliplevel( medim , 0.5 ) ; 00283 fprintf(stderr,"++ Autoclip value = %g\n",clip_val) ; 00284 } 00285 00286 if( clip_val > 0.0 ){ 00287 nkeep = 0 ; 00288 for( ii=0 ; ii < nvox ; ii++ ) 00289 if( medar[ii] >= clip_val ) nkeep++ ; 00290 if( nkeep < 9 ){ 00291 fprintf(stderr,"*** %d voxels survive the clip: can't continue\n",nkeep) ; 00292 exit(1) ; 00293 } else if( nkeep == nvox ){ /* no clipping */ 00294 fprintf(stderr,"++ All %d voxels survive the clip\n",nvox) ; 00295 mkeep = medar ; 00296 } else { /* nkeep < nvox */ 00297 keep = (int *) malloc( sizeof(int) *nkeep ) ; 00298 mkeep = (float *) malloc( sizeof(float)*nkeep ) ; 00299 tkeep = (float *) malloc( sizeof(float)*nkeep ) ; 00300 for( jj=ii=0 ; ii < nvox ; ii++ ) 00301 if( medar[ii] >= clip_val ){ 00302 keep[jj] = ii ; mkeep[jj] = medar[ii] ; jj++ ; 00303 } 00304 mri_free(medim) ; 00305 fprintf(stderr,"++ %d out of %d voxels survive the clip\n",nkeep,nvox) ; 00306 } 00307 } else { 00308 nkeep = nvox ; mkeep = medar ; /* no clipping */ 00309 } 00310 00311 switch( method ){ 00312 default: 00313 case SPEARMAN: 00314 rv = spearman_rank_prepare( nkeep , mkeep ) ; break ; 00315 00316 case QUADRANT: 00317 rv = quadrant_corr_prepare( nkeep , mkeep ) ; break ; 00318 } 00319 00320 /*-- loop over input bricks, and correlate with median --*/ 00321 00322 corr = (float *) malloc( sizeof(float)*nvals ) ; 00323 00324 for( iv=0 ; iv < nvals ; iv++ ){ 00325 00326 /*- get sub-brick -*/ 00327 00328 tsim = THD_extract_float_brick( iv , dset ) ; 00329 if( tsim == NULL ) exit(1) ; 00330 var = MRI_FLOAT_PTR(tsim) ; 00331 00332 if( nkeep < nvox ){ 00333 for( jj=0 ; jj < nkeep ; jj++ ) tkeep[jj] = var[keep[jj]] ; 00334 } else { 00335 tkeep = var ; 00336 } 00337 00338 /*- compute correlation -*/ 00339 00340 switch( method ){ 00341 default: 00342 case SPEARMAN: 00343 corr[iv] = 1.0-spearman_rank_corr(nkeep,tkeep,rv,mkeep) ; break ; 00344 00345 case QUADRANT: 00346 corr[iv] = 1.0-quadrant_corr(nkeep,tkeep,rv,mkeep) ; break ; 00347 } 00348 00349 mri_free(tsim) ; DSET_unload_one(dset,iv) ; 00350 00351 if( !do_range ) printf( "%g\n" , corr[iv] ) ; 00352 00353 } /* end of loop over sub-bricks */ 00354 00355 /*-- now compute median and MAD of corr[] --*/ 00356 00357 qmedmad_float( nvals,corr , &cmed,&cmad ) ; 00358 00359 cbot = cmed - 3.5*cmad ; 00360 ctop = cmed + 3.5*cmad ; if( cbot < 0.0 ) cbot = 0.0 ; 00361 00362 /*-- now print results (if not already out) --*/ 00363 00364 if( do_range ){ 00365 for( iv=0 ; iv < nvals ; iv++ ) 00366 printf( "%g %g %g\n", corr[iv] , cbot,ctop ) ; 00367 } 00368 00369 fprintf(stderr,"++ Median=%g Median-3.5*MAD=%g Median+3.5*MAD=%g\n", 00370 cmed,cbot,ctop) ; 00371 exit(0) ; 00372 } |
|
Definition at line 120 of file 3dTqual.c. References quadrant_corr_prepare(), and r.
00121 { 00122 register int ii ; 00123 register float ss ; float xv ; 00124 00125 xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ; 00126 00127 for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ; 00128 00129 return ( ss/sqrt(rv*xv) ) ; 00130 } |
|
Definition at line 81 of file 3dTqual.c. References a, and rank_order_float().
00082 { 00083 register int ii ; 00084 register float rb , rs ; 00085 00086 rank_order_float( n , a ) ; 00087 00088 rb = 0.5*(n-1) ; rs=0.0 ; 00089 for( ii=0 ; ii < n ; ii++ ){ 00090 a[ii] = (a[ii] > rb) ? 1.0 00091 : (a[ii] < rb) ? -1.0 : 0.0 ; 00092 rs += a[ii]*a[ii] ; 00093 } 00094 00095 return rs ; 00096 } |
|
Definition at line 10 of file 3dTqual.c. References a, c, free, malloc, n1, and qsort_floatint().
00011 { 00012 register int ii , ns , n1 , ib ; 00013 static int nb = 0 ; 00014 static int * b = NULL ; /* workspaces */ 00015 static float * c = NULL ; 00016 float cs ; 00017 00018 /*- handle special cases -*/ 00019 00020 if( a == NULL ){ 00021 if( b != NULL ){ free(b); free(c); b=NULL ; c=NULL; nb=0; } /* free workspaces */ 00022 return ; 00023 } 00024 00025 if( n < 1 ) return ; /* meaningless input */ 00026 if( n == 1 ){ a[0] = 0.0 ; return ; } /* only one point!? */ 00027 00028 /*- make workspaces, if needed -*/ 00029 00030 if( nb < n ){ 00031 if( b != NULL ){ free(b); free(c); } 00032 b = (int *) malloc(sizeof(int )*n) ; 00033 c = (float *) malloc(sizeof(float)*n) ; 00034 nb = n ; 00035 } 00036 00037 for( ii=0 ; ii < n ; ii++ ) c[ii] = b[ii] = ii ; 00038 00039 /*- sort input, carrying b along -*/ 00040 00041 qsort_floatint( n , a , b ) ; /* see cs_sort_fi.c */ 00042 00043 /* compute ranks into c[] */ 00044 00045 n1 = n-1 ; 00046 for( ii=0 ; ii < n1 ; ii++ ){ 00047 if( a[ii] == a[ii+1] ){ /* handle ties */ 00048 cs = 2*ii+1 ; ns = 2 ; ib=ii ; ii++ ; 00049 while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; } 00050 for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ; 00051 } 00052 } 00053 00054 for( ii=0 ; ii < n ; ii++ ) a[b[ii]] = c[ii] ; 00055 00056 return ; 00057 } |
|
Definition at line 106 of file 3dTqual.c. References r, and spearman_rank_prepare().
00107 { 00108 register int ii ; 00109 register float ss ; float xv ; 00110 00111 xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0 ) return 0.0 ; 00112 00113 for( ii=0,ss=0.0 ; ii < n ; ii++ ) ss += x[ii] * r[ii] ; 00114 00115 return ( ss/sqrt(rv*xv) ) ; 00116 } |
|
Definition at line 63 of file 3dTqual.c. References a, and rank_order_float().
00064 { 00065 register int ii ; 00066 register float rb , rs ; 00067 00068 rank_order_float( n , a ) ; 00069 00070 rb = 0.5*(n-1) ; rs=0.0 ; 00071 for( ii=0 ; ii < n ; ii++ ){ 00072 a[ii] -= rb ; 00073 rs += a[ii]*a[ii] ; 00074 } 00075 00076 return rs ; 00077 } |