diff --git a/assets/img/project_fbp_reconstruction.png b/assets/img/project_fbp_reconstruction.png new file mode 100644 index 0000000000000000000000000000000000000000..d444e9f0a2d89243030a803ea5ff91cd55b3bbfe Binary files /dev/null and b/assets/img/project_fbp_reconstruction.png differ diff --git a/project.md b/project.md index 79bebb2377a5d36fcafed176048714e6a0b33b2f..57b9fc8e3ae52a533d94aee0f3f65f92b573494a 100644 --- a/project.md +++ b/project.md @@ -32,6 +32,7 @@ projections at different angles. ## Fourier Reconstruction + ### Fourier Slice Theorem Let $$f : \mathbb{R}^2 \rightarrow \mathbb{R}$$ be a density function. If we map density to @@ -111,7 +112,7 @@ algorithms out there that come closer to the theoretically ideal The popular [gridrec](https://www.ncbi.nlm.nih.gov/pubmed/23093766) method is one. -### Implementation +### Results 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 @@ -148,6 +149,7 @@ It would be nice if we could avoid the interpolation required for Fourier recons simplest way of doing so, and still the most popular method of performing image reconstruction, is called filtered back projection. + ### Theory Let's hop back to the continuous theory for a moment. The Fourier reconstruction technique is based @@ -189,6 +191,11 @@ $$ f(x, y) = \int_0^\pi \mathcal{F}^{-1}[q_\theta](x \cos \theta + y \sin \theta) \mathrm{d} \theta $$ +Since multiplication in the frequency domain is equivalent to convolution in the spatial domain, +$$\mathcal{F}^{-1}[q_\theta]$$ is simply a filtered version of $$p_\theta$$. Hence the name filtered +back projection. + + ### Discretization We'll replace the continuous Fourier transforms with DFTs as before. @@ -210,6 +217,13 @@ them directly for quadrature. For $$r$$, on the other hand, the regular samples line up with the values $$x \cos \theta + y \sin \theta$$ that we want. So we'll still have to do some interpolation. But now it's a simple 1D interpolation problem that's much easier to do without introducing as much error. In particular, we only have to do interpolation in the spatial domain, so -errors will only accumulate locally. [Lanczos -resampling](https://en.wikipedia.org/wiki/Lanczos_resampling) would probably be the best, but I'll -just use Catmull-Rom again since I already have it implemented. +errors will only accumulate locally. + + +### Results + +I again use [FFTW](http://www.fftw.org/) for Fourier transforms. For interpolation I just take the +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. + +