Description
Introduction
In class, we discussed some fundamental concepts in image processing. In this project, you will
apply this newfound knowledge to a real volume scan – a CT scan of the thorax made available by
Mulrecon. It is the first scan on their webpage:
http://www.castlemountain.dk/atlas/index.php?page=mulrecon.
The concepts tested here will be: 1) Noise, 2) Blurring, 3) Filtering, and 4) Edge detection.
You have 2 weeks to complete this project. It should be submitted as a zip file emailed to me
(b.frost@columbia.edu) containing the MATLAB or Python files and the documentation as a PDF.
The documentation should contain the relevant figures and analyses asked for in each section below,
and doesn’t need to be too fancy.
MATLAB functions that may be useful: rgb2gray, squeeze, imfilter, imgaussfilt, medfilt2,
edge.
I. Slicing
The stack is given in the form of 237 JPG images. Each image is a slice through the thorax at
a fixed z. Call the horizontal (left-to-right) x, the vertical (up-to-down) y and the direction of
increasing image index z (this is a right-handed coordinate system). This format is clunky – you
cannot naturally see other cross-sections of the scan (fixed x or fixed y) looking at these images.
Create a 3D array containing “the volume” – that is, a stack of the images so that the first index
corresponds to the x pixel, the second to the y pixel and the third to the z pixel. The volume should
be grayscale. If you have done this correctly, the slice indexed by (:,:,p) should correspond to
the pth jpeg image.
In the documentation, present one slice in each of the three dimensions near the middle of the
subject. Label the axes with the corresponding dimensions, and label the pixel at which the slice
is taken (for example, title a figure x = 215). Comment on the units of the axes, how they relate
to meters, and what that means about the aspect ratio of your images.
II. Salt-N-Pepa
Salt and pepper noise or speckle noise is common in CT. This means that noise appears in the form
of spurious very high- or very low-intensity pixels. In small neighborhoods, this could skew the
mean downward or upward due to these statistical outliers. However, the spurious pixels should
not affect the median.
In larger neighborhood, the statistical outliers probably “average out.” We might intuitively
expect the mean and the median to behave similarly in large neighborhoods. However, filtering
1
over large neighborhoods has a larger smoothing effect, which is not helpful if we want to recognize
finer features.
For one of the slices you presented in section I, show the signal median-filtered and mean-filtered
in a 3×3 and 5×5 neighborhood (4 images total). In the text, describe the effects of either filter,
compare the median with the mean filter at each neighborhood size.
For many cross-sections, you may see that both filters are very similar. To exaggerate this effect
further, you can insert salt and pepper noise yourself! Write a function that replaces a pixel with
a white pixel with probability p, a black pixel with probability p, or leaves the pixel unaffected
with probability 1 − 2p. Show the effects of the mean and median filter with neighborhood size 3×3
for three reasonable p values (in your report, show both the noisy and filtered images, so 6 images
total), and describe the difference.
III. Gauss’ Wrath
We will explore the effects of Gaussian noise, Gaussian filtering and Gaussian models of blurring,
one by one.
Gaussian Noise – Gaussian noise is present in most measurements due in part to thermal
effects on the hardware side. Traditionally, local averaging through Gaussian filters is the best
way to deal with this sort of noise. Gaussian filters are low-pass filters, so they have undesirable
smoothing effects if they are made too large. This is very similar to the mean/median filters
above. Add Gaussian noise to the signal with different standard deviations (pick two reasonable
but different values). Improve the signal quality by using Gaussian filters on these images. For
these filters, just pick the standard deviation which gives you the qualitatively nicest results. Report
these values and how they relate to the noise standard deviation, and discuss if or why this makes
sense.
Gaussian Blurring – Blurring in imaging is usually caused by the image being out of focus.
However, all images are out of focus (by the diffraction limit), so there will always be some blurring
which we model as convolution with the PSF (impulse response) of the system. The PSF is usually
an Airy pattern in 2D or 3D, but this is well-modeled by a multidimensional Gaussian pattern.
Assume that the CT scan is not a scan at all, but instead a perfect representation of our sample.
That is, the CT intensity map is precisely the reflectivity of the sample. We will now model the
imaging of this sample by a device with an anisotropic point spread function.
Say the resolution as quantified by the full width half max (FWHM) value of the PSF is δx in the
x and y directions (lateral) and δz in the z direction (axial). Model this PSF by a three-dimensional
Gaussian pattern in space – what is the form of this PSF (normalized so that it integrates to one)?
It might help to consider the one-dimensional Gaussians in each direction – what standard deviation
would give them the correct full-width value at half of the maximum? Once you have those three
standard deviations you can write the formula quickly up to the normalization constant, which you
can simply compute in MATLAB. Present this 3D formula in the documentation.
For each slice you presented in part I, present the slice blurred in the proper way with δx = 2
mm, δz = 1 mm. This is non-trivial, as the pixels in these images are rectangular as well. You
may want to use the 2D functional form applied to a meshgrid, and one of convolution or an FFT
multiplication equivalent. Compare these to the original scans.
2
III. The Plural of Foramen is Apparently Foramina
The spaces between your vertebrae are called the interverbral foramina. Look at slice y = 180 –
you can make out some features of the subject’s spine, including several intervertebral disks and
foramina. The disks appear as subtler intensity variations, with the foramina being more starkly
contrasted from the vertabrae.
Use MATLAB’s built-in Canny edge finder (you should recall this method from class) on the
raw slice. You should be able to make out the intervertebral disks and foramena as well as many
other edges. Do you see an issue with the result here?
Use this Canny edge detection method on four variations on the slice: a) a 3×3 median-filter is
applied, b) the object is imaged by the system as described in part II (fancy 2D Gaussian filter),
c) the image is impacted by one of your Gaussian noise paradigms from part II, and d) the image
is filtered by a Gaussian filter with standard deviation 3 pixels in both directions.
Describe the effects on the ability to discriminate the intervertebral foramena and disks. What
about other seemingly spurious edges?
Finally, for your own sake, look at how pretty the Canny edge detected images are for the
z-slices!
3
