1 #include "mrilib.h"
2 #include "thd.h"
3
4 /*---------------------------------------------------------------------------
5 This program catenates multiple 3D datasets in the slice direction.
6 -- RWCox - 08 Aug 2001
7 ----------------------------------------------------------------------------*/
8
9 /*-------------------------- global data --------------------------*/
10
11 static THD_3dim_dataset_array * ZCAT_dsar = NULL ; /* input datasets */
12 static int ZCAT_nxy = -1 ; /* # voxels/slice */
13 static int ZCAT_nbrik = -1 ; /* # bricks */
14 static int ZCAT_verb = 0 ; /* verbose? */
15 static int ZCAT_type = -1 ; /* dataset type */
16 static int ZCAT_datum = -1 ; /* dataset datum */
17 static int ZCAT_fscale = 0 ;
18 static int ZCAT_gscale = 0 ;
19 static int ZCAT_nscale = 0 ;
20 static int ZCAT_frugal = 0 ; /* 05 Apr 2006 */
21
22 #define DSUB(id) DSET_IN_3DARR(ZCAT_dsar,(id))
23
24 static char ZCAT_output_prefix[THD_MAX_PREFIX] = "zcat" ;
25 static RwcBoolean write_output = False; /* 21 Jun 2006 [dg] -force rewrite as in 3drefit by rickr */
26 static RwcBoolean NIFTI_mode = False; /* saving NIFTI output */
27 static int cmode = COMPRESS_NOFILE; /* check compression mode for NIFTI separately */
28
29 /*--------------------------- prototypes ---------------------------*/
30
31 void ZCAT_read_opts( int , char ** ) ;
32 void ZCAT_Syntax(void) ;
33
34 /*--------------------------------------------------------------------
35 read the arguments, load the global variables
36 ----------------------------------------------------------------------*/
37
ZCAT_read_opts(int argc,char * argv[])38 void ZCAT_read_opts( int argc , char * argv[] )
39 {
40 int nopt = 1 , ii ;
41 THD_3dim_dataset * dset ;
42
43 INIT_3DARR(ZCAT_dsar) ; /* array of datasets */
44
45 while( nopt < argc ){
46
47 /**** -frugal [05 Apr 2006] ****/
48
49 if( strncmp(argv[nopt],"-frugal",4) == 0 ){
50 ZCAT_frugal = 1 ; nopt++ ; continue ;
51 }
52
53 /**** -nscale ****/
54
55 if( strncmp(argv[nopt],"-nscale",6) == 0 ){
56 ZCAT_gscale = ZCAT_fscale = 0 ; ZCAT_nscale = 1 ; nopt++ ; continue ;
57 }
58
59 /**** -fscale ****/
60
61 if( strncmp(argv[nopt],"-fscale",6) == 0 ){
62 ZCAT_fscale = 1 ; ZCAT_nscale = 0 ; nopt++ ; continue ;
63 }
64
65 #if 0
66 /**** -gscale ****/
67
68 if( strncmp(argv[nopt],"-gscale",6) == 0 ){
69 ZCAT_gscale = ZCAT_fscale = 1 ; ZCAT_nscale = 0 ; nopt++ ; continue ;
70 }
71 #endif
72
73 /**** -prefix prefix ****/
74
75 if( strncmp(argv[nopt],"-prefix",6) == 0 ){
76 nopt++ ;
77 if( nopt >= argc ){
78 ERROR_exit("need argument after -prefix!\n") ;
79 }
80 MCW_strncpy( ZCAT_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
81
82 if( strstr(ZCAT_output_prefix,".nii") != NULL ) {
83 write_output = True;
84 NIFTI_mode = True;
85 if( strstr(ZCAT_output_prefix,".nii.gz") != NULL ) {
86 cmode = 0; /* force gzip compression (actually zlib from nifti library)*/
87 }
88 }
89 else if( !THD_filename_ok(ZCAT_output_prefix) )
90 ERROR_exit("Illegal character in -prefix '%s'",ZCAT_output_prefix) ;
91 continue ;
92 }
93
94 /**** -verb ****/
95
96 if( strncmp(argv[nopt],"-verb",5) == 0 ){
97 ZCAT_verb++ ; nopt++ ; continue ;
98 }
99
100 /**** -datum type ****/
101
102 if( strncmp(argv[nopt],"-datum",6) == 0 ){
103 if( ++nopt >= argc ){
104 ERROR_exit("need an argument after -datum!\n");
105 }
106 if( strcmp(argv[nopt],"short") == 0 ) ZCAT_datum = MRI_short ;
107 else if( strcmp(argv[nopt],"float") == 0 ) ZCAT_datum = MRI_float ;
108 else if( strcmp(argv[nopt],"byte" ) == 0 ) ZCAT_datum = MRI_byte ;
109 #if 0
110 else if( strcmp(argv[nopt],"complex")== 0) ZCAT_datum = MRI_complex ;
111 #endif
112 else {
113 ERROR_exit("-datum %s not supported in 3dZcat!\n",
114 argv[nopt] ) ;
115 }
116 nopt++ ; continue ; /* go to next arg */
117 }
118
119 /**** Garbage ****/
120
121 if( strcmp(argv[nopt],"-input") == 0 ){
122 nopt++ ;
123 } else if( argv[nopt][0] == '-' ){
124 ERROR_exit("Unknown option: %s\n",argv[nopt]) ;
125 }
126
127 /**** read dataset ****/
128
129 if( ZCAT_verb )
130 fprintf(stderr,"++ Opening dataset %s\n",argv[nopt]) ;
131
132 dset = THD_open_dataset( argv[nopt++] ) ;
133 if( dset == NULL ){
134 ERROR_exit("Can't open dataset %s\n",argv[nopt-1]) ;
135 }
136 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
137 THD_make_cardinal(dset); /* deoblique 21 Oct, 2011 [rickr] */
138
139 if( ZCAT_type < 0 ) ZCAT_type = dset->type ;
140
141 /* check if voxel counts match, etc. */
142
143 ii = dset->daxes->nxx * dset->daxes->nyy ;
144 if( ZCAT_nxy < 0 ){
145 ZCAT_nxy = ii ;
146 ZCAT_nbrik = DSET_NVALS(dset) ;
147 } else if( ii != ZCAT_nxy ){
148 ERROR_exit("** Dataset %s differs in slice size from others\n",argv[nopt-1]);
149 } else if ( DSET_NVALS(dset) != ZCAT_nbrik ){
150 ERROR_exit("** Dataset %s has different number of sub-bricks\n",argv[nopt-1]) ;
151 } else if ( DSET_BRICK_TYPE(dset,0) == MRI_complex ){
152 ERROR_exit("** Dataset %s is complex-valued -- ILLEGAL\n",argv[nopt-1]) ;
153 }
154
155 ADDTO_3DARR(ZCAT_dsar,dset) ; /* list of datasets */
156
157 } /* end of loop over command line arguments */
158
159 if(NIFTI_mode && ZCAT_frugal){
160 ERROR_message("Frugality and NIFTI output do not mix.\n");
161 ERROR_exit("Try without -frugal or with different output type.\n");
162 }
163 return ;
164 }
165
166 /*-------------------------------------------------------------------------*/
167
ZCAT_Syntax(void)168 void ZCAT_Syntax(void)
169 {
170 printf(
171 "Usage: 3dZcat [options] dataset dataset ...\n"
172 "Concatenates datasets in the slice (z) direction. Each input\n"
173 "dataset must have the same number of voxels in each slice, and\n"
174 "must have the same number of sub-bricks.\n"
175 "\n"
176 "Options:\n"
177 " -prefix pname = Use 'pname' for the output dataset prefix name.\n"
178 " [default='zcat']\n"
179 " -datum type = Coerce the output data to be stored as the given\n"
180 " type, which may be byte, short, or float.\n"
181 " -fscale = Force scaling of the output to the maximum integer\n"
182 " range. This only has effect if the output datum\n"
183 " is byte or short (either forced or defaulted).\n"
184 " This option is sometimes necessary to eliminate\n"
185 " unpleasant truncation artifacts.\n"
186 " -nscale = Don't do any scaling on output to byte or short datasets.\n"
187 " This may be especially useful when operating on mask\n"
188 " datasets whose output values are only 0's and 1's.\n"
189 " -verb = Print out some verbositiness as the program proceeds.\n"
190 " -frugal = Be 'frugal' in the use of memory, at the cost of I/O time.\n"
191 " Only needed if the program runs out of memory.\n"
192 " Note frugality cannot be combined with NIFTI output\n"
193 "\n"
194 "Command line arguments after the above are taken as input datasets.\n"
195 "\n"
196 "Notes:\n"
197 "* You can use the '3dinfo' program to see how many slices a\n"
198 " dataset comprises.\n"
199 "* There must be at least two datasets input (otherwise, the\n"
200 " program doesn't make much sense, does it?).\n"
201 "* Each input dataset must have the same number of voxels in each\n"
202 " slice, and must have the same number of sub-bricks.\n"
203 "* This program does not deal with complex-valued datasets.\n"
204 "* See the output of '3dZcutup -help' for a C shell script that\n"
205 " can be used to take a dataset apart into single slice datasets,\n"
206 " analyze them separately, and then assemble the results into\n"
207 " new 3D datasets.\n"
208 "* Also see program 3dXYZcat for a version that can catenate along\n"
209 " the x and y axes as well (with some limitations).\n"
210 ) ;
211
212 PRINT_COMPILE_DATE ; exit(0) ;
213 }
214
215 /*-------------------------------------------------------------------------*/
216
main(int argc,char * argv[])217 int main( int argc , char * argv[] )
218 {
219 int ninp , ids , iv ,kz,new_nz, nx,ny,nz,nxy,nxyz ;
220 THD_3dim_dataset * new_dset=NULL , * dset=NULL ;
221 THD_ivec3 iv_nxyz ;
222 float * fvol , *ffac ;
223 void * svol ;
224 int fscale ; FILE * data_file = NULL ;
225 MRI_IMAGE *svol_im;
226 MRI_IMARR *im_array=NULL;
227
228 /*** read input options ***/
229
230 if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) ZCAT_Syntax() ;
231
232 /*-- addto the arglist, if user wants to --*/
233
234 { int new_argc ; char ** new_argv ;
235 addto_args( argc , argv , &new_argc, &new_argv ) ;
236 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
237 }
238
239 mainENTRY("3dZcat main") ; machdep() ; AFNI_logger("3dZcat",argc,argv) ;
240 PRINT_VERSION("3dZcat") ;
241
242 ZCAT_read_opts( argc , argv ) ;
243
244 /*** create new dataset (empty) ***/
245
246 ninp = ZCAT_dsar->num ;
247 if( ninp < 2 ){
248 ERROR_exit("Must have at least 2 input datasets!\n");
249 }
250
251 /* compute total number of z-slices in output */
252 new_nz = 0 ;
253 for( ids=0 ; ids < ninp ; ids++ ) new_nz += DSET_NZ(DSUB(ids)) ;
254
255 if( ZCAT_verb )
256 fprintf(stderr,"++ Output will have %d slices\n",new_nz) ;
257
258 /** find 1st dataset that is time dependent, if any **/
259
260 for( ids=0 ; ids < ninp ; ids++ ){
261 dset = DSUB(ids) ;
262 if( DSET_TIMESTEP(dset) > 0.0 ) break ;
263 }
264 if( ids == ninp ) dset = DSUB(0) ; /* fallback position */
265 if( ZCAT_verb )
266 fprintf(stderr,"++ Using %s as 'master' dataset\n",DSET_HEADNAME(dset)) ;
267
268 new_dset = EDIT_empty_copy( dset ) ; /* make a copy of its header */
269
270 tross_Copy_History( dset , new_dset ) ;
271 tross_Make_History( "3dZcat" , argc,argv , new_dset ) ;
272
273 /* modify its header */
274
275 nx = DSET_NX(dset) ;
276 ny = DSET_NY(dset) ; nxy = nx*ny ; nxyz = nxy*new_nz ;
277
278 LOAD_IVEC3( iv_nxyz , nx , ny , new_nz ) ;
279
280 /* if data type has not been selected by user, get the data type from master */
281 if( ZCAT_datum < 0 ) ZCAT_datum = DSET_BRICK_TYPE(dset,0) ;
282
283 /*-- open output BRIK file --*/
284 /* default cmode is COMPRESS_NOFILE unless nifti.gzip output (.nii.gz) */
285 if (cmode == COMPRESS_NOFILE) { /* ignore compression for NIFTI - do in write
286 automatically later */
287 cmode = THD_get_write_compression() ; /* check env. variable for compression*/
288 #if 0
289 if(NIFTI_mode && (cmode!=0)) /* have to compress this NIFTI data, add .gz to prefix */
290 cmode = COMPRESS_NOFILE;
291 sprintf(ZCAT_output_prefix, "%s.gz", ZCAT_output_prefix);
292 #endif
293 }
294
295 EDIT_dset_items( new_dset ,
296 ADN_prefix , ZCAT_output_prefix ,
297 ADN_type , ZCAT_type ,
298 ADN_nxyz , iv_nxyz ,
299 ADN_datum_all , ZCAT_datum ,
300 ADN_nsl , 0 , /* kill time offsets */
301 ADN_none ) ;
302
303 /* can't re-write existing dataset */
304 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(new_dset)) ){
305 ERROR_exit("Fatal error: dataset %s already exists!\n",
306 DSET_HEADNAME(new_dset) ) ;
307 }
308
309 THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
310
311 if( !ZCAT_frugal ){ /* 05 Apr 2006 */
312 for( ids=0 ; ids < ninp ; ids++ ){
313 dset = DSUB(ids) ; DSET_load(dset) ; CHECK_LOAD_ERROR(dset) ;
314 }
315 }
316 else {
317 data_file = COMPRESS_fopen_write( DSET_BRIKNAME(new_dset) , cmode ) ;
318 if( data_file == NULL ){
319 ERROR_exit(
320 "\a\n*** cannot open output file %s\n",DSET_BRIKNAME(new_dset)) ;
321 }
322 }
323
324 /* make space for sub-brick scaling factors */
325
326 ffac = (float *) malloc(sizeof(float)*DSET_NVALS(new_dset)) ;
327
328 if (!ZCAT_frugal)
329 INIT_IMARR(im_array);
330
331
332 /*** Loop over output sub-bricks ***/
333
334 for( iv=0 ; iv < DSET_NVALS(new_dset) ; iv++ ){
335
336 if( ZCAT_verb )
337 fprintf(stderr,"++ Computing output sub-brick #%d\n",iv) ;
338
339 fscale = ZCAT_fscale ; /* local for this brick */
340
341 /* make temporary holding space */
342
343 fvol = (float *) malloc(sizeof(float)*nxyz) ;
344
345 /* loop over input datasets */
346 /* kz = slice index in output that this input dataset starts at */
347
348 for( kz=ids=0 ; ids < ninp ; ids++ ){
349
350 dset = DSUB(ids) ; nz = DSET_NZ(dset) ;
351 if( ZCAT_frugal ){
352 DSET_load(dset) ;
353 if( ! DSET_LOADED(dset) ){
354 ERROR_message("Fatal error: can't load data from %s\n",
355 DSET_BRIKNAME(dset)) ;
356 COMPRESS_fclose(data_file) ; remove(DSET_BRIKNAME(new_dset)) ;
357 exit(1) ;
358 }
359 }
360 if( ZCAT_verb == 2 )
361 fprintf(stderr," + processing input %s\n",DSET_BRIKNAME(dset)) ;
362
363 /* copy data from input to holding space, converting to float */
364
365 switch( DSET_BRICK_TYPE(dset,iv) ){
366
367 default:
368 ERROR_exit(
369 "** Illegal input brick type=%s in dataset %s\n",
370 MRI_TYPE_name[DSET_BRICK_TYPE(dset,iv)] ,
371 DSET_HEADNAME(dset) ) ;
372
373 case MRI_short:{
374 register int ii , itop = DSET_NVOX(dset) ;
375 register short *dar = DSET_ARRAY(dset,iv) ;
376 register float fac = DSET_BRICK_FACTOR(dset,iv) ;
377 register float *var = fvol + kz*nxy ;
378 if( fac == 0.0 ) fac = 1.0 ;
379 for( ii=0 ; ii < itop ; ii++ ) var[ii] = fac*dar[ii] ;
380 if( fac != 1.0 ) fscale = 1 ;
381 }
382 break ;
383
384 case MRI_byte:{
385 register int ii , itop = DSET_NVOX(dset) ;
386 register byte *dar = DSET_ARRAY(dset,iv) ;
387 register float fac = DSET_BRICK_FACTOR(dset,iv) ;
388 register float *var = fvol + kz*nxy ;
389 if( fac == 0.0 ) fac = 1.0 ;
390 for( ii=0 ; ii < itop ; ii++ ) var[ii] = fac*dar[ii] ;
391 if( fac != 1.0 ) fscale = 1 ;
392 }
393 break ;
394
395 case MRI_float:{
396 register int ii , itop = DSET_NVOX(dset) ;
397 register float *dar = DSET_ARRAY(dset,iv) ;
398 register float fac = DSET_BRICK_FACTOR(dset,iv) ;
399 register float *var = fvol + kz*nxy ;
400 if( fac == 0.0 ) fac = 1.0 ;
401 for( ii=0 ; ii < itop ; ii++ ) var[ii] = fac*dar[ii] ;
402 }
403 break ;
404
405 } /* end of switch over input brick type */
406
407 kz += nz ;
408 if( ZCAT_frugal ) DSET_unload(dset) ;
409 else DSET_unload_one(dset,iv) ;
410
411 } /* end of loop over input datasets */
412
413 /* at this point, fvol holds the correctly
414 scaled values for the output dataset;
415 now, must store it in the correct datum into svol;
416 this code is lifted/adapted from 3dcalc.c */
417
418 switch( ZCAT_datum ){
419
420 default: /* should never transpire */
421 ERROR_exit(
422 "** Fatal Error **\n"
423 "** Somehow ended up with -datum = %d\n",ZCAT_datum) ;
424
425 case MRI_float: /* the trivial case */
426 svol = (void *) fvol ;
427 ffac[iv] = 0.0 ; /* don't need no stinking factor */
428 break ;
429
430 case MRI_byte: /* harder cases: */
431 case MRI_short:{ /* must create svol and scale to it */
432 float gtop , fimfac ;
433
434 /* find largest value in fvol */
435
436 gtop = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, fvol ) ;
437 if( gtop == 0.0 )
438 fprintf(stderr,"++ Warning: output sub-brick %d is all zeros!\n",iv) ;
439
440 /* compute scaling factor */
441
442 if( fscale ){
443 fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[ZCAT_datum] / gtop : 0.0 ;
444 } else if( !ZCAT_nscale ){
445 fimfac = (gtop > MRI_TYPE_maxval[ZCAT_datum] || (gtop > 0.0 && gtop <= 1.0) )
446 ? MRI_TYPE_maxval[ZCAT_datum]/ gtop : 0.0 ;
447 } else {
448 fimfac = 0.0 ;
449 }
450
451 if( ZCAT_verb == 2 ){
452 if( fimfac != 0.0 )
453 fprintf(stderr," + Output sub-brick %d scale factor = %f\n",iv,fimfac) ;
454 else
455 fprintf(stderr,"++ Output sub-brick %d: no scale factor\n" ,iv) ;
456 }
457
458 /* now scale fvol into svol */
459
460 svol = (void *) malloc( mri_datum_size(ZCAT_datum) * nxyz ) ;
461
462 EDIT_coerce_scale_type( nxyz , fimfac ,
463 MRI_float, fvol , ZCAT_datum,svol ) ;
464
465 free(fvol) ;
466 ffac[iv] = (fimfac > 0.0 && fimfac != 1.0) ? 1.0/fimfac : 0.0 ;
467 }
468 break ;
469 } /* end of switch on output datum type */
470
471 /* if being frugal with memory (why?),
472 save the data to disk for each output sub-brick */
473 if(ZCAT_frugal) {
474 /*-- now save svol to disk --*/
475 fwrite( svol , mri_datum_size(ZCAT_datum) , nxyz , data_file ) ;
476 free(svol) ;
477 }
478 /* otherwise, save them up and write all sub-bricks all at once */
479 else {
480 /* copy the sub-brick volume, svol, into an MRI_IMAGE structure and
481 then append that to the image array */
482 svol_im = mri_new_vol_empty( nx, ny, kz, ZCAT_datum) ;
483 mri_fix_data_pointer( svol , svol_im ) ;
484 ADDTO_IMARR(im_array, svol_im);
485 }
486
487 } /* end of loop over output sub-bricks */
488
489
490 /*-- cleanup, and write dataset header (and brick if not being frugal) --*/
491
492 if( ZCAT_verb )
493 fprintf(stderr,"++ Computing output sub-brick statistics\n") ;
494
495 for( ids=0 ; ids < ninp ; ids++ ) /* remove inputs from memory */
496 DSET_delete( DSUB(ids) ) ;
497
498 /* update output dataset image array with pointer to all of the data */
499 if(!ZCAT_frugal)
500 new_dset->dblk->brick = im_array; /* update pointer to data */
501 else
502 COMPRESS_fclose(data_file) ; /* close output BRIK file */
503
504
505 EDIT_dset_items( new_dset , /* set sub-brick factors */
506 ADN_brick_fac,ffac ,
507 ADN_none ) ;
508 free(ffac) ; /* don't need ffac no more */
509
510 /* read all the data back in if being frugal with memory */
511 if(ZCAT_frugal)
512 DSET_load(new_dset) ; /* read new dataset from disk */
513 else
514 write_output = 1 ; /* have to write everything out */
515
516 THD_load_statistics(new_dset) ; /* compute sub-brick statistics */
517
518 THD_write_3dim_dataset(NULL,NULL,new_dset,write_output); /* (re)write output file */
519 INFO_message("output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
520
521 exit(0) ; /* stage left, pursued by a bear */
522 }
523