+ Module for computing digitally reconstructed radiographs
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
DRR
+
DRR is a PyTorch module that compues differentiable digitally reconstructed radiographs. The viewing angle for the DRR (known generally in computer graphics as the camera pose) is parameterized by the following parameters:
+
+
SDR : Source-to-Detector radius (half of the source-to-detector distance)
+
\(\mathbf R \in \mathrm{SO}(3)\) : a rotation
+
\(\mathbf t \in \mathbb R^3\) : a translation
+
+
+
+
+
+
+
+Tip
+
+
+
+
DiffDRR can take a rotation parameterized in any of the following forms to move the detector plane:
+
+
axis_angle
+
euler_angles (note: also need to specify the convention for the Euler angles)
(bx, by, bz) are translational parameters and (alpha, beta, gamma) are rotational parameters. The rotational parameters are detailed in Spherical Coordiantes Tutorial.
PyTorch module that computes differentiable digitally reconstructed radiographs.
+
+
+
+
+
+
+
+
+
+
+
Type
+
Default
+
Details
+
+
+
+
+
volume
+
np.ndarray
+
+
CT volume
+
+
+
spacing
+
np.ndarray
+
+
Dimensions of voxels in the CT volume
+
+
+
sdr
+
float
+
+
Source-to-detector radius for the C-arm (half of the source-to-detector distance)
+
+
+
height
+
int
+
+
Height of the rendered DRR
+
+
+
delx
+
float
+
+
X-axis pixel size
+
+
+
width
+
int | None
+
None
+
Width of the rendered DRR (default to height)
+
+
+
dely
+
float | None
+
None
+
Y-axis pixel size (if not provided, set to delx)
+
+
+
x0
+
float
+
0.0
+
Principal point X-offset
+
+
+
y0
+
float
+
0.0
+
Principal point Y-offset
+
+
+
p_subsample
+
float | None
+
None
+
Proportion of pixels to randomly subsample
+
+
+
reshape
+
bool
+
True
+
Return DRR with shape (b, 1, h, w)
+
+
+
reverse_x_axis
+
bool
+
False
+
If pose includes reflection (in E(3) not SE(3)), reverse x-axis
+
+
+
patch_size
+
int | None
+
None
+
Render patches of the DRR in series
+
+
+
bone_attenuation_multiplier
+
float
+
1.0
+
Contrast ratio of bone to soft tissue
+
+
+
renderer
+
str
+
siddon
+
Rendering backend, either “siddon” or “trilinear”
+
+
+
renderer_kwargs
+
+
+
+
+
+
+
The forward pass of the DRR module generated DRRs from the input CT volume. The pose parameters (i.e., viewing angles) from which the DRRs are generated are passed to the forward call.
Generate DRR with rotational and translational parameters.
+
+
+
+
Registration
+
The Registration module uses the DRR module to perform differentiable 2D-to-3D registration. Initial guesses for the pose parameters are as stored as nn.Parameters of the module. This allows the pose parameters to be optimized with any PyTorch optimizer. Furthermore, this design choice allows DRR to be used purely as a differentiable renderer.
We represent rigid transforms as \(4 \times 4\) matrices
+
\[\begin{equation}
+\begin{bmatrix}
+ \mathbf R & \mathbf t \\
+ \mathbf 0 & 1
+\end{bmatrix}
+\in \mathbf{SE}(3) \,,
+\end{equation}\]
+
where \(\mathbf R \in \mathbf{SO}(3)\) is a rotation matrix and \(\mathbf t\in \mathbb R^3\) represents a translation.
+
Note that since rotation matrices are orthogonal (\(\mathbf R \mathbf R^T = \mathbf R^T \mathbf R = \mathbf I_3\)), we have a simple closed-form equation for the inverse: \[\begin{equation}
+\begin{bmatrix}
+ \mathbf R & \mathbf t \\
+ \mathbf 0 & 1
+\end{bmatrix}^{-1} =
+\begin{bmatrix}
+ \mathbf R^T & -\mathbf R^T \mathbf t \\
+ \mathbf 0 & 1
+\end{bmatrix} \,.
+\end{equation}\]
Applies rigid transforms in SE(3) to point clouds. Can handle batched rigid transforms, composition of transforms, closed-form inversion, and conversions to various representations of SE(3).
Convert a 10-vector in the quaternion adjugate, a symmetric matrix whose eigenvector corresponding to the eigenvalue of maximum modulus is the (unnormalized) quaternion. Uses a fast method to solve for the eigenvector without explicity computing the eigendecomposition.
+
Source: https://arxiv.org/abs/2205.09116
+
+
+
PyTorch3D conversions port
+
PyTorch3D has many useful conversion functions for transforming between multiple parameterizations of \(\mathbf{SO}(3)\) and \(\mathbf{SE}(3)\). However, installing PyTorch3D can be annoying for users not on Linux. We include the required conversion functions for PyTorch3D below. The original LICENSE from PyTorch3D is also included:
+
BSD License
+
+For PyTorch3D software
+
+Copyright (c) Meta Platforms, Inc. and affiliates. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+
+ * Redistributions of source code must retain the above copyright notice, this
+ list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+
+ * Neither the name Meta nor the names of its contributors may be used to
+ endorse or promote products derived from this software without specific
+ prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
+ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
DRRs are generated by modeling the geometry of an idealized projectional radiography system. Let \(\mathbf s \in \mathbb R^3\) be the X-ray source and \(\mathbf p \in \mathbb R^3\) be a target pixel on the detector plane. Then, \(R(\alpha) = \mathbf s + \alpha (\mathbf p - \mathbf s)\) is a ray that originates from \(\mathbf s\) (\(\alpha=0\)), passes through the imaged volume, and hits the detector plane at \(\mathbf p\) (\(\alpha=1\)). The proportion of energy attenuation experienced by the X-ray at the time it reaches pixel \(\mathbf p\) is given by the following line integral:
+
\[\begin{equation}
+ E(R) = \|\mathbf p - \mathbf s\|_2 \int_0^1 \mathbf V \left( \mathbf s + \alpha (\mathbf p - \mathbf s) \right) \, \mathrm d\alpha \,,
+\end{equation}\]
+
where \(\mathbf V : \mathbb R^3 \mapsto \mathbb R\) is the imaged volume. The units term \(\|\mathbf p - \mathbf s\|_2\) serves to cancel out the units of \(\mathbf V(\cdot)\), reciprocal length, such that the final proportion \(E\) is unitless. For DRR synthesis, \(\mathbf V\) is approximated by a discrete 3D CT volume, and Eq. (1) becomes
+
\[\begin{equation}
+ E(R) = \|\mathbf p - \mathbf s\|_2 \sum_{m=1}^{M-1} (\alpha_{m+1} - \alpha_m) \mathbf V \left[ \mathbf s + \frac{\alpha_{m+1} + \alpha_m}{2} (\mathbf p - \mathbf s) \right] \,,
+\end{equation}\]
+
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.
+
+
+
+
+
+
+Note
+
+
+
+
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.
Modules can also contain other Modules, allowing to nest them in a tree structure. You can assign the submodules as regular attributes::
+
import torch.nn as nn
+import torch.nn.functional as F
+
+class Model(nn.Module):
+ def __init__(self):
+ super().__init__()
+ self.conv1 = nn.Conv2d(1, 20, 5)
+ self.conv2 = nn.Conv2d(20, 20, 5)
+
+ def forward(self, x):
+ x = F.relu(self.conv1(x))
+ return F.relu(self.conv2(x))
+
Submodules assigned in this way will be registered, and will have their parameters converted too when you call :meth:to, etc.
+
.. note:: As per the example above, an __init__() call to the parent class must be made before assignment on the child.
+
:ivar training: Boolean represents whether this module is in training or evaluation mode. :vartype training: bool
+
Siddon’s method 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 \[\begin{equation}
+ \label{eqn:x-intersect}
+ \alpha_x(i) = \frac{b_x + i \Delta X - \mathbf s_x}{\mathbf p_x - \mathbf s_x} ,
+\end{equation}\] with analogous expressions for \(\alpha_y(\cdot)\) and \(\alpha_z(\cdot)\).
+
We can use Eq. (3) to compute the values \(\mathbf \alpha_x\) for all the intersections between \(R\) and the planes in the \(x\)-direction: \[\begin{equation*}
+ \mathbf\alpha_x = \{ \alpha_x(i_{\min}), \dots, \alpha_x(i_{\max}) \} ,
+\end{equation*}\] where \(i_{\min}\) and \(i_{\max}\) denote the first and last intersections of \(R\) with the \(x\)-direction planes.
+
Defining \(\mathbf\alpha_y\) and \(\mathbf\alpha_z\) analogously, we construct the array \[\begin{equation}
+ \label{eqn:alphas}
+ \mathbf\alpha = \mathrm{sort}(\mathbf\alpha_x, \mathbf\alpha_y, \mathbf\alpha_z) ,
+\end{equation}\] 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 \(\mathbf\alpha\) into Eq. (2) to evaluate \(E(R)\), which corresponds to the intensity of pixel \(\mathbf p\) in the synthesized DRR.
+
+
+
+
Trilinear interpolation
+
Instead of computing the exact line integral over the voxel grid (i.e., Siddon’s method), we can sample colors at points along the each ray using trilinear interpolation.
df is a pandas.DataFrame with columns ["alpha", "beta", "gamma", "bx", "by", "bz"]. Each row in df is an iteration of optimization with the updated values for that timestep.
+
+
+
+
3D Visualization
+
Uses pyvista and trame to interactively visualize DRR geometry in 3D.
For a given pose (not batched), turn the camera and detector into a mesh. Additionally, render the DRR for the pose. Convert into a texture that can be applied to the detector mesh.
Auto-differentiable DRR synthesis and optimization in PyTorch
+
+
+
DiffDRR is a PyTorch-based digitally reconstructed radiograph (DRR) generator that provides
+
+
Auto-differentiable X-ray rendering
+
GPU-accelerated DRR synthesis
+
A pure Python implementation
+
+
Most importantly, DiffDRR implements DRR rendering as a PyTorch module, making it interoperable with deep learning pipelines.
+
Below is a comparison of DiffDRR to a real X-ray (X-rays and CTs from the DeepFluoro dataset):
+
+
+
+DiffDRR rendered from the same camera pose as a real X-ray.
+
+
+
+
Install
+
To install DiffDRR from PyPI:
+
pip install diffdrr
+
+
+
Usage
+
The following minimal example specifies the geometry of the projectional radiograph imaging system and traces rays through a CT volume:
+
+
import matplotlib.pyplot as plt
+import torch
+
+from diffdrr.drr import DRR
+from diffdrr.data import load_example_ct
+from diffdrr.visualization import plot_drr
+
+# Read in the volume and get the isocenter
+volume, spacing = load_example_ct()
+bx, by, bz = torch.tensor(volume.shape) * torch.tensor(spacing) /2
+
+# Initialize the DRR module for generating synthetic X-rays
+device ="cuda"if torch.cuda.is_available() else"cpu"
+drr = DRR(
+ volume, # The CT volume as a numpy array
+ spacing, # Voxel dimensions of the CT
+ sdr=300.0, # Source-to-detector radius (half of the source-to-detector distance)
+ height=200, # Height of the DRR (if width is not seperately provided, the generated image is square)
+ delx=4.0, # Pixel spacing (in mm)
+).to(device)
+
+# Set the camera pose with rotation (yaw, pitch, roll) and translation (x, y, z)
+rotation = torch.tensor([[torch.pi, 0.0, torch.pi /2]], device=device)
+translation = torch.tensor([[bx, by, bz]], device=device)
+
+# 📸 Also note that DiffDRR can take many representations of SO(3) 📸
+# For example, quaternions, rotation matrix, axis-angle, etc...
+img = drr(rotation, translation, parameterization="euler_angles", convention="ZYX")
+plot_drr(img, ticks=False)
+plt.show()
+
+
+
+
+
+
+
+
+
On a single NVIDIA RTX 2080 Ti GPU, producing such an image takes
+
+
+
+
33.6 ms ± 27.8 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
We demonstrate the utility of our auto-differentiable DRR generator by solving the 2D/3D registration problem with gradient-based optimization. Here, we generate two DRRs:
+
+
A fixed DRR from a set of ground truth parameters
+
A moving DRR from randomly initialized parameters
+
+
To solve the registration problem, we use gradient descent to maximize an image loss similarity metric between the two DRRs. This produces optimization runs like this:
DiffDRR reformulates Siddon’s method,1 the canonical algorithm for calculating the radiologic path of an X-ray through a volume, as a series of vectorized tensor operations. This version of the algorithm is easily implemented in tensor algebra libraries like PyTorch to achieve a fast auto-differentiable DRR generator.
