AFNI program: 3dcalc

Output of -help

++ 3dcalc: AFNI version=AFNI_24.3.06 (Oct 31 2024) [64-bit]
++ Authored by: A cast of thousands
Program: 3dcalc                                                         
Author:  RW Cox et al                                                   

3dcalc - AFNI's calculator program ~1~                                  

     This program does voxel-by-voxel arithmetic on 3D datasets         
     (only limited inter-voxel computations are possible).              

     The program assumes that the voxel-by-voxel computations are being 
     performed on datasets that occupy the same space and have the same 
     orientations.                                                      

     3dcalc has a lot of input options, as its capabilities have grown  
     over the years.  So this 'help' output has gotten kind of long.    

     For simple voxel-wise averaging of datasets:    cf. 3dMean         
     For averaging along the time axis:              cf. 3dTstat        
     For smoothing in time:                          cf. 3dTsmooth      
     For statistics from a region around each voxel: cf. 3dLocalstat    

------------------------------------------------------------------------
Usage: ~1~                                                              
-----                                                                   
       3dcalc -a dsetA [-b dsetB...] \                                 
              -expr EXPRESSION       \                                 
              [options]                                                 

Examples: ~1~                                                           
--------                                                                
1. Average datasets together, on a voxel-by-voxel basis:                

     3dcalc -a fred+tlrc -b ethel+tlrc -c lucy+tlrc \
            -expr '(a+b+c)/3' -prefix subjects_mean                     

   Averaging datasets can also be done by programs 3dMean and 3dmerge.  
   Use 3dTstat to averaging across sub-bricks in a single dataset.      

2. Perform arithmetic calculations between the sub-bricks of a single   
   dataset by noting the sub-brick number on the command line:          

     3dcalc -a 'func+orig[2]' -b 'func+orig[4]' -expr 'sqrt(a*b)'       

3. Create a simple mask that consists only of values in sub-brick #0    
   that are greater than 3.14159:                                       

     3dcalc -a 'func+orig[0]' -expr 'ispositive(a-3.14159)' \
            -prefix mask                                                

4. Normalize subjects' time series datasets to percent change values in 
   preparation for group analysis:                                      

   Voxel-by-voxel, the example below divides each intensity value in    
   the time series (epi_r1+orig) with the voxel's mean value (mean+orig)
   to get a percent change value. The 'ispositive' command will ignore  
   voxels with mean values less than 167 (i.e., they are labeled as     
  'zero' in the output file 'percent_change+orig') and are most likely  
   background/noncortical voxels.                                       

     3dcalc -a epi_run1+orig -b mean+orig     \
            -expr '100 * a/b * ispositive(b-167)' -prefix percent_chng  

5. Create a compound mask from a statistical dataset, where 3 stimuli   
   show activation.                                                     
      NOTE: 'step' and 'ispositive' are identical expressions that can  
            be used interchangeably:                                    

     3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]' \
            -expr 'step(a-4.2)*step(b-2.9)*step(c-3.1)'              \
            -prefix compound_mask                                       

   In this example, all 3 statistical criteria must be met at once for  
   a voxel to be selected (value of 1) in this mask.                    

6. Same as example #5, but this time create a mask of 8 different values
   showing all combinations of activations (i.e., not only where        
   everything is active, but also each stimulus individually, and all   
   combinations).  The output mask dataset labels voxel values as such: 

        0 = none active    1 = A only active    2 = B only active       
        3 = A and B only   4 = C only active    5 = A and C only        
        6 = B and C only   7 = all A, B, and C active                   

     3dcalc -a 'func+orig[12]' -b 'func+orig[15]' -c 'func+orig[18]' \
            -expr 'step(a-4.2)+2*step(b-2.9)+4*step(c-3.1)'          \
            -prefix mask_8                                              

   In displaying such a binary-encoded mask in AFNI, you would probably 
   set the color display to have 8 discrete levels (the '#' menu).      

