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
If we set
x = 0, we can plot it as function of
y.
The projection of this function is given in
IPMB; Homework Problem 36 asks the reader to derive it.
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.
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
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)
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
To find the inverse Fourier transform, we need
This integral appears in Abramowitz and Stegun (Page 487, Equation 11.4.37)
After some simplification (which I leave to you), the filtered projection becomes
Below is a plot of the filtered projection, which you should compare with 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
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
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 with the bottom panel of Fig. 12.22.
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
which is 4
a times the complete
elliptic integral of the second kind
For |
y| greater than
a, you get the more complicated expression
which is 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
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