1 #include "mrilib.h"
2
3 /*==============================================================================*/
4
5 #define SPEARMAN 1
6 #define QUADRANT 2
7
main(int argc,char * argv[])8 int main( int argc , char *argv[] )
9 {
10 THD_3dim_dataset *dset ;
11 int nopt=1 , method=SPEARMAN , do_range=0 , do_autoclip=0 , do_mask=0;
12 int nvox , nvals , ii,iv,jj ;
13 float *medar, *var , rv , *corr , cmed,cmad,cbot,ctop , clip_val=0.0 ;
14 MRI_IMAGE *tsim , *medim ;
15 int nkeep=0 , *keep=NULL ;
16 float *mkeep=NULL , *tkeep=NULL ;
17 bytevec *bvec=NULL; /* vars for do_mask 15 Nov 2012 [rickr] */
18 byte *mask=NULL;
19 int nmask=-1, nmaskset=0;
20
21 /*----*/
22
23 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
24 printf("Usage: 3dTqual [options] dataset\n"
25 "Computes a `quality index' for each sub-brick in a 3D+time dataset.\n"
26 "The output is a 1D time series with the index for each sub-brick.\n"
27 "The results are written to stdout.\n"
28 "\n"
29 "Note that small values of the index are 'good', indicating that\n"
30 "the sub-brick is not very different from the norm. The purpose\n"
31 "of this program is to provide a crude way of screening FMRI\n"
32 "time series for sporadic abnormal images, such as might be\n"
33 "caused by large subject head motion or scanner glitches.\n"
34 "\n"
35 "Do not take the results of this program too literally. It\n"
36 "is intended as a GUIDE to help you find data problems, and no\n"
37 "more. It is not an assurance that the dataset is good, and\n"
38 "it may indicate problems where nothing is wrong.\n"
39 "\n"
40 "Sub-bricks with index values much higher than others should be\n"
41 "examined for problems. How you determine what 'much higher' means\n"
42 "is mostly up to you. I suggest graphical inspection of the indexes\n"
43 "(cf. EXAMPLE, infra). As a guide, the program will print (stderr)\n"
44 "the median quality index and the range median-3.5*MAD .. median+3.5*MAD\n"
45 "(MAD=Median Absolute Deviation). Values well outside this range might\n"
46 "be considered suspect; if the quality index were normally distributed,\n"
47 "then values outside this range would occur only about 1%% of the time.\n"
48 "\n"
49 "OPTIONS:\n"
50 " -spearman = Quality index is 1 minus the Spearman (rank)\n"
51 " correlation coefficient of each sub-brick\n"
52 " with the median sub-brick.\n"
53 " [This is the default method.]\n"
54 " -quadrant = Similar to -spearman, but using 1 minus the\n"
55 " quadrant correlation coefficient as the\n"
56 " quality index.\n"
57 "\n"
58 " -autoclip = Clip off low-intensity regions in the median sub-brick,\n"
59 " -automask = so that the correlation is only computed between\n"
60 " high-intensity (presumably brain) voxels. The\n"
61 " intensity level is determined the same way that\n"
62 " 3dClipLevel works. This prevents the vast number\n"
63 " of nearly 0 voxels outside the brain from biasing\n"
64 " the correlation coefficient calculations.\n"
65 #if 1
66 "\n"
67 " -clip val = Clip off values below 'val' in the median sub-brick.\n"
68 #endif
69 "\n"
70 " -mask MSET = Compute correlation only across masked voxels.\n"
71 "\n"
72 " -range = Print the median-3.5*MAD and median+3.5*MAD values\n"
73 " out with EACH quality index, so that they\n"
74 " can be plotted (cf. Example, infra).\n"
75 " Notes: * These values are printed to stderr in any case.\n"
76 " * This is only useful for plotting with 1dplot.\n"
77 " * The lower value median-3.5*MAD is never allowed\n"
78 " to go below 0.\n"
79 "\n"
80 "EXAMPLE:\n"
81 " 3dTqual -range -automask fred+orig | 1dplot -one -stdin\n"
82 "will calculate the time series of quality indexes and plot them\n"
83 "to an X11 window, along with the median+/-3.5*MAD bands.\n"
84 "\n"
85 "NOTE: cf. program 3dToutcount for a somewhat different quality check.\n"
86 "\n"
87 "-- RWCox - Aug 2001\n"
88 ) ;
89 PRINT_COMPILE_DATE ; exit(0) ;
90 }
91
92 mainENTRY("3dTqual main"); machdep(); AFNI_logger("3dTqual",argc,argv);
93 PRINT_VERSION("3dTqual");
94
95 /*-- option processing --*/
96
97 while( nopt < argc && argv[nopt][0] == '-' ){
98
99 if( strcmp(argv[nopt],"-autoclip") == 0 ||
100 strcmp(argv[nopt],"-automask") == 0 ){
101
102 do_autoclip = 1 ; nopt++ ; continue ;
103 }
104
105 #if 1
106 if( strcmp(argv[nopt],"-clip") == 0 ){
107 do_autoclip = 0 ;
108 clip_val = strtod(argv[++nopt],NULL) ;
109 if( clip_val <= 0.0 ){
110 fprintf(stderr,"** value after -clip is illegal!\n"); exit(1);
111 }
112 nopt++ ;continue ;
113 }
114 #endif
115
116 if( strcmp(argv[nopt],"-mask") == 0 ){ /* set mask and nmask */
117 if(++nopt >= argc) ERROR_exit("Need argument after '%s'",argv[nopt-1]);
118 bvec = THD_create_mask_from_string(argv[nopt]);
119 if(! bvec) ERROR_exit("Can't create mask from '-mask' option from %s",
120 argv[nopt]);
121 mask = bvec->ar; nmask = bvec->nar;
122
123 nmaskset = THD_countmask(nmask, mask);
124 INFO_message("%d voxels in -mask dataset", nmaskset);
125
126 do_mask = 1 ; nopt++ ; continue ;
127 }
128
129 if( strcmp(argv[nopt],"-range") == 0 ){
130 do_range = 1 ; nopt++ ; continue ;
131 }
132
133 if( strcmp(argv[nopt],"-spearman") == 0 ){
134 method = SPEARMAN ; nopt++ ; continue ;
135 }
136
137 if( strcmp(argv[nopt],"-quadrant") == 0 ){
138 method = QUADRANT ; nopt++ ; continue ;
139 }
140
141 fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]); exit(1);
142 }
143
144 /*-- open dataset --*/
145
146 if( nopt >= argc ){
147 fprintf(stderr,"*** No dataset on command line!?\n"); exit(1);
148 }
149
150 if( do_mask && do_autoclip )
151 ERROR_exit("Cannot use both mask and clip");
152
153 dset = THD_open_dataset( argv[nopt] ) ;
154 if( dset == NULL ){
155 fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1);
156 }
157 if( DSET_NUM_TIMES(dset) < 2 ){
158 fprintf(stderr,"*** Input dataset is not 3D+time\n"); exit(1);
159 }
160 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
161
162 /*-- compute median brick, and order it for correlation --*/
163
164 nvox = DSET_NVOX(dset) ;
165 nvals = DSET_NVALS(dset) ;
166
167 medim = THD_median_brick( dset ) ; if( medim == NULL ) exit(1) ;
168 medar = MRI_FLOAT_PTR(medim) ;
169
170 if( do_autoclip ){
171 clip_val = THD_cliplevel( medim , 0.5 ) ;
172 fprintf(stderr,"++ Autoclip value = %g\n",clip_val) ;
173 }
174
175 if( clip_val > 0.0 ){
176 nkeep = 0 ;
177 for( ii=0 ; ii < nvox ; ii++ )
178 if( medar[ii] >= clip_val ) nkeep++ ;
179 if( nkeep < 9 ){
180 fprintf(stderr,"*** %d voxels survive the clip: can't continue\n",nkeep) ;
181 exit(1) ;
182 } else if( nkeep == nvox ){ /* no clipping */
183 fprintf(stderr,"++ All %d voxels survive the clip\n",nvox) ;
184 mkeep = medar ;
185 } else { /* nkeep < nvox */
186 keep = (int *) malloc( sizeof(int) *nkeep ) ;
187 mkeep = (float *) malloc( sizeof(float)*nkeep ) ;
188 tkeep = (float *) malloc( sizeof(float)*nkeep ) ;
189 for( jj=ii=0 ; ii < nvox ; ii++ )
190 if( medar[ii] >= clip_val ){
191 keep[jj] = ii ; mkeep[jj] = medar[ii] ; jj++ ;
192 }
193 mri_free(medim) ;
194 fprintf(stderr,"++ %d out of %d voxels survive the clip\n",nkeep,nvox) ;
195 }
196 } else if ( mask ) { /* -mask option for evangelou 15 Nov 2012 [rickr] */
197 if( nmask != nvox )
198 ERROR_exit("dset has %d voxels, but mask has %d\n", nvox, nmask);
199 if( nmaskset == nvox ) {
200 nkeep = nvox ; mkeep = medar ; /* no masking */
201 } else if( nmaskset < 9 ){
202 ERROR_exit("only %d voxels in mask: can't continue\n",nmaskset) ;
203 } else { /* nkeep < nvox */
204 nkeep = nmaskset;
205 keep = (int *) malloc( sizeof(int) *nkeep ) ;
206 mkeep = (float *) malloc( sizeof(float)*nkeep ) ;
207 tkeep = (float *) malloc( sizeof(float)*nkeep ) ;
208 for( jj=ii=0 ; ii < nvox ; ii++ )
209 if( mask[ii] ){
210 keep[jj] = ii ; mkeep[jj] = medar[ii] ; jj++ ;
211 }
212 mri_free(medim) ;
213 KILL_bytevec(bvec) ;
214 fprintf(stderr,"++ %d out of %d voxels in mask\n",nkeep,nvox) ;
215 }
216 } else {
217 nkeep = nvox ; mkeep = medar ; /* no clipping */
218 }
219
220 switch( method ){
221 default:
222 case SPEARMAN:
223 rv = spearman_rank_prepare( nkeep , mkeep ) ; break ;
224
225 case QUADRANT:
226 rv = quadrant_corr_prepare( nkeep , mkeep ) ; break ;
227 }
228
229 /*-- loop over input bricks, and correlate with median --*/
230
231 corr = (float *) malloc( sizeof(float)*nvals ) ;
232
233 for( iv=0 ; iv < nvals ; iv++ ){
234
235 /*- get sub-brick -*/
236
237 tsim = THD_extract_float_brick( iv , dset ) ;
238 if( tsim == NULL ) exit(1) ;
239 var = MRI_FLOAT_PTR(tsim) ;
240
241 if( nkeep < nvox ){
242 for( jj=0 ; jj < nkeep ; jj++ ) tkeep[jj] = var[keep[jj]] ;
243 } else {
244 tkeep = var ;
245 }
246
247 /*- compute correlation -*/
248
249 switch( method ){
250 default:
251 case SPEARMAN:
252 corr[iv] = 1.0-spearman_rank_corr(nkeep,tkeep,rv,mkeep) ; break ;
253
254 case QUADRANT:
255 corr[iv] = 1.0-quadrant_corr(nkeep,tkeep,rv,mkeep) ; break ;
256 }
257
258 mri_free(tsim) ; DSET_unload_one(dset,iv) ;
259
260 if( !do_range ) printf( "%g\n" , corr[iv] ) ;
261
262 } /* end of loop over sub-bricks */
263
264 /*-- now compute median and MAD of corr[] --*/
265
266 qmedmad_float( nvals,corr , &cmed,&cmad ) ;
267
268 cbot = cmed - 3.5*cmad ;
269 ctop = cmed + 3.5*cmad ; if( cbot < 0.0 ) cbot = 0.0 ;
270
271 /*-- now print results (if not already out) --*/
272
273 if( do_range ){
274 for( iv=0 ; iv < nvals ; iv++ )
275 printf( "%g %g %g\n", corr[iv] , cbot,ctop ) ;
276 }
277
278 fprintf(stderr,"++ Median=%g Median-3.5*MAD=%g Median+3.5*MAD=%g\n",
279 cmed,cbot,ctop) ;
280 exit(0) ;
281 }
282