Flesh out CT theoretical discussion

project.md

@@ -2,17 +2,41 @@
 title: Final Project
 ---
 
-# CT Imaging from Scratch
+# Computed Tomography
+
 
 ## Background
 
 - [Radon Transform](http://www-math.mit.edu/~helgason/Radonbook.pdf) by Sigurdur Helgason
 
-## 2d Reconstruction
 
-In two dimensions, the theory of image reconstruction from projections is pretty simple. Assume some
-density function $$f : \mathbb{R}^2 \rightarrow \mathbb{R}$$ (with compact support). The projection
-of this density function to the x axis is
+## Introduction
+
+Computed Tomography (CT) can turn 2d projections of a 3d shape like these (TODO insert image)
+
+into a 3d model like this (TODO insert image).
+
+This is pretty cool. In particular, voxels in the computed model can tell you the density at
+specific locations *inside* the scanned object, whereas pixels in the projections can only tell you
+the average density along lines that pass *through* the object. So we're getting a lot more out of
+those images than meets the eye.
+
+How does it work? That's what I'd like to understand. In particular, I hope to clearly explain the
+basic principles of CT, use them to implement a basic reconstruction from scratch, and experiment
+with existing solutions to get a sense of the state of the art.
+
+
+## Theory
+
+Let's start in 2d, since it will make everything simpler. In this case we have an image like this
+(TODO insert image), and we'd like to reconstruct it only knowing its 1d projections at different
+angles.
+
+We can describe the image as a density function $$f : \mathbb{R}^2 \rightarrow \mathbb{R}$$. Think
+of $$f(x, y)$$ as the brightness of the image at $$(x, y)$$. We'll assume that this function is
+defined everywhere, but is always zero outside our image.
+
+The projection of this density function to the x axis is obtained by integrating along y:
 
 $$
 p(x) = \int_\mathbb{R} f(x, y) dy
@@ -25,7 +49,7 @@ $$
 = \int_\mathbb{R} \int_\mathbb{R} f(x, y) e^{-2 \pi i (u x + v y)} dx dy
 $$
 
-Note that the slice along the $$u$$ axis in frequency space is described by
+Now comes the key insight. The slice along the $$u$$ axis in frequency space is
 
 $$
 \begin{align*}
@@ -36,3 +60,27 @@ $$
 &= \hat{p}(u)
 \end{align*}
 $$
+
+So the Fourier transform of the 1d projection is a 1d slice through the 2d Fourier transform of the
+image. Since the x axis is arbitrary (we can rotate the image however we want), this works for other
+angles as well. That is, if you take the projection onto the line that makes an angle $$\theta$$
+with the x axis, its Fourier transform is equal to the slice along the same line in frequency space.
+
+Conceptually this tells us everything we need to know about the reconstruction. First we take the
+Fourier transform of each projection individually. Then we combine them by arranging them radially.
+Finally we take the inverse Fourier transform of the resulting 2d function. We'll end up with the
+image.
+
+It also tells us how to generate the projections, given that we don't have a 1d x-ray machine. First
+we take the Fourier transform of the image. Then we extract radial slices from it. Finally we take
+the inverse Fourier transform of each slice. These are the projections. This will come in handy for
+generating testing data.
+
+
+## Discretization
+
+Of course the clean math of the theory has be modified a bit to make room for reality. First of all,
+our projections will be represented as discrete pixels, not perfect functions. So we'll want our
+Fourier transforms to be DFTs. And when we combine the transforms of the projections, we won't
+naturally end up with a nice 2d grid of values on which to perform an inverse DFT. So we'll have to
+interpolate.