sLecture

↳ Topic 2: Tomographic Reconstruction

↳ Intro

↳ CT

↳ PET

↳ Co-ordinate Rotation

↳ Radon Transform

↳ Fourier Slice Theorem

↳ Convolution Back Projection

Excerpt from Prof. Bouman's Lecture

Accompanying Lecture Notes

introduction Convolution back projection (CBP) or forward back projection. do not need any calculation in frequency domain so it is computationally much faster.

definition? the back projection is like the adjoint or the transpose of the projection operator.

Summary of CBP algorithm

1. Measure projections $ p_{\theta}(r) $.
2. Filter the projections to obtain $ g_{\theta}(r) = h(r)*p_{\theta}(r) $.
3. Back project filtered projections using the following
   

$ f(x,y) = \int_0^{\pi}g_{\theta}(x\cos(\theta)+x\sin(\theta))d\theta $

describe with images the concept of back projection and smearing.

Back projection "smears" values of $ g(r) $ back over the image, and then adds smeared images for each angle. Then you add the back projections of (a) and (b). the theory says that if I take an infinite number of filtered back projections, I can get perfect reconstruction of the original image given that the original is band limited. But of course this is not possible in practice. In general, if you have an $ n $ by $ n $ image, you need to perform roughly $ n $ or more back projections.

MATLAB demo and images from the example at the end of the video clip on this page.

Derivation of the filter $ h(r) $

In order to compute the inverse CSFT of $ F(u,v) $ in polar coordinates, we must use the Jacobian of the polar coordinate transformation.
$ du dv = |\rho|d\theta d\rho \ $
where $ u = \rho\cos(\theta) $ and $ v = \rho\sin(\theta) $

This results in the expression
$ f(x,y) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)e^{2\pi j(xu+yv)}dudv $

Using the Fourier slice theorem, we can replace F(x,y) with its representation in polar coordinates, $ P_{\theta}(\rho) $
$ \begin{align} \Rightarrow f(x,y) &= \int_{-\infty}^{\infty}\int_0^{\pi}P_{\theta}(r)e^{2\pi j(x\rho\cos(\theta) +y\rho\sin(\theta))}dudv \\ &= \int_0^{\pi}\underbrace{[\int_{-\infty}^{\infty}|\rho|P_{\theta}e^{2\pi j\rho(x\cos(\theta) +y\sin(\theta))}d\rho]}_{g_{\theta}(x\cos(\theta) + y\sin(\theta))}d\theta \end{align} $

Then $ g_{\theta}(t) $ is given by
$ \begin{align} g_{\theta}(t) &= \int_{\infty}^{\infty}|\rho|P_{\theta}(\rho)e^{2\pi j\rho t}d\rho \\ &= CTFT^{-1}\{|\rho|P_{\theta}(\rho)\} \\ &= h(t) * p_{\theta}(r) \end{align} $

where $ h(t) = CTFT^{-1}\{|\rho|\} $, and
$ f(x,y) = \int_0^{\pi}g_{\theta}(x\cos(\theta) + y\sin(\theta))d\theta $

A closer look at Projection Filter

At each angle, projections are filtered to produce $ g){\theta}(r) $:
$ g_{\theta}(r) = h(r)*p_{\theta}(r) \ $.

The frequency response of the filter is given by
$ H(\rho) = |\rho| \ $
As you can see from figure 3, $ h(r) $ is like a high pass filter.

Fig 3: Frequency response of $ H(\rho) = |\rho| $

The frequency response depends on the absolute value of $ \rho $. Recall that $ \rho $ is the radial distance from the origin so this filter is essentially weighting the input by the magnitude of $ \rho $. This is because as you get further from the origin, the shaded area shown in figure 4 grows as $ \rho $ gets larger.

Fig 4: Area of a slice through circular rings.

The filter could theoretically have unlimited bandwidth but if the image is bandlimited, then the real filter does not need to have infinite bandwidth. The frequencies can be clipped at $ |\rho| $ ≤ $ f_c $ for some cutoff frequency $ f_c $.

The frequency response in figure 3 resembles a rect function minus a tri function. So
$ \begin{align} H(\rho) &= f_c[rect(f/(2f_c))-\wedge(f/f_c)] \\ h(r) &= f_c^2[2sinc(t2f_c)-sinc^2(tf_c)] \end{align} $

A Closer Look at Back Projection

Back projection function is given by
$ f(x,y) = \int_0^{\pi}b_{\theta}(x,y)d\theta $

where
$ b_{\theta}(x,y) = g_{\theta}(x\cos(\theta)+y\sin(\theta)) $

Consider the set of points $ (x,y) $ such that
$ r = x\cos(\theta)+y\sin(\theta) $

This set looks like

Fig 5: Graphical representation of the set $ \{(x,y):r=x\cos(\theta)+y\sin(\theta)\} $

Along this line, $ b_{\theta}(x,y) = g_{\theta}(r) $.

For each angle $ \theta $ back projection is constant along the lines of angle $ \theta $ and take on value $ g_{\theta}(r) $. As stated earlier, this is basically like the adjoint of the projection operator. In the forward projection, you sum together all the quantities along the solid line shown in figure 5. In the back projection,you take some value and proportionately put it in each location in proportion to the weight you used to add it up in the first place so they're like transpose operation. The domain of the forward projection operator becomes the range of the back projection operator and vice versa. you assign quantities along the line Complete back projection is formed by integrating (summing) back projections for angles ranging from $ 0 $ to $ \pi $.
$ \begin{align} f(x,y) &= \int_0^{\pi}b_{\theta}(x,y)d\theta \\ &\approx \frac{\pi}{M} \sum_{m=0}^{M-1} b_{\frac{m\pi}{M}}(x,y) \end{align} $

Fig 1: $ P_{\theta}(\rho) $ is $ F(u,v) $ in polar coordinates

References

• C. A. Bouman. ECE 637. Class Lecture. Digital Image Processing I. Faculty of Electrical Engineering, Purdue University. Spring 2013.

Questions and comments

If you have any questions, comments, etc. please post them on this page

Back to the "Bouman Lectures on Image Processing" by Maliha Hossain