From d6628e55c1abc55c026c4c91e65b3c7ad8a1515f Mon Sep 17 00:00:00 2001 From: Erik Strand <erik.strand@cba.mit.edu> Date: Tue, 30 Apr 2019 20:47:06 -0400 Subject: [PATCH] Flesh out CT theoretical discussion --- project.md | 60 ++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 54 insertions(+), 6 deletions(-) diff --git a/project.md b/project.md index fdb1a63..f867fe7 100644 --- a/project.md +++ b/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. -- GitLab