Slice-projection theorem and Radon transforms
Ever did a CT-Scan ? Well the maths behind the reconstruction of your inside business really are fascinating. This post covers the Slice-projection theorem that links projections, Fourier transforms and slices. Then we will proceed to reconstruct an image like it is done everyday in CT-Scans or X-rays.
But first, 2D Fourier transforms.
Fourier transforms, in 2D
One may be fully comfortable with 1-dimensional Fourier transforms, but it's less common to work in higher dimensions. No worry, the concept stays the same.
In 1-dimension, the Fourier coefficient for frequency $f$ is the result of the integral of the signal multiplied with a complex phasor of frequency $f$. The scalar product between the signal and a given phasor can be thought of a similarity measure between two signals: if they match, the greater the Fourier coefficient. Actually there's more than that: a phase difference is also computed by doing the scalar product (which is complex), but the idea holds still.
$$\mathcal{F}(s, f)= \int_{-\infty}^{\infty} s(t)\ e^{-2j\pi f t}dt,\quad \forall\ f \in \mathbb R.$$
You can compute the Fourier transform for any frequency $f$.
In 2-dimensions, your signal $s(x, y)$ is 2 dimensional. It could be an image for example, and the basis is also 2-dimensional. This means we have now two frequencies along which the integral can be evaluated, one for each axis. We will name those two frequencies $u$ and $v$ for the axis $x$ and $y$ respectively.
In order to compute the 2D Fourier transform, we now integrate our signal (2D) multiplied by this new 2D basis function of $u$ and $v$ with respect to $x$ and $y$.
$$\mathcal{F}(s, u, v)= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}s(x, y)\ e^{-2j\pi(ux+vy)}dxdy,\quad \forall\ u, v \in \mathbb R^2.$$
But what does $e^{-2j\pi(ux+vy)}$ looks like ? Below is a plot of its real part for a given $x, y$ domain for various couples of $u$ and $v$ with this profile:
Remember, only the real part of the complex value is here plotted. The frequencies coordinates $(u, v)$ are plotted on the right graph. Notice how the spatial frequencies of the sinusoidal pattern on the left changes as well as its direction.
If we want to compute the (discrete) 2D Fourier transform from here we could, very crudely, multiply our image with each basis (every $u, v$ couple) and sum. In real life, don't compute 2D DFT like this, use well optimized Faster Fourier Transform implementations (like from numpy or scipy).
So, what does a 2D Fourier transform looks like ? Let me take an example image first (a Shepp-logan phantom):
I will apply a 2D FFT and plot the absolute value (logscale) of the (complex) result:
Every point $u, v$ is the result of an integral with a pattern that looks like the figure $2$. Points close to the center represent low frequencies, while the borders are occupied by higher frequencies up to Shannon limit (half the sampling rate). Yes, there is also a sampling rate even in images and it's related to the resolution.
If we think about this image as being a sampled version of a continuous real image, then the resolution dictates the spacing between each pixel, or sampling points. Because each pixel is akin to a sample, each pixel spacing is actually the sampling period. The unit do not have a clear signification here, since there is no real physical phenomenon behind my example image.
The Shannon limit here is however easy to grasp by following the 1D analogy: given an image resolution, the couple $u, v$ shouldn't yield faster spatial frequencies than the spacing between two pixels can resolve. In other words, the sinusoidal pattern I showed earlier shouldn't go above one half cycle between two pixels, otherwise the frequency would be higher than half the sampling frequency.
Fun fact: the center point in actually the integral of the whole signal (because $e^{0} = 1$). This observation can be used in various ways to our advantage.
Now that 2D FFT are introduced, let's dive in the Slice-Projection theorem!
Slice-Projection theorem
This beautiful theorem states that the two following things are equivalent:
- A slice of the 2D Fourier transform passing going through the origin.
- The 1D FFT of the projection of the original image onto that slice.
At first, it's really not clear why that should be the case. Let me clarify things a bit.
By projection we mean the line integral (which we'll later call Radon transform) of the original image onto a segment. Again, an animation is worth thousands of words:
As you can see on the animation, the right yellow curve is the projection of the original image (the two white square) onto the yellow line along the white lines. Each white line is one integration line that yield one x value on the right. By varying the angle, another projection is obtained. This is what I mean by projection or line integral.
And what do I mean by slices of 2D Fourier transform ? Well, like if you were to slice the 2D FFT image and look at the value along that slice:
Well, to rephrase, the Projection-Slice theorem states that the Fourier slice you are observing on the right of previous figure is the same thing as the 1D FFT of the projection on the same slice. Neat.
Why is that neat you may ask? Well, you probably can see it coming: in CT-scans you have this kind of projection. Because you are trying to diagnose something without killing the human under-test, you can't just open him in two part, do your own business in order to check what is wrong, and fix back the body.
Instead, you can shoot many X-rays, or high frequency electromagnetic radiations (order of peta to hexahertz) from one side and check what you get on the other side. Because biological tissues absorb (and not all in the same way) those rays, you indeed get as a result an integration of the human body along lines (X-rays). Those individual projections are of no use if you can't get back to the actual body composition but, because the Fourier-Slice theorem exists, I will show how one can recover the full picture.
But before, can we develop an intuition on why is the Projection-Theorem true? My first question when I encountered this theorem was where the projection could arise from such a formulation? Remember, we are actually trying to prove that a slice into the result of the following equation reduces down to the Fourier transform of the projection onto that slice:
$$\mathcal{F}(s, u, v)= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}s(x, y)\ e^{-2j\pi(ux+vy)}dxdy,\quad \forall\ u, v \in \mathbb R^2.$$
Before rigorous mathematical proof, I will explain how I now intuitively understand it: for any $u, v$ couple, the pattern of $e^{-2j\pi(ux+vy)}$ is constant for only one given direction (which is perpendicular to the sinusoidal direction). When thinking about integration, we can either consider integration along $x$ and $y$, but also change variables and integrate along the invariant direction and its perpendicular (the slice). The integral along the invariant direction really is the product of the image and the constant slice. The constant slice can then be factorized out and what is left is the integral of the signal along that line, or the projection.
This is where the projection arises, the sinusoidal pattern can be factorized out for a specific direction, which is perpendicular to the slice of interest.
Enough talking, let's prove the idea:
$$\mathcal{F}(s, u, v)= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}s(x, y)\ e^{-2j\pi(ux+vy)}dxdy,\quad \forall\ u, v \in \mathbb R^2.$$
We can parametrize a slice (or line) of angle $\theta$ by substituting $u$ and $v$ by $k\cos\theta$ and $k\sin\theta$.
$$\mathcal{F}(s, k\cos\theta, k\sin\theta)= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}s(x, y)\ e^{-2j\pi(k\cos\theta x+k\sin\theta y)}dxdy$$
We will then introduce the following change of variables:
$$ \begin{aligned} \alpha &= x\cos\theta + y\sin\theta\cr \beta &= -x\sin\theta + y\cos\theta\cr \end{aligned} $$
Which is just really a rotation of angle $\theta$. In order to rewrite the double integral, we need to compute the jacobian $J(\alpha, \beta)$ (or more precisely its determinant).
$$ \int \int s(x, y),dxdy = \int \int s(h(\alpha, \beta), g(\alpha, \beta)),\mathopen|J(\alpha, \beta)\mathclose| d\alpha d\beta $$
Because our change of variable is a rotation, we already know that its determinant will be equal to 1 (rotations do not alter areas), but nonetheless:
$$ J(\alpha, \beta)= \begin{bmatrix} \frac{\partial x}{\partial \alpha} & \frac{\partial x}{\partial \beta} \cr \frac{\partial y}{\partial \alpha} & \frac{\partial y}{\partial \beta} \cr \end{bmatrix}= \begin{bmatrix} \cos\theta & -\sin\theta \cr \sin\theta & \cos\theta \cr \end{bmatrix} $$
And indeed,
$$ \mathopen|J(\alpha, \beta)\mathclose| = \cos^2\theta + \sin^2\theta=1$$
Rewriting, we get:
$$\mathcal{F}(s, k\cos\theta, k\sin\theta)= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}s(\alpha\cos\theta - \beta\sin\theta, \alpha\cos\theta + \beta\sin\theta)\ e^{-2j\pi k\alpha}d\alpha d\beta$$
Now, notice that $e^{-2j\pi k\alpha}$ do not depend on $\beta$, which mean we can rearrange terms:
$$\int_{-\infty}^{\infty} \left[\int_{-\infty}^{\infty}s(\alpha\cos\theta - \beta\sin\theta, \alpha\cos\theta + \beta\sin\theta),d\beta\right] e^{-2j\pi k\alpha}d\alpha$$
The middle term is the X-ray projection $p(\theta)$ we talked about earlier onto the slice of angle $\theta$! The signal $s$ is simply rotated in order to align with the axis $\beta$ so that we can integrate it as a projection. This leads us to:
$$\boxed{\mathcal{F}(s, k\cos\theta, k\sin\theta)= \int_{-\infty}^{\infty}p(\theta) e^{-2j\pi k\alpha}d\alpha}$$
Which is the definition of the Fourier transform in 1 dimension.
So here is theorem proved: the slice of angle $\theta$ of the 2D Fourier transform $\mathcal{F}(s)$ is the 1D Fourier transform of the projection $p(\theta)$.
The Radon transform
Building things gradually, I will now develop the Radon transformation. In fact, you already know what is does: see Figure. $5$ ? Swipe all angles from $0$ to $\pi$ and stack all the projections in one big matrix. This is the Radon transform. Let's to it for our phantom ( from $-\pi/4$ to $\pi/4$ for now, I will explain why later):
The bottom right plot is the result of the Radon transform for every angle. This representation is also called a sinogram, since an individual point in an input image would result in a sinusoidal pattern at the output.
Now that all projections were computed, how do we get back to the original image ? Well remember the Fourier-Slice theorem, each of the rows of this sinogram is actually one slice of the 2D Fourier transform. This mean that if we take the FFT (row-wise) and place each transformed row on a rotating slice, we would be reconstructing (with errors in precision) the original 2D Fourier result.
See, each FFT row on the right (which is now the result of the 1D FFT of each sinogram row) has to be mapped to the correct location in the polar coordinates on the left. This is simply a polar conversion, nothing too fancy here.
Let's see how it looks on our real sinogram:
You can see that we are missing data in the four corner of the 2D FFT, this is because my projection has the same number of rays for all orientations where the actually slice is actually changing in length. This is not an issue for my explanations purposes, we are just missing out high frequency details. In synthesis, what have we done until there ?
- Projected the original image onto orientated slices in order to compute the Radon transform.
- Taken the FFT row wise of the projections
- Placed those FFT in a polar fashion.
The last piece missing consists in taking the inverse 2D Fourier transform of our polar representation. If everything is well and good, the Projection-Slice theorem should allow us to reconstruct our original image.
First, we need to interpolate this data that does not belong to a grid (remember, individual points are on polar coordinates). Once the interpolation on the grid is done, the 2D inverse FFT is computed. I show below the results of the inverse FFT result on the left and the original image on the right.
Success! This does indeed looks like our phantom. Many artifact are however visible, and for various reasons:
- We are a finite number of rays on a finite number of orientations. We are loosing information here.
- We are missing out the corner for reasons explained before.
- The last step involves an interpolation of the polar FFT slices on a regular grid.
Better attention to details in implementation as well as better rays / angles resolution would make this result better. Actually, while playing with the concept a few days ago, I found my Python implementation of the Radon transform too slow.
The transformation is not complex per se, but it involves mainly for loops, which Python is really bad at. Before doing that I implemented algorithmic optimizations, like pre-computing valid rays indices but the bottleneck here was really the multiple iterations. This is why I switched language and reimplemented my algorithm in Rust, where I gained approximatively a factor of 1000 with the same algorithm. If you would like to fizzle around with it, I made the Python and Rust implementations (as well as Python bindings for the Rust implementation) available on Github.
I hope you found the Projection-Slice theorem interesting and that your curiosity was tickled by those concepts. See you next time.
