The Filtered Backprojection Method

  

Figure 2.1: X rays penetrating an object in a Computerised Tomography set up.

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:

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

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

$$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$$ (2.24)

The l=0 case of this formula is:

$$M(k, 0) = \int_{-\infty}^\infty \int_{-\infty}^\infty \mu(x,y) e^{-2\pi i k x} \,\textrm{d} x\,\textrm{d} y$$

$$\quad = \int_{-\infty}^\infty \left( \int_{-\infty}^\infty \mu(x,y) \,\textrm{d} y \right) e^{-2\pi i k x}\,\textrm{d} x$$

$$\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:

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

Then:

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

and

$$P(\theta, k') \int_{-\infty}^\infty p(\theta, x') e^{-2\pi i k' x'}\,\textrm{d} x'$$ =

$$\int_{-\infty}^\infty \int_{-\infty}^\infty \mu(x',y') e^{-2\pi i k' x'} \,\textrm{d} y' \,\textrm{d} x'$$ =

$$\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)$:

$$\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$$ (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.,

$$p(\theta, x') = p(\theta+\pi, L-x')$$ (2.29)

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

$$P(\theta, k') = P(\theta+\pi, -k')$$ (2.30)

and therefore:

 

$$\mu(x,y) \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)

$$\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)

$$\int_0^{\pi} C(\theta, x')\,\textrm{d}\theta$$ = (2.33)

where

$$C(\theta, x') = \int_{-\infty}^\infty P(\theta, k') e^{2\pi i k' x'} \vert k'\vert \,\textrm{d} k'$$ (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:

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

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

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

The function in the x' space that B(k') corresponds to is b(x'):

$$\int_{-\infty}^\infty B(k') e^{2\pi i k' x'} \,\textrm{d} k'$$ b(x') =

$$\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:

 

$$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'$$

$$\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:

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

In other words,

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