Another Analytical Example of Filtered Back Projection

Some people collect stamps, aged bottles of wine, or fine art. I collect analytical examples of tomographic reconstruction: determining a function f(x,y) from its projections F(θ,x'). Today I’ll share my latest acquisition: an example of filtered back projection. I discussed a similar example in a previous post, but you can’t have too many of these. This analysis illustrates the process that Russ Hobbie and I describe in Section 12.5 of Intermediate Physics for Medicine and Biology.

The method has two steps: filtering the projection (that is, taking its Fourier transform, multiplying by a filter function, and doing an inverse Fourier transform) and then back projecting. We start with a projection F(θ,x'), which in the clinic would be the output of your tomography machine.

This projection is independent of the angle θ, implying that the function f(x,y) looks the same in all directions. A plot of F(θ,x') as a function of x' is shown below.

I suggest you pause for a minute and guess f(x,y) (after all, our goal is to build intuition). Once you make your guess, continue reading.

Step 1a: Fourier transform

The Fourier transform of the projection F(θ,x') is

with S_f(θ, k) = 0 (Russ and I divide the Fourier transform into a cosine part C and a sine part S, see Eq. 11.57 in IPMB). C_f(θ, k) is a function of k, the wave number (also known as the spatial frequency). I won’t show all calculations; to gain the most from this post, fill in the missing steps. C_f(θ, k) is plotted below as a function of k.

Step 1b: Filter

To filter, multiply C_f(θ, k) by |k|/2π (Eq. 12.39 in IPMB) to get the Fourier transform of the filtered projection C_g(θ, k)

with S_g(θ, k) = 0. The result is shown below. Notice that filtering removes the dc contribution at k = 0 and causes the function to fall off more slowly at large |k| (it's a high-pass filter).

Step 1c: Inverse Fourier transform

Next use Eq. 11.57 in IPMB to calculate the inverse Fourier transform of C_g(θ, k). Initially I didn’t think the required integral could be solved analytically. I even checked the best integral table in the world, with no luck. However, when I typed the integral into the WolframAlpha website, it gave me the answer

The expression contains the cosine integral, Ci(z). After evaluating it at the integral’s end points, using limiting expressions for Ci(z) at small z, and expending much effort, I derived the filtered projection G(θ,x')

A plot of G(θ,x') is shown below.

Step 2: Back projection

The back projection (Eq. 12.30 in IPMB) of G(θ,x') requires some care. Because f(x,y) does not depend on direction, you can evaluate the back projection along any radial line. I chose y = 0, which means the rotated coordinate x' = x cosθ + y sinθ in the back projection is simply x' = x cosθ. You must examine two cases.

Case 1: |x| less than a

In this case, integrate using only the section of G(θ,x') for |x'| less than a.

Instead of integrating from θ = 0 to π, use symmetry, multiply by two, and integrate from 0 to π/2.

At this point, I was sure I could not integrate such a complicated function analytically, but again WolframAlpha came to the rescue.

(Beware: Wolfram assumes you integrate over x, so in the above screenshot x is my θ and b is my x.) The solution contains inverse hyperbolic tangent functions, but once you evaluate them at the end points you get a simple expression

Case 2: |x| greater than a

When |x| is greater than a, angles around θ = 0 use the section of G(θ,x') for |x'| greater than a and angles around θ = π/2 use the section for |x'| less than a.

The angle where you switch from one to the other is θ = cos^-1(a/x). (To see why, analyze the right triangle in the above figure with side of length a and hypotenuse of length x.) The back projection integral becomes

If you evaluate this integral, you get zero.


Final Result

Putting all this together, and remembering that f(x,y) doesn’t depend on direction so you can replace x by the radial distance, you find

We did it! We solved each step of filtered back projection analytically, and found f(x,y).

I’ll end with a few observations.

1. Most clinical tomography devices use discrete data and computer computation. You rarely see the analysis done analytically. Yet, I think the analytical process builds insight. Plus, it’s fun.
2. Want to check our result? Calculate the projection of f(x,y). You get the function F(θ,x') that we started with. To learn how to project, click here.
3. Regular readers of this blog might remember that I analyzed this function in a previous post, where you can see what you get if you don’t filter before back projecting. 
4. I have a newfound respect for WolframAlpha. It solved integrals analytically that I thought were hopeless. In addition, it’s online, free, and open to all.
5. Most filtered back projection algorithms don’t filter using Fourier transforms. Instead, they use a convolution. I think Fourier analysis provides more insight, but that may be a matter of taste. 
6. My bucket list includes finding an analytical example of filtered back projection when f(x,y) depends on direction. Wouldn’t that be cool! 
7. Remember, there is another method for doing tomography: the Fourier method (see Section 12.4 in IPMB). Homework Problems 26 and 27 in Chapter 12 provide analytical examples of that technique.