*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

*J*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.

_{1}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 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

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 to 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