sLecture

Topic 2: Tomographic Reconstruction
Intro
CT
PET
Co-ordinate Rotation
Radon Transform
Fourier Slice Theorem
Convolution Back Projection


The Bouman Lectures on Image Processing

A Slecture by Maliha Hossain

Subtopic 6: Fourier Slice Theorem

© 2013




Excerpt from Prof. Bouman's Lecture

Click here to view video directly on Youtube


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} x \\ y \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.



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

Alumni Liaison

Basic linear algebra uncovers and clarifies very important geometry and algebra.

Dr. Paul Garrett