Line 67: | Line 67: | ||
\end{align} | \end{align} | ||
</math><br /> | </math><br /> | ||
− | Pulling out <math>g_{\theta}(t)</math>,<br /> | + | Pulling out <math>g_{\theta}(t) \ </math>,<br /> |
<math>\begin{align} | <math>\begin{align} | ||
− | g_{\theta}(t) &= \int_{\infty}^{\infty}|\rho | P_{\theta}(\rho )e^{2\pi j\rho t} d\rho \\ | + | g_{\theta}(t) |
+ | &= \int_{\infty}^{\infty}|\rho | P_{\theta}(\rho )e^{2\pi j\rho t} d\rho \\ | ||
&= CTFT^{-1}\{|\rho |P_{\theta}(\rho)\} ", by comparison with the inverse CTFT"\\ | &= CTFT^{-1}\{|\rho |P_{\theta}(\rho)\} ", by comparison with the inverse CTFT"\\ | ||
&= h(t) * p_{\theta}(r) \\ | &= h(t) * p_{\theta}(r) \\ | ||
Line 76: | Line 77: | ||
=Projection Filter Analysis= | =Projection Filter Analysis= | ||
Now let's focus on <math>g_{\theta}(r)</math>.<br /> | Now let's focus on <math>g_{\theta}(r)</math>.<br /> | ||
− | <math>g_{\theta}(r) = h(r) * p_{\theta}(r)</math><br /> | + | <math>g_{\theta}(r) = h(r) * p_{\theta}(r) \ </math><br /> |
− | The frequency response of the filter is | + | The frequency response of the filter is<br /> |
− | <math>H(\rho) = CTFT{h(r)} = |\rho|</math><br /> | + | <math>H(\rho) = CTFT{(h(r))} = |\rho| \ </math><br /> |
After graphing the frequency response, it is apparent that the filter is a high pass filter. <br /> | After graphing the frequency response, it is apparent that the filter is a high pass filter. <br /> | ||
Revision as of 13:17, 21 December 2014
Convolution/Fourier Back Projection Algorithm
A slecture by ECE student Sahil Sanghani
Partly based on the ECE 637 material of Professor Bouman.
Contents
Introduction
Convolution Back Projection (CBP) offers a reconstruction method that is not computationally expensive. Although the method is based on the Fourier Slice Theorem, there's never a transformation to the frequency domain. Theoretically CBP involves extruding every projection back through the origin and then summing the results. This operation is called back projection. The projections in Figure 1 were all assumed to be the same regardless of $ \theta $. In Figure 1a, the extrusion is demonstrated. In Figure 1b, the summing is demonstrated.
Summary
There are 3 steps to reconstructing an object from its projections using CBP.
1. Measure the projections $ p_{\theta}(r) $. (i.e. using CT)
2. Filter the projections $ g_{\theta}(r) = h(r) * p_{\theta}(r) $
3. Back project filtered projections
$ f(x,y) = \int_0^{\pi}g_{\theta}(x\cos(\theta)+x\sin(\theta))d\theta $
An infinite number of filtered back projections will result in a perfect reconstruction of the original image of the object given it is band limited. Since it is impossible to take that many projections, a good practice is to take at least $ n $ back projections for a $ n $ by $ n $ image.
Derivation
From the Fourier Slice Theorem, we get the following relationships.
$ \begin{align} u &= \rho\cos(\theta) \\ v &= \rho\sin(\theta) \end{align} $
Now let's calculate the Jacobian of the polar coordinate transformation.
$ dudv = \left | \frac{\partial (u,v)}{\partial (\theta,\rho)} \right \vert d\theta d\rho $
$ \begin{align} \left | \frac{\partial (u,v)}{\partial (\theta,\rho)} \right \vert &= det\begin{bmatrix} \frac{\partial u}{\partial \theta} & \frac{\partial u}{\partial \rho} \\ \frac{\partial v}{\partial \theta} & \frac{\partial u}{\partial \rho} \end{bmatrix} \\ &= det\begin{bmatrix} \frac{\partial (\rho\cos(\theta))}{\partial \theta} & \frac{\partial (\rho\cos(\theta))}{\partial \rho} \\ \frac{\partial (\rho\sin(\theta))}{\partial \theta} & \frac{\partial (\rho\sin(\theta))}{\partial \rho} \end{bmatrix} \\ &= det\begin{bmatrix} -\rho\sin(\theta) & \cos(\theta) \\ \rho\cos(\theta) & \sin(\theta) \end{bmatrix} \\ &= |-\rho\sin^{2}(\theta) - \rho\cos^{2}(\theta)| \\ &= |-\rho(\sin^{2}(\theta) + \cos^{2}(\theta))| \\ &= |-\rho| \\ &= |\rho| \\ \end{align} $
$ \Rightarrow dudv = |\rho|d\theta d\rho $
Starting with the formula for the inverse CSFT and using the Fourier Slice theorem, we will end up with
$ \begin{align} f(x,y) &= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)e^{2\pi j(xu+yv)}dudv \\ &= \int_{-\infty}^{\infty}\int_{0}^{\pi}F(\rho\cos(\theta),\rho\sin(\theta))e^{2\pi j(x\rho\cos(\theta) +y\rho\sin(\theta))}|\rho|d\theta d\rho \\ &= \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} $
Pulling out $ g_{\theta}(t) \ $,
$ \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)\} ", by comparison with the inverse CTFT"\\ &= h(t) * p_{\theta}(r) \\ \end{align} $
Projection Filter Analysis
Now let's focus on $ g_{\theta}(r) $.
$ g_{\theta}(r) = h(r) * p_{\theta}(r) \ $
The frequency response of the filter is
$ H(\rho) = CTFT{(h(r))} = |\rho| \ $
After graphing the frequency response, it is apparent that the filter is a high pass filter.
This filter can be represented by a rect function minus a triangle function, so its equation is
$ \begin{align} H(\rho) &= f_c[rect(f/(2f_c))-\wedge(f/f_c)] \\ h(r) &= f_c^2[2sinc(2f_c t)-sinc^2(tf_c)] \end{align} $
Back Projection Analysis
References:
[1] C. A. Bouman. ECE 637. Class Lecture. Digital Image Processing I. Faculty of Electrical Engineering, Purdue University. Spring 2013.