1 #include "mrilib.h"
2 #include "thd.h"
3 
4 /*---------------------------------------------------------------------------
5   This program catenates multiple 3D datasets in various directions.
6 ----------------------------------------------------------------------------*/
7 
8 /*-------------------------- global data --------------------------*/
9 
10 static THD_3dim_dataset_array *XCAT_dsar  = NULL ;  /* input datasets */
11 static int                     XCAT_dir   = 1 ;     /* 1=x, 2=y, 3=z */
12 static int                     XCAT_nbrik = -1 ;    /* # bricks */
13 static int                     XCAT_verb  = 0 ;     /* verbose? */
14 static int                     XCAT_datum = -1 ;    /* dataset datum */
15 
16 #define DSUB(id) DSET_IN_3DARR(XCAT_dsar,(id))
17 
18 static char XCAT_output_prefix[THD_MAX_PREFIX] = "xyzcat" ;
19 
20 /*--------------------------- prototypes ---------------------------*/
21 
22 void XCAT_read_opts( int , char ** ) ;
23 void XCAT_Syntax(void) ;
24 
25 /*--------------------------------------------------------------------
26    read the arguments, load the global variables
27 ----------------------------------------------------------------------*/
28 
XCAT_read_opts(int argc,char * argv[])29 void XCAT_read_opts( int argc , char *argv[] )
30 {
31    int nopt=1 , ii , nerr=0 ;
32    THD_3dim_dataset * dset ;
33 
34    INIT_3DARR(XCAT_dsar) ;  /* array of datasets */
35 
36    while( nopt < argc ){
37 
38       /**** -prefix prefix ****/
39 
40       if( strncmp(argv[nopt],"-prefix",6) == 0 ){
41         if( ++nopt >= argc ) ERROR_exit("need argument after -prefix!\n") ;
42 
43         MCW_strncpy( XCAT_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
44 
45         if( !THD_filename_ok(XCAT_output_prefix) )
46           ERROR_exit("Illegal character in -prefix '%s'",XCAT_output_prefix) ;
47 
48         continue ;
49       }
50 
51       /**** -verb ****/
52 
53       if( strncmp(argv[nopt],"-verb",5) == 0 ){
54         XCAT_verb++ ; nopt++ ; continue ;
55       }
56 
57       /*** -dir ***/
58 
59       if( strncmp(argv[nopt],"-dir",4) == 0 ){
60         nopt++ ;
61         if( nopt >= argc ) ERROR_exit("need argument after -dir!") ;
62         switch( toupper(argv[nopt][0]) ){
63           default: ERROR_exit("don't understand '-dir %s'",argv[nopt]) ;
64           case 'X': case 'I': XCAT_dir = 1 ; break ;
65           case 'Y': case 'J': XCAT_dir = 2 ; break ;
66           case 'Z': case 'K': XCAT_dir = 3 ; break ;
67         }
68         nopt++ ; continue ;
69       }
70 
71       /**** Garbage of various flavors ****/
72 
73       if( strcmp(argv[nopt],"-input") == 0 ){
74         nopt++ ;  /* and fall thru */
75       } else if( argv[nopt][0] == '-' ){
76         ERROR_exit("Unknown option: %s\n",argv[nopt]) ;
77       }
78 
79       /**** read dataset ****/
80 
81       if( XCAT_verb ) INFO_message("Opening dataset %s",argv[nopt]) ;
82 
83       dset = THD_open_dataset( argv[nopt++] ) ;
84       if( dset == NULL ){
85         ERROR_message("Can't open dataset %s",argv[nopt-1]) ; nerr++ ; continue ;
86       }
87       THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
88 
89       if( XCAT_datum < 0 ) XCAT_datum = DSET_BRICK_TYPE(dset,0) ;
90       if( XCAT_nbrik < 0 ) XCAT_nbrik = DSET_NVALS(dset) ;
91 
92       /* check stuff */
93 
94       if( !DSET_datum_constant(dset) ){
95         ERROR_message("Dataset %s doesn't have a constant data type :-(",argv[nopt-1]) ;
96         nerr++ ;
97       }
98       if( DSET_NVALS(dset) != XCAT_nbrik ){
99         ERROR_message("Dataset %s has different number of sub-bricks",argv[nopt-1]) ;
100         nerr++ ;
101       }
102       if( DSET_BRICK_TYPE(dset,0) != XCAT_datum ){
103         ERROR_message("Dataset %s has a different data type",argv[nopt-1]) ;
104         nerr++ ;
105       }
106       if( THD_need_brick_factor(dset) ){
107         ERROR_message("Can't catenated dataset %s that uses brick factors - sorry!",argv[nopt-1]) ;
108         nerr++ ;
109       }
110 
111       ADDTO_3DARR(XCAT_dsar,dset) ;  /* list of datasets */
112 
113    }  /* end of loop over command line arguments */
114 
115    if( nerr > 0 ) ERROR_exit("Can't continue after above ERROR message%s",
116                              (nerr==1) ? "\0" : "s" ) ;
117    return ;
118 }
119 
120 /*-------------------------------------------------------------------------*/
121 
XCAT_Syntax(void)122 void XCAT_Syntax(void)
123 {
124    printf(
125     "Usage: 3dXYZcat [options] dataset dataset ...\n"
126     "* Catenates datasets spatially (for time cat-ing, cf. 3dTcat).\n"
127     "* The input datasets must match, in the sense that the pieces\n"
128     "  fit together properly (spatially and in time).\n"
129     "* Unlike in 3dZcat, all input datasets must be stored with the\n"
130     "  same data type (e.g., shorts, floats, ...); also, sub-brick scale\n"
131     "  factors are not allowed.  If you need to spatially catenate scaled\n"
132     "  short datasets, for example, convert them to float format using\n"
133     "  '3dcalc -float', then catenate THOSE datasets.\n"
134     "\n"
135     "Options:\n"
136     "--------\n"
137     "  -prefix pname = Use 'pname' for the output dataset prefix name.\n"
138     "                    [default prefix = 'xyzcat']\n"
139     "  -verb         = Print out some verbositiness as the program proceeds.\n"
140     "  -dir Q        = Catenate along direction 'Q', which is one of\n"
141     "                    X or Y or Z (synonyms are I or J or K)\n"
142     "                  which are the STORAGE directions (not DICOM) of the\n"
143     "                  3D grid of the input datasets.\n"
144     "                    [default direction = 'X', for no good reason]\n"
145     "\n"
146     "Command line arguments after the above are taken as input datasets.\n"
147     "\n"
148     "Notes:\n"
149     "------\n"
150     "* If the i-th input dataset has dimensions nx[i] X ny[i] X nz[i], then\n"
151     "    case Q = X | I ==> all ny[i] and nz[i] must be the same;\n"
152     "                       the output dataset has nx = sum{ nx[i] }\n"
153     "    case Q = Y | J ==> all nx[i] and nz[i] must be the same;\n"
154     "                       the output dataset has ny = sum{ ny[i] }\n"
155     "    case Q = Z | K ==> all nx[i] and ny[i] must be the same;\n"
156     "                       the output dataset has nz = sum{ nz[i] }\n"
157     "* In all cases, the input datasets must have the same number of\n"
158     "  sub-bricks (time points) and the same data storage type.\n"
159     "* You can use the '3dinfo' program to see the orientation and\n"
160     "    grid size of a dataset, to help you decide how to glue your\n"
161     "    inputs together.\n"
162     "* There must be at least two datasets input (otherwise, the\n"
163     "    program doesn't make much sense, now does it?).\n"
164     "* This is mostly useful for making side-by-side pictures from\n"
165     "    multiple datasets, for edification and elucidation.\n"
166     "* If you have some other use for 3dXYZcat, let me know!\n"
167     "\n"
168     "** Author: RW Cox [Dec 2010] **\n"
169    ) ;
170 
171    PRINT_COMPILE_DATE ; exit(0) ;
172 }
173 
174 /*-------------------------------------------------------------------------*/
175 
main(int argc,char * argv[])176 int main( int argc , char *argv[] )
177 {
178    int ninp , ids , iv ,kz,new_nz, nx,ny,nz,nxy,nxyz ;
179    THD_3dim_dataset *new_dset=NULL , *dset=NULL ;
180    MRI_IMAGE *outim;
181    MRI_IMARR *im_array;
182 
183    /*** read input options ***/
184 
185    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) XCAT_Syntax() ;
186 
187    /*-- addto the arglist, if user wants to --*/
188 
189    { int new_argc ; char ** new_argv ;
190      addto_args( argc , argv , &new_argc, &new_argv ) ;
191      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
192    }
193 
194    mainENTRY("3dZcat main") ; machdep() ; AFNI_logger("3dZcat",argc,argv) ;
195    PRINT_VERSION("3dZcat") ;
196 
197    XCAT_read_opts( argc , argv ) ;
198 
199    /*** create new dataset (empty) ***/
200 
201    ninp = XCAT_dsar->num ;
202    if( ninp < 2 )
203      ERROR_exit("Must have at least 2 input datasets!") ;
204 
205    for( ids=0 ; ids < ninp ; ids++ ){
206      dset = DSUB(ids) ; if( DSET_TIMESTEP(dset) > 0.0f ) break ;
207    }
208    if( ids == ninp ) dset = DSUB(0) ; /* fallback position */
209    if( XCAT_verb )
210      INFO_message("Using %s as 'master' dataset",DSET_HEADNAME(dset)) ;
211 
212    new_dset = EDIT_empty_copy( dset ) ; /* make a copy of its header */
213 
214    tross_Copy_History( dset , new_dset ) ;
215    tross_Make_History( "3dXYZcat" , argc,argv , new_dset ) ;
216 
217    EDIT_dset_items( new_dset ,
218                       ADN_prefix    , XCAT_output_prefix ,
219                       ADN_datum_all , XCAT_datum ,
220                       ADN_nsl       , 0 , /* kill time offsets  */
221                     ADN_none ) ;
222 
223    /* can't re-write existing dataset */
224 
225    if( THD_deathcon() && THD_is_file(DSET_HEADNAME(new_dset)) )
226      ERROR_exit("Output dataset %s already exists!",DSET_HEADNAME(new_dset) ) ;
227 
228    THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
229 
230    if( XCAT_verb ) INFO_message("Loading input datasets") ;
231    for( ids=0 ; ids < ninp ; ids++ ){
232      dset = DSUB(ids) ; DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
233    }
234 
235    /*** Loop over output sub-bricks ***/
236 
237    for( iv=0 ; iv < DSET_NVALS(new_dset) ; iv++ ){
238 
239       if( XCAT_verb ) ININFO_message("Computing output sub-brick #%d",iv) ;
240 
241       INIT_IMARR(im_array) ;
242       for( ids=0 ; ids < ninp ; ids++ )
243         ADDTO_IMARR( im_array , DSET_BRICK(DSUB(ids),iv) ) ;
244 
245       outim = mri_catvol_1D(im_array,XCAT_dir) ;
246       if( outim == NULL )
247         ERROR_exit("Can't catenate! Dataset dimensions mismatch? :-(") ;
248 
249       FREE_IMARR(im_array) ;
250       for( ids=0 ; ids < ninp ; ids++ ) DSET_unload_one(DSUB(ids),iv) ;
251 
252       if( iv == 0 ){  /* mangle sub-brick xyz dimensions */
253         THD_ivec3 iv_nxyz ;
254         LOAD_IVEC3( iv_nxyz , outim->nx , outim->ny , outim->nz ) ;
255         EDIT_dset_items( new_dset , ADN_nxyz , iv_nxyz , ADN_none ) ;
256       }
257 
258       EDIT_substitute_brick( new_dset,iv, XCAT_datum,mri_data_pointer(outim) ) ;
259       mri_clear_data_pointer(outim) ; mri_free(outim) ;
260 
261    } /* end of loop over output sub-bricks */
262 
263    /*-- write dataset --*/
264 
265    for( ids=0 ; ids < ninp ; ids++ ) DSET_delete( DSUB(ids) ) ;
266    DSET_write(new_dset) ;
267    WROTE_DSET(new_dset) ;
268    exit(0) ;
269 }
270