A Discrete Formulation of Filtered Backprojection

Our first step is to replace

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

with

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

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

$$x' = x \cos{\theta} + y\sin{\theta}$$ (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

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

with

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

But the index k has only meaning for x'₀ 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:

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

where we have restricted k to run between -(K - 1) and (K - 1), because, for example, for x'_n = x₀ 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'₀. 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'₀, x'_K-1].

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

The discrete values of b_k are going to be:

$$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}$$ (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

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

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

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

2.

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