Fourier Series

One of the most important mathematical techniques for a physicist is the Fourier series. I discussed Joseph Fourier, the inventor of this method, previously in this blog. In the 4th edition of Intermediate Physics for Medicine and Biology, Russ Hobbie and I discuss the Fourier series in Sections 11.4 and 11.5.

The classic example of a Fourier series is the representation of a periodic square wave: y(t) = 1 for t between 0 and T/2, and y(t) = -1 for t between T/2 and T, where T is the period. The Fourier series represents this function as a sum of sines and cosines, with frequencies of k/T, where k is an integer, k = 0, 1, 2, …. The square wave function y(t) is odd, so the contributions of the cosine functions vanish. The sine functions contribute for half the frequencies, those with odd values of k. The amplitude of each non-zero frequency is 4/πk (Eq. 11.34 in IPMB), so the very high frequency terms (large k) don’t contribute much.

Being able to calculate the Fourier series is nice, but much more important is being able to visualize it. When I teach my Medical Physics class (PHY 326), based on the last half of IPMB, I stress that students should “think before you calculate.” One ought to be able to predict qualitatively the Fourier coefficients by inspection. Being able to understand a mathematical calculation in pictures and in physical terms is crucially important for a physicist. The Wikipedia article about a square wave has a nice animation of the square wave being built up by adding more and more frequencies to the series. I always insist that students draw figures showing better and better approximations to a function as more terms are added, at least for the first three non-zero Fourier components. You can also find a nice discussion of the square wave at the Wolfram website. However, the best visualization of the Fourier series that I have seen was brought to my attention by one of the PHY 326 students, Melvin Kucway. He found this lovely site, which shows the different Fourier components as little spinning wheels attached to wheels attached to wheels, each with the correct radius and spinning frequency so that their sum traces out the square wave. Watch this animation carefully. Notice how the larger wheels rotate at a lower frequency, while the smaller wheels spin around at higher frequencies. This picture reminds me of the pre-Copernican view of the rotation of planets based on epicycles proposed by Ptolemy.

What is unique about the development of Fourier series in IPMB? Our approach, which I rarely, if ever, see elsewhere, is to derive the Fourier coefficients using a least-squares approach. This may not be the simplest or most elegant route to the coefficients, but in my opinion it is the most intuitive. Also, we emphasize the Fourier series written in terms of sines and cosines, rather than complex exponentials. Why? Understanding Fourier series on an intuitive level is hard enough with trigonometric functions; it becomes harder still when you add in complex numbers. I admit, the math appears in a more compact expression using complex exponentials, but for me it is more difficult to visualize.

If you want a nice introduction to Fourier series, click here or here (in the second site, scroll down to the bottom on the left). If you prefer listening to reading, click here for an MIT Open Courseware lecture about the Fourier series. The two subsequent lectures are also useful: see here and here. The last of these lectures examines the square wave specifically.

One of the fascinating things about the Fourier representation of the square wave is the Gibbs phenomenon. But, I have discussed that in the blog before, so I won’t repeat myself.

What is the Fourier series used for? In IPMB, the main application is in medical imaging. In particular, computed tomography (Chapter 12) and magnetic resonance imaging (Chapter 18) are both difficult to understand quantitatively without using Fourier methods.

As a new year’s resolution, I suggest you master the Fourier series, with a focus on understanding it on a graphical and intuitive level. What is my new year’s resolution for 2014? It is for Russ and I to finish and submit the 5th edition of IPMB to our publishers. With luck, you will be able to purchase a copy before the end of 2015.

Joseph Fourier

The August 2010 issue of Physics Today, published by the American Institute of Physics, contains an article by T. N. Narasimhan about “Thermal Conductivity Through the 19th Century.” A large part of the article deals with Joseph Fourier (17687–1830), the French physicist and mathematician. Russ Hobbie and I discuss Fourier’s mathematical technique of representing a periodic function as a sum of sines and cosines of different frequencies in Chapter 11 of the 4th edition of Intermediate Physics for Medicine and Biology. Interestingly, this far-reaching mathematical idea grew out of Fourier’s study of heat conduction and thermal conductivity. Russ and I introduce thermal conductivity in Homework Problem 15 of Chapter 4 about diffusion. This is not as odd as it sounds because, as shown in the problem, heat conduction and diffusion are both governed by the same partial differential equation, typically called the diffusion equation (Eq. 4.24). The concept of heat conduction is crucial when developing the bioheat equation (Chapter 14), which has important medical applications in tissue heating and ablation.

