Mathematical Theory
Introduction
Computer tomography (CT) entails the reconstruction of a function f from line integrals of f. For for our project, the original function was a 2-D matrix in MATLAB that represented an image we created. The basic CODE for our pictures can be found on this link. Numerous projections, a vector of sums of the values of the image points taken along a specific direction, must be taken and manipulated in order to prepare them for backprojection. This process involves filtering and taking the image's transform, the most common transform used being the Radon Transform, in order to reconstruct our original image.
Radon Transform
CT uses two primary mathematical tools: the Radon Transform and the backprojection algorithm. It is based around the properties that exist in a space of n-dimensions of real numbers. Our process focuses around the Radon Transform, i.e. projection operator, to create our projections. The Radon Transform of a function is defined as the line integral along a line at an angle, phi, from the y-axis and at a distance s from the origin. Mathematically we write it as follows:
The delta functions allows us to implicitly define the line of integration due to the properties of the delta function where it is only equal to one when the statement inside is equal to zero. As can be seen from Figure 1.1, once our original image is passed through the radon transform, we have a collection of projections gathered from N different angles around our image. We will denote our projections as g(s, phi) Again, we had to turn to MATLAB in order to simulate this reconstruction. Our CODE for finding the projections can be found here.
Projection-Slice Theorem & Backprojection
The reconstruction of f(x,y) from our set of projections is based on the Fourier Slice Theorem. This theorem describes the idea that the 2-D Fourier transform of f(x,y) can be found from the set of 1-D Fourier transforms of our projections, g(s, phi). As we have discussed all semester in Elec 301, the simple Fourier transform of our projections, which we will denote as G(phi, w), is as follows:
Next, we can perform some tedious calculus that involves plugging the Radon Transform we described above into the previous equations. Then after simplifying we see the defining idea behind the slice theorem. We see that G(phi, w) is just the 2-D Fourier transform, F(u,v), of f(x,y) evaluated at u = w*sin(phi) and v = -w*cos(phi) - this is simply the line that our projections were taken on.
Next, we simply take the 2-D inverse Fourier transform and make a change of variables so that we are now in polar coordinates. Then once we substitute in G(phi, w), we arrive at the following defining equation:
Now, after dealing with the Radon transform and the basic Fourier transform we discussed in class, we have developed a relationship between our projections taken above and the image that we wish to reconstruct. One thing we must note from this equation is that our projections are multiplied by |w| in the Fourier domain. What this amount to is applying a filter to each projection the Fourier domain before we move back into the time domain. This idea of applying a filter will be discussed further in the report. We will represent this filtered, back-projection as follows:
Now that we have a relationship between our projections, which we already have found a way to find and have written the MATLAB code for, we are ready to create a new program that will use our projections and, through filtering, will reconstruct our original image. In order to implement this, we used three programs. The CODE for all three of these can be found at this page. The first program, pojfilter.m, simply applied the various filter we created to the projections. Then, back.m performed the heart of the algorithm. Finally, we created go.m in order to be able to go from start to finish with one simple command in MATLAB.