Mental Calculation
By taking a lot - really, a lot -
of X-ray pictures we can in effect even see
into objects... Bill Casselman
University of British Columbia, Vancouver, Canada
cass
at math.ubc.ca
IntroductionAmong the most fascinating applications of mathematics in the modern
world are those to the science of imaging in medicine
and biological research. Some that I myself have recently
experienced first hand are the use of the Doppler effect and
ultrasound to visualize the motion of fluids inside the body,
and that of optical coherence tomography (OCT)
to view sections of the retina
from the exterior of the eye. The visual effects are
spectacular as well as informative,
and as close to wizardry as I hope to come.
In a different category is the use of similar
techniques to track the brain as it performs various tasks.
Among the tasks so far examined are those involved
in at least elementary arithmetical computations.
One paper on the subject (by Dehaene et al.) opens with this quote
from Hadamard's book:
Will it ever happen that mathematicians will know enough of that
subject of the physiology of the brain
and that neurophysiologists know enough of mathematical discovery
for efficient cooperation to be possible?
but so far the promise has not been fulfilled.
It is nonetheless curious that in this business, at least potentially,
mathematicians are both producers and consumers.
The original example of this sort of technology, and the ancestor
of many of these technologies, is what is now
called computed tomography,
for which Allan Cormack, a
physicist whose research became more and more mathematical as time went on,
laid down the theoretical foundations around 1960.
He shared the 1979 Nobel prize in medicine
for his work in this field.
In fact the basic idea of tomography had been
discovered for purely theoretical reasons in 1917
by the Austrian mathematician Johann Radon,
and it had been rediscovered several times since by others,
but Cormack was not to know this until much later than his own
independent discovery.
The problem he solved is this: Suppose we know
all the line integrals through a body of varying density.
Can we reconsruct the body itself? The answer, perhaps
surprisingly, is that we can, and furthermore we can do so
constructively. In practical terms, we know that a single
X-ray picture can give only limited information because certain
things are obscured by other, heavier things. We might take
more X-ray pictures in the hope
that we can somehow see behind the obscuring objects,
but it is not at all obvious that by taking a lot - really, a lot -
of X-ray pictures we can in effect even see
into objects, which is what Radon tells us,
at least in principle.
Making Radon's theorem into a practical tool was not
a trivial matter.
Projection
There are several different kinds of tomography. The simplest one,
which is what I'll look at here, is one in which
several beams are shot through an object,
in parallel series from several angles.
As always with the transmission of electromagnetic waves,
a ray is continually attenuated by a certain factor as it passes through
a substance. The strength I of the ray
at any point at distance s from the source satisfies a differential equation
I'(s) = -F(s) I(s)
where F(s) is the attenuation factor,
which means that in the plane of the detector the intensity I
is the exponential of the line integral of -F(s).
Radon's problem is therefore this:
Given all the line integrals
of a function f(x,y), can we determine f itself?
Explicitly, lines L in the plane can be parametrized by pairs (A, h)
where A is an angle
and h the oriented height of the line above the origin, as in this picture:
The parametrization of points on the line
is also easy to see, as shown in the figure:
(x,y) = h(- sin A, cos A) + s(cos A, sin A)
Since a line has two possible orientations,
the line with parameters (A, h) is the same
as that with (A +  , -h).
For a body with
attenuation factor f(x,y) the `X-ray' projection gives us therefore
a function g(L) and an associated function g(A, h) with g(A, h) = g(A + , -h) equal to the line integral of f along L.
The explicit formula for the projection of f is
If the function f is radially symmetric around the origin,
say f(x, y) = F(r),
its projection doesn't depend on the projection angle A. In that case
g(A, h) = g(h) and we have
For example, the projection of a uniform unit disk gives this picture:
In any case,
Radon 's problem is to find a formula for f
in terms of g. In the radially symmetric case,
this reduces directly to solving Abel's integral equation,
and this connection played an important role in
the discovery of results.
Radon's inversion formula
Radon's statement of how to recover f
from g is this:
Suppose we are given the function g(L). For a point P = (x,y)
and r > 0
let FP(r) be the average of the value of g over all lines
that are tangent to the circle of radius r centred at P. In particular
FP(0) is the average over all lines that pass through
P, which is called the back-projection of g.
Then
The formula is invariant under translations,
so it suffices to prove it for P equal to the origin.
In this case he reduces it to solving Abel's integral equation,
at least in favourable cases. But he goes on to
give a second, more enlightening derivation.
Modern proofs of the formula use the Fourier transform,
and none of them seems to me quite as natural
as Radon's original.
There is a problem even in making sense of
Radon's formula, because of the singularity 1/r
near 0. To understand
better what's going on, the modern version of the statement,
if not of the proof, is better in both theory and practice.
Starting with g(L), start with the point P; for each angle A
calculate R = R(A) to be the projection of P
perpendicular to the direction of A;
define an auxiliary function H and then continue:
The equivalence of this with Radon's formula is essentially a matter of interchanging
the order of integration.
We still have the singularity, but the integral can now be defined as a
principal value. If g is smooth enough,
the positive and negative terms in the integral
on either side cancel out. This idea manifests itself
even more concretely, since under the assumption that functions are of compact support
we can write very generally
This makes perfect sense since ln|x| is integrable near 0.
It even turns out to be practical.
Practical reconstruction
Mathematically, the process looks fairly simple,
but making the mathematics practical is not all
that simple.
First of all the projection can be obtained for only
a finite number of values of A and h.
This causes intrinsic problems, since in reconstructing
values at points (x,y) we must have access to values for other h.
The standard method used in the field
is one called Filtered Back Projection (FBP), and it uses
linear interpolation for this purpose.
If we are working with grids of
linear size N, the process takes time proportional to O(N3),
which is reasonably good but not spectacular.
It is better than might be expected, however, because due to
noise very large N are probably not feasible.
Recently, practical methods of time
O(N2 log N) have been proposed (for example in the paper
by Brandt et al.).
The F in FBP is not conceptually simple,
and I was therefore pleased to see that
recently a more direct method of inverse
computation has been proposed first in the paper by
la Rivière, and again in a paper by
Fokas, Iserles, and Marinakis who incorporate this
in working with a more complicated kind of tomography. It is
not fast (at least not yet),
being of order O(N4) in the straightforward
implementation. But since it is direct, one might hope
that errors are easier to understand.
In this direct calculation for each
angle Ai one calculates the cubic spline in h
for the function g(Ai,h), and then calculates the principal values directly
for these splines.
Incidentally, the papers I mention above
claim to use a formula for the inverse calculation of
the derivative of principal value integrals
substantially more complicated than the one
I mention above.
What can go wrong ...
This procedure is conceptually simple, but does it work? Programming the
direct method turned out to be simpler than that of other methods.
Of course I ran some simple tests
to check my work. Below, I
take a solid disk of radius 3/4, i.e. mathematically the
characteristic function of that disk, and first calculate
the projection, then the inverse projection,
and then see what there is to see.
Here are the reconstructed images I got with
N = 32, 64, and 128:
Something peculiar is going on in the backgrounds.
Now, this is not quite a straightforward
record of the data I computed.
I anticipated a few negative values in the reconstructed
function, which of course have no physical
meaning, so in drawing the picture I assigned them -
somewhat arbitrarily - to be 0. This
is of course the color of the background. Since the
background doesn't seem to be
very uniform, it looked like a good idea to investigate
the negative values more carefully. First I did a frequency analysis
of all the pixel values I calculated. Here is what the
graph looks like for N=128:
What we get appears to be a great deal of variability clustered
around 0, which I hadn't expected at all, as well
as a reasonably accurate representation of the values around 1.
These values clustered around 0
appear to be spread in lines,
so it seems that we are looking at an artefact of the discretization of
theoretically continous data. I don't know a thorough explanation
of what's going on,
but you can see with a little thought that some problem
is likely. The inverse transform,
like the original,
has to be of compact support,
but in the computation the sweep through angles A
involves accessing non-zero values. There has to be
some kind of high-powered cancellation going on,
which suggests rapid spatial oscillations. The discretization
doesn't handle this well.
Final remarks
There are several standard models for testing out ways
to analyze X-ray data. One is an artificial
model of a head with some tumours in it,
first proposed in a classic paper
by Logan and Shepp. The model itself is a collection of ellipses
in an elliptical skull of finite and variable thickness:
The light blue areas represent slices of tumors, the oblong
regions the ventricles, the rest gray matter. Here is what
my program does with this for N = 32, 64, 128:
For these images, as for many
medical problems, the main problem
is to separate the regions whose densities are spread closely
around that of water. It was therefore quite reasonable
to manipulate the data automatically
to enhance the distinctions between these densities.
In the literature there is much discussion about the thin
artificial band just inside the skull. Aside from this band,
the results seem satisfactory, especially
considering how unsophisticated this has been.
References
The literature is vast. The topic is one
involving life and death,
and among other rewards for clever work
are fame and fortune.
The following are items I have found interesting
and useful. One of the themes followed here is
not only that of mathematics
as a producer of images,
but of mathematics as a process in the brain
that can be tracked with imaging, in this
case with emission tomography and functional
magnetic resonance imaging (fMRI)
-
J. R. Anderson, Y. Qin, M.-H. Sohn,
V. A. Stenger, and C. Carter,
`An information-processing model
of the BOLD response in symbol manipulation tasks',
Psychonomic Bulletin and Review 10 (2003), 241-261.
Raises the intriguing possibility that some day we will be able to follow
in detail the thought
processes involved in doing mathematics.
-
A. Brandt, J. Mann, M. Brodski,
and M. Galun,
`A fast and accurate multilevel
inversion of
the Radon transform',
Siam Journal of Applied Mathematics 60 (1999),
437-462. This is the latest in a long series of proposals
to replace the standard technology by something much
faster, here of time O(N2 log N).
-
Allan M. Cormack,
`Early two-dimensional reconstruction and recent topics stemming from it',
Nobel Prize lecture, 1979. Available at
http://nobelprize.org/medicine/laureates/1979/cormack-lecture.html.
An enjoyable account of early days in tomography.
-
S. Dehaene, E. Spelke, P. Pinel,
R. Stanescu, S. Tsivkin,
`Sources of mathematical thinking',
Science 284 (7 May 1999), 970-974.
One of the first reports of attempts to see (literally) what
goes on in the process of doing ... well, the authors call it
mathematics, but to a mathematician it
will look rather elementary.
-
A. S. Fokas and L.-Y. Sung,
`Generalized Fourier transforms, their nonlinearization
and the imagine of the brain',
to appear in the NOTICES of the American Mathematical Society,
December, 2005.
-
A. S. Fokas, A. Iserles, and V. Marinakis,
`Reconstruction algorithm for single
photon emission computed tomography
and its numerical implementation',
Journal of the Royal Society Interface,
http://www.pubs.royalsoc.ac.uk/interface.shtml,
DOI: 10.1098/rsif.2005.0061. This is
where I found the idea I interpret here,
although I think it is fair to
say that in so far as classical
tomography is concerned most credit goes to
la Rivière and Pan.
But this paper also contains
an explicit algorithm for inversion of emission tomography,
which up to now has been
a notoriously inaccessible process.
This is the technology which traces brain activity by reading output
from a chemical inside an object that's stimulated by
activity of some kind - such as mathematical thought.
Jacques Hadamard ,
An essay on the psychology of invention in
the mathematical field,
Princeton University Press, 1945.
Reprinted by Dover.
-
Sigurdur Helgason,
The Radon Transform,
Birkhaüser, 1980.
-
F. Natterer,
`Numerical methods in tomography',
Acta Numerica (1999), 1-35. Also found at
http://wwwmath1.uni-muenster.de/num/Preprints/1998/natterer_1/.
Standard summary of the state of the art.
-
Ivars Peterson,
`Brainy figuring', in Math Trek (May 31 1999),
http://www.maa.org/mathland/mathtrek_5_31_99.html.
Tony Phillips, `The mathematical CAT scan',
the December 1999 column in the
archive
for this monthly AMS feature. A very different but also very
enjoyable look at the
discrete Radon transform.
-
W. H. Press et al., Numerical Recipes, Cambridge University Press, (various editions).
Section 3.3 is on cubic spline interpolation, and this can be found on the Internet at
Cornell's NR web site.
-
Johann Radon,
`Über die Bestimmung von Funktionen
durch ihre Integralwerte längs gewisser Mannigfaltigkeiten',
Berichte Sächsische Akademie der Wissenschaften
Math.-Phys. Kl. 69 (1917), 262-267.
Reprinted in facsimile as an appendix to Helgason's book.
A surprisingly short and uninformative
biography of Radon can be found at
the St. Andrews' history site. It does not mention the Radon transform,
which is surely his most famous contribution to mathematics,
if not the most difficult.
-
P. J. la Rivière and X. Pan,
`Spline-based inverse Radon transform
in two and three dimensions',
IEEE Transactions on Nuclear Science 45 (1998), 2224-2231.
This is the first place I am aware of that suggests using
cubic splines for interpolation in
Radon's formula.
-
Lawrence Shepp, a lecture on tomography, at
http://www.pitt.edu/~super1/lecture/lec13611/001.htm.
-
L. A. Shepp and B. F. Logan,
`The Fourier reconstruction of a head section',
IEEE Transactions on Nuclear Science NS-21 (June 1974),
21-43. A classic.
Bill Casselman
University of British Columbia, Vancouver, Canada
cass
at math.ubc.ca
NOTE: Those who can access JSTOR can find some of the
papers mentioned above there. For those with access, the American Mathematical
Society's MathSciNet can be used to get
additional bibliographic information and reviews of some these materials. Some of the
items above can be accessed via the ACM
Portal, which also provides bibliographic services. |
