00001 #include "mrilib.h"
00002 
00003 /*---------------------------------------------------------------------------
00004   This program catenates multiple 3D datasets in the slice direction.
00005   -- RWCox - 08 Aug 2001
00006 ----------------------------------------------------------------------------*/
00007 
00008 #ifndef myXtFree
00009 #   define myXtFree(xp) (XtFree((char *)(xp)) , (xp)=NULL)
00010 #endif
00011 
00012 /*-------------------------- global data --------------------------*/
00013 
00014 static THD_3dim_dataset_array * ZCAT_dsar  = NULL ;  /* input datasets */
00015 static int                      ZCAT_nxy   = -1 ;    /* # voxels/slice */
00016 static int                      ZCAT_nbrik = -1 ;    /* # bricks */
00017 static int                      ZCAT_verb  = 0 ;     /* verbose? */
00018 static int                      ZCAT_type  = -1 ;    /* dataset type */
00019 static int                      ZCAT_datum = -1 ;    /* dataset datum */
00020 static int                      ZCAT_fscale = 0 ;
00021 static int                      ZCAT_gscale = 0 ;
00022 static int                      ZCAT_nscale = 0 ;
00023 
00024 #define DSUB(id) DSET_IN_3DARR(ZCAT_dsar,(id))
00025 
00026 static char ZCAT_output_prefix[THD_MAX_PREFIX] = "zcat" ;
00027 
00028 /*--------------------------- prototypes ---------------------------*/
00029 
00030 void ZCAT_read_opts( int , char ** ) ;
00031 void ZCAT_Syntax(void) ;
00032 
00033 /*--------------------------------------------------------------------
00034    read the arguments, load the global variables
00035 ----------------------------------------------------------------------*/
00036 
00037 void ZCAT_read_opts( int argc , char * argv[] )
00038 {
00039    int nopt = 1 , ii ;
00040    THD_3dim_dataset * dset ;
00041 
00042    INIT_3DARR(ZCAT_dsar) ;  /* array of datasets */
00043 
00044    while( nopt < argc ){
00045 
00046       /**** -nscale ****/
00047 
00048       if( strncmp(argv[nopt],"-nscale",6) == 0 ){
00049          ZCAT_gscale = ZCAT_fscale = 0 ; ZCAT_nscale = 1 ;
00050          nopt++ ; continue ;
00051       }
00052 
00053       /**** -fscale ****/
00054 
00055       if( strncmp(argv[nopt],"-fscale",6) == 0 ){
00056          ZCAT_fscale = 1 ; ZCAT_nscale = 0 ;
00057          nopt++ ; continue ;
00058       }
00059 
00060 #if 0
00061       /**** -gscale ****/
00062 
00063       if( strncmp(argv[nopt],"-gscale",6) == 0 ){
00064          ZCAT_gscale = ZCAT_fscale = 1 ; ZCAT_nscale = 0 ;
00065          nopt++ ; continue ;
00066       }
00067 #endif
00068 
00069       /**** -prefix prefix ****/
00070 
00071       if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00072          nopt++ ;
00073          if( nopt >= argc ){
00074             fprintf(stderr,"*** need argument after -prefix!\n") ; exit(1) ;
00075          }
00076          MCW_strncpy( ZCAT_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00077          continue ;
00078       }
00079 
00080       /**** -verb ****/
00081 
00082       if( strncmp(argv[nopt],"-verb",5) == 0 ){
00083          ZCAT_verb++ ;
00084          nopt++ ; continue ;
00085       }
00086 
00087       /**** -datum type ****/
00088 
00089       if( strncmp(argv[nopt],"-datum",6) == 0 ){
00090          if( ++nopt >= argc ){
00091             fprintf(stderr,"*** need an argument after -datum!\n"); exit(1);
00092          }
00093               if( strcmp(argv[nopt],"short") == 0 ) ZCAT_datum = MRI_short ;
00094          else if( strcmp(argv[nopt],"float") == 0 ) ZCAT_datum = MRI_float ;
00095          else if( strcmp(argv[nopt],"byte" ) == 0 ) ZCAT_datum = MRI_byte  ;
00096 #if 0
00097          else if( strcmp(argv[nopt],"complex")== 0) ZCAT_datum = MRI_complex ;
00098 #endif
00099          else {
00100             fprintf(stderr,"*** -datum %s not supported in 3dZcat!\n",
00101                     argv[nopt] ) ;
00102             exit(1) ;
00103          }
00104          nopt++ ; continue ;  /* go to next arg */
00105       }
00106 
00107       /**** Garbage ****/
00108 
00109       if( argv[nopt][0] == '-' ){
00110          fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00111       }
00112 
00113       /**** read dataset ****/
00114 
00115       if( ZCAT_verb )
00116          fprintf(stderr,"+++ Opening dataset %s\n",argv[nopt]) ;
00117 
00118       dset = THD_open_dataset( argv[nopt++] ) ;
00119       if( dset == NULL ){
00120          fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt-1]) ; exit(1) ;
00121       }
00122       THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00123 
00124       if( ZCAT_type < 0 ) ZCAT_type = dset->type ;
00125 
00126       /* check if voxel counts match, etc. */
00127 
00128       ii = dset->daxes->nxx * dset->daxes->nyy ;
00129       if( ZCAT_nxy < 0 ){
00130          ZCAT_nxy   = ii ;
00131          ZCAT_nbrik = DSET_NVALS(dset) ;
00132       } else if( ii != ZCAT_nxy ){
00133          fprintf(stderr,"*** Dataset %s differs in slice size from others\n",argv[nopt-1]);
00134          exit(1) ;
00135       } else if ( DSET_NVALS(dset) != ZCAT_nbrik ){
00136          fprintf(stderr,"*** Dataset %s has different number of sub-bricks\n",argv[nopt-1]) ;
00137          exit(1) ;
00138       } else if ( DSET_BRICK_TYPE(dset,0) == MRI_complex ){
00139          fprintf(stderr,"*** Dataset %s is complex-valued -- ILLEGAL\n",argv[nopt-1]) ;
00140          exit(1) ;
00141       }
00142 
00143       ADDTO_3DARR(ZCAT_dsar,dset) ;  /* list of datasets */
00144 
00145    }  /* end of loop over command line arguments */
00146 
00147    return ;
00148 }
00149 
00150 /*-------------------------------------------------------------------------*/
00151 
00152 void ZCAT_Syntax(void)
00153 {
00154    printf(
00155     "Usage: 3dZcat [options] dataset dataset ...\n"
00156     "Concatenates datasets in the slice (z) direction.  Each input\n"
00157     "dataset must have the same number of voxels in each slice, and\n"
00158     "must have the same number of sub-bricks.\n"
00159     "\n"
00160     "Options:\n"
00161     "  -prefix pname = Use 'pname' for the output dataset prefix name.\n"
00162     "                    [default='zcat']\n"
00163     "  -datum type   = Coerce the output data to be stored as the given\n"
00164     "                    type, which may be byte, short, or float.\n"
00165     "  -fscale     = Force scaling of the output to the maximum integer\n"
00166     "                  range.  This only has effect if the output datum\n"
00167     "                  is byte or short (either forced or defaulted).\n"
00168     "                  This option is sometimes necessary to eliminate\n"
00169     "                  unpleasant truncation artifacts.\n"
00170     "  -nscale     = Don't do any scaling on output to byte or short datasets.\n"
00171     "                   This may be especially useful when operating on mask\n"
00172     "                   datasets whose output values are only 0's and 1's.\n"
00173     "  -verb         = Print out some verbositiness as the program\n"
00174     "                    proceeds.\n"
00175     "\n"
00176     "Command line arguments after the above are taken as input datasets.\n"
00177     "A dataset is specified using one of these forms:\n"
00178     "   'prefix+view', 'prefix+view.HEAD', or 'prefix+view.BRIK'.\n"
00179     "\n"
00180 
00181     MASTER_SHORTHELP_STRING
00182 
00183     "\n"
00184     "Notes:\n"
00185     "* You can use the '3dinfo' program to see how many slices a\n"
00186     "    dataset comprises.\n"
00187     "* There must be at least two datasets input (otherwise, the\n"
00188     "    program doesn't make much sense, does it?).\n"
00189     "* Each input dataset must have the same number of voxels in each\n"
00190     "    slice, and must have the same number of sub-bricks.\n"
00191     "* This program does not deal with complex-valued datasets.\n"
00192     "* See the output of '3dZcutup -help' for a C shell script that\n"
00193     "    can be used to take a dataset apart into single slice datasets,\n"
00194     "    analyze them separately, and then assemble the results into\n"
00195     "    new 3D datasets.\n"
00196    ) ;
00197 
00198    exit(0) ;
00199 }
00200 
00201 /*-------------------------------------------------------------------------*/
00202 
00203 int main( int argc , char * argv[] )
00204 {
00205    int ninp , ids , iv , iz,kz,new_nz, nx,ny,nz,nxy,nxyz ;
00206    THD_3dim_dataset * new_dset=NULL , * dset ;
00207    THD_ivec3 iv_nxyz ;
00208    float * fvol , *ffac ;
00209    void  * svol ;
00210    int cmode , fscale ; FILE * far ;
00211 
00212    /*** read input options ***/
00213 
00214    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) ZCAT_Syntax() ;
00215 
00216    /*-- addto the arglist, if user wants to --*/
00217 
00218    { int new_argc ; char ** new_argv ;
00219      addto_args( argc , argv , &new_argc , &new_argv ) ;
00220      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00221    }
00222 
00223    mainENTRY("3dZcat main") ; machdep() ; AFNI_logger("3dZcat",argc,argv) ;
00224    PRINT_VERSION("3dZcat") ;
00225 
00226    ZCAT_read_opts( argc , argv ) ;
00227 
00228    /*** create new dataset (empty) ***/
00229 
00230    ninp = ZCAT_dsar->num ;
00231    if( ninp < 2 ){
00232       fprintf(stderr,"*** Must have at least 2 input datasets!\n") ; exit(1) ;
00233    }
00234 
00235    new_nz = 0 ;
00236    for( ids=0 ; ids < ninp ; ids++ ) new_nz += DSET_NZ(DSUB(ids)) ;
00237 
00238    if( ZCAT_verb )
00239       fprintf(stderr,"+++ Output will have %d slices\n",new_nz) ;
00240 
00241    /** find 1st dataset that is time dependent, if any **/
00242 
00243    for( ids=0 ; ids < ninp ; ids++ ){
00244       dset = DSUB(ids) ;
00245       if( DSET_TIMESTEP(dset) > 0.0 ) break ;
00246    }
00247    if( ids == ninp ) dset = DSUB(0) ; /* fallback position */
00248    if( ZCAT_verb )
00249       fprintf(stderr,"+++ Using %s as 'master' dataset\n",DSET_HEADNAME(dset)) ;
00250 
00251    new_dset = EDIT_empty_copy( dset ) ; /* make a copy of its header */
00252 
00253    tross_Copy_History( dset , new_dset ) ;
00254    tross_Make_History( "3dZcat" , argc,argv , new_dset ) ;
00255 
00256    /* modify its header */
00257 
00258    nx = DSET_NX(dset) ;
00259    ny = DSET_NY(dset) ; nxy = nx*ny ; nxyz = nxy*new_nz ;
00260 
00261    LOAD_IVEC3( iv_nxyz , nx , ny , new_nz ) ;
00262 
00263    if( ZCAT_datum < 0 ) ZCAT_datum = DSET_BRICK_TYPE(dset,0) ;
00264 
00265    EDIT_dset_items( new_dset ,
00266                       ADN_prefix    , ZCAT_output_prefix ,
00267                       ADN_type      , ZCAT_type ,
00268                       ADN_nxyz      , iv_nxyz ,
00269                       ADN_datum_all , ZCAT_datum ,
00270                       ADN_nsl       , 0 , /* kill time offsets  */
00271                     ADN_none ) ;
00272 
00273    /* can't re-write existing dataset */
00274 
00275    if( THD_is_file(DSET_HEADNAME(new_dset)) ){
00276      fprintf(stderr,"*** Fatal error: dataset %s already exists!\n",
00277              DSET_HEADNAME(new_dset) ) ;
00278      exit(1) ;
00279    }
00280 
00281    THD_force_malloc_type( new_dset->dblk , DATABLOCK_MEM_MALLOC ) ;
00282 
00283    /*-- open output BRIK file --*/
00284 
00285    cmode = THD_get_write_compression() ;
00286    far = COMPRESS_fopen_write( DSET_BRICKNAME(new_dset) , cmode ) ;
00287    if( far == NULL ){
00288       fprintf(stderr,
00289         "\a\n*** cannot open output file %s\n",DSET_BRICKNAME(new_dset)) ;
00290       exit(1) ;
00291    }
00292 
00293    /* make space for sub-brick scaling factors */
00294 
00295    ffac = (float *) malloc(sizeof(float)*DSET_NVALS(new_dset)) ;
00296 
00297    /*** Loop over output sub-bricks ***/
00298 
00299    for( iv=0 ; iv < DSET_NVALS(new_dset) ; iv++ ){
00300 
00301       if( ZCAT_verb )
00302          fprintf(stderr,"+++ Computing output sub-brick #%d\n",iv) ;
00303 
00304       fscale = ZCAT_fscale ;  /* local for this brick */
00305 
00306       /*  make temporary holding space */
00307 
00308       fvol = (float *) malloc(sizeof(float)*nxyz) ;
00309 
00310       /* loop over input datasets */
00311       /* kz = slice index in output that this input dataset starts at */
00312 
00313       for( kz=ids=0 ; ids < ninp ; ids++ ){
00314 
00315          dset = DSUB(ids) ; nz = DSET_NZ(dset) ; DSET_load(dset) ;
00316          if( ! DSET_LOADED(dset) ){
00317             fprintf(stderr,"*** Fatal error: can't load data from %s\n",
00318                     DSET_HEADNAME(dset)) ;
00319             exit(1) ;
00320          }
00321 
00322          if( ZCAT_verb == 2 )
00323             fprintf(stderr," ++ processing input %s\n",DSET_HEADNAME(dset)) ;
00324 
00325          /* copy data from input to holding space, converting to float */
00326 
00327          switch( DSET_BRICK_TYPE(dset,iv) ){
00328 
00329             default:
00330                fprintf(stderr,
00331                        "*** Illegal input brick type=%s in dataset %s\n",
00332                        MRI_TYPE_name[DSET_BRICK_TYPE(dset,iv)] ,
00333                        DSET_HEADNAME(dset) ) ;
00334             exit(1) ;
00335 
00336             case MRI_short:{
00337                register int ii , itop = DSET_NVOX(dset) ;
00338                register short *dar = DSET_ARRAY(dset,iv) ;
00339                register float fac  = DSET_BRICK_FACTOR(dset,iv) ;
00340                register float *var = fvol + kz*nxy ;
00341                if( fac == 0.0 ) fac = 1.0 ;
00342                for( ii=0 ; ii < itop ; ii++ ) var[ii] = fac*dar[ii] ;
00343                if( fac != 1.0 ) fscale = 1 ;
00344             }
00345             break ;
00346 
00347             case MRI_byte:{
00348                register int ii , itop = DSET_NVOX(dset) ;
00349                register byte *dar  = DSET_ARRAY(dset,iv) ;
00350                register float fac  = DSET_BRICK_FACTOR(dset,iv) ;
00351                register float *var = fvol + kz*nxy ;
00352                if( fac == 0.0 ) fac = 1.0 ;
00353                for( ii=0 ; ii < itop ; ii++ ) var[ii] = fac*dar[ii] ;
00354                if( fac != 1.0 ) fscale = 1 ;
00355             }
00356             break ;
00357 
00358             case MRI_float:{
00359                register int ii , itop = DSET_NVOX(dset) ;
00360                register float *dar = DSET_ARRAY(dset,iv) ;
00361                register float fac  = DSET_BRICK_FACTOR(dset,iv) ;
00362                register float *var = fvol + kz*nxy ;
00363                if( fac == 0.0 ) fac = 1.0 ;
00364                for( ii=0 ; ii < itop ; ii++ ) var[ii] = fac*dar[ii] ;
00365             }
00366             break ;
00367 
00368          } /* end of switch over input brick type */
00369 
00370          kz += nz ; DSET_unload(dset) ;
00371 
00372       } /* end of loop over input datasets */
00373 
00374       /* at this point, fvol holds the correctly
00375          scaled values for the output dataset;
00376          now, must store it in the correct datum into svol;
00377          this code is lifted/adapted from 3dcalc.c          */
00378 
00379       switch( ZCAT_datum ){
00380 
00381          default:                   /* should never transpire */
00382             fprintf(stderr,
00383                     "*** Fatal Error ***\n"
00384                     "*** Somehow ended up with -datum = %d\n",ZCAT_datum) ;
00385          exit(1) ;
00386 
00387          case MRI_float:            /* the trivial case */
00388             svol = (void *) fvol ;
00389             ffac[iv] = 0.0 ;        /* don't need no stinking factor */
00390          break ;
00391 
00392          case MRI_byte:             /* harder cases: */
00393          case MRI_short:{           /* must create svol and scale to it */
00394             float gtop , fimfac ;
00395 
00396             /* find largest value in fvol */
00397 
00398             gtop = MCW_vol_amax( nxyz , 1 , 1 , MRI_float, fvol ) ;
00399             if( gtop == 0.0 )
00400               fprintf(stderr,"+++ Warning: output sub-brick %d is all zeros!\n",iv) ;
00401 
00402             /* compute scaling factor */
00403 
00404             if( fscale ){
00405                fimfac = (gtop > 0.0) ? MRI_TYPE_maxval[ZCAT_datum] / gtop : 0.0 ;
00406             } else if( !ZCAT_nscale ){
00407                fimfac = (gtop > MRI_TYPE_maxval[ZCAT_datum] || (gtop > 0.0 && gtop <= 1.0) )
00408                         ? MRI_TYPE_maxval[ZCAT_datum]/ gtop : 0.0 ;
00409             } else {
00410                fimfac = 0.0 ;
00411             }
00412 
00413             if( ZCAT_verb == 2 ){
00414                if( fimfac != 0.0 )
00415                   fprintf(stderr," ++ Output sub-brick %d scale factor = %f\n",iv,fimfac) ;
00416                else
00417                   fprintf(stderr,"+++ Output sub-brick %d: no scale factor\n" ,iv) ;
00418             }
00419 
00420             /* now scale fvol into svol */
00421 
00422             svol = (void *) malloc( mri_datum_size(ZCAT_datum) * nxyz ) ;
00423 
00424             EDIT_coerce_scale_type( nxyz , fimfac ,
00425                                     MRI_float, fvol , ZCAT_datum,svol ) ;
00426 
00427             free(fvol) ;
00428             ffac[iv] = (fimfac > 0.0 && fimfac != 1.0) ? 1.0/fimfac : 0.0 ;
00429          }
00430          break ;
00431       } /* end of switch on output datum type */
00432 
00433       /*-- now save svol to disk --*/
00434 
00435       fwrite( svol , mri_datum_size(ZCAT_datum) , nxyz , far ) ;
00436       free(svol) ;
00437 
00438    } /* end of loop over output sub-bricks */
00439 
00440    /*-- cleanup, and write dataset header --*/
00441 
00442    if( ZCAT_verb )
00443       fprintf(stderr,"+++ Computing output sub-brick statistics\n") ;
00444 
00445    for( ids=0 ; ids < ninp ; ids++ )       /* remove inputs from memory */
00446       DSET_delete( DSUB(ids) ) ;
00447 
00448    COMPRESS_fclose(far) ;                  /* close output BRIK file */
00449    EDIT_dset_items( new_dset ,             /* set sub-brick factors */
00450                       ADN_brick_fac,ffac ,
00451                     ADN_none ) ;
00452    free(ffac) ;                            /* don't need ffac no more */
00453    DSET_load(new_dset) ;                   /* read new dataset from disk */
00454    THD_load_statistics(new_dset) ;         /* compute sub-brick statistics */
00455    DSET_write_header(new_dset) ;           /* write output HEAD */
00456    fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00457 
00458    exit(0) ;                               /* stage left */
00459 }