Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • develop
  • pset7
  • pset6
  • pset3
5 results
Blame Permalink
project.md 3.37 KiB 
title: Final Project

Computed Tomography

Background

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

Meanwhile, the Fourier transform of f is

\hat{f}(u, v) = \int_\mathbb{R} \int_\mathbb{R} f(x, y) e^{-2 \pi i (u x + v y)} dx dy

Now comes the key insight. The slice along the u axis in frequency space is

\begin{align*} \hat{f}(u, 0) &= \int_\mathbb{R} \int_\mathbb{R} f(x, y) e^{-2 \pi i u x} dx dy \\ &= \int_\mathbb{R} \left( \int_\mathbb{R} f(x, y) dy \right) e^{-2 \pi i u x} dx \\ &= \int_\mathbb{R} p(x) e^{-2 \pi i u x} dx \\ &= \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.