Nov 14,2008 Compiled on September 7, 2023 at 11:48pm

Digital signal processing plays an increasingly important rule in the ﬁeld of medical imaging. This report examines in details the mathematical formulation behind an important medical imaging method called Computed Tomography or CT.

The mathematics of CT are outlined showing the central role played by spatial Discrete Fourier Transform (DFT) and the 2D Inverse DFT in the formulation of the method. The theory of reconstruction of a 2D medical image from a sequence of 1D projections taken at diﬀerent angles between zero and \(\pi \) is described.

Projections are generated by applying the Radon transform to the original image at diﬀerent angles. Only the parallel rays case is discussed in this report.

To reconstruct the 2D image from the sequence of projections, Filtered Backprojection (FB) method is used. Filtering is performed the frequency domain. In the report only the RAM-LAK ﬁlter is considered, but other ﬁlters are also possible.

The use of RAM-LAK ﬁlter allows for much improved 2D image reconstruction. The use of FFT in the implementation of the DFT and IDFT algorithms makes this medical image method practical

1 Introduction

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

2.1 Solving the problem algebraically

2.2 Solving the problem using the Fourier central Slice theorem

2.3 Solving the problem using ﬁltered backprojection

3 Radon transform introduction

4 Image reconstruction

5 Conclusion

6 References

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

2.1 Solving the problem algebraically

2.2 Solving the problem using the Fourier central Slice theorem

2.3 Solving the problem using ﬁltered backprojection

3 Radon transform introduction

4 Image reconstruction

5 Conclusion

6 References

We are familiar with the standard x-ray imaging, the type one encounters at a doctor or a dentist oﬃce. Brieﬂy, 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 diﬀerently, and will do so in proportion to the type and amount of tissue it passes though. This attenuation occurs due to the fact that diﬀerent types of tissue (for example, bone vs. muscle vs. soft tissue) absorbs diﬀerent amount of radiation.

The recording of the varying intensities of the x-rays leaving the body (the x-ray ﬂux hitting the detector
surface) gives a 2D image which reﬂects the content of the section of the body that the x-rays passed
through.^{1}

Due to overlapping and shadowing between the internal tissues in the body, the ﬁnal 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 speciﬁc angle \(\theta \). When the beams hits the detector on the other side of the body, the ﬂux 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 ﬁrst 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 ﬁrst phase is mathematically called the radon transform. The operation which accomplishes the second phase is the ﬁltered backprojection. We start by reviewing the ﬁrst phase, showing how the projection is ﬁrst generated. The equation of the projection is then used in the second phase.

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 simpliﬁcation 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

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 inﬁnity many solutions

Hence the best we can do to obtain a solution using the above approach, is to ﬁnd 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.

The following diagram illustrates the main idea behind this approach

The following diagram illustrates the main idea behind this approach. In ﬁltered backprojection, the projection spectrum is ﬁrst ﬁltered, then backprojected directly into a 2D space, and additional projections and backprojections are accumulated on top of this, which results in the ﬁnal 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 ﬁltered backprojection operation.

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

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 speciﬁc \(p\) and speciﬁc \(\theta \) satisﬁes \(x\cos \theta +y\sin \theta =p\)

Assume now there exist a function \(f\left ( x,y\right ) \) deﬁned 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 diﬀerential 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 speciﬁc \(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 diﬀerent \(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 ﬂux 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 \]

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 diﬀerent such sequence for diﬀerent 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 deﬁnition 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 ﬁltered 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 diﬀerentials 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 ﬁltered 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 ﬁlter whose spectrum given by \(\left \vert f\right \vert \) is called RAM-LAK

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

We see the reason for this from the shape of the ﬁlter in the frequency domain. Low pass frequencies are attenuated, this is what reduces the blurring eﬀect, while higher frequencies are ampliﬁed, 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 ﬁlter, and one was done with the use of the above ﬁlter. We see that the one that used the ﬁlter is sharper, while the other without the ﬁlter is more blurred.

We see that without the use of ﬁltered 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 ﬁltered image. However, the eye can tolerate these streaking eﬀect and the advantages of ﬁltered backprojection remain clear over the unﬁltered image.

In the above derivation, the RAM-LAK ﬁlter 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 ﬁlters that can be used to replace the RAM-LAK ﬁlter with. The simulation software supports other types of ﬁlters.

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

We have obtained an expression for ﬁltered 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 diﬀerent 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 speciﬁcally.

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

^{1}Thanks for the Image goes to http://labspace.open.ac.uk