Fast Auto-Differentiable Digitally Reconstructed Radiographs for Solving Inverse Problems in Intraoperative Imaging

The process of generating a DRR models the geometry of an idealized projectional radiography system (Fig. 1a). Let $s\in R^3$ be the X-ray source, and $p\in R^3$ be a target pixel on the detector plane. Then $R\left(\alpha \right)=s+\alpha (p-s)$ is a ray that originates from $s$ ($\alpha =0$), passes through the imaged volume, and hits the detector plane at $p$ ($\alpha =1$). The total energy attenuation experienced by the X-ray by the time it reaches pixel $p$ is given by the following line integral:

where $V:R^3\mapsto R$ is the imaged volume. The term $\parallel p-s{\parallel }_2$ endows the unit-free $d\alpha$ with the physical unit of length. For DRR synthesis, $V$ is approximated by a discrete 3D CT volume, and Eq. (1) becomes

where ${\alpha }_m$ parameterizes the locations where ray $R$ intersects one of the orthogonal planes comprising the CT volume, and $M$ is the number of such intersections (Fig. 1b). Note that this model does not account for patterns of reflection and scattering that are present in real X-ray systems. While these simplifications preclude synthesis of realistic X-rays, the model in Eq. (2) has been widely and successfully used in slice-to-volume registration [13]. Additionally, our approach of vectorizing DRR generation might also be interoperable with more sophisticated image synthesis models, an extension we examine further in the Discussion.

Siddon’s method [10] provides a parametric method to identify the plane intersections $\{{\alpha }_m{\}}_{m=1}^M$. Let $\Delta X$ be the CT voxel size in the $x$-direction and $b_x$ be the location of the $0$-th plane in this direction. Then the intersection of ray $R$ with the $i$-th plane in the $x$-direction is given by

with analogous expressions for ${\alpha }_y(\cdot )$ and ${\alpha }_z(\cdot )$. We can use Eq. (3) to compute the values ${\alpha }_x$ for all the intersections between $R$ and the planes in the $x$-direction:

where $i_{min}$ and $i_{max}$ denote the first and last intersections of $R$ with the $x$-direction planes. Defining ${\alpha }_y$ and ${\alpha }_z$ analogously, we construct the array

which contains $M$ values of $\alpha$ parameterizing the intersections between $R$ and the orthogonal $x$-, $y$-, and $z$-directional planes. We substitute values in the sorted set $\alpha$ into Eq. (2) to evaluate $E\left(R\right)$, which corresponds to the intensity of pixel $p$ in the synthesized DRR.

A faster variant determines consecutive intersecting planes iteratively [4]. For example, the value of $\alpha$ at the second plane intersected by $R$ is given by ${\alpha }_2=min\left\{{\alpha }_x\right(i_{min}+1),{\alpha }_y(j_{min}+1),{\alpha }_z(k_{min}+1\left)\right\}$. The algorithm iteratively finds the next value of $\alpha$ until we reach the edge of the CT volume, making this approach more memory efficient by requiring fewer intermediate values to be stored. This modified algorithm, known as Siddon-Jacobs’ method, is commonly implemented in CUDA and C++ to create multi-threaded GPU-accelerated DRR generators that exploit data parallelism by assigning each thread to trace an independent ray intersecting the detector plane [1, 5, 7].

While Siddon-Jacobs’ method is more memory efficient, the iterative loop it relies on is not amenable to implementations in vectorized tensor algebra libraries. Thus we vectorize the original Siddon’s method as follows. Let $P\in R^{H\times W\times 3}$ contain the 3D pixel coordinates of a DRR with dimension $H\times W$. We compute the $\alpha$ values for intersections with all of the $x$-, $y$-, and $z$-planes for all $p\in P$ in parallel:

