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 creates AFNI "bucket" type datasets.
10
11 File: 3dbucket.c
12 Author: R. W. Cox
13 Date: 17 December 1997
14
15
16 Mod: Changes to implement "-glueto" command option.
17 Also, modified output to preserve sub-brick labels.
18 Author: B. D. Ward
19 Date: 04 February 1998
20
21 Mod: If more than one input dataset, copy command line history from the
22 first input dataset to the output bucket dataset.
23 Date: 14 March 2002
24
25 Mod: When verifying view type extension for -glueto dataset, scan
26 from the end (in case there are extra '+' characters).
27 Author: R. C. Reynolds
28 Date: 30 October 2003
29
30 */
31 /*---------------------------------------------------------------------------*/
32
33
34 #define PROGRAM_NAME "3dbucket" /* name of this program */
35 #define LAST_MOD_DATE "30 October 2003" /* date of last program mod */
36
37 #include "mrilib.h"
38
39 /*-------------------------- global data --------------------------*/
40
41 static THD_3dim_dataset_array * BUCK_dsar = NULL ;
42 static RwcPointer_array * BUCK_subv = NULL ;
43 static int BUCK_nvox = -1 ;
44 static int BUCK_dry = 0 ;
45 static int BUCK_verb = 0 ;
46 static int BUCK_type = -1 ;
47 static int BUCK_glue = 0 ;
48 static int BUCK_ccode = COMPRESS_NONE ; /* 16 Mar 2010 */
49
50 static char BUCK_output_prefix[THD_MAX_PREFIX] = "buck" ;
51 static char BUCK_session[THD_MAX_NAME] = "./" ;
52
53 #define NSUBV(id) ( ((int *)BUCK_subv->ar[(id)])[0] )
54 #define SUBV(id,jj) ( ((int *)BUCK_subv->ar[(id)])[(jj)+1] )
55 #define DSUB(id) DSET_IN_3DARR(BUCK_dsar,(id))
56
57 /*--------------------------- prototypes ---------------------------*/
58
59 void BUCK_read_opts( int , char ** ) ;
60 void BUCK_Syntax(int) ;
61 int * BUCK_get_subv( int , char * ) ;
62
63 /*--------------------------------------------------------------------
64 read the arguments, load the global variables
65 ----------------------------------------------------------------------*/
66
BUCK_read_opts(int argc,char * argv[])67 void BUCK_read_opts( int argc , char * argv[] )
68 {
69 int nopt = 1 , ii ;
70 char dname[THD_MAX_NAME] ;
71 char subv[THD_MAX_NAME] ;
72 char *cpt ;
73 THD_3dim_dataset *dset , *fset=NULL ;
74 int *svar ;
75 char *str;
76 int ok, ilen, nlen;
77
78 INIT_3DARR(BUCK_dsar) ;
79 INIT_XTARR(BUCK_subv) ;
80
81 while( nopt < argc ){
82 if( strcmp(argv[nopt],"-help") == 0 ||
83 strcmp(argv[nopt],"-h") == 0) {
84 BUCK_Syntax(strlen(argv[nopt])>3?2:1) ;
85 exit(0);
86 }
87 /**** -prefix prefix ****/
88
89 if( strncmp(argv[nopt],"-prefix",6) == 0 ||
90 strncmp(argv[nopt],"-output",6) == 0 ){
91 if (BUCK_glue){
92 fprintf(stderr,"-prefix and -glueto options are not compatible\n");
93 exit(1) ;
94 }
95 nopt++ ;
96 if( nopt >= argc ){
97 fprintf(stderr,"need argument after -prefix!\n") ; exit(1) ;
98 }
99 MCW_strncpy( BUCK_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
100 continue ;
101 }
102
103 /**** -session directory ****/
104
105 if( strncmp(argv[nopt],"-session",6) == 0 ){
106 if (BUCK_glue){
107 fprintf(stderr,
108 "-session and -glueto options are not compatible\n");
109 exit(1) ;
110 }
111 nopt++ ;
112 if( nopt >= argc ){
113 fprintf(stderr,"need argument after -session!\n") ; exit(1) ;
114 }
115 MCW_strncpy( BUCK_session , argv[nopt++] , THD_MAX_NAME ) ;
116 continue ;
117 }
118
119 if( strncmp(argv[nopt],"-dry",3) == 0 ){
120 BUCK_dry = BUCK_verb = 1 ;
121 nopt++ ; continue ;
122 }
123
124 if( strncmp(argv[nopt],"-fbuc",4) == 0 ){
125 BUCK_type = HEAD_FUNC_TYPE ;
126 nopt++ ; continue ;
127 }
128
129 if( strncmp(argv[nopt],"-abuc",4) == 0 ){
130 BUCK_type = HEAD_ANAT_TYPE ;
131 nopt++ ; continue ;
132 }
133
134 if( strncmp(argv[nopt],"-verb",5) == 0 ){
135 BUCK_verb = 1 ;
136 nopt++ ; continue ;
137 }
138
139 if( strncmp(argv[nopt],"-glueto",5) == 0 ||
140 strncmp(argv[nopt],"-aglueto",5) == 0){ /* ZSS April 16 2010 */
141 if( strncmp(BUCK_output_prefix, "buck", 5) != 0 ){
142 fprintf(stderr,
143 "-prefix, -glueto, and -aglueto options are not compatible\n");
144 exit(1) ;
145 }
146 if( strncmp(BUCK_session, "./", 5) != 0 ){
147 fprintf(stderr,
148 "-session, -glueto, and -aglueto options are not compatible\n");
149 exit(1) ;
150 }
151 nopt++ ;
152 if( nopt >= argc ){
153 fprintf(stderr,"need argument after -glueto or -aglueto!\n") ;
154 exit(1) ;
155 }
156 if (strncmp(argv[nopt-1],"-aglueto",5) == 0) { /* ZSS April 16 2010 */
157 THD_3dim_dataset *ddd = THD_open_dataset( argv[nopt] );
158 if( !ISVALID_DSET(ddd) ){
159 /* treat as -prefix */
160 MCW_strncpy( BUCK_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
161 continue ;
162 } else {
163 /* go on as standard -glueto option */
164 }
165 }
166 BUCK_glue = 1 ;
167
168 /*----- Verify that file name ends in View Type -----*/
169 ok = 1;
170 nlen = strlen(argv[nopt]);
171 if (nlen <= 5) ok = 0;
172
173 if (ok)
174 {
175 #if 0 /* old code - scan from end, instead */
176
177 for (ilen = 0; ilen < nlen; ilen++)
178 {
179 str = argv[nopt] + ilen;
180 if (str[0] == '+') break;
181 }
182 if (ilen == nlen) ok = 0;
183 #endif
184
185 /* scan from end for view type extension, require one char */
186 /* 30 Oct 2003 [rickr] */
187 for (ilen = nlen - 1; ilen > 0; ilen--)
188 {
189 str = argv[nopt] + ilen;
190 if (str[0] == '+') break;
191 }
192 if (ilen == 0) ok = 0;
193 }
194
195 if (ok)
196 {
197 str = argv[nopt] + ilen + 1;
198
199 for (ii=FIRST_VIEW_TYPE ; ii <= LAST_VIEW_TYPE ; ii++)
200 if (! strncmp(str,VIEW_codestr[ii],4)) break ;
201
202 if( ii > LAST_VIEW_TYPE ) ok = 0;
203 }
204
205 if (! ok)
206 {
207 fprintf(stderr,
208 "File name must end in +orig, +acpc, or +tlrc after -glueto\n"
209 "(consider: 3dbucket -prefix dsetA -overwrite dsetA dsetB ...)\n");
210 exit(1);
211 }
212
213 /*----- Remove View Type from string to make output prefix -----*/
214 MCW_strncpy( BUCK_output_prefix , argv[nopt] , ilen+1) ;
215
216 /*----- Note: no "continue" statement here. File name will now
217 be processed as an input dataset -----*/
218 }
219
220 if( strncmp(argv[nopt],"-aglueto",5) == 0 ){
221 if( strncmp(BUCK_output_prefix, "buck", 5) != 0 ){
222 fprintf(stderr,"-prefix and -aglueto options are not compatible.\n"
223 "Make sure you do not have two -agluto options on command.\n");
224 exit(1) ;
225 }
226 if( strncmp(BUCK_session, "./", 5) != 0 ){
227 fprintf(stderr,
228 "-session and -aglueto options are not compatible\n");
229 exit(1) ;
230 }
231 BUCK_glue = 1 ;
232 nopt++ ;
233 if( nopt >= argc ){
234 fprintf(stderr,"need argument after -glueto!\n") ; exit(1) ;
235 }
236
237 /*----- Verify that file name ends in View Type -----*/
238 ok = 1;
239 nlen = strlen(argv[nopt]);
240 if (nlen <= 5) ok = 0;
241
242 if (ok)
243 {
244 #if 0 /* old code - scan from end, instead */
245
246 for (ilen = 0; ilen < nlen; ilen++)
247 {
248 str = argv[nopt] + ilen;
249 if (str[0] == '+') break;
250 }
251 if (ilen == nlen) ok = 0;
252 #endif
253
254 /* scan from end for view type extension, require one char */
255 /* 30 Oct 2003 [rickr] */
256 for (ilen = nlen - 1; ilen > 0; ilen--)
257 {
258 str = argv[nopt] + ilen;
259 if (str[0] == '+') break;
260 }
261 if (ilen == 0) ok = 0;
262 }
263
264 if (ok)
265 {
266 str = argv[nopt] + ilen + 1;
267
268 for (ii=FIRST_VIEW_TYPE ; ii <= LAST_VIEW_TYPE ; ii++)
269 if (! strncmp(str,VIEW_codestr[ii],4)) break ;
270
271 if( ii > LAST_VIEW_TYPE ) ok = 0;
272 }
273
274 if (! ok)
275 {
276 fprintf(stderr,
277 "File name must end in +orig, +acpc, or +tlrc after -glueto\n"
278 "(consider: 3dbucket -prefix dsetA -overwrite dsetA dsetB ...)\n");
279 exit(1);
280 }
281
282 /*----- Remove View Type from string to make output prefix -----*/
283 MCW_strncpy( BUCK_output_prefix , argv[nopt] , ilen+1) ;
284
285 /*----- Note: no "continue" statement here. File name will now
286 be processed as an input dataset -----*/
287 }
288
289 if( argv[nopt][0] == '-' ){
290 fprintf(stderr,"Unknown option: %s\n",argv[nopt]) ;
291 suggest_best_prog_option(argv[0], argv[nopt]);
292 exit(1) ;
293 }
294
295 /**** read dataset ****/
296
297 cpt = strstr(argv[nopt],"[") ;
298 if( cpt == NULL ){
299 if (strlen(argv[nopt]) > THD_MAX_NAME-1) {
300 ERROR_exit( "Too long a filename for '%s'\n"
301 "Maximum limit is %d\n", argv[nopt], THD_MAX_NAME-1);
302 }
303 strcpy(dname,argv[nopt]) ;
304 subv[0] = '\0' ; /* make sure subv is reset ZSS Nov. 2010*/
305 } else if( cpt == argv[nopt] ){
306 fprintf(stderr,"illegal dataset specifier: %s\n",argv[nopt]) ;
307 exit(1) ;
308 } else {
309 ii = cpt - argv[nopt] ;
310 if (ii > THD_MAX_NAME-1) {
311 ERROR_exit( "Too long a filename for '%s'\n"
312 "Maximum character limit is %d, have %d\\n",
313 argv[nopt], THD_MAX_NAME-1, ii);
314 }
315 if (strlen(argv[nopt])-ii > THD_MAX_NAME-1) {
316 ERROR_exit( "Too long a sub-brick selection for '%s'\n"
317 "Maximum limit is %d, have %d\n"
318 "Consider using '[1dcat FF.1D]' or '[count ...]' methods\n"
319 "for sub-brick selection. See 3dTcat -help for details.\n",
320 argv[nopt], THD_MAX_NAME-1, strlen(argv[nopt])-ii);
321 }
322 memcpy(dname,argv[nopt],ii) ; dname[ii] = '\0' ;
323 strcpy(subv,cpt) ;
324 }
325 nopt++ ;
326
327 dset = THD_open_one_dataset( dname ) ;
328 /* akin to 3dTcat, handle failure in open_one
329 (todo: remove open_one from these programs) 18 Apr 2016 [rickr] */
330 if ( dset == NULL ) {
331 dset = THD_open_dataset( argv[nopt-1] ) ;
332 subv[0] = '\0'; /* use selectors via THD_O_D */
333 if( dset && BUCK_verb )
334 INFO_message("failure w/open_one_dset, success w/open_dset");
335 }
336 if( dset == NULL ){
337 fprintf(stderr,"can't open dataset %s\n",dname) ; exit(1) ;
338 }
339 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
340
341 if( BUCK_type < 0 ) BUCK_type = dset->type ;
342
343 BUCK_ccode = COMPRESS_filecode(dset->dblk->diskptr->brick_name) ;
344 /* 16 Mar 2010 */
345
346 ii = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz ;
347 if( BUCK_nvox < 0 ){
348 BUCK_nvox = ii ; fset = dset ;
349 } else if( ii != BUCK_nvox ){
350 ERROR_exit("Dataset %s differs in size from first one",dname);
351 } else if( !EQUIV_GRIDS(dset,fset) ){
352 WARNING_message("Dataset %s grid differs from first one",dname) ;
353 }
354 ADDTO_3DARR(BUCK_dsar,dset) ;
355 if (subv == NULL || subv[0] == '\0') { /* lazy way for 3dbucket special */
356 svar = BUCK_get_subv( DSET_NVALS(dset) , subv ) ; /* ZSS Dec 09 */
357 } else {
358 svar = MCW_get_thd_intlist (dset, subv); /* ZSS Dec 09 */
359 }
360 if( svar == NULL || svar[0] <= 0 ){
361 fprintf(stderr,"can't decipher index codes from %s%s\n",dname,subv) ;
362 exit(1) ;
363 }
364 ADDTO_XTARR(BUCK_subv,svar) ;
365
366 } /* end of loop over command line arguments */
367
368 if( argc < 2) {
369 ERROR_message("Too few options");
370 BUCK_Syntax(0) ;
371 exit(1);
372 }
373 return ;
374 }
375
376 /*------------------------------------------------------------------*/
377
BUCK_get_subv(int nvals,char * str)378 int * BUCK_get_subv( int nvals , char * str )
379 {
380 int * subv = NULL ;
381 int ii , ipos , nout , slen ;
382 int ibot,itop,istep , nused ;
383 char * cpt ;
384
385 /* Meaningless input? */
386
387 if( nvals < 1 ) return NULL ;
388
389 /* No selection list ==> select it all */
390
391 if( str == NULL || str[0] == '\0' ){
392 subv = (int *) RwcMalloc( sizeof(int) * (nvals+1) ) ;
393 subv[0] = nvals ;
394 for( ii=0 ; ii < nvals ; ii++ ) subv[ii+1] = ii ;
395 return subv ;
396 }
397
398 /* skip initial '[' */
399
400 subv = (int *) RwcMalloc( sizeof(int) * 2 ) ;
401 subv[0] = nout = 0 ;
402
403 ipos = 0 ;
404 if( str[ipos] == '[' ) ipos++ ;
405
406 /*** loop through each sub-selector until end of input ***/
407
408 slen = strlen(str) ;
409 while( ipos < slen && str[ipos] != ']' ){
410
411 /** get starting value **/
412
413 if( str[ipos] == '$' ){ /* special case */
414 ibot = nvals-1 ; ipos++ ;
415 } else { /* decode an integer */
416 ibot = strtol( str+ipos , &cpt , 10 ) ;
417 if( ibot < 0 ){ myRwcFree(subv) ; return NULL ; }
418 if( ibot >= nvals ) ibot = nvals-1 ;
419 nused = (cpt-(str+ipos)) ;
420 if( ibot == 0 && nused == 0 ){ myRwcFree(subv) ; return NULL ; }
421 ipos += nused ;
422 }
423
424 /** if that's it for this sub-selector, add one value to list **/
425
426 if( str[ipos] == ',' || str[ipos] == '\0' || str[ipos] == ']' ){
427 nout++ ;
428 subv = (int *) RwcRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
429 subv[0] = nout ;
430 subv[nout] = ibot ;
431 ipos++ ; continue ; /* re-start loop at next sub-selector */
432 }
433
434 /** otherwise, must have '..' or '-' as next inputs **/
435
436 if( str[ipos] == '-' ){
437 ipos++ ;
438 } else if( str[ipos] == '.' && str[ipos+1] == '.' ){
439 ipos++ ; ipos++ ;
440 } else {
441 myRwcFree(subv) ; return NULL ;
442 }
443
444 /** get ending value for loop now **/
445
446 if( str[ipos] == '$' ){ /* special case */
447 itop = nvals-1 ; ipos++ ;
448 } else { /* decode an integer */
449 itop = strtol( str+ipos , &cpt , 10 ) ;
450 if( itop < 0 ){ myRwcFree(subv) ; return NULL ; }
451 if( itop >= nvals ) itop = nvals-1 ;
452 nused = (cpt-(str+ipos)) ;
453 if( itop == 0 && nused == 0 ){ myRwcFree(subv) ; return NULL ; }
454 ipos += nused ;
455 }
456
457 /** set default loop step **/
458
459 istep = (ibot <= itop) ? 1 : -1 ;
460
461 /** check if we have a non-default loop step **/
462
463 if( str[ipos] == '(' ){ /* decode an integer */
464 ipos++ ;
465 istep = strtol( str+ipos , &cpt , 10 ) ;
466 if( istep == 0 ){ myRwcFree(subv) ; return NULL ; }
467 nused = (cpt-(str+ipos)) ;
468 ipos += nused ;
469 if( str[ipos] == ')' ) ipos++ ;
470 }
471
472 /** add values to output **/
473
474 for( ii=ibot ; (ii-itop)*istep <= 0 ; ii += istep ){
475 nout++ ;
476 subv = (int *) RwcRealloc( (char *)subv , sizeof(int) * (nout+1) ) ;
477 subv[0] = nout ;
478 subv[nout] = ii ;
479 }
480
481 /** check if we have a comma to skip over **/
482
483 if( str[ipos] == ',' ) ipos++ ;
484
485 } /* end of loop through selector string */
486
487 return subv ;
488 }
489
490 /*------------------------------------------------------------------*/
491
BUCK_Syntax(int detail)492 void BUCK_Syntax(int detail)
493 {
494 printf(
495 "Concatenate sub-bricks from input datasets into one big 'bucket' dataset. ~1~\n"
496 "Usage: 3dbucket options\n"
497 "where the options are: ~1~\n"
498 ) ;
499
500 printf(
501 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
502 " OR -output pname [default='buck']\n"
503 "\n"
504 " -session dir = Use 'dir' for the output dataset session directory.\n"
505 " [default='./'=current working directory]\n"
506 " -glueto fname = Append bricks to the end of the 'fname' dataset.\n"
507 " This command is an alternative to the -prefix \n"
508 " and -session commands.\n"
509 " * Note that fname should include the view, as in\n"
510 " 3dbucket -glueto newset+orig oldset+orig'[7]'\n"
511 " -aglueto fname= If fname dset does not exist, create it (like -prefix).\n"
512 " Otherwise append to fname (like -glueto).\n"
513 " This option is useful when appending in a loop.\n"
514 " * As with -glueto, fname should include the view, e.g.\n"
515 " 3dbucket -aglueto newset+orig oldset+orig'[7]'\n"
516 " -dry = Execute a 'dry run'; that is, only print out\n"
517 " what would be done. This is useful when\n"
518 " combining sub-bricks from multiple inputs.\n"
519 " -verb = Print out some verbose output as the program\n"
520 " proceeds (-dry implies -verb).\n"
521 " -fbuc = Create a functional bucket.\n"
522 " -abuc = Create an anatomical bucket. If neither of\n"
523 " these options is given, the output type is\n"
524 " determined from the first input type.\n"
525 "\n"
526 "Command line arguments after the above are taken as input datasets.\n"
527 "A dataset is specified using one of these forms:\n"
528 " 'prefix+view', 'prefix+view.HEAD', or 'prefix+view.BRIK'.\n"
529 "You can also add a sub-brick selection list after the end of the\n"
530 "dataset name. This allows only a subset of the sub-bricks to be\n"
531 "included into the output (by default, all of the input dataset\n"
532 "is copied into the output). A sub-brick selection list looks like\n"
533 "one of the following forms:\n"
534 " fred+orig[5] ==> use only sub-brick #5\n"
535 " fred+orig[5,9,17] ==> use #5, #9, and #17\n"
536 " fred+orig[5..8] or [5-8] ==> use #5, #6, #7, and #8\n"
537 " fred+orig[5..13(2)] or [5-13(2)] ==> use #5, #7, #9, #11, and #13\n"
538 "Sub-brick indexes start at 0. You can use the character '$'\n"
539 "to indicate the last sub-brick in a dataset; for example, you\n"
540 "can select every third sub-brick by using the selection list\n"
541 " fred+orig[0..$(3)]\n"
542 "\n"
543 "Notes: ~1~\n"
544 "N.B.: The sub-bricks are output in the order specified, which may\n"
545 " not be the order in the original datasets. For example, using\n"
546 " fred+orig[0..$(2),1..$(2)]\n"
547 " will cause the sub-bricks in fred+orig to be output into the\n"
548 " new dataset in an interleaved fashion. Using\n"
549 " fred+orig[$..0]\n"
550 " will reverse the order of the sub-bricks in the output.\n"
551 "\n"
552 "N.B.: Bucket datasets have multiple sub-bricks, but do NOT have\n"
553 " a time dimension. You can input sub-bricks from a 3D+time dataset\n"
554 " into a bucket dataset. You can use the '3dinfo' program to see\n"
555 " how many sub-bricks a 3D+time or a bucket dataset contains.\n"
556 "\n"
557 "N.B.: The '$', '(', ')', '[', and ']' characters are special to\n"
558 " the shell, so you will have to escape them. This is most easily\n"
559 " done by putting the entire dataset plus selection list inside\n"
560 " single quotes, as in 'fred+orig[5..7,9]'.\n"
561 "\n"
562 "N.B.: In non-bucket functional datasets (like the 'fico' datasets\n"
563 " output by FIM, or the 'fitt' datasets output by 3dttest), sub-brick\n"
564 " [0] is the 'intensity' and sub-brick [1] is the statistical parameter\n"
565 " used as a threshold. Thus, to create a bucket dataset using the\n"
566 " intensity from dataset A and the threshold from dataset B, and\n"
567 " calling the output dataset C, you would type\n"
568 " 3dbucket -prefix C -fbuc 'A+orig[0]' -fbuc 'B+orig[1]'\n"
569 "\n"
570 "WARNING: ~1~\n"
571 " Using this program, it is possible to create a dataset that\n"
572 " has different basic datum types for different sub-bricks\n"
573 " (e.g., shorts for brick 0, floats for brick 1).\n"
574 " Do NOT do this! Very few AFNI programs will work correctly\n"
575 " with such datasets!\n"
576 ) ;
577
578 PRINT_COMPILE_DATE ; return ;
579 }
580
581 /*------------------------------------------------------------------*/
582
main(int argc,char * argv[])583 int main( int argc , char * argv[] )
584 {
585 int ninp , ids , nv , iv,jv,kv , ivout , new_nvals , have_fdr = 0, nfdr = 0 ;
586 THD_3dim_dataset * new_dset=NULL , * dset ;
587 char buf[256] ;
588 double angle;
589
590 /*----- identify program -----*/
591 #if 0
592 printf ("\n\nProgram %s \n", PROGRAM_NAME);
593 printf ("Last revision: %s \n\n", LAST_MOD_DATE);
594 #endif
595
596 /*** read input options ***/
597
598
599 mainENTRY("3dbucket main"); machdep(); PRINT_VERSION("3dbucket") ;
600 set_obliquity_report(0); /* silence obliquity */
601
602 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
603
604 { int new_argc ; char ** new_argv ;
605 addto_args( argc , argv , &new_argc , &new_argv ) ;
606 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
607 }
608
609 AFNI_logger("3dbucket",argc,argv) ;
610
611 BUCK_read_opts( argc , argv ) ;
612
613 /*** create new dataset (empty) ***/
614 ninp = BUCK_dsar->num ;
615 if( ninp < 1 ){
616 fprintf(stderr,"*** No input datasets?\n") ; exit(1) ;
617 }
618
619 new_nvals = 0 ;
620 for( ids=0 ; ids < ninp ; ids++ ) new_nvals += NSUBV(ids) ;
621
622 if( BUCK_verb ) printf("-verb: output will have %d sub-bricks\n",new_nvals) ;
623
624 new_dset = EDIT_empty_copy( DSUB(0) ) ;
625
626 /* 23 May 2005: check for axis consistency */
627 /* 06 Feb 2008: and see if there are fdrcurves to perpetuate */
628
629 if( DSUB(0)->dblk->brick_fdrcurve ) have_fdr = 1 ;
630 for( iv=1 ; iv < ninp ; iv++ ){
631 if( !EQUIV_DATAXES(new_dset->daxes,DSUB(iv)->daxes) )
632 fprintf(stderr,"++ WARNING: %s grid mismatch with %s\n",
633 DSET_BRIKNAME(DSUB(0)) , DSET_BRIKNAME(DSUB(iv)) ) ;
634 if( DSUB(iv)->dblk->brick_fdrcurve ) have_fdr = 1 ;
635 /* allow for small angle differences 22 May 2015 [rickr] */
636 angle = dset_obliquity_angle_diff(new_dset, DSUB(iv), OBLIQ_ANGLE_THRESH);
637 if (angle > 0.0) {
638 WARNING_message(
639 "dataset %s has an obliquity difference of %f degress with %s\n",
640 new_dset ,
641 angle, DSUB(iv) );
642 }
643 }
644
645 /* if( ninp == 1 ) */ tross_Copy_History( DSUB(0) , new_dset ) ;
646 tross_Make_History( "3dbucket" , argc,argv , new_dset ) ;
647
648 EDIT_dset_items( new_dset ,
649 ADN_prefix , BUCK_output_prefix ,
650 ADN_directory_name, BUCK_session ,
651 ADN_type , BUCK_type ,
652 ADN_func_type , ISANATTYPE(BUCK_type) ? ANAT_BUCK_TYPE
653 : FUNC_BUCK_TYPE,
654 ADN_ntt , 0 ,
655 ADN_nvals , new_nvals ,
656 ADN_none ) ;
657
658 /* can't re-write existing dataset, unless glueing is used */
659
660 if (! BUCK_glue){
661 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(new_dset)) ){
662 fprintf(stderr,"*** Fatal error: file %s already exists!\n",
663 DSET_HEADNAME(new_dset) ) ;
664 exit(1) ;
665 }
666 } else { /* if glueing is used, make the 'new'
667 dataset have the same idcode as the old one */
668
669 new_dset->idcode = DSUB(0) -> idcode ; /* copy the struct */
670 }
671
672 THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
673
674 /* if there are fdr curves, allocate space 06 Feb 2008 [rickr] */
675 if( have_fdr ){
676 new_dset->dblk->brick_fdrcurve = (floatvec **)calloc(sizeof(floatvec *),
677 new_nvals) ;
678 if( !new_dset->dblk->brick_fdrcurve ){
679 fprintf(stderr,"** failed to alloc %d fdrcurves\n",new_nvals);
680 exit(1);
681 }
682 if( BUCK_verb ) printf("-verb: adding fdrcurve list\n");
683
684 new_dset->dblk->brick_mdfcurve = (floatvec **)calloc(sizeof(floatvec *),
685 /* 22 Oct 2008 */ new_nvals) ;
686 }
687
688 /*** loop over input datasets ***/
689
690 if( ninp > 1 ) myRwcFree( new_dset->keywords ) ;
691
692 ivout = 0 ;
693 for( ids=0 ; ids < ninp ; ids++ ){
694 dset = DSUB(ids) ;
695 nv = NSUBV(ids) ;
696
697 if( ! BUCK_dry ){
698 DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
699 }
700 /** loop over sub-bricks to output **/
701
702 for( iv=0 ; iv < nv ; iv++ ){
703 jv = SUBV(ids,iv) ; /* which sub-brick to use */
704
705 if( ! BUCK_dry ){
706 EDIT_substitute_brick( new_dset , ivout ,
707 DSET_BRICK_TYPE(dset,jv) , DSET_ARRAY(dset,jv) ) ;
708
709 /*----- preserve label when one exists --- Modified March 2010 ZSS*/
710 if (DSET_HAS_LABEL(dset, jv) )
711 sprintf (buf, "%s", DSET_BRICK_LABEL(dset,jv));
712 else
713 sprintf(buf,"%.12s[%d]",DSET_PREFIX(dset),jv) ;
714 EDIT_dset_items( new_dset , ADN_brick_label_one+ivout, buf , ADN_none ) ;
715
716 #if 0
717 sprintf(buf,"%s[%d]",DSET_FILECODE(dset),jv) ;
718 EDIT_dset_items(
719 new_dset, ADN_brick_keywords_replace_one+ivout, buf, ADN_none ) ;
720 #endif
721
722 EDIT_dset_items(
723 new_dset ,
724 ADN_brick_fac_one +ivout, DSET_BRICK_FACTOR(dset,jv),
725 #if 0
726 ADN_brick_keywords_append_one+ivout, DSET_BRICK_KEYWORDS(dset,jv) ,
727 #endif
728 ADN_none ) ;
729
730 /** possibly write statistical parameters for this sub-brick **/
731
732 kv = DSET_BRICK_STATCODE(dset,jv) ;
733
734 if( FUNC_IS_STAT(kv) ){ /* input sub-brick has stat params */
735
736 int npar = FUNC_need_stat_aux[kv] , lv ;
737 float * par = (float *) malloc( sizeof(float) * (npar+2) ) ;
738 float * sax = DSET_BRICK_STATAUX(dset,jv) ;
739 par[0] = kv ;
740 par[1] = npar ;
741 for( lv=0 ; lv < npar ; lv++ )
742 par[lv+2] = (sax != NULL) ? sax[lv] : 0.0 ;
743
744 EDIT_dset_items(new_dset ,
745 ADN_brick_stataux_one+ivout , par ,
746 ADN_none ) ;
747 free(par) ;
748
749 /* 2: if the input dataset has statistical parameters */
750
751 } else if( ISFUNC(dset) && /* dset has stat */
752 FUNC_IS_STAT(dset->func_type) && /* params */
753 jv == FUNC_ival_thr[dset->func_type] ){ /* thr sub-brick */
754
755 int npar , lv ;
756 float * par , * sax ;
757 kv = dset->func_type ;
758 npar = FUNC_need_stat_aux[kv] ;
759 par = (float *) malloc( sizeof(float) * (npar+2) ) ;
760 sax = dset->stat_aux ;
761 par[0] = kv ;
762 par[1] = npar ;
763 for( lv=0 ; lv < npar ; lv++ )
764 par[lv+2] = (sax != NULL) ? sax[lv] : 0.0 ;
765
766 EDIT_dset_items(new_dset ,
767 ADN_brick_stataux_one+ivout , par ,
768 ADN_none ) ;
769 free(par) ;
770 }
771
772 /** append any fdrcurve **/
773 if( have_fdr ){
774 /* fixed iv->jv (ick!), noticed by dglen 16 Mar 2010 [rickr] */
775 if(dset->dblk->brick_fdrcurve && dset->dblk->brick_fdrcurve[jv]){
776 COPY_floatvec(new_dset->dblk->brick_fdrcurve[ivout],
777 dset->dblk->brick_fdrcurve[jv]) ;
778 nfdr++;
779 }
780 else new_dset->dblk->brick_fdrcurve[ivout] = NULL ;
781
782 if(dset->dblk->brick_mdfcurve && dset->dblk->brick_mdfcurve[jv]){
783 COPY_floatvec(new_dset->dblk->brick_mdfcurve[ivout],
784 dset->dblk->brick_mdfcurve[jv]) ;
785 }
786 else new_dset->dblk->brick_mdfcurve[ivout] = NULL ;
787 }
788
789 /** print a message? **/
790
791 if( BUCK_verb ) printf("-verb: copied %s[%d] into %s[%d]\n" ,
792 DSET_FILECODE(dset) , jv ,
793 DSET_FILECODE(new_dset) , ivout ) ;
794 } else {
795 printf("-dry: would copy %s[%d] into %s[%d]\n" ,
796 DSET_FILECODE(dset) , jv ,
797 DSET_FILECODE(new_dset) , ivout ) ;
798 }
799
800 ivout++ ;
801 }
802
803 /** loop over all bricks in input dataset and
804 unload them if they aren't going into the output
805 (not required, but is done to economize on memory) **/
806
807 if( ! BUCK_dry && nv < DSET_NVALS(dset) ){
808
809 for( kv=0 ; kv < DSET_NVALS(dset) ; kv++ ){ /* all input sub-bricks */
810 for( iv=0 ; iv < nv ; iv++ ){ /* all output sub-bricks */
811 jv = SUBV(ids,iv) ;
812 if( jv == kv ) break ; /* input matches output */
813 }
814 if( iv == nv ){
815 mri_free( DSET_BRICK(dset,kv) ) ;
816 #if 0
817 if( BUCK_verb ) printf("-verb: unloaded unused %s[%d]\n" ,
818 DSET_FILECODE(dset) , kv ) ;
819 #endif
820 }
821 }
822 }
823
824 } /* end of loop over input datasets */
825
826 if( ! BUCK_dry ){
827 if( BUCK_verb ){
828 if( have_fdr ) fprintf(stderr,"-verb: added %d of %d fdr curves\n",
829 nfdr, new_nvals);
830 fprintf(stderr,"-verb: loading statistics\n") ;
831 }
832 THD_load_statistics( new_dset ) ;
833 if( BUCK_glue ) putenv("AFNI_DECONFLICT=OVERWRITE") ;
834 if( BUCK_glue && BUCK_ccode >= 0 )
835 THD_set_write_compression(BUCK_ccode) ; /* 16 Mar 2010 */
836 if ( ! THD_write_3dim_dataset( NULL,NULL , new_dset , True ) )
837 exit(1) ;
838 if( BUCK_verb ) fprintf(stderr,"-verb: wrote output: %s\n",DSET_BRIKNAME(new_dset)) ;
839 }
840
841 exit(0) ;
842 }
843