Project 3 Wavelets, compression and compressed sensing

Wavelets, compression and compressed sensing in
magnetic resonance imaging

Introduction
In class, we discussed the theory of wavelets, and its application to compression and compressed
sensing. In this project, you will explore these ideas practically using real MRI data. The data we
will use is on OpenNeuro, found at this url:
https://openneuro.org/datasets/ds001499/versions/1.3.0.
As can be read in the ReadMe on this website, there are four patients: CSI1, 2, 3 and 4.
Most of the sessions here were for functional recordings (fMRI), but we are interested in anatomical
recordings (MRI). For patients CSI1-3, the anatomical recordings are from session 16, and for patient
CSI4, the anatomical recordings are from session 10. The images are either T1- or T2-weighted.
The datasets are titled, for example, “sub-CSI1 ses-16 run-01 T1w.nii.gz.” Extracted, these .nii
files are called NIfTI files, which is a standardized file type used in MRI, fMRI, PET and other
functional imaging modalities. They can be read into MATLAB via niftiread, and Python by
using load in the nibabel library. These contain three-dimensional 16-bit integer arrays, in the
spatial domain. That means the IFT has already been applied by the machine.
As usual, you are expected to write a report in which you display and explain all of your results,
along with answering any specific questions asked in this document. While I suggest the use of
MATLAB or Python, any software package may be used so long as you give me your code and
instructions on how to run it. You will need to take wavelet transforms, so MATLAB’s Wavelet
toolbox, or Pythons pywt will be helpful.
I. Navigating the volumes
To start, consider Patient CSI1’s anatomical session, and load in both the T1- and T2-weighted
data. Create a figure with 6 subplots that displays, in grayscale and scaled magnitude (imagesc
in MATLAB), sagittal, coronal and axial cross-sections from each volume. Label each in a format
like “Axial slice 150, T1-weighted.” Make note of how the dimensions of the array corresponds to
anatomical axes.
In your report, use the following table of T1 and T2 values to describe which structures are
shown most brightly. Does T1/T2 weighting reward large or small T values?
Structure T1 (ms) T2 (ms)
White Matter 510 67
Grey Matter 760 77
Cerebrospinal Fluid 2650 280
1
Look at the corresponding json files in the run’s folder. Explain how the inversion time, echo
time and repetition time relate to the weighting. Explain by appealing to either conceptual NMR
principles, or the imaging equations. They don’t give units for the pulse sequence times – what do
you think the units are, guided by the values in the table above?
For the next several problems, we will want a “guiding image,” which will serve as a visual
marker for the success/failure of the methods we will attempt here. I suggest the center-index
fixed-x image.
II. Auxiliary Functions
Although you will probably write several other functions as you complete this project, I think you
will do well to make sure you have these “right off the bat.” This isn’t a requirement but rather a
suggestion.
1. A function that takes an array and a number s between 0 and 100. The returned array is the
same size as the input, but the lowest-magnitude s% values in the input array are 0 in the
output. This is a simple compression.
2. A function that takes two arrays, one of which is considered a ground truth, and computes
the mean square error.
3. A function that takes an array and a name of a wavelet domain, and plots the approximation
and detail images for the given wavelet transform at two levels. They can be shown either
concatenated on one set of axes, or in a single figure in many subplots (former preferred).
This is trivial in Python but at the cost of clarity. In MATLAB, this requires understanding
the wavedec2, appcoef2 and detcoef2 functions. An easy example can be found here.
III. Fourier Domain Compression?
Show the (scaled) magnitude of your guiding image in the Fourier domain. Make sure you are
taking a 2D FT, and not column-wise 1D FTs. Apply fftshift, as well. Show a histogram of this
magnitude with a few hundred bins, and comment on the sparsity. When you take histograms, be
sure your data is reshaped to 1D.
Note that this is the domain in which our images are acquired. If many of these points are lowmagnitude, maybe we can just sparsely sample this domain and not have to do much else! Consider
your guiding image as ground truth, and observe the reconstructed image after zeroing the smallest
Fourier coefficients with several s values. Note also the mean-squared errors. Comment on the
results here in a detailed fashion, showing comparisons between the original and reconstructed
images labeled by error and what percentage of coefficients were zeroed.
You should find that you can’t achieve significant compression on account of the visual quality
reconstructed of the image, even if the mean-squared error appears low. This is because the phase
of the Fourier domain information is very important, even in the very low-magnitude components.
Thereby even though we have sparsity in the Fourier domain, we cannot use this domain for image
compression!
2
IV. Wavelet Domains
For this problem, we will consider three wavelet bases – Haar, Daubuchies 4 and Coiflet 3. Familiarize yourself with the workings of the wavelet transform functions in the toolbox you are using.
To prove to yourself you know what is doing on, take a 2-level 2D wavelet decomposition of a
standard MATLAB image (for example, ‘cameraman’) and display the approximation image, and
the 6 detail coefficient images, using the function you wrote above.
Try plotting the magnitude of the 2D, 2-level wavelet transforms in these three bases to a
standard basis element in 256×256 matrix space (a standard element is one which is 0 everywhere
except one pixel at which it is 1).How do they look? Plot their Fourier transform magnitudes. How
are these elements localized in space/frequency?
Use either MATLAB or Google to see canonical spatial representations of these wavelets and
see if your commentary above matches what you read, if you’d like.
Now display the 2-level decompositions of the guiding image in these three bases. Is the image
sparse in these bases? Show histograms for the total set of coefficients in these three bases after
1-level, 2-level and 3-level decomposition, for a total of 9 histograms. Which domain looks like the
best candidate for compression/compressed sensing?
Now just to be sure you know what you’re doing, pick the domain you like best and reconstruct
the image where the lowest 10% of the coefficients are zeroed. Display this alongside the original
image, and the mean-squared error.
V. How Sparse?
Say the percent of coefficients being deleted is s. Plot the mean squared error between the reconstructed signal and the true signal as a function of s for the guide image for each basis. Do you
see diminishing returns? Plot a few reconstructed images with each basis. Toy around with the
number of levels, and see if it works for images other than your guide image, including one from
the T2 dataset.
Make qualitative statements, and also say how many coefficients we can toss if we can accept
10% error. Use a for loop to see what the max error among all images would be with this number
of coefficients kept.
VI. Compressed Sensing
Write a gradient descent algorithm that simulates the compressed sensing of MRI images based
on sparsity in the wavelet domain. Pick the wavelet domain based on what you see works best
qualitatively/quantitatively in Section V. Ground truth is the guiding image in the k-domain, i.e.
the FT of the guiding image. This is undersampled by a random sensing matrix (you generate
a Bernoulli random vector to perform the undersampling with probability p). The domain of
sparseness is the wavelet transform of the IFT of the k-domain. Show the average convergence of
this algorithm in objective function as a function of iteration for 1000 sensing matrices with the
same p. Do this for many values of p and give qualitative/quantitative arguments for a good p to
use. I am being purposefully vague because I want you to explore the space. Try different values of
the sparsifying parameter λ early on and settle on a good one before you do the rest of your tests.
3
Suggested Completion Timeline
By Week 12 (after Thanksgiving) complete Section I above, write functions (1) and (2) from section
II and consider Section III. By Week 13, Complete Sections III and IV. By Week 14, Complete V
and start VI. By Week 15 be done with the project, spending the last week exploring section VI.
Update your report each week, and come to class with questions each week. You can hand in the
report by the last day of the semester.
4