7. Create a region-of-interest mask comprised of a 3-dimensional sphere.
   Values within the ROI sphere will be labeled as '1' while values     
   outside the mask will be labeled as '0'. Statistical analyses can    
   then be done on the voxels within the ROI sphere.                    

   The example below puts a solid ball (sphere) of radius 3=sqrt(9)     
   about the point with coordinates (x,y,z)=(20,30,70):                 

     3dcalc -a anat+tlrc                                              \
            -expr 'step(9-(x-20)*(x-20)-(y-30)*(y-30)-(z-70)*(z-70))' \
            -prefix ball                                                

   The spatial meaning of (x,y,z) is discussed in the 'COORDINATES'     
   section of this help listing (far below).                            

8. Some datasets are 'short' (16 bit) integers with a scalar attached,   
   which allow them to be smaller than float datasets and to contain    
   fractional values.                                                   

   Dataset 'a' is always used as a template for the output dataset. For 
   the examples below, assume that datasets d1+orig and d2+orig consist 
   of small integers.                                                   

   a) When dividing 'a' by 'b', the result should be scaled, so that a  
      value of 2.4 is not truncated to '2'. To avoid this truncation,   
      force scaling with the -fscale option:                            

        3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -fscale   

   b) If it is preferable that the result is of type 'float', then set  
      the output data type (datum) to float:                            

        3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot \
                -datum float                                            

   c) Perhaps an integral division is desired, so that 9/4=2, not 2.24. 
      Force the results not to be scaled (opposite of example 8a) using 
      the -nscale option:                                               

        3dcalc -a d1+orig -b d2+orig -expr 'a/b' -prefix quot -nscale   

9. Compare the left and right amygdala between the Talairach atlas,     
   and the CA_N27_ML atlas.  The result will be 1 if TT only, 2 if CA   
   only, and 3 where they overlap.                                      

     3dcalc -a 'TT_Daemon::amygdala' -b 'CA_N27_ML::amygdala' \
            -expr 'step(a)+2*step(b)'  -prefix compare.maps             

   (see 'whereami_afni -help' for more information on atlases)               

10. Convert a dataset from AFNI short format storage to NIfTI-1 floating
    point (perhaps for input to an non-AFNI program that requires this):

      3dcalc -a zork+orig -prefix zfloat.nii -datum float -expr 'a'     

    This operation could also be performed with program 3dAFNItoNIFTI.  

11. Compute the edge voxels of a mask dataset.  An edge voxel is one    
    that shares some face with a non-masked voxel.  This computation    
    assumes 'a' is a binary mask (particularly for 'amongst').          

      3dcalc -a mask+orig -prefix edge                     \
             -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k     \
             -expr 'a*amongst(0,b,c,d,e,f,g)'                           

    consider similar erode or dilate operations:                        
        erosion:  -expr 'a*(1-amongst(0,b,c,d,e,f,g))'                  
        dilation: -expr 'amongst(1,a,b,c,d,e,f,g)'                      

------------------------------------------------------------------------
ARGUMENTS for 3dcalc (must be included on command line): ~1~            
---------                                                               

 -a dname    = Read dataset 'dname' and call the voxel values 'a' in the
               expression (-expr) that is input below. Up to 26 dnames  
               (-a, -b, -c, ... -z) can be included in a single 3dcalc  
               calculation/expression.                                  
               ** 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.                            
               ** You can use the subscript '[]' method                 
                  to select sub-bricks of datasets, as in               
                     -b dname+orig'[3]'                                 
               ** If you just want to test some 3dcalc expression,      
                  you can supply a dataset 'name' of the form           
                    jRandomDataset:64,64,16,40                          
                  to have the program create and use a dataset          
                  with a 3D 64x64x16 grid, with 40 time points,         
                  filled with random numbers (uniform on [-1,1]).       

 -expr       = Apply the expression - within quotes - to the input      
               datasets (dnames), one voxel at time, to produce the     
               output dataset.                                          
               ** You must use 1 and only 1 '-expr' option!             

 NOTE: If you want to average or sum up a lot of datasets, programs     
       3dTstat and/or 3dMean and/or 3dmerge are better suited for these 
       purposes.  A common request is to increase the number of input   
       datasets beyond 26, but in almost all cases such users simply    
       want to do simple addition!                                      

 NOTE: If you want to include shell variables in the expression (or in  
       the dataset sub-brick selection), then you should use double     
       "quotes" and the '$' notation for the shell variables; this    
       example uses csh notation to set the shell variable 'z':         

         set z = 3.5                                                    
         3dcalc -a moose.nii -prefix goose.nii -expr "a*$z"

       The shell will not expand variables inside single 'quotes',      
       and 3dcalc's parser will not understand the '$' character.       

 NOTE: You can use the ccalc program to play with the expression        
       evaluator, in order to get a feel for how it works and           
       what it accepts.                                                 

