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 /*---------------------------------------------------------------------------*/
8 /*
9 This program calculates a functional image from an AFNI 3D+time data set,
10 by calculating the correlation between the image time series and one or
11 more reference ("ideal") time series.
12
13 File: 3dfim.c
14 Date: 06 September 1996
15
16 Incorporated "ort" time series into 3dfim program.
17 BDW 12 December 1996
18
19 Added -percent option for calculating relative effect of reference waveform
20 upon observed time series.
21 BDW 19 May 1997
22
23 Corrected reference to "ort" time series data structure.
24 BDW 05 Sept 1997
25
26 Print a more explicit error message when ideal or "ort" time series are not
27 of sufficient length.
28 BDW 14 January 1998
29
30 Mod: Added changes for incorporating History notes.
31 Date: 10 September 1999
32
33 Mod: Added call to AFNI_logger.
34 Date: 15 August 2001
35
36 */
37
38 /*---------------------------------------------------------------------------*/
39
40 #define PROGRAM_NAME "3dfim" /* name of this program */
41 #define PROGRAM_AUTHOR "R. W. Cox and B. D. Ward" /* program author */
42 #define PROGRAM_INITIAL "06 Sept 1996" /* date of initial program release */
43 #define PROGRAM_LATEST "15 August 2001" /* date of latest program revision */
44
45 /*---------------------------------------------------------------------------*/
46
47 #define SO_BIG 33333
48
49 #include "mrilib.h"
50 #include "ts.h"
51 #include "afni_pcor.c"
52
53 /*-------------------------------------------------------------------
54 Lots of CPU work in here!
55 ---------------------------------------------------------------------*/
56
57 #undef FIM_THR
58 #define FIM_THR 0.0999
59
60
fim3d_fimmer_compute(THD_3dim_dataset * dset_time,time_series_array * ref_ts,time_series_array * ort_ts,int itbot,char * new_prefix,float max_percent)61 THD_3dim_dataset * fim3d_fimmer_compute ( THD_3dim_dataset * dset_time ,
62 time_series_array * ref_ts , time_series_array * ort_ts ,
63 int itbot, char * new_prefix,
64 float max_percent /* 19 May 1997 */ )
65 {
66 THD_3dim_dataset * new_dset ;
67 int ifim , it,iv , nvox=0 , ngood_ref , ntime , it1 , dtyp , nxyz;
68 float * vval , * tsar , * aval , * rbest , * abest ;
69 int * indx=NULL ;
70 short * bar ;
71 void * ptr ;
72 float stataux[MAX_STAT_AUX];
73 float fthr , topval ;
74 int nx_ref , ny_ref , ivec , nnow ;
75 PCOR_references ** pc_ref ;
76 PCOR_voxel_corr ** pc_vc ;
77 int save_resam ;
78
79 int fim_nref , nx_ort , ny_ort=0 , internal_ort ; /* 10 Dec 1996 */
80 static float * ref_vec = NULL ;
81 static int nref_vec = -666 ;
82
83 float * ref_ts_min = NULL,
84 * ref_ts_max = NULL,
85 * baseline = NULL; /* 19 May 1997 */
86
87 int i;
88
89 int nupdt = 0 , /* number of updates done yet */
90 min_updt = 5 ; /* min number needed for display */
91
92
93 /*--- check for legal inputs ---*/ /* 14 Jan 1998 */
94
95 if (!DSET_GRAPHABLE(dset_time))
96 {
97 fprintf (stderr, "Error: Invalid 3d+time input data file \n");
98 RETURN (NULL);
99 }
100
101 if (ref_ts == NULL)
102 {
103 fprintf (stderr, "Error: No ideal time series \n");
104 RETURN (NULL);
105 }
106
107 for (i = 0; i < ref_ts->num; i++)
108 if (ref_ts->tsarr[i]->len < DSET_NUM_TIMES(dset_time))
109 {
110 fprintf (stderr,
111 "Error: ideal time series is too short: ntime=%d num_ts=%d \n",
112 DSET_NUM_TIMES(dset_time),
113 ref_ts->tsarr[i]->len);
114 RETURN (NULL) ;
115 }
116
117
118 /** 10 Dec 1996: allow for orts **/
119
120 if( ort_ts->num > 0 ) /** 05 Sept 1997 **/
121 {
122 internal_ort = 0;
123 ny_ort = ort_ts->num;
124 for (i = 0; i < ny_ort; i++)
125 {
126 nx_ort = ort_ts->tsarr[i]->len ;
127 if (nx_ort < DSET_NUM_TIMES(dset_time)) /* 14 Jan 1998 */
128 {
129 fprintf (stderr,
130 "Error: ort time series is too short: ntime=%d ort_ts=%d \n",
131 DSET_NUM_TIMES(dset_time),
132 ort_ts->tsarr[i]->len);
133 RETURN (NULL) ;
134 }
135 }
136 }
137 else
138 {
139 internal_ort = 1 ;
140 }
141 fim_nref = (internal_ort) ? 3 : (ny_ort+3) ;
142
143 if( nref_vec < fim_nref )
144 {
145 ref_vec = (float *) malloc (sizeof(float)*fim_nref) ;
146 nref_vec = fim_nref;
147 }
148
149
150 /* arrays to store maximum change in the ideal time series */
151 if (max_percent > 0.0) /* 19 May 1997 */
152 {
153 ref_ts_max = (float *) malloc (sizeof(float) * (ref_ts->num));
154 ref_ts_min = (float *) malloc (sizeof(float) * (ref_ts->num));
155 }
156
157
158 nx_ref = ref_ts->tsarr[0]->len;
159 ny_ref = ref_ts->num;
160 ntime = DSET_NUM_TIMES(dset_time) ;
161 ngood_ref = 0 ;
162 it1 = -1 ;
163 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
164 tsar = ref_ts->tsarr[ivec]->ts;
165 ifim = 0 ;
166
167 if (max_percent > 0.0) /* 19 May 1997 */
168 {
169 ref_ts_min[ivec] = (float) SO_BIG;
170 ref_ts_max[ivec] = - (float) SO_BIG;
171 }
172
173 for( it=itbot ; it < ntime ; it++ )
174 {
175 if( tsar[it] < SO_BIG )
176 {
177 ifim++ ;
178 if( it1 < 0 ) it1 = it ;
179
180 if (max_percent > 0.0) /* 19 May 1997 */
181 {
182 if (tsar[it] > ref_ts_max[ivec]) ref_ts_max[ivec] = tsar[it];
183 if (tsar[it] < ref_ts_min[ivec]) ref_ts_min[ivec] = tsar[it];
184 }
185 }
186 }
187
188 if( ifim < min_updt ){
189 STATUS("ref_ts has too few good entries!") ;
190 RETURN(NULL) ;
191 }
192
193 ngood_ref = MAX( ifim , ngood_ref ) ;
194 }
195
196 /** at this point, ngood_ref = max number of good reference points,
197 and it1 = index of first point used in first reference **/
198
199 dtyp = DSET_BRICK_TYPE(dset_time,it1) ;
200 if( ! AFNI_GOOD_FUNC_DTYPE(dtyp) ){
201 STATUS("illegal input data type!") ;
202 RETURN(NULL) ;
203 }
204
205
206 #ifdef AFNI_DEBUG
207 { char str[256] ;
208 sprintf(str,"new prefix = %s",new_prefix) ; STATUS(str) ; }
209 #endif
210
211 /*--- FIM: find values above threshold to fim ---*/
212
213 DSET_load(dset_time); CHECK_LOAD_ERROR(dset_time);
214
215 nxyz = dset_time->dblk->diskptr->dimsizes[0]
216 * dset_time->dblk->diskptr->dimsizes[1]
217 * dset_time->dblk->diskptr->dimsizes[2] ;
218
219 /** find the mean of the first array,
220 compute the threshold (fthr) from it,
221 make indx[i] be the 3D index of the i-th voxel above threshold **/
222
223 switch( dtyp ){
224
225 case MRI_short:{
226 short * dar = (short *) DSET_ARRAY(dset_time,it1) ;
227 for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += abs(dar[iv]) ;
228 fthr = FIM_THR * fthr / nxyz ;
229 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
230 if( abs(dar[iv]) > fthr ) nvox++ ;
231 indx = (int *) malloc( sizeof(int) * nvox ) ;
232 if( indx == NULL ){
233 fprintf(stderr,"\n*** indx malloc failure in fim3d_fimmer_compute\n") ;
234 RETURN(NULL) ;
235 }
236 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
237 if( abs(dar[iv]) > fthr ) indx[nvox++] = iv ;
238 }
239 break ;
240
241 case MRI_float:{
242 float * dar = (float *) DSET_ARRAY(dset_time,it1) ;
243 for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += fabs(dar[iv]) ;
244 fthr = FIM_THR * fthr / nxyz ;
245 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
246 if( fabs(dar[iv]) > fthr ) nvox++ ;
247 indx = (int *) malloc( sizeof(int) * nvox ) ;
248 if( indx == NULL ){
249 fprintf(stderr,"\n*** indx malloc failure in fim3d_fimmer_compute\n") ;
250 RETURN(NULL) ;
251 }
252 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
253 if( fabs(dar[iv]) > fthr ) indx[nvox++] = iv ;
254 }
255 break ;
256
257 case MRI_byte:{
258 byte * dar = (byte *) DSET_ARRAY(dset_time,it1) ;
259 for( iv=0,fthr=0.0 ; iv < nxyz ; iv++ ) fthr += dar[iv] ;
260 fthr = FIM_THR * fthr / nxyz ;
261 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
262 if( dar[iv] > fthr ) nvox++ ;
263 indx = (int *) malloc( sizeof(int) * nvox ) ;
264 if( indx == NULL ){
265 fprintf(stderr,"\n*** indx malloc failure in fim3d_fimmer_compute\n") ;
266 RETURN(NULL) ;
267 }
268 for( iv=0,nvox=0 ; iv < nxyz ; iv++ )
269 if( dar[iv] > fthr ) indx[nvox++] = iv ;
270 }
271 break ;
272 }
273
274 /** allocate space for voxel values **/
275
276 vval = (float *) malloc( sizeof(float) * nvox) ;
277 if( vval == NULL ){
278 fprintf(stderr,"\n*** vval malloc failure in fim3d_fimmer_compute\n") ;
279 free(indx) ; RETURN(NULL) ;
280 }
281
282
283 /*----- allocate space for baseline values -----*/
284 if (max_percent > 0.0) /* 19 May 1997 */
285 {
286 baseline = (float *) malloc (sizeof(float) * nvox);
287 if (baseline == NULL)
288 {
289 fprintf(stderr,
290 "\n*** baseline malloc failure in fim3d_fimmer_compute\n") ;
291 free(indx) ; free(vval); RETURN(NULL) ;
292 }
293 else /* initialize baseline values to zero */
294 for (iv = 0; iv < nvox; iv++)
295 baseline[iv] = 0.0;
296 }
297
298
299 /** allocate extra space for comparing results from multiple ref vectors **/
300
301 if( ny_ref > 1 ){
302 aval = (float *) malloc( sizeof(float) * nvox) ;
303 rbest = (float *) malloc( sizeof(float) * nvox) ;
304 abest = (float *) malloc( sizeof(float) * nvox) ;
305 if( aval==NULL || rbest==NULL || abest==NULL ){
306 fprintf(stderr,"\n*** abest malloc failure in fim3d_fimmer_compute\n") ;
307 free(vval) ; free(indx) ;
308 if( aval != NULL ) free(aval) ;
309 if( rbest != NULL ) free(rbest) ;
310 if( abest != NULL ) free(abest) ;
311 RETURN(NULL) ;
312 }
313 } else {
314 aval = rbest = abest = NULL ;
315 }
316
317 #ifdef AFNI_DEBUG
318 { char str[256] ;
319 sprintf(str,"nxyz = %d nvox = %d",nxyz,nvox) ; STATUS(str) ; }
320 #endif
321
322 /*--- FIM: initialize recursive updates ---*/
323
324 pc_ref = (PCOR_references **) malloc( sizeof(PCOR_references *) * ny_ref ) ;
325 pc_vc = (PCOR_voxel_corr **) malloc( sizeof(PCOR_voxel_corr *) * ny_ref ) ;
326
327 if( pc_ref == NULL || pc_vc == NULL ){
328 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
329 if( aval != NULL ) free(aval) ;
330 if( rbest != NULL ) free(rbest) ;
331 if( abest != NULL ) free(abest) ;
332 fprintf(stderr,"\n*** FIM initialization fails in fim3d_fimmer_compute\n") ;
333 RETURN(NULL) ;
334 }
335
336 ifim = 0 ;
337 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
338 pc_ref[ivec] = new_PCOR_references( fim_nref ) ;
339 pc_vc[ivec] = new_PCOR_voxel_corr( nvox , fim_nref ) ;
340 if( pc_ref[ivec] == NULL || pc_vc[ivec] == NULL ) ifim++ ;
341 }
342
343 if( ifim > 0 ){
344 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
345 free_PCOR_references(pc_ref[ivec]) ;
346 free_PCOR_voxel_corr(pc_vc[ivec]) ;
347 }
348 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
349 if( aval != NULL ) free(aval) ;
350 if( rbest != NULL ) free(rbest) ;
351 if( abest != NULL ) free(abest) ;
352 fprintf(stderr,"\n*** FIM initialization fails in fim3d_fimmer_compute\n") ;
353 RETURN(NULL) ;
354 }
355
356 /*--- Make a new dataset to hold the output ---*/
357
358 new_dset = EDIT_empty_copy( dset_time ) ;
359
360 it = EDIT_dset_items( new_dset ,
361 ADN_prefix , new_prefix ,
362 ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
363 ADN_type , ISHEAD(dset_time)
364 ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
365 ADN_func_type , FUNC_COR_TYPE ,
366 ADN_nvals , FUNC_nvals[FUNC_COR_TYPE] ,
367 ADN_datum_all , MRI_short ,
368 ADN_ntt , 0 ,
369 ADN_none ) ;
370
371 if( it > 0 ){
372 fprintf(stderr,
373 "\n*** EDIT_dset_items error %d in fim3d_fimmer_compute\n",it) ;
374 THD_delete_3dim_dataset( new_dset , False ) ;
375 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
376 free_PCOR_references(pc_ref[ivec]) ;
377 free_PCOR_voxel_corr(pc_vc[ivec]) ;
378 }
379 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
380 if( aval != NULL ) free(aval) ;
381 if( rbest != NULL ) free(rbest) ;
382 if( abest != NULL ) free(abest) ;
383 RETURN(NULL) ;
384 }
385
386 for( iv=0 ; iv < new_dset->dblk->nvals ; iv++ ){
387 ptr = malloc( DSET_BRICK_BYTES(new_dset,iv) ) ;
388 mri_fix_data_pointer( ptr , DSET_BRICK(new_dset,iv) ) ;
389 }
390
391 if( THD_count_databricks(new_dset->dblk) < new_dset->dblk->nvals ){
392 fprintf(stderr,
393 "\n*** failure to malloc new bricks in fim3d_fimmer_compute\n") ;
394 THD_delete_3dim_dataset( new_dset , False ) ;
395 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
396 free_PCOR_references(pc_ref[ivec]) ;
397 free_PCOR_voxel_corr(pc_vc[ivec]) ;
398 }
399 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
400 if( aval != NULL ) free(aval) ;
401 if( rbest != NULL ) free(rbest) ;
402 if( abest != NULL ) free(abest) ;
403 RETURN(NULL) ;
404 }
405
406
407 /*--- FIM: do recursive updates ---*/
408
409 for( it=itbot ; it < ntime ; it++ ){
410
411 nnow = 0 ;
412 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
413 tsar = ref_ts->tsarr[ivec]->ts ;
414 if( tsar[it] >= SO_BIG ) continue ; /* skip this */
415
416 ref_vec[0] = 1.0 ; /* we always supply orts */
417 ref_vec[1] = (float) it ; /* for mean and linear trend */
418
419 if (internal_ort) /* 10 Dec 1996 */
420 {
421 ref_vec[2] = tsar[it] ;
422 }
423 else
424 {
425 for( iv=0 ; iv < ny_ort ; iv++ )
426 ref_vec[iv+2] = ort_ts->tsarr[iv]->ts[it];
427 ref_vec[ny_ort+2] = tsar[it] ;
428 }
429
430
431 #ifdef AFNI_DEBUG
432 { char str[256] ;
433 sprintf(str,"time index=%d ideal[%d]=%f" , it,ivec,tsar[it] ) ;
434 if (ivec == 0) STATUS(str) ; }
435 #endif
436
437
438 update_PCOR_references( ref_vec , pc_ref[ivec] ) ;
439
440 switch( dtyp ){
441 case MRI_short:{
442 short * dar = (short *) DSET_ARRAY(dset_time,it) ;
443 for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
444 }
445 break ;
446
447 case MRI_float:{
448 float * dar = (float *) DSET_ARRAY(dset_time,it) ;
449 for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
450 }
451 break ;
452
453 case MRI_byte:{
454 byte * dar = (byte *) DSET_ARRAY(dset_time,it) ;
455 for( iv=0 ; iv < nvox ; iv++ ) vval[iv] = (float) dar[indx[iv]] ;
456 }
457 break ;
458 }
459
460 PCOR_update_float( vval , pc_ref[ivec] , pc_vc[ivec] ) ;
461 nnow++ ;
462
463 /*----- update baseline value calculation -----*/
464 if (max_percent > 0.0) /* 19 May 1997 */
465 if (ivec == 0)
466 for (iv = 0; iv < nvox; iv++)
467 baseline[iv] += vval[iv] / ngood_ref;
468
469 }
470 if( nnow > 0 ) nupdt++ ;
471
472
473 /*--- Load results into the dataset and redisplay it ---*/
474
475 if( nupdt == ngood_ref )
476 {
477 /*--- set the statistical parameters ---*/
478
479 stataux[0] = nupdt ; /* number of points used */
480 stataux[1] = (ny_ref==1) ? 1 : 2 ; /* number of references */
481 stataux[2] = fim_nref - 1 ; /* number of orts */ /* 12 Dec 96 */
482 for( iv=3 ; iv < MAX_STAT_AUX ; iv++ ) stataux[iv] = 0.0 ;
483
484 STATUS("setting statistical parameters") ;
485
486 (void) EDIT_dset_items( new_dset ,
487 ADN_stat_aux , stataux ,
488 ADN_none ) ;
489
490 /*** Compute brick arrays for new dataset ***/
491
492 if( ny_ref == 1 ){
493
494 /*** Just 1 ref vector --> load values directly into dataset ***/
495
496 /*--- get alpha (coef) into vval,
497 find max value, scale into brick array ---*/
498
499 STATUS("getting 1 ref alpha") ;
500
501 PCOR_get_coef( pc_ref[0] , pc_vc[0] , vval ) ;
502
503 /*--- replace alpha with percentage change, if so requested ---*/
504 if (max_percent > 0.0) /* 19 May 1997 */
505 {
506 for (iv = 0; iv < nvox; iv++)
507 {
508 vval[iv] *= 100.0 * (ref_ts_max[0] - ref_ts_min[0]);
509 if (fabs(vval[iv]) < max_percent * fabs(baseline[iv]))
510 vval[iv] = fabs( vval[iv] / baseline[iv] );
511 else
512 vval[iv] = max_percent;
513 }
514 topval = max_percent;
515 }
516 else
517 {
518 topval = 0.0 ;
519 for( iv=0 ; iv < nvox ; iv++ )
520 if( fabs(vval[iv]) > topval ) topval = fabs(vval[iv]) ;
521 }
522
523 bar = DSET_ARRAY( new_dset , FUNC_ival_fim[FUNC_COR_TYPE] ) ;
524 memset( bar , 0 , sizeof(short)*nxyz ) ;
525
526 if( topval > 0.0 ){
527 topval = MRI_TYPE_maxval[MRI_short] / topval ;
528 for( iv=0 ; iv < nvox ; iv++ )
529 bar[indx[iv]] = (short)(topval * vval[iv] + 0.499) ;
530
531 stataux[0] = 1.0/topval ;
532 } else {
533 stataux[0] = 0.0 ;
534 }
535
536 /*--- get correlation coefficient (pcor) into vval,
537 scale into brick array (with fixed scaling factor) ---*/
538
539 STATUS("getting 1 ref pcor") ;
540
541 PCOR_get_pcor( pc_ref[0] , pc_vc[0] , vval ) ;
542
543 bar = DSET_ARRAY( new_dset , FUNC_ival_thr[FUNC_COR_TYPE] ) ;
544 memset( bar , 0 , sizeof(short)*nxyz ) ;
545
546 for( iv=0 ; iv < nvox ; iv++ )
547 bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * vval[iv] + 0.499) ;
548
549 stataux[1] = 1.0 / FUNC_COR_SCALE_SHORT ;
550
551 } else {
552
553 /*** Multiple references --> find best correlation at each voxel ***/
554
555 /*--- get first ref results into abest and rbest (best so far) ---*/
556
557 PCOR_get_coef( pc_ref[0] , pc_vc[0] , abest ) ;
558
559 /*--- modify alpha for percentage change calculation ---*/
560 if (max_percent > 0.0) /* 19 May 1997 */
561 for (iv = 0; iv < nvox; iv++)
562 abest[iv] *= 100.0 * (ref_ts_max[0] - ref_ts_min[0]);
563
564 PCOR_get_pcor( pc_ref[0] , pc_vc[0] , rbest ) ;
565
566 /*--- for each succeeding ref vector,
567 get results into aval and vval,
568 if |vval| > |rbest|, then use that result instead ---*/
569
570 for( ivec=1 ; ivec < ny_ref ; ivec++ ){
571
572 PCOR_get_coef( pc_ref[ivec] , pc_vc[ivec] , aval ) ;
573
574 PCOR_get_pcor( pc_ref[ivec] , pc_vc[ivec] , vval ) ;
575
576 for( iv=0 ; iv < nvox ; iv++ ){
577 if( fabs(vval[iv]) > fabs(rbest[iv]) ){
578 rbest[iv] = vval[iv] ;
579 abest[iv] = aval[iv] ;
580
581 /*--- modify alpha for percentage change calculation ---*/
582 if (max_percent > 0.0) /* 19 May 1997 */
583 abest[iv] *= 100.0 *
584 (ref_ts_max[ivec] - ref_ts_min[ivec]);
585
586 }
587 }
588
589 }
590
591 /*--- at this point, abest and rbest are the best
592 results, so scale them into the dataset bricks ---*/
593
594 /*--- finish percentage change calculation, if so requested ---*/
595 if (max_percent > 0.0) /* 19 May 1997 */
596 {
597 for (iv = 0; iv < nvox; iv++)
598 {
599 if (fabs(abest[iv]) < max_percent * fabs(baseline[iv]))
600 abest[iv] = fabs( abest[iv] / baseline[iv] );
601 else
602 abest[iv] = max_percent;
603 }
604 topval = max_percent;
605 }
606 else
607 {
608 topval = 0.0 ;
609 for( iv=0 ; iv < nvox ; iv++ )
610 if( fabs(abest[iv]) > topval ) topval = fabs(abest[iv]) ;
611 }
612
613 bar = DSET_ARRAY( new_dset , FUNC_ival_fim[FUNC_COR_TYPE] ) ;
614 memset( bar , 0 , sizeof(short)*nxyz ) ;
615
616 if( topval > 0.0 ){
617 topval = MRI_TYPE_maxval[MRI_short] / topval ;
618 for( iv=0 ; iv < nvox ; iv++ )
619 bar[indx[iv]] = (short)(topval * abest[iv] + 0.499) ;
620
621 stataux[0] = 1.0/topval ;
622 } else {
623 stataux[0] = 0.0 ;
624 }
625
626 bar = DSET_ARRAY( new_dset , FUNC_ival_thr[FUNC_COR_TYPE] ) ;
627 memset( bar , 0 , sizeof(short)*nxyz ) ;
628
629 for( iv=0 ; iv < nvox ; iv++ )
630 bar[indx[iv]] = (short)(FUNC_COR_SCALE_SHORT * rbest[iv] + 0.499) ;
631
632 stataux[1] = 1.0 / FUNC_COR_SCALE_SHORT ;
633
634 }
635
636 STATUS("setting brick_fac") ;
637
638 (void) EDIT_dset_items( new_dset ,
639 ADN_brick_fac , stataux ,
640 ADN_none ) ;
641
642 }
643 }
644
645
646 /*--- End of recursive updates; now free temporary workspaces ---*/
647
648 for( ivec=0 ; ivec < ny_ref ; ivec++ ){
649 free_PCOR_references(pc_ref[ivec]) ;
650 free_PCOR_voxel_corr(pc_vc[ivec]) ;
651 }
652 free(vval) ; free(indx) ; free(pc_ref) ; free(pc_vc) ;
653 if( aval != NULL ) free(aval) ;
654 if( rbest != NULL ) free(rbest) ;
655 if( abest != NULL ) free(abest) ;
656
657 if (ref_ts_min != NULL) free (ref_ts_min); /* 19 May 1997 */
658 if (ref_ts_max != NULL) free (ref_ts_max);
659 if (baseline != NULL) free (baseline);
660
661
662 /* --- load the statistics --- */
663 THD_load_statistics (new_dset);
664
665 /*--- Return new dataset ---*/
666
667 RETURN(new_dset) ;
668 }
669
670 /*--------------------------------------------------------------------------*/
671
672 /** structure type to hold results from scanning the command line options **/
673
674 typedef struct
675 {
676 char prefix_name[THD_MAX_NAME]; /* filename root for output */
677 THD_3dim_dataset * dset; /* input 3d+time data set */
678 time_series_array * idts; /* array of ideal time series */
679 time_series_array * ortts; /* array of ort time series */
680 int first_im; /* first time series point to be used */
681 float max_percent; /* max. percentage change due to ideal
682 time series, 19 May 1997 */
683 } line_opt ;
684
685
686 /*** internal prototype ***/
687
688 void Syntax( char * ) ;
689
690 /*--------------------------------------------------------------------------*/
691
692 /*** read and interpret command line arguments ***/
693
694
695 #define OPENERS "[{/%"
696 #define CLOSERS "]}/%"
697
698 #define IS_OPENER(sss) (strlen((sss))==1 && strstr(OPENERS,(sss))!=NULL)
699 #define IS_CLOSER(sss) (strlen((sss))==1 && strstr(CLOSERS,(sss))!=NULL)
700
701
get_line_opt(int argc,char * argv[],line_opt * opt)702 void get_line_opt( int argc , char *argv[] , line_opt *opt )
703 {
704
705 int nopt ;
706 float ftemp ;
707 MRI_IMAGE *im ;
708 time_series *ideal;
709 time_series *ort;
710 char err_note[128];
711 char input_filename[THD_MAX_NAME];
712
713
714 /* --- give help if needed --- */
715 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) Syntax(NULL) ;
716
717 /*----- add to program log -----*/
718 AFNI_logger (PROGRAM_NAME,argc,argv);
719
720 /* --- setup opt with defaults --- */
721 strcpy (opt->prefix_name, ""); /* no prefix name yet */
722 opt->first_im = 0 ; /* first image to actually use */
723 opt->max_percent = 0.0;
724
725 /* --- initialize array of time series data --- */
726 INIT_TSARR(opt->idts) ;
727
728 /* --- initialize array of ort time series data --- */ /* 10 Dec 1996 */
729 INIT_TSARR (opt->ortts);
730
731 /* --- set defaults in local variables --- */
732 ideal = NULL ; /* time_series of ideal response vector */
733 strcpy (input_filename, ""); /* no input file name yet */
734
735
736 /* --- read command line arguments and process them:
737 coding technique is to examine each argv, and if it matches
738 something we expect (e.g., -ideal), process it, then skip to
739 the next loop through ('continue' statements in each strncmp if).
740 Note that the 'while' will increment the argument index (nopt) --- */
741
742 nopt = 1 ;
743
744 do{
745
746 #ifdef DEBUG
747 # define DBOPT(str) fprintf(stderr,"at option %s: %s\n",argv[nopt],str)
748 #else
749 # define DBOPT(str) /* nothing */
750 #endif
751
752 /* --- -im1 # ==> index (1...) of first image to actually use --- */
753 if( strncmp(argv[nopt],"-im1",4) == 0 )
754 {
755 DBOPT("-im1") ;
756 if( ++nopt >= argc ) Syntax("-im1 needs an argument") ;
757 ftemp = strtod( argv[nopt] , NULL ) ;
758 if( ftemp < 1.0 )
759 {
760 sprintf( err_note , "illegal -im1 value: %f" , ftemp ) ;
761 Syntax(err_note) ;
762 }
763 opt->first_im = (int)(ftemp+0.499) - 1;
764 continue ;
765 }
766
767
768 /* --- -ideal file ==> ideal response vector time series --- */
769 if( strncmp(argv[nopt],"-ideal",5) == 0 )
770 {
771 DBOPT("-ideal") ;
772 if( ++nopt >= argc ) Syntax("-ideal needs a filename") ;
773
774 /** July 1995: read a bunch of ideals using [ a b c ... ] format **/
775
776 if( ! IS_OPENER(argv[nopt]) )
777 { /* --- one file --- */
778 ideal = RWC_read_time_series( argv[nopt] ) ;
779 if( ideal == NULL ) Syntax ("Unable to read ideal time series.");
780 ADDTO_TSARR( opt->idts , ideal ) ;
781 }
782 else
783 { /* --- many files --- */
784 for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ )
785 {
786 ideal = RWC_read_time_series( argv[nopt] ) ;
787 if( ideal == NULL ) Syntax ("Unable to read ideal time series.");
788 ADDTO_TSARR( opt->idts , ideal ) ;
789 }
790 if( nopt == argc ) Syntax("-ideal never finishes!") ;
791 }
792 continue ;
793 }
794
795
796 /*----- -percent p Calculate percentage change from baseline -----*/
797 if( strncmp(argv[nopt],"-percent",8) == 0 )
798 {
799 DBOPT("-percent") ;
800 if( ++nopt >= argc ) Syntax("-percent needs an argument") ;
801 ftemp = strtod( argv[nopt] , NULL ) ;
802 if( ftemp <= 0.0 )
803 {
804 sprintf( err_note , "illegal -percent value: %f" , ftemp ) ;
805 Syntax(err_note) ;
806 }
807 opt->max_percent = ftemp;
808 continue ;
809 }
810
811
812 /* --- -ort file ==> ort vector time series --- */ /* 10 Dec 1996 */
813 if( strncmp(argv[nopt],"-ort",4) == 0 )
814 {
815 DBOPT("-ort") ;
816 if( ++nopt >= argc ) Syntax("-ort needs a filename") ;
817
818 /** read a bunch of orts using [ a b c ... ] format **/
819
820 if( ! IS_OPENER(argv[nopt]) )
821 { /* --- one file --- */
822 ort = RWC_read_time_series( argv[nopt] ) ;
823 if( ort == NULL ) Syntax ("Unable to read ort time series.");
824 ADDTO_TSARR( opt->ortts , ort ) ;
825 }
826 else
827 { /* --- many files --- */
828 for( nopt++ ; !IS_CLOSER(argv[nopt]) && nopt<argc ; nopt++ )
829 {
830 ort = RWC_read_time_series( argv[nopt] ) ;
831 if( ort == NULL ) Syntax ("Unable to read ort time series.");
832 ADDTO_TSARR( opt->ortts , ort ) ;
833 }
834 if( nopt == argc ) Syntax("-ort never finishes!") ;
835 }
836 continue ;
837 }
838
839
840 /* --- -prefix name ==> prefix name of output file --- */
841 if( strncmp(argv[nopt], "-prefix", 5) == 0 ){
842 DBOPT("-prefix") ;
843 if( ++nopt >= argc ) Syntax("-prefix needs a name") ;
844 if( strcmp(opt->prefix_name, "") )
845 Syntax("Cannot have two prefix output names!" ) ;
846 strcpy (opt->prefix_name, argv[nopt]) ;
847 DBOPT("stored as prefix") ;
848 continue ;
849 }
850
851 /* --- -input fname ==> file name of input 3d+time data --- */
852
853 if( strncmp(argv[nopt], "-input", 5) == 0 )
854 {
855 DBOPT("-input") ;
856 if( ++nopt >= argc ) Syntax("-input needs a name") ;
857 if( strcmp(input_filename, "") )
858 Syntax ("Cannot have two input names!" ) ;
859 strcpy (input_filename, argv[nopt] );
860 /* open 3d data set */
861 opt->dset = THD_open_one_dataset( argv[nopt] ) ;
862 if( opt->dset == NULL )
863 {
864 fprintf(stderr, "\n** Unable to open 3d+time data file: %s\n",
865 argv[nopt]);
866 Syntax ("input dataset open failure");
867 }
868 continue ;
869 }
870
871 /* unidentified option */
872 sprintf( err_note , "Unknown option %s" , argv[nopt] ) ;
873 Syntax(err_note) ;
874
875 } while( ++nopt < argc ) ; /* loop until all args are read, or we break */
876
877
878 /* --- check for valid user inputs --- */
879 if (!strcmp(input_filename,"")) Syntax ("No input file name was given.");
880 if( opt->idts->num == 0 ) Syntax("No ideal response vector was defined!") ;
881 if (!strcmp(opt->prefix_name,"")) Syntax ("No prefix name was specified.");
882
883
884 /* --- Done! --- */
885
886 return ;
887 }
888
889 /*----------------------------------------------------------------------------*/
890
891 /* --- give some help and exit --- */
892
Syntax(char * note)893 void Syntax( char *note )
894 {
895 static char *mesg[] = {
896 "Program: 3dfim \n\n"
897 "Purpose: Calculate functional image from 3d+time data file. \n"
898 "Usage: 3dfim [-im1 num] -input fname -prefix name \n"
899 " -ideal fname [-ideal fname] [-ort fname] \n"
900 " " ,
901 "options are:",
902 "-im1 num num = index of first image to be used in time series ",
903 " correlation; default is 1 ",
904 " ",
905 "-input fname fname = filename of 3d + time data file for input",
906 " " ,
907 "-prefix name name = prefix of filename for saving functional data",
908 " ",
909 "-ideal fname fname = filename of a time series to which the image data",
910 " is to be correlated. ",
911 " ",
912 "-percent p Calculate percentage change due to the ideal time series ",
913 " p = maximum allowed percentage change from baseline ",
914 " Note: values greater than p are set equal to p. ",
915 " ",
916 "-ort fname fname = filename of a time series to which the image data",
917 " is to be orthogonalized ",
918 " ",
919 " N.B.: It is possible to specify more than",
920 " one ideal time series file. Each one is separately correlated",
921 " with the image time series and the one most highly correlated",
922 " is selected for each pixel. Multiple ideals are specified",
923 " using more than one '-ideal fname' option, or by using the",
924 " form '-ideal [ fname1 fname2 ... ]' -- this latter method",
925 " allows the use of wildcarded ideal filenames.",
926 " The '[' character that indicates the start of a group of",
927 " ideals can actually be any ONE of these: " OPENERS ,
928 " and the ']' that ends the group can be: " CLOSERS ,
929 " ",
930 " [Format of ideal time series files:",
931 " ASCII; one number per line;",
932 " Same number of lines as images in the time series;",
933 " Value over 33333 --> don't use this image in the analysis]",
934 " ",
935 " N.B.: It is also possible to specify more than",
936 " one ort time series file. The image time series is ",
937 " orthogonalized to each ort time series. Multiple orts are ",
938 " specified by using more than one '-ort fname' option, ",
939 " or by using the form '-ort [ fname1 fname2 ... ]'. This ",
940 " latter method allows the use of wildcarded ort filenames.",
941 " The '[' character that indicates the start of a group of",
942 " ideals can actually be any ONE of these: " OPENERS ,
943 " and the ']' that ends the group can be: " CLOSERS ,
944 " ",
945 " [Format of ort time series files:",
946 " ASCII; one number per line;",
947 " At least same number of lines as images in the time series]",
948 " ",
949 " ",
950 } ;
951
952 int nsize , ii ;
953
954 if( note != NULL && (nsize=strlen(note)) > 0 ){
955 fprintf(stderr,"\n") ;
956 for(ii=0;ii<nsize+9;ii++) fprintf(stderr,"*") ;
957 fprintf(stderr,"\a\n Error: %s\n", note ) ;
958 for(ii=0;ii<nsize+9;ii++) fprintf(stderr,"*") ;
959 fprintf(stderr,"\nTry 3dfim -help for instructions\a\n") ;
960 exit(-1) ;
961 }
962
963 for( ii=0 ; ii < sizeof(mesg)/sizeof(char *) ; ii++ ){
964 printf( " %s\n" , mesg[ii] );
965 }
966 PRINT_COMPILE_DATE ; exit(0) ;
967 }
968
969
970 /*----------------------------------------------------------------------------*/
971
main(int argc,char * argv[])972 int main( int argc , char *argv[] )
973 {
974 line_opt opt ; /* holds data constructed from command line */
975 THD_3dim_dataset * new_dset; /* functional data set to be calculated */
976 RwcBoolean ok; /* true if 3d write is successful */
977
978
979 /*----- Identify software -----*/
980 #if 0
981 printf ("\n\n");
982 printf ("Program: %s \n", PROGRAM_NAME);
983 printf ("Author: %s \n", PROGRAM_AUTHOR);
984 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
985 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
986 printf ("\n");
987 #endif
988
989 PRINT_VERSION("3dfim") ; AUTHOR(PROGRAM_AUTHOR) ;
990 mainENTRY("3dfim main") ; machdep() ;
991
992 WARNING_message("This program (3dfim) is very old, may not be useful, and will not be maintained.") ;
993
994 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
995
996 { int new_argc ; char ** new_argv ;
997 addto_args( argc , argv , &new_argc , &new_argv ) ;
998 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
999 }
1000
1001
1002 /* --- read command line --- */
1003 get_line_opt( argc , argv , &opt ) ;
1004
1005 /* --- calculate functional image --- */
1006 new_dset = fim3d_fimmer_compute (opt.dset, opt.idts, opt.ortts,
1007 opt.first_im, opt.prefix_name,
1008 opt.max_percent); /* 19 May 1997 */
1009
1010 /*----- Record history of dataset -----*/
1011 tross_Copy_History( opt.dset , new_dset ) ;
1012 tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
1013
1014 /* --- write 3d functional image data --- */
1015 ok = THD_write_3dim_dataset ("", opt.prefix_name, new_dset, TRUE);
1016 if (!ok) Syntax ("Failure to write functional data set.");
1017
1018 /* --- cleanup --- */
1019 THD_delete_3dim_dataset (new_dset, FALSE);
1020
1021 return (0);
1022 }
1023
1024 /*---------------------------------------------------------------------------*/
1025