This program calculates the depth of ROIs, masks and 'background', using
the fun Euclidean Distance Transform (EDT).

Basically, this means calculating the Euclidean distance of each
voxel's centroid to the nearest boundary  with a separate ROI (well, to be
brutally technical, to centroid of the nearest voxel in a neighboring ROI.
The input dataset should be a map of ROIs (so, integer-valued). The
EDT values are calculated throughout the entire FOV by default,
even in the zero/background regions (there is an option to control this).

written by: PA Taylor and P Lauren (SSCC, NIMH, NIH)


This code calculates the Euclidean Distance Transform (EDT) for 3D
volumes following this nice, efficient algorithm, by Felzenszwalb
and Huttenlocher (2012;  FH2012):

   Felzenszwalb PF, Huttenlocher DP (2012). Distance Transforms of
   Sampled Functions. Theory of Computing 8:415-428.

Thanks to C. Rorden for pointing this paper out and discussing it.

The current code here extends/tweaks the FH2012 algorithm to a more
general case of having several different ROIs present, for running
in 3D (trivial extension), and for having voxels of non-unity and
non-isotropic lengths.  It does this by utilizing the fact that at
its very heart, the FH2012 algorithm works line by line and can even
be thought of as working boundary-by-boundary.

Here, the zero-valued 'background' is also just treated like an ROI,
with one difference.  At a FOV boundary, the zero-valued
ROI/background is treated as open, so that the EDT value at each
'zero' voxel is always to one of the shapes within the FOV.  For
nonzero ROIs, one can treat the FOV boundary *either* as an ROI edge
(EDT value there will be 1 edge length) *or* as being open.

Command usage and option list

    3dDepthMap [options] -prefix PREF -input DSET


  -input DSET      :(req) input dataset

  -prefix PREF     :(req) output prefix name

  -mask  MASK      :mask dataset.  NB: this mask is only applied *after*
                    the EDT has been calculated.  Therefore, the boundaries
                    of this mask have no affect on the calculated distance
                    values, except for potentially zeroing some out at the

  -dist_sq         :by default, the output EDT volume contains distance
                    values.  By using this option, the output values are

 -ignore_voxdims   :this EDT algorithm works in terms of physical distance
                    and uses the voxel dimension info in each direction, by
                    default.  However, using this option will ignore voxel
                    size, producing outputs as if each voxel dimension was

  -rimify RIM      :instead of outputting a depthmap for each ROI, output
                    a map of each ROI's 'rim' voxels---that is, the boundary
                    layer or periphery up to thickness RIM---if RIM>0.
                    +  Note that RIM is applied to whatever kind of depth
                    information you are calculating: if you use '-dist_sq'
                    then the voxel's distance-squared value to the ROI edge
                    is compared with RIM; if using '-ignore_voxdims', then
                    the number-of-voxels to the edge is compared with RIM.
                    The depthmap thresholding is applied as:
                    +  When using this opt, any labeltable/atlastable
                    from the original should be passed along, as well.
                    +  A negative RIM value inverts the check, and the
                    output is kept if the depth info is:
                    NB: with a negative RIM value, it is possible an ROI
                    could disappear!

  -zeros_are_zero  :by default, EDT values are output for the full FOV,
                    even zero-valued regions.  If this option is used, EDT
                    values are only reported within the nonzero locations
                    of the input dataset.

  -zeros_are_neg   :if this option is used, EDT in the zero/background
                    of the input will be negative (def: they are positive).
                    This opt cannot be used if '-zeros_are_zero' is.

  -nz_are_neg      :if this option is used, EDT in the nonzero ROI regions
                    of the input will be negative (def: they are positive).

  -bounds_are_not_zero :this flag affects how FOV boundaries are treated for
                    nonzero ROIs: by default, they are viewed as ROI
                    boundaries (so the FOV is a closed boundary for an ROI,
                    as if the FOV were padded by an extra layer of zeros);
                    but when this option is used, the ROI behaves as if it
                    continued 'infinitely' at the FOV boundary (so it is
                    an open boundary).  Zero-valued ROIs (= background)
                    are not affected by this option.

  -only2D   SLI    :instead of running full 3D EDT, run just in 2D, per.
                    plane.  Provide the slice plane you want to run along
                    as the single argument SLI:
                       "axi"  -> for axial slice
                       "cor"  -> for coronal slice
                       "sag"  -> for sagittal slice

  -binary_only     :if the input is a binary mask or should be treated as
                    one (all nonzero voxels -> 1; all zeros stay 0), then
                    using this option will speed up the calculation.  See
                    Notes below for more explanation of this. NOT ON YET!

 -verb V           :manage verbosity when running code (def: 1).
                    Providing a V of 0 means to run quietly.


Depth and the Euclidean Distance Transform

The original EDT algorithm of FH2012 was developed for a simple binary
mask input (and actually for homogeneous data grids of spacing=1). This
program, however, was built to handle more generalized cases of inputs,
namely ROI maps (and arbitrary voxel dimensions).

The tradeoff of the expansion to handling ROI maps is an increase in
processing time---the original binary-mask algorithm is *very* efficient,
and the generalized one is still pretty quick but less so.

So, if you know that your input should be treated as a binary mask, then
you can use the '-binary_only' option to utilize the more efficient
(and less generalized) algorithm.  The output dataset should be the same
in either case---this option flag is purely about speed of computation.

All other options about outputting dist**2 or negative values/etc. can be
used in conjunction with the '-binary_only', too.


1) Basic case:
   3dDepthMap                                                      \
       -input  roi_map.nii.gz                                      \
       -prefix roi_map_EDT.nii.gz

2) Same as above, but only output distances within nonzero regions/ROIs:
   3dDepthMap                                                      \
       -zeros_are_zero                                             \
       -input  roi_map.nii.gz                                      \
       -prefix roi_map_EDT_ZZ.nii.gz

3) Output distance-squared at each voxel:
   3dDepthMap                                                      \
       -dist_sq                                                    \
       -input  mask.nii.gz                                         \
       -prefix mask_EDT_SQ.nii.gz

4) Distinguish ROIs from nonzero background by making the former have
   negative distance values in output:
   3dDepthMap                                                      \
       -nz_are_neg                                                 \
       -input  roi_map.nii.gz                                      \
       -prefix roi_map_EDT_NZNEG.nii.gz

5) Have output voxel values represent (number of vox)**2 from a boundary;
   voxel dimensions are ignored here:
   3dDepthMap                                                      \
       -ignore_voxdims                                             \
       -dist_sq                                                    \
       -input  roi_map.nii.gz                                      \
       -prefix roi_map_EDT_SQ_VOX.nii.gz

6) Basic case, with option for speed-up because the input is a binary mask
   (i.e., only ones and zeros); any of the other above options can
   be combined with this, too:
   3dDepthMap                                                      \
       -binary_only                                                \
       -input  roi_mask.nii.gz                                     \
       -prefix roi_mask_EDT.nii.gz

7) Instead of outputting ROI depth, output a map of the ROI rims, keeping:
   the part of the ROI up to where depth is >=1.6mm
   3dDepthMap                                                      \
       -input  roi_map.nii.gz                                      \
       -rimify 1.6                                                 \
       -prefix roi_map_rim.nii.gz