------------------------------------------------------------------------
 OPTIONS for 3dcalc: ~1~                                                
 -------                                                                

  -help      = Show this help.

  -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]                 
  -float }                                                              
  -short }   = Alternative options to specify output data format.       
  -byte  }                                                              

  -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 behavior,
                   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.                       
            ** N.B.: -usetemp and -gscale are incompatible!!            

  -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.       
                  ** Only use this option if you are sure you           
                     want the output dataset to be integer-valued!      

  -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]              
                  You can also include the output directory in the      
                  'pname' parameter to the -prefix option.              

  -usetemp      = With this option, a temporary file will be created to 
                  hold intermediate results.  This will make the program
                  run slower, but can be useful when creating huge      
                  datasets that won't all fit in memory at once.        
                * The program prints out the name of the temporary      
                  file; if 3dcalc crashes, you might have to delete     
                  this file manually.                                   
               ** N.B.: -usetemp and -gscale are incompatible!!         

  -dt tstep     = Use 'tstep' as the TR for "manufactured" 3D+time    
    *OR*          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  
  -taxis N:tstep: this 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 pre-   
                  defined 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.                                             

  -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.)                      

  -cx2r METHOD  = For complex input datasets, the 2 channels must be    
                  converted to 1 real number for calculation.  The      
                  methods available are:  REAL  IMAG  ABS  PHASE        
                * The default method is ABS = sqrt(REAL^2+IMAG^2)       
                * PHASE = atan2(IMAG,REAL)                              
                * Multiple '-cx2r' options can be given:                
                    when a complex dataset is given on the command line,
                    the most recent previous method will govern.        
                    This also means that for -cx2r to affect a variable 
                    it must precede it. For example, to compute the     
                    phase of data in 'a' you should use                 
                 3dcalc -cx2r PHASE -a dft.lh.TS.niml.dset -expr 'a'    
                    However, the -cx2r option will have no effect in    
                 3dcalc -a dft.lh.TS.niml.dset -cx2r PHASE -expr 'a'    
                    which will produce the default ABS of 'a'           
                    The -cx2r option in the latter example only applies 
                    to variables that will be defined after it.         
                    When in doubt, check your output.                   
                * If a complex dataset is used in a differential        
                    subscript, then the most recent previous -cx2r      
                    method applies to the extraction; for example       
                      -cx2r REAL -a cx+orig -cx2r IMAG -b 'a[0,0,0,0]'  
                    means that variable 'a' refers to the real part     
                    of the input dataset and variable 'b' to the        
                    imaginary part of the input dataset.                
                * 3dcalc cannot be used to CREATE a complex dataset!    
                    [See program 3dTwotoComplex for that purpose.]      

  -sort         = Sort each output brick separately, before output:     
  -SORT           'sort' ==> increasing order, 'SORT' ==> decreasing.   
                  [This is useful only under unusual circumstances!]    
                  [Sorting is done in spatial indexes, not in time.]    
                  [Program 3dTsort will sort voxels along time axis]    

  -isola        = After computation, remove isolated non-zero voxels.   
                  This option can be repeated to iterate the process;   
                  each copy of '-isola' will cause the isola removal    
                  process to be repeated one more time.