where $(n_x,n_y,n_z)$ are the dimensions of the CT volume $V$, $(i,j,k)$ are the CT voxel indices, $(\Delta X,\Delta Y,\Delta Z)$ are the CT voxel sizes, and $\otimes$ and $\oslash$ are the Hadamard product and division operators, respectively. Rather than explicitly compute the indices $(i_{min},i_{max})$, $(j_{min},j_{max})$, and $(k_{min},k_{max})$ for each ray, as is done in Siddon’s original method, we instead compute the minimum and maximum values of $\alpha$, corresponding to when each ray enters and exits the volume:

where ${\alpha }_{min},{\alpha }_{max}\in R^{H\times W}$. We filter $A$ to include only values in the range $[{\alpha }_{min},{\alpha }_{max}]$ and sort each row $A(h,w,\cdot )$ for $h\in \{1,\dots ,H\},w\in \{1,\dots ,W\}$. Finally, we evaluate Eq. (2) with this sorted tensor to compute the intensity for each pixel in the DRR, completing a chain of vectorized tensor operations.

Because we reformulated the original Siddon’s method as a series of tensor operations, our vectorized version benefits from the mature GPU compilers and memory allocators developed for optimizing large-scale deep learning models. For empirical evaluation of our method, we also implemented a partially-vectorized version of Siddon-Jacobs’ method in which the $\alpha$ updates are still computed iteratively (i.e., with a loop), but the updates are applied in a vectorized form to every target pixel in the detector plane.

We specify the 3D position of the X-ray source and detector plane relative to the CT volume with the following seven geometric parameters: radius $\rho$ that acts as a scaling factor; three rotational degrees of freedom (DoF) $(\theta ,\phi ,\gamma )$; and three translational DoF $(b_x,b_y,b_z)$. Using spherical coordinates, we express the position of the X-ray source as $s=(\rho ,\theta ,\phi )$, where $\rho$ is half of the source-to-detector distance, and $\theta$ and $\phi$ are the azimuthal and polar angles, respectively. We assume the detector plane is tangent to this implied sphere at the point opposite $s$. The orientation of this plane is determined by a rotation about the $x$-axis by the angle $\gamma$. We add the translation $(b_x,b_y,b_z)$ to the coordinates of the X-ray source and detector plane to create a reference frame wherein the patient is not perfectly centered relative to the X-ray scanner.

Since every step of our pipeline, from the generation of the pixels on the detector plane to the computation of the pixel intensities, is performed in PyTorch’s tensor framework, the resulting DRRs are differentiable with respect to the parameters described above. That is, the gradient ${\nabla }_{\eta }L\left(I\right(\eta \left)\right)$ of the loss function $L(\cdot )$ evaluated for the DRR $I\left(\eta \right)$ with respect to the parameters $\eta =(\rho ,\theta ,\phi ,\gamma ,b_x,b_y,b_z)$ is obtained automatically for any differentiable $L(\cdot )$.

We compare our vectorized GPU version of Siddon’s method (VGS) to two baseline approaches: a widely used CPU implementation in the Plastimatch package (CP) [9], and our vectorized GPU implementation of Siddon-Jacobs’ method (VGSJ) which we described in Section 2.3. Note that the CUDA-accelerated DRR generator in Plastimatch is not working at the time of publication. GPU benchmarks are run on an NVIDIA GeForce RTX 2080 Ti, and CPU benchmarks are run on an 18-core Linux computer with Intel(R) Xeon(R) CPU E5-2697 v4 @ 2.30GHz processors.

We use our auto-differentiable DRR generator to implement slice-to-volume registration with synthetic DRRs. Specifically, we generate a fixed DRR from a set of ground truth parameters ${\eta }^{\ast }=(\theta ,\phi ,\gamma ,b_x,b_y,b_z)$, and generate a second moving DRR from a set of random initial parameters ${\eta }_0$. We use basic gradient descent to minimize the negative Zero-Normalized Cross-Correlation (ZNCC) between the fixed DRR and the moving DRR.