00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "parser.h"
00008 #include "mrilib.h"
00009 
00010 /*-------------------------- global data --------------------------*/
00011 
00012 static int                CALC_datum = -1 ;
00013 static int                CALC_nvox  = -1 ;
00014 static PARSER_code *      CALC_code  = NULL ;
00015 
00016 static int CALC_nx=-1 , CALC_ny=-1 ;  /* 27 Jul 2000 */
00017 
00018 static MRI_IMAGE * CALC_im[26] ;
00019 static int         CALC_type[26] ;
00020 static byte *      CALC_byte[26] ;
00021 static short *     CALC_short[26] ;
00022 static float *     CALC_float[26] ;
00023 
00024 static char CALC_output_name[256] = "imcalc.out" ;
00025 
00026 /*--------------------------- prototypes ---------------------------*/
00027 void CALC_read_opts( int , char ** ) ;
00028 void CALC_Syntax(void) ;
00029 
00030 /*--------------------------------------------------------------------
00031    read the arguments, load the global variables
00032 ----------------------------------------------------------------------*/
00033 
00034 void CALC_read_opts( int argc , char * argv[] )
00035 {
00036    int nopt = 1 ;
00037    int ids ;
00038 
00039    for( ids=0 ; ids < 26 ; ids++ ){
00040       CALC_type[ids] = -1 ;
00041    }
00042 
00043    while( nopt < argc && argv[nopt][0] == '-' ){
00044 
00045       /**** -nxy nx ny ****/
00046 
00047       if( strncmp(argv[nopt],"-nxy",4) == 0 ){
00048          if( ++nopt >= argc-1 ){
00049             fprintf(stderr,"need 2 arguments after -nxy!\n"); exit(1);
00050          }
00051          CALC_nx = strtol(argv[nopt++],NULL,10) ;
00052          CALC_ny = strtol(argv[nopt++],NULL,10) ;
00053          continue ;  /* go to next arg */
00054       }
00055 
00056       /**** -datum type ****/
00057 
00058       if( strncmp(argv[nopt],"-datum",6) == 0 ){
00059          if( ++nopt >= argc ){
00060             fprintf(stderr,"need an argument after -datum!\n") ; exit(1) ;
00061          }
00062          if( strcmp(argv[nopt],"short") == 0 ){
00063             CALC_datum = MRI_short ;
00064          } else if( strcmp(argv[nopt],"float") == 0 ){
00065             CALC_datum = MRI_float ;
00066          } else if( strcmp(argv[nopt],"byte") == 0 ){
00067             CALC_datum = MRI_byte ;
00068          } else {
00069             fprintf(stderr,"-datum of type '%s' is not supported in 3dmerge!\n",
00070                     argv[nopt] ) ;
00071             exit(1) ;
00072          }
00073          nopt++ ; continue ;  /* go to next arg */
00074       }
00075 
00076       /**** -output name ****/
00077 
00078       if( strncmp(argv[nopt],"-output",6) == 0 ){
00079          nopt++ ;
00080          if( nopt >= argc ){
00081             fprintf(stderr,"need argument after -output!\n") ; exit(1) ;
00082          }
00083          strncpy( CALC_output_name , argv[nopt++] , 255 ) ; CALC_output_name[255] = '\0' ;
00084          continue ;
00085       }
00086 
00087       /**** -expr expression ****/
00088 
00089       if( strncmp(argv[nopt],"-expr",4) == 0 ){
00090          if( CALC_code != NULL ){
00091             fprintf(stderr,"cannot have 2 -expr options!\n") ; exit(1) ;
00092          }
00093          nopt++ ;
00094          if( nopt >= argc ){
00095             fprintf(stderr,"need argument after -expr!\n") ; exit(1) ;
00096          }
00097          CALC_code = PARSER_generate_code( argv[nopt++] ) ;
00098          if( CALC_code == NULL ){
00099             fprintf(stderr,"illegal expression!\n") ; exit(1) ;
00100          }
00101          continue ;
00102       }
00103 
00104       /**** -[a-z] filename ****/
00105 
00106       ids = strlen(argv[nopt]) ;
00107 
00108       if( (argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z') && (ids == 2) ){
00109 
00110          int ival , nxyz ;
00111          MRI_IMAGE * im ;
00112 
00113          ival = argv[nopt][1] - 'a' ;
00114          if( CALC_im[ival] != NULL ){
00115             fprintf(stderr,"can't open duplicate %s images!\n",argv[nopt]) ;
00116             exit(1) ;
00117          }
00118 
00119          nopt++ ;
00120          if( nopt >= argc ){
00121             fprintf(stderr,"need argument after %s!\n",argv[nopt-1]) ; exit(1) ;
00122          }
00123 
00124          im = mri_read_just_one( argv[nopt++] ) ;
00125          if( im == NULL ){
00126             fprintf(stderr,"can't open image %s\n",argv[nopt-1]) ; exit(1) ;
00127          }
00128 
00129          nxyz = im->nvox ;
00130          if( CALC_nvox < 0 ){
00131             CALC_nvox = nxyz ;
00132          } else if( nxyz != CALC_nvox ){
00133             fprintf(stderr,"image %s differs in size from others\n",argv[nopt-1]);
00134             exit(1) ;
00135          }
00136 
00137          CALC_type[ival] = im->kind ;
00138          CALC_im[ival]   = im ;
00139 
00140          switch( CALC_type[ival] ){
00141             default: fprintf(stderr,"illegal datum in image %s\n",argv[nopt-1]);
00142                      exit(1) ;
00143 
00144             case MRI_short: CALC_short[ival] = MRI_SHORT_PTR(im) ; break ;
00145             case MRI_float: CALC_float[ival] = MRI_FLOAT_PTR(im) ; break ;
00146             case MRI_byte : CALC_byte [ival] = MRI_BYTE_PTR(im)  ; break ;
00147          }
00148 
00149          if( CALC_datum < 0 ) CALC_datum = CALC_type[ival] ;
00150          continue ;
00151       }
00152 
00153       fprintf(stderr,"Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00154 
00155    }  /* end of loop over options */
00156 
00157    /*** cleanup ***/
00158 
00159    if( nopt < argc ){
00160      fprintf(stderr,"Extra command line arguments puzzle me! %s ...\n",argv[nopt]) ;
00161      exit(1) ;
00162    }
00163 
00164    for( ids=0 ; ids < 26 ; ids++ ) if( CALC_im[ids] != NULL ) break ;
00165 
00166    if( ids == 26 && (CALC_nx < 2 || CALC_ny < 2) ){
00167       fprintf(stderr,"No input images given!\n") ; exit(1) ;
00168    }
00169 
00170    if( CALC_code == NULL ){
00171       fprintf(stderr,"No expression given!\n") ; exit(1) ;
00172    }
00173 
00174    return ;
00175 }
00176 
00177 /*------------------------------------------------------------------*/
00178 
00179 void CALC_Syntax(void)
00180 {
00181    printf(
00182     "Do arithmetic on 2D images, pixel-by-pixel.\n"
00183     "Usage: imcalc options\n"
00184     "where the options are:\n"
00185    ) ;
00186 
00187    printf(
00188     "  -datum type = Coerce the output data to be stored as the given type,\n"
00189     "                  which may be byte, short, or float.\n"
00190     "                  [default = datum of first input image]\n"
00191     "  -a dname    = Read image 'dname' and call the voxel values 'a'\n"
00192     "                  in the expression.  'a' may be any letter from 'a' to 'z'.\n"
00193     "               ** If some letter name is used in the expression, but not\n"
00194     "                  present in one of the image options here, then that\n"
00195     "                  variable is set to 0.\n"
00196     "  -expr \"expression\"\n"
00197     "                Apply the expression within quotes to the input images,\n"
00198     "                  one voxel at a time, to produce the output image.\n"
00199     "                  (\"sqrt(a*b)\" to compute the geometric mean, for example)\n"
00200     "  -output name = Use 'name' for the output image filename.\n"
00201     "                  [default='imcalc.out']\n"
00202     "\n"
00203     "See the output of '3dcalc -help' for details on what kinds of expressions\n"
00204     "are possible.  Note that complex-valued images cannot be processed (byte,\n"
00205     "short, and float are OK).\n"
00206    ) ;
00207    exit(0) ;
00208 }
00209 
00210 /*------------------------------------------------------------------*/
00211 
00212 int main( int argc , char * argv[] )
00213 {
00214    double atoz[26] ;
00215    int ii , ids ;
00216    MRI_IMAGE * new_im ;
00217    float * fnew ;
00218 
00219    /*** read input options ***/
00220 
00221    if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) CALC_Syntax() ;
00222 
00223    machdep() ;
00224 
00225    CALC_read_opts( argc , argv ) ;
00226 
00227    /*** make output image (always float datum) ***/
00228 
00229    if( mri_filesize(CALC_output_name) > 0 ){
00230       fprintf(stderr,
00231               "*** Output file %s already exists -- cannot continue!\n",
00232               CALC_output_name ) ;
00233       exit(1) ;
00234    }
00235 
00236    for( ids=0 ; ids < 26 ; ids++ ) if( CALC_im[ids] != NULL ) break ;
00237    if( ids < 26 ){
00238       new_im = mri_new_conforming( CALC_im[ids] , MRI_float ) ;
00239    } else if( CALC_nx > 1 && CALC_ny > 1 ){
00240       new_im = mri_new( CALC_nx , CALC_ny , MRI_float ) ;
00241       CALC_nvox = new_im->nvox ;
00242    }
00243    fnew = MRI_FLOAT_PTR(new_im) ;
00244 
00245    for( ids=0 ; ids < 26 ; ids++ ) atoz[ids] = 0.0 ;
00246 
00247    /*** loop over voxels and compute result ***/
00248 
00249    for( ii=0 ; ii < CALC_nvox ; ii++ ){
00250 
00251       for( ids=0 ; ids < 26 ; ids++ ){
00252          switch( CALC_type[ids] ){
00253             case MRI_short: atoz[ids] = CALC_short[ids][ii] ; break;
00254             case MRI_float: atoz[ids] = CALC_float[ids][ii] ; break;
00255             case MRI_byte : atoz[ids] = CALC_byte [ids][ii] ; break;
00256          }
00257       }
00258 
00259       fnew[ii] = PARSER_evaluate_one( CALC_code , atoz ) ;
00260    }
00261 
00262    /*** toss input images ***/
00263 
00264    for( ids=0 ; ids < 26 ; ids++ )
00265       if( CALC_im[ids] != NULL ) mri_free( CALC_im[ids] ) ;
00266 
00267    /*** scale to output image, if needed ***/
00268 
00269    switch( CALC_datum ){
00270 
00271       case MRI_short:{
00272          float top ;
00273          MRI_IMAGE * shim ;
00274 
00275          top = mri_maxabs( new_im ) ;
00276          if( top < 32767.0 ) shim = mri_to_short( 1.0 , new_im ) ;
00277          else                shim = mri_to_short_scl( 0.0 , 10000.0 , new_im ) ;
00278 
00279          mri_free(new_im) ; new_im = shim ;
00280       }
00281       break ;
00282 
00283       case MRI_byte:{
00284          float top ;
00285          MRI_IMAGE * bim ;
00286 
00287          top = mri_maxabs( new_im ) ;
00288          if( top < 255.0 ) bim = mri_to_byte_scl( 1.0 ,   0.0 , new_im ) ;
00289          else              bim = mri_to_byte_scl( 0.0 , 255.0 , new_im ) ;
00290 
00291          mri_free(new_im) ; new_im = bim ;
00292       }
00293       break ;
00294    }
00295 
00296    /*** done ***/
00297 
00298    mri_write( CALC_output_name , new_im ) ;
00299    exit(0) ;
00300 }