Siddon’s Method
DRR Generation
The process of generating a DRR models 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 by 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. 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_raycast
siddon_raycast (source:torch.Tensor, target:torch.Tensor, volume:torch.Tensor, spacing:torch.Tensor, eps:float=1e-08)
Compute 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.