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 "mrilib.h"
8
9 /*---------------------------------------------------------------------------
10 This program catenates multiple 3D+time datasets to create one
11 super-dataset. Adapted from 3dbucket.c -- RWCox 21 Sep 1998.
12 ----------------------------------------------------------------------------*/
13
14 /*-------------------------- global data --------------------------*/
15
16 static THD_3dim_dataset_array *TCAT_dsar = NULL ; /* input datasets */
17 static RwcPointer_array *TCAT_subv = NULL ; /* sub-brick selectors */
18 static int TCAT_nvox = -1 ; /* # voxels */
19 static int TCAT_dry = 0 ; /* dry run? */
20 static int TCAT_verb = 0 ; /* verbose? */
21 static int TCAT_type = -1 ; /* dataset type */
22 static int TCAT_glue = 0 ; /* glueing run? */
23 static int TCAT_ccode = COMPRESS_NONE ; /* 16 Mar 2010 */
24 static int TCAT_rlt = 0 ; /* remove linear trend? */
25
26 static int TCAT_rqt = 0 ; /* 15 Nov 1999 */
27 static int TCAT_rct = 0 ; /* 15 Nov 1999 */
28
29 static int TCAT_relabel = 0 ; /* 03 Nov 2010 */
30
31 static char * TCAT_tpattern = NULL;/* 08 Mar 2013 [rickr] */
32 static float TCAT_tr = 0.0 ; /* 08 Mar 2013 [rickr] */
33
34 static char TCAT_output_prefix[THD_MAX_PREFIX] = "tcat" ;
35 static char TCAT_session[THD_MAX_NAME] = "./" ;
36
37 /* macros to get
38 DSUB = a particular input dataset
39 NSUBV = the number of sub-bricks selected from a dataset
40 SUBV = a particular sub-brick selected from a dataset */
41
42 #define DSUB(id) DSET_IN_3DARR(TCAT_dsar,(id))
43 #define NSUBV(id) ( ((int *)TCAT_subv->ar[(id)])[0] )
44 #define SUBV(id,jj) ( ((int *)TCAT_subv->ar[(id)])[(jj)+1] )
45
46 /*--------------------------- prototypes ---------------------------*/
47
48 void TCAT_read_opts( int , char ** ) ;
49 void TCAT_Syntax ( void ) ;
50 int * TCAT_get_subv ( int , char * ) ;
51
52 /*--------------------------------------------------------------------
53 read the arguments, load the global variables
54 ----------------------------------------------------------------------*/
55
TCAT_read_opts(int argc,char * argv[])56 void TCAT_read_opts( int argc , char *argv[] )
57 {
58 int nopt = 1 , ii ;
59 char dname[THD_MAX_NAME] ;
60 char *subv=NULL ; /* no sub-brick length limit 17 Jun 2010 [rickr] */
61 char *cpt ;
62 THD_3dim_dataset *dset, *fset=NULL ;
63 int *svar ;
64 char *str;
65 int ok, ilen=0, nlen , max_nsub=0 ;
66
67 INIT_3DARR(TCAT_dsar) ; /* array of datasets */
68 INIT_XTARR(TCAT_subv) ; /* array of sub-brick selector arrays */
69
70 PUTENV("AFNI_GLOB_SELECTORS","YES") ; /* 09 Mar 2011 */
71
72 while( nopt < argc ){
73
74 /**** -relabel ****/
75
76 if( strcmp(argv[nopt],"-relabel") == 0 ){ /* 03 Nov 2010 */
77 TCAT_relabel++ ; nopt++ ; continue ;
78 }
79
80 /**** -prefix prefix ****/
81
82 if( strncmp(argv[nopt],"-prefix",6) == 0 ||
83 strncmp(argv[nopt],"-output",6) == 0 ){
84 if(TCAT_glue)
85 ERROR_exit("-prefix and -glueto options are not compatible");
86 nopt++ ;
87 if( nopt >= argc )
88 ERROR_exit("need argument after -prefix!") ;
89 MCW_strncpy( TCAT_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
90 continue ;
91 }
92
93 /**** -tpattern PATTERN ****/
94
95 if( strncmp(argv[nopt],"-tpat",5) == 0 ) {
96 nopt++ ;
97 if( nopt >= argc )
98 ERROR_exit("need argument after -tpattern!") ;
99 TCAT_tpattern = argv[nopt] ;
100 /* allow @tpattern 02 Sep 2014 [rickr] */
101 nopt++ ; continue ;
102 }
103
104 /**** -tr TR ****/
105
106 if( strcmp(argv[nopt],"-tr") == 0 || strcmp(argv[nopt],"-TR") == 0 ){
107 nopt++ ;
108 if( nopt >= argc )
109 ERROR_exit("need argument after -tr!") ;
110 TCAT_tr = atof(argv[nopt]) ;
111 if( TCAT_tr <= 0.0 ) ERROR_exit("illegal -tr value: %s", argv[nopt]);
112 if( TCAT_tr >= 100.0 )
113 WARNING_message("-tr is in seconds, so %g seems big", TCAT_tr);
114 nopt++ ; continue ;
115 }
116
117 /**** -session directory ****/
118
119 if( strncmp(argv[nopt],"-session",6) == 0 ){
120 if(TCAT_glue)
121 ERROR_exit("-session and -glueto options are not compatible");
122 nopt++ ;
123 if( nopt >= argc )
124 ERROR_exit("need argument after -session!") ;
125 MCW_strncpy( TCAT_session , argv[nopt++] , THD_MAX_NAME ) ;
126 continue ;
127 }
128
129 /**** -dry ****/
130
131 if( strncmp(argv[nopt],"-dry",3) == 0 ){
132 TCAT_dry = 1 ; TCAT_verb++ ;
133 nopt++ ; continue ;
134 }
135
136 /**** -verb ****/
137
138 if( strncmp(argv[nopt],"-verb",5) == 0 ){
139 TCAT_verb++ ;
140 nopt++ ; continue ;
141 }
142
143 /**** -rlt ****/
144
145 if( strcmp(argv[nopt],"-rlt") == 0 ){
146 TCAT_rlt = 1 ;
147 nopt++ ; continue ;
148 }
149
150 /**** -rlt+ [16 Sep 1999] ****/
151
152 if( strcmp(argv[nopt],"-rlt+") == 0 ){ /* 16 Sep 1999 */
153 TCAT_rlt = 2 ;
154 nopt++ ; continue ;
155 }
156
157 /**** -rlt++ [16 Sep 1999] ****/
158
159 if( strcmp(argv[nopt],"-rlt++") == 0 ){ /* 16 Sep 1999 */
160 TCAT_rlt = 3 ;
161 nopt++ ; continue ;
162 }
163
164 /**** -rqt [15 Nov 1999] ****/
165
166 if( strcmp(argv[nopt],"-rqt") == 0 ){
167 TCAT_rqt = 1 ;
168 nopt++ ; continue ;
169 }
170
171 /**** -rct [15 Nov 1999] ****/
172
173 if( strcmp(argv[nopt],"-rct") == 0 ){
174 TCAT_rct = 1 ;
175 nopt++ ; continue ;
176 }
177
178 /**** -glueto fname ****/
179
180 if( strncmp(argv[nopt],"-glueto",5) == 0 ){
181 if( strncmp(TCAT_output_prefix, "tcat", 5) != 0 )
182 ERROR_exit("-prefix and -glueto options are not compatible");
183 if( strncmp(TCAT_session, "./", 5) != 0 )
184 ERROR_exit("-session and -glueto options are not compatible");
185 TCAT_glue = 1 ;
186 nopt++ ;
187 if( nopt >= argc )
188 ERROR_exit("need argument after -glueto!") ;
189
190 /*----- Verify that file name ends in View Type -----*/
191 ok = 1;
192 nlen = strlen(argv[nopt]);
193 if (nlen <= 5) ok = 0;
194
195 #define BACKASS /* 03 Oct 2002 -- RWCox */
196
197 if (ok)
198 {
199 #ifndef BACKASS
200 for (ilen = 0; ilen < nlen; ilen++) /* BDW: scan forward */
201 #else
202 for( ilen=nlen-3 ; ilen >= 0 ; ilen-- ) /* RWC: scan backward */
203 #endif
204 {
205 str = argv[nopt] + ilen;
206 if (str[0] == '+') break;
207 }
208 #ifndef BACKASS
209 if (ilen == nlen) ok = 0;
210 #else
211 if (ilen <= 0 ) ok = 0;
212 #endif
213 }
214
215 if (ok)
216 {
217 str = argv[nopt] + ilen + 1;
218
219 for (ii=FIRST_VIEW_TYPE ; ii <= LAST_VIEW_TYPE ; ii++)
220 if (! strncmp(str,VIEW_codestr[ii],4)) break ;
221
222 if( ii > LAST_VIEW_TYPE ) ok = 0;
223 }
224
225 if (! ok)
226 ERROR_exit(
227 "File name must end in +orig, +acpc, or +tlrc after -glueto");
228
229 /*----- Remove View Type from string to make output prefix -----*/
230 MCW_strncpy( TCAT_output_prefix , argv[nopt] , ilen+1) ;
231
232 /*----- Note: no "continue" statement here. File name
233 will now be processed as an input dataset -----*/
234
235 } /* end of -glueto prefatory actions */
236
237 /*-- does this look like another option? that's bad! --*/
238
239 if( argv[nopt][0] == '-' ) ERROR_exit("Unknown option: %s",argv[nopt]) ;
240
241 /*-- 09 Mar 2011: do filename expansion --*/
242
243 { int nexp,ee,ebad=0 ; char **fexp ;
244
245 MCW_file_expand( 1 , argv+nopt , &nexp , &fexp ) ;
246
247 if( nexp == 0 ){ ebad = nexp = 1 ; fexp = argv+nopt ; }
248 nopt++ ;
249
250 if( TCAT_verb )
251 INFO_message("3dTcat, file_expand nexp = %d, fexp[0] = %s\n",
252 nexp, *fexp);
253
254 for( ee=0 ; ee < nexp ; ee++ ){
255
256 /**** read dataset ****/
257
258 cpt = strstr(fexp[ee],"[") ; /* look for the sub-brick selector */
259
260 subv = NULL; /* Need to resest this variable ZSS Nov. 24 2010 */
261 if( cpt == NULL ){ /* no selector */
262 strcpy(dname,fexp[ee]) ;
263 } else if( cpt == fexp[ee] ){ /* can't be at start!*/
264 ERROR_exit("Illegal dataset specifier: %s",fexp[ee]) ;
265 } else { /* found selector */
266 ii = cpt - fexp[ee] ;
267 memcpy(dname,fexp[ee],ii) ; dname[ii] = '\0' ;
268 subv = cpt; /* no length limit 17 Jun 2010 [rickr] */
269 }
270
271 if( TCAT_verb > 1 )
272 INFO_message("opening one tcat dset #%d, %s", ee, dname);
273
274 dset = THD_open_one_dataset( dname ) ;
275 /* rather than failing, try new-fangled open 4 Apr 2016 [rickr] */
276 /* NOTE: let THD_open_dataset parse sub-brick selectors */
277 if( dset == NULL ) {
278 if( TCAT_verb > 1 )
279 INFO_message("opening tcat dset #%d, %s", ee, fexp[ee]);
280
281 dset = THD_open_dataset( fexp[ee] ) ;
282 subv = NULL;
283 }
284 if( dset == NULL ) ERROR_exit("Can't open dataset %s",dname) ;
285 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
286
287 if( TCAT_type < 0 ) TCAT_type = dset->type ;
288
289 /* 16 Mar 2010 */
290 TCAT_ccode = COMPRESS_filecode(dset->dblk->diskptr->brick_name) ;
291
292 /* check if voxel counts match */
293
294 ii = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
295 if( TCAT_nvox < 0 ){
296 TCAT_nvox = ii ; fset = dset ;
297 } else if( ii != TCAT_nvox ){
298 ERROR_exit("Dataset %s differs in size from first one!",dname);
299 } else if( !EQUIV_GRIDS(dset,fset) ){
300 WARNING_message("Dataset %s grid differs from first one!",dname);
301 }
302 ADDTO_3DARR(TCAT_dsar,dset) ; /* list of datasets */
303
304 if (subv == NULL || subv[0] == '\0') { /* lazy way for 3dTcat special */
305 svar = TCAT_get_subv( DSET_NVALS(dset) , subv ) ; /* ZSS March 2010 */
306 } else {
307 svar = MCW_get_thd_intlist (dset, subv); /* ZSS March 2010 */
308 }
309 if( svar == NULL || svar[0] <= 0 ){
310 ERROR_exit("can't decipher index codes from %s%s\n",dname,subv) ;
311 }
312 ADDTO_XTARR(TCAT_subv,svar) ; /* list of sub-brick selectors */
313
314 max_nsub = MAX( max_nsub , svar[0] ) ;
315
316 if( TCAT_rlt == 3 && svar[0] < 3 ) /* 16 Sep 1999 */
317 WARNING_message(
318 "-rlt++ option won't work properly with"
319 " less than 3 sub-bricks per input dataset!") ;
320 } /* end of loop over filename expansion*/
321
322 if( !ebad ) MCW_free_expand( nexp , fexp ) ;
323
324 } /* end of filename expansion stuff */
325
326 } /* end of while loop over command line arguments */
327
328 /*--- final sanity checks ---*/
329
330 if( max_nsub < 3 && TCAT_rlt ){
331 WARNING_message("can't apply -rlt option -- "
332 "Not enough points per input dataset." ) ;
333 TCAT_rlt = 0 ;
334 }
335
336 if( TCAT_rlt && TCAT_dry ){
337 WARNING_message("-rlt option does nothing with -dry!") ;
338 TCAT_rlt = 0 ;
339 }
340
341 return ;
342 }
343
344 /*-----------------------------------------------------------------------
345 Decode a string like [1..3,5..9(2)] into an array of integers.
346 -------------------------------------------------------------------------*/
347
TCAT_get_subv(int nvals,char * str)348 int * TCAT_get_subv( int nvals , char *str )
349 {
350 int *subv = NULL ;
351 int ii , ipos , nout , slen ;
352 int ibot,itop,istep , nused ;
353 char * cpt ;
354
355 /* Meaningless input? */
356
357 if( nvals < 1 ) return NULL ;
358
359 /* No selection list ==> select it all */
360
361 if( str == NULL || str[0] == '\0' ){
362 subv = (int *) RwcMalloc( sizeof(int) * (nvals+1) ) ;
363 subv[0] = nvals ;
364 for( ii=0 ; ii < nvals ; ii++ ) subv[ii+1] = ii ;
365 return subv ;
366 }
367
368 /* skip initial '[' */
369
370 subv = (int *) RwcMalloc( sizeof(int) * 2 ) ;
371 subv[0] = nout = 0 ;
372
373 ipos = 0 ;
374 if( str[ipos] == '[' ) ipos++ ;
375
376 /* do we have a count string in there ? */
377 if (strstr(str,"count ")) {
378 return(get_count_intlist (str, &ii, nvals-1));
379 }
380
381
382 /*** loop through each sub-selector until end of input ***/
383
384 slen = strlen(str) ;
385 while( ipos < slen && str[ipos] != ']' ){
386
387 /** get starting value **/
388
389 if( str[ipos] == '$' ){ /* special case */
390 ibot = nvals-1 ; ipos++ ;
391 } else { /* decode an integer */
392 ibot = strtol( str+ipos , &cpt , 10 ) ;
393 if( ibot < 0 ){ myRwcFree(subv) ; return NULL ; }
394 if( ibot >= nvals ) ibot = nvals-1 ;
395 nused = (cpt-(str+ipos)) ;
396 if( ibot == 0 && nused == 0 ){ myRwcFree(subv) ; return NULL ; }
397 ipos += nused ;
398 }
399
400 /** if that's it for this sub-selector, add one value to list **/
401
402 if( str[ipos] == ',' || str[ipos] == '\0' || str[ipos] == ']' ){
403 nout++ ;
404 subv = (int *) RwcRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
405 subv[0] = nout ;
406 subv[nout] = ibot ;
407 ipos++ ; continue ; /* re-start loop at next sub-selector */
408 }
409
410 /** otherwise, must have '..' or '-' as next inputs **/
411
412 if( str[ipos] == '-' ){
413 ipos++ ;
414 } else if( str[ipos] == '.' && str[ipos+1] == '.' ){
415 ipos++ ; ipos++ ;
416 } else {
417 myRwcFree(subv) ; return NULL ;
418 }
419
420 /** get ending value for loop now **/
421
422 if( str[ipos] == '$' ){ /* special case */
423 itop = nvals-1 ; ipos++ ;
424 } else { /* decode an integer */
425 itop = strtol( str+ipos , &cpt , 10 ) ;
426 if( itop < 0 ){ myRwcFree(subv) ; return NULL ; }
427 if( itop >= nvals ) itop = nvals-1 ;
428 nused = (cpt-(str+ipos)) ;
429 if( itop == 0 && nused == 0 ){ myRwcFree(subv) ; return NULL ; }
430 ipos += nused ;
431 }
432
433 /** set default loop step **/
434
435 istep = (ibot <= itop) ? 1 : -1 ;
436
437 /** check if we have a non-default loop step **/
438
439 if( str[ipos] == '(' ){ /* decode an integer */
440 ipos++ ;
441 istep = strtol( str+ipos , &cpt , 10 ) ;
442 if( istep == 0 ){ myRwcFree(subv) ; return NULL ; }
443 nused = (cpt-(str+ipos)) ;
444 ipos += nused ;
445 if( str[ipos] == ')' ) ipos++ ;
446 }
447
448 /** add values to output **/
449
450 for( ii=ibot ; (ii-itop)*istep <= 0 ; ii += istep ){
451 nout++ ;
452 subv = (int *) RwcRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
453 subv[0] = nout ;
454 subv[nout] = ii ;
455 }
456
457 /** check if we have a comma to skip over **/
458
459 if( str[ipos] == ',' ) ipos++ ;
460
461 } /* end of loop through selector string */
462
463 return subv ;
464 }
465
466 /*-------------------------------------------------------------------------*/
467
TCAT_Syntax(void)468 void TCAT_Syntax(void)
469 {
470 printf(
471 "Concatenate sub-bricks from input datasets into one big 3D+time dataset.\n"
472 "Usage: 3dTcat options\n"
473 "where the options are:\n"
474 ) ;
475
476 printf(
477 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
478 " OR -output pname [default='tcat']\n"
479 "\n"
480 " -session dir = Use 'dir' for the output dataset session directory.\n"
481 " [default='./'=current working directory]\n"
482 " -glueto fname = Append bricks to the end of the 'fname' dataset.\n"
483 " This command is an alternative to the -prefix \n"
484 " and -session commands. \n"
485 " -dry = Execute a 'dry run'; that is, only print out\n"
486 " what would be done. This is useful when\n"
487 " combining sub-bricks from multiple inputs.\n"
488 " -verb = Print out some verbose output as the program\n"
489 " proceeds (-dry implies -verb).\n"
490 " Using -verb twice results in quite lengthy output.\n"
491 " -rlt = Remove linear trends in each voxel time series loaded\n"
492 " from each input dataset, SEPARATELY. That is, the\n"
493 " data from each dataset is detrended separately.\n"
494 " At least 3 sub-bricks from a dataset must be input\n"
495 " for this option to apply.\n"
496 " Notes: (1) -rlt removes the least squares fit of 'a+b*t'\n"
497 " to each voxel time series; this means that\n"
498 " the mean is removed as well as the trend.\n"
499 " This effect makes it impractical to compute\n"
500 " the %% Change using AFNI's internal FIM.\n"
501 " (2) To have the mean of each dataset time series added\n"
502 " back in, use this option in the form '-rlt+'.\n"
503 " In this case, only the slope 'b*t' is removed.\n"
504 " (3) To have the overall mean of all dataset time\n"
505 " series added back in, use this option in the\n"
506 " form '-rlt++'. In this case, 'a+b*t' is removed\n"
507 " from each input dataset separately, and the\n"
508 " mean of all input datasets is added back in at\n"
509 " the end. (This option will work properly only\n"
510 " if all input datasets use at least 3 sub-bricks!)\n"
511 " (4) -rlt can be used on datasets that contain shorts\n"
512 " or floats, but not on complex- or byte-valued\n"
513 " datasets.\n"
514 " -relabel = Replace any sub-brick labels in an input dataset\n"
515 " with the input dataset name -- this might help\n"
516 " identify the sub-bricks in the output dataset.\n"
517 "\n"
518 " -tpattern PATTERN = Specify the timing pattern for the output\n"
519 " dataset, using patterns described in the\n"
520 " 'to3d -help' output (alt+z, seq, alt-z2, etc).\n"
521 "\n"
522 " -tr TR = Specify the TR (in seconds) for the output dataset.\n"
523 "\n"
524 " -DAFNI_GLOB_SELECTORS=YES\n"
525 " Setting the environment variable AFNI_GLOB_SELECTORS\n"
526 " to YES (as done temporarily with this option) means\n"
527 " that sub-brick selectors '[..]' will not be used\n"
528 " as wildcards. For example:\n"
529 " 3dTcat -DAFNI_GLOB_SELECTORS=YES -relabel -prefix EPIzero 'rest_*+tlrc.HEAD[0]'\n"
530 " will work to make a dataset with the #0 sub-brick\n"
531 " from each of a number of 3D+time datasets.\n"
532 " ** Note that the entire dataset specification is in quotes\n"
533 " to prevent the shell from doing the '*' wildcard expansion\n"
534 " -- it will be done inside the program itself, after the\n"
535 " sub-brick selector is temporarily detached from the string\n"
536 " -- and then a copy of the selector is re-attached to each\n"
537 " expanded filename.\n"
538 " ** Very few other AFNI '3d' programs do internal\n"
539 " wildcard expansion -- most of them rely on the shell.\n"
540 "\n"
541 "Command line arguments after the above are taken as input datasets.\n"
542 "A dataset is specified using one of these forms:\n"
543 " prefix+view\n"
544 " prefix+view.HEAD\n"
545 " prefix+view.BRIK\n"
546 " prefix.nii\n"
547 " prefix.nii.gz\n"
548 "\n"
549 "SUB-BRICK SELECTION:\n"
550 "--------------------\n"
551 "You can also add a sub-brick selection list after the end of the\n"
552 "dataset name. This allows only a subset of the sub-bricks to be\n"
553 "included into the output (by default, all of the input dataset\n"
554 "is copied into the output). A sub-brick selection list looks like\n"
555 "one of the following forms:\n"
556 " fred+orig[5] ==> use only sub-brick #5\n"
557 " fred+orig[5,9,17] ==> use #5, #9, and #17\n"
558 " fred+orig[5..8] or [5-8] ==> use #5, #6, #7, and #8\n"
559 " fred+orig[5..13(2)] or [5-13(2)] ==> use #5, #7, #9, #11, and #13\n"
560 "Sub-brick indexes start at 0. You can use the character '$'\n"
561 "to indicate the last sub-brick in a dataset; for example, you\n"
562 "can select every third sub-brick by using the selection list\n"
563 " fred+orig[0..$(3)]\n"
564 "You can reverse the order of sub-bricks with a list like\n"
565 " fred+origh[$..0(-1)]\n"
566 "(Exactly WHY you might want to time-reverse a dataset is a mystery.)\n"
567 "\n"
568 "You can also use a syntax based on the usage of the program count.\n"
569 "This would be most useful when randomizing (shuffling) the order of\n"
570 "the sub-bricks. Example:\n"
571 " fred+orig[count -seed 2 5 11 s] is equivalent to something like:\n"
572 " fred+orig[ 6, 5, 11, 10, 9, 8, 7] \n"
573 "You could also do: fred+orig[`count -seed 2 -digits 1 -suffix ',' 5 11 s`]\n"
574 "but if you have lots of numbers, the command line would get too\n"
575 "long for the shell to process it properly. Omit the seed option if\n"
576 "you want the code to generate a seed automatically.\n"
577 "You cannot mix and match count syntax with other selection gimmicks.\n"
578 "\n"
579 "If you have a lot of bricks to select in a particular order, you will\n"
580 "also run into name length problems. One solution is to put the indices\n"
581 "in a .1D file then use the following syntax. For example, say you have\n"
582 "the selection in file reorder.1D. You can extract the sub-bricks with:\n"
583 " fred+orig'[1dcat reorder.1D]' \n"
584 "As with count, you cannot mix and match 1dcat syntax with other \n"
585 "selection gimmicks.\n"
586 "\n"
587 "NOTES:\n"
588 "------\n"
589 "You can also add a sub-brick selection list after the end of the\n"
590 "* The TR and other time-axis properties are taken from the\n"
591 " first input dataset that is itself 3D+time. If no input\n"
592 " datasets contain such information, then TR is set to 1.0.\n"
593 " This can be altered later using the 3drefit program.\n"
594 "\n"
595 "* The sub-bricks are output in the order specified, which may\n"
596 " not be the order in the original datasets. For example, using\n"
597 " fred+orig[0..$(2),1..$(2)]\n"
598 " will cause the sub-bricks in fred+orig to be output into the\n"
599 " new dataset in an interleaved fashion. Using\n"
600 " fred+orig[$..0]\n"
601 " will reverse the order of the sub-bricks in the output.\n"
602 " If the -rlt option is used, the sub-bricks selected from each\n"
603 " input dataset will be re-ordered into the output dataset, and\n"
604 " then this sequence will be detrended.\n"
605 "\n"
606 "* You can use the '3dinfo' program to see how many sub-bricks\n"
607 " a 3D+time or a bucket dataset contains.\n"
608 "\n"
609 "* The '$', '(', ')', '[', and ']' characters are special to\n"
610 " the shell, so you will have to escape them. This is most easily\n"
611 " done by putting the entire dataset plus selection list inside\n"
612 " single quotes, as in 'fred+orig[5..7,9]'.\n"
613 "\n"
614 "* You may wish/need to use the 3drefit program on the output\n"
615 " dataset to modify some of the .HEAD file parameters.\n"
616 "\n"
617 "* The program does internal wildcard expansion on the filenames\n"
618 " provided to define the datasets. The software first strips the\n"
619 " sub-brick selector string '[...]' off the end of each filename\n"
620 " BEFORE wildcard expansion, then re-appends it to the results\n"
621 " AFTER the expansion; for example, '*+orig.HEAD[4..7]' might\n"
622 " expand to 'fred+orig.HEAD[4..7]' and 'wilma+orig.HEAD[4..7]'.\n"
623 " ++ However, the '[...]' construct is also a shell wildcard,\n"
624 " It is not practicable to use this feature for filename\n"
625 " selection with 3dTcat if you are also using sub-brick\n"
626 " selectors.\n"
627 " ++ Since wildcard expansion looks for whole filenames, you must\n"
628 " use wildcard expansion in the form (e.g.) of '*+orig.HEAD',\n"
629 " NOT '*+orig' -- since the latter form doesn't match filenames.\n"
630 " ++ Don't use '*+orig.*' since that will match both the .BRIK and\n"
631 " .HEAD files, and each dataset will end up being read in twice!\n"
632 " ++ If you want to see the filename expansion results, run 3dTcat\n"
633 " with the option '-DAFNI_GLOB_DEBUG=YES'\n"
634 ) ;
635
636 PRINT_COMPILE_DATE ; exit(0) ;
637 }
638
639 /*-------------------------------------------------------------------------*/
640
main(int argc,char * argv[])641 int main( int argc , char *argv[] )
642 {
643 int ninp , ids , nv , iv,jv,kv , ivout , new_nvals , ivbot,ivtop ;
644 THD_3dim_dataset *new_dset=NULL , * dset=NULL ;
645 char buf[256] ;
646 float *rlt0=NULL , *rlt1=NULL ;
647 float *rltsum=NULL ; /* 16 Sep 1999 */
648 int nrltsum ;
649 float dTR , nTR;
650 double angle=0.0 ;
651
652 /*** read input options ***/
653
654 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) TCAT_Syntax() ;
655
656 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
657
658 mainENTRY("3dTcat main"); machdep() ; PRINT_VERSION("3dTcat") ;
659 set_obliquity_report(0); /* silence obliquity */
660 (void)COX_clock_time() ;
661
662 { int new_argc ; char ** new_argv ;
663 addto_args( argc , argv , &new_argc , &new_argv ) ;
664 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
665 }
666
667 AFNI_logger("3dTcat",argc,argv) ;
668
669 TCAT_read_opts( argc , argv ) ;
670
671 /*** create new dataset (empty) ***/
672
673 ninp = TCAT_dsar->num ;
674 if( ninp < 1 ) ERROR_exit("No input datasets?") ;
675
676 new_nvals = 0 ;
677 for( ids=0 ; ids < ninp ; ids++ ) new_nvals += NSUBV(ids) ;
678
679 /* 14 Oct 2009: allow creation of single volume dataset [rickr] */
680 if( new_nvals < 1 )
681 ERROR_exit("Can't create 3D+time dataset with only %d sub-bricks!",
682 new_nvals) ;
683
684 if( TCAT_verb ) printf("-verb: output will have %d sub-bricks\n",new_nvals) ;
685
686 /** find 1st dataset that is time dependent **/
687
688 for( ids=0 ; ids < ninp ; ids++ ){
689 dset = DSUB(ids) ;
690 if( DSET_TIMESTEP(dset) > 0.0 ) break ;
691 }
692 if( ids == ninp ){ ids = 0 ; dset = DSUB(0) ; }
693
694 new_dset = EDIT_empty_copy( dset ) ; /* make a copy of its header */
695
696 /* 23 May 2005: check for axis consistency */
697
698 for( iv=0 ; iv < ninp ; iv++ ){
699 if( iv != ids ) {
700 if (!EQUIV_DATAXES(new_dset->daxes,DSUB(iv)->daxes) ) {
701 WARNING_message("%s grid mismatch with %s",
702 DSET_BRIKNAME(dset) , DSET_BRIKNAME(DSUB(iv)) ) ;
703 }
704 /* allow for small angle differences 22 May 2015 [rickr] */
705 angle = dset_obliquity_angle_diff(new_dset, DSUB(iv),OBLIQ_ANGLE_THRESH);
706 if (angle > 0.0) {
707 WARNING_message(
708 "dataset %s has an obliquity difference of %f degress with %s\n",
709 new_dset ,
710 angle, DSUB(iv) );
711 }
712 }
713 }
714
715 tross_Make_History( "3dTcat" , argc,argv , new_dset ) ;
716
717 /* modify its header */
718
719 EDIT_dset_items( new_dset ,
720 ADN_prefix , TCAT_output_prefix ,
721 ADN_directory_name, TCAT_session ,
722 ADN_type , TCAT_type ,
723 ADN_func_type , ISANATTYPE(TCAT_type) ? ANAT_EPI_TYPE
724 : FUNC_BUCK_TYPE ,
725 ADN_ntt , new_nvals , /* both ntt and nvals */
726 ADN_nvals , new_nvals , /* must be altered */
727 ADN_none ) ;
728
729 /* check if we have a valid time axis; if not, make one up */
730 /* (do this as an initialization, even if -TCAT_tr/pattern) */
731
732 if ( DSET_TIMESTEP(new_dset) <= 0.0 ){
733 float TR = 1.0 ;
734 float torg = 0.0 , tdur = 0.0 ;
735 int tunits = UNITS_SEC_TYPE ;
736
737 #if 0
738 for( ids=0 ; ids < ninp ; ids++ ){
739 dset = DSUB(ids) ;
740 if( DSET_TIMESTEP(dset) > 0.0 ){
741 TR = DSET_TIMESTEP(dset) ; tunits = DSET_TIMEUNITS(dset) ;
742 torg = DSET_TIMEORIGIN(dset) ; tdur = DSET_TIMEDURATION(dset) ;
743 break ;
744 }
745 }
746 #endif
747
748 EDIT_dset_items( new_dset ,
749 ADN_tunits , tunits ,
750 ADN_ttdel , TR ,
751 ADN_ttorg , torg ,
752 ADN_ttdur , tdur ,
753 ADN_none ) ;
754 if (DSET_NVALS(new_dset) > 1 && TCAT_tr <= 0.0) {
755 WARNING_message("Set TR of output dataset to 1.0 s") ;
756 }
757 }
758
759 /* 8 Mar 2013 [rickr] : */
760 /* now that time axis is set, adjust it for TR or tpat options */
761 if( TCAT_tr > 0.0 ) {
762 EDIT_dset_items( new_dset ,
763 ADN_tunits , UNITS_SEC_TYPE ,
764 ADN_ttdel , TCAT_tr ,
765 ADN_none ) ;
766 if( TCAT_verb ) printf("-verb: setting TR to %g\n", TCAT_tr);
767 }
768
769 if( TCAT_tpattern ) {
770 int nz = DSET_NZ(new_dset);
771 float tr = DSET_TIMESTEP(new_dset);
772 float * tpat = TS_parse_tpattern(nz, tr, TCAT_tpattern);
773 if( !tpat ) ERROR_exit("cannot get timing for nz=%d, tr=%g, tpat=%s",
774 nz, tr, TCAT_tpattern);
775
776 EDIT_dset_items( new_dset ,
777 ADN_nsl , nz ,
778 ADN_toff_sl , tpat ,
779 ADN_none ) ;
780 if( TCAT_verb ) {
781 printf("-verb: setting timing to :");
782 for( kv=0; kv < nz; kv++ ) printf(" %g", tpat[kv]);
783 putchar('\n');
784 }
785
786 free(tpat);
787 }
788
789 /* 10 Dec 2007: check if time steps are coherent */
790
791 nTR = DSET_TIMESTEP(new_dset) ;
792 for( ids=0 ; ids < ninp ; ids++ ){
793 dset = DSUB(ids) ; dTR = DSET_TIMESTEP(dset) ;
794 if( dTR > 0.0f && fabsf(dTR-nTR) > 0.001f &&
795 DSET_NVALS(new_dset) > 1)
796 WARNING_message("TR=%g in dataset %s; differs from output TR=%g",
797 dTR , DSET_HEADNAME(dset) , nTR ) ;
798 }
799
800 /* can't re-write existing dataset, unless glueing is used */
801
802 if (! TCAT_glue){
803 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(new_dset)) ){
804 ERROR_exit("file %s already exists!", DSET_HEADNAME(new_dset) ) ;
805 }
806 } else { /* if glueing is used, make the 'new'
807 dataset have the same idcode as the old one */
808
809 new_dset->idcode = DSUB(0) -> idcode ; /* copy the struct */
810 }
811
812 THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
813
814 /*** if needed, create space for detrending ***/
815
816 if( TCAT_rlt ){
817 rlt0 = (float *) malloc( sizeof(float) * TCAT_nvox ) ;
818 rlt1 = (float *) malloc( sizeof(float) * TCAT_nvox ) ;
819 if( rlt0 == NULL || rlt1 == NULL )
820 ERROR_exit("can't malloc memory for detrending!") ;
821
822 if( TCAT_rlt == 3 ){
823 rltsum = (float *) malloc( sizeof(float) * TCAT_nvox ) ;
824 if( rltsum == NULL )
825 ERROR_exit("can't malloc memory for detrending!") ;
826
827 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rltsum[iv] = 0.0 ;
828 nrltsum = 0 ;
829 }
830 }
831
832 /*** loop over input datasets ***/
833
834 if( ninp > 1 ) myRwcFree( new_dset->keywords ) ;
835
836 ivout = 0 ;
837 for( ids=0 ; ids < ninp ; ids++ ){
838 dset = DSUB(ids) ;
839 nv = NSUBV(ids) ;
840
841 if( ! TCAT_dry ){
842 if( TCAT_verb ) printf("-verb: loading %s\n",DSET_FILECODE(dset)) ;
843 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
844 }
845
846 /** loop over sub-bricks to output **/
847
848 ivbot = ivout ; /* save this for later */
849 for( iv=0 ; iv < nv ; iv++ ){
850 jv = SUBV(ids,iv) ; /* which sub-brick to use */
851
852 if( ! TCAT_dry ){
853 EDIT_substitute_brick( new_dset , ivout ,
854 DSET_BRICK_TYPE(dset,jv) , DSET_ARRAY(dset,jv) ) ;
855
856 /*----- If this sub-brick is from a bucket dataset,
857 preserve the label for this sub-brick -----*/
858
859 if( !TCAT_relabel && DSET_HAS_LABEL(dset,jv) )
860 sprintf (buf, "%s", DSET_BRICK_LABEL(dset,jv));
861 else
862 sprintf(buf,"%.16s[%d]",DSET_PREFIX(dset),jv) ;
863 EDIT_dset_items( new_dset, ADN_brick_label_one+ivout, buf, ADN_none );
864
865 EDIT_dset_items( new_dset ,
866 ADN_brick_fac_one+ivout, DSET_BRICK_FACTOR(dset,jv),
867 ADN_none ) ;
868
869 /** possibly write statistical parameters for this sub-brick **/
870
871 kv = DSET_BRICK_STATCODE(dset,jv) ;
872
873 if( FUNC_IS_STAT(kv) ){ /* input sub-brick has stat params */
874 int npar = FUNC_need_stat_aux[kv] , lv ;
875 float *par = (float *) malloc( sizeof(float) * (npar+2) ) ;
876 float *sax = DSET_BRICK_STATAUX(dset,jv) ;
877 par[0] = kv ;
878 par[1] = npar ;
879 for( lv=0 ; lv < npar ; lv++ )
880 par[lv+2] = (sax != NULL) ? sax[lv] : 0.0 ;
881
882 EDIT_dset_items(new_dset ,
883 ADN_brick_stataux_one+ivout , par ,
884 ADN_none ) ;
885 free(par) ;
886 }
887
888 /** print a message? **/
889
890 if( TCAT_verb > 1 ) printf("-verb: copied %s[%d] into %s[%d]\n" ,
891 DSET_FILECODE(dset) , jv ,
892 DSET_FILECODE(new_dset) , ivout ) ;
893 } else {
894 printf("-dry: would copy %s[%d] into %s[%d]\n" ,
895 DSET_FILECODE(dset) , jv ,
896 DSET_FILECODE(new_dset) , ivout ) ;
897 }
898
899 ivout++ ;
900 }
901 ivtop = ivout ; /* new_dset[ivbot..ivtop-1] are from the current dataset */
902
903 /** loop over all bricks in input dataset and
904 unload them if they aren't going into the output
905 (not required, but is done to economize on memory) **/
906
907 if( ! TCAT_dry && nv < DSET_NVALS(dset) ){
908
909 for( kv=0 ; kv < DSET_NVALS(dset) ; kv++ ){ /* all input sub-bricks */
910 for( iv=0 ; iv < nv ; iv++ ){ /* all output sub-bricks */
911 jv = SUBV(ids,iv) ;
912 if( jv == kv ) break ; /* input matches output */
913 }
914 if( iv == nv ) mri_free( DSET_BRICK(dset,kv) ) ;
915 }
916 }
917
918 /*** remove linear trend? ***/
919
920 if( TCAT_rlt ){
921
922 /* have enough data? */
923
924 if( ivtop-ivbot < 3 ){
925 if( TCAT_verb )
926 printf("-verb: skipping -rlt for %s\n",DSET_FILECODE(dset)) ;
927
928 } else {
929 float c0,c1,c2 , det , a0,a1,a2 , qq ;
930 int iv , ns , kk , err=0 ;
931
932 if( TCAT_verb )
933 printf("-verb: applying -rlt to data from %s\n",DSET_FILECODE(dset)) ;
934
935 /* compute weighting coefficients */
936
937 ns = ivtop - ivbot ; /* number of sub-bricks */
938 c0 = ns ; /* sum[ 1 ] */
939 c1 = 0.5 * ns * (ns-1) ; /* sum[ qq ] */
940 c2 = 0.16666667 * ns * (ns-1) * (2*ns-1) ; /* sum[ qq*qq ] */
941 det = c0*c2 - c1*c1 ;
942 a0 = c2 / det ; /* -1 */
943 a1 = -c1 / det ; /* [ c0 c1 ] */
944 a2 = c0 / det ; /* [ c1 c2 ] */
945
946 /* set voxel sums to 0 */
947
948 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rlt0[iv] = rlt1[iv] = 0.0 ;
949
950 /* compute voxel sums */
951
952 for( kk=ivbot ; kk < ivtop ; kk++ ){
953 qq = kk - ivbot ;
954 switch( DSET_BRICK_TYPE(new_dset,kk) ){
955 default:
956 err = 1 ;
957 WARNING_message(
958 "Warning: -rlt can't use datum type %s from %s",
959 MRI_TYPE_name[DSET_BRICK_TYPE(new_dset,kk)] ,
960 DSET_FILECODE(dset) ) ;
961 break ;
962
963 case MRI_short:{
964 short * bar = (short *) DSET_ARRAY(new_dset,kk) ;
965 float fac = DSET_BRICK_FACTOR(new_dset,kk) ;
966
967 if( fac == 0.0 ) fac = 1.0 ;
968 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
969 rlt0[iv] += (fac * bar[iv]) ; /* sum of voxel */
970 rlt1[iv] += (fac * bar[iv]) * qq ; /* sum of voxel*qq */
971 }
972 }
973 break ;
974
975 case MRI_float:{
976 float * bar = (float *) DSET_ARRAY(new_dset,kk) ;
977 float fac = DSET_BRICK_FACTOR(new_dset,kk) ;
978
979 if( fac == 0.0 ) fac = 1.0 ;
980 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
981 rlt0[iv] += (fac * bar[iv]) ;
982 rlt1[iv] += (fac * bar[iv]) * qq ;
983 }
984 }
985 break ;
986 }
987 if( err ) break ;
988 } /* end of loop over sub-bricks */
989
990 /* only do the detrending if no errors happened */
991
992 if( !err ){
993 float qmid = 0.0 ; /* 16 Sep 1999 */
994
995 for( iv=0 ; iv < TCAT_nvox ; iv++ ){ /* transform voxel sums */
996 c0 = a0 * rlt0[iv] + a1 * rlt1[iv] ;
997 c1 = a1 * rlt0[iv] + a2 * rlt1[iv] ;
998 rlt0[iv] = c0 ; rlt1[iv] = c1 ;
999 }
1000
1001 if( TCAT_rlt == 2 ){ /* 16 Sep 1999 */
1002 qmid = 0.5 * (ns-1) ;
1003 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rlt0[iv] = 0.0 ;
1004 } else if( TCAT_rlt == 3 ){
1005 nrltsum += ns ;
1006 for( iv=0 ; iv < TCAT_nvox ; iv++ )
1007 rltsum[iv] += (rlt0[iv] + (0.5*ns)*rlt1[iv])*ns ;
1008 }
1009
1010 for( kk=ivbot ; kk < ivtop ; kk++ ){ /* detrend */
1011 qq = kk - ivbot ;
1012 switch( DSET_BRICK_TYPE(new_dset,kk) ){
1013 default: break ; /* should not happen, I hope */
1014
1015 #undef ROUND
1016 #define ROUND(qq) ((short)rint((qq)+0.00001))
1017
1018 case MRI_short:{
1019 short * bar = (short *) DSET_ARRAY(new_dset,kk) ;
1020 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
1021
1022 if( fac == 0.0 ) fac = 1.0 ;
1023 finv = 1.0 / fac ;
1024 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
1025 val = fac*bar[iv] - rlt0[iv] - rlt1[iv]*(qq-qmid) ;
1026 bar[iv] = ROUND(finv*val) ;
1027 }
1028 }
1029 break ;
1030
1031 case MRI_float:{
1032 float * bar = (float *) DSET_ARRAY(new_dset,kk) ;
1033 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
1034
1035 if( fac == 0.0 ) fac = 1.0 ;
1036 finv = 1.0 / fac ;
1037 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
1038 val = fac*bar[iv] - rlt0[iv] - rlt1[iv]*(qq-qmid) ;
1039 bar[iv] = (finv*val) ;
1040 }
1041 }
1042 break ;
1043 }
1044 }
1045 }
1046 }
1047 } /* end of -rlt */
1048
1049 } /* end of loop over input datasets */
1050
1051 /* 16 Sep 1999: add overall average back in */
1052
1053 if( TCAT_rlt == 3 && rltsum != NULL && nrltsum > 0 ){
1054 float scl = 1.0/nrltsum ; int kk ;
1055
1056 for( iv=0 ; iv < TCAT_nvox ; iv++ ) rltsum[iv] *= scl ;
1057
1058 for( kk=0 ; kk < new_nvals ; kk++ ){
1059 switch( DSET_BRICK_TYPE(new_dset,kk) ){
1060 default: break ; /* should not happen, I hope */
1061 case MRI_short:{
1062 short * bar = (short *) DSET_ARRAY(new_dset,kk) ;
1063 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
1064
1065 if( fac == 0.0 ) fac = 1.0 ;
1066 finv = 1.0 / fac ;
1067 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
1068 val = fac*bar[iv] + rltsum[iv] ; bar[iv] = ROUND(finv*val) ;
1069 }
1070 }
1071 break ;
1072
1073 case MRI_float:{
1074 float * bar = (float *) DSET_ARRAY(new_dset,kk) ;
1075 float fac = DSET_BRICK_FACTOR(new_dset,kk) , val,finv ;
1076
1077 if( fac == 0.0 ) fac = 1.0 ;
1078 finv = 1.0 / fac ;
1079 for( iv=0 ; iv < TCAT_nvox ; iv++ ){
1080 val = fac*bar[iv] + rltsum[iv] ; bar[iv] = (finv*val) ;
1081 }
1082 }
1083 break ;
1084 }
1085 }
1086 }
1087
1088 if( TCAT_rlt ){ free(rlt0); free(rlt1); if(rltsum!=NULL)free(rltsum); }
1089
1090 if( ! TCAT_dry ){
1091 if( TCAT_verb ) INFO_message("-verb: computing sub-brick statistics") ;
1092 THD_load_statistics( new_dset ) ;
1093 if( TCAT_glue ) putenv("AFNI_DECONFLICT=OVERWRITE") ;
1094 if( TCAT_glue && TCAT_ccode >= 0 )
1095 THD_set_write_compression(TCAT_ccode) ; /* 16 Mar 2010 */
1096 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
1097 if( TCAT_verb ) INFO_message("-verb: Wrote output to %s",
1098 DSET_BRIKNAME(new_dset) ) ;
1099 INFO_message("elapsed time = %.1f s",COX_clock_time()) ;
1100 }
1101
1102 exit(0) ;
1103 }
1104