Narasimhan’s article provides some interesting insights into Fourier and his times.

In 1802, upon his return to France from Napoleon’s Egyptian campaign, Fourier was appointed perfect of the department of Isere. Despite heavy administrative responsibilities, Fourier found time to study heat diffusion. He was inspired by deep curiosity about Earth and such phenomena as the attenuation of seasonal temperature variations in Earth’s subsurface, oceanic and atmospheric variations in Earth’s subsurface, oceanic and atmospheric circulation driven by solar heat, and the background temperature of deep space…

Thermal conductivity, appropriate for characterizing the internal conduction, was defined by Fourier as the quantity of heat per unit time passing through a unit cross-section divided by the temperature difference of two constant-temperature surfaces separated by unit distance… Fourier presented his ideas in an unpublished 1807 paper submitted to the Institut de France.

Fourier was not satisfied with the 1807 work. It took him an additional three years to go beyond the discrete finite-difference description of flow between constant-temperature surfaces and to express heat flow across an infinitesimally thin surface segment in terms of the temperature gradient.

When Fourier presented his mathematical theory, the nature of heat was unknown… Fourier considered mathematical laws governing the effects of heat to be independent of all hypotheses about the nature of heat… No method was available to measure flowing heat. Consequently, in order to demonstrate that his mathematical theory was physically credible, Fourier had to devise suitable experiments and methods to measure thermal conductivity.

It is not widely recognized that in his unpublished 1807 manuscript and in the prize essay he submitted to the Institut de France in 1811, Fourier provided results from transient and steady-state experiments and outlined methods to invert exponential data to estimate thermal conductivity. For some reason, he decided to restrict his 1822 masterpiece, The Analytical Theory of Heat, to mathematics and omit experimental results.

For more insight on Fourier’s life and times, see Keston’s article “Jospeh Fourier: Policitian and Scientist.” It begins

The life of Baron Jean Baptiste Joseph Fourier (1768–1830) the mathematical physicist has to be seen in the context of the French Revolution and its reverberations. One might say his career followed the peaks and troughs of the political wave. He was in turns: a teacher; a secret policeman; a political prisoner; governor of Egypt; prefect of Isère and Rhône; friend of Napoleon; and secretary of the Académie des Sciences. His major work, The Analytic Theory of Heat, (Théorie analytique de la chaleur) changed the way scientists think about functions and successfully stated the equations governing heat transfer in solids. His life spanned the eruption and aftermath of the Revolution; Napoleon's rise to power, defeat and brief return (the so-called Hundred Days); and the Restoration of the Bourbon Kings.

The Fast Fourier Transform

In Chapter 11 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I discuss the fast Fourier transform.

The calculation of the Fourier coefficients using our equations involves N evaluations of the sine or cosine, N multiplications, and N additions for each coefficient. There are N coefficients, so that there must be N² evaluations of the sines and cosines, which uses a lot of computer time. Cooley and Tukey (1965) showed that it is possible to group the data in such a way that the number of multiplications is about (N/2)log₂N instead of N² and the sines and cosines need to be evaluated only once, a technique known as the fast Fourier transform (FFT).

Additional analysis of the FFT is found in the homework problems at the end of the chapter.

Problem 17. This problem provides some insight into the fast Fourier transform. Start with the expression for an N-point Fourier transform in complex notation, Y_k in Eq. 11.29a. Show that Y_k can be written as the sum of two N/2-point Fourier transforms: Y_k = ½[Y_k^e + W^k Y_k^o], where W = exp(-i2π/N), superscript e stands for even values of j, and o stands for odd values.

Numerical Recipes:
The Art of Scientific Computing,
by Press et al.

The FFT is a famous algorithm in the field of numerical methods. Below is how Press et al. describe it in one of my favorite books, Numerical Recipes.

