Three-dimensional multiscale discrete Radon and John transforms

TL;DR

The paper introduces two forward transforms for discrete 3D data: a multiscale discrete Radon transform (DRT) of planes and a multiscale discrete John transform (DJT) of lines, both operating with integer arithmetic and no interpolation. Each transform employs a recursive, dyadic subdivision (dyadic cubes for planes, dodecant-like partitions for lines) that yields linearithmic complexity and fast, exact adjoint-based inverses suitable for iterative refinement. The authors provide detailed mapping equations, output shapes, and implementation strategies, including a Halide-based, architecture-agnostic pipeline and a comprehensive comparison to Fourier-based Radon methods. They demonstrate practical applicability through a depth-map/scene-scanning example and discuss inversion stability, convergence via multigrid schemes, and potential for real-time processing on standard hardware. The results indicate substantial speed advantages for forward paths and integers-only computation, with implications for tomography, computational photography, and acoustic localization.

Abstract

Two algorithms are introduced for the computation of discrete integral transforms with a multiscale approach operating in discrete three-dimensional (3D) volumes while considering its real-time implementation. The first algorithm, referred to as 3D discrete Radon transform (DRT) of planes, will compute the summation set of values lying in discrete planes in a cube that imitates, in discrete data, the integrals on two-dimensional planes in a 3D volume similar to the continuous Radon transform. The normals of these planes, equispaced in ascents, cover a quadrilateralized hemisphere and comprise 12 dodecants. The second proposed algorithm, referred to as the 3D discrete John transform (DJT) of lines, will sum elements lying on discrete 3D lines while imitating the behavior of the John or X-ray continuous transform on 3D volumes. These discrete integral transforms do not perform interpolation on input or intermediate data, and they can be computed using only integer arithmetics with linearithmic complexity; thus, outperforming the methods based on the Fourier slice-projection theorem for real-time applications. We briefly prove that these transforms have fast inversion algorithms that are exact for discrete inputs.

Figures (16)

• Figure 1: Top left: Depiction of the relative sizes and shapes of an input and its 2D DRT for the quadrant covering $0^\circ$ to $45^\circ$. Top right: Set of displacements through an $8\times8$ domain for digital lines with slope 3. Bottom: DRT of the four quadrants covering $180^\circ$ merged together.
• Figure 2: Discrete plane $z = l^5_{4}(\mathbf{x}) + l^5_{12}(\mathbf{y})$, and ray $(x,y,z) = ( \lambda(\mathbf{u}), \, 26 - l^5_{23}(\mathbf{u}), \, 4 + l^5_{7}(\mathbf{u})\, )$
• Figure 3: Left: Two unique alternatives for a discrete line segment crossing the black dot in 2D at the lowest scale. Middle: Alternatives for a 3D line crossing through the black dot at the lowest scale. Right: Four unique alternatives for a discrete plane patch crossing the black dot in 3D at the lowest scale.
• Figure 4: Right: Depiction of the parameters of the basic algorithm of 3D DRT. Left: Boundary x-y coordinates of plane normals covered by the basic algorithm, as projected into positive and negative hemispheres cut through $z=0$.
• Figure 5: Depiction of the four dodecants that compute the entire set of planes near the $z-$normal: null plane (face in yellow), origin point $\mathbf{p_0}$ (black circle), direction of ascents (black arrow), direction of neighbors points $\mathbf{p_1}$ and $\mathbf{p_2}$ (red and blue arrows), and location of boundary $x-y$ coordinates for plane normals (dots) within positive and negative $z-$hemispheres (dashed circles). Left: Quadrilateralized spherical caps emerge from the z-faces, inspired from Chan & O'Neill (1975).