up

PDF (letter size)

Contents

1 Introduction

We are familiar with the standard x-ray imaging, the type one encounters at a doctor or a dentist office. Briefly, this technique of medical imaging works as follows: A source is used to emit x-rays which traverse through the body and is then detected and recorded as an image on the detector surface located behind the body.

As x-rays travel throughout the body, its intensity attenuates. Each x-ray beam will attenuate differently, and will do so in proportion to the type and amount of tissue it passes though. This attenuation occurs due to the fact that different types of tissue (for example, bone vs. muscle vs. soft tissue) absorbs different amount of radiation.

The recording of the varying intensities of the x-rays leaving the body (the x-ray flux hitting the detector surface) gives a 2D image which reflects the content of the section of the body that the x-rays passed through.¹

Due to overlapping and shadowing between the internal tissues in the body, the final image recorded will not give a clear picture of the section of the body being examined.

CT uses x-rays as well, however in parallel x-ray CT, thin x-ray parallel beams are transmitted across a section of the body at a specific angle \(\theta \). When the beams hits the detector on the other side of the body, the flux is recorded which represents the projection of the cross section at the angle \(\theta .\) The angle is incremented to \(\theta +\delta \theta \) and another projection is obtained of the same body section.

By repeating this process we will obtain a sequence of projections. This sequence of projections is used to reconstruct a 2D image of that section of the body. This report will show the mathematical derivation of the reconstructed image using Filtered Back projection method (FBP) and the central role played by the spatial Fourier transform in this process. The following diagram illustrates the above discussion.

Therefore, we see that there are two main phases in CT. The first phase is one in which a sequence of projections are generated, and the second phase is one in which these images are combined to reconstruct an approximation to the original 2D image.

The operation which accomplishes the first phase is mathematically called the radon transform. The operation which accomplishes the second phase is the filtered backprojection. We start by reviewing the first phase, showing how the projection is first generated. The equation of the projection is then used in the second phase.

2 Illustration of the problem of image reconstruction by use of projections

Before going further into the discussion, it is useful to spend a little more time to illustrate what we mean by a projection, and to make sure we understand the problem we are about to solve.  To do this, we will use an example of a simple 2D image of 4 pixels of some random values as follows

The above 2D image represent a simplification of the cross section of the body which the xrays in CT scan would go through.

Now, let us obtain 2 projections (we can obtain more projections if we want to, but for now, we will assume we have only 2 projections).

We point the xray source at 2 angles and generate the projections as follows

Notice that the projection is an integral operation along the path of the ray. In other words, it sums the pixel values it encounters. Now that we have obtained 2 projections, the problem is how to determine the original image from the knowledge of just these projections and the angles they were obtained at

We see that this is an inverse problem with many solutions. In other words, one can come up with more than one image whose projections are those given.

There are number of approaches to solve this problem, as illustrated by the following diagram

2.1 Solving the problem algebraically

We can try to algebraically solve the above problem by setting 4 equations as follows. First let the image pixels be variables as follows

Hence we have the following equations

\begin {align*} A+B & =5\\ C+D & =8\\ C+A & =6\\ D+B & =7 \end {align*}

Or in matrix form\[\begin {pmatrix} 1 & 1 & 0 & 0\\ 0 & 0 & 1 & 1\\ 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 1 \end {pmatrix}\begin {pmatrix} A\\ B\\ C\\ D \end {pmatrix} =\begin {pmatrix} 5\\ 8\\ 6\\ 7 \end {pmatrix} \] The above is in the form \(Qx=b\).  The \(Q\) matrix above is not invertible since its determinant is zero. Hence there does not exist one unique solution to the problem. These are 3 possible solutions among infinity many solutions

Hence the best we can do to obtain a solution using the above approach, is to find the solution with the least square error. The more projections we obtain, the less this error will becomes.

Solving this inverse problem algebraically is not done in practice due to the large number of equations.

2.2 Solving the problem using the Fourier central Slice theorem

The following diagram illustrates the main idea behind this approach

2.3 Solving the problem using filtered backprojection

The following diagram illustrates the main idea behind this approach. In filtered backprojection, the projection spectrum is first filtered, then backprojected directly into a 2D space, and additional projections and backprojections are accumulated on top of this, which results in the final 2D reconstruction. This approach is the one used in practice in place of the direct approach of using the inverse 2D FFT as shown in the above diagram since it leads to a better quality image.

Now we begin a more mathematical analysis of the above solution, showing the derivation of the filtered backprojection operation.  

We start with discussion of the radon transform, which is the mathematical basis of the projection operation.

3 Radon transform introduction