The discrete Fourier transform can, in fact, be computed in O(Nlog₂N) operations with an algorithm called the fast Fourier transform, or FFT. The difference between Nlog₂N and N² is immense. With N = 10⁶, for example, it is the difference between, roughly, 30 seconds of CUP time and 2 weeks of CPU time on a microsecond cycle time computer. The existence of an FFT algorithm became generally known only in the mid-1960s, from the work of J. W. Cooley and J. W. Tukey. Retrospectively, we now know…that efficient methods for computing the DFT [discrete Fourier transform] had been independently discovered, and in some cases implemented, by as many as a dozen individuals, starting with Gauss in 1805!

One “rediscovery” of the FFT, that of Danielson and Lanczos in 1942, provides one of the clearest derivations of the algorithm. Danielson and Lanczos showed that a discrete Fourier transform of length N can be rewritten as the sum of two discrete Fourier transforms, each of length N/2. One of the two is formed from the even-numbered points of the original N, the other from the odd-numbered points…

The wonderful thing about the Danielson-Lanczos Lemma is that it can be used recursively. Having reduced the problem of computing F_k to that of computing F_k^e and F_k^o, we can do the same reduction of F_k^e to the problem of the transform of its N/4 even-numbered input data and N/4 odd-numbered data…

Although there are ways of treating other cases, by far the easiest case is the one in which the original N is an integer power of 2…With this restriction on N, it is evident that we can continue applying the Danielson-Lanczos Lemma until we have subdivided the data all the way down to transforms of length 1…The points as given are the one-point transforms. We combine adjacent pairs to get two-point transforms, then combine adjacent pairs of pairs to get 4-point transforms, and so on, until the first and second halves of the whole data set are combined into the final transform. Each combination takes on order N operations, and there are evidently log₂N combinations, so the whole algorithm is of order Nlog₂N.

This process, called decimation-in-time, is summarized in this lovely butterfly diagram.

A butterfly diagram of a decimation-in-time
Fast Fourier Transform.

The Fourier Series of the Cotangent Function

In Section 11.5 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I describe the Fourier series. I’m always looking for illustrative examples of Fourier series to assign as homework (see Problems 10.20 and 10.21), to explain in class, or to include on an exam. Not every function will work; it must be well behaved (in technical jargon, it must obey the Dirichlet conditions). Sometimes, however, I like to examine a case that does not satisfy these conditions, just to see what happens.

Consider the Fourier series for the cotangent function, cot(x) = cos(x)/sin(x).

The cotangent function, from Schaum's Outlines: Mathematical Handbook of Formulas and Tables.

The function is periodic with period π, but at zero it asymptotically approaches infinity. Its Fourier series is defined as

where

The cotangent is odd, implying that only sines contribute to the sum and a₀ = a_n = 0. Because the product of two odd functions is even, we can change the lower limit of the integral for b_n to zero and multiply the integral by two

To evaluate this integral, I looked in the best integral table in the world (Gradshteyn and Ryzhik) and found

implying that b_n = 2, independent of n. The Fourier series of the cotangent is therefore

When I teach Fourier series, I require that students plot the function using just a few terms in the sum, so they can gain intuition about how the function is built from several frequencies. The first plot shows only the first term (red). It's not a good approximation to the cotangent (black), but what can you expect from a single frequency?

The cotangent function approximated by a single frequency.

The second plot shows the first term (green, solid), the second term (green dashed), and their sum (red). It’s better, but still has a long ways to go.

The cotangent function approximated by two frequencies.

If you add lots of frequencies the fit resembles the third plot (red, first ten terms). The oscillations don’t seem to converge to the function and their amplitude remains large.

The cotangent function approximated by ten frequencies.

The Youtube video below shows that the oscillation amplitude never dies down. It is like the Gibbs phenomenon on steroids; instead of narrow spikes near a discontinuity you get large oscillations everywhere.

The bottom line: the Fourier method fails for the cotangent; its Fourier series doesn’t converge. High frequencies contribute as much as low ones, and there are more of them (infinitely more). Nevertheless, we do gain insight by analyzing this case. The method fails in a benign enough way to be instructive.

