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 <stdlib.h>
8 #include <stdio.h>
9 #include <string.h>
10 #include <math.h>
11 #include <ctype.h>
12
13 #include "mrilib.h"
14
15 #define MAX_NAME 64
16
17 typedef struct {
18 int count ;
19 char name[MAX_NAME] ;
20 MRI_IMAGE * avim , * sdim ;
21 } SF_interval ;
22
23 /*** global data!!! ***/
24
25 static int SF_numint = 0 ;
26 static SF_interval * SF_int = NULL ;
27 static char SF_bname[MAX_NAME] = "rest" ;
28 static char SF_pname[MAX_NAME] = "sfim." ;
29 static char SF_iname[MAX_NAME] = "sfint" ;
30 static MRI_IMARR * SF_imts = NULL ;
31 static int SF_localbase = 0 ;
32
33 /*** prototypes ***/
34
35 void SFIM_getopts( int , char * argv[] ) ;
36 SF_interval * SFIM_load_intervals( char * ) ;
37 void SFIM_syntax( char * ) ;
38 void SFIM_write_avs() ;
39
40 /**********************************************************************/
41
main(int argc,char * argv[])42 int main( int argc , char * argv[] )
43 {
44 int lin , kim , kbot,ktop , nx,ny , npix , ii ,
45 lbase , lup,ldown ;
46 MRI_IMAGE ** stat_ret ;
47 MRI_IMAGE * imb=NULL ;
48 float * bar , * bav ;
49
50 printf(
51 "MCW SFIM: Stepwise Functional IMages, by RW Cox\n") ;
52
53 if( argc < 2 ) SFIM_syntax("type sfim -help for usage details") ;
54 else if( strcmp(argv[1],"-help") == 0 ) SFIM_syntax(NULL) ;
55
56 machdep() ;
57
58 SFIM_getopts( argc , argv ) ;
59
60 /*----- average over each interval -----*/
61
62 nx = SF_imts->imarr[0]->nx ;
63 ny = SF_imts->imarr[0]->ny ; npix = nx * ny ;
64
65 lin = 0 ; kbot = 0 ;
66 do {
67 ktop = kbot + SF_int[lin].count ;
68 if( ktop > SF_imts->num ) ktop = SF_imts->num ;
69
70 if( isalpha(SF_int[lin].name[0]) ){ /* average if a good name */
71 for( kim=kbot ; kim < ktop ; kim++ )
72 (void) mri_stat_seq( SF_imts->imarr[kim] ) ;
73
74 stat_ret = mri_stat_seq( NULL ) ;
75 SF_int[lin].avim = stat_ret[0] ;
76 SF_int[lin].sdim = stat_ret[1] ;
77 }
78
79 kbot = ktop ; lin ++ ;
80 } while( SF_int[lin].count > 0 && kbot < SF_imts->num ) ;
81 SF_numint = lin ;
82
83 /*----- find the number of base intervals -----*/
84
85 lbase = 0 ;
86 for( lin=0 ; lin < SF_numint ; lin++ )
87 if( strcmp(SF_int[lin].name,SF_bname) == 0 ) lbase++ ;
88
89 /* no bases --> write averages out now and quit */
90
91 if( lbase <= 0 ){
92 printf("** no 'base' intervals --> task means not adjusted\n") ;
93 SFIM_write_avs() ;
94 exit(0) ;
95 }
96
97 /* bases yes, but not localbase --> compute global average of bases */
98
99 if( ! SF_localbase ){
100 int knum = 0 ;
101
102 kbot = 0 ;
103 for( lin=0 ; lin < SF_numint ; lin++ ){
104 ktop = kbot + SF_int[lin].count ;
105 if( ktop > SF_imts->num ) ktop = SF_imts->num ;
106 if( strcmp(SF_int[lin].name,SF_bname) == 0 ){ /* if a base */
107 for( kim=kbot ; kim < ktop ; kim++ ){
108 (void) mri_stat_seq( SF_imts->imarr[kim] ) ; /* average in */
109 knum ++ ;
110 }
111 }
112 kbot = ktop ;
113 }
114 stat_ret = mri_stat_seq( NULL ) ;
115 imb = stat_ret[0] ; /* average of all bases */
116 mri_free( stat_ret[1] ) ; /* don't keep st. dev. */
117
118 printf("** global base = average of %d images\n",knum) ;
119 }
120
121 /*----- for each non-base interval,
122 subtract the relevant base average -----*/
123
124 for( lin=0 ; lin < SF_numint ; lin++ ){
125 int free_imb = 0 ;
126
127 if( !isalpha(SF_int[lin].name[0]) ||
128 strcmp(SF_int[lin].name,SF_bname) == 0 ) continue ; /* skip this */
129
130 if( SF_localbase ){
131 for( lup=lin+1 ; lup < SF_numint ; lup++ ) /* look for a base above */
132 if( strcmp(SF_int[lup].name,SF_bname) == 0 ) break ;
133
134 for( ldown=lin-1 ; ldown >=0 ; ldown-- ) /* look for a base below */
135 if( strcmp(SF_int[ldown].name,SF_bname) == 0 ) break ;
136
137 if( ldown < 0 && lup >= SF_numint ){ /* no base? an error! */
138 fprintf(stderr,"*** can't find base above or below at lin=%d\n",lin) ;
139 SFIM_syntax("INTERNAL ERROR -- should not occur!") ;
140 }
141
142 /* if only have one neighbor, use it, otherwise make average */
143
144 if( ldown < 0 ){
145 imb = SF_int[lup].avim ; free_imb = 0 ;
146 printf("** local base for %s = average of %d images above\n",
147 SF_int[lin].name , SF_int[lup].count ) ;
148 }
149 else if( lup >= SF_numint ){
150 imb = SF_int[ldown].avim ; free_imb = 0 ;
151 printf("** local base for %s = average of %d images below\n",
152 SF_int[lin].name , SF_int[ldown].count ) ;
153 }
154 else {
155 float * bup , * bdown ;
156 bup = mri_data_pointer( SF_int[lup].avim ) ;
157 bdown = mri_data_pointer( SF_int[ldown].avim ) ;
158 imb = mri_new( nx , ny , MRI_float ) ; free_imb = 1 ;
159 bar = mri_data_pointer( imb ) ;
160 for( ii=0 ; ii < npix ; ii++ )
161 bar[ii] = 0.5 * ( bup[ii] + bdown[ii] ) ;
162
163 printf("** local base for %s = average of %d below, %d above\n",
164 SF_int[lin].name, SF_int[ldown].count, SF_int[lup].count ) ;
165 }
166 }
167
168 /* subtract imb (base average) from current interval average */
169
170 bar = mri_data_pointer( imb ) ;
171 bav = mri_data_pointer( SF_int[lin].avim ) ;
172 for( ii=0 ; ii < npix ; ii++ ) bav[ii] -= bar[ii] ;
173
174 if( SF_localbase && free_imb ) mri_free( imb ) ;
175 }
176
177 /*----- now write the averages out -----*/
178
179 SFIM_write_avs() ;
180 exit(0) ;
181 }
182
183 /**********************************************************************/
184
SFIM_write_avs()185 void SFIM_write_avs()
186 {
187 char name[256] ;
188 int lin , ldown , tnum ;
189
190 for( lin=0 ; lin < SF_numint ; lin++ ){
191 if( !isalpha(SF_int[lin].name[0]) ||
192 SF_int[lin].avim == NULL ) continue ; /* skip */
193
194 tnum = 1 ;
195 for( ldown=lin-1 ; ldown >=0 ; ldown-- )
196 if( strcmp(SF_int[lin].name,SF_int[ldown].name) == 0 ) tnum++ ;
197
198 sprintf( name , "%s%s.%04d" , SF_pname , SF_int[lin].name , tnum ) ;
199
200 printf("-- writing file %s\n",name) ;
201 mri_write( name , SF_int[lin].avim ) ;
202 }
203 }
204
205 /**********************************************************************/
206
SFIM_getopts(int argc,char * argv[])207 void SFIM_getopts( int argc , char * argv[] )
208 {
209 int nopt = 1 ;
210 int kk , nx , ny , kount ;
211
212 /*--- do options ---*/
213
214 while( nopt < argc && argv[nopt][0] == '-' ){
215
216 /** -sfint iname **/
217
218 if( strncmp(argv[nopt],"-sfint",5) == 0 ){
219 if( ++nopt >= argc ) SFIM_syntax("-sfint needs a name after it!") ;
220 strcpy(SF_iname,argv[nopt]) ;
221 nopt++ ; continue ;
222 }
223
224 /** -base bname **/
225
226 if( strncmp(argv[nopt],"-base",5) == 0 ){
227 if( ++nopt >= argc ) SFIM_syntax("-base needs a name after it!") ;
228 strcpy(SF_bname,argv[nopt]) ;
229 nopt++ ; continue ;
230 }
231
232 /** -prefix pname **/
233
234 if( strncmp(argv[nopt],"-prefix",5) == 0 ){
235 if( ++nopt >= argc ) SFIM_syntax("-prefix needs a name after it!") ;
236 strcpy(SF_pname,argv[nopt]) ;
237 kk = strlen(SF_pname) ;
238 if( SF_pname[kk-1] != '.' ){
239 SF_pname[kk] = '.' ;
240 SF_pname[kk+1] = '\0' ;
241 }
242 nopt++ ; continue ;
243 }
244
245 /** -localbase **/
246
247 if( strncmp(argv[nopt],"-localbase",5) == 0 ){
248 SF_localbase = 1 ;
249 nopt++ ; continue ;
250 }
251
252 /** illegal option **/
253
254 fprintf(stderr,"*** illegal command line option: %s\n",argv[nopt]) ;
255 SFIM_syntax("type sfim -help for more details") ;
256
257 }
258
259 /*--- do images ---*/
260
261 SF_imts = mri_read_many_files( argc-nopt , argv+nopt ) ;
262 if( SF_imts == NULL ) SFIM_syntax("cannot continue without input images!") ;
263
264 /* check images for consistency */
265
266 nx = SF_imts->imarr[0]->nx ;
267 ny = SF_imts->imarr[0]->ny ;
268
269 for( kk=1 ; kk < SF_imts->num ; kk++ ){
270 if( nx != SF_imts->imarr[kk]->nx || ny != SF_imts->imarr[kk]->ny ){
271 fprintf(stderr,"*** image %d not conformant to image 0\n",kk) ;
272 SFIM_syntax("cannot continue with images whose sizes differ!") ;
273 }
274 }
275
276 /*--- read intervals ---*/
277
278 SF_int = SFIM_load_intervals( SF_iname ) ;
279 if( SF_int[0].count <= 0 ) SFIM_syntax("cannot read interval definitions") ;
280
281 kount = kk = 0 ;
282 do {
283 kount += SF_int[kk++].count ;
284 } while( SF_int[kk].count > 0 ) ;
285
286 if( kount < SF_imts->num ){
287 fprintf(stderr,"*** # images = %d but only %d images in intervals!\n",
288 SF_imts->num , kount ) ;
289 SFIM_syntax("must have at least as many in intervals as actual images!") ;
290 }
291
292 return ;
293 }
294
295 /**********************************************************************/
296
SFIM_syntax(char * str)297 void SFIM_syntax( char * str )
298 {
299 if( str != NULL ){
300 fprintf(stderr,"*** %s\n",str) ;
301 exit(-1) ;
302 }
303
304 printf(
305 "\n"
306 "Usage: sfim [options] image_files ...\n"
307 "\n"
308 " + image_files are in the same format AFNI accepts\n"
309 " + options are from the following:\n"
310 "\n"
311 " -sfint iname: 'iname' is the name of a file which has\n"
312 " the interval definitions; an example is\n"
313 " 3*# 5*rest 4*A 5*rest 4*B 5*rest 4*A 5*rest\n"
314 " which says:\n"
315 " - ignore the 1st 3 images\n"
316 " - take the next 5 as being in task state 'rest'\n"
317 " - take the next 4 as being in task state 'A'\n"
318 " and so on;\n"
319 " task names that start with a nonalphabetic character\n"
320 " are like the '#' above and mean 'ignore'.\n"
321 " *** the default 'iname' is 'sfint'\n"
322 "\n"
323 " -base bname: 'bname' is the task state name to use as the\n"
324 " baseline; other task states will have the mean\n"
325 " baseline state subtracted; if there are no task\n"
326 " states from 'iname' that match 'bname', this\n"
327 " subtraction will not occur.\n"
328 " *** the default 'bname' is 'rest'\n"
329 "\n"
330 " -localbase: if this option is present, then each non-base\n"
331 " task state interval has the mean of the two\n"
332 " nearest base intervals subtracted instead of the\n"
333 " grand mean of all the base task intervals.\n"
334 "\n"
335 " -prefix pname: 'pname' is the prefix for output image filenames for\n"
336 " all states: the i'th interval with task state name\n"
337 " 'fred' will be writen to file 'pname.fred.i'.\n"
338 " *** the default 'pname' is 'sfim'\n"
339 "\n"
340 " Output files are the base-mean-removed averages for each non-base\n"
341 " task interval, and simply the mean for each base task interval.\n"
342 " Output images are in the 'flim' (floating pt. image) format, and\n"
343 " may be converted to 16 bit shorts using the program 'ftosh'.\n"
344 ) ;
345
346 exit(0) ;
347 }
348
349 /**********************************************************************/
350
SFIM_load_intervals(char * fname)351 SF_interval * SFIM_load_intervals( char * fname )
352 {
353 SF_interval * sfint ;
354 FILE * fd ;
355 int num_int , num_all , ic , cc ;
356 char nn[MAX_NAME] ;
357 #define INC_SFINT 8
358
359 if( fname == NULL || strlen(fname) == 0 ){
360 fprintf(stderr,"SFIM_load_intervals illegal filename\n") ;
361 exit(-1) ;
362 }
363 fd = fopen( fname , "r" ) ;
364 if( fd == NULL ){
365 fprintf(stderr,"SFIM_load_intervals cannot open filename %s\n",fname) ;
366 exit(-1) ;
367 }
368
369 /* malloc -> calloc 13 Feb 2009 [lesstif patrol] */
370 sfint = (SF_interval *) calloc( INC_SFINT, sizeof(SF_interval) ) ;
371 if( sfint == NULL ){
372 fprintf(stderr,"SFIM_load_intervals fails to malloc!\n") ;
373 exit(-1) ;
374 }
375 num_int = 0 ;
376 num_all = INC_SFINT ;
377
378 do {
379 ic = fscanf( fd , "%d*%s" , &cc , nn ) ;
380 if( num_int == num_all ){
381 num_all += INC_SFINT ;
382 sfint = (SF_interval *) realloc( (void *) sfint ,
383 sizeof(SF_interval) * num_all ) ;
384 if( sfint == NULL ){
385 fprintf(stderr,"SFIM_load_intervals fails to realloc!\n") ;
386 exit(-1) ;
387 }
388 }
389 if( ic != 2 ){ /* bad scanning */
390 sfint[num_int].count = -1 ;
391 break ;
392 }
393 sfint[num_int].count = cc ; /* put results into interval */
394 sfint[num_int].avim = NULL ;
395 sfint[num_int].sdim = NULL ;
396 strcpy( sfint[num_int].name , nn ) ;
397
398 num_int ++ ;
399 } while( cc > 0 ) ;
400
401 fclose( fd ) ;
402 return sfint ;
403 }
404