Do arithmetic on 3D datasets, voxel-by-voxel [no inter-voxel computation]. Usage: 3dcalc [options] where the options are: -verbose = Makes the program print out various information as it progresses. -datum type = Coerce the output data to be stored as the given type, which may be byte, short, or float. [default = datum of first input dataset] -fscale = Force scaling of the output to the maximum integer range. This only has effect if the output datum is byte or short (either forced or defaulted). This option is often necessary to eliminate unpleasant truncation artifacts. [The default is to scale only if the computed values seem to need it -- are all <= 1.0 or there is at least one value beyond the integer upper limit.] ** In earlier versions of 3dcalc, scaling (if used) was applied to all sub-bricks equally -- a common scale factor was used. This would cause trouble if the values in different sub-bricks were in vastly different scales. In this version, each sub-brick gets its own scale factor. To override this behaviour, use the '-gscale' option. -gscale = Same as '-fscale', but also forces each output sub-brick to get the same scaling factor. This may be desirable for 3D+time datasets, for example. -nscale = Don't do any scaling on output to byte or short datasets. This may be especially useful when operating on mask datasets whose output values are only 0's and 1's. -a dname = Read dataset 'dname' and call the voxel values 'a' in the expression that is input below. 'a' may be any single letter from 'a' to 'z'. ** If some letter name is used in the expression, but not present in one of the dataset options here, then that variable is set to 0. ** If the letter is followed by a number, then that number is used to select the sub-brick of the dataset which will be used in the calculations. For example, '-b3 dname' specifies that the variable 'b' refers to sub-brick 3 of the dataset (indexes start at 0). N.B.: Another way to achieve the effect of '-b3' is described below in the 'INPUT DATASET SPECIFICATION' section. -expr "expression" Apply the expression within quotes to the input datasets, one voxel at a time, to produce the output dataset. ("sqrt(a*b)" to compute the geometric mean, for example). -prefix pname = Use 'pname' for the output dataset prefix name. [default='calc'] -session dir = Use 'dir' for the output dataset session directory. [default='./'=current working directory] -dt tstep = Use 'tstep' as the TR for manufactured 3D+time datasets. -TR tstep If not given, defaults to 1 second. -taxis N = If only 3D datasets are input (no 3D+time or .1D files), *OR* = then normally only a 3D dataset is calculated. With this -taxis N:tstep= option, you can force the creation of a time axis of length 'N', optionally using time step 'tstep'. In such a case, you will probably want to use the predefined time variables 't' and/or 'k' in your expression, or each resulting sub-brick will be identical. For example, '-taxis 121:0.1' will produce 121 points in time, spaced with TR 0.1. N.B.: You can also specify the TR using the -dt option. N.B.: You can specify 1D input datasets using the '1D:n@val,n@val' notation to get a similar effect; for example: -dt 0.1 -w '1D:121@0' will have pretty much the same effect as -taxis 121:0.1 N.B.: For both '-dt' and '-taxis', the 'tstep' value is in seconds. You can suffix it with 'ms' to specify that the value is in milliseconds instead; e.g., '-dt 2000ms'. -rgbfac A B C = For RGB input datasets, the 3 channels (r,g,b) are collapsed to one for the purposes of 3dcalc, using the formula value = A*r + B*g + C*b The default values are A=0.299 B=0.587 C=0.114, which gives the grayscale intensity. To pick out the Green channel only, use '-rgbfac 0 1 0', for example. Note that each channel in an RGB dataset is a byte in the range 0..255. Thus, '-rgbfac 0.001173 0.002302 0.000447' will compute the intensity rescaled to the range 0..1.0 (i.e., 0.001173=0.299/255, etc.) 3D+TIME DATASETS: This version of 3dcalc can operate on 3D+time datasets. Each input dataset will be in one of these conditions: (A) Is a regular 3D (no time) dataset; or (B) Is a 3D+time dataset with a sub-brick index specified ('-b3'); or (C) Is a 3D+time dataset with no sub-brick index specified ('-b'). If there is at least one case (C) dataset, then the output dataset will also be 3D+time; otherwise it will be a 3D dataset with one sub-brick. When producing a 3D+time dataset, datasets in case (A) or (B) will be treated as if the particular brick being used has the same value at each point in time. Multi-brick 'bucket' datasets may also be used. Note that if multi-brick (bucket or 3D+time) datasets are used, the lowest letter dataset will serve as the template for the output; that is, '-b fred+tlrc' takes precedence over '-c wilma+tlrc'. (The program 3drefit can be used to alter the .HEAD parameters of the output dataset, if desired.) INPUT DATASET NAMES ------------------- An input dataset is specified using one of these forms: 'prefix+view', 'prefix+view.HEAD', or 'prefix+view.BRIK'. You can also add a sub-brick selection list after the end of the dataset name. This allows only a subset of the sub-bricks to be read in (by default, all of a dataset's sub-bricks are input). A sub-brick selection list looks like one of the following forms: fred+orig[5] ==> use only sub-brick #5 fred+orig[5,9,17] ==> use #5, #9, and #12 fred+orig[5..8] or [5-8] ==> use #5, #6, #7, and #8 fred+orig[5..13(2)] or [5-13(2)] ==> use #5, #7, #9, #11, and #13 Sub-brick indexes start at 0. You can use the character '$' to indicate the last sub-brick in a dataset; for example, you can select every third sub-brick by using the selection list fred+orig[0..$(3)] N.B.: The sub-bricks are read in the order specified, which may not be the order in the original dataset. For example, using fred+orig[0..$(2),1..$(2)] will cause the sub-bricks in fred+orig to be input into memory in an interleaved fashion. Using fred+orig[$..0] will reverse the order of the sub-bricks. N.B.: You may also use the syntaxafter the name of an input dataset to restrict the range of values read in to the numerical values in a..b, inclusive. For example, fred+orig[5..7]<100..200> creates a 3 sub-brick dataset with values less than 100 or greater than 200 from the original set to zero. If you use the <> sub-range selection without the [] sub-brick selection, it is the same as if you had put [1..$] in front of the sub-range selection. N.B.: Datasets using sub-brick/sub-range selectors are treated as: - 3D+time if the dataset is 3D+time and more than 1 brick is chosen - otherwise, as bucket datasets (-abuc or -fbuc) (in particular, fico, fitt, etc datasets are converted to fbuc!) N.B.: The characters '$ ( ) [ ] < >' are special to the shell, so you will have to escape them. This is most easily done by putting the entire dataset plus selection list inside forward single quotes, as in 'fred+orig[5..7,9]', or double quotes "x". ** WARNING: you cannot combine sub-brick selection of the form -b3 bambam+orig (the old method) with sub-brick selection of the form -b 'bambam+orig[3]' (the new method) If you try, the Doom of Mandos will fall upon you! 1D TIME SERIES: You can also input a '*.1D' time series file in place of a dataset. In this case, the value at each spatial voxel at time index n will be the same, and will be the n-th value from the time series file. At least one true dataset must be input. If all the input datasets are 3D (single sub-brick) or are single sub-bricks from multi-brick datasets, then the output will be a 'manufactured' 3D+time dataset For example, suppose that 'a3D+orig' is a 3D dataset: 3dcalc -a a3D+orig -b b.1D -expr "a*b" The output dataset will 3D+time with the value at (x,y,z,t) being computed by a3D(x,y,z)*b(t). The TR for this dataset will be set to 'tstep' seconds -- this could be altered later with program 3drefit. Another method to set up the correct timing would be to input an unused 3D+time dataset -- 3dcalc will then copy that dataset's time information, but simply do not use that dataset's letter in -expr. If the *.1D file has multiple columns, only the first read will be used in this program. You can select a column to be the first by using a sub-vector selection of the form 'b.1D[3]', which will choose the 4th column (since counting starts at 0). '{...}' row selectors can also be used - see the output of '1dcat -help' for more details on these. Note that if multiple timeseries or 3D+time or 3D bucket datasets are input, they must all have the same number of points along the 'time' dimension. '1D:' INPUT: You can input a 1D time series 'dataset' directly on the command line, without an external file. The 'filename for such input takes the general format '1D:n_1@val_1,n_2@val_2,n_3@val_3,...' where each 'n_i' is an integer and each 'val_i' is a float. For example -a '1D:5@0,10@1,5@0,10@1,5@0' specifies that variable 'a' be assigned to a 1D time series of 35, alternating in blocks between values 0 and value 1. COORDINATES and PREDEFINED VALUES: If you don't use '-x', '-y', or '-z' for a dataset, then the voxel spatial coordinates will be loaded into those variables. For example, the expression 'a*step(x*x+y*y+z*z-100)' will zero out all the voxels inside a 10 mm radius of the origin x=y=z=0. Similarly, the '-t' value, if not otherwise used by a dataset or *.1D input, will be loaded with the voxel time coordinate, as determined from the header file created for the OUTPUT. Please note that the units of this are variable; they might be in milliseconds, seconds, or Hertz. In addition, slices of the dataset might be offset in time from one another, and this is allowed for in the computation of 't'. Use program 3dinfo to find out the structure of your datasets, if you are not sure. If no input datasets are 3D+time, then the effective value of TR is tstep in the output dataset, with t=0 at the first sub-brick. Similarly, the '-i', '-j', and '-k' values, if not otherwise used, will be loaded with the voxel spatial index coordinates. The '-l' (letter 'ell') value will be loaded with the temporal index coordinate. Otherwise undefined letters will be set to zero. In the future, new default values for other letters may be added. DIFFERENTIAL SUBSCRIPTS [22 Nov 1999]: Normal calculations with 3dcalc are strictly on a per-voxel basis: there is no 'cross-talk' between spatial or temporal locations. The differential subscript feature allows you to specify variables that refer to different locations, relative to the base voxel. For example, -a fred+orig -b 'a[1,0,0,0]' -c 'a[0,-1,0,0]' -d 'a[0,0,2,0]' means: symbol 'a' refers to a voxel in dataset fred+orig, symbol 'b' refers to the following voxel in the x-direction, symbol 'c' refers to the previous voxel in the y-direction symbol 'd' refers to the 2nd following voxel in the z-direction To use this feature, you must define the base dataset (e.g., 'a') first. Then the differentially subscripted symbols are defined using the base dataset symbol followed by 4 integer subscripts, which are the shifts in the x-, y-, z-, and t- (or sub-brick index) directions. For example, -a fred+orig -b 'a[0,0,0,1]' -c 'a[0,0,0,-1]' -expr 'median(a,b,c)' will produce a temporal median smoothing of a 3D+time dataset (this can be done more efficiently with program 3dTsmooth). Note that the physical directions of the x-, y-, and z-axes depend on how the dataset was acquired or constructed. See the output of program 3dinfo to determine what direction corresponds to what axis. For convenience, the following abbreviations may be used in place of some common subscript combinations: [1,0,0,0] == +i [-1, 0, 0, 0] == -i [0,1,0,0] == +j [ 0,-1, 0, 0] == -j [0,0,1,0] == +k [ 0, 0,-1, 0] == -k [0,0,0,1] == +l [ 0, 0, 0,-1] == -l The median smoothing example can thus be abbreviated as -a fred+orig -b a+l -c a-l -expr 'median(a,b,c)' When a shift calls for a voxel that is outside of the dataset range, one of three things can happen: STOP => shifting stops at the edge of the dataset WRAP => shifting wraps back to the opposite edge of the dataset ZERO => the voxel value is returned as zero Which one applies depends on the setting of the shifting mode at the time the symbol using differential subscripting is defined. The mode is set by one of the switches '-dsSTOP', '-dsWRAP', or '-dsZERO'. The default mode is STOP. Suppose that a dataset has range 0..99 in the x-direction. Then when voxel 101 is called for, the value returned is STOP => value from voxel 99 [didn't shift past edge of dataset] WRAP => value from voxel 1 [wrapped back through opposite edge] ZERO => the number 0.0 You can set the shifting mode more than once - the most recent setting on the command line applies when a differential subscript symbol is encountered. PROBLEMS: * Complex-valued datasets cannot be processed. * This program is not very efficient (but is faster than it once was). * Differential subscripts slow the program down even more. EXPRESSIONS: Arithmetic expressions are allowed, using + - * / ** and parentheses. As noted above, datasets are referred to by single letter variable names. At this time, C relational, boolean, and conditional expressions are NOT implemented. Built in functions include: sin , cos , tan , asin , acos , atan , atan2, sinh , cosh , tanh , asinh , acosh , atanh , exp , log , log10, abs , int , sqrt , max , min , J0 , J1 , Y0 , Y1 , erf , erfc , qginv, qg , rect , step , astep, bool , and , or , mofn , sind , cosd , tand , median, lmode , hmode , mad , gran , uran , iran , eran , lran , orstat, where * qg(x) = reversed cdf of a standard normal distribution * qginv(x) = inverse function to qg * min, max, atan2 each take 2 arguments ONLY * J0, J1, Y0, Y1 are Bessel functions (see Watson) * erf, erfc are the error and complementary error functions * sind, cosd, tand take arguments in degrees (vs. radians) * median(a,b,c,...) computes the median of its arguments * mad(a,b,c,...) computes the MAD of its arguments * orstat(n,a,b,c,...) computes the n-th order statistic of {a,b,c,...} - that is, the n-th value in size, starting at the bottom (e.g., orstat(1,a,b,c) is the minimum) * lmode(a,b,c,...) and hmode(a,b,c,...) compute the mode of their arguments - lmode breaks ties by choosing the smallest value with the maximal count, hmode breaks ties by choosing the largest value with the maximal count [median,lmode,hmode take a variable number of arguments] * gran(m,s) returns a Gaussian deviate with mean=m, stdev=s * uran(r) returns a uniform deviate in the range [0,r] * iran(t) returns a random integer in the range [0..t] * eran(s) returns an exponentially distributed deviate * lran(t) returns a logistically distributed deviate You may use the symbol 'PI' to refer to the constant of that name. This is the only 2 letter symbol defined; all input files are referred to by 1 letter symbols. The case of the expression is ignored (in fact, it is converted to uppercase as the first step in the parsing algorithm). The following functions are designed to help implement logical functions, such as masking of 3D volumes against some criterion: step(x) = {1 if x>0 , 0 if x<=0}, astep(x,y) = {1 if abs(x) > y , 0 otherwise} = step(abs(x)-y) rect(x) = {1 if abs(x)<=0.5, 0 if abs(x)>0.5}, bool(x) = {1 if x != 0.0 , 0 if x == 0.0}, and(a,b,...,c) = {1 if all arguments are nonzero, 0 if any are zero} or(a,b,...,c) = {1 if any arguments are nonzero, 0 if all are zero} mofn(m,a,...,c) = {1 if at least 'm' arguments are nonzero, 0 otherwise} argmax(a,b,...) = index of largest argument; = 0 if all args are 0 argnum(a,b,...) = number of nonzero arguments [These last 5 functions take a variable number of arguments.] The following 27 new [Mar 1999] functions are used for statistical conversions, as in the program 'cdf': fico_t2p(t,a,b,c), fico_p2t(p,a,b,c), fico_t2z(t,a,b,c), fitt_t2p(t,a) , fitt_p2t(p,a) , fitt_t2z(t,a) , fift_t2p(t,a,b) , fift_p2t(p,a,b) , fift_t2z(t,a,b) , fizt_t2p(t) , fizt_p2t(p) , fizt_t2z(t) , fict_t2p(t,a) , fict_p2t(p,a) , fict_t2z(t,a) , fibt_t2p(t,a,b) , fibt_p2t(p,a,b) , fibt_t2z(t,a,b) , fibn_t2p(t,a,b) , fibn_p2t(p,a,b) , fibn_t2z(t,a,b) , figt_t2p(t,a,b) , figt_p2t(p,a,b) , figt_t2z(t,a,b) , fipt_t2p(t,a) , fipt_p2t(p,a) , fipt_t2z(t,a) . See the output of 'cdf -help' for documentation on the meanings of and arguments to these functions. (After using one of these, you may wish to use program '3drefit' to modify the dataset statistical auxiliary parameters.) Computations are carried out in double precision before being truncated to the final output 'datum'. Note that the quotes around the expression are needed so the shell doesn't try to expand * characters, or interpret parentheses. (Try the 'ccalc' program to see how the expression evaluator works. The arithmetic parser and evaluator is written in Fortran-77 and is derived from a program written long ago by RW Cox to facilitate compiling on an array processor hooked up to a VAX. It's a mess, but it works - somewhat slowly.)