Consider a 2D object in the xy reference frame. Let L be an x-ray beam going through this object at some angle \(\theta \) as follows

The equation of the line \(L\) can be written in 2 ways. The standard way is \[ y=mx+b \] where \(m\) is the slope and \(b\) is the intercept. It can also be written in terms of the parameters \(p\) and \(\theta \) as \[ L\left ( p,\theta \right ) =\left \{ \left ( x,y\right ) :x\cos \theta +y\sin \theta =p\right \} \] Any point \(\left ( x,y\right ) \) on the line \(L\) with specific \(p\) and specific \(\theta \) satisfies \(x\cos \theta +y\sin \theta =p\)

Assume now there exist a function \(f\left ( x,y\right ) \) defined over the region shown above. The integral of this function over the line \(L\left ( p,\theta \right ) \) is\[ I=\int _{L}f\left ( x,y\right ) ds \] where \(ds\) is a differential element of the line

It would be simpler to express the above integral in terms of \(x\) and \(y\). To do that, a trick is used with the help of the delta function. The above integral can be written as\[ I={\displaystyle \int \limits _{y=-\infty }^{\infty }} {\displaystyle \int \limits _{x=-\infty }^{\infty }} f\left ( x,y\right ) \ \delta \left ( x\cos \theta +y\sin \theta -p\right ) \ dxdy \] Hence for a specific \(p,\theta \) the above will integrate \(f\left ( x,y\right ) \) over the line \(L\left ( p,\theta \right ) \). The above is the radon transform of \(f\left ( x,y\right ) \) over the line \(L\left ( p,\theta \right ) \).

The result of the above radon transform is one numerical value. It is the line integral value. We can imagine a projection line into which we accumulate the result of these line integrals as follows