------------------------------------------------------------------------
DATASET TYPES: ~1~                                                      
-------------                                                           

 The most common AFNI dataset types are 'byte', 'short', and 'float'.   

 A byte value is an 8-bit signed integer (0..255), a short value ia a   
 16-bit signed integer (-32768..32767), and a float value is a 32-bit   
 real number.  A byte value has almost 3 decimals of accuracy, a short  
 has almost 5, and a float has approximately 7 (from a 23+1 bit         
 mantissa).                                                             

 Datasets can also have a scalar attached to each sub-brick. The main   
 use of this is allowing a short type dataset to take on non-integral   
 values, while being half the size of a float dataset.                  

 As an example, consider a short dataset with a scalar of 0.0001. This  
 could represent values between -32.768 and +32.767, at a resolution of 
 0.001.  One could represent the difference between 4.916 and 4.917, for 
 instance, but not 4.9165. Each number has 15 bits of accuracy, plus a  
 sign bit, which gives 4-5 decimal places of accuracy. If this is not   
 enough, then it makes sense to use the larger type, float.             

------------------------------------------------------------------------
3D+TIME DATASETS: ~1~                                                   
----------------                                                        

 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 ('[3]'); 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 #17
   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 syntax <a..b> after 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 [0..$] 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".

CATENATED AND WILDCARD DATASET NAMES
------------------------------------
 Datasets may also be catenated or combined in memory, as if one first
 ran 3dTcat or 3dbucket.
 
 An input with space-separated elements will be read as a concatenated
 dataset, as with 'dset1+tlrc dset2+tlrc dset3+tlrc', or with paths,
 'dir/dset1+tlrc dir/dset2+tlrc dir/dset3+tlrc'.
 The datasets will be combined (as if by 3dTcat) and then treated as a
 single input dataset.  Note that the quotes are required to specify
 them as a single argument.
 
 Sub-brick selection using '[]' works with space separated dataset
 names.  If the selector is at the end, it is considered global and
 applies to all inputs.  Otherwise, it applies to the adjacent input.
 For example:
    local:  'dset1+tlrc[2,3] dset2+tlrc[7,0,1] dset3+tlrc[5,0,$]'
    global: 'dset1+tlrc dset2+tlrc dset3+tlrc[5,6]'
 
 N.B. If AFNI_PATH_SPACES_OK is set to Yes, will be considered as part
 of the dataset name, and not as a separator between them.
 
 Similar treatment applies when specifying datasets using a wildcard
 pattern, using '*' or '?', as in: 'dset*+tlrc.HEAD'.  Any sub-brick
 selectors would apply to all matching datasets, as with:
    'dset*+tlrc.HEAD[2,5,3]'
 
 N.B.: complete filenames are required when using wildcard matching,
 or no files will exist to match, e.g. 'dset*+tlrc' would not work.
 
 N.B.: '[]' are processed as sub-brick or time point selectors.  They
 are therefore not allowed as wildcard characters in this context.
 
 Space and wildcard catenation can be put together.  In such a case,
 spaces divide the input into wildcard pieces, which are processed
 individually.
 
 Examples (each is processed as a single, combined dataset):
 
    'dset1+tlrc dset2+tlrc dset3+tlrc'
    'dset1+tlrc dset2+tlrc dset3+tlrc[2,5,3]'
    'dset1+tlrc[3] dset2+tlrc[0,1] dset3+tlrc[3,0,1]'
 
    'dset*+tlrc.HEAD'
    'dset*+tlrc.HEAD[2,5,3]'
    'dset1*+tlrc.HEAD[0,1] dset2*+tlrc.HEAD[7,8]'
 
    'group.*/subj.*/stats*+tlrc.HEAD[7]'