I hope this analysis of a function that does not have a Fourier series helps you understand better functions that do. Enjoy!

Another Analytical Example of Filtered Back Projection

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 the answer

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 to π, 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.

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

The Central Slice Theorem: An Example

The central slice theorem is key to understanding tomography. In Intermediate Physics for Medicine and Biology, Russ Hobbie and I ask the reader to prove the central slice theorem in a homework problem. Proofs are useful for their generality, but I often understand a theorem better by working an example. In this post, I present a new homework problem that guides you through every step needed to verify the central slice theorem. This example contains a lot of math, but once you get past the calculation details you will find it provides much insight.

The central slice theorem states that taking a one-dimensional Fourier transform of a projection is equivalent to taking the two-dimensional Fourier transform and evaluating it along one direction in frequency space. Our “object” will be a mathematical function (representing, say, the x-ray attenuation coefficient as a function of position). Here is a summary of the process, cast as a homework problem.

Section 12.4 

Problem 21½. Verify the central slice theorem for the object

(a) Calculate the projection of the object using Eq. 12.29,

Then take a one-dimensional Fourier transform of the projection using Eq. 11.59,

 

(b) Calculate the two-dimensional Fourier transform of the object using Eq. 12.11a,

Then transform F̂(k_x,k_y ) to F̂(θ,k) by converting from Cartesian to polar coordinates in frequency space.

(c) Compare your answers to parts (a) and (b). Are they the same?

I’ll outline the solution to this problem, and leave it to the reader to fill in the missing steps. 

 

Fig. 12.12 from IPMB, showing how to do a projection.

The Projection 

Figure 12.12 shows that the projection is an integral of the object along various lines in the direction θ, as a function of displacement perpendicular to each line, x'. The integral becomes

Note that you must replace x and y by the rotated coordinates x' and y'

You can verify that x² + y²= x'² + y'².

After some algebra, you’re left with integrals involving e^−by'2 (Gaussian integrals) such as those analyzed in Appendix K of IPMB. The three you’ll need are

The resulting projection is

Think of the projection as a function of x', with the angle θ being a parameter.

 

The One-Dimensional Fourier Transform

The next step is to evaluate the one-dimensional Fourier transform of the projection

The variable k is the spatial frequency. This integral isn’t as difficult as it appears. The trick is to complete the square of the exponent

Then make a variable substitution u = x' + ik⁄2b. Finally, use those Gaussian integrals again. You get

This is our big result: the one-dimensional Fourier transform of the projection. Our next goal is to show that it’s equal to the two-dimensional Fourier transform of the object evaluated in the direction θ.

Two-Dimensional Fourier Transform

To calculate the two-dimensional Fourier transform, we must evaluate the double integral

The variables k_x and k_y are again spatial frequencies, and they make up a two-dimensional domain we call frequency space.

You can separate this double integral into the product of an integral over x and an integral over y. Solving these requires—guess what—a lot of algebra, completing the square, and Gaussian integrals. But the process is straightforward, and you get

Select One Direction in Frequency Space

If we want to focus on one direction in frequency space, we must convert to polar coordinates: k_x = k cosθ and k_y = k sinθ. The result is 

This is exactly the result we found before! In other words, we can take the one-dimensional Fourier transform of the projection, or the two-dimensional Fourier transform of the object evaluated in the direction θ in frequency space, and we get the same result. The central slice theorem works.

I admit, the steps I left out involve a lot of calculations, and not everyone enjoys math (why not?!). But in the end you verify the central slice theorem for a specific example. I hope this helps clarify the process, and provides insight into what the central slice theorem is telling us.

Dr. Euler’s Fabulous Formula Cures Many Mathematical Ills

Dr. Euler's Fabulous Formula,
by Paul Nahin.

I’ve decided that Leonhard Euler (1707-1783) is my favorite mathematician, and I’m mad at myself for not getting his name mentioned in the 5th edition of Intermediate Physics for Medicine and Biology. I’ve discussed Euler before in this blog, when I reviewed William Dunham’s book Euler: The Master of Us All. I’ve just finished another book, this one by Paul Nahin, about Dr. Euler’s Fabulous Formula. Of course, Nahin means Eq. 11.28 in IPMB,
e^iθ = cosθ + i sinθ .

