*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*= 0 (Russ and I divide the Fourier transform into a cosine part

_{f}(θ, k)*C*and a sine part

*S*, see Eq. 11.57 in

*IPMB*).

*C*is a function of

_{f}(θ, k)*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*is plotted below as a function of

_{f}(θ, k)*k*.

**Step 1b: Filter**

To filter, multiply

*C*by |

_{f}(θ, k)*k*|/2

*π*(Eq. 12.39 in

*IPMB*) to get the Fourier transform of the filtered projection

*C*

_{g}(θ, k)with

*S*. The result is shown below. Notice that filtering removes the dc contribution at

_{g}(θ, k) = 0*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*. 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

_{g}(θ, k)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 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 - 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.

## No comments:

## Post a Comment