------------------------------------------------------------------------
1D TIME SERIES: ~1~                                                     
--------------                                                          

 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.                                     

 N.B.: To perform calculations ONLY on .1D files, use program 1deval.   
       3dcalc takes .1D files for use in combination with 3D datasets!  

 N.B.: If you auto-transpose a .1D file on the command line, (by ending 
       the filename with \'), then 3dcalc will NOT treat it as the     
       special case described above, but instead will treat it as       
       a normal dataset, where each row in the transposed input is a    
       'voxel' time series.  This would allow you to do differential    
       subscripts on 1D time series, which program 1deval does not      
       implement.  For example:                                         

        3dcalc -a '1D: 3 4 5 6'\' -b a+l -expr 'sqrt(a+b)' -prefix -   

       This technique allows expression evaluation on multi-column      
       .1D files, which 1deval also does not implement.  For example:   

        3dcalc -a '1D: 3 4 5 | 1 2 3'\' -expr 'cbrt(a)' -prefix -      

------------------------------------------------------------------------
'1D:' INPUT: ~1~                                                        
-----------                                                             

 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.                    

 You can combine 3dUndump with 3dcalc to create an all zero 3D+time     
 dataset from 'thin air', as in the commands                            
    3dUndump -dimen 128 128 32 -prefix AllZero_A -datum float           
    3dcalc -a AllZero_A+orig -b '1D: 100@' -expr 0 -prefix AllZero_B    
 If you replace the '0' expression with 'gran(0,1)', you'd get a        
 random 3D+time dataset, which might be useful for testing purposes.    

------------------------------------------------------------------------
'I:*.1D' and 'J:*.1D' and 'K:*.1D' INPUT: ~1~                           
----------------------------------------                                

 You can input a 1D time series 'dataset' to be defined as spatially    
 dependent instead of time dependent using a syntax like:               

   -c I:fred.1D                                                         

 This indicates that the n-th value from file fred.1D is to be associated
 with the spatial voxel index i=n (respectively j=n and k=n for 'J: and 
 K: input dataset names).  This technique can be useful if you want to  
 scale each slice by a fixed constant; for example:                     

   -a dset+orig -b K:slicefactor.1D -expr 'a*b'                         

 In this example, the '-b' value only varies in the k-index spatial     
 direction.                                                             

------------------------------------------------------------------------
COORDINATES and PREDEFINED VALUES: ~1~                                  
---------------------------------                                       

 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.

 The '-n' value, if not otherwise used, will be loaded with the overall 
 voxel 1D index.  For a 3D dataset, n = i + j*NX + k*NX*NY, where       
 NX, NY, NZ are the array dimensions of the 3D grid.  [29 Jul 2010]     

 Otherwise undefined letters will be set to zero.  In the future, new   
 default values for other letters may be added.                         

 NOTE WELL: By default, the coordinate order of (x,y,z) is the order in 
 *********  which the data array is stored on disk; this order is output
            by 3dinfo.  The options below control can change this order:

 -dicom }= Sets the coordinates to appear in DICOM standard (RAI) order,
 -RAI   }= (the AFNI standard), so that -x=Right, -y=Anterior , -z=Inferior,
                                        +x=Left , +y=Posterior, +z=Superior.

 -SPM   }= Sets the coordinates to appear in SPM (LPI) order,           
 -LPI   }=                      so that -x=Left , -y=Posterior, -z=Inferior,
                                        +x=Right, +y=Anterior , +z=Superior.

 The -LPI/-RAI behavior can also be achieved via the AFNI_ORIENT        
 environment variable (27 Aug, 2014).                                   
------------------------------------------------------------------------
DIFFERENTIAL SUBSCRIPTS [22 Nov 1999]: ~1~                              
-----------------------                                                 

 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.

N.B.: You can also use program 3dLocalstat to process data from a
      spatial neighborhood of each voxel; for example, to compute
      the maximum over a sphere of radius 9 mm placed around
      each voxel:
        3dLocalstat -nbhd 'SPHERE(9)' -stat max -prefix Amax9 A+orig

------------------------------------------------------------------------
ISSUES: ~1~
------ 

 * Complex-valued datasets cannot be processed, except via '-cx2r'.
 * This program is not very efficient (but is faster than it once was).
 * Differential subscripts slow the program down even more.

