1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include <string.h>
8 #include "mrilib.h"
9 #include <stdlib.h>
10 #include <ctype.h>
11
12 /*-------------------------- global data --------------------------*/
13
14 /** inputs **/
15
16 static THD_3dim_dataset * UC_dset = NULL ; /* dataset */
17
18 static char UC_prefix[THD_MAX_PREFIX] = "uuu" ;
19
20 static int UC_be_quiet = 1 ;
21
22 static byte * UC_mask = NULL ;
23 static int UC_mask_nvox = 0 ;
24 static int UC_mask_hits = 0 ;
25
26 static int UC_nvec = 0 ; /* # of vectors to use from dataset */
27 static int UC_vdim = 0 ; /* length of each vector */
28
29 static float ** UC_vec = NULL ; /* UC_vec[k][i] is the i-th component */
30 /* of the k-th vector 0 <= i < UC_vdim */
31 /* 0 <= k < UC_nvec */
32
33 static int * UC_iv = NULL ; /* UC_vec[k] comes from voxel # UC_iv[k] */
34
35 static float UC_ptail = 0.0001 ;
36
37 void UC_syntax(char * msg) ;
38
39 #include "uuu.c"
40
41 /*-------------------------------------------------------------------
42 detrend: routine to remove unwanted components from time series
43 ---------------------------------------------------------------------*/
44
detrend(int n,float vec[])45 void detrend( int n , float vec[] )
46 {
47 register int ii ;
48 register float sum0 , sum1 , cf , lf ;
49 float sum2 , det ;
50
51 static int nold = -1 ; /* initialization flag */
52 static float cf0,cf1 , lf0,lf1 ; /* to be initialized */
53
54 /*** initialize coefficients for detrending ***/
55
56 if( n != nold ){
57 nold = n ; sum0 = sum1 = sum2 = 0.0 ;
58 for( ii=0 ; ii < n ; ii++ ){
59 sum0 += 1.0 ; sum1 += ii ; sum2 += ii*ii ;
60 }
61 det = sum0 * sum2 - sum1 * sum1 ;
62 cf0 = sum2 / det ; /* constant factor for sum0 */
63 cf1 = -sum1 / det ; /* constant factor for sum1 */
64 lf0 = cf1 ; /* linear factor for sum0 */
65 lf1 = sum0 / det ; /* linear factor for sum1 */
66 }
67
68 /*** remove mean and linear trend ***/
69
70 sum0 = sum1 = 0.0 ;
71 for( ii=0 ; ii < n ; ii++ ){
72 sum0 += vec[ii] ; sum1 += vec[ii] * ii ;
73 }
74
75 cf = cf0 * sum0 + cf1 * sum1 ;
76 lf = lf0 * sum0 + lf1 * sum1 ;
77 for( ii=0 ; ii < n ; ii++ ) vec[ii] -= cf + ii*lf ;
78 }
79
80 /*----------------------------------------------------------
81 normalize: routine to scale a time series to unit vector
82 ------------------------------------------------------------*/
83
normalize(int n,float vec[])84 void normalize( int n , float vec[] )
85 {
86 register int ii ;
87 register float sqsum ;
88
89 detrend( n , vec ) ;
90
91 sqsum = 0.0 ;
92 for( ii=0 ; ii < n ; ii++ ) sqsum += vec[ii] * vec[ii] ;
93
94 if( sqsum < 1.e-10 ){
95 for( ii=0 ; ii < n ; ii++ ) vec[ii] = 0.0 ;
96 } else {
97 sqsum = 1.0 / sqrt(sqsum) ;
98 for( ii=0 ; ii < n ; ii++ ) vec[ii] *= sqsum ;
99 }
100 }
101
102 /*------------------------------------------------------------------*/
103
UC_read_opts(int argc,char * argv[])104 void UC_read_opts( int argc , char * argv[] )
105 {
106 int nopt = 1 ;
107 float val ;
108 int kk, nxyz, mm,nn ;
109 float * vv , * bb ;
110
111 while( nopt < argc && argv[nopt][0] == '-' ){
112
113 /**** -verbose ****/
114
115 if( strncmp(argv[nopt],"-verbose",5) == 0 ){
116 UC_be_quiet = 0 ;
117 nopt++ ; continue ;
118 }
119
120 /**** -prefix prefix ****/
121
122 if( strncmp(argv[nopt],"-prefix",6) == 0 ){
123 nopt++ ;
124 if( nopt >= argc ) UC_syntax("-prefix needs an argument!") ;
125 MCW_strncpy( UC_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
126 continue ;
127 }
128
129 /**** -mask mset ****/
130
131 if( strncmp(argv[nopt],"-mask",5) == 0 ){
132 THD_3dim_dataset * mset ; int ii,nn ;
133 nopt++ ;
134 if( nopt >= argc ) UC_syntax("need arguments after -mask!") ;
135 mset = THD_open_dataset( argv[nopt] ) ;
136 if( mset == NULL ) UC_syntax("can't open -mask dataset!") ;
137 UC_mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
138 UC_mask_nvox = DSET_NVOX(mset) ;
139 DSET_delete(mset) ;
140 if( UC_mask == NULL ) UC_syntax("can't use -mask dataset!") ;
141 UC_mask_hits = THD_countmask( UC_mask_nvox , UC_mask ) ;
142 if( UC_mask_hits == 0 ) UC_syntax("mask is all zeros!") ;
143 if( !UC_be_quiet ) printf("--- %d voxels in mask\n",UC_mask_hits) ;
144
145 UC_iv = (int *) malloc( sizeof(int) * UC_mask_hits ) ;
146 for( nn=ii=0 ; ii < UC_mask_nvox ; ii++ )
147 if( UC_mask[ii] ) UC_iv[nn++] = ii ;
148
149 nopt++ ; continue ;
150 }
151
152 /**** -ptail p ****/
153
154 if( strcmp(argv[nopt],"-ptail") == 0 ){
155 if( ++nopt >= argc ) UC_syntax("-ptail needs an argument!") ;
156 UC_ptail = strtod( argv[nopt] , NULL ) ;
157 if( UC_ptail <= 0.0 || UC_ptail >= 0.499 )
158 UC_syntax("value after -ptail is illegal!") ;
159 nopt++ ; continue ;
160 }
161
162 /**** unknown switch ****/
163
164 fprintf(stderr,"\n*** unrecognized option %s\n",argv[nopt]) ;
165 exit(1) ;
166
167 } /* end of loop over options */
168
169 /*--- a simple consistency check ---*/
170
171 /*--- last input is dataset name ---*/
172
173 if( nopt >= argc ) UC_syntax("no input dataset name?") ;
174
175 UC_dset = THD_open_dataset( argv[nopt] ) ;
176 if( !ISVALID_3DIM_DATASET(UC_dset) ){
177 fprintf(stderr,"\n*** can't open dataset file %s\n",argv[nopt]) ;
178 exit(1) ;
179 }
180
181 nxyz = DSET_NVOX(UC_dset) ;
182 if( UC_mask != NULL && nxyz != UC_mask_nvox )
183 UC_syntax("mask and input dataset size mismatch!") ;
184
185 /*--- load vectors ---*/
186
187 UC_nvec = (UC_mask_hits > 0) ? UC_mask_hits : nxyz ;
188 UC_vdim = DSET_NVALS(UC_dset) ;
189 if( UC_vdim < 4 )
190 UC_syntax("input dataset needs at least 4 sub-bricks!") ;
191
192 vv = (float *) malloc( sizeof(float) * UC_nvec * UC_vdim ) ;
193 UC_vec = (float **) malloc( sizeof(float *) * UC_nvec ) ;
194 for( kk=0 ; kk < UC_nvec ; kk++ ) UC_vec[kk] = vv + (kk*UC_vdim) ;
195
196 if( !UC_be_quiet ) printf("--- reading input dataset\n") ;
197 DSET_load(UC_dset) ; CHECK_LOAD_ERROR(UC_dset) ;
198
199 /* copy brick data into float storage */
200
201 if( !UC_be_quiet ) printf("--- loading vectors\n") ;
202
203 bb = (float *) malloc( sizeof(float) * nxyz ) ;
204 for( mm=0 ; mm < UC_vdim ; mm++ ){
205
206 EDIT_coerce_type( nxyz ,
207 DSET_BRICK_TYPE(UC_dset,mm) , DSET_ARRAY(UC_dset,mm) ,
208 MRI_float , bb ) ;
209
210 DSET_unload_one( UC_dset , mm ) ;
211
212 if( UC_mask == NULL ){
213 for( kk=0 ; kk < nxyz ; kk++ ) UC_vec[kk][mm] = bb[kk] ;
214 } else {
215 for( nn=kk=0 ; kk < nxyz ; kk++ )
216 if( UC_mask[kk] ) UC_vec[nn++][mm] = bb[kk] ;
217 }
218 }
219 free(bb) ; DSET_unload( UC_dset ) ;
220
221 /* detrend and normalize vectors */
222
223 if( !UC_be_quiet ) printf("--- normalizing vectors\n") ;
224
225 for( kk=0 ; kk < UC_nvec ; kk++ )
226 normalize( UC_vdim , UC_vec[kk] ) ;
227
228 return ;
229 }
230
231 /*-----------------------------------------------------------------------
232 Compute the unusuality index of a reference vector in a set of vectors
233 -------------------------------------------------------------------------*/
234
UC_unusuality(int ndim,float * ref,int nvec,float ** vec)235 float UC_unusuality( int ndim , float * ref , int nvec , float ** vec )
236 {
237 register int ii , kk ;
238 register float psum , * vv ;
239
240 static int nvold=-1 ;
241 static float * zval =NULL ;
242
243 if( ndim < 4 || nvec < 4 || ref == NULL || vec == NULL ) return 0.0 ;
244
245 /* initialize if number of vectors has changed */
246
247 if( nvold != nvec ){
248 if( zval != NULL ) free(zval) ;
249 zval = (float *) malloc(sizeof(float)*nvec) ;
250 nvold = nvec ;
251 }
252
253 /* compute dot products */
254
255 for( kk=0 ; kk < nvec ; kk++ ){
256 psum = 0.0 ; vv = vec[kk] ;
257 for( ii=0 ; ii < ndim ; ii++ ) psum += ref[ii] * vv[ii] ;
258 zval[kk] = psum ;
259 }
260
261 psum = unusuality( nvec, zval ) ;
262
263 return psum ;
264 }
265
266 /*---------------------------------------------------------------------------*/
267
UC_syntax(char * msg)268 void UC_syntax(char * msg)
269 {
270 if( msg != NULL ){ fprintf(stderr,"\n*** %s\n",msg) ; exit(1) ; }
271
272 printf(
273 "Usage: 3duuu [options] dataset ...\n"
274 "\n"
275 "The input dataset may have a sub-brick selector list.\n"
276 "Otherwise, all sub-bricks from a dataset will be used.\n"
277 "\n"
278 "OPTIONS:\n"
279 " -prefix pname \n"
280 " -verbose\n"
281 " -mask mset\n"
282 " -ptail p\n"
283 ) ;
284
285 exit(0) ;
286 }
287
288 /*------------------------------------------------------------------*/
289
main(int argc,char * argv[])290 int main( int argc , char * argv[] )
291 {
292 int kk , nvox , ii ;
293 THD_3dim_dataset * oset ;
294 float * far ;
295
296 /*-- read command line arguments --*/
297
298 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) UC_syntax(NULL) ;
299
300 (void) my_getenv("junk") ;
301
302 UC_read_opts( argc , argv ) ;
303 set_unusuality_tail( UC_ptail ) ;
304
305 oset = EDIT_empty_copy( UC_dset ) ;
306 EDIT_dset_items( oset ,
307 ADN_prefix , UC_prefix ,
308 ADN_ntt , 0 ,
309 ADN_nvals , 1 ,
310 ADN_datum_all , MRI_float ,
311 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
312 ADN_none ) ;
313
314 nvox = DSET_NVOX(oset) ;
315 far = (float *) malloc( sizeof(float) * nvox ) ;
316 for( kk=0 ; kk < nvox ; kk++ ) far[kk] = 0.0 ;
317 EDIT_substitute_brick( oset , 0 , MRI_float , far ) ;
318
319 if( !UC_be_quiet ){ printf("--- computing u") ; fflush(stdout) ; }
320
321 for( kk=0 ; kk < UC_nvec ; kk++ ){
322 ii = (UC_iv == NULL) ? kk : UC_iv[kk] ;
323 far[ii] = UC_unusuality( UC_vdim, UC_vec[kk] , UC_nvec, UC_vec ) ;
324 if( !UC_be_quiet && kk%1000==999 ){ printf(".");fflush(stdout); }
325 }
326 if( !UC_be_quiet ) printf("\n--- writing output\n") ;
327
328 DSET_write(oset) ;
329 exit(0) ;
330 }
331