Now suppose we have many parallel lines and we perform the same line integral of \(f\left ( x,y\right ) \) over each one of these lines (since all these lines are parallel to line \(L\), then all of them will have the same \(\theta \), but each they will have different \(p.\) This will result in many radon transform integrals as  shown in the following diagram

The projection shown above is a discrete function. It is a function of \(p\) and \(\theta \). For the purpose of the discussion that follows, it is assumed that the projections generated by CT are continuous function (each projection is the flux of the x-ray on the plate) and is then written as\[ g\left ( p,\theta \right ) ={\displaystyle \int \limits _{y=-\infty }^{\infty }} {\displaystyle \int \limits _{x=-\infty }^{\infty }} f\left ( x,y\right ) \ \delta \left ( x\cos \theta +y\sin \theta -p\right ) \ dxdy \]

4 Image reconstruction

The input to CT 2D reconstruction is the continuous function \(g\left ( p,\theta \right ) .\)Before we go any further, we need to make sure we keep track of the angle \(\theta \) being used in the above expression. Recall that this is the angle that each projection is generated at. This angle goes from \(0\) to \(\pi .\) (Since going beyond that will duplicate measurements). Let the number of projections to be generated be \(M\). Hence \[ \Delta \theta =\frac {\pi }{M}\] Hence \(g\left ( p,\theta _{m}\right ) \) means the projection at angle \(\theta _{m}\), which is the angle \(m\frac {\pi }{M}\). For example, \(\theta _{0}=0^{0}\), \(\theta _{1}=\frac {\pi }{M}\) and \(\theta _{M}=M\frac {\pi }{M}=\pi \) and so on.

The next step is to sample \(g\left ( p,\theta _{m}\right ) \) and then to obtain the DFT of the sampled sequence. To sample \(g\left ( p,\theta _{m}\right ) \), we assume that the smallest distance between 2 adjacent repeating cycles of intensities is \(\tau \) cm (or in millimeter).  What this means is that largest spatial frequency present in the projection data is \(W=\frac {1}{\tau }\). Therefore, the projection needs to be sampled at frequency no less than \(2W\) (Nyquist sampling theorem). Hence we will us as the sampling frequency\[ f_{s}=2W \] The above means that \(g\left ( p,\theta _{m}\right ) \) is sampled with an interval of width \(\frac {\tau }{2}\) cm. Let the number of samples obtained from the above sampling frequency be \(N\) (and we assume it to be even as we can always reduce the sampling interval to obtain the next even value of \(N\))

The result of sampling \(g\left ( p,\theta _{m}\right ) \) will be the sequence of numbers \(g\left ( \frac {n}{2W},\theta _{m}\right ) \) where \[ n=-\frac {N}{2},\cdots ,0,1,\cdots ,\frac {N}{2}-1 \] There will be different such sequence for different angle \(\theta _{m}.\)

The following diagram helps illustrate the above formulation

Once a projection is sampled, its DFT is obtained as follows\begin {equation} S\left ( k,\theta _{m}\right ) ={\displaystyle \sum \limits _{n=-\frac {N}{2}}^{\frac {N}{2}-1}} g\left ( \frac {n}{2W},\theta _{m}\right ) e^{-j\frac {2\pi }{N}kn}\ \ \ \ \ k=0,1,\cdots ,N-1 \tag {1} \end {equation} Since the sampling frequency is \(f_{s}=2W\), hence the frequency axis will extend from \(-W\) to \(W\) and \(\Delta f=\frac {2W}{N}\)as shown in the following diagram

Let us now review what we have done so far. We sampled a projection at some angle \(\theta _{m}\) and obtained the discrete Fourier transform of the sampled projection. Now to make more progress, we must resort to the Fourier central slice theorem. This theorem tells us that the Fourier transform of a projection taken at angle \(\theta _{m}\) is equal to the values found along a slice in the 2D Fourier transform of the original image itself, as long as this line goes through the origin of the 2D Fourier transform plane and has the same angle \(\theta _{m}\). The following diagram illustrates the central slice theorem

Let the 2D Fourier transform of \(f\left ( x,y\right ) \) be \(F\left ( u,v\right ) \), hence from the definition of the inverse 2D Fourier transform we write\begin {equation} h\left ( x,y\right ) ={\displaystyle \int \limits _{-\infty }^{\infty }} {\displaystyle \int \limits _{-\infty }^{\infty }} F\left ( u,v\right ) e^{j2\pi xu}e^{j2\pi yv}dudv \tag {2} \end {equation} Where \(h\left ( x,y\right ) \) is an approximation to the original image \(f\left ( x,y\right ) \).

\(h\left ( x,y\right ) \) is what we call the filtered backprojection image. We now need to rewrite (2) in polar coordinates (to allow us to later substitute the Fourier transform of the projection that we obtained earlier into the above double integral).  Using the following diagram, we obtain the substitution needed

In the above diagram, \(f\) is radial distance from the original (the DC frequency) and \(\theta \) is the angle measured counter clock wise from the \(u\) axis.

Hence \(u\left ( f,\theta \right ) =f\cos \theta \) and \(v\left ( f,\theta \right ) =f\sin \theta \), and the Jacobian of transformation is the determinant of \(J=\) \(\left ( \begin {array} [c]{cc}\frac {\partial u}{\partial f} & \frac {\partial u}{\partial \theta }\\ \frac {\partial v}{\partial f} & \frac {\partial v}{\partial \theta }\end {array} \right ) \), hence \begin {align*} \left \vert J\right \vert & =\left \vert \begin {array} [c]{cc}\frac {\partial u}{\partial f} & \frac {\partial u}{\partial \theta }\\ \frac {\partial v}{\partial f} & \frac {\partial v}{\partial \theta }\end {array} \right \vert =\left \vert \begin {array} [c]{cc}\cos \theta & -f\sin \theta \\ \sin \theta & f\cos \theta \end {array} \right \vert \\ & =f\cos ^{2}\theta +f\sin ^{2}\theta \\ & =f \end {align*}

Hence the integration differentials can now be changed as \begin {align} du\ dv & =\left \vert J\right \vert dfd\theta \nonumber \\ & =f\ dfd\theta \tag {3} \end {align}

Using (3) in (2), we obtain, after changing the limit of integration\begin {align} h\left ( x,y\right ) & ={\displaystyle \int \limits _{\theta =0}^{2\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta \right ) e^{j2\pi xf\cos \theta }e^{j2\pi yf\sin \theta }\ f\ dfd\theta \nonumber \\ & ={\displaystyle \int \limits _{\theta =0}^{2\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\ f\ dfd\theta \tag {4} \end {align}

The above is the 2D inverse Fourier transform in polar coordinates.

Split the outside integral in (4) which goes from \(0\) to \(2\pi \) into 2 halves as follows\begin {align} h\left ( x,y\right ) & ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\ f\ dfd\theta \nonumber \\ & +{\displaystyle \int \limits _{\theta =\pi }^{2\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\ f\ dfd\theta \tag {5} \end {align}

The above 2 integrals can be combined into one integral using the property of the symmetry of the Fourier Transform \(F\left ( f,\theta \right ) \). To do that, we manipulate the second term in the above as follows. First, let the second term be called \(\Delta \) for now, i.e. \[ \Delta ={\displaystyle \int \limits _{\theta =\pi }^{2\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\ f\ dfd\theta \] Then\[ \Delta ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta +\pi \right ) e^{j2\pi f(x\cos \left ( \theta +\pi \right ) +y\sin \left ( \theta +\pi \right ) )}\ f\ dfd\theta \] and since the object \(f\left ( x,y\right ) \) is real valued, the 2D Fourier transform \(F\left ( f,\theta \right ) \) is symmetric. Hence \(F\left ( f,\theta +\pi \right ) =F\left ( -f,\theta \right ) \), then the above becomes\[ \Delta ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( -f,\theta \right ) e^{j2\pi f(x\cos \left ( \theta +\pi \right ) +y\sin \left ( \theta +\pi \right ) )}f\ dfd\theta \] and \(\cos \left ( \theta +\pi \right ) =-\cos \theta \) and \(\sin \left ( \theta +\pi \right ) =-\sin \theta \), hence the above becomes\[ \Delta ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( -f,\theta \right ) e^{j2\pi \left ( -f\right ) (x\cos \theta +y\sin \theta )}f\ dfd\theta \] Let \(\alpha =-f\), the above becomes\begin {align*} \Delta & ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{\alpha =0}^{-\infty }} F\left ( \alpha ,\theta \right ) e^{j2\pi \left ( \alpha \right ) (x\cos \theta +y\sin \theta )}\left ( -\alpha \right ) \ \left ( -d\alpha \right ) d\theta \\ & ={\displaystyle \int \limits _{\theta =0}^{\pi }} \left [ -{\displaystyle \int \limits _{\alpha =0}^{-\infty }} F\left ( \alpha ,\theta \right ) e^{j2\pi \left ( \alpha \right ) (x\cos \theta +y\sin \theta )}\left ( -\alpha \right ) \ d\alpha d\theta \right ] \\ & ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{\alpha =-\infty }^{0}} F\left ( \alpha ,\theta \right ) e^{j2\pi \left ( \alpha \right ) (x\cos \theta +y\sin \theta )}\left ( -\alpha \right ) \ d\alpha d\theta \end {align*}

Since \(\alpha \) is free variable, we call rename it as we wish. Let it be called \(f\) again, then the above becomes\begin {equation} \Delta ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=-\infty }^{0}} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\left ( -f\right ) \ dfd\theta \tag {6} \end {equation} The above completes the manipulation needed in the second term, replace the above equation (6) into equation (5), then the 2D inverse polar Fourier transform becomes\begin {align*} h\left ( x,y\right ) & ={\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=0}^{\infty }} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\ f\ dfd\theta \\ & +{\displaystyle \int \limits _{\theta =0}^{\pi }} {\displaystyle \int \limits _{f=-\infty }^{0}} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\left ( -f\right ) \ df d\theta \end {align*}

Using \(\left \vert f\right \vert \) we can combine the above 2 terms into one\begin {equation} h\left ( x,y\right ) ={\displaystyle \int \limits _{\theta =0}^{\pi }} \left [ {\displaystyle \int \limits _{f=-\infty }^{\infty }} F\left ( f,\theta \right ) e^{j2\pi f(x\cos \theta +y\sin \theta )}\ \left \vert f\right \vert \ df\right ] d\theta \tag {7} \end {equation} Converting the outside integral to Riemann sum, and recalling that we have \(M\) projections (or equivalently \(M\) angles), we obtain\begin {equation} h(x,y) =\frac {\pi }{M}{\displaystyle \sum \limits _{m=0}^{M-1}} \left [ {\displaystyle \int \limits _{f=-\infty }^{\infty }} F\left (f,m\frac {\pi }{M}\right ) e^{j 2\pi f \left ( x\cos \left ( m\frac {\pi }{M}\right ) +y \sin \left ( m\frac {\pi }{M}\right ) \right )}\ \left \vert f\right \vert \ df\right ] \tag {8} \end {equation} Now we are ready to use the central slice theorem as follows. In equation (4) above, the function \(F\left ( f,\theta \right ) \) is the 2D Fourier transform along a radial line at angle \(\theta \). But this is the same as the Fourier transform of a projection taken at the same angle \(\theta .\) Hence we can replace \(F\left ( f,m\frac {\pi }{M}\right ) \) in (8)  with the Fourier transform of \(g\left ( p,\theta \right ) \), which we already found its DFT in (1) as \(S\left ( k,\theta _{m}\right ) .\)

Recalling that \(\Delta f=\frac {2W}{N}\) equation (8) above becomes\begin {align} h(x,y) & =\frac {\pi }{M}{\sum \limits _{m=0}^{M-1}}\frac {2W}{N}{\sum \limits _{n=-\frac {N}{2}}^{\frac {N}{2}-1}}S\left ( n\frac {2W}{N},m\frac {\pi }{M}\right ) e^{j2\pi \left ( i\frac {2W}{N}\right ) \left ( x\cos \left ( m\frac {\pi }{M}\right ) +y\sin \left ( m\frac {\pi }{M}\right ) \right ) }\,\left \vert n\frac {2W}{N}\right \vert \,\ \nonumber \\ & =\frac {2W\pi }{M}{\sum \limits _{m=0}^{M-1}}\left ( \frac {1}{N}{\sum \limits _{n=-\frac {N}{2}}^{\frac {N}{2}-1}}\overset {\text {filtered projection}}{\overbrace {S\left ( n\frac {2W}{N},m\frac {\pi }{M}\right ) }}\left \vert n\frac {2W}{N}\right \vert e^{nj\frac {4W\pi }{N}\left ( x\cos \left ( m\frac {\pi }{M}\right ) +y\sin \left ( m\frac {\pi }{M}\right ) \right ) }\right ) \tag {9} \end {align}

And now we see why \(h\left ( x,y\right ) \) is called a filtered backprojection. The spectrum \(S\left ( k,\theta _{m}\right ) \) of each projection is being multiplied by \(\left \vert n\frac {2W}{N}\right \vert \) or \(\left \vert f\right \vert \) before the 2D inverse Fourier transform is obtained. The filter whose spectrum given by \(\left \vert f\right \vert \) is called RAM-LAK

The use of this filter before performing backprojecting, causes the resulting image to be much sharper and details shown more clearly than without this filter.

We see the reason for this from the shape of the filter in the frequency domain. Low pass frequencies are attenuated, this is what reduces the blurring effect, while higher frequencies are amplified, which corresponds to sharpening those parts of the image where sudden changes occur, causing boundaries between portions of the image to become more apparent. The following diagram, shows a 2D reconstruction generated by simulation using Matlab, where one backprojection was done without the use of the above filter, and one was done with the use of the above filter. We see that the one that used the filter is sharper, while the other without the filter is more blurred.

We see that without the use of filtered backprojection, computed tomography will not be practical as images reconstructed could not be useful for medical images purposes.  We see that with RAM-LAK, there is now streak lines across the image, which are not there in the non filtered image. However, the eye can tolerate these streaking effect and the advantages of filtered backprojection remain clear over the unfiltered image.

In the above derivation, the RAM-LAK filter came out naturally from the use of the Jacobian of transformation. The above derivation follows that given by Kak and Slaney book. However, in this report more details and steps are shown and explained.

We note that in (9) above, that the outside loop accumulates each generated 2D image after each projections 2D inverse Fourier transform has been found. The 2D Inverse Fourier Transform can be implemented using 2D IFFT for speed.

There are other filters that can be used to replace the RAM-LAK filter with. The simulation software supports other types of filters.

The Matlab simulation also shows how the final backprojection images converges to the original image as more projections are added. The following is a screen shot of the application GUI

5 Conclusion

We have obtained an expression for filtered backprojection for parallel lines CT. Which is the following\[ h\left ( x,y\right ) =\frac {2W\pi }{MN}{\displaystyle \sum \limits _{m=0}^{M-1}} \left [ {\displaystyle \sum \limits _{n=-\frac {N}{2}}^{\frac {N}{2}-1}} S\left ( n\frac {2W}{N},m\frac {\pi }{M}\right ) e^{nj\frac {4W\pi }{N}(x\cos \left ( m\frac {\pi }{M}\right ) +y\sin \left ( m\frac {\pi }{M}\right ) )}\ \left \vert n\frac {2W}{N}\right \vert \right ] \] Where \(M\) is the number of projections, \(N\) is the number of samples per projection (All projections are sampled using the same number of samples).

Hence, given an \(x\) and \(y\) coordinates, we can obtain the pixel value at that location using the above expression. \(W\) is the largest spatial frequency in all the projections, and is used to determine the Nyquist sampling frequency. \(S\left ( k,\theta _{m}\right ) \) is the DFT of the projection at angle \(\theta _{m}\) and is found using FFT for speed.

We have shown above the steps needed to obtain a 2D reconstruction from a set of projections taken at different angles. In the above we used a number of digital signal processing techniques such as FFT, Central Slice theorem and Nyquist theorem. This shows the importance of signal processing in image processing and in medical imaging specifically.

6 References

1. Principles of computerized tomographic imaging by Kak and Slaney
2. Matlab radon and iradon help.
3. http://rad.usuhs.mil/rad/home/comptom.html
4. One image was obtained from http://labspace.open.ac.uk
5. My own reports written during summer 2008 for final project for Msc Applied Mathematics at CSUF Mathematics department.