00001
00002
00003
00004
00005
00006
00007 #include "parser.h"
00008 #include <math.h>
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <string.h>
00012 #include "mrilib.h"
00013
00014 #ifndef PI
00015 # define PI 3.14159265358979323846
00016 #endif
00017
00018 double ztone( double x ) ;
00019 double waveform( double t ) ;
00020 double waveform_GAM( double t ) ;
00021 double waveform_WAV( double t ) ;
00022 void Process_Options( int c , char * a[] ) ;
00023 void Syntax(void) ;
00024
00025 double waveform_EXPR( double t ) ;
00026
00027
00028
00029
00030
00031 #define ZT_FAC 0.50212657
00032 #define ZT_ADD 0.99576486
00033
00034 double ztone( double x )
00035 {
00036 double y , z ;
00037
00038 if( x <= 0.0 ) return 0.0 ;
00039 if( x >= 1.0 ) return 1.0 ;
00040
00041 y = (0.5*PI) * ( 1.6 * x - 0.8 ) ;
00042 z = ZT_FAC * ( tanh(tan(y)) + ZT_ADD ) ;
00043 return z ;
00044 }
00045
00046
00047
00048
00049
00050 #define WAV_TYPE 1
00051 #define GAM_TYPE 2
00052 #define EXPR_TYPE 3
00053
00054 static int waveform_type = WAV_TYPE ;
00055
00056 static double WAV_delay_time = 2.0 ,
00057 WAV_rise_time = 4.0 ,
00058 WAV_fall_time = 6.0 ,
00059 WAV_undershoot = 0.2 ,
00060 WAV_restore_time = 2.0 ;
00061
00062 static double WAV_rise_start = -666.0 ,
00063 WAV_fall_start = -666.0 ,
00064 WAV_fall_end = -666.0 ,
00065 WAV_restore_end = -666.0 ;
00066
00067 static double GAM_power = 8.6 ;
00068 static double GAM_time = 0.547 ;
00069 static double GAM_ampl = 0.0 ;
00070 static double GAM_delay_time = 0.0 ;
00071
00072 static PARSER_code * EXPR_pcode = NULL ;
00073 static double EXPR_fac = 1.0 ;
00074
00075 double waveform( double t )
00076 {
00077 switch( waveform_type ){
00078
00079 default:
00080 case WAV_TYPE: return waveform_WAV(t) ;
00081
00082 case GAM_TYPE: return waveform_GAM(t) ;
00083
00084 case EXPR_TYPE: return waveform_EXPR(t);
00085 }
00086 return 0.0 ;
00087 }
00088
00089 #define TT 19
00090 double waveform_EXPR( double t )
00091 {
00092 static int first=1 ;
00093 static double atoz[26] ;
00094
00095 if( t < 0.0 ) return 0.0 ;
00096
00097 if( first ){
00098 int ii ;
00099 for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ;
00100 first =0 ;
00101 }
00102 atoz[TT] = t ;
00103 return (EXPR_fac * PARSER_evaluate_one(EXPR_pcode,atoz) ) ;
00104 }
00105
00106 double waveform_GAM( double t )
00107 {
00108 if( GAM_ampl == 0.0 ){
00109 GAM_ampl = exp(GAM_power) / pow(GAM_power*GAM_time,GAM_power) ;
00110 }
00111
00112 if( t-GAM_delay_time <= 0.0 ) return 0.0 ;
00113
00114 return GAM_ampl * pow((t-GAM_delay_time),GAM_power) * exp(-(t-GAM_delay_time)/GAM_time) ;
00115 }
00116
00117 double waveform_WAV( double t )
00118 {
00119 if( WAV_rise_start < 0.0 ){
00120 WAV_rise_start = WAV_delay_time ;
00121 WAV_fall_start = WAV_rise_start + WAV_rise_time ;
00122 WAV_fall_end = WAV_fall_start + WAV_fall_time ;
00123 WAV_restore_end = WAV_fall_end + WAV_restore_time ;
00124 }
00125
00126 if( t < WAV_rise_start )
00127 return 0.0 ;
00128
00129 if( t < WAV_fall_start )
00130 return ztone( (t-WAV_rise_start)/WAV_rise_time ) ;
00131
00132 if( t < WAV_fall_end )
00133 return (1.0+WAV_undershoot) * ztone( (WAV_fall_end-t)/WAV_fall_time )
00134 - WAV_undershoot ;
00135
00136 if( t < WAV_restore_end )
00137 return -WAV_undershoot * ztone( (WAV_restore_end-t)/WAV_restore_time ) ;
00138
00139 return 0.0 ;
00140 }
00141
00142 static double WAV_peak = 100.0 ;
00143 static double WAV_dt = 0.1 ;
00144
00145 static int IN_npts = -666 ;
00146 static double * IN_ts = NULL ;
00147
00148 static int WAV_duration = -666.0 ;
00149 static int WAV_npts = -666 ;
00150 static double * WAV_ts = NULL ;
00151
00152 static int OUT_xy = 0 ;
00153 static int OUT_npts = -666 ;
00154 static double * OUT_ts = NULL ;
00155
00156 static int OUT_numout = -666 ;
00157
00158 static int IN_num_tstim = -666 ;
00159 static double IN_top_tstim = 0.0 ;
00160 static double * IN_tstim_a = NULL ;
00161 static double * IN_tstim_b = NULL ;
00162 static double * IN_tstim_c = NULL ;
00163
00164 int main( int argc , char * argv[] )
00165 {
00166 int ii , jj ;
00167 double val ;
00168
00169 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) Syntax() ;
00170
00171 machdep(); if(argc > 1) AFNI_logger("waver",argc,argv) ;
00172
00173 Process_Options( argc , argv ) ;
00174
00175
00176
00177 switch( waveform_type ){
00178
00179
00180
00181 default:
00182 case WAV_TYPE:
00183 WAV_duration = WAV_delay_time + WAV_rise_time + WAV_fall_time
00184 + ( (WAV_undershoot != 0.0) ? WAV_restore_time : 0.0 ) ;
00185 break ;
00186
00187
00188
00189 case GAM_TYPE:{
00190 double bal = 5.5/GAM_power + 1.0 ;
00191 double al = bal ;
00192 double lal = log(al) ;
00193
00194 while( al < bal + lal ){ al = bal + 1.1*lal ; lal = log(al) ; }
00195
00196 WAV_duration = al * GAM_power * GAM_time ;
00197 }
00198 break ;
00199
00200
00201
00202 #define FPASS 500
00203 #define SPASS 10000
00204 #define STHR 5
00205 case EXPR_TYPE:{
00206 double val , vtop=0.0 , vthr ;
00207 int itop=-1 , icount ;
00208 for( ii=0 ; ii < FPASS ; ii++ ){
00209 val = waveform_EXPR( WAV_dt * ii ) ; val = fabs(val) ;
00210 if( val > vtop ){ vtop = val ; itop = ii ; }
00211 }
00212 if( itop < 0 ){
00213 fprintf(stderr,"** -EXPR is 0 for 1st %d points!\n",FPASS);
00214 exit(1) ;
00215 }
00216 vthr = 0.01 * vtop ;
00217 for( icount=0,ii=itop+1 ; ii < SPASS && icount < STHR ; ii++ ){
00218 val = waveform_EXPR( WAV_dt * ii ) ; val = fabs(val) ;
00219 if( val <= vthr ) icount++ ;
00220 else icount=0 ;
00221 }
00222 if( ii == SPASS && icount < STHR ){
00223 fprintf(stderr,"** -EXPR doesn't decay away in %d points!\n",SPASS);
00224 exit(1) ;
00225 }
00226 WAV_duration = WAV_dt * ii ;
00227
00228 if( WAV_peak != 0.0 ) EXPR_fac = WAV_peak / vtop ;
00229 WAV_peak = 1.0 ;
00230 }
00231 break ;
00232 }
00233
00234
00235
00236 WAV_npts = (int)(1 + ceil(WAV_duration/WAV_dt)) ;
00237 WAV_ts = (double *) malloc( sizeof(double) * WAV_npts ) ;
00238
00239 for( ii=0 ; ii < WAV_npts ; ii++ )
00240 WAV_ts[ii] = WAV_peak * waveform( WAV_dt * ii ) ;
00241
00242
00243
00244 if( IN_npts < 1 && IN_num_tstim < 1 ){
00245 int top = WAV_npts ;
00246 if( OUT_numout > 0 ) top = OUT_numout ;
00247 if( OUT_xy ){
00248 for( ii=0 ; ii < top ; ii++ )
00249 printf( "%g %g\n" , WAV_dt * ii , WAV_ts[ii] ) ;
00250 } else {
00251 for( ii=0 ; ii < top ; ii++ )
00252 printf( "%g\n" , WAV_ts[ii] ) ;
00253 }
00254 exit(0) ;
00255 }
00256
00257
00258
00259 if( IN_npts > 0 ){
00260 OUT_npts = IN_npts + WAV_npts ;
00261 if( OUT_numout > 0 ) OUT_npts = OUT_numout ;
00262 OUT_ts = (double *) malloc( sizeof(double) * OUT_npts ) ;
00263 for( ii=0 ; ii < OUT_npts ; ii++ ) OUT_ts[ii] = 0.0 ;
00264
00265 for( jj=0 ; jj < IN_npts ; jj++ ){
00266 val = IN_ts[jj] ;
00267 if( val == 0.0 || fabs(val) >= 33333.0 ) continue ;
00268 for( ii=0 ; ii < WAV_npts && ii+jj < OUT_npts ; ii++ )
00269 OUT_ts[ii+jj] += val * WAV_ts[ii] ;
00270 }
00271
00272 for( jj=0 ; jj < IN_npts ; jj++ ){
00273 val = IN_ts[jj] ;
00274 if( fabs(val) >= 33333.0 ) OUT_ts[jj] = 99999.0 ;
00275 }
00276
00277 } else if( IN_num_tstim > 0 ){
00278 #undef TSTEP
00279 #define TSTEP 10
00280 int ibot,itop , kk ;
00281 int nts ;
00282 double dts = WAV_dt/TSTEP , dur , aa ;
00283 double *tst , *ast ;
00284
00285
00286
00287 OUT_npts = ceil(IN_top_tstim/WAV_dt) + WAV_npts ;
00288 if( OUT_numout > 0 ) OUT_npts = OUT_numout ;
00289 OUT_ts = (double *) malloc( sizeof(double) * OUT_npts ) ;
00290 for( ii=0 ; ii < OUT_npts ; ii++ ) OUT_ts[ii] = 0.0 ;
00291
00292
00293
00294 nts = 0 ;
00295 for( jj=0 ; jj < IN_num_tstim ; jj++ ){
00296 dur = IN_tstim_b[jj] - IN_tstim_a[jj] ; dur = MAX(dur,0.0) ;
00297 ii = ((int)ceil(dur/dts)) + 1 ;
00298 nts += ii ;
00299 }
00300
00301
00302
00303 tst = (double *) malloc( sizeof(double) * nts ) ;
00304 ast = (double *) malloc( sizeof(double) * nts ) ;
00305 nts = 0 ;
00306 for( jj=0 ; jj < IN_num_tstim ; jj++ ){
00307 dur = IN_tstim_b[jj] - IN_tstim_a[jj] ; dur = MAX(dur,0.0) ;
00308 ii = ((int)ceil(dur/dts)) + 1 ;
00309 if( ii == 1 ){
00310 tst[nts] = IN_tstim_a[jj] ; ast[nts] = IN_tstim_c[jj] ;
00311 } else {
00312 aa = dur/(WAV_dt*ii) * IN_tstim_c[jj];
00313 dur = dur / (ii-1) ;
00314 for( kk=0 ; kk < ii ; kk++ ){
00315 tst[nts+kk] = IN_tstim_a[jj]+kk*dur ;
00316 ast[nts+kk] = aa ;
00317 }
00318 }
00319 nts += ii ;
00320 }
00321
00322
00323
00324 for( jj=0 ; jj < nts ; jj++ ){
00325 ibot = (int) (tst[jj]/WAV_dt) ;
00326 itop = ibot + WAV_npts ;
00327 if( itop > OUT_npts ) itop = OUT_npts ;
00328 for( ii=ibot ; ii < itop ; ii++ ){
00329 val = WAV_peak * ast[jj] * waveform( WAV_dt * ii - tst[jj] ) ;
00330 OUT_ts[ii] += val ;
00331 }
00332 }
00333 free(ast); free(tst);
00334 }
00335
00336 if( OUT_xy ){
00337 for( ii=0 ; ii < OUT_npts ; ii++ )
00338 printf( "%g %g\n" , WAV_dt * ii , OUT_ts[ii] ) ;
00339 } else {
00340 for( ii=0 ; ii < OUT_npts ; ii++ )
00341 printf( "%g\n" , OUT_ts[ii] ) ;
00342 }
00343
00344 exit(0) ;
00345 }
00346
00347 void Syntax(void)
00348 {
00349 printf(
00350 "Usage: waver [options] > output_filename\n"
00351 "Creates an ideal waveform timeseries file.\n"
00352 "The output goes to stdout, and normally would be redirected to a file.\n"
00353 "\n"
00354 "Options: (# refers to a number; [xx] is the default value)\n"
00355 " -WAV = Sets waveform to Cox special [default]\n"
00356 " (cf. AFNI FAQ list for formulas)\n"
00357 " -GAM = Sets waveform to form t^b * exp(-t/c)\n"
00358 " (cf. Mark Cohen)\n"
00359 "\n"
00360 " -EXPR \"expression\" = Sets waveform to the expression given,\n"
00361 " which should depend on the variable 't'.\n"
00362 " e.g.: -EXPR \"step(t-2)*step(12-t)*(t-2)*(12-t)\"\n"
00363 " N.B.: The peak value of the expression on the '-dt' grid will\n"
00364 " be scaled to the value given by '-peak'; if this is not\n"
00365 " desired, set '-peak 0', and the 'natural' peak value of\n"
00366 " the expression will be used.\n"
00367 "\n"
00368 "These options set parameters for the -WAV waveform.\n"
00369 " -delaytime # = Sets delay time to # seconds [2]\n"
00370 " -risetime # = Sets rise time to # seconds [4]\n"
00371 " -falltime # = Sets fall time to # seconds [6]\n"
00372 " -undershoot # = Sets undershoot to # times the peak [0.2]\n"
00373 " (this should be a nonnegative factor)\n"
00374 " -restoretime # = Sets time to restore from undershoot [2]\n"
00375 "\n"
00376 "These options set parameters for the -GAM waveform:\n"
00377 " -gamb # = Sets the parameter 'b' to # [8.6]\n"
00378 " -gamc # = Sets the parameter 'c' to # [0.547]\n"
00379 " -gamd # = Sets the delay time to # seconds [0.0]\n"
00380 "\n"
00381 "These options apply to all waveform types:\n"
00382 " -peak # = Sets peak value to # [100]\n"
00383 " -dt # = Sets time step of output AND input [0.1]\n"
00384 " -TR # = '-TR' is equivalent to '-dt'\n"
00385 "\n"
00386 "The default is just to output the waveform defined by the parameters\n"
00387 "above. If an input file is specified by one the options below, then\n"
00388 "the timeseries defined by that file will be convolved with the ideal\n"
00389 "waveform defined above -- that is, each nonzero point in the input\n"
00390 "timeseries will generate a copy of the waveform starting at that point\n"
00391 "in time, with the amplitude scaled by the input timeseries value.\n"
00392 "\n"
00393 " -xyout = Output data in 2 columns:\n"
00394 " 1=time 2=waveform (useful for graphing)\n"
00395 " [default is 1 column=waveform]\n"
00396 "\n"
00397 " -input infile = Read timeseries from *.1D formatted 'infile';\n"
00398 " convolve with waveform to produce output\n"
00399 " N.B.: you can use a sub-vector selector to choose\n"
00400 " a particular column of infile, as in\n"
00401 " -input 'fred.1D[3]'\n"
00402 "\n"
00403 " -inline DATA = Read timeseries from command line DATA;\n"
00404 " convolve with waveform to produce output\n"
00405 " DATA is in the form of numbers and\n"
00406 " count@value, as in\n"
00407 " -inline 20@0.0 5@1.0 30@0.0 1.0 20@0.0 2.0\n"
00408 " which means a timeseries with 20 zeros, then 5 ones, then 30 zeros,\n"
00409 " a single 1, 20 more zeros, and a final 2.\n"
00410 " [The '@' character may actually be any of: '@', '*', 'x', 'X'.\n"
00411 " Note that * must be typed as \\* to prevent the shell from\n"
00412 " trying to interpret it as a filename wildcard.]\n"
00413 "\n"
00414 " -tstim DATA = Read discrete stimulation times from the command line\n"
00415 " and convolve the waveform with delta-functions at\n"
00416 " those times. In this input format, the times do\n"
00417 " NOT have to be at intervals of '-dt'. For example\n"
00418 " -dt 2.0 -tstim 5.6 9.3 13.7 16.4\n"
00419 " specifies a TR of 2 s and stimuli at 4 times\n"
00420 " (5.6 s, etc.) that do not correspond to integer\n"
00421 " multiples of TR. DATA values cannot be negative.\n"
00422 " If the DATA is stored in a file, you can read it\n"
00423 " onto the command line using something like\n"
00424 " -tstim `cat filename`\n"
00425 " where using the backward-single-quote operator\n"
00426 " of the usual Unix shells.\n"
00427 " ** 12 May 2003: The times after '-tstim' can now also be specified\n"
00428 " in the format 'a:b', indicating a continuous ON\n"
00429 " period from time 'a' to time 'b'. For example,\n"
00430 " -dt 2.0 -tstim 13.2:15.7 20.3:25.3\n"
00431 " The amplitude of a response of duration equal to\n"
00432 " 'dt' is equal the the amplitude of a single impulse\n"
00433 " response (which is the special case a=b). N.B.: This\n"
00434 " means that something like '5:5.01' is very different\n"
00435 " from '5' (='5:5'). The former will have a small amplitude\n"
00436 " because of the small duration, but the latter will have\n"
00437 " a large amplitude because the case of an instantaneous\n"
00438 " input is special. It is probably best NOT to mix the\n"
00439 " two types of input to '-tstim' for this reason.\n"
00440 " Compare the graphs from the 2 commands below:\n"
00441 " waver -dt 1.0 -tstim 5:5.1 | 1dplot -stdin\n"
00442 " waver -dt 1.0 -tstim 5 | 1dplot -stdin\n"
00443 " If you prefer, you can use the form 'a%%c' to indicate\n"
00444 " an ON interval from time=a to time=a+c.\n"
00445 " ** 13 May 2005: You can now add an amplitude to each response individually.\n"
00446 " For example\n"
00447 " waver -dt 1.0 -peak 1.0 -tstim 3.2 17.9x2.0 23.1x-0.5\n"
00448 " puts the default response amplitude at time 3.2,\n"
00449 " 2.0 times the default at time 17.9, and -0.5 times\n"
00450 " the default at time 23.1.\n"
00451 "\n"
00452 " -when DATA = Read time blocks when stimulus is 'on' (=1) from the\n"
00453 " command line and convolve the waveform with with\n"
00454 " a zero-one input. For example:\n"
00455 " -when 20..40 60..80\n"
00456 " means that the stimulus function is 1.0 for time\n"
00457 " steps number 20 to 40, and 60 to 80 (inclusive),\n"
00458 " and zero otherwise. (The first time step is\n"
00459 " numbered 0.)\n"
00460 "\n"
00461 " -numout NN = Output a timeseries with NN points; if this option\n"
00462 " is not given, then enough points are output to\n"
00463 " let the result tail back down to zero.\n"
00464 "\n"
00465 "At least one option is required, or the program will just print this message\n"
00466 "to stdout. Only one of the 3 timeseries input options above can be used.\n"
00467 "\n"
00468 "If you have the 'xmgr' graphing program, then a useful way to preview the\n"
00469 "results of this program is through a command pipe like\n"
00470 " waver -dt 0.25 -xyout -inline 16@1 40@0 16@1 40@0 | xmgr -source stdin\n"
00471 "Using the cruder AFNI package program 1dplot, you can do something like:\n"
00472 " waver -GAM -tstim 0 7.7 | 1dplot -stdin\n"
00473 "\n"
00474 "If a square wave is desired, see the 'sqwave' program.\n"
00475 ) ;
00476 exit(0) ;
00477 }
00478
00479 #define ERROR \
00480 do{fprintf(stderr,"Illegal %s option\n",argv[nopt]);Syntax();}while(0)
00481
00482 void Process_Options( int argc , char * argv[] )
00483 {
00484 int nopt = 1 ;
00485
00486 while( nopt < argc ){
00487
00488 if( strncmp(argv[nopt],"-GAM",4) == 0 ){
00489 waveform_type = GAM_TYPE ;
00490 nopt++ ; continue ;
00491 }
00492
00493 if( strncmp(argv[nopt],"-WAV",4) == 0 ){
00494 waveform_type = WAV_TYPE ;
00495 nopt++ ; continue ;
00496 }
00497
00498 if( strncmp(argv[nopt],"-EXPR",4) == 0 ){
00499 waveform_type = EXPR_TYPE ;
00500 if( EXPR_pcode != NULL ){
00501 fprintf(stderr,"** Can't have 2 -EXPR options!\n") ;
00502 exit(1) ;
00503 }
00504 nopt++ ;
00505 if( nopt >= argc ){
00506 fprintf(stderr,"** -EXPR needs an argument!\n") ; exit(1) ;
00507 }
00508 EXPR_pcode = PARSER_generate_code( argv[nopt] ) ;
00509 if( EXPR_pcode == NULL ){
00510 fprintf(stderr,"** Illegal -EXPR expression!\n") ;
00511 exit(1) ;
00512 }
00513 if( !PARSER_has_symbol("T",EXPR_pcode) ){
00514 fprintf(stderr,"** -EXPR expression doesn't use variable 't'!\n");
00515 exit(1) ;
00516 }
00517 nopt++ ; continue ;
00518 }
00519
00520 if( strncmp(argv[nopt],"-gamb",5) == 0 ){
00521 if( nopt+1 >= argc ) ERROR ;
00522 GAM_power = strtod(argv[nopt+1],NULL) ;
00523 if( GAM_power <= 0.0 ) ERROR ;
00524 waveform_type = GAM_TYPE ;
00525 nopt++ ; nopt++ ; continue ;
00526 }
00527
00528 if( strncmp(argv[nopt],"-gamc",5) == 0 ){
00529 if( nopt+1 >= argc ) ERROR ;
00530 GAM_time = strtod(argv[nopt+1],NULL) ;
00531 if( GAM_time <= 0.0 ) ERROR ;
00532 waveform_type = GAM_TYPE ;
00533 nopt++ ; nopt++ ; continue ;
00534 }
00535
00536 if( strncmp(argv[nopt],"-gamd",5) == 0 ){
00537 if( nopt+1 >= argc ) ERROR ;
00538 GAM_delay_time = strtod(argv[nopt+1],NULL) ;
00539
00540 waveform_type = GAM_TYPE ;
00541 nopt++ ; nopt++ ; continue ;
00542 }
00543
00544 if( strncmp(argv[nopt],"-del",4) == 0 ){
00545 if( nopt+1 >= argc ) ERROR ;
00546 WAV_delay_time = strtod(argv[nopt+1],NULL) ;
00547 if( WAV_delay_time < 0.0 ) ERROR ;
00548 waveform_type = WAV_TYPE ;
00549 nopt++ ; nopt++ ; continue ;
00550 }
00551
00552 if( strncmp(argv[nopt],"-ris",4) == 0 ){
00553 if( nopt+1 >= argc ) ERROR ;
00554 WAV_rise_time = strtod(argv[nopt+1],NULL) ;
00555 if( WAV_rise_time <= 0.0 ) ERROR ;
00556 waveform_type = WAV_TYPE ;
00557 nopt++ ; nopt++ ; continue ;
00558 }
00559
00560 if( strncmp(argv[nopt],"-fal",4) == 0 ){
00561 if( nopt+1 >= argc ) ERROR ;
00562 WAV_fall_time = strtod(argv[nopt+1],NULL) ;
00563 if( WAV_fall_time <= 0.0 ) ERROR ;
00564 waveform_type = WAV_TYPE ;
00565 nopt++ ; nopt++ ; continue ;
00566 }
00567
00568 if( strncmp(argv[nopt],"-und",4) == 0 ){
00569 if( nopt+1 >= argc ) ERROR ;
00570 WAV_undershoot = strtod(argv[nopt+1],NULL) ;
00571 waveform_type = WAV_TYPE ;
00572 nopt++ ; nopt++ ; continue ;
00573 }
00574
00575 if( strncmp(argv[nopt],"-pea",4) == 0 ){
00576 if( nopt+1 >= argc ) ERROR ;
00577 WAV_peak = strtod(argv[nopt+1],NULL) ;
00578 nopt++ ; nopt++ ; continue ;
00579 }
00580
00581 if( strncmp(argv[nopt],"-res",4) == 0 ){
00582 if( nopt+1 >= argc ) ERROR ;
00583 WAV_restore_time = strtod(argv[nopt+1],NULL) ;
00584 if( WAV_restore_time <= 0.0 ) ERROR ;
00585 waveform_type = WAV_TYPE ;
00586 nopt++ ; nopt++ ; continue ;
00587 }
00588
00589 if( strncmp(argv[nopt],"-dt",3) == 0 || strncmp(argv[nopt],"-TR",3) == 0 ){
00590 if( nopt+1 >= argc ) ERROR ;
00591 WAV_dt = strtod(argv[nopt+1],NULL) ;
00592 if( WAV_dt <= 0.0 ) ERROR ;
00593 nopt++ ; nopt++ ; continue ;
00594 }
00595
00596 if( strncmp(argv[nopt],"-xyo",4) == 0 ){
00597 OUT_xy = 1 ;
00598 nopt++ ; continue ;
00599 }
00600
00601 if( strncmp(argv[nopt],"-inp",4) == 0 ){
00602 MRI_IMAGE * tsim ;
00603 float * tsar ;
00604 int ii ;
00605
00606 if( IN_npts > 0 || IN_num_tstim > 0 ){
00607 fprintf(stderr,"Cannot input two timeseries!\n") ;
00608 exit(1) ;
00609 }
00610
00611 if( nopt+1 >= argc ) ERROR ;
00612 tsim = mri_read_1D( argv[nopt+1] ) ;
00613 if( tsim == NULL ) ERROR ;
00614
00615 IN_npts = tsim->nx ;
00616 IN_ts = (double *) malloc( sizeof(double) * IN_npts ) ;
00617 tsar = MRI_FLOAT_PTR(tsim) ;
00618 for( ii=0 ; ii < IN_npts ; ii++ ) IN_ts[ii] = tsar[ii] ;
00619 mri_free(tsim) ;
00620 nopt++ ; nopt++ ; continue ;
00621 }
00622
00623 if( strcmp(argv[nopt],"-tstim") == 0 ){
00624 int iopt , nnn , zero_valc=0 ;
00625 float value , valb , valc ;
00626 char *cpt , *dpt ;
00627
00628 if( IN_num_tstim > 0 || IN_npts > 0 ){
00629 fprintf(stderr,"Cannot input two timeseries!\n"); exit(1);
00630 }
00631 if( nopt+1 >= argc ) ERROR ;
00632
00633 iopt = nopt+1 ;
00634 IN_num_tstim = 0 ;
00635 IN_tstim_a = (double *) malloc( sizeof(double) ) ;
00636 IN_tstim_b = (double *) malloc( sizeof(double) ) ;
00637 IN_tstim_c = (double *) malloc( sizeof(double) ) ;
00638
00639
00640
00641 while( iopt < argc && argv[iopt][0] != '-' ){
00642
00643 if( isspace(argv[iopt][0]) ){
00644 fprintf(stderr,
00645 "** Skipping -tstim value #%d that starts with whitespace!\n",
00646 IN_num_tstim ) ;
00647 iopt++; continue;
00648 }
00649
00650
00651
00652
00653
00654
00655 valb = 0.0 ; valc = 1.0 ;
00656 if( strchr(argv[iopt],'%') != NULL ){
00657 nnn = sscanf( argv[iopt] , "%f%%%f" , &value , &valb ) ;
00658 if( nnn == 2 ) valb += value ;
00659 } else if( strchr(argv[iopt],':') != NULL ){
00660 nnn = sscanf( argv[iopt] , "%f:%f" , &value , &valb ) ;
00661 } else {
00662 nnn = sscanf( argv[iopt] , "%f" , &value ) ;
00663 }
00664 if( nnn < 1 || value < 0.0 ){
00665 fprintf(stderr,"** Weird value after -tstim: argv='%s'\n",argv[iopt] ) ;
00666 fprintf(stderr,"** previous argv='%s'\n",argv[iopt-1]) ;
00667 fprintf(stderr,"** ==> Skipping this value!\n") ;
00668 iopt++; continue;
00669 }
00670 if( nnn == 1 || valb < value ) valb = value ;
00671
00672
00673
00674 cpt = strchr(argv[iopt],'x') ;
00675 if( cpt != NULL ){
00676 cpt++ ; valc = strtod( cpt , &dpt ) ;
00677 if( valc == 0.0 && dpt == cpt ) valc = 1.0 ;
00678 if( valc == 0.0 ) zero_valc++ ;
00679 }
00680
00681 IN_tstim_a = (double *)realloc(IN_tstim_a,sizeof(double)*(IN_num_tstim+1));
00682 IN_tstim_b = (double *)realloc(IN_tstim_b,sizeof(double)*(IN_num_tstim+1));
00683 IN_tstim_c = (double *)realloc(IN_tstim_c,sizeof(double)*(IN_num_tstim+1));
00684 IN_tstim_a[IN_num_tstim] = value ;
00685 IN_tstim_b[IN_num_tstim] = valb ;
00686 IN_tstim_c[IN_num_tstim] = valc ;
00687 IN_num_tstim++ ;
00688 if( valb > IN_top_tstim ) IN_top_tstim = valb ;
00689 iopt++ ;
00690 }
00691 if( zero_valc == IN_num_tstim )
00692 fprintf(stderr,
00693 "** WARNING: all '/' amplitudes in 'waver -tstim' are zero!\n") ;
00694
00695 nopt = iopt ; continue ;
00696 }
00697
00698 if( strncmp(argv[nopt],"-inl",4) == 0 ){
00699 int iopt , count , nnn ;
00700 float value ;
00701 char sep ;
00702
00703 if( IN_npts > 0 || IN_num_tstim > 0 ){
00704 fprintf(stderr,"Cannot input two timeseries!\n") ;
00705 exit(1) ;
00706 }
00707
00708 if( nopt+1 >= argc ) ERROR ;
00709 iopt = nopt+1 ;
00710 IN_npts = 0 ;
00711 IN_ts = (double *) malloc( sizeof(double) ) ;
00712 while( iopt < argc && argv[iopt][0] != '-' ){
00713
00714 if( strstr(argv[iopt],"@") != NULL ||
00715 strstr(argv[iopt],"x") != NULL ||
00716 strstr(argv[iopt],"X") != NULL ||
00717 strstr(argv[iopt],"*") != NULL ){
00718
00719 nnn = sscanf( argv[iopt] , "%d%c%f" , &count , &sep , &value ) ;
00720 if( nnn != 3 || count < 1 ){
00721 fprintf(stderr,"Illegal value after -inline: %s\n",argv[iopt]) ;
00722 exit(1) ;
00723 }
00724
00725 } else {
00726 count = 1 ;
00727 nnn = sscanf( argv[iopt] , "%f" , &value ) ;
00728 if( nnn != 1 ){
00729 fprintf(stderr,"Illegal value after -inline: %s\n",argv[iopt]) ;
00730 exit(1) ;
00731 }
00732 }
00733
00734 IN_ts = (double *) realloc( IN_ts , sizeof(double) * (IN_npts+count) ) ;
00735 for( nnn=0 ; nnn < count ; nnn++ )
00736 IN_ts[nnn+IN_npts] = value ;
00737
00738 IN_npts += count ; iopt++ ;
00739 }
00740 nopt = iopt ; continue ;
00741 }
00742
00743 if( strcmp(argv[nopt],"-when") == 0 ){
00744 int iopt , bot,top , nn , nbt,*bt , count=0 , ii,kk ;
00745 float value ;
00746 char sep ;
00747
00748 if( IN_npts > 0 || IN_num_tstim > 0 ){
00749 fprintf(stderr,"Cannot input two timeseries!\n") ;
00750 exit(1) ;
00751 }
00752
00753 if( nopt+1 >= argc ) ERROR ;
00754 iopt = nopt+1 ;
00755 nbt = 0 ;
00756 bt = (int *) malloc(sizeof(int)) ;
00757 while( iopt < argc && argv[iopt][0] != '-' ){
00758
00759
00760
00761 bot = top = -1 ;
00762 nn = sscanf( argv[iopt],"%d..%d",&bot,&top) ;
00763 if( nn < 1 || bot < 0 ){
00764 fprintf(stderr,"Illegal value after -when: %s\n",argv[iopt]) ;
00765 exit(1) ;
00766 }
00767 if( nn == 1 ){
00768 top = bot ;
00769 } else if( top < bot ){
00770 fprintf(stderr,"Illegal value after -when: %s\n",argv[iopt]) ;
00771 exit(1) ;
00772 }
00773
00774
00775
00776 bt = (int *) realloc( bt , sizeof(int)*(nbt+1)*2 ) ;
00777 bt[2*nbt ] = bot ;
00778 bt[2*nbt+1] = top ; nbt++ ;
00779 if( count < top ) count = top ;
00780 iopt++ ;
00781 }
00782
00783 if( nbt < 1 ){
00784 fprintf(stderr,"No ranges after -when?\n") ; exit(1) ;
00785 }
00786
00787 IN_npts = count+1 ;
00788 IN_ts = (double *) malloc( sizeof(double) * IN_npts ) ;
00789 for( ii=0 ; ii < IN_npts ; ii++ ) IN_ts[ii] = 0.0 ;
00790
00791 for( kk=0 ; kk < nbt ; kk++ ){
00792 bot = bt[2*kk] ; top = bt[2*kk+1] ;
00793 for( ii=bot ; ii <= top ; ii++ ) IN_ts[ii] = 1.0 ;
00794 }
00795
00796 free(bt) ; nopt = iopt ; continue ;
00797 }
00798
00799 if( strcmp(argv[nopt],"-numout") == 0 ){
00800 int val = -1 ;
00801 if( nopt+1 >= argc ) ERROR ;
00802 sscanf(argv[nopt+1],"%d",&val) ;
00803 if( val <= 1 ){
00804 fprintf(stderr,"Illegal value after -numout: %s\n",argv[nopt]);
00805 exit(1);
00806 }
00807 OUT_numout = val ;
00808 nopt++ ; nopt++ ; continue ;
00809 }
00810
00811 ERROR ;
00812 }
00813
00814 if( WAV_peak == 0.0 && waveform_type != EXPR_TYPE ){
00815 fprintf(stderr,"** Illegal -peak 0 for non-EXPR waveform type!\n") ;
00816 exit(1) ;
00817 }
00818
00819 return ;
00820 }