I liked the book, in part because Nahin and I seem to have similar tastes: we both favor the illustrations of Norman Rockwell over the paintings of Jackson Pollock, we both like to quote Winston Churchill, and we both love limericks:

I used to think math was no fun,
‘Cause I couldn’t see how it was done.
    Now Euler’s my hero
    For I now see why zero,
Equals e^πi + 1.

Nahin’s book contains a lot of math, and I admit I didn’t go through it all in detail. A large chunk of the text talks about the Fourier series, which Russ Hobbie and I develop in Chapter 11 of IPMB. Nahin motivates the study of the Fourier series as a tool to solve the wave equation. We discuss the wave equation in Chapter 13 of IPMB, but never make the connection between the Fourier series and this equation, perhaps because biomedical applications don’t rely on such an analysis as heavily as, say, predicting how a plucked string vibrates.

Nahin delights in showing how interesting mathematical relationships arise from Fourier analysis. I will provide one example, closely related to a calculation in IPMB. In Section 11.5, we show that the Fourier series of the square wave (y(t) = 1 for t from 0 to T/2 and equal to -1 for t from T/2 to T) is

y(t) = Σ b_k cos(k2πt/T)

where the sum is over all odd values of k (k = 1, 3, 5, ....) and b_k = 4/(π k). Evaluate both expressions for y(t) at t = T/4. You get

π/4 = 1 – 1/3 + 1/5 – 1/7 +…

This lovely result is hidden in IPMB’s Eq. 11.36. Warning: this is not a particularly useful algorithm for calculating π, as it converges slowly; including ten terms in the sum gives π = 3.04, which is still over 3% off.

In Figure 11.17, Russ and I discuss the Gibbs phenomenon: spikes that occur in y(t) at discontinuities when the Fourier series includes only a finite number of terms. Nahin makes the same point with the periodic function y(t) = (π – t)/2 for t from 0 to 2π. He describes the history of the Gibbs phenomena, which arises from a series of published letters between Josiah Gibbs, Albert Michelson, A. E. H. Love, and Henri Poincare. Interestingly, the Gibbs phenomenon was discovered long before Gibbs by the English mathematician Henry Wilbraham.

Fourier series did not originate with Joseph Fourier. Euler, for example, was known to write such trigonometric series. Fourier transforms (the extension of Fourier series to nonperiodic functions), on the other hand, were first presented by Fourier. Nahin discusses many of the same topics that Russ and I cover, including the Dirac delta function, Parseval’s theorem, convolutions, and the autocorrelation.

Nahin concludes with a section about Euler the man and mathematical physicist. I found an interesting connection to biology and medicine: when hired in 1727 by the Imperial Russian Academy of Sciences, it was as a professor of physiology. Euler spent several months before he left for Russia studying physiology, so he would not be totally ignorant of the subject when he arrived in Saint Petersburg!

I will end with a funny story of my own. I was working at Vanderbilt University just as Nashville was enticing a professional football team to move there. One franchise that was looking to move was the Houston Oilers. Once the deal was done, folks in Nashville began debating what to call their new team. They wanted a name that would respect the team’s history, but would also be fitting for its new home. Nashville has always prided itself as the home of many colleges and universities, so a name out of academia seemed appropriate. Some professors in Vanderbilt’s Department of Mathematics came up with what I thought was the perfect choice: call the team the Nashville Eulers. Alas, the name didn’t catch on, but at least I never again was uncertain about how to pronounce Euler.

Central Slice Theorem and Ronald Bracewell

Chapter 12 of the 4th edition of Intermediate Physics for Medicine and Biology deals with images and tomography. One of the key ideas in tomography is the “central slice theorem.” Russ Hobbie and I write in Section 12.4 that

The Fourier transform of the projection at angle θ is equal to the two-dimensional Fourier transform of the object, evaluated in the direction θ in Fourier transform space. This result is known as the projection theorem or the central slice theorem (Problem 17). The transforms of a set of projections at many different angles provide values of C and S [the cosine and sine parts of the 2-d Fourier transform] throughout the k_xk_y plane [frequency space] that can be used in Eq. 12.9a [the definition of the 2-d Fourier transform] to calculate f(x,y).

