MATLAB Code Overview
This section provides a brief overview of each of the MATLAB files and a description of how they all work together. We attempt to include as much information as possible, so that others can copy the programs, duplicate our work, and expand on it. With that said, this section does not substitute for a person actually working with the code and the theory behind it. We also request that all copied files maintain the author and Rice University information in the comments at the beginning of each M file. Please arrange the origination statement such that the "help m_file_name" command in MATLAB displays the statement along with the help file description.
All of the code was written using MATLAB version 5.X.X. Most of the code will run using the basic package; however, we may have unknowingly used commands from some of the toolboxes. If this occurs, more than likely, the upsetting command came from the Signal Processing toolbox. MATLAB has announced that it will bundle the Signal Processing, Symbolic, and Controls toolbox in a student package in August, 1999. This bundle should be at a greatly reduced price, and is suppose to include the complete toolboxes without limitation. At the same time, the student version of MATLAB will be released without the matrix size limitation. Please let us know if you have any problems.
Main.m
Main.m is not a MATLAB function, but instead a MATLAB script that
calls the other functions. It sets the main parameters. Explicitly it requires
the number of angles that will be used to reconstruct the 2-D image. Also, the
length of the detector array and the number of evenly spaced detectors on that
array are needed. The inverse of the resolution determines the number of pixels
in the reconstructed image. Implicitly, the horizontal shift, heart growth, and
lung growth, as described in the Thoracic Cavity Phantom
section, are determined in calcproj.m as inputs into the
project.m function. Finally, the "filtertype" argument in backproj.m selects
between fftfilter.m or convfilter.m
, if set to one. All reconstructed images at the beginning of the Results section were run with fftfilter.m. The convfilter.m was a successful experiment shown later in the Results section.
Project.m
Project.m calculates the integral of the densities on a line
from point (X1,Y1) to (X2,Y2) of a synthetic 2-D slice model of a human thoracic
cavity (chest). The (0,0) point is in the lower, left corner. The "shift"
argument allows a horizontal translation of the entire structure, while "grow_hrt"
and "grow_lung" allows a dilation of the "heart", and "lungs & chest wall", respectively.
The model is contained within a square box with vertices at (0,0) and (1,1).
Any coordinate set can be given for the line, including inside the model.
"Shift" is limited to +/- 0.04, "grow_hrt" is limited to 0 to 0.04, and "
grow_lng" is limited from 0 to 0.08. By looping on project, multiple projection
lines can be built with various angles, as needed for the simulation.
Modgraph.m
Modgraph.m displays a synthetic 2-D slice model
of a human chest using some of the input arguments needed by the
project.m function. "Shift" allows a horizontal translation
while "grow_hrt" allows a dilation of the "heart", and "grow_lng" allows
a dilation of the "lungs" and "chest wall". The model is contained within
a square box with vertices at (0,0) and (1,1). "Shift" is limited to +/-0.04,
"grow_hrt" is limited from 0 to 0.04, and "grow_lng" is limited from
0 to 0.08. The inverse of "res" determines the number of pixels in both the x
and y directions. A suggestion for "res" is 0.004. Due to MATLAB rounding errors,
modgraph.m does not give a completely accurate image at
low resolutions, but it generates the image relatively fast compared to
projgraf.m. Also, the code in modgraph.m
provides a clear description of the equations for the seven ellipses used
to generate the image.
Projgraf.m
Projgraf.m displays and returns a synthetic 2-D slice model
of a human chest using some of the input arguments needed by project.m.
It works by inputting the same x and y coordinates for the coordinate pair, used
in the project.m function, to identify a point
within the image. The returned value from project.m is the density
at that point. Projgraf.m requires a
relatively long time to loop through all points in the grid, but the output
is an image which should be used for comparison with the output from the back-projection
algorithm. The arguments for projgraf.m are "shift" which
allows a horizontal translation, "grow_hrt" which allows a dilation of the "heart", and
"grow_lng" which allows a dilation of the "lungs" and "chest wall". The model is
contained within a square box with vertices at (0,0) and (1,1). "Shift" is limited to +/-0.04,
"grow_hrt" is limited from 0 to 0.04, and "grow_lng" is limited from 0 to 0.08. "Res" is the rectangular sampling grid spacing. The inverse of "res" determines the number of pixels for both the
x and y directions. A suggestion for "res" is 0.004.
Calcproj.m
Calcproj.m calculates the projections at the specified angles with
the specified detector structure using project.m.
Backproj.m
Backproj.m is the main back-projection function that
reconstructs the image pixel by pixel. It uses either fftfilter.m or convfilter.m to filter the projections. To back-project, it uses bproj.m.
Fftfilter.m
Fftfilter.m is a filtering function that is used
by backproj.m to filter the projections before back-projection.
The filter is linearly increasing as frequency increases. At the highest frequencies,
1/16 of the the spectrum is zeroed. Fftfilter.m performs
the filtering operation in the frequency domain.
Convfilter.m
Convfilter.m is a filtering function that can be used
by backproj.m to filter the projections before back-projection, if the "filtertype" argument is set to one.
It performs the filtering operation in the time domain, and provides improved performance
over fftfilter.m.
Bproj.m
Bproj.m is the function used to calculated the
samples of the filtered projections needed for back-projection. Each pixel
is then calculated by integrating the interpolated values from each
filtered projection. This function is called by backproj.m.
Zeroneg.m
Zeroneg.m is a function that moves all negative values
in an image to zero. Typically, a negative value in an image has no realistic meaning.
Within the context of this work, a negative value implies that a negative amount of x-rays
were present. The minimum amount possible is zero. This function in not part of the
main.m script. It was used as one of the post-processing experiments.
Shiftfix.m
Shiftfix.m can be used to repair the artifacts caused by a horizontal shift of the phantom during the acquisition of the projections. The function works by recentering all projections to an estimated midpoint. The estimated midpoint is a sinusoid centered on the length of detectors with amplitude proportional to the offset of the first projection.
Enhance.m
Enhance.m calculates the squared error between an original tomographic image and six reconstructed images that have gone through one or more filtering processes. It also calculates the total squared error from the results.
Trans.m
Trans.m takes in an input image and outputs the normalized sums of adjacent row elements, then differences between adjacent row elements. This is part of the foundation for a Haar wavelet transform.
Itrans.m
Itrans.m takes in an input image and outputs the inverse Haar wavelet transform.
Rpeak.m
Rpeak.m takes an input matrix, scans over the entire image in 1x3 blocks, and returns only local peaks into the output matrix while zeroing out everything else.