**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 *x*_{0} with the detector *x*_{1}, equation 14.1 gives rise to the integral

where *dx* is the restriction to *L* of the Lebesgue measure in **R**^{n}. The task is to compute *a* in a domain.W Í **R**^{n}.from the values of equation 14.2 where *x*_{0} and *x*_{1} 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 **R**^{n} ; i.e., q Î *S*^{n-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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 201

For *n* = 3, the relevant integral transform is the x-ray transform, *P,* given by

where q Î *S*^{n}^{-1} and *x* Î q^{^}. *P* admits a similar inversion formula as * R* :

with *K* very similar to equation 14.6 and

where *E*q 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 *x*_{0}*, x*_{l} Î ¶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 *x*_{0} on a circle surrounding supp(a) and *x*_{1} Î *S*^{n}^{-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 **R**^{n} (n = 2 or 3), and *x*_{0} Î ¶W. Since no radiation comes in from outside we have

where *v*_{x} * *is the exterior normal on ¶W at *x* Î ¶W. Equation 14.1 is now replaced by

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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, *h***is 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 Í **R**^{n} for each * q* Î *S*^{n}^{-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 - t*q *,* 0 __<__ *t* < ¥. Thus, equation 14.13 leads to the integral equation

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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} *= E**j* _{ 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* *= (g*_{1}*,. . . ,g*_{n}*)* of measurements is considered a realization of the random variable g = (g_{1},. . . *,*g_{n}), where gi is the number of events detected in detector *i.* The model is determined by the (*n,m*)-matrix *A* = *(a*_{ij} *)* whose elements are

*a*_{ij} * ***=** P(event in pixel/voxel *j* is detected in detector *i),*

where *P* denotes probability. Normalizing the * a*_{ij} 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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 Í **R**^{n} 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*^{-ikt}*u(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 ~ u*q in the integral in equation 14.18. With this approximation, equation 14.17 reads

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 *x*_{0} a source outside the object, consider the initial conditions

We want to determine *c* inside the object from knowing

for many sources *x*_{0} and receivers *x*_{1} 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

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 Î *S*^{n-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*?

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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* - *x*_{0}), where x_{0} Î ¶W is a source point. *u* is considered stationary, too. A second possibility is the light flash *f =* *d*(*x - x*_{0})*d**(t).* Finally, one can also use time harmonic illumination, in which case *f =* *d*(*x - x*_{0}) *e*i*w* * t*. This case reduces to the stationary case with a replaced by *a + i**w* * .* 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} *, b*_{ij} *, r*_{ij} *, e*_{ij} 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 * P*_{ij} *,*_{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 *f*_{ij} *, b*_{ij} *, r*_{ij} *, l*_{ij} from the values of * P*_{ij} *,*_{i' j'} for all boundary pixels. Preliminary

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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* about_{}the direction of that field. The frequency of this precession is the Larmor frequency g*H* where g is the gyromagnetic ratio. By making the magnetic field *H* space-dependent in a controlled way, the local magnetization *M*0*(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, *M*_{i} and * M*_{i}^{0} are the *i*-th components of *M* and *M*^{0}, respectively, and *e*_{i} is the *i*-th unit vector *i* = 1,2,3. The significance of *T* _{ 1}*, T* _{2}, and *M*_{0}

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 210

becomes apparent if equation 14.29 is solved with the static field *H* = *H*_{0}*e*_{3} and initial values * M(x,0)* = *M*^{0}*(x).* Setting *w* _{ 0} = g*H*_{0} * *leads to

Thus the magnetization rotates in the (*x*_{l}, *x*_{2})-plane with Larmor frequency *w* _{ 0} and returns to the equilibrium condition (0, 0, *M*_{0}) with speed controlled by *T* _{2} in the (*x*_{l}, *x*_{2})-plane and by *T* _{1}_{}in the *x*_{3}-direction.

In an MRI scanner one generates a field

where *G* and *H*_{1}_{}are under control. In the language of MRI, *H*_{0}*e*_{3} is the static field, *G* the gradient, and *H*_{1}_{}the radio-frequency (RF) field. The input *G, H*_{1} produces in the detecting system of the scanner the output signal

where *B* characterizes the detecting system. Depending on the choice of *H*_{1}*,* various approximations to *S* can be derived, two of which are detailed here.

Choosing *G* constant for t**£ ***t* £ *r+T* and zero otherwise, we get for *T «* *T* _{ 2}

where is the three-dimensional Fourier transform of *M*_{0}*.*

From here we can proceed in two ways. We can use equation 14.32 to determine the three-dimensional Fourier transform of *M*_{0} * *and compute

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 211

*M*_{0} by an inverse three-dimensional Fourier transform. This requires to be known on a cartesian grid, which can be achieved by a proper choice of the gradients or by interpolation. Alternatively, one can invoke the central slice theorem to obtain the three-dimensional Radon transform *RM*_{0} of * M*_{0} by a series of one-dimensional Fourier transforms. *M*_{0} is then recovered by inverting the three-dimensional Radon transform.

**Case 2, shaped pulse ** In this case, *H*_{1} is the shaped pulse

where *f* is a smooth positive function supported in [0, t]. Then, with *x', G'* the first two components of *x, G,* respectively, we have

where

with a function *Q* essentially supported in a small neighborhood of the origin. Equation 14.33 is the two-dimensional analog of equation 14.32, and again we face the choice between Fourier imaging (i.e., computing the two-dimensional Fourier transform from 14.33 and doing an inverse two-dimensional Fourier transform) and projection imaging (i.e., doing a series of one-dimensional Fourier transforms on equation 14.33 and inverting the two-dimensional Radon transform).

Some of the mathematical problems of MRI, e.g., interpolation in Fourier space, are common to other techniques in medical imaging and are discussed in section 14.5. An interesting mathematical problem occurs if the magnets do not produce sufficiently homogeneous fields. This then calls for the reconstruction of a function from integrals over slightly curved manifolds. Even though this is a problem of classical integral geometry, a satisfactory solution has not yet been found.

##### 14.1.7 Vector Tomography

If the domain under consideration contains a moving fluid, then the Doppler shift can be used to measure the velocity * u*(*x*) of motion. Assume the time-harmonic signal * e*^{i}^{w}^{0t} is transmitted along the oriented line *L.* This signal is

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 212

reflected by particles traveling with speed *v* in the direction of *L* according to *e*^{i(w0-kv)t}*,* where *k* = 2w_{0}/*c* and *c* is the speed of the probing signal; thus, *kv* is the Doppler shift. Let * S(L, v)* be the Lebesgue measure of those particles on *L* that move with speed less than *v;* i.e., *u(x)* × *e*_{L} *< v,* where * e*_{L} is the tangent vector on *L.* Then the total response is

Thus *S* can be recovered from *g* by a Fourier transform. The problem is to recover *u* from *S*.

Not much is known about uniqueness of the solution. However, the first moment of *S,*

is similar to the Radon transform from section 14.1.1. Inverting *R* is a problem of vector tomography. Mathematically, this problem belongs to the recently developed field of integral geometry of tensor fields; see, for example, the paper by Sharafutdinov listed in section 14.8. One can show that curl * u* can be computed from *Ru,* and an inversion formula similar to the Radon inversion formula exists. Numerical simulations of Doppler tomography are given in the work by Sparr et al. listed in section 14.8.

##### 14.1.8 Tensor Tomography

As an immediate extension of transmission CT to anisotropic media, consider a matrix-valued attenuation *a(x) = (a*_{ij}*(x)), i, j =* 1*,. . . ,n.* We then solve the vector differential equation

for the vector-valued function *u(t)* = [*u*_{i}*(t)*]_{i=l}*, . . . ,n .* Let * a* be defined in a convex domain W, and let *x*_{0}, *x*_{1} Î ¶W*.* Then, *u(x*_{1}*) = U(x*_{0}*,x*_{1}*)u(x*_{0} *)* with a nonlinear map * U(x*_{0}*,x*_{1}*)* depending on *a.* The problem is to recover *a* in W from the knowledge of *U(x*_{0}*,x*_{1}*)* for *x*_{0}, *x*_{l}_{}Î ¶W. For *n* = 1 we regain equation 14.1. Applications of equation 14.35 for *n* > 1 have become known in photoelasticity (see Aben's book, listed in the suggested reading), but applications to medicine are not totally out of the question.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 213

In a further extension, let a depend on the direction x *=* (*x*_{1} - *x*_{0})/½*x*_{1} - *x*_{0}½. Such situations occur, for instance, in the polarization of harmonic electromagnetic and elastic waves in anisotropic media. In linearized form these problems give rise to the transverse x-ray transform

and to the longitudinal x-ray transform

Equation 14.36 can easily be reduced to the (*n-* 1)-dimensional x-ray transform in the plane *H*_{w} _{, s} = {*y* : *y* × w* * = *s*}*.* We only have to introduce the function *a*_{w} (*y*) *=* *w*^{T} *a*(*y*)*w**.* Then, 14.36 provides for x Î * H* _{ w,}_{s} * *all the line integrals of *a*_{w}, on *H*_{w},_{s}. For equation 14.37, the situation is not so easy. We decompose * a* into its solenoidal and potential parts; i.e.,

It can be shown that *a*_{1} can be recovered uniquely from 14.37, but *a*_{2} is completely undetermined. This is reminiscent of vector tomography as discussed in subsection 14.1.7.

##### 14.1.9 Magnetic Source Imaging

In magnetic source imaging (MSI) one wants to find the electric currents inside the body from measurements of the induced magnetic field outside the patient. These fields may be caused by epileptic fits or by cardiac infarction. The appropriate mathematical model relating the electric field *J* inside the body W to the magnetic field *B* is the Biot-Savart law (compare section 10.2)

Here, *J = J* ^{i} * *+ s*E* with *J* ^{i} the impressed source current and s*E* the ohmic current. Since s*E* depends on *J* ^{i}*,* 14.38 is a nonlinear equation for *J* ^{ i}*.*

In general the relation between *J* ^{ i} and s*E* is quite complicated and requires extensive computations (see the work by Hämäläinen and Sarvas listed in section 14.8). If the body is modeled as a spherically symmetric

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 214

conductor, this relation can be made explicit. In that case, 14.38 reduces to the linear integral equation

for the dipole moment *Q*. Obviously, equation 14.39 is not uniquely solvable for *Q*. However, it makes sense to compute the generalized solution defined by 14.47 below, which is called the lead field solution in the biomagnetism literature.

Besides non-uniqueness, the main problem in MSI is again instability. One tries to overcome this difficulty either by assuming that there are only a few (1-3) magnetic dipoles or by invoking a strong regularization scheme for a continuous dipole distribution.

##### 14.1.10 Electrical Source Imaging

With electrical source imaging (ESI) one wants to find the electrical potential *u* generated by electrical sources inside the body (typically in the heart) on a surface G_{l} (''epicardial") close to the sources from the potential on the surface G_{0} of the body W. Thus we have (see also section 8.2)

The problem is to recover *u* on G_{1} from the values of * u* on a part of G_{0}. The degree of ill-posedness of this inverse problem depends on how close G_{1} is to the sources. In an idealized setting this can be analyzed quantitatively as was done by Lavrentiev et al. (see suggested reading list). Let G_{1}, G_{0} be concentric circles with radius *rl* < *ro,* respectively, and let all the sources lie inside a third concentric circle with radius *r*_{2} < *r*1. Then, the reconstruction error *e*(d) * *on G_{1} for a data error d on G_{0} is

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 215

In the terminology of section 14.4, the problem is only modestly ill posed if *c* is not too small, i.e., if the sources are not too close to G_{l}.

Unfortunately, in practice G_{l} is always very close to the sources; i.e., *c* is small. Thus the inverse problem 14.40 behaves like a severely ill posed problem, and regularization is mandatory. Let *A : L*_{2}(G_{l}) ® *L*_{2}(G_{0}) be the forward operator, which associates with each potential *f* on G_{l} the function *g* = *u* on G_{0} where *u* is the solution of equation 14.40 with *u* = *f* on G_{l}. Then, the Tikhonov-Phillips regularization is obtained by minimizing

A possible choice for *f* _{ 0} is *f* _{ 0} = 0. However, since measurements are done for many time instants, one can also use an approximation to * f* computed from earlier time frames.

#### 14.2 Forward Problems

Many of the equations described above are derived as approximations to complex physical phenomena. A good formulation in terms of the tissue parameters relative to each modality is the basis of any attack on the corresponding inverse problem. Considerable effort should be spent producing increasingly realistic models that allow for effective numerical simulation and validation in terms of real data.

Numerical methods for solving inverse problems are usually iterative in nature-solving the forward problem over and over again. Thus it is important to have fast and reliable algorithms for the forward problems. Such algorithms are available for impedance tomography, since algorithms for Poissontype equations are a well-established field in numerical analysis (multigrid methods, domain decomposition, preconditioning). For ultrasound, one has to solve Helmholtz-type equations at high spatial frequencies, which still is a challenge. Likewise, the development of codes for transport problems is still an area of active research.

#### 14.3 Inverse Problems

Most of the mathematical problems in medical imaging belong to the class of inverse problems. In an abstract framework these problems can be described as follows. Let *L*(*a*) * *be an operator depending on the parameter *a*. Assume

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 216

the equation

to be solvable for each *f.* Equation 14.42 is called the forward problem in this context. Now assume we are given an "observation operator" *B,* and let

be known. The inverse problem calls for computing *a* from the overdetermined system 14.42-14.43. A typical example is the ultrasound problem discussed in section 14.1.3, in which * L*(*a*) is the Helmholtz operator D*+k*^{2}(l*+a*) or the hyperbolic operator ¶^{2}*u/*¶*t*^{2} - *c*^{2}D*u.* In transmission and optical tomography, *L(a)* is a transport operator. For impedance tomography the operator is elliptic, *L(a)* = div(*a*Ñ*u*). Parabolic operators show up in the diffusion approximation to optical tomography.

Very often the inverse problem 14.42-14.43 can be reduced either exactly or approximately to a linear integral equation of the first kind,

Depending on the smoothness of the kernel *K,* 14.44 is more or less ill posed. We come back to this point in section 14.4.

Since the inverse operator *L*^{-1}(*a*) exists, the inverse problem 14.42-14.43 is equivalent to the single nonlinear equation

In principle, every method for solving nonlinear ill-posed problems (such as a properly regularized Newton method or the Levenberg-Marquardt method of section 14.4) can be used for 14.42-14.43. However, it is clear that the structure of the problem at hand has to be exploited. For instance, in all the cases discussed here, *L*(*a*) depends in an affine-linear way on *a*. Also, *L*(*a*) * *decomposes in many cases naturally into lower-dimensional operators *L*_{j}(*a*)*, j =* 1, . . . , *p,* of identical structure. For instance, *j* may refer to the direction of incidence in x-ray/ultrasound tomography or to the position of the electrode in impedance tomography.

The mathematical theory of inverse problems is still in its infancy. Uniqueness theorems have been obtained only recently. However, the most serious problem is instability. Integral equations of the first kind with smooth kernels, such as 14.44, tend to be severely ill posed; see the next section. The

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 217

few exceptions to this rule include inverse problems for hyperbolic equations. In particular in the one-dimensional case, the Gelfand-Levitan procedure not only yields uniqueness and stability, but also provides efficient numerical methods (background can be found in the work by Burridge listed in section 14.8).

#### 14.4 Ill-Posedness and Regularization

All problems in medical imaging are ill posed, even to very different degrees. We describe ill-posedness for linear equations such as 14.44 and discuss some widely used regularization techniques. The results are extended to nonlinear ill-posed problems. Discussion of the use of nonlinear and stochastic side information is postponed to section 14.6.

Let *A* be a map from a Hilbert space *X* into a Hilbert space *Y* . The equation

is called well posed if

(a) the equation is solvable for each *g,*

(b) the solution *f* is uniquely determined, and

(c) the solution *f* depends continuously on *g*.

Otherwise, equation 14.46 is ill posed. The (Moore-Penrose) generalized solution *f* ^{ +} of 14.46 is defined in the following way: among all the least squares solutions, i.e., the minimizers of úï*Af - g* úï, select the one for which úï*f* úï is minimal. If *A* is linear and bounded, *f* ^{+} can be characterized by

*f* ^{+} is uniquely determined. In many cases in medical imaging, *f* ^{+} is a very useful substitute for the exact solution of equation 14.46.

In medical imaging, the failure to satisfy condition (c) is quite common, which leads to instability. In that case, a small data error *d* of *g* can lead to a very large error in the computed solution *f.*

Stability can be restored by side information. By this we mean information on *f* that is not contained in 14.46, such as * f* Î *M* where *M* is a subset of *X*. Often one can find a function e with e(*d*) ® 0 as *d* ® 0 such that

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 218

the following (conditional) stability estimate holds: let * g*^{d} Î *Y* such that úï*g* - *g*^{d}úï £ d*,* and let *Af* ^{d} ^{}*= g*^{d}*.* Let *f* Î *M.* Then,

The behavior of e(d) as d ® 0 reflects the degree of ill-posedness. If e(d) is of the order d, then there is no ill-posedness at all. On the other hand, if e(d) tends to zero very slowly, e.g., e(d) = 1/log(1/d) (logarithmic stability), then the problem is severely ill posed. For such problems, reducing the reconstruction error by a factor of 2 requires doubling the number of correct digits in the data. Such problems are extremely difficult to handle. As a rule, only very poor approximations to the true solution can be computed. The only problems in medical imaging that have been solved successfully so far belong to the class of modestly ill-posed problems, for which e(d) is of the order d^{c} with 0 < *c* < 1 (Holder stability). For this class, *m* digit accuracy in the data leads to *cm* correct digits in the reconstruction. This loss in accuracy seems bearable if *c* is not too small.

A regularization method is any method that computes from * g*^{d} and from a knowledge of *M * an element *f* ^{d}^{}such that 14.48 holds. Some common methods are explored below.

##### 14.4.1 The Tikhonov-Phillips Method

With the Tikhonov-Phillips method of regularization (see also section 8.2.2), *M =* {*f* Î *X* : úï*f - f* _{0}úï £ *r*} for some *f* _{0} Î *X*, *r* > 0. *F*^{}^{d}^{}is obtained by minimizing

in X. If *A* is linear, then *f* ^{d} can be computed by solving the linear system

or, equivalently,

where *A*^{*} is the adjoint of *A.* The choice of the regularization parameter ca is critical. Of course it depends primarily on d : a = a(d).

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 219

##### 14.4.2 The Truncated Singular Value Decomposition

*A* is said to admit a singular value decomposition (SVD) if

where *(f* _{k}*), (g*_{k}*)* are orthonormal systems in *X*, *Y,* respectively, and s_{k} * *are positive numbers, the singular values of *A.* The spectrum of *A*^{*}*A* consists of the numbers and possibly the eigenvalue 0. The decay of the singular values is a measure for the degree of ill-posedness. Exponential decay corresponds to severe ill-posedness.

The SVD provides more information on the character of the ill-posedness than does any other method. It tells us precisely which features of *f* cannot be recovered in a stable way, namely those that resemble functions *f*_{k} that have small singular values s_{k}. A good example is the limited angle problem in transmission CT where the SVD has been computed analytically.

With the help of SVD one can compute an approximation * f*^{d} * *to *f* by

where s = s(d) plays the role of a regularization parameter. In most cases in medical imaging, the SVD is too expensive to provide an efficient practical reconstruction method. But for analyzing a problem it is a very valuable tool.

##### 14.4.3 Iterative Methods

Let *f* _{k+l} *= B* _{k} *f* _{k} *+ C* _{k} *g* be any convergent iterative method for the solution of *Af = g.* When applied to *g*^{d} instead of *g*, the sequence will in general not converge: even if the first few iterates provide a good approximation to *f,* they will deteriorate as the iteration goes on. This semi-convergence phenomenon has been observed in many iterative methods in medical imaging, most prominently in the EM algorithm in emission tomography (checkerboard effect, see page 206).

Semi-convergence forces one to stop the iteration after a certain number *k(**d* * )* of steps. We put

The optimal choice of *k(**d* *)* is the subject of many papers.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

Page 220

##### 14.4.4 Regularization by Discretization

Let *A*_{h} be a discrete approximation to *A* with step size *A.* For *h* ® 0 we assume *A*_{h} *®** A* in some sense. Since *A*^{-1} is not continuous, (if it exists) tends to infinity in some sense. Thus,

with some discrete approximation *g*_{h} to *g*, will not converge, certainly not when the procedure is applied to * g*^{d} * * 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

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 x_{j} 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 x_{j} are the points of a grid in polar coordinates. This occurs also in some instances of MRI. In other instances, the x_{j} 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!)

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 *M*_{1}*, . . . , M*_{m}*,* and put *M*_{0} *=* {*f* : * *úï*Af - g*^{d} ÷ï £ *d* }*.* Then it is natural to define a sequence *f* _{k} of approximations to * f * by successively projecting on *M*_{0}*, M*_{1}*, . . . , M*_{m}*.* More specifically, let *P*_{j} *, j* = 0, . . . , *m* be the orthogonal projections of *X* onto * M*_{j} *.* 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 *M*_{j} the *j*-th hyperplane in the system 14.46.

Typical examples of the sets *M*_{j} are *M*_{1}* =* {*f* : a £ *f* £ *b*}*, M*_{2}_{}= {*f: f* = 0 on some subset}, and *M*_{3}_{}= {*f* : úï*D*^{b } *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 *g*0 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

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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,

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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 **R**^{3} 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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

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.

**Suggested Citation:**"14 A CROSS-CUTTING LOOK AT THE MATHEMATICS OF EMERGING BIOMEDICAL IMAGING." National Research Council. 1996.

*Mathematics and Physics of Emerging Biomedical Imaging*. Washington, DC: The National Academies Press. doi: 10.17226/5066.

There was a problem loading page 230.