Page 199

A Cross-Cutting Look at the Mathematics of Emerging Biomedical Imaging

This chapter describes the mathematical models used in medical imaging, their limitations, and the pertinent mathematical methods and problems. In many cases these problems have not yet been solved satisfactorily.

This description of mathematical models is restricted to those features that are important for the mathematical scientist. More about the physical details, the assumptions made, and the limitations can be found in Chapters 2-11.

Transmission computed tomography (CT) is the original and simplest case of CT. In transmission tomography one probes an object with non-diffracting radiation, e.g., x-rays for the human body. If *I* _{0} is the intensity of the source, *a(x)* the linear attenuation coefficient of the object at point *x, L* the ray along which the radiation propagates, and *I* the intensity past the

Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.

Intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.

OCR for page 199

Page 199
Chapter 14 A Cross-Cutting Look at the Mathematics of Emerging Biomedical Imaging This chapter describes the mathematical models used in medical imaging, their limitations, and the pertinent mathematical methods and problems. In many cases these problems have not yet been solved satisfactorily.
14.1 Mathematical Models for Particular Imaging Modalities This description of mathematical models is restricted to those features that are important for the mathematical scientist. More about the physical details, the assumptions made, and the limitations can be found in Chapters 2-11.
14.1.1 Transmission Computed Tomography Transmission computed tomography (CT) is the original and simplest case of CT. In transmission tomography one probes an object with non-diffracting radiation, e.g., x-rays for the human body. If I 0 is the intensity of the source, a(x) the linear attenuation coefficient of the object at point x, L the ray along which the radiation propagates, and I the intensity past the

OCR for page 199

Page 200
object, then
(compare to section 3.1.2). In the simplest case the ray L may be thought of as a straight line. Modeling L as a strip or cone, possibly with a weight factor to account for detector inhomogeneities, may be more appropriate. Equation 14.1 neglects the dependence of a with the energy (beam hardening effect) and other nonlinear phenomena (e.g., partial volume effect); see section 3.2.3.
The mathematical problem in transmission tomography is to determine a from measurements of I for a large set of rays L. If L is simply the straight line connecting the source x0 with the detector x1, equation 14.1 gives rise to the integral
where dx is the restriction to L of the Lebesgue measure in Rn. The task is to compute a in a domain.W Í Rn.from the values of equation 14.2 where x0 and x1 run through certain subjects of ¶W.
For n = 2, equation 14.2 is simply a reparametrization of the Radon transform R. The operator R is defined to be
where is a unit vector in Rn ; i.e., q Î Sn-1 , and s Î R. Thus a is in principle found through Radon's inversion formula for R,
R* is given explicitly by
and the operator K is given by
where H is the Hilbert transform. In fact the numerical implementation of equation 14.4 leads to the filtered backprojection algorithm, which is the standard algorithm in commercial CT scanners.

OCR for page 199