I consider the central slice theorem to be one of the most important concepts in medical imaging. How was this fundamental idea first developed? The answer to that question provides a fascinating example of how physics and engineering can contribute to medicine.

Ronald Bracewell first developed the central slice theorem while working in the field of radio astronomy. His 2007 New York Times obituary states

Ronald N. Bracewell, an astronomer and engineer who used radio telescopes to make early images of the Sun’s surface, in work that also led to advances in medical imaging, died on Aug. 12 at his home in Stanford, Calif. He was 86…

With his colleagues at Stanford University in the 1950s, Dr. Bracewell designed a specialized radio telescope, called a spectroheliograph, to receive and evaluate microwaves emitted by the Sun…

Later, in the 1970s, the techniques and a formula devised by Dr. Bracewell were applied by other scientists in developing X-ray imaging of tumors, called tomography, and other forms of medical imaging that scan electromagnetic and radio waves. Dr. Bracewell advised researchers at Stanford and other institutions, but did not conduct laboratory research in the field.

The Fourier Transform
and Its Applications,
by Ronald Bracewell.

Russ and I cite Bracewell’s 1990 paper “Numerical Transforms” (Science, Volume 248, Pages 697–704). The central slice theorem was published in 1956 in the Australian Journal of Physics (Volume 9, Pages 198–217). Early in this career Bracewell published a lot in that journal, which is now defunct but maintains a website with free access to all the papers. Bracewell also wrote a marvelous book: The Fourier Transform and Its Applications (originally published in 1965, the revised 2nd edition is published by McGraw-Hill, New York, 1986). When writing this blog entry, I checked this book out of Kresge Library here at Oakland University. Once I opened it, I realized it is an old friend. I am sure I read this book in graduate school. It contains many pictures that allow the student to gain an intuition about the Fourier transform; an extraordinarily valuable skill to develop. The introduction states

The present work began as a pictorial guide to Fourier transforms to complement the standard lists of pairs of transforms expressed mathematically. It quickly became apparent that the commentary would far outweigh the pictorial list in value, but the pictorial dictionary of transforms is nevertheless important, for a study of the entries reinforces the intuition, and many valuable and common types of function are included which, because of their awkwardness when expressed algebraically, do not occur in other lists.

The text also does a fine job describing convolutions.

Convolution is used a lot here. Experience shows that it is a fairly tricky concept when it is presented bluntly under its integral definition, but it becomes easy if the concept of a functional is first understood.

Many of the ideas that Russ and I present in Chapter 11 of Intermediate Physics for Medicine and Biology are examined in more detail in Bracewell’s book. I recommend it as a reference to keep at your side as your plow through the mathematics of Fourier analysis.

Finally, Bracewell’s view of homework problems, as stated in his Preface to the second edition, mirrors my own.

A good problem assigned at the right stage can be extremely valuable for the student, but a good problem is hard to compose. Among the collection of supplementary problems now included at the end of the book are several that go beyond being mathematical exercises by inclusion of technical background or by asking for opinions.

The Gibbs Phenomenon

In chapter 11 of the 4th edition of Intermediate Physics for Medicine and Biology, Russ Hobbie and I discuss Fourier analysis, a fascinating but very mathematical subject. One of the most surprising results of Fourier analysis is the Gibbs phenomenon, which we describe at the end of Sec. 11.5 (Fourier Series for a Periodic Function).

Table 11.4 shows the first few coefficients for the Fourier series representing the square wave, obtained from Eq. 11.34… Figure 11.16 shows the fits for n = 3 and n = 39. As the number of terms in the fit is increased, the value of Q [measuring the least squares fit between the function at its Fourier series] decreases. However, spikes of constant height (about 18% of the amplitude of the square wave or 9% of the discontinuity in y) remain. These are seen in Fig. 11.16. These spikes appear whenever there is a discontinuity in y and are called the Gibbs phenomenon.