------------------------------------------------------------------------
------------------------------------------------------------------------
EXPRESSIONS: ~1~
----------- 

 As noted above, datasets are referred to by single letter variable names.
 Arithmetic expressions are allowed, using + - * / ** ^ and parentheses.
 C relational, boolean, and conditional expressions are NOT implemented!
* Note that the expression evaluator is designed not to fail;  illegal  *
* operations like 'sqrt(-1)' are changed to legal ones to avoid crashes.*
 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, mod  ,       
    mean , stdev, sem  , Pleg  , cbrt  , rhddc2, hrfbk4,hrfbk5
    minabove, maxbelow, extreme, absextreme    , acfwxm
    gamp , gampq

 where some of the less obvious functions are:
 * 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 the holy book: Watson)
 * Pleg(m,x) is the m'th Legendre polynomial evaluated at x
 * 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
 * mean(a,b,c,...) computes the mean of its arguments
 * stdev(a,b,c,...) computes the standard deviation of its arguments
 * sem(a,b,c,...) computes standard error of the mean of its arguments,
                  where sem(n arguments) = stdev(same)/sqrt(n)
 * 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)
 * minabove(X,a,b,c,...) computes the smallest value amongst {a,b,c,...}
    that is LARGER than the first argument X; if all values are smaller
    than X, then X will be returned
 * maxbelow(X,a,b,c,...) similarly returns the largest value amongst
    {a,b,c,...} that is SMALLER than the first argument X.
 * extreme(a,b,c,...) finds the largest absolute value amongst
    {a,b,c,...} returning one of the original a,b,c,... values.
 * absextreme(a,b,c,...) finds the largest absolute value amongst
    {a,b,c,...} returning the maximum absolute value of a,b,c,... values.
 * 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
    ["a,b,c,..." indicates 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
               with parameter s; mean=s
 * lran(t)   returns a logistically distributed deviate
               with parameter t; mean=0, stdev=t*1.814
 * mod(a,b)  returns (a modulo b) = a - b*int(a/b)
 * hrfbk4(t,L) and hrfbk5(t,L) are the BLOCK4 and BLOCK5 hemodynamic
    response functions from 3dDeconvolve (L=stimulus duration in sec,
    and t is the time in sec since start of stimulus); for example:
 1deval -del 0.1 -num 400 -expr 'hrfbk5(t-2,20)' | 1dplot -stdin -del 0.1
    These HRF functions are scaled to return values in the range [0..1]

 * ACFWXM(a,b,c,x) returns the Full Width at X Maximum for the mixed
   model ACF function
     f(r) = a*expr(-r*r/(2*b*b))+(1-a)*exp(-r/c)
   for X between 0 and 1 (not inclusive).  This is the model function
   estimated in program 3dFWHMx.
 * gamp(peak,fwhm) returns the parameter p in the formula
      g(t) = (t/(p*q))^p * exp(p-t/q)
   that gives the peak value of g(t) occurring at t=peak when the
   FWHM of g(t) is given by fwhm; gamq(peak,fwhm) gives the q parameter.
   These functions are largely used for creating FMRI hemodynamic shapes.

 You may use the symbol 'PI' to refer to the constant of that name.
 This is the only 2 letter symbol defined; all variables 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},
       posval(x)  = {x if x>0           , 0 if x<=0},
       astep(x,y) = {1 if abs(x) > y    , 0 otherwise} = step(abs(x)-y)
  within(x,MI,MX) = {1 if MI <= x <= MX , 0 otherwise},
       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},
    notzero(x)    = bool(x),
     iszero(x)    = 1-bool(x) = { 0 if x != 0.0, 1 if x == 0.0 },
        not(x)    = same as iszero(x)
     equals(x,y)  = 1-bool(x-y) = { 1 if x == y , 0 if x != y },
   ispositive(x)  = { 1 if x > 0; 0 if x <= 0 },
   isnegative(x)  = { 1 if x < 0; 0 if x >= 0 },
   ifelse(x,t,f)  = { t if x != 0; f if x == 0 },
        not(x)    = same as iszero(x) = Boolean negation
   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, else 0 }
  argmax(a,b,...) = index of largest argument; = 0 if all args are 0
  argnum(a,b,...) = number of nonzero arguments
  pairmax(a,b,...)= finds the 'paired' argument that corresponds to the
                    maximum of the first half of the input arguments;
                    for example, pairmax(a,b,c,p,q,r) determines which
                    of {a,b,c} is the max, then returns corresponding
                    value from {p,q,r}; requires even number of args.
  pairmin(a,b,...)= Similar to pairmax, but for minimum; for example,
                    pairmin(a,b,c,p,q,r} finds the minimum of {a,b,c}
                    and returns the corresponding value from {p,q,r};
                      pairmin(3,2,7,5,-1,-2,-3,-4) = -2
                    (The 'pair' functions are Lukas Pezawas specials!)
  amongst(a,b,...)= Return value is 1 if any of the b,c,... values
                    equals the a value; otherwise, return value is 0.
 choose(n,a,b,...)= chooses the n-th value from the a,b,... values.
                    (e.g., choose(2,a,b,c) is b)
  isprime(n)      = 1 if n is a positive integer and a prime number
                    0 if n is a positive integer and not a prime number
                   -1 if n is not a positive integer
                   or if n is bigger than 2^31-1

  [These last 9 functions take a variable number of arguments.]

 The following 27 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.  The two functions below use the
 NIfTI-1 statistical codes to map between statistical values and
 cumulative distribution values:
   cdf2stat(val,code,p1,p2,p3) -- val is between 0 and 1
   stat2cdf(val,code,p1,p2,p3) -- val is legal for the given distribution
 where code is
   2 = correlation statistic     p1 = DOF
   3 = t statistic (central)     p1 = DOF
   4 = F statistic (central)     p1 = num DOF, p2 = den DOF
   5 = N(0,1) statistic          no parameters (p1=p2=p3=0)
   6 = Chi-squared (central)     p1 = DOF
   7 = Beta variable (central)   p1 = a , p2 = b
   8 = Binomial variable         p1 = #trials, p2 = prob per trial
   9 = Gamma distribution        p1 = shape, p2 = scale
  10 = Poisson distribution      p1 = mean
  11 = N(mu,variance) normal     p1 = mean, p2 = scale
  12 = noncentral F statistic    p1 = num DOF, p2 = den DOF, p3 = noncen
  13 = noncentral chi-squared    p1 = DOF, p2 = noncentrality parameter
  14 = Logistic distribution     p1 = mean, p2 = scale
  15 = Laplace distribution      p1 = mean, p2 = scale
  16 = Uniform distribution      p1 = min, p2 = max
  17 = noncentral t statistic    p1 = DOF, p2 = noncentrality parameter
  18 = Weibull distribution      p1 = location, p2 = scale, p3 = power
  19 = Chi statistic (central)   p1 = DOF
  20 = inverse Gaussian variable p1 = mu, p2 = lambda
  21 = Extreme value type I      p1 = location, p2 = scale
  22 = 'p-value'                 no parameters
  23 = -ln(p)                    no parameters
  24 = -log10(p)                 no parameters
When fewer than 3 parameters are needed, the values for later parameters
are still required, but will be ignored.  An extreme case is code=5,
where the correct call is (e.g.) cdf2stat(p,5,0,0,0)

Finally, note that the expression evaluator is designed not to crash, or
to return NaN or Infinity.  Illegal operations, such as division by 0,
logarithm of negative value, etc., are intercepted and something else
(usually 0) will be returned.  To find out what that 'something else'
is in any specific case, you should play with the ccalc program.

** If you modify a statistical sub-brick, you may want 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 - but hey, computers are fast these days.)

++ Compile date = Oct 31 2024 {AFNI_24.3.06:linux_ubuntu_24_64}
This page auto-generated on Thu Oct 31 09:41:28 PM EDT 2024