2D Distance transforms

u-Segment3D categorizes the distance transforms according to whether the limit or attractor of propagating points using gradient descent over an infinite number of steps is implicitly or explicitly defined (Extended Data Fig.1). Explicitly defined transforms are further categorized by the type of attractor: a single fixed point source or comprises a set of points.

u-Segment3D implements distance transforms, Φ that solve the Eikonal equation (∥∇Φ∥² = 1, which gives the shortest geodesic solution) or Poisson’s equation (∇²Φ = −1, which gives a smooth harmonic solution), for the cell interior using numerically stable methods. The Eikonal equation finds the shortest time of propagation for a point. Poisson’s equation can also be viewed as solving the shortest time of propagation but with the additional constraint of minimizing curvature, yielding smoother solutions.

With only the boundary condition Φ = 0, the Eikonal and Poisson equation conceptually propagates a wave inwards symmetrically from the cell boundaries. The limit solution is the definition of the medial axis skeleton, the locus of the centers of all inscribed maximal spheres of the object where these spheres touch the boundary at more than one point^77,124,125.

Solves the Eikonal equation using fast image morphological operations. u-Segment3D uses the memory and speed optimized implementation in the Python edt package released by the Seung Lab (https://github.com/seung-lab/euclidean-distance-transform-3d).

Solves the Poisson equation using LU decomposition (Python Scipy, scipy.sparse.linalg.spsolve) for each cell in an image. It is solved in parallel in u-Segment3D using the Python Dask library.

The implicit attractor solves the equations everywhere in the cell interior. The explicit attractor variants modifies the equations to have different source terms (right hand side of equation) in different parts of the cell interior. For the Eikonal equation, Φ = 0 at the cell boundary and outside, non-source points obey ∥∇Φ∥² = 1 and source points act as obstacles with vanishing speed, so that ∥∇Φ∥² = 0. For the Poisson equation, Φ = 0 at the cell boundary and outside, non-source points obey the Laplace equation, ∇²Φ = 0 and source points obey ∇²Φ = −1.

A single interior point is designated as a point source. u-Segment3D finds the interior point with Euclidean distance transform value greater than the percentile threshold (default: 10^th percentile) nearest the median coordinate of all points.

At the interior point, ∥∇Φ∥² = 0. The modified equations are solved using the Fast Marching Method (FMM)⁸⁰, with the constraint enforced using a masked array by the Python scikit-fmm library. Central first order differences are used to compute the unit normalized 2D gradient.

Only at the interior point, ∇²Φ = −1. The modified equations are solved using LU decomposition as before. To apply power transformation with exponent p > 0, the minimum is first subtracted from Φ to ensure positivity, Φ^p ≔ (Φ − Φ_min)^p. Central first order differences are used to compute the respective unit normalized 2D gradient.

Any number of interior points are designated as point sources. u-Segment3D computes the 2D medial axis skeleton as the point set. The binary skeleton is computed from the binary cell image by iteratively removing border pixels over multiple image passes¹²⁶ (Python Scikit-image, skimage.morphology.skeletonize). This raw result can often produce skeletons that have extraneous branches that may be too close to a neighboring, contacted cell. To improve the skeleton quality, the binary image is Gaussian filtered with σ = 3 pixels, rebinarized by mean value thresholding and reskeletonized.

For all points part of the skeleton, ∥∇Φ∥² = 0. The modified equations are then solved using the Fast Marching Method (FMM)⁸⁰ as above with central first order differences for computing the unit normalized 2D gradient. The gradients for all points part of the skeleton is set to zero to enforce the limiting behavior under gradient descent.

For all points part of the skeleton, ∇²Φ = −1. The modified equations are solved using LU decomposition as above with central first order differences for computing the unit normalized 2D gradient. The gradients for all points part of the skeleton is set to zero to enforce the limiting behavior under gradient descent.