You have to be amazed by the Gibbs phenomenon. Think about it: as you add terms in the sum, the fit between the function and its Fourier series gets better and better, but the overshoot in amplitude does not get any smaller. Instead, the region containing ringing near the discontinuity gets narrower and narrower. If you want to see a figure like our Fig. 11.16 presented as a neat animation, take a look at http://www.sosmath.com/fourier/fourier3/gibbs.html. Also, check out http://ocw.mit.edu/ans7870/18/18.06/javademo/Gibbs/ for an interactive demo that will let you include up to 200 terms in the Fourier series.

The Gibbs phenomenon is important in medical imaging. The entry for the Gibbs phenomenon from the Encyclopedia of Medical Imaging is reproduced below.

Gibbs phenomenon, (J. Willard Gibbs, 1839-1903, American physicist), phenomenon occurring whenever a “curve” with sharp edges is subject to Fourier analysis. The Gibbs phenomenon is relevant in MR imaging, where it is responsible for so-called Gibbs artefacts. Consider a signal intensity profile across the skull, where at the edge of the brain the signal intensity changes from virtually zero to a finite value. In MR imaging the measurement process is a breakdown of such intensity profiles into their Fourier harmonics. i.e. sine and cosine functions. Representation of the profiles measured with a limited number of Fourier harmonics is imperfect, resulting in high frequency oscillations at the edges, and the image can therefore exhibit some noticeable signal intensity variations at intensity boundaries: the Gibbs phenomenon, overshoot artefacts, or “ringing.” The artefacts can be suppressed by filtering the images. However, filtering can in turn reduce spatial resolution.

Figures 12.24 and 12.25 of our book show a CT scan with ringing inside the skull and its removal by filtering, an example of the Gibbs phenomenon.

Josiah Willard Gibbs was a leading American physicist from the 19th century. He is particularly well known for his contributions to thermodynamics. Gibbs appears at several places in Intermediate Physics for Medicine and Biology. Section 3.17 discusses the Gibbs free energy, a quantity that provides a simple way to keep track of the changes in total entropy when a system is in contact with a reservoir at constant temperature and pressure. A footnote on page 68 addresses the Gibbs paradox (which deserves an entire blog entry of its own), and Problem 47 in Chapter 3 introduces the Gibbs factor (similar to the Boltzmann factor but including the chemical potential).

Selected Papers of Great
American Physicists.

The preface to Gibbs’ book on statistical mechanics is reproduced in Selected Papers of Great American Physicists: The Bicentennial Commemorative Volume of the American Physical Society 1976, edited by Spencer Weart. I recall being quite impressed by this book when in graduate school at Vanderbilt University. Below is a quote from Weart’s biographical notes about Gibbs.

Gibbs, son of a Yale professor of sacred literature, descended from a long line of New England college graduates. He studied at Yale, received his Ph.D. there in 1863—one of the first doctorates granted in the United States—tutored Latin and natural philosophy there, and then left for three decisive years in Europe. Up to that time, Gibbs had shown interest in both mathematics and engineering, which he combined in his dissertation “On the Form of the Teeth in Wheels in Spur Gearing.” The lectures he attended in Paris, Berlin and Heidelberg, given by some of the greatest men of the day, changed him once and for all. In 1871, two years after his return from Europe, he became Yale’s first Professor of Mathematical Physics. He had not yet published any papers on this subject. For nine years he held the position without pay, living on the comfortable inheritance his father had left; only when Johns Hopkins University offered Gibbs a post did Yale give him a small salary.

Gibbs never married. He lived out a calm and uneventful life in the house where he grew up, which he shared with his sisters. He was a gentle and considerate man, well-liked by those who knew him, but he tended to avoid society and was little known even in New Haven. Nor was he known to more than a few of the world’s scientists—partly because his writings were extremely compact, abstract and difficult. As one of Gibb’s European colleagues wrote, “Having once condensed a truth into a concise and very general formula, he would not think of churning out the endless succession of specific cases that were implied by the general proposition; his intelligence, like his character, was of a retiring disposition.” The Europeans paid for their failure to read Gibbs: A large part of the work they did in thermodynamics before the turn of the century could have been found already in his published work.

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

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.

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 4a 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.