diff --git a/assets/img/project_brain.png b/assets/img/project_brain.png new file mode 100644 index 0000000000000000000000000000000000000000..59eb3c3e87e7e77150e248f7de3f31ee488638ac Binary files /dev/null and b/assets/img/project_brain.png differ diff --git a/assets/img/project_fourier_reconstruction.png b/assets/img/project_fourier_reconstruction.png new file mode 100644 index 0000000000000000000000000000000000000000..059d29c5b63cbce3f8f8b768d7320d7f5d45f07d Binary files /dev/null and b/assets/img/project_fourier_reconstruction.png differ diff --git a/assets/img/project_sinogram.png b/assets/img/project_sinogram.png new file mode 100644 index 0000000000000000000000000000000000000000..906b7c28f3ca9e29c4e9346039bca951e7dff82c Binary files /dev/null and b/assets/img/project_sinogram.png differ diff --git a/project.md b/project.md index e77cf2806f5779db0d880e273591611a805708b2..79bebb2377a5d36fcafed176048714e6a0b33b2f 100644 --- a/project.md +++ b/project.md @@ -105,13 +105,43 @@ appropriately spaced samples. You can read more about these requirements 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 -interpolate. This step is tricky and tends to introduce a lot of error. I'll rather arbitrarily use -[Catmull-Rom interpolation](http://entropymine.com/imageworsener/bicubic/), but there are better +interpolate. This step is tricky and tends to introduce a lot of error. , but there are better algorithms out there that come closer to the theoretically ideal [sinc interpolation](https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula). The popular [gridrec](https://www.ncbi.nlm.nih.gov/pubmed/23093766) method is one. +### Implementation + +I used [FFTW](http://www.fftw.org/) to compute Fourier transforms. It's written in C and is very +[fast](http://www.fftw.org/benchfft/). I implemented my own polar resampling routine. It uses a +[Catmull-Rom interpolation](http://entropymine.com/imageworsener/bicubic/) kernel. + +I started with this image of a +[brain](https://commons.wikimedia.org/wiki/File:FMRI_coronal_scan.jpg) from Wikimedia Commons. + + + +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. + + + +The reconstruction, however, isn't so clean. + + + +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 +frequency content looks ok, but there's a lot of error in low frequency content. This is encoded in +the middle of the Fourier transform, which is what is most distorted by the polar resampling. I +could implement a resampling routine specifically designed for polar resampling, or indeed +specifically for polar resampling for CT reconstruction, but there are better algorithms out there +anyway so I'll move on. + + ## Filtered Back Projection It would be nice if we could avoid the interpolation required for Fourier reconstruction. The