Friday, May 27, 2016

An Analytical Example of Filtered Back Projection

One of my hobbies is to find tomography problems that can be solved analytically. I know this is artificial—all tomography for medical imaging uses numerical computation—but as a learning tool analytical analysis helps build insight. I have some nice analytical examples using the Fourier method to solve the tomography problem (see homework problems 26 and 27 in chapter 12 of Intermediate Physics for Medicine and Biology), but I don't have a complete analytical example to illustrate the filtered back projection method (see a previous post for a partial example). Russ Hobbie and I do include a numerical example in section 12.6 of IPMB. I have always wondered if I can do that example analytically. Guess what. I can! Well, almost.

Start with a top-hat function for your object
A mathematical function of the top-hat function, which is part of an analytical example of filtered back projection.
If we set x = 0, we can plot it as function of y.
A plot of the top-hat function, which is part of an analytical example of filtered back projection.
The projection of this function is given in IPMB; Homework Problem 36 asks the reader to derive it.
A mathematical expression for the projection of the top-hat function, which is part of an analytical example of filtered back projection.
Because the object looks the same from all directions, the projection is independent of the angle. Below is a plot of the projection as a function of x'. It is identical to the top panel of IPMB's Figure 12.22.
A plot of the projection of the top-hat function, which is part of an analytical example of filtered back projection.
The next step is to filter the projection, which means we have to take its Fourier transform, multiply the transform by a high-pass filter, and then do the inverse Fourier transform. The Fourier transform of the projection is
A mathematical expression for the Fourier transform of the projection of the top-hat function, which is part of an analytical example of filtered back projection.
This integral is not trivial, but Abramowitz and Stegun’s Handbook of Mathematical Functions With Formulas, Graphs and Mathematical Tables contains (Page 360, Equation 9.1.20)
An integral expression for a Bessel function.
where J1 is a first-order Bessel function (see Homework Problem 10). Because the projection is an even function, the sine part of the Fourier transform vanishes.

Filtering is easy; multiply by |k|/2π. The result is
A mathematical expression for the Fourier transform of the filtered projection of the top-hat function, which is part of an analytical example of filtered back projection.
To find the inverse Fourier transform, we need
An integral needed to calculation the filtered projection of the top-hat function, which is part of an analytical example of filtered back projection.
This integral appears in Abramowitz and Stegun (Page 487, Equation 11.4.37)
The filtered projection of the top-hat function, which is part of an analytical example of filtered back projection.
After some simplification (which I leave to you), the filtered projection becomes
A mathematical expression of the filtered projection of the top-hat function, which is part of an analytical example of filtered back projection.
Below is a plot of the filtered projection, which you should compare to the middle panel of Fig. 12.22. It looks the same as the plot in IPMB, except in the numerical calculation there is some ringing near the discontinuity that is not present in the analytical solution
A plot of the filtered projection of the top-hat function, which is part of an analytical example of filtered back projection.
The final step is back projection. Because the projection is independent of the angle, we can calculate the back projection along any radial line, such as along the y axis
A mathematical expression for the process of backprojection.
If |y| is less than a, the back projection is easy: you just get 1. Thus, the filtered back projection is the same as the object, as it should be. If |y| is greater than a, the result should be zero. This is where I get stuck; I cannot do the integral. If any reader can solve this integral (and presumably show that it gives zero), I would greatly appreciate hearing about it. Below is a plot of the result; the part in red is what I have not proven yet. Compare this plot to the bottom panel of Fig. 12.22.
A plot of the filtered back projection of the top-hat function, which is part of an analytical example of filtered back projection.

What happens if you do the back projection without filtering? You end up with a blurry image of the object. I can solve this case analytically too. For |y| less than a, the back projection without filtering is
A mathematical expression for the back projection of the top-hat function without filtering, which is part of an analytical example of filtered back projection.
which is 4a times the complete elliptic integral of the second kind
The definition of an elliptic integral.
For |y| greater than a, you get the more complicated expression
A mathematical expression for the back projection of the top-hat function without filtering, which is part of an analytical example of filtered back projection.
which is the incomplete elliptic integral of the second kind
The definition of the incomplete elliptic integral of the second kind.
The trickiest part of the calculation is determining the upper limit on the integral, which arises because for some angles the projection is zero (you run into the same situation in homework problem 35, which I highly recommend). Readers who are on the ball may worry that the elliptic integral is tabulated only for kappa less than one, but there are ways around this (see Abramowitz and Stegun, Page 593, Equation 17.4.16). When I plot the result, I get
A plot of the back projection of the top-hat function without filtering, which is part of an analytical example of filtered back projection.
which looks like Fig. 12.23 in IPMB.

So, now you have an analytical example that illustrates the entire process of filtered back projection. It even shows what happens if you forget to filter before back projecting. For people like me, the Bessel functions and elliptic integrals in this calculation are a source of joy and beauty. I know that for others they may be less appealing. To each his own.

I’ll rely on you readers to fill in the one missing step: show that the filtered back projection is zero outside the top hat.

No comments:

Post a Comment