next up previous index
Next: FFTPACK Up: Computerised Tomography Previous: The Filtered Backprojection Method

A Discrete Formulation of Filtered Backprojection

Our first step is to replace

\begin{displaymath}\mu(x,y) = \int_0^\pi C(\theta, x')\,\textrm{d}\theta
\end{displaymath}

with

\begin{displaymath}\mu(x,y) \approx \frac{\pi}{N}\sum_{m=0}^{N-1}C\left(\theta_m, x'\right)
\end{displaymath} (2.38)

where a single pair $(\theta_m, x'_n)$ yields such (x,y) that

\begin{displaymath}x' = x \cos{\theta} + y\sin{\theta}
\end{displaymath} (2.39)

so that going through all measured values of $(\theta_m, x'_n)$ will eventually cover the whole plane $x\times y$.

Some comments:

1.
Points for small values of k' are much denser than points for large values of k', hence the importance of the filtering function. Filtering $M(\theta, k')$ with $\vert k'\vert$ weights the data to compensate for the sparser fill of the frequency domain at larger $\vert k'\vert$.
2.
Filter B(k') does the same job, whereas its relative b(x') is equivalent, but operates in the spatial domain.

The next step is to discretise

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

with

\begin{displaymath}C(\theta_m, x'_n) = \sum_{k = -\infty}^\infty p(\theta_m, x'_k)
b(x'_n - x'_k) \Delta x'
\end{displaymath} (2.40)

But the index k has only meaning for x'0 through x'K-1, and the positions for which we have any measurements of $p(\theta, x')$ are $x'_k = x'_0 + k\Delta x'$.

Convolution is a symmetric operation, i.e., $p \star b = b \star p$. This can be proven simply by subsitution. This means that we can move x'n - x'k to p getting:

\begin{displaymath}C(\theta_m, x'_n) = \Delta x'
\sum_{k=-(K-1)}^{K-1}p(\theta_m, x'_n - x'_k) b(x'_k)
\end{displaymath} (2.41)

where we have restricted k to run between -(K - 1) and (K - 1), because, for example, for x'n = x0 we can subtract a negative x'-(K-1) and still stay within the measured region, and for x'n = x'K-1 we can subtract x'K-1 and go back to x'0. But we are still going to hit some combinations of x'n and x'k that are going to push us beyond the physically meaningful region [x'0, x'K-1].

What then? We are simply going to pad our vector p with zeros.

The discrete values of bk are going to be:

\begin{displaymath}b_k = b(x'_k) = \Bigg\{
\begin{array}{ll}
1/(2\Delta x')^2,...
...
-1/(k\pi\Delta x')^2 &\qquad \hbox{for odd}\> k
\end{array}\end{displaymath} (2.42)

In summary:

1.
to compute $C(\theta_m, x'_n)$ for a given nwe have to evaluate a scalar product of p and b, both of length 2K-1
2.
to compute $C(\theta_m, x'_n)$ for all n will have to perform $K \times (2K-1) \approx {\cal O}(K^2)$floating point multiplications.

This is going to be expensive.

Because of the low cost of FFT it is cheaper to

1.
Take FFT of a padded data function $p = \left( p_{-(K-1)}, \ldots,
p_{K-1}\right)$ to form $P(\theta_m, k')$
2.
Discretise

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

This algorithm is going to be cheaper, because its major cost is associated with FFTs and these will run like $(2K-1){\cal O}(2K - 1)$.

What next?

1.
For a given angle $\theta_m$ the vector $C(\theta_m, x'_n)$ gives values of C at points (x,y) such that

\begin{displaymath}x_{m,n}\cos{\theta_m} + y_{m,n}\sin{\theta_m} = x'_n
\end{displaymath} (2.43)

2.
The filtered projection values are mapped onto their closest grid box (xg, yg), and then summed within that box to produce density $\mu(x_g, y_g)$.


next up previous index
Next: FFTPACK Up: Computerised Tomography Previous: The Filtered Backprojection Method
Zdzislaw Meglicki
2001-02-26