Projections and Back Projections

Tomography is one of the most important contributions of mathematics to medicine. In the 4th edition of Intermediate Physics for Medicine and Biology, Russ Hobbie and I describe two methods to solve the problem of tomographic reconstruction.

The reconstruction problem can be stated as follows. A function f(x,y) exists in two dimensions. Measurements are made that give projections: the integrals of f(x,y) along various lines as a function of displacement perpendicular to each line…F(θ,x'), where x' is the distance along the axis at angle θ with the x axis. The problem is to reconstruct f(x,y) from the set of functions F(θ,x'). Several different techniques can be used… We will consider two of these techniques: reconstruction by Fourier transform […] and filtered back projection…. The projection at angle θ is integrated along the line y':

[where x = x' cosθ – y' sinθ and y = x' sinθ + y' cosθ]… The definition of the back projection is

where x' is determined for each projection using Eq. 12.27 [x' = x cosθ + y sinθ].

In IPMB, Homework Problem 32 asks you can take the function (the “object”)

and calculate the projection using Eq. 12.29, and then calculate the back projection using Eq. 12.30. The object and the back projection are different. The moral of the story is that you cannot solve the tomography problem by back-projection alone. Before you back-project, you must filter.*

I like having two homework problems that illustrate the same point, one that I can do in class and another that I can assign to the students. IPMB contains only one example of projecting and then back-projecting, but recently I have found another. So, dear reader, here is a new homework problem; do this one in class, and then assign Problem 32 as homework.

Problem 32 ½ Consider the object f(x,y) = √(a² − x² − y²)/a for √(x² +y²) less than a, and 0 otherwise.

(a) Plot f(x,0) vs. x 

(b) Calculate the projection F(θ,x'). Plot F(0,x') vs. x'.

(c) Use the projection from part (b) to calculate the back projection f_b(x,y). Plot f_b(x,0) vs. x.

(d) Compare the object and the back projection. Explain qualitatively how they differ.

The nice thing about this function

(as well as the function in Problem 32) is that f(x,y) does not depend on direction, so F(θ,x') is independent of θ; you can make you life easier and solve it for θ = 0. Similarly, you can calculate the back projection along any line through the origin, such as y = 0. I won’t solve Problem 32½ here in detail, but let me outline the solution.

Below is a plot of f(x,0) versus x

To take the projection in part (b) use θ = 0, so x' = x and y' = y. If |x| is greater than a, you integrate over a function that is always zero, so F(0,x') = 0. If |x| is less than a, you must do more work. The limits of the integral over y become −√(a² – x²) to +√(a² – x²). The integral is not too difficult (it’s in my mathematical handbook), and you get F(0,x) = π(a² – x²)/2a. Because the projection is independent of θ, you can generalize to

The plot, an inverted parabola, is

In part (c) you need to find the back-projection. I suggest calculating it for the line y = 0. Once you have it, you can find it for any point by replacing x by √(x² +y²). The back-projection for |x| less than a is easy. The integral in Eq. 12.30 gives another inverted parabola. For |x| is greater than a, the calculation is more complicated because some angles give zero and some don’t. A little geometry will convince you that the integral should range from cos^-1(a/x) to π - cos^-1(a/x). Because the function is even around π/2, you can make your life easier by multiplying by two and integrating from cos^-1(a/x) to π/2. The only way I know to show how you get cos^-1(a/x) is to draw a picture of a line through the point x greater than a, y = 0 that is tangent to the circle x² + y² = a², and then do some trigonometry. When you then evaluate the integral you get

A plot of this complicated looking function is

To answer part (d), compare your plots in (a) and (c). The object in (a) is confined entirely inside the circle x² + y² = a², but the back-projection in (c) spreads over all space. One could say the back-projection produced a smeared-out version of the object. Why are they not the same? We didn’t filter before we back projected.

If you really want to have fun, show that the limit of the back-projection goes to zero as x goes to infinity. This is surprisingly difficult—that last term doesn’t look like it’ll go to zero—and you need to keep more than one term in your Taylor series to make it work.

The back projection shown above looks similar to the back projection shown in Fig. 12.23 of IPMB. My original goal was to calculate the back projection in Fig. 12.23 exactly, but I got stuck trying to evaluate the integral

I sometimes ask myself: why do I assign these analytical examples of projections and back-projections? Anyone solving the tomography problem in the clinic will certainly use a computer. Is anything gained by doing analytical calculations? Yes. The goal in doing these problems is to gain insight. A computer program is a black box: you put projections in and out comes an image, but you don’t know what happens in between. Analytical calculations force you to work through each step. Please don’t skip the plots in parts (a) and (c), and the comparison of the plots in part (d), otherwise you defeat the purpose; all you have done is an exercise in calculus and you learn nothing about tomography.

I wish I had come up with this example six months ago, so Russ and I could have included it in the 5th edition of IPMB. But it’s too late, as the page proofs have already been corrected and the book will be printed soon. This new homework problem will have to wait for the 6th edition!


*Incidentally, Prob. 12.32 has a typo in the second line: "|x|" should really be "√(x² +y²)".

Consider an Impervious Plane Containing A Circular Disk

When I teach, one of my goals is to help students analyze a mathematical equation to uncover the underlying physics. Equations are not merely expressions you plug numbers into to get other numbers; equations tell a story. Students often begin the semester knowing the required math and physics, but have difficulty connecting the two. One way students learn this connection is by taking limits of mathematical equations, so they can show that the equation behaves as it should when some key variable is large or small. If the equation is too complicated, taking the limit may be so difficult that the student gets lost in the algebra (see the February 20 blog entry about “the ugliest equation”). If the equation is too simple, taking the limit may be trivial. The best examples are of intermediate difficulty: students must do some work to take the limit, but their result provides physical insight.

One of my favorite examples is Problem 29 of Chapter 4 in the 4th edition of Intermediate Physics for Medicine and Biology. I should be more specific: parts (b)-(d) of that problem are just the sort of practice the students need (part (a) of the problem is for masochists).

Problem 29 Consider an impervious plane at z = 0 containing a circular disk of radius a having a concentration C₀. The concentration at large z goes to zero. Carslaw and Jaeger (1959) show that the steady-state solution to the diffusion equation is 

(a) (optional) Verify that C(r, z) satisfies ∇ ²C = 0. The calculation is quite involved, and you may wish to use a computer algebra program such as Mathematica or Maple.

(b) Show that for z = 0, C = C₀ if r is less than a. 

(c) Show that for z = 0, dC/dz = 0 if r is greater than a. 

(d) Integrate J_z over the disk (z = 0, r between 0 and a) and show that i = 4DaC₀.

I’m amazed that the convoluted expression for the concentration, C(r, z), containing an inverse sine function with a complicated argument does indeed obey Laplace’s equation. I once proved that this is true, but have no intention of doing so again—it’s a mess. Before analyzing parts (b)-(d) of the problem, let’s examine the claim that “the concentration at large z goes to zero.” Exactly how does it go to zero? Far from the disc (z is much greater than a) the concentration should fall off as the inverse distance (just as the electrical potential of a point charge falls off as 1/r). When z is large, the argument of the inverse sine function is small, and it can be approximated using the first term of its Taylor series: sin^-1(x) ~ x. In that case, the concentration becomes 2aC₀/πz, falling off with the inverse distance as expected.

The behavior at small z is more puzzling and the solution has the astonishing property that it changes its behavior at r = a. For r less than a the concentration is constant, but for r greater than a the derivative of the concentration is zero. How can such a relatively simple analytical expression display such a discontinuous behavior?

To answer that question, let’s analyze the equation for the concentration. The key step in our analysis is when we factor out (r – a)² from the first square root in the denominator of the argument. This factor is always positive because the quantity is squared. But when removed from inside the square root, it becomes |r – a|, with the absolute value needed to ensure that the expression remains positive. There are two cases: r less than a, when we replace |r – a| by (a – r); and r greater than a, when we replace |r – a| by (r – a). The rest is an exercise in Taylor series resulting from the limit as z goes to zero. Use the expansion √(1 + x²) = 1 + x²/2 for the two square roots, simplify the denominator as much as possible, and then use the expansion 1/(1 - x) = 1 + x to find that for r less than a

and for r greater than a

For all the gory details see the solution manual, available free-of-charge to instructors using IPMB in their class; email Russ Hobbie or me (roth@oakland.edu) for a copy.

Now we can do parts (b)-(d) of the homework problem. For part (b), when r is less than a set z = 0 and use that sin⁻¹(1) = π/2 to find that C(r,0) = C₀, independent of r. For part (c), when r is greater than a take the derivative of the inverse sine function with respect to z, and then set z = 0. Because the argument of the inverse sine is even in z, at z = 0 the derivative dC/dz vanishes.

Part (d) is tricky. Using the same argument as above, we might expect that because the argument of the inverse sine is even for r less than a the derivative dC/dz must be zero. But no! The derivative of sin⁻¹(u(x)) is 1/√(1−u²) du/dx. When we take the derivative of the expression for C(r,z) for r less than a and then set z = 0, we get not only zero in the numerator (because du/dx = 0; the function is even) but also zero in the denominator (because √(1−u²) = 0 in this case). The result is dC/dz = − 2C₀/π√(a² − r²).

To find the total number of particles per unit time i coming out of the disk, we need to integrate the particle current density J_z over the plane z = 0. The current density is J_z = – D dC/dz, where D is the diffusion constant. Because dC/dz is zero for r greater than a, we need to integrate the current density over only the circular area r less than a. You can look this integral up and find that i = 4C₀Da.

There are other aspects of the solution that are not emphasized in this problem, but are equally interesting. For instance, J_z goes to infinity at r = a. In other words, the number of particles leaving the disk is largest at the edge of the disk. This result is analogous to the concentration of electrical current near the edge of a disk electrode, discussed in Problem 22 of Chapter 6 in IPMB; J_z changes abruptly from an extremely large value for r just less than a to zero for r just greater than a. Another interesting result is that the total number of particles leaving the disk per unit time is proportional to the disk radius, not the disk area as you might naively expect. This result also has an analogy in electrostatics: the capacitance of a disk electrode relative to a distant ground is proportional to the radius of the disk, not the radius squared, which surprised me when I first encountered it (see Section 6.19 in IPMB).

In retrospect, maybe this problem is slightly too complicated to be the perfect example of taking the limits of a mathematical expression to gain insight into the physical behavior. But for those whose mathematical skills are up to the job, it provides an excellent case study.

Mathematical equations tell a physical story. You must practice extracting that story from the math.

If You Want Healthy Cows Feed Them Magnets

I saw an article on the internet claiming “If you want healthy cows feed them magnets” and I thought “oh no, not more biomagnetism nonsense.” First magnets in shoes to relieve foot pain, then magnetic bracelets for arthritis, and finally “biomagnetic therapy” for all sorts of disorders; I thought it couldn’t get worse, but feeding magnets to heifers? Really? Sounds like bull to me.

A drawing of a cow
with a magnet in its stomach.

Well, I can’t vouch for the accuracy of this story or the effectiveness of the treatment, but at least the mechanism underlying the feeding of magnets to cows is plausible. Cattle swallow a lot of junk while eating, including some that is magnetic (for example, wires and nails...yikes!). The article says

That's where magnets come in. A magnet about the size and shape of a finger is placed inside a bolus gun, essentially a long tube that ensures the magnet goes down the cow's throat. Then it settles in the reticulum, collecting any stray pieces of metal. The magnets, which cost a few bucks a pop, can also be placed preventatively. To check if a cow already has a magnet, farmers use a compass.

Apparently the “bolus gun” is inserted through the mouth; I wasn’t so sure. Wikipedia has a page about cow magnets, titled “hardware disease.” Companies make money selling cow magnets (these are big magnets, about four inches long). But even though calves eat magnets, kids should not (note the plural: the problems arise when magnets interact).

Consider a spherical cow.

The 4th edition of Intermediate Physics for Medicine and Biology has an entire chapter about biomagnetism, but no mention of magnets in bovine stomachs. What is wrong with Russ and me? The only place we mention cattle at all is in Homework Problem 30 in Chapter 4, where we analyze the temperature distribution throughout a spherical cow. A small-scale analogy of magnets in steers’ stomachs are rows of magnetosomes in magnetotactic bacteria (see Fig. 8.25 in IPMB), but I doubt the bacteria use them to collect nails before they can puncture their membrane. Yet, could we misunderstand the biological purpose of magnetosomes?

Finally, I have some good news and bad news about the 5th edition of IPMB. The good news: we submitted the page proofs and the book should be published in the next few months. The bad news: no more mention of livestock in the revised edition.

If you want healthy cows feed them magnets.

A Mathematical Model of Agonist-Induced Propagation of Calcium Waves in Astrocytes

When I was working at the National Institutes of Health in the mid-1990s, I spent most of my time studying transcranial magnetic stimulation and theoretical cardiac electrophysiology. But also I collaborated with James Russell to study calcium waves in astrocytes (a type of glial cell in the brain), and we published a paper in the journal Cell Calcium describing “A Mathematical Model of Agonist-Induced Propagation of Calcium Waves in Astrocytes” (Volume 17, Pages 53–64, 1995). The introduction is reproduced below (with citations removed):

Recent experiments have clearly shown that astroglia in brain participate in long distance signaling together with neurons. Such signalling in astrocytes is supported by intracellular calcium oscillations induced by neurotransmitters that are propagated as waves through the cytoplasm of individual cells and through astrocyte networks. These calcium oscillations generally are triggered by activation of metabotropic receptors which are coupled to inositol 1,4,5-trisphosphate (IP₃) generation and intracellular calcium release through IP₃-gated calcium channels on the endoplasmic reticulum (ER) membrane. Yagodin et al. have shown that, in astrocytes, wave propagation is saltatory, with discrete loci of nonlinear wave amplification separated by regions through which passive diffusion of calcium occurs. These wave amplification loci appear to be intracellular specializations that remain invariant and support a qualitatively characteristic response pattern in any given cell. The loci may each have different intrinsic oscillatory frequencies, resulting in complex spatio-temporal dynamics, with wave collisions and annihilations.

Several mathematical models have been presented that describe the temporal characteristics of agonist-induced calcium oscillations in different types of cells. A few of these models address the spatial characteristics of wave propagation, but none have addressed the complex wave dynamics observed in different types of cells including astrocytes. The purpose of this paper is to extend a previously developed model of calcium oscillations so that it includes spatial diffusion of calcium in a cell with discrete active loci of wave amplification. This model is then used to analyze experimental data and to gain insight into the mechanism of wave collisions and annihilations.

As you might expect, my contribution to this paper was developing the mathematical model, while Russell and his team provided the experimental data as well as the biological knowledge and insight. The model was based on a paper by Li and Rinzel (“Equations for InsP3 Receptor-Mediated [Ca2+]i Oscillations Derived from a Detailed Kinetic Model: A Hodgkin-Huxley Like Formalism,” Journal of Theoretical Biology, Volume 166, Pages 461–473, 1994). At that time, John Rinzel was at NIH, heading the Mathematical Research Branch of the National Institute of Diabetes and Digestive and Kidney Diseases. John, now at the Center for Neural Science at New York University, has contributed much to theoretical biology, but I remember him best for his work on bursting of pancreatic beta cells. He is this year’s winner of the Society of Mathematical Biology’s Arthur T. Winfree Prize “for his elegant work on the analysis of dynamical behavior in models of neural activity and the contributions that work has made in the neurobiological community to the understanding of a host of phenomena (including simple excitability as well as bursting) in single neurons, small populations of neurons, and other excitable cells.”

Russell remains at NIH with the Microscopy and Imaging Core of the Eunice Kennedy Shriver National Institute of Child Health and Development. He leads a multi-user research facility providing training and instrumentation for high resolution microscopy and image processing.

As so often happens, an echo of my work on calcium wave modeling with Russell appears in the 4th edition of Intermediate Physics for Medicine and Biology. Homework Problem 24 in Chapter 4 contains a simplified model of calcium waves. This system is a classic “reaction-diffusion” system: calcium diffuses down the cell, triggering calcium-induced calcium release, which produces more diffusion, triggering more calcium release, resulting in positive feedback and a calcium wave. The process is analogous to action potential propagation along a nerve.