Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
1deval.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "parser.h"
00008 #include <ctype.h>
00009 #include <stdlib.h>
00010 #include "mrilib.h"
00011
00012 int main( int argc , char * argv[] )
00013 {
00014 PARSER_code * pcode = NULL ;
00015 char sym[4] ;
00016 double atoz[26] , value , del=1.0 ;
00017 int ii,jj , kvar , nopt , qvar , num=-1 , verbose=0 ;
00018 MRI_IMAGE * inim[26] ;
00019 float * inar[26] ;
00020 MRI_IMAGE *dindex_im = NULL;
00021 float *dindex = NULL;
00022 char abet[] = "abcdefghijklmnopqrstuvwxyz" ;
00023
00024
00025
00026 if( argc < 3 ){
00027 printf("Usage: 1deval [options] -expr 'expression'\n"
00028 "Evaluates the expression and writes the result to stdout.\n"
00029 "Any single letter from a-z can be used as the independent\n"
00030 "variable in the expression.\n"
00031 "\n"
00032 "Options:\n"
00033 " -del d = Use 'd' as the step for the variable in the\n"
00034 " expression [default = 1.0]\n"
00035 " -num n = Evaluate the expression 'n' times.\n"
00036 " If -num is not used, then the length of an\n"
00037 " input time series is used. If there are no\n"
00038 " time series input, then -num is required.\n"
00039 " -a q.1D = Read time series file q.1D and assign it\n"
00040 " to the symbol 'a' (as in 3dcalc).\n"
00041 " -index i.1D = Read index column from file i.1D and\n"
00042 " write it out as 1st column of output.\n"
00043 " This option is useful when working with\n"
00044 " surface data.\n"
00045 "Examples:\n"
00046 " 1deval -expr 'sin(2*PI*t)' -del 0.01 -num 101 > sin.1D\n"
00047 " 1deval -expr 'a*b*x' -a fred.1D -b ethel.1D > x.1D\n"
00048
00049 "\n"
00050 TS_HELP_STRING
00051 ) ;
00052 exit(0) ;
00053 }
00054
00055 machdep() ;
00056
00057
00058
00059 for( ii=0 ; ii < 26 ; ii++ ){
00060 atoz[ii] = 0.0 ; inim[ii] = NULL ; inar[ii] = NULL ;
00061 }
00062
00063
00064
00065 nopt = 1 ;
00066 while( nopt < argc ){
00067
00068 if( strcmp(argv[nopt],"-verb") == 0 ){
00069 verbose++ ;
00070 nopt++ ; continue ;
00071 }
00072
00073 if( strlen(argv[nopt]) == 2 && argv[nopt][0] == '-' &&
00074 argv[nopt][1] >= 'a' && argv[nopt][1] <= 'z' ){
00075
00076 int ival = argv[nopt][1] - 'a' ;
00077
00078 if( inim[ival] != NULL ){
00079 fprintf(stderr,"** Can't define symbol %c twice!\n",argv[nopt][1]);
00080 exit(1) ;
00081 }
00082
00083 nopt++ ;
00084 if( nopt >= argc ){
00085 fprintf(stderr,"** -%c needs an argument!\n",argv[nopt][1]); exit(1);
00086 }
00087
00088 inim[ival] = mri_read_1D( argv[nopt] ) ;
00089 if( inim[ival] == NULL ){
00090 fprintf(stderr,"** Can't read time series file %s\n",argv[nopt]);
00091 exit(1) ;
00092 }
00093
00094 inar[ival] = MRI_FLOAT_PTR(inim[ival]) ;
00095 nopt++ ; continue ;
00096 }
00097
00098 if( strcmp(argv[nopt],"-expr") == 0 ){
00099 if( pcode != NULL ){
00100 fprintf(stderr,"** Can't have 2 -expr options!\n") ;
00101 exit(1) ;
00102 }
00103 nopt++ ;
00104 if( nopt >= argc ){
00105 fprintf(stderr,"** -expr needs an argument!\n") ; exit(1) ;
00106 }
00107 pcode = PARSER_generate_code( argv[nopt] ) ;
00108 if( pcode == NULL ){
00109 fprintf(stderr,"** Illegal expression!\n") ;
00110 exit(1) ;
00111 }
00112 nopt++ ; continue ;
00113 }
00114
00115 if( strcmp(argv[nopt],"-del") == 0 ){
00116 nopt++ ;
00117 if( nopt >= argc ){
00118 fprintf(stderr,"** -del needs an argument!\n") ;
00119 exit(1) ;
00120 }
00121 del = strtod( argv[nopt] , NULL ) ;
00122 if( del == 0 ){
00123 fprintf(stderr,"** -del value must not be zero!\n") ;
00124 exit(1) ;
00125 }
00126 if( verbose ) fprintf(stderr,"del set to %g\n",del) ;
00127 nopt++ ; continue ;
00128 }
00129
00130 if( strcmp(argv[nopt],"-index") == 0 ){
00131 nopt++ ;
00132 if( nopt >= argc ){
00133 fprintf(stderr,"** -index needs an argument!\n") ;
00134 exit(1) ;
00135 }
00136
00137 dindex_im = mri_read_1D( argv[nopt] ) ;
00138 if( dindex_im == NULL ){
00139 fprintf(stderr,"** Can't read time series file %s\n",argv[nopt]);
00140 exit(1) ;
00141 }
00142 if (dindex_im->ny != 1) {
00143 fprintf(stderr,"** Only one column allowed for indexing.\n Found %d columns\n", dindex_im->ny);
00144 exit(1) ;
00145 }
00146 dindex = MRI_FLOAT_PTR(dindex_im) ;
00147 nopt++ ; continue ;
00148
00149 }
00150
00151
00152 if( strcmp(argv[nopt],"-num") == 0 ){
00153 nopt++ ;
00154 if( nopt >= argc ){
00155 fprintf(stderr,"** -num needs an argument!\n") ;
00156 exit(1) ;
00157 }
00158 num = strtol( argv[nopt] , NULL , 10 ) ;
00159 if( num <= 0 ){
00160 fprintf(stderr,"** -num value must be positive!\n") ;
00161 exit(1) ;
00162 }
00163 nopt++ ; continue ;
00164 }
00165
00166 fprintf(stderr,"** %s = unknown command line option!\n",argv[nopt]) ;
00167 exit(1) ;
00168 }
00169
00170 if( num <= 0 ){
00171 for( ii=0 ; ii < 26 ; ii++ ){
00172 if( inim[ii] != NULL ){ num = inim[ii]->nx ; break ; }
00173 }
00174 if( num > 0 ){
00175 if( verbose )
00176 fprintf(stderr,"++ Set num = %d from input time series\n",num);
00177 } else {
00178 fprintf(stderr,"** Need to supply -num on command line!\n"); exit(1);
00179 }
00180 }
00181
00182 for( qvar=ii=0 ; ii < 26 ; ii++ ){
00183 if( inim[ii] != NULL && inim[ii]->nx < num ){
00184 fprintf(stderr,"** Time series file %s is too short!\n",inim[ii]->name);
00185 qvar++ ;
00186 }
00187 }
00188
00189 if (dindex) {
00190 if ( dindex_im->nx != num) {
00191 fprintf(stderr,"** Number of values in index column (%d) \n not equal to number of values in data (%d)\n", dindex_im->nx, num );
00192 exit(1) ;
00193 }
00194 }
00195
00196 if( qvar > 0 ) exit(1) ;
00197
00198 if( pcode == NULL ){
00199 fprintf(stderr,"** -expr is missing!\n") ; exit(1) ;
00200 }
00201
00202 qvar = 0 ; kvar = -1 ;
00203 for( ii=0 ; ii < 26 ; ii++ ){
00204 sym[0] = 'A' + ii ; sym[1] = '\0' ;
00205 if( PARSER_has_symbol(sym,pcode) ){
00206 if( inim[ii] == NULL ){
00207 qvar++ ; if( kvar < 0 ) kvar = ii ;
00208 }
00209 } else if( inim[ii] != NULL ){
00210 fprintf(stderr,"++ Symbol %c defined but not used!\n",abet[ii]) ;
00211 }
00212 }
00213 if( qvar > 1 ){
00214 fprintf(stderr,"++ Found %d indeterminate symbols in expression,\n"
00215 "++ but there should only be 0 or 1.\n"
00216 "++ Will use symbol %c as the variable.\n" ,
00217 qvar , abet[kvar] ) ;
00218 }
00219
00220
00221
00222 for( ii=0 ; ii < num ; ii++ ){
00223 for( jj=0 ; jj < 26 ; jj++ )
00224 if( inar[jj] != NULL ) atoz[jj] = inar[jj][ii] ;
00225 if( kvar >= 0 ) atoz[kvar] = ii * del ;
00226 value = PARSER_evaluate_one( pcode , atoz ) ;
00227 if (dindex) printf(" %d\t%g\n", (int)dindex[ii], value) ;
00228 else printf(" %g\n",value) ;
00229 }
00230 exit(0) ;
00231 }