Siddon’s Method
DRR Generation
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 total 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 term \(\|\mathbf p - \mathbf s\|_2\) endows the unit-free \(\mathrm d \alpha\) with the physical unit of length. 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 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.
siddon_raycast
siddon_raycast (source:torch.Tensor, target:torch.Tensor, volume:torch.Tensor, spacing:torch.Tensor, eps:float=1e-08)
An auto-differentiable implementation of the raycasting algorithm known as Siddon’s method.
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.