Skip to content
Snippets Groups Projects
Commit dae2f5a3 authored by Erik Strand's avatar Erik Strand
Browse files

Fix links; add link to code

parent bf4f5277
No related branches found
No related tags found
No related merge requests found
......@@ -9,10 +9,19 @@ title: Final Project
- [Radon Transform](http://www-math.mit.edu/~helgason/Radonbook.pdf) by Sigurdur Helgason
- [Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency
Information](../assets/pdf/2004_Candes_Romberg_Tao.pdf) by Candes, Romberg, and Tao
- [Total Variation and Tomographic Imaging from Projections](../assets/pdf/2011_hansen_jorgensen.pdf)
Information](assets/pdf/2004_Candes_Romberg_Tao.pdf) by Candes, Romberg, and Tao
- [Total Variation and Tomographic Imaging from Projections](assets/pdf/2011_Hansen_Jorgensen.pdf)
by Hansen and Jørgensen
## Code
All my code is in this [repo](https://gitlab.cba.mit.edu/erik/funky_ct/tree/develop), named in honor
of the [Radon transform](https://en.wikipedia.org/wiki/Radon_transform)'s more eccentric
[cousin](https://en.wikipedia.org/wiki/Funk_transform). Building requires a modern C++ compiler and
cmake. [Libpng](http://www.libpng.org/pub/png/libpng.html) must be installed as well. All told it's
about 1,800 lines (excluding third party code), but there's some unnecessary duplication in there
since I've been favoring velocity over hygiene.
## Fourier Reconstruction
......@@ -30,7 +39,7 @@ $$
p(x) = \int_\mathbb{R} f(x, y) dy
$$
Meanwhile, the [Fourier transform](../notes/fourier_transform.html) of $$f$$ is
Meanwhile, the [Fourier transform](notes/fourier_transform.html) of $$f$$ is
$$
\hat{f}(u, v)
......@@ -86,7 +95,7 @@ First, we'll want our Fourier transforms to be discrete Fourier transforms (DFTs
continuous and discrete Fourier transforms are effectively interchangeable, as long as the functions
we work with are mostly spatially and bandwidth limited, and we take an appropriate number of
appropriately spaced samples. You can read more about these requirements
[here](../notes/fourier_series.html).
[here](notes/fourier_series.html).
Second, since we combine the DFTs of the projections radially, we'll end up with samples (of the 2D
Fourier transform of our image) on a polar grid rather than a cartesian one. So we'll have to
......@@ -105,18 +114,18 @@ I used [FFTW](http://www.fftw.org/) to compute Fourier transforms. It's written
I started with this image of a
[brain](https://commons.wikimedia.org/wiki/File:FMRI_coronal_scan.jpg) from Wikimedia Commons.
![brain](../assets/img/project_brain.png)
![brain](assets/img/project_brain.png)
Fourier reconstruction produces a nice [sinogram](https://en.wikipedia.org/wiki/Radon_transform).
Each row is one projection, with angle going from 0 at the top to $$\pi$$ at the bottom, and $$r =
0$$ down the middle column. You can clearly see the skull (a roughly circular feature) unwrapped to
a line on the left side of the sinogram.
![sinogram](../assets/img/project_sinogram.png)
![sinogram](assets/img/project_sinogram.png)
The reconstruction, however, isn't so clean.
![brain](../assets/img/project_fourier_reconstruction.png)
![brain](assets/img/project_fourier_reconstruction.png)
The Fourier library I'm using is solid (and indeed it reproduces images very well even after
repeated applications), so the error must be coming from my interpolation code. Indeed, the high
......@@ -131,12 +140,12 @@ anyway so I'll move on.
I got a number of interesting failures before getting my code to work correctly.
![outtake](../assets/img/project_cool_1.png)
![outtake](../assets/img/project_cool_2.png)
![outtake](../assets/img/project_cool_3.png)
![outtake](../assets/img/project_cool_4.png)
![outtake](assets/img/project_cool_1.png)
![outtake](assets/img/project_cool_2.png)
![outtake](assets/img/project_cool_3.png)
![outtake](assets/img/project_cool_4.png)
These all started as attempts to project and reconstruct a [disk](../notes/fourier_examples.html),
These all started as attempts to project and reconstruct a [disk](notes/fourier_examples.html),
though I did experiment once interesting mistakes started happening. They mostly result from
indexing errors and an issue with my integration with FFTW. The latter problem relates to the
periodicity of discrete Fourier transforms and the resulting ambiguity in frequency interpretation.
......@@ -148,13 +157,13 @@ For this application it's very important that the center of the polar coordinate
resampling is right at the DC sample. So though the raw Fourier transform of the image looks like
this (real and imaginary parts shown separately),
![transform](../assets/img/project_transform_swapped_re.png)
![transform](../assets/img/project_transform_swapped_im.png)
![transform](assets/img/project_transform_swapped_re.png)
![transform](assets/img/project_transform_swapped_im.png)
we want to permute the quadrants so that it looks like this:
![transform](../assets/img/project_transform_re.png)
![transform](../assets/img/project_transform_im.png)
![transform](assets/img/project_transform_re.png)
![transform](assets/img/project_transform_im.png)
Then the center of the image can be the center of the polar coordinate system. (Note: to make these
I linearly mapped the full range of each image to [1, e], then applied the natural logarithm. So
......@@ -250,17 +259,17 @@ I again use [FFTW](http://www.fftw.org/) for Fourier transforms. For interpolati
sample to the left of the desired location, and for integration I just use a left Riemann sum. Even
with these lazy techniques, the reconstruction looks quite good.
![brain](../assets/img/project_fbp_reconstruction.png)
![brain](assets/img/project_fbp_reconstruction.png)
## Total Variation / Compressed Sensing
Unfortunately this is still a work in progress. I started implementing a [modern
approach](../assets/pdf/2011_hansen_jorgensen.pdf), but found I didn't have enough background to
approach](assets/pdf/2011_hansen_jorgensen.pdf), but found I didn't have enough background to
make it work. They gloss over some details that I tried to get around with brute force, only to find
that the resulting problem was computationally intractable. So yesterday I started over,
re-implementing results from one of the original compressed sensing
[papers](../assets/pdf/2004_Candes_Romberg_Tao.pdf). However I did learn some things from my failed
[papers](assets/pdf/2004_Candes_Romberg_Tao.pdf). However I did learn some things from my failed
attempts...
......@@ -278,24 +287,24 @@ of the Fourier and filtered back projection algorithms implemented as matrices.
and create a sinogram with a single matrix multiplication. As mentioned the matrices are huge, so
here I'll work with a 32 by 32 brain image.
![dft](../assets/img/project_tiny_brain.png)
![dft](assets/img/project_tiny_brain.png)
The DFTs were easy. Very similar to the cosine transform we used in the last problem set.
![dft](../assets/img/project_matrix_dft_re.png)
![dft](../assets/img/project_matrix_dft_im.png)
![dft](assets/img/project_matrix_dft_re.png)
![dft](assets/img/project_matrix_dft_im.png)
To construct the polar interpolation matrix, I rewrote my interpolation routine to drop the weights
it calculated in the appropriate entries in a matrix rather than summing things up as it goes. These
weights depend on the type of interpolation used and the sizes of the images involved.
![dft](../assets/img/project_matrix_polar_re.png)
![dft](../assets/img/project_matrix_polar_im.png)
![dft](assets/img/project_matrix_polar_re.png)
![dft](assets/img/project_matrix_polar_im.png)
Finally it's just some (inverse) DFTs to get the sinogram. They can all be expressed simultaneously
in one matrix.
![dft](../assets/img/project_matrix_sinogram.png)
![dft](assets/img/project_matrix_sinogram.png)
I bump up the pixel count for the polar projections so that the resampling doesn't lose as much
information. Thus the final matrix that performs all three of these operations at once has
......@@ -316,11 +325,11 @@ I scanned a 3d print that I made in [How to Make (almost)
Anything](http://fab.cba.mit.edu/classes/863.18/CBA/people/erik/04_3d_printing_scanning/) last
semester.
![3d print](../assets/img/project_3d_print.jpg)
![3d print](assets/img/project_3d_print.jpg)
I reconstructed it using TomoPy's gridrec implementation.
![3d print](../assets/img/project_3d_print_reconstruction.jpg)
![3d print](assets/img/project_3d_print_reconstruction.jpg)
The results aren't that great, so evidently the settings will require some tweaking. Ultimately I'd
like to set up an easy to use TomoPy/ASTRA toolchain to use with CBA's CT scanner; this is just a
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment