# Siddon’s Method

Mapping rays to pixel intensities

## 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

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.

source

### 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.