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 an answer.

The above screenshot above 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 t

*o* *Ï€*, 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.

- 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.
- 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.
- 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.
- 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.
- 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.
- My bucket list includes finding an analytical example of filtered back projection when
*f(x,y)* depends on direction. Wouldn’t that be cool!
- 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.