renderers
Siddon’s Method
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 the first equation 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 Siddon’s method has been widely and successfully used in slice-to-volume registration.
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} \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 this equation 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} \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 the first equation to evaluate \(E(R)\), which corresponds to the intensity of pixel \(\mathbf p\) in the synthesized DRR.
Siddon
Siddon (mode:str='nearest', stop_gradients_through_grid_sample:bool=False, filter_intersections_outside_volume:bool=True, reducefn:str='sum', eps:float=1e-08)
Differentiable X-ray renderer implemented with Siddon’s method for exact raytracing.
| Type | Default | Details | |
|---|---|---|---|
| mode | str | nearest | Interpolation mode for grid_sample |
| stop_gradients_through_grid_sample | bool | False | Apply torch.no_grad when calling grid_sample |
| filter_intersections_outside_volume | bool | True | Use alphamin/max to filter the intersections |
| reducefn | str | sum | Function for combining samples along each ray |
| eps | float | 1e-08 | Small constant to avoid div by zero errors |
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.
Now, the rendering equation is \[\begin{equation} E(R) = \|\mathbf p - \mathbf s\|_2\frac{\alpha_{\max} - \alpha_{\min}}{M-1} \sum_{m=1}^{M} \mathbf V \left[ \mathbf s + \alpha_m (\mathbf p - \mathbf s) \right] \,, \end{equation}\] where \(\mathbf V[\cdot]\) is the trilinear interpolation function and \(M\) is the number of points sampled per ray.
Trilinear
Trilinear (mode:str='bilinear', reducefn:str='sum', eps:float=1e-08)
Differentiable X-ray renderer implemented with trilinear interpolation.
| Type | Default | Details | |
|---|---|---|---|
| mode | str | bilinear | Interpolation mode for grid_sample |
| reducefn | str | sum | Function for combining samples along each ray |
| eps | float | 1e-08 | Small constant to avoid div by zero errors |