Digital signal processing plays an increasingly important rule in the field 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 different angles between zero and is described.
Projections are generated by applying the Radon transform to the original image at different 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 filter is considered, but other filters are also possible.
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.
X-ray projections in CT. Thanks for the Image goes to http://labspace.open.ac.uk
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 . 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.
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.
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
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
Or in matrix form
The above is in the form . The 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.
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.
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.
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
The equation of the line can be written in 2 ways. The standard way is
where is the slope and is the intercept. It can also be written in terms of the parameters and as
Any point on the line with specific and specific satisfies
Assume now there exist a function defined over the region shown above. The integral of this function over the line is
where is a differential element of the line
It would be simpler to express the above integral in terms of and . To do that, a trick is used with the help of the delta function. The above integral can be written as
Hence for a specific the above will integrate over the line . The above is the radon transform of over the line .
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 over each one of these lines (since all these lines are parallel to line , then all of them will have the same , but each they will have different 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 and . 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
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 to (Since going beyond that will duplicate measurements). Let the number of projections to be generated be . Hence
Hence means the projection at angle , which is the angle . For example, , and and so on.
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 above means that is sampled with an interval of width cm. Let the number of samples obtained from the above sampling frequency be (and we assume it to be even as we can always reduce the sampling interval to obtain the next even value of )
The result of sampling will be the sequence of numbers where
There will be different such sequence for different angle
The following diagram helps illustrate the above formulation
Once a projection is sampled, its DFT is obtained as follows
Since the sampling frequency is , hence the frequency axis will extend from to and as shown in the following diagram
Let us now review what we have done so far. We sampled a projection at some angle 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 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 . The following diagram illustrates the central slice theorem
Let the 2D Fourier transform of be , hence from the definition of the inverse 2D Fourier transform we write
Where is an approximation to the original image .
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, is radial distance from the original (the DC frequency) and is the angle measured counter clock wise from the axis.
Hence and , and the Jacobian of transformation is the determinant of , hence
Hence the integration differentials can now be changed as
Using (3) in (2), we obtain, after changing the limit of integration
The above is the 2D inverse Fourier transform in polar coordinates.
Split the outside integral in (4) which goes from to into 2 halves as follows
The above 2 integrals can be combined into one integral using the property of the symmetry of the Fourier Transform . To do that, we manipulate the second term in the above as follows. First, let the second term be called for now, i.e.
and since the object is real valued, the 2D Fourier transform is symmetric. Hence , then the above becomes
and and , hence the above becomes
Let , the above becomes
Since is free variable, we call rename it as we wish. Let it be called again, then the above becomes
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
Using we can combine the above 2 terms into one
Converting the outside integral to Riemann sum, and recalling that we have projections (or equivalently angles), we obtain
Now we are ready to use the central slice theorem as follows. In equation (4) above, the function is the 2D Fourier transform along a radial line at angle . But this is the same as the Fourier transform of a projection taken at the same angle Hence we can replace in (8) with the Fourier transform of , which we already found its DFT in (1) as
Recalling that equation (8) above becomes
And now we see why is called a filtered backprojection. The spectrum of each projection is being multiplied by or before the 2D inverse Fourier transform is obtained. The filter whose spectrum given by 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
We have obtained an expression for filtered backprojection for parallel lines CT. Which is the following
Where is the number of projections, is the number of samples per projection (All projections are sampled using the same number of samples).
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.