by Nasser M. Abbasi, Nov 14,2008
We outline the mathematics of CT showing the central role played by spatial
Discrete Fourier Transform (DFT) and the 2D Inverse DFT in the formulation of
the method. We outline the theory of reconstruction of a 2D medical image from
a sequence of 1D projections taken at different angles between 0 and .
Projections are generated by applying the Radon transform to the original
image at different angles. We will only examine the parallel rays case.
To reconstruct the 2D image from the sequence of projections, the Filtered
Backprojection (FB) method is used. Filtering is performed the frequency
domain. In the report we consider only the RAM-LAK filter, but other filters
The use of RAM-LAK filter 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.
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.
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 . 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 The angle is incremented to 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.
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.
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
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
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
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 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.
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 as follows
Assume now there exist a function defined over the region shown above. The integral of this function over the line is
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
The input to CT 2D reconstruction is the continuous function Before we go any further, we need to make sure we keep track of the angle being used in the above expression. Recall that this is the angle that each projection is generated at. This angle goes from 0 to (Since going beyond that will duplicate measurements). Let the number of projections to be generated be . Hence
The next step is to sample and then to obtain the DFT of the sampled sequence. To sample , we assume that the smallest distance between 2 adjacent repeating cycles of intensities is cm (or in millimeter). What this means is that largest spatial frequency present in the projection data is . Therefore, the projection needs to be sampled at frequency no less than (Nyquist sampling theorem). Hence we will us as the sampling frequency
The result of sampling will be the sequence of numbers where
The following diagram helps illustrate the above formulation
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
Hence and , and the Jacobian of transformation is the determinant of , hence
Split the outside integral in (4) which goes from 0 to into 2 halves as follows
Recalling that equation (8) above becomes
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.
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
We have obtained an expression for filtered backprojection for parallel lines CT. Which is the following
Hence, given an and coordinates, we can obtain the pixel value at that location using the above expression. is the largest spatial frequency in all the projections, and is used to determine the Nyquist sampling frequency. is the DFT of the projection at angle 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.