Skip to content
Snippets Groups Projects
Select Git revision
  • master default protected
  • develop
  • pset7
  • pset6
  • pset3
5 results

project.md

Blame
    • project.md 4.54 KiB 
    title: Final Project

    Computed Tomography

    Background

    Introduction

    Computed Tomography (CT) is used to 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 2d 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.

    I'm going to describe the algorithms 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.

    Fourier Reconstruction

    Fourier Slice Theorem

    Let f : \mathbb{R}^2 \rightarrow \mathbb{R} be a density function. If we map density to brightness we can view f as describing an image. We'll assume that this function is defined everywhere, but is always zero outside some finite neighborhood of the origin (say, the bounds of the image).

    The projection of f 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 project f onto the line that makes an angle \theta with the x axis, its Fourier transform is equal to the slice along the corresponding line in frequency space. This result is known as the Fourier Slice Theorem, and is the foundation of most reconstruction techniques.

    Conceptually this tells us everything we need to know about the reconstruction. First we take the 1D Fourier transform of each projection. 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 reconstructed 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

    Naturally the clean math of the theory has be modified a bit to make room for reality. In particular, we only have discrete samples of f (i.e. pixel values) rather than the full continuous function (which in our current formalism may contain infinite information). This has two important implications.

    First, we'll want our Fourier transforms to be discrete Fourier transforms (DFTs). Luckily the continuous and discrete Fourier transforms are effectively interchangeable, as long as the functions we work with are mostly spatially and bandwidth limited, and we take an appropriate number of appropriately spaced samples. You can read more about these requirements here.

    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, but there are better algorithms out there that come closer to the theoretically ideal sinc interpolation. The popular gridrec method is one.