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

Definition

Let
$ \begin{align} P_{\theta}(\rho) &= CTFT \{p_\theta(r)\} \\ F(u,v) &= CSFT\{f(x,y)\} \end{align} $

where $ \rho $ is the frequency variable corresponding to $ r $ just as $ u $ and $ v $ are the frequency variables corresponding to $ x $ and $ y $ respectively.

Then
$ P_{\theta}(\rho) = F(\rho\cos(\theta),\rho\sin(\theta)) \ $

Recall that $ p_{\theta}(r) $ is the projection of image $ f(x,y) $ at angle $ \theta $. $ P_{\theta}(\rho) $ is its 1-D Fourier transform. $ F(u,v) $ on the other hand, is the 2-D Fourier transform of image $ f(x,y) $.

So essentially, the theorem tells us that $ P_{\theta}(\rho) $ is $ F(u,v) $ in polar coordinates.

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

Another way to think about it is that the projection $ p_{\theta}(r) $ at angle $ \theta $ is equal to the slice of $ F(u,v) $ that passes through the origin at angle $ \theta $. Let us look at the simple example in figure 2 to illustrate this relationship. For the given $ f(x,y) $, the projection for $ \theta = 0 $° is a 1-D rect function. Let us assume it is of unit area.

Fig 2: $ p_{\theta}(r) $ at angle $ \theta = 0 $°

So we have that
$ \begin{align} p_{\theta}(r) |_{\theta = 0} &= rect(r) \\ \Rightarrow P_{\theta}(\rho) |_{\theta = 0} &= CTFT\{p_{\theta}(r) \}|_{\theta = 0} \\ &= sinc(\rho) \end{align} $

But $ f(x,y) $ is a 2-D rect so
$ \begin{align} f(x,y) &= rect(x)rect(y) \\ \Rightarrow F(u,v) &= CSFT\{f(x,y)\} \\ &= sinc(u)sinc(v) \end{align} $

Now lets convert $ F $ from cartesian coordinates $ u,v $ to polar coordinates $ (\rho,\theta) $ where $ \theta = 0 $°
$ \begin{align} F(u,v) &= sinc(u)sinc(v) \\ &= sinc(\rho \cos\theta)sinc(\rho \sin\theta)|_{\theta = 0} \\ &= sinc(\rho) \\ &= P_{\theta}(\rho)|_{\theta = 0} \end{align} $

So we see that the theorem holds for the above scenario. But it does not prove the theorem, only verifies it for a particular situation. The theorem can be proved in different ways. Two of the methods are presented below.

Proof

Method 1

The first method relies on plugging variable substitutions into the defintion.

By definition,
$ \begin{align} P_{\theta}(\rho) &= CTFT\{p_{\theta}(r)\} \\ &= \int_{-\infty}^{\infty}p_{\theta}(r)e^{-j2\pi\rho r}dr \\ &=\int_{-\infty}^{\infty} [\int_{-\infty}^{\infty}f(\mathbf{A_{\theta}} \begin{bmatrix} r \\ z \end{bmatrix}) dz]e^{-j2\pi\rho r}dr \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}f(\mathbf{A_{\theta}} \begin{bmatrix} r \\ z \end{bmatrix})e^{-j2\pi\rho r}dzdr \end{align} $

Next we make the following change of variables
$ \begin{bmatrix} r \\ z \end{bmatrix} = \mathbf{A_{-\theta}}\begin{bmatrix} r \\ z \end{bmatrix} $

Notice that the Jacobian is $ |{A_{-\theta}}| $$ =1 $ since
$ \begin{align} |\frac{\partial (r,z)}{\partial (x,y)}| &= det \begin{bmatrix} \frac{\partial r}{\partial x} & \frac{\partial r}{\partial y} \\ \frac{\partial z}{\partial x} & \frac{\partial z}{\partial y} \end{bmatrix} \\ &= det \begin{bmatrix} \frac{\partial (x\cos(\theta)+y\sin(\theta))}{\partial x} & \frac{\partial (x\cos(\theta)+y\sin(\theta))}{\partial y} \\ \frac{\partial (-x\sin(\theta)+y\cos(\theta))}{\partial x} & \frac{\partial (-x\sin(\theta)+y\cos(\theta))}{\partial y} \end{bmatrix} \\ &= det \begin{bmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{bmatrix} \\ &= \cos^2\theta + \sin^2\theta \\ &=1 \end{align} $

Then,
$ drdz = |\frac{\partial(r,z)}{\partial(x,y)}|dxdy = dxdy $

Also notice that $ r=x\cos(\theta)+y\sin(\theta) $. So finally we have that
$ \begin{align} P_{\theta}(\rho) &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)e^{-j2\pi\rho[x\cos(\theta)+y\sin(\theta)]}dxdy \\ &= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)e^{-j2\pi[x\rho\cos(\theta)+y\rho\sin(\theta)]}dxdy \\ &= F(\rho\cos(\theta),\rho\sin(\theta))_{\blacksquare} \end{align} $

Method 2

The second method is slightly more compact than the first.

First, let $ \theta = 0 $, then
$ p_0(r) = \int_{-\infty}^{\infty}f(r,y)dy $
Note that when $ \theta = 0 $, the $ y $ and $ z $ axes are the same, as are the $ x $ and $ r $ axes, and in the above integral, $ r $ and $ y $ are just dummy variables. Next, taking the CTFT of both sides, we get
$ \begin{align} \Rightarrow P_0(\rho) &= \int_{-\infty}^{\infty}p_0(r)e^{-2\pi jr\rho}dr \\ &= \int_{-\infty}^{\infty}[\int_{-\infty}^{\infty}f(r,y)dy]e^{-2\pi jr\rho}dr \\ &= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(r,y)e^{(-2\pi j(r\rho+y0)}drdy \\ &= F(\rho,0) \end{align} $

Recall the rotation property of the CSFT; if $ CSFT\{f(x,y)\} = F(u,v) $, then
$ CSFT\{f(\mathbf{A}\begin{bmatrix}x \\ y\end{bmatrix})\} = |\mathbf{A}|^{-1}F([u,v]\mathbf{A}^{-1}) $

In polar coordinates, this becomes
$ CSFT\{f(r,\theta+\alpha)\} = F(\rho,\phi+\alpha) \ $

So we see that if it holds under $ \theta = 0 $, it must hold for any $ \theta $ since the CSFT is invariant under rotation.

Drawbacks of the Fourier Slice Theorem

The Fourier slice theorem is beautiful because it allows you to reconstruct $ f(x,y) $ but unfortunately, there are many technical challenges behind implementing this in reality.

Fig 3: Inverse Radon Transform

The physical system measures $ p_{\theta}(r) $. From these, we can calculate $ P_{\theta}(\rho) = CTFT\{p_{\theta}(r)\} $, which is $ F(u,v) $ in polar coordinates. So we can plot $ F(u,v) $ from projections but using this method, data points for $ F(u,v) $ are available only along the radial lines shown figure 3. To implement the inverse Fourier transform to obtain $ f(x,y) $, we need data points on a rectangular grid. This grid does not match the radial grid. We could interpolate data points but this is computationally expensive, and far from the origin, the data is sparse so it is harder to interpolate. So while the Fourier slice theorem illustrates a simple and beautiful relationship between the image and its projections, we cannot put it to use in practical implementation. Instead, convolution back projection is the most commonly used method to recover the image and this will be the topic of discussion in the next section.

References

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

• "Table of Continuous-space (CS) Fourier Transform Pairs and Properties." From ProjectRhea's. Collective Table of Formulas. https://www.projectrhea.org/rhea/index.php/Continuous_Space_Fourier_Transform_(frequences_in_hertz). April 23rd, 2013. [May 27th, 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