Page 201
For n = 3, the relevant integral transform is the x-ray transform, P, given by
where q Î Sn-1 and x Î q^. P admits a similar inversion formula as R :
with K very similar to equation 14.6 and
where Eq is the orthogonal projection onto q^. Unfortunately, equation 14.7 is not as useful as equation 14.4. The reason is that equation 14.7 requires g for all q and y Î q^; i.e., equation 14.2 has to be available for all x0, xl Î ¶W. This is not practical. Also, it is not necessary for unique reconstruction of a. In fact it can be shown that a can be recovered uniquely from equation 14.2 with sources x0 on a circle surrounding supp(a) and x1 Î Sn-1. Unfortunately, as section 14.1.3 makes clear, the determination of a in such an arrangement is, though uniquely possible, highly unstable. The condition of stability is the following: each plane meeting supp(a) must contain at least one source. This condition is obviously violated for sources on a circle. Cases in which the condition is satisfied include the helix and a pair of orthogonal circles. A variety of inversion formulas have been derived for these cases.
If scatter is to be included, a transport model is more appropriate. Let u(x, q ) be the density of the particles at x traveling (with speed 1) in direction q. Then,
Here, h(x, q , q') is the probability that a particle at x traveling in direction q is scattered in direction q'. Again we neglect dependence on energy. d is the Dirac d-function modeling a source of unit strength. Equation 14.8 holds in a domain W of Rn (n = 2 or 3), and x0 Î ¶W. Since no radiation comes in from outside we have
where vx is the exterior normal on ¶W at x Î ¶W. Equation 14.1 is now replaced by

OCR for page 199

Page 202
The problem of recovering a from equations 14.8-14.10 is much harder. An explicit formula for a such as equation 14.4 has not been found and is unlikely to exist. Nevertheless, numerical methods have been developed for special choices of h; see the papers by Anikonov et al. and Bondarenko and Antyufeev listed in the suggested reading at the end of this chapter. The situation gets even more difficult if one takes into account that, strictly speaking, his object dependent and hence not known in advance. Equations 14.8-14.10 are a typical example of an inverse problem for a partial differential equation. In an inverse problem one has to determine the differential equation—in this case a, h—from information about the solution—in this case equation 14.10. More on inverse problems can be found in section 14.3.
14.1.2 Emission Computed Tomography In emission computed tomography one determines the distribution f of radiating sources in the interior of an object by measuring the radiation outside the object in a tomographic fashion. Let u(x,q ) again be the density of particles at x traveling in direction qwith speed 1, and let a be the attenuation distribution of the object. (This is the quantity that is sought in transmission CT.) Then,
This equation holds in the object region W Í Rn for each q Î Sn-1. Again there exists no incoming radiation; i.e.,
while the outgoing radiation
is measured and hence known. Equations 14.11-14.13 again constitute an inverse problem for a transport equation. If a is known, 14.11-14.12 are readily solved to yield
The exterior integral is over the ray x - tq , 0 < t < ¥. Thus, equation 14.13 leads to the integral equation

OCR for page 199

Page 203
for f; compare to section 5.4. Apart from the exponential factor, equation 14.14 is identical—up to notation—to the integral equation in transmission CT. Except for very special cases—e.g., a constant in a known domain (see works by Tretiak and Metz and by Palamodov in the suggested reading list)—no explicit inversion formulas are available. Numerical techniques have been developed but are considered to be slow. Again the situation becomes worse if scatter is taken into account. This can be done by simply adding the scattering integral in equation 14.8 to the right-hand side of equation 14.11.
What we have described so far is called single photon emission CT (SPECT). In positron emission tomography (PET) the sources eject the particles pairwise in opposite directions. They are detected in coincidence mode; i.e., only events with two particles arriving at opposite detectors at the same time are counted. Equation 14.14 has to be replaced (see section 6.1.3) by
Thus PET is even closer to the case of transmission CT. If a is known, we simply have to invert the x-ray transform. Inversion formulas are available that can make use of the data collected in PET, and their numerical implementations are currently under consideration.
The main problems in emission CT are (compare to sections 5.4.2 and 6.2.4) unknown attenuation, noise, and scatter. For the attenuation problem, the ideal mathematical solution would be a method for determining f and a in equations 14.11-14.13 simultaneously. Under strong assumptions on a—a constant in a known region (Hertle's paper in section 14.8), a an affine distortion of a prototype (Natterer's 1993 paper on tissue attenuation in section 14.8), or a close to a known distribution (Bronnikov's paper in section 14.8)—encouraging results have been obtained. Theoretical results based on the transport formulation have been derived, even for models including scatter (see Romanov's work listed in section 14.8). But a clinically useful way of determining a from the emission data has not yet been found.
Noise and scatter are stochastic phenomena. Thus, in addition to models using integral equations, stochastic models have been set up for emission tomography. These models are completely discrete. The reconstruction region

OCR for page 199

Page 204
is subdivided into m pixels or voxels. The number of events in pixel/voxel j is a Poisson random variable j whose mathematical expectation f j = Ej j is a measure for the activity in pixel/voxel j. The vector f = (f 1,. . ., f n ) is the sought-for quantity. The vector g = (g1,. . . ,gn) of measurements is considered a realization of the random variable g = (g1,. . . ,gn), where gi is the number of events detected in detector i. The model is determined by the (n,m)-matrix A = (aij ) whose elements are
aij = P(event in pixel/voxel j is detected in detector i),
where P denotes probability. Normalizing the aij such that = 1 leads to E(g) = Af. f is determined from g by the maximum likelihood method. A numerical method for doing this is the expectation maximization (EM) algorithm. In its basic form it reads (compare section 6.2.4)
where division and multiplication are to be understood componentwise. The problem with the EM-algorithm is that it is only semi-convergent (compare section 14.4); i.e., noise is amplified at high iteration numbers. This is known as the checkerboard effect. Various suggestions have been made for getting rid of this effect. The most exciting and interesting ones use prior information and attempt to maximize "posterior likelihood." Thus f is assumed to have a prior probability distribution, called a Gibbs-Markov random field p( f ), which gives preference to certain functions f. Most priors p simply add a penalty term to the likelihood function to account for correction between neighboring pixels and do not use biological information. However, if p is carefully chosen so that piecewise constant functions f with smooth boundaries forming the region of constancy are preferred, then the noise amplification at high iteration numbers can be avoided. The question remains as to whether this conclusion will remain valid for functions f that are assigned low probability by p—or, more to the point, whether "real" emission densities f will be well-resolved by this Bayesian method. A double-blind study in which radiologists tried to find lesions from images produced by two different algorithms concluded that maximum likelihood methods were superior to the filtered backprojection algorithm in certain clinical applications. The same type of study is needed to determine whether or not Gibbs priors will improve the maximum likelihood reconstruction (stopped short of convergence to avoid noise amplification) on real data.

OCR for page 199

Page 205
14.1.3 Ultrasound Computed Tomography X-rays travel along straight lines. For other sources of radiation, such as ultrasound and microwaves, this is not the case. The paths are not straight, and their exact shape depends on the internal structure of the object. Simple projections and linear integral equations will not suffice, and more sophisticated nonlinear models have to be used.
In the following we consider an object W Í Rn with refractive index g. Assume g= 1 outside the object. The object is probed by a plane wave
with wave number , l the wave length, traveling in the direction q. The resulting wave e-iktu(x) satisfies the reduced wave equation
plus suitable boundary conditions at infinity. The inverse problem to be solved is now the following. Assume that
is known outside W. Determine f inside W!
The uniqueness and stability of the inverse problem equations 14.16 14.17 have recently been settled. However, stability is only logarithmic; i.e., a data error of size d results in a reconstruction error of 1/log(1/d ) (see section 14.4). Numerical algorithms did not emerge from this work.
Numerical methods for equations 14.16-14.17 are based mostly on linearizations, such as the Born and Rytov approximations. In order to derive the Born approximation, one rewrites 14.16 as
where G is an appropriate Green's function. For n = 3, the Green's function is
The Born approximation is now obtained by assuming u ~ uq in the integral in equation 14.18. With this approximation, equation 14.17 reads

OCR for page 199

Page 206
This is a linear integral equation for f, valid for all x outside the object and for all measured directions q.
Numerical methods based on equation 14.20—and a similar equation for the Rytov approximation—have become known as diffraction tomography. Unfortunately, the assumptions underlying the Born and Rytov approximations are not satisfied in medical imaging. Thus, the reconstructions of f obtained from equation 14.20 are very poor. However, we may use equation 14.20 to get some encouraging information about stability. For | x | large, equation 14.20 assumes the form
with the Fourier transform of f . Equation 14.21 determines within a ball of radius from the data in equation 14.17 in a completely stable way. Therefore, the stability of the inverse problem 14.16-14.17 is much better than logarithmic. If the resolution is restricted to spatial frequencies below —which is perfectly reasonable from a physical point of view—then we can expect equations 14.16-14.17 to be perfectly stable.
So far we have considered plane wave irradiation at fixed frequency, and have worked in the frequency domain. Time domain methods are conceivable as well. We start with the wave equation
with the propagation speed c assumed to be 1 outside the object. With x0 a source outside the object, consider the initial conditions
We want to determine c inside the object from knowing
for many sources x0 and receivers x1 outside the object. In the one-dimensional case, the inverse problem described by equations 14.22-14.24 can be solved by the famous Gelfand-Levitan method in a stable way. It is not clear

OCR for page 199

Page 207
how Gelfand-Levitan can be extended to dimensions two and three. The standard methods use sources and receivers on all of the boundary of the object. This is not practical in medical imaging. However, for reduced data sets, comparable to those in three-dimensional x-ray tomography (compare equation 14.1), it is not known how to use Gelfand-Levitan, nor is anything known about stability.
Of course one can always solve the nonlinear problem of equations 14.2214.24 by a Newton-type method. Such methods have been developed (see papers by Gutman and Klibanov and by Kleinman and van den Berg in section 14.8). They suffer from the requirement for excessive computing time and from their apparent inability to handle large wave numbers k.
14.1.4 Optical Tomography With optical tomography one uses near-infrared lasers for the illumination of the body (see Chapter 11). The process is described by the transport equation
for the density u(x, q ,t) of the photons at x Î W flying in direction q Î Sn-1 at time t. a and b are the sought-for tissue parameters. In the simplest case the scattering kernel h is assumed to be known. The source term f is under the control of the experimenter. Together with the initial and boundary conditions
equation 14.25 has a unique solution under natural conditions on a, b, h, and f. As in equation 14.1, we pose the inverse problem: assume that the outward radiation is known,
can we determine one or both of the quantities a, b?

OCR for page 199

Page 208
There are essentially three methods for illuminating the object, i.e., for choosing the source term f in equation 14.25. In the stationary case one posits f = d(x - x0), where x0 Î ¶W is a source point. u is considered stationary, too. A second possibility is the light flash f = d(x - x0)d(t). Finally, one can also use time harmonic illumination, in which case f = d(x - x0) eiw t. This case reduces to the stationary case with a replaced by a + iw . In all three cases, the data function g of equation 14.27 is measured at x Î ¶W, possibly averaged over one or both of the variables q, t.
Light tomography is essentially an absorption and scattering phenomenon. This means that the scattering integral in equation 14.25 is essential; it can no longer be treated merely as a perturbation, as in x-ray CT. Thus the mathematical analysis and the numerical methods are expected to be quite different from what is seen in other types of tomography.
The mathematical theory of the inverse problem 14.25-14.27 is in a deplorable state. There exist some Russian papers on uniqueness (see Anikonov in section 14.8). General methods have been developed, too, but apparently they have been applied to one-dimensional problems only (see Sanchez and McCormick in section 14.8). Nothing seems to be known about stability. The numerical methods that have become known are of the Newton type, either applied directly to the transport equation or to the so-called diffusion approximation (e.g., Arridge et al., listed in section 14.8). The diffusion approximation is an approximation to the transport equation by a parabolic differential equation. Since inverse problems for parabolic equations are severely ill posed, this approach is questionable. Higher-order approximations (see, e.g., Gratton et al., in section 14.8) are hyperbolic, making the inverse problem much more stable.
As an alternative to the transport equation one can also model optical tomography by a discrete stochastic model. We consider only a very simple model of the two-dimensional case. Break up the object into a rectangular arrangement of pixels labeled by indices i, j with a £ i £ b and c £ j £ d. Attach to each pixel the quantities f ij , bij , rij , eij denoting the probability (respectively) of a forward, backward, rightward, or leftward transition out of the pixel {i, j} with respect to the direction used to get into this pixel. For each pair of boundary pixels {i, j} and {i', j'}, let Pij ,i' j' be the probability that a particle that enters the object at pixel {i, j} will eventually leave the object at pixel {i', j'}. The problem is to determine the quantities fij , bij , rij , lij from the values of Pij ,i' j' for all boundary pixels. Preliminary

OCR for page 199

Page 209
numerical tests show that this is possible, at least in principle. However, the computations are very time-consuming. More seriously, they reveal a very high degree of instability.
14.1.5 Electrical Impedance Tomography Here, the sought-for quantity is the electrical impedance s of an object W. Voltages are applied via electrodes on ¶W, and the resulting currents at these electrodes are measured. With u the potential in W, we have
Knowing many voltage-current pairs (g, f) on ¶W, the task is to determine s from equation 14.28.
Uniqueness for the inverse problem 14.28 has recently been settled. Unfortunately, the stability properties are very bad. Numerical methods based on Newton's method, linearization, simple backprojection, and layer stripping have been tried. All these methods suffer from the severe ill-posedness of the problem. There seems to be no way to improve stability by purely mathematical means.
14.1.6 Magnetic Resonance Imaging The physical phenomenon exploited in magnetic resonance imaging (MRI) is the precession of the spin of a proton in a magnetic field of strength H aboutthe direction of that field. The frequency of this precession is the Larmor frequency gH where g is the gyromagnetic ratio. By making the magnetic field H space-dependent in a controlled way, the local magnetization M0(x) (together with the relaxation times T 1(x) and T 2(x)) can be imaged. In the following we derive the imaging equations; compare to section 4.1.
The magnetization M(x,t) caused by a magnetic field H(x,t) satisfies the Bloch equation
Here, Mi and Mi0 are the i-th components of M and M0, respectively, and ei is the i-th unit vector i = 1,2,3. The significance of T 1, T 2, and M0

OCR for page 199

Page 220
14.4.4 Regularization by Discretization Let Ah be a discrete approximation to A with step size A. For h ® 0 we assume Ah ® A in some sense. Since A-1 is not continuous, (if it exists) tends to infinity in some sense. Thus,
with some discrete approximation gh to g, will not converge, certainly not when the procedure is applied to gd instead of g. However, if we choose h = h(d) suitably, then, subject to certain assumptions,
will satisfy equation 14.48. Often the discretization is done by projecting on finite dimensional subspaces. In the medical imaging community such methods have become known under the names "method of sieves," "generalized series model" (see paper by Liang et al. in section 14.8), "finite series expansion methods" (see Censor's work in the suggested reading list), and others.
In each of these regularization techniques, the choice of the regularization parameter—be it a a , s, k, or h—is crucial. For most problems in medical imaging, choices based on simulations, trial and error, and experience suffice. However, there exist mathematically justified methods for choosing these parameters. Some are completely deterministic, such as the discrepancy principle and the L-curve. Others, such as generalized cross-validation (GCV), are stochastic in nature.
14.4.5 Maximum Entropy In the maximum entropy (ME) method one determines a non-negative solution of the (underdetermined) system 14.46 by maximizing a measure for the entropy of f, e.g., the Shannon entropy
or a discrete version of it. Many algorithms have become known for maximizing entropy, most notably the MART algorithm in CT (see the paper by Herman and Lent listed in section 14.8).
For the nonlinear problem

OCR for page 199

Page 221
with a nonlinear operator A, the well-known Newton method reads
Here, A' denotes the derivative of A. If equation 14.46 is ill posed, equation 14.56 has to be regularized, e.g., by the Tikhonov-Phillips method. Other iterative methods for solving equation 14.55 are the conjugate gradient algorithm and the Levenberg-Marquardt algorithm.
14.5 Sampling All the reconstruction algorithms in medical imaging have to be implemented in a discrete setting. Questions of how densely and in which way the data have to be sampled in order to get a certain image quality are of paramount interest. We distinguish between sampling in real space and sampling in Fourier space.
14.5.1 Sampling in Real Space Sampling in real space occurs, for example, in transmission CT. In principle it can be dealt with using classical sampling theory (such as is associated with the names Shannon, Petersen-Middleton, and Beurling). The two-dimensional case is fairly well understood. Assume that the reconstruction region is the circle with radius r, and that we want to reconstruct reliably functions that are essentially w-band-limited, i.e., whose Fourier transform is negligible in some sense outside the circle of radius w. According to the Shannon sampling theorem, this restriction on the Fourier transform corresponds to a spatial resolution of 2p/w. One of the classical results of transmission CT is that one needs pieces of data to reconstruct such a function reliably. An appropriate sampling geometry is the interlaced parallel geometry suggested as early as 1978 by Cormack (whose 1978 paper is listed in section 14.8). Fan beam sampling requires only slightly more pieces of data (see Natterer's 1993 paper on fan beam sampling listed in the suggested reading). Algorithms are available that actually achieve the resolution 2p/w for these data sets (see, e.g., Faridani in section 14.8).
Not much is known about sampling in three dimensions, not even in fairly well developed fields such as helical scanning in transmission CT or in three-dimensional PET. It seems quite possible that the techniques used in two dimensions lead to substantial data reduction in three dimensions.

OCR for page 199

Page 222
14.5.2 Sampling in Fourier Space Sampling in Fourier space is quite different from sampling in real space and is much less well understood. The problem of recovering a function f from scattered values of its Fourier transform arises in many fields of medical imaging, e.g., in two- and three-dimensional transmission tomography and in MRI; compare to section 4.3. If the xj are the points of a regular cartesian grid the problem can be easily solved by the fast Fourier transform (FFT). However, in two-dimensional transmission tomography the xj are the points of a grid in polar coordinates. This occurs also in some instances of MRI. In other instances, the xj are sitting on a spiral, or they are distributed in a non-symmetric way with respect to the origin. Many methods have been developed for dealing with this situation: the gridding method, a method based on the principle of the stationary phase, methods that estimate the phase of f based on small samples, localized polynomial approximation, autoregressive moving averages (ARMAs), projection on convex subsets (POCS), and maximum entropy (compare section 14.4 and see the paper by Liang et al. in section 14.8). Sometimes these methods work quite well, but clearly better methods are needed and are likely to exist.
14.6 Priors and Side Information The regularization techniques mentioned in section 14.4 permit one to incorporate non-specific side information, such as smoothness and size. In many cases, quite specific side information is at hand, and much effort has been made to make use of such information. An early example is the theory of optimal recovery (Michelli and Rivlin in suggested reading list). An optimal recovery scheme S for the problem 14.46 with side information f Î M is defined as the minimizer of the worst-case error
among all operators S :Y ® X. While this approach looks natural to a mathematician, it obviously has not been of much use. Another natural approach is to minimize
in M. This constraint optimization approach has been used with some success. But it is easy to construct examples in which ignoring the (correct!)

OCR for page 199

Page 223
side information f Î M yields better results. The conclusion from these examples is that the obvious things have already been tried.
We now describe some more recent approaches, beginning with POCS. Assume that the set M is the intersection of convex subsets M1, . . . , Mm, and put M0 = {f : úïAf - gd ÷ï £ d }. Then it is natural to define a sequence f k of approximations to f by successively projecting on M0, M1, . . . , Mm. More specifically, let Pj , j = 0, . . . , m be the orthogonal projections of X onto Mj . Then,
The well-known algebraic reconstruction technique (ART) or Kaczmarz method, which was used in the first commercially available CT-scanner, is a special case of 14.57 with Mj the j-th hyperplane in the system 14.46.
Typical examples of the sets Mj are M1 = {f : a £ f £ b}, M2= {f: f = 0 on some subset}, and M3= {f : úïDb f úï £ r} with a derivative of order b.
POCS permits the use of any side information that can be expressed by convex sets. Even though it has been used successfully for various problems, the mathematical theory is not well developed, and algorithms for POCS tend to be slow. Special methods have been developed for exploiting nonnegativity (e.g., by Janssen; see section 14.8).
Rather than using deterministic side information that can be expressed by some subset M, one can also use stochastic side information. An example is the well-known Bayesian method. Here we think of f and g as families of random variables that are jointly normally distributed with mean values f 0 and g0 and covariance operators F and G, respectively; i.e.,
with E the mathematical expectation, and correspondingly for G. Equation 14.46 is replaced by the linear stochastic model
where n is a family of random variables independent of f that is normally distributed with mean value 0 and covariance å. The Bayesian estimate for f given g is

OCR for page 199

Page 224
The similarity to equation 14.51 is obvious. In fact Bayesian methods very often yield results very similar to those of the Tikhonov-Phillips methods, but provide more flexibility in the modeling of the noise and the properties of f.
A plethora of numerical methods for computing f B have been developed. Most of them are iterative. In the case of circular symmetry, direct methods can be used, too.
In medical imaging, the random variables are usually not normally distributed. Then, the Bayesian estimate is no longer a linear function of the data, and one has to resort to other techniques. One of them is expectation maximization (EM) already mentioned on page 206. In any case, Bayesian methods require a prior F, which models the interpixel correlation of the image f. Various choices for F can be found in the paper by Geman and McClure in the suggested reading list.
14.7 Research Opportunities · Investigation of the trade-offs of stability versus resolution for the inverse problem for the Helmholtz equation, as applied to ultrasound and microwave imaging.
· Development of the mathematical theory and algorithms for inverse problems for the transport equation. There are applications not only to light tomography but also to scatter correction in CT and emission CT. Also needed are investigations of diffusion and diffusion-wave approximations.
· Development of numerical methods for parabolic inverse problems for application to light tomography (the diffusion approximation). Methods using adjoint fields look promising.
· Investigation of Gelfand-Levitan theory for multi-dimensional hyperbolic inverse problems. The one-dimensional Gelfand-Levitan theory is the backbone for one-dimensional inverse scattering; extension to three dimensions would solve the mathematical and numerical problems in ultrasound and microwave imaging.
· Development of a general-purpose algorithm for bilinear inverse problems. The inverse problems of medical imaging frequently have a bilinear structure, irrespective of the type of underlying equation (elliptic,

OCR for page 199

Page 225
parabolic, hyperbolic, or transport). A general-purpose algorithm for discretized problems is conceivable and would save a lot of work.
· Development of methods of scatter correction through transport models for transmission and emission CT. Preliminary results are already available. This challenge should be easier than the inverse problem in light tomography, at least for cases in which scatter is not too large.
· Creation of reconstruction algorithms for three-dimensional CT and efficient algorithms for cone-beam and helical scanning.
· Classification of three-dimensional scanning geometries according to the stability of the inversion problem.
· Development of reconstruction algorithms, possibly of Fourier type, for three-dimensional PET. Recent progress has been made by use of the stationary phase principle.
· Development of faster methods for computing maximum likelihood estimates with priors, more efficient iterative methods, and methods that exploit symmetries of the scanning geometries through efficient numerical algorithms such as the FFT.
· Investigation of preconditioning for nonlinear iterations such as expectation maximization (EM).
· Construction of good priors for SPECT and PET computations.
· Creation of mathematical attenuation corrections for emission CT, i.e., determination of the attenuation map from the emission data (without transmission measurements). Encouraging mathematical results for two source distributions are available, and simulations with templates have been performed.
· Creation of specialized regularization methods, particularly for regularization in time, which would have applications to magnetic and electrical source imaging.
· Reconstruction of functions under global shape information using templates. A typical application is the reconstruction of attenuation maps for emission CT and magnetic and electrical source imaging.

OCR for page 199

Page 226
· Constraint reconstruction, i.e., reconstruction of a function from transmission or emission CT or MRI data with side conditions such as nonnegativity. There are applications to limited-angle CT and to MRI scans with insufficient sampling.
· Reconstruction of a function from irregularly spaced samples of its Fourier transform, which has applications to CT and MRI.
· Removal of artifacts caused by opaque objects, such as hip joints in CT. Many papers are available, but there is no satisfactory solution.
· Reconstruction of a function in R3 from integrals over almost-planar surfaces. There are applications to MRI data collected with imperfect magnets.
· Reconstruction of vector fields. Are there ways to determine more than the curl? If not, what can be concluded from the curl? There are applications to Doppler imaging.
14.8 Suggested Reading 1. Aben, H.K. Integrated Photoelasticity, Valgus, Talin, 1975 (in Russian).
2. Anikonov, D.S., Uniqueness of simultaneous determination of two coefficients of the transport equation, Sov. Math. Dokl. 30 (1984), 149-151.
3. Anikonov, D.S., Prokhorov, I.V., and Kovtanyuk, E.E., Investigation of scattering and absorbing media by the methods of x-ray tomography, J. Invest. Ill-Posed Probl. 1 (1993), 259-281.
4. Arridge, S.R., Van der Zee, P., Cope, M., and Delpy, D.T., Reconstruction methods for infrared absorption imaging, Proc. SPIE 1431 (1991), 204-215.
5. Bondarenko, A., and Antyufeev, V., X-ray Tomography in Scattering Media, Technical report, Institute of Mathematics, Novosibirsk, Russia (1990).
6. Bronnikov, A.V., Degradation transform in tomography, Pattern Recognition Letters 15 (1994), 527-592.

OCR for page 199

Page 227
7. Burridge, R., The Gelfand-Levitan, the Marchenko, and the GopinathSondhi integral equations of inverse scattering theory, regarded in the context of inverse impulse-response problems, Wave Motion 2 (1980), 305-323.
8. Censor, Y., Finite series-expansion reconstruction methods, Proc. IEEE 71 (1983), 409-419.
9. Conference on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, June 23-25, 1993, Snowbird, Utah, Phys. Med. Biol. 39 (1994).
10. Cormack, A.M., Sampling the Radon transform with beams of finite width, Phys. Med. Biol. 23 (1978), 1141-1148.
11. Faridani, A., Reconstructing from efficiently sampled data in parallelbeam computed tomography, in Inverse Problems and Imaging, G.F. Roach, ed., Pitman Res. Notes Math. 245, John Wiley & Sons, New York, 68-102, 1991.
12. Geman, S., and McClure, D., Statistical methods for tomographic image reconstruction, Bull. Int. Stat. Inst. LII (1987), 5-21.
13. Gratton, E., et al., A novel approach to laser tomography, Bioimaging 1 (1993), 40-46.
14. Gutman, S., and Klibanov, M.V., Regularized quasi-Newton method for inverse scattering problems, Math. Comput. Modeling 18 (1993), 5-31.
15. Hamalainen, M.S., and Sarvas, J., Realistic conductivity model of the human head for interpretation of neuromagnetic data, IEEE Trans. Biomed. Eng. BME-36 (1989), 165-171.
16. Herman, G., Image Reconstruction from Projections: The Fundamentals of Computerized Tomography, Academic Press, New York, 1980.
17. Herman, G.T., and Lent, A., Iterative reconstruction algorithms, Comput. Biol. Med. 6 (1976), 273-294.
18. Hertle, A., The identification problem for the constant attenuation Radon transform, Math. Z. 197 (1988), 9-13.

OCR for page 199

Page 228
19. Hinshaw, W.S., and Lent, A.H., An introduction to NMR imaging: From the Bloch equation to the imaging equation, Proc. IEEE 71 (1983), 338-350.
20. Isakov, V., Uniqueness and stability in multi-dimensional inverse problems, Inverse Problems 9 (1993), 579-617.
21. Janssen, A.J.E.M., Frequency-domain bounds for non-negative bandlimited functions, Philips J. Res. 45 (1990), 325-366.
22. Kak, A.C., and Slaney, M., Principles of Computerized Tomography Imaging, IEEE Press, Los Alamidas, Calif., 1987.
23. Kleinman, R.E., and van den Berg, P.M., A modified gradient method for two-dimensional problems in tomography, J. Comput. Appl. Math. 42 (1992), 17-35.
24. Lavrentiev, M.M., Romanov, V.G., and Shishatskij, S.P., Ill-Posed Problems of Mathematical Physics and Analysis, Translations of Mathematical Monographs, vol. 64, American Mathematical Society, Providence, R.I., 1986.
25. Liang, Z.-P., Boada, F.E., Constable, R.T., Haacke, E.M., Lauterbur, P.C., and Smith, M.R., Constrained reconstruction methods in MR imaging, Rev. Magn. Reson. Med. 4 (1992), 67-185.
26. Louis, A.K., Medical imaging: State of the art and future development, Inverse Problems 8 (1992), 709-738.
27. Michelli, C.D., and Rivlin, T.J., eds., Optimal Estimation in Approximation Theory, Plenum Press, New York, 1977.
28. Natterer, F., Determination of tissue attenuation in emission tomography of optically dense media, Inverse Problems 9 (1993), 731-736.
29. Natterer, F., Sampling in Fan Beam Tomography, SIAM J. Appl. Math. 53 (1993), 358-380.
30. Natterer, F., The Mathematics of Computerized Tomography, John Wiley & Sons, New York, 1986.
31. Palamodov, V., An inversion method for attenuated x-ray transform in space, to appear in Inverse Problems.

OCR for page 199

Page 229
32. Quinto, E.T., Cheney, M., and Kuchment, P., eds., Tomography, Impedance Imaging, and Integral Geometry, Proceedings of the AMSSIAM Summer Seminar, June 7-18, 1993, Mount Holyoke College, Massachusetts, American Mathematical Society, Providence, R.I., 1993.
33. Romanov, V.G., Conditional stability estimates for the two-dimensional problem of restoring the right-hand side and absorption in the transport equation (in Russian), Sib. Math. J. 35 (1994), 1184-1201.
34. Sanchez, R., and McCormick, N.J., General solutions to inverse transport problems, J. Math. Phys. 22 (1981), 847-855.
35. Sharafutdinov, V.A., Integral Geometry of Tensor Fields, VSP, Utrecht, The Netherlands, 1994.
36. Shepp, L.A., and Vardi, Y., Maximum likelihood reconstruction for emission tomography, IEEE Trans. Med. Imaging 1 (1982), 115-122.
37. Singer, J.R., Griinbaum, F.A., Kohn, P., and Zubelli, J.R., Image reconstruction of the interior of bodies that diffuse radiation, Science 248 (1990), 990-993.
38. Smith, K.T., Solmon, D.C., and Wagner, S.L., Practical and mathematical aspects of the problem of reconstructing objects from radiographs, Bull. Am. Math. Soc. 83 (1977), 1227-1270.
39. Sparr, G., Stråklén, K., Lindstrom, K., and Persson, W., Doppler tomography for vector fields, Inverse Problems 11 (1995), 1051-1061.
40. Tretiak, O.J., and Metz, C., The exponential Radon transform, SIAM J. Appl. Math. 39 (1980), 341-354.

OCR for page 199