next up previous index
Next: A Discrete Formulation of Up: Computerised Tomography Previous: Computerised Tomography

The Filtered Backprojection Method


  
Figure 2.1: X rays penetrating an object in a Computerised Tomography set up.
\begin{figure}
\par\begin{center}
\par\begin{picture}
(8,8)(0,0)
\par % \put(2,6...
...,0)(0.5, 0.0){17}{\vector(0,1){8}}
\end{picture}\par\end{center}\par\end{figure}

Assume a computerised tomography set up with rays parallel to the y-axis directed upwards, see Figure 2.1. Assume absorbing density in the plane to be $\mu(x,y)$, where function $\mu$ has a compact support. Signal intensity registered on the x axis, $p(\theta=0, x)$ will then be:

\begin{displaymath}p(\theta=0, x) = \int_{-\infty}^{\infty}\mu(x,y)\,\textrm{d}y
\end{displaymath} (2.23)

Now consider a two-dimensional Fourier transform of function $\mu(x,y)$:

\begin{displaymath}M(k, l) = \int_{-\infty}^\infty \int_{-\infty}^{\infty} \mu(x,y)
e^{-2\pi i (k x + l y)} \,\textrm{d} x\,\textrm{d} y
\end{displaymath} (2.24)

The l=0 case of this formula is:
    $\displaystyle M(k, 0) = \int_{-\infty}^\infty \int_{-\infty}^\infty
\mu(x,y) e^{-2\pi i k x} \,\textrm{d} x\,\textrm{d} y$  
    $\displaystyle \quad = \int_{-\infty}^\infty
\left(
\int_{-\infty}^\infty \mu(x,y) \,\textrm{d} y
\right)
e^{-2\pi i k x}\,\textrm{d} x$  
    $\displaystyle \quad = \int_{-\infty}^\infty
p(\theta=0, x)
e^{-2\pi i k x}\,\textrm{d} x$ (2.25)

In other words:
Taking the one-dimensional Fourier Transform of the projection $p(\theta=0, x)$ along the y axis gives us the (k,0) line in the Fourier space of $\mu(x,y)$.

We can now rotate our set-up and obtain lines in the Fourier space of $\mu(x,y)$ under any angle by calculating one-dimensional Fourier transforms of measured parallel projections, and in this way effectively reconstructing the whole M(k, l). Once we have M(k, l) we can then obtain $\mu(x,y)$ by taking the inverse Fourier Transform of M(k, l).

It is convenient in this case to write it all down in terms of the rotation angle $\theta$.

Consider rotating the system of coordinates as follows:

\begin{displaymath}\left(
\begin{array}{rr}
\cos{\theta} & \sin{\theta} \\
-...
...t)
=
\left(
\begin{array}{c}
x' \\
y'
\end{array}\right)
\end{displaymath} (2.26)

Then:

\begin{displaymath}p(\theta, x') = \int_{-\infty}^\infty \mu(x',y')\,\textrm{d} y'
\end{displaymath} (2.27)

and
$\displaystyle P(\theta, k')$ = $\displaystyle \int_{-\infty}^\infty p(\theta, x')
e^{-2\pi i k' x'}\,\textrm{d} x'$  
  = $\displaystyle \int_{-\infty}^\infty \int_{-\infty}^\infty
\mu(x',y') e^{-2\pi i k' x'} \,\textrm{d} y'
\,\textrm{d} x'$  
  = $\displaystyle \int_{-\infty}^\infty \int_{-\infty}^\infty
\mu(x',y') e^{-2\pi i k' (x\cos{\theta} + y\sin{\theta})}
\,\textrm{d} x \,\textrm{d} y = M(\theta, k')$  

The inverse transform of $P(\theta, k')$ yields $\mu(x,y)$:

\begin{displaymath}\mu(x,y) = \int_0^{2\pi}\int_0^\infty P(\theta, k')
e^{2\pi ...
...\theta}+ y\sin{\theta})} k'\,\textrm{d} k'
\,\textrm{d}\theta
\end{displaymath} (2.28)

A closer inspection shows that it is not necessary to rotate the apparatus through the whole $[0, 2\pi]$, because swapping the source and the detector does not change the resulting attenuation, i.e.,

\begin{displaymath}p(\theta, x') = p(\theta+\pi, L-x')
\end{displaymath} (2.29)

where L is the length of the detector. This, in turn, implies that

\begin{displaymath}P(\theta, k') = P(\theta+\pi, -k')
\end{displaymath} (2.30)

and therefore:
 
$\displaystyle \mu(x,y)$ = $\displaystyle \int_0^{\pi}\int_{-\infty}^\infty P(\theta, k')
e^{2\pi i k' (x\cos{\theta}+ y\sin{\theta})} \vert k'\vert
\,\textrm{d} k'\,\textrm{d}\theta$ (2.31)
  = $\displaystyle \int_0^{\pi}\int_{-\infty}^\infty P(\theta, k')
e^{2\pi i k' x'} \vert k'\vert
\,\textrm{d} k'\,\textrm{d}\theta$ (2.32)
  = $\displaystyle \int_0^{\pi}
C(\theta, x')\,\textrm{d}\theta$ (2.33)

where

\begin{displaymath}C(\theta, x') = \int_{-\infty}^\infty P(\theta, k')
e^{2\pi i k' x'} \vert k'\vert
\,\textrm{d} k'
\end{displaymath} (2.34)

is called a filtered projection, and $\vert k'\vert$ plays the role of an inverse filter: multiplying $P(\theta, k')$ by $\vert k'\vert$ increases the influence of $P(\theta, k')$ at high frequencies.
Integral (2.36) defines the map of density $\mu(x,y)$in the $x\times y$ plane by accumulating all of the filtered projections for all angles $\theta$ from 0 to $\pi$. Each filtered projection $C(\theta, x')$ contributes to the density along the line of constant $x' = x\cos{\theta} + y\sin{\theta}$for a particular value of $\theta$. And so each filtered projection is backprojected into the $x\times y$ plane.

In order to avoid aliasing problems associated with the Nyquist critical frequency $k' = 1/(2\Delta x')$, we shall introduce a new filter function defined as follows:

\begin{displaymath}B(k') = \Bigg\{
\begin{array}{ll}
\vert k'\vert, & \qquad ...
...eq 1/(2\Delta x')\\
0, & \qquad \hbox{otherwise}
\end{array}\end{displaymath} (2.35)

and convolve it with $P(\theta, k')$ in the expression that defines the filtered projection:

\begin{displaymath}C(\theta, x') = \int_{-\infty}^\infty P(\theta, k') B(k')
e^{2\pi i k' x'}
\,\textrm{d} k'
\end{displaymath} (2.36)

The function in the x' space that B(k') corresponds to is b(x'):
b(x') = $\displaystyle \int_{-\infty}^\infty B(k') e^{2\pi i k' x'} \,\textrm{d} k'$  
  = $\displaystyle \frac{1}{2(\Delta x')^2}
\frac{\sin\left(2\pi x' / (2\Delta x')\r...
...t(
\frac{\sin \left(\pi x'/(2\Delta x')\right)}
{\pi x'/(2\Delta x')}
\right)^2$  

Now, remember that a Fourier transform of a convolution is equal to a product of Fourier transforms. Since $P(\theta, k')$ and B(k') are Fourier transforms of $p(\theta, x')$ and b(x'), their product is equal to Fourier transform of a convolution of the latter two functions:
 
    $\displaystyle P(\theta, k') B(k') = \int_{-\infty}^\infty p(\theta, x') e^{-2\p...
...xtrm{d} x'
\times
\int_{-\infty}^\infty b(x') e^{-2\pi i k' x'}
\,\textrm{d} x'$  
    $\displaystyle \quad = \int_{-\infty}^\infty
\left(
\int_{-\infty}^\infty p(\theta, x) b(x' - x) \,\textrm{d} x
\right)
e^{-2\pi i k' x'} \,\textrm{d} x'$  

Since $C(\theta, x')$ is the inverse Fourier Transform of $P(\theta, k') B(k')$, the Fourier Transform in the equation (2.42) gets undone and we're simply left with:

\begin{displaymath}C(\theta, x') = \int_{-\infty}^\infty p(\theta, x) b(x' - x) \,\textrm{d} x
\end{displaymath} (2.37)

In other words,
the filtered projection $C(\theta, x')$ is the convolution of $p(\theta, x')$ and b(x').


next up previous index
Next: A Discrete Formulation of Up: Computerised Tomography Previous: Computerised Tomography
Zdzislaw Meglicki
2001-02-26