Digitally Reconstructed Radiographs

The CT volume consists of a 3D array of x-ray attenuation values, captured by scanning a patient. CT scanners produce attenuation values in Hounsfield units, but we will assume that we have converted this to attenuation values in units of "fraction of x-ray absorbed per unit length", denoted $\mu$. Neither of these units is perfect, as they do not take into account that x-ray absorption varies with x-ray energy.

The CT volume has a known position and orientation (given by the transformation, $T$, from the previous lecture) in the coordinate system of the fluoro.

The output of the DRR algorithm is a simulated image, $I_\textrm{proj}$, on the image plane.

Recall the Beer-Lambert law from the lecture on x-ray imaging:

The received x-ray energy, $I_\textrm{out}$, is

$\displaystyle I_\textrm{out} = I_\textrm{in} \; e^{- \int_0^D \mu(t) \; dt }$

where $\mu(t)$ is the fraction of x-rays absorbed per unit length at position $t$ on the ray.

For each pixel of the image, we draw a ray, $\ell(t)$, from the fluoro origin to that pixel, going through the CT volume. We determine the point at which the ray enters the volume ($t = 0$) and the point at which the ray exits the volume ($t = D$). Note that the length, $D$, of the ray inside the volume, can be different for each pixel.

Given the entry and exit points, we can evalutate $e^{-\int_0^D \mu(t) \; dt}$ between those points, plug that result into the equation above, and get a value for the pixel.

The ray segment $\ell(0)$ to $\ell(D)$ is cut into small pieces, each of width $\Delta t$:

where $t_i = i \; \Delta t$ and $n = {D \over \Delta t}$.

Then the exponential over the whole segment can be broken into a product of exponentials over the intervals:

$\begin{array}{rcl} \displaystyle e^{-\int_0^D \mu(t) \; dt} \;\; & = & \displaystyle e^{-\int_{t_0}^{t_1} \mu(t) \; dt} \;\;\times\;\; e^{-\int_{t_1}^{t_2} \mu(t) \; dt} \;\;\times\;\; \cdots \;\;\times\;\; e^{-\int_{t_{n-1}}^{t_n} \mu(t) \; dt} \\ & = & \displaystyle \prod_{i=0}^{n-1} e^{-\int_{t_i}^{t_{i+1}} \mu(t) \; dt} \end{array}$

If we assume that $\mu(t)$ is constant on each interval, $[t_i,t_{i+1}]$, then each exponential can be written as

$\begin{array}{rcl} \displaystyle e^{-\int_{t_i}^{t_{i+1}} \mu(t) \; dt} & = & \displaystyle e^{\; -\int_{t_i}^{t_{i+1}} \mu_i \; dt} \\ & = & \displaystyle e^{\; -\mu_i \int_{t_i}^{t_{i+1}} dt} \\ & = & \displaystyle e^{\; -\mu_i \Delta t} \end{array}$

where $\mu_i$ is the constant $\mu(t)$ in the interval $[t_i,t_{i+1}]$.

Further simplifying: The Taylor Series of $e^{-\alpha}$ is $1 - \alpha + {\alpha^2 \over 2} + \cdots$, so we can approximate the exponential with the first two terms of its Taylor Series:

$ e^{\; -\mu_i \Delta t} \approx 1 - \mu_i \; \Delta t$

Substituting that approximation back into the product, we get the following very important formula:

$\displaystyle e^{-\int_0^D \mu(t) \; dt} \;\; \approx \;\; \prod_{i=0}^{n-1} \; (1 - \mu_i \Delta t)$

Note that $\mu_i \Delta t$ is the "fraction absorbed per unit length" times "length of the interval", which is the fraction of light absorbed over the interval. So $\mu_i \Delta t$ is the opacity of the interval, and $1 - \mu_i \Delta t$ is the transparency of the interval.

The algorithm to compute one pixel of the image, $I_\textrm{proj}$, is:

1. For the ray from the fluoro origin to the pixel, determine its entry point and exit point on the CT volume.
2. From the entry and exit points, determine $D$ and $n$.
3. Set $p = 1$.
4. Step through the volume:
   
   For $i$ = $0$ to $n-1$:
   1. Determine $\mu_i$ by looking up the attenuation in the CT volume at position $\ell( \frac12 (t_i + t_{i+1}) )$, which is in the middle of the interval. This requires a trilinear interpolation of the discretely sampled values in the CT volume.
   2. Set $p = p \times (1 - \mu_i \Delta t)$.
   
   
5. Then the pixel intensity is $I_\textrm{out} = p \; I_\textrm{in}$.

Note that this is only an approximation and can be improved by

• Taking smaller steps, $\Delta t$.
• Using more terms of the Taylor Series in approximating $e^{-\mu_i \Delta t}$.

This algorithm can be implemented in software, but it is now usually implemented in hardware on the Graphics Processing Unit (GPU).