- Email: [email protected]
1D inverse problem in diffusion based optical tomography using iteratively regularized Gauss–Newton algorithm
Applied Mathematics and Computation 161 (2005) 149–170 www.elsevier.com/locate/amc
1D inverse problem in diffusion based optical tomography using iteratively regularized Gauss–Newton algorithm q T. Khan
a,*
, A. Smirnova
b
a
b
Department of Mathematical Sciences, O-201 Martin Hall, P.O. Box 340975, Clemson, SC 29634-0975, USA Department of Mathematics and Statistics, Georgia State University, Atlanta, GA 30303, USA
Abstract In this paper, we investigate an one-dimensional inverse problem in diffusion based optical tomography using iteratively regularized Gauss–Newton (IRGN) algorithm for ill-posed nonlinear problems. We devise a stable reconstruction algorithm for the inverse problem using iterative regularization with Armijo–Goldstein–Wolf (AGW) type line search strategy. We demonstrate the efficacy of the IRGN combined with AGW by reconstructing the scattering parameter relevant to the inverse problem in optical tomography. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Inverse problems; Nonlinear ill-posed; Iteratively regularized Gauss–Newton; Biomedical imaging; Reconstruction algorithms
1. Introduction In the last decade research in the area of biomedical optics has flourished. In particular there has been considerable new development in biomedical imaging
q
This work is supported by NSF grant (award number DMS-0207050). Corresponding author. E-mail addresses: [email protected] (T. Khan), [email protected] (A. Smirnova).
*
0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.019
150
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
using optical tomography [2,8,10,18]. Optical tomography is a way to probe highly scattering media using low-energy visible or near infra-red light (NIR) and then to reconstruct images of these media. Light in the near-infrared range (wavelength from 700 to 1200 nm) penetrates tissue and interacts with it. The predominant effects are absorption and scattering [6,11,14]. We assume: given a set of measurements of transmitted light between pairs of points on the surface of an object, there exists a unique distribution of internal scatters and absorbers which would yield that set. The formation of an image for the optical properties of the tissue from a series of boundary measurements is the inverse problem of optical absorption and scattering tomography (OAST) or optical tomography (OT) for short. The widely accepted photon transport model is the radiative transfer equation (RTE). The transport equation is an integro-differential equation for the radiance and has spatially dependent diffusion and absorption parameters as coefficients which are a priori unknown. Hence the problem is to infer from the measurements of some function of the photon density on the boundary, the coefficients of absorption and diffusion in the tissue. A low order diffusion approximation to the RTE has been derived and studied in the last several years. This is an approximation to the RTE by a parabolic differential equation in the time domain and by an elliptic differential equation in the frequency domain [2]. The diffusion approximation to the transport equation has been widely used to calculate photon migration in biological tissues [15]. The existing computational methods for the inverse problem for photon migration in biological tissues are almost exclusively based on the diffusion approximation [7]. It is well known that the diffusion based inverse problem in optical tomography is exponentially ill-posed or unstable [2,23]. In order to understand the effect of iteratively regularized Gauss–Newton (IRGN) algorithm on the stability of the inverse problem, we investigated the one-dimensional inverse problem in diffusion based optical tomography. The outline of the paper is as follows. In Section 2, we discuss the transport model and its diffusion approximation. In Section 3, we describe system approximation using finite element method (FEM). In Section 4, we describe the iterative Gauss–Newton algorithm for the inverse problem. In Section 5, we discuss our results and work in progress.
2. Photon transport model In optical imaging, low-energy visible light is used to illuminate the biological tissue. The illumination of the tissue can be modelled as a photon transport phenomenon. The process is described by the most widely applied equation in optical imaging, the radiative transfer or transport equation (RTE)
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
151
[9,17]. Let X Rn , n P 2 be bounded with C 2 ––smooth boundary oX, mðxÞ denote the outward unit normal to oX at the point x 2 oX, and C is defined as C :¼ ðx; s; tÞ 2 oX S n1 ½0; T ; mðxÞ s > 0 : S n1 ½0; T Þ and f 2 L1 ðX Let a, b 2 L1 ðXÞ with a > b > 0, u 2 H 1 ðX S ½0; T Þ. Then the RTE is given by n1
1 ou ðx; s; tÞ þ s ruðx; s; tÞ þ aðxÞuðx; s; tÞ c ot Z bðxÞ Hðs s0 Þuðx; s0 ; tÞ ds0 ¼ f ðx; s; tÞ
ð1Þ
S n1
together with the initial and boundary conditions uðx; s; 0Þ ¼ 0 in X S n1 ;
ð2Þ
uðx; s; tÞ ¼ 0
ð3Þ
on C ;
where uðx; s; tÞ describes the density function of particles (photons) which travel in X at time t though the point x in the direction s 2 S n1 , unit sphere in Rn . The parameter a is the total cross section or attenuation, and b is the scattering cross section. The difference l :¼ a b has the physical meaning of an absorption cross section. The parameters a, b and l are assumed to be real, nonnegative and bounded functions of the position. The parameters a and b are the sought-for tissue parameters and c is the velocity of light. The function H is the scattering phase function characterizing the intensity of a wave incident in direction s0 scattered in the direction s. It is assumed to be a real, nonnegative function and is normalized to one. The inverse problem is to recover a and b from measurements of some given functionals of the outgoing density uj jCþ at the boundary oC for ms different set of source distribution fj , j ¼ 1; . . . ; ms . An example of a measured quantity in optical tomography is Z 1 mðxÞ suðx; s; tÞ ds; x 2 oX; t P 0; ð4Þ gðx; tÞ ¼ 4p S 2 where g is known for many incoming fluxes fj . Simpler deterministic models can be derived from RTE for a constant refractive index by expanding the density u and source f in spherical harmonics and retaining a limited number of terms [2,3,5,20]. Due to the prevalence of scattering, the flux is essentially isotropic a small distance away from the sources; i.e., it depends only linearly on s. Thus we may describe the process adequately by the first two moments
152
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
Z 1 uðx; s; tÞ ds; 4p S 2 Z 1 suðx; s; tÞ ds u1 ðx; tÞ ¼ 4p S 2
u0 ðx; tÞ ¼
of u. We assume that u depends only linearly on s, i.e., uðx; s; tÞ ¼ c0 u0 ðx; tÞ þ c1 s u1 ðx; tÞ;
ð5Þ
where the constants c0 ¼ 1 and c1 ¼ 3 are easily obtained by computing the is the mean scattering cosine moments of order 0 and 1. If we let H Z ¼ 1 H s0 sHðs s0 Þ ds0 4p S 2 which does not depend on s and if we introduce the reduced scattering cross and section b0 ¼ ð1 HÞ D¼
1 3ðl þ b0 Þ
which is the diffusion coefficient. We get diffusion approximation P0 (DA): 1 ou0 ðx; tÞ r DðxÞru0 ðx; tÞ þ la ðxÞu0 ðx; tÞ ¼ f0 ðx; tÞ c ot
ð6Þ
together with initial condition u0 ðx; 0Þ ¼ 0
ð7Þ
in X;
and the boundary condition u0 ðx; tÞ þ 2DðxÞ
ou0 ðx; tÞ ¼ 0; om
x 2 oX:
ð8Þ
Similarly the equation for the measurement reads gðx; tÞ ¼ mðxÞ u1 ðx; tÞ ¼ DðxÞmðxÞ ru0 ðx; tÞ ¼ DðxÞ
ou0 : om
ð9Þ
Frequency domain diffusion approximation can easily be obtained by Fourier transforming the time-domain equation. The frequency domain analogue of Eq. (6) is given by ix r DðxÞru0 ðx; xÞ þ la ðxÞ þ ð10Þ u0 ðx; xÞ ¼ f0 ðx; xÞ; c where we made use of the relation o ix: ot
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
153
The frequency domain DE is elliptic where the time-domain DE is parabolic, an important distinction for numerical solutions. Furthermore, the diffusion approximation to the radiative transfer model can be written in the time independent (dc) case as in [2] r DðxÞru0 ðxÞ þ la ðxÞu0 ðxÞ ¼ f0 ðxÞ:
ð11Þ
The associated boundary condition is u0 ðxÞ þ 2DðxÞ
ou0 ðxÞ ¼ 0; om
x 2 oX:
ð12Þ
If we let X be the domain under consideration with surface oX, the weak forward problem corresponding to Eq. (11) is to find uðxÞ 2 H 1 ðXÞ such that for all vðxÞ 2 H 1 ðXÞ, the following variational equation is satisfied Z Z Z Z 1 vu ds ¼ rv Dru dX þ vla u dX þ vf dX: ð13Þ X X oX 2 X Now, we can define the forward problem as: given sources fj in X and q in Q, a vector of model parameters, for example the coefficient of diffusion D and the T coefficient of absorption la (i.e. q ¼ ðD; la Þ ) that belongs to a parameter set Q, find the data u on oX and the inverse problem as: given data z on oX find q. We can recast the forward problem in an abstract setting as the following parameter dependent equation: r DðxÞruðx; qÞ þ la ðxÞuðx; qÞ ¼ f ðxÞ;
ð14Þ
where x is in X, uðqÞ is in an appropriate abstract space H , and f represents a source or a forcing distribution. In general, measurement of uðqÞ may not be possible, only some observable part CuðqÞ of the actual state uðqÞ may be measured. In this abstract setting, the objective of the inverse or parameter estimation problem is to choose a parameter q in Q, that minimizes an error e Q criterion or cost functional J ðuðqÞ; CuðqÞ; qÞ over all possible q in Q subject to uðqÞ satisfying the diffusion approximation. A typical observation operator is m m ou 1 uðxi ; q; f Þ Cf uðqÞ ¼ D ðxi ; q; f Þ ¼ ; ð15Þ om 2 i¼1 i¼1 where xi is in oX, m is the number of measurements, and the second equality comes from the boundary condition (12). For example, in one dimension with X ¼ ½0; L , we have m 2 f1; 2g and x1 ¼ 0, x2 ¼ L. A typical cost functional Jk is given as ms X m
2 1X
fj f
2 ð16Þ Jk ðqÞ ¼
Ci uðqÞ zi j þ kkq ^qk ; 2 j¼1 i¼1
154
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 f
where zi j is the measured data at the boundary for a given source fj , k is the regularization parameter, and ^ q is an a priori choice, for example could be just the nominal values for the parameter. Now composing uðqÞ and CuðqÞ we obtain the parameter-to-output mapping: T ½q ¼ Cu. This is the nonlinear mapping of diffusion based optical tomography in abstract setting. 2.1. A simple example In what follows, we illustrate the theoretical concepts with a simple example in one dimension. Let X ¼ ½0; L , the diffusion approximation with constant background is the Strum–Louiville equation: r2 u þ q2 u ¼
f ; D
ð17Þ
where q2 ¼ la =D is constant, with the Rubin boundary condition: uð0Þ 2Druð0Þ ¼ 0;
ð18Þ
uðLÞ þ 2DruðLÞ ¼ 0:
The inverse problem is to estimate the scalar q from the data z measured at x ¼ 0 or x ¼ L. For example, the GreenÕs function (solution for a delta distribution source f ðxÞ ¼ dðx gÞ where xs ¼ g is the location of the source) for the problem (17) and (18) can be solved analytically and the solution for x < g is: uðx; g; qÞ ¼
ðeqg ceqg Þðeqx beqx Þ ; 2qDðb cÞ
ð19Þ
where c¼
e2Lq ; b
b¼
1 2Dq : 1 þ 2Dq
If we measure Dru m at x ¼ 0, then the inverse problem is to estimate q from the parameter to output map which is given by T ½q ¼ Cuðx; g; qÞ ¼ D
ou ðeqg ceqg Þð1 þ bÞ ð0; g; qÞ ¼ ; om 2ðb cÞ
ð20Þ
which is a nonlinear function of the parameter q as discussed. In Fig. 1, we plot Jk ðqÞ for the 1D diffusion approximation with Rubin boundary conditions for a homogeneous background medium with la ¼ 0:012 mm1 and D ¼ 0:33 mm. This is the simulation offfi the 1D diffusion approximation on the interval pffiffiffiffiffiffiffiffiffiffi ð0; 43:0Þ with q ¼ la =D. We computed cost functional Jk ðqÞ ¼ jTq z1 j2 þ kkqk2 for q0 ¼ 0:1907 mm1 (corresponding to z1 ¼ Tq0 , x1 ¼ 0 la ¼ 0:012 mm1 , and D ¼ 0:33 mm) over a range of q starting from 0:14 to 0:4. The solid
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
7
x 10
155
–6
cost functional Jλ(q)
6 5 4 3 2 1 0 0.15
0.2
0.25 0.3 parameter value q
0.35
0.4
Fig. 1. Parameter estimation and regularization.
curve represents J without regularization (k ¼ 0) and the broken curve represents Jk with regularization parameter k ¼ 105 . From Fig. 1, it is clear that without regularization ðk ¼ 0Þ the function J is rather insensitive to parameter ffi pffiffiffiffiffiffiffiffiffiffi q ¼ la =D (i.e. the numerical method starting with an overestimate of the true parameter is bound to fail). But with regularization ðk ¼ 105 Þ, Jk is more convex. We note here that the regularization has changed the problem so that we are solving for a minimum qk that is no longer the same as the problem without regularization, mainly q0 . This simple example for a constant background illustrates the complexity of a nonlinear ill-posed inverse problem. In the rest of the paper we investigate further the one-dimensional inverse problem in general with spatially varying parameters D and la .
3. System approximation The general analytic method for solving ODEs or PDEs containing a delta distribution source is the GreenÕs function method. However for complex geometries, the analytic solution is intractable. Therefore one requires numerical solutions. The FEM is more versatile than other methods including the finite difference method because of its ease in complex geometries and modelling boundary effects. The FEM is a variational method used to approximate the solution by a family of finite dimensional basis functions. Then the forward problem is reduced to one of linear algebra and one computes the approximate solution using FEM codes.
156
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
If we let X ¼ ½0; L be the 1D domain under consideration with surface oX ¼ f0; Lg, the weak formulation of the problem corresponding to Eq. (11) is to find uðxÞ 2 H 1 ðXÞ such that for all vðxÞ 2 H 1 ðXÞ, the following variational equation is satisfied Z Z 1 1 rv Dru dX þ la vu dX þ vð0Þuð0Þ þ vðLÞuðLÞ 2 2 X X Z ¼ vf dX: ð21Þ X
The FEM is derived by projecting the weak form of (21) onto a finite dimensional function space. The finite dimensional function space is the set of continuous and twice differentiable, piecewise cubic polynomials on subintervals of the domain X. We approximate u with the functions uNp ðxÞ ¼
N p þ1 X
cfi Bi ðxÞ;
ð22Þ
i¼1
where Bi ðxÞ ¼ B
x x i
h
with partitions x1 < x0 ¼ 0 < x1 < < xNp 1 < xNp ¼ L < xNp þ1 ; and stepsize h ¼ xi xi1
for 0 6 i < Np þ 1;
and BðxÞ is the cubic B-spline. By using the test functions v ¼ Bj on the weak form of the PDE in (21), we obtain the Galerkin approximation for all 1 6 j 6 Np þ 1 Z Z 1 1 rBj DruNp dX þ la Bj uNp dX þ Bj ð0ÞuNp ð0Þ þ Bj ðLÞuNp ðLÞ 2 2 X X Z ¼ Bj f dX: ð23Þ X
From the above we obtain a system of equations for cfi where the coefficients of the system equations are given by ij þ Mij þ Pij ; Aij ¼ D where the diffusion inner product Z ij ¼ D rBi DrBj dX; X
ð24Þ
ð25Þ
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
the absorption inner product Z Bi la Bj dX; Mij ¼
157
ð26Þ
X
and the pointwise inner product is 1 ð27Þ Pij ¼ ðBi ð0ÞBj ð0Þ þ Bi ðLÞBj ðLÞÞ: 2 and M using the Gaussian quadrature. We computed the integrals for D Therefore we obtain the system ACf ¼ F;
ð28Þ
where the forcing vector is Z Bi f dX; Fi ¼ X
and Cf is the ðNp þ 3Þ-vector containing the unknown coefficients cfi : Cf ¼ ðcf1 ; cf0 ; . . . ; cfNp þ1 ÞT :
ð29Þ
This linear system then can be solved using any one of the number of techniques. 3.1. Forward simulation To test our forward solver, we first computed the approximate solutions in two cases: (1) constant values of D and la , using f ¼ ex (see Fig. 2, and compare to solution in Fig. 3 where f is shown in Fig. 4 and absolute error is shown in Fig. 5) 3 computed solution 2.5 2 1.5 1 0.5 0 0
10
20
30
40
50
Fig. 2. Constant values of D and la , f ¼ expðxÞ.
158
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 3
exact solution
2.5 2 1.5 1 0.5 0
0
10
20
30
40
50
Fig. 3. Constant values of D and la , f ¼ expðxÞ.
1 forcing distribution
0.8 0.6 0.4 0.2 0
0
10
20
30
40
50
Fig. 4. Constant values of D and la , f ¼ expðxÞ.
1
x 10
–3
0.5
absolute error
0
–0.5
1
0
10
20
30
40
50
Fig. 5. Constant values of D and la , f ¼ expðxÞ.
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
159
Computed solution for approx delta distrib source 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
50 jxxs j2 22
A ffi Fig. 6. Constant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
.
jxx j2 A s f ¼ pffiffiffiffiffiffiffiffiffi e 22 ; 2p2
ð30Þ
with A ¼ 1 and ¼ 0:1 where xs ¼ 42:5 is the source location (see Fig. 8). Here this particular f is an approximation to the dðx xs Þ source for which the GreenÕs function solution is available (see Fig. 6 for approximate solution and compare it to analytic solution in Fig. 7. The absolute error is given in Fig. 9) Fig. 13 illustrates numerical simulation with nonconstant D and la (see Figs. 10 and 11), and f as in (30) with A ¼ 1, ¼ 0:1, and xs ¼ 42:5 (Fig. 12).
4. Inversion using iterative regularized Gauss–Newton algorithm In the last section, we outlined the numerical method for the forward problem. In this section, we will describe the computational method for the
Exact solution for delta distributed source 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
50 jxxs j2 22
A ffi Fig. 7. Constant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
.
160
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 Approximated delta distributed source 3.5 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
Fig. 8. Constant values of D and la , f ¼
50 jxxs j2 22
.
jxxs j2 22
.
A ffi pffiffiffiffiffiffi e 2p2
Absolute error 0.05
0
0.05
0.1
0.15
0
10
20
30
40
50
A ffi Fig. 9. Constant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
The coefficient of absorption mua –3
12
x 10
11 10 9 8 7 6
0
10
20
30
40
50 jxxs j2 22
A ffi Fig. 10. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
.
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 The coefficient of diffusion D 0.6
0.5
0.4
0.3 0
10
20
30
40
50 jxxs j2 22
.
jxxs j2 22
.
jxxs j2 22
.
A ffi Fig. 11. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
Approximated delta distributed source f 3.5 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
50
A ffi Fig. 12. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
Computed solution for variable D and mu a 2.5 2 1.5 1 0.5 0 0
10
20
30
40
50
A ffi Fig. 13. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
161
162
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
inverse problem. We approximate the infinite dimensional parameters D and la with the functions DðxÞ ¼
N X
^ k ðxÞ; dk B
ð31Þ
k¼1
la ðxÞ ¼
N X
^ k ðxÞ; ek B
ð32Þ
k¼1
^ k are the basis functions may be on a different set of grid points than where B the grid points for the basis functions Bk . The model parameter vector q is approximated by T
qN ¼ ðd1 ; . . . ; dN ; e1 ; . . . ; eN Þ ;
ð33Þ
where qN is to be identified from the minimization of the cost functional Jk ms X m
2 1X
fj Np N f
2 Jk qN ¼ ð34Þ
Ci u q zi j þ kkqN ^qN k ; 2 j¼1 i¼1 where T ½qN ¼ CuNp ðqN Þ and k is the regularization parameter. 4.1. Frechet derivative Any numerical method for solving a nonlinear inverse problem requires computation of the Frechet derivative and it is adjoint via linearization. T ½q
associates with q ¼ ðD; la ÞT say the detector response of the form T ½q ¼ CuðqÞ ¼ D
ou ; om
ð35Þ
where u is the solution to Eq. (11) with the forcing distribution f . The inverse problem is to determine D; la from a set of f
f
Ci j uðqÞ ¼ zi j
ð36Þ
detector outputs for a set of sources fj , j ¼ 1; . . . ; ms . Therefore we derive the linearization of T ½q . Let us assume that approximations D0 and l0 for D and la is available such that D ¼ D0 þ h1 , la ¼ l0 þ h2 with h1 and h2 small, and for simplicity, D ¼ D0 near oX. Let v0 be the solution of r D0 ðxÞrv0 ðxÞ þ l0 ðxÞv0 ðxÞ ¼ f0 ðxÞ
ð37Þ
with boundary condition v0 ðxÞ þ 2D0 ðxÞ
ov0 ðxÞ ¼ 0; om
x 2 oX:
ð38Þ
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
163
Let u ¼ v0 þ v where u solves Eqs. (11) and (12) and subtracting Eqs. (37) and (38) from Eqs. (11) and (12), we obtain ignoring second-order terms in h1 , h2 , and v r D0 ðxÞrvðxÞ þ l0 ðxÞvðxÞ ¼ r h1 ðxÞrv0 ðxÞ h2 ðxÞv0 ðxÞ
ð39Þ
with boundary condition vðxÞ þ 2DðxÞ
ov ðxÞ ¼ 0; om
x 2 oX:
The linearization T 0 ½q is given by ov h 0 T ½q 1 ¼ D ; h2 om
ð40Þ
ð41Þ
where v is the solution to Eq. (39). The adjoint of this linearized equation (39) can be derived using the GreenÕs formula Z ð r Drv þ la vÞ w vðr Drw þ la w Þ dX ð42Þ X ! Z ow ov v w D ¼ ds ð43Þ om om oX and the adjoint of T 0 ½q is h rv0 rw T 0 ½q 1 r ¼ ; h2 v0 w
ð44Þ
where w is a solution of the adjoint equation r DðxÞrw ðxÞ þ la ðxÞw ðxÞ ¼ 0
ð45Þ
with the boundary condition w ðxÞ þ 2DðxÞ
ow ðxÞ ¼ r ; om
x 2 oX:
ð46Þ
4.2. Discrete sensitivity relations Now we will derive the discrete sensitivity relations which will be used for the Jacobian calculation for our numerical scheme below. If we plug in the expression for DðxÞ and lðxÞ Eqs. (31) and (32) into the equations for the system matrices D and M Eqs. (25) and (26) respectively, we get ij ¼ D
N X k¼1
k;ij ; dk D
ð47Þ
164
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
where Z
k;ij ¼ D
^ k rBi rBj dX; B
ð48Þ
ek Mk;ij ;
ð49Þ
^ k Bi Bj dX: B
ð50Þ
X
and similarly Mij ¼
N X k¼1
where Mk;ij ¼
Z X
Now if we differentiate Eq. (28) with respect to the parameters qk and if we assume F is independent of qk we get A
oCf oA þ Cf ¼ 0; oqk oqk
ð51Þ
and if we solve for the partial derivatives with respect to the parameters qk oCf þ M þ P Þ1 Dk Cf ; ¼ ðD oqk
1 6 k 6 N;
ð52Þ
oCf þ M þ P Þ1 Mk Cf ; ¼ ðD oqk
N þ 1 6 k 6 2N ;
ð53Þ
þ M þ P, D k and Mk are the ðNp þ 3Þ ðNp þ 3Þ matrices with where A ¼ D k;ij and Mk;ij respectively. elements D 4.3. Numerical scheme Now we describe our numerical scheme for the inverse problem. The inverse problem in Eq. (34) can be casted as a finite dimensional nonlinear minimization problem mainly ms X mtot m
2 1 1X 1X
fj Np N f
2 min J qN ¼ Ci u q zi j ¼ kR qN k ¼ ri2 qN
N 2N 2 j¼1 i¼1 2 2 i¼1 q 2R T
with zero or nonzero residual, where, RðqN Þ ¼ ðr1 ðqN Þ; r2 ðqN Þ; . . . ; rmtot ðqN ÞÞ f f with mtot ¼ m ms , rn ðqN Þ ¼ Ci j uNp ðqN Þ zi j with n ¼ i þ ðj 1Þm, and the N ðqN Þ ¼ K ðqN ÞK ðqN Þ, where K ðqN Þ norm kRk is called the residual at q . Let H N is the Jacobian matrix of Rðq Þ, orn ðqN Þ ð54Þ K qN ¼ ; n ¼ 1; . . . ; mtot ; k ¼ 1; . . . ; N : oqk
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
The gradient is given by rJ qN ¼ K qN R qN ;
165
ð55Þ
and the Hessian is given by qN þ A qN ; ð56Þ r2 J qN ¼ H P m tot where AðqN Þ ¼ i¼1 ri r2 ri ðqN Þ. Thus, if the Jacobian matrix K ðqN Þ is avail of the Hessian r2 J ðqN Þ without requiring any able, we know the first part H second-order information. For example, if we use an observation operator as in Eq. (15) we get Np orn ðqN Þ o ou 1 ouNp N ðxi ; q ; fj Þ ¼ ¼ D ðxi ; qN ; fj Þ: ð57Þ oqk oqk 2 oqk om If we now plug the expression for uNp from Eq. (22) into Eq. (57) we get Np þ1 f orn ðqN Þ 1 X ocl j T oCfj ; ¼ Bl ðxi Þ ¼ B oqk 2 l¼1 oqk oqk
ð58Þ
where ¼ 1 B1 ðxi Þ; B0 ðxi Þ; . . . ; BN ðxi Þ; BN þ1 ðxi Þ T ; B p p 2
ð59Þ
and Cfj is given by Eq. (29). Therefore using the sensitivity equations (52) and (53) for the expression oCfj =oqk in Eq. (58) we can compute the Jacobian matrix KqN for our inverse calculation. One can also use the adjoint formulation to derive this Jacobian matrix which can be beneficial in two and three dimensions. However in one dimension speed of the forward sensitivity relations is adequate to compute the Jacobian of the inverse problem. Given the Jacobian of the inverse problem, we can define the IRGN method as follows:
qNkþ1 ¼ qNk ck ½K qNk K qNk þ kk I 1 K qNk R qNk þ kk ðqNk ^qN Þ ; where kk is a regularizing sequence that satisfies the conditions: kk > 0;
16
kk 6 r; kkþ1
lim kk ¼ 0:
k!1
The reader may consult [1,4,16] for the convergence analysis of IRGN under different types of conditions on the initial guess q0 . The best convergence rate was achieved when we used the iteratively regularized analog of a line search procedure suggested in [22] for the well-posed case. Namely, a backtracking strategy with ck ¼ 1; 1=2; 1=4, until the following two requirements were simultaneously fulfilled: Jkk ðqkþ1 Þ 6 Jkk ðqk Þ þ dchrJkk ðqk Þ; sk i
ð60Þ
166
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
which is the Armijo–Goldstein strategy and krJkk ðqkþ1 Þk 6 qkrJkk ðqk Þk
ð61Þ
which is the Wolfe type strategy. Here sk is the solution to ðK ðqk ÞKðqk Þ þ kk IÞsk ¼ ðK ðqk ÞRðqk Þ þ kk qk Þ. Our code uses the parameters d ¼ 0:0001 and q ¼ 0:99 as in [13]. 4.4. Simulation results The optimization of the functional Jk was performed using the IRGN iterative method described above. After 20 iterations using equally spaced grid points with Np ¼ 40 and N ¼ 20 with k0 ¼ 1:0E 6 with the choice of cubic Bsplines satisfying Rubin boundary condition for the solution and Neumann boundary condition for the parameter, we reconstructed the D profile for a set of 10 sources. We modelled the sources as in Eq. (30) with xs ¼ 0:1, 0.45, 0.8, 1.15, 1.5, 41.0, 41.3750, 41.75, 41.1250, 42.5. The reconstructed profile is depicted in Fig. 14. The line corresponds to true D, the circles correspond to the optimal profile for the space of N -dimensional cubic B-spline basis functions for D, and the stars correspond to the constructed D using IRGN. Since the problem is exponentially ill-posed the proper choice of the initial value of kk was extremely important for the numerical simulations. In Figs. 15– 18 the dependence of the accuracy of the reconstructed solutions on the choice of k0 is illustrated. For a posteriori stopping of IRGN iterations the generalized discrepancy principle (see [4]) was implemented, i.e. the iterations were stopped at the first index ^k, for which the residual kRðq^Nk Þk is less than or equal to sr, s > 1: 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5
exaact D ****** computed D oooo optimal D initial D
0.4 0.3 0.2
0
10
20
30
Fig. 14. lk ¼ 106 /k.
40
50
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
167
1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 15. lk ¼ 103 /k.
1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 16. lk ¼ 104 /k.
1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 17. lk ¼ 105 /k.
kRðq^Nk Þk2 6 sr < kRðqNk Þk2 ;
0 6 k < ^k;
where r is the error of the simulated data in L2 -norm.
ð62Þ
168
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 18. lk ¼ 106 /k.
Table 1 Effect of regularization parameter k0
IRGN iteration
Inner iterations
Relative error
Residual
1.0E)1 1.0E)2 1.0E)3 1.0E)4 1.0E)5 1.0E)6 1.0E)7
10 10 10 10 4 6 1
56 56 56 36 11 28 1
0.12677 0.11652 0.09710 0.08179 0.07483 0.07024 div
7.5060E)4 4.2413E)4 2.3323E)4 1.10069E)4 3.5921E)5 4.5446E)5 div
In Table 1 one can see the summary of measured accuracy and work done for values of k0 in the range from 1.0E)4 to 1.0E)7. For every k0 IRGN approximations were calculated until either (62) was satisfied with sr ¼ 1:0E)4 or the number of iterations k became greater than 10. One can notice that the required accuracy of the residual was attained for two values of k0 :1.0E)5 and 1.0E)6. For k0 ¼ 1:0E)4 the residual is very close to 1.0E)4 but still greater. For other regularizing sequences kk , for example kk ¼ k0 =k n , 0:5 6 n, and kk ¼ k0 = expðkÞ the optimal values of the parameters can also be found, however when kk ¼ k0 =k the accuracy of reconstruction is higher.
5. Conclusion In this paper, we investigated an one-dimensional inverse problem in diffusion based optical tomography. We devised a stable reconstruction algorithm for the inverse problem using IRGN with Armijo–Goldstein–Wolf (AGW) type line search. The algorithm worked very well for the ill-posed inverse problem, for example in comparison to the MATLAB Levenberg–Marquardt trust region nonlinear least square routine LSQNONLIN [12,19,21] with
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 1.2 1.1
1.2 1.1
1
1
0.9
0.9
0.8
0.8
0.7 0.6
0.7 0.6
0.5
0.5
exact D
****** computed D oooo optimal D
0.4 0.3
(a)
169
initial D
0.3
0.2 0
10
20
30
40
exact D
****** computed D oooo optimal D
0.4
50
(b)
0.2
initial D
0
10
20
30
40
50
Fig. 19. D reconstructed using (a) IRGN with AGW algorithm, (b) Levenberg–Marquardt trust region algorithm.
Np ¼ 40 and N ¼ 20, the reconstruction D profile is depicted in Fig. 19(b). The line corresponds to true D, circles correspond to the optimal D, and stars correspond to the constructed D. The comparison of Fig. 19(a) and (b) demonstrates that the Levenberg–Marquardt trust region did not converge to the optimal solution where IRGN/AGW did converge accurately. Therefore, we have demonstrated the efficacy of our algorithm with accurate reconstruction of the scattering parameter relevant to optical tomography. We are currently investigating theoretical convergence and implementation of the proposed IRGN–AGW algorithm in 2D and 3D cases with clinical data for early detection of breast cancer at the Biomedical Optics Laboratory at Clemson University.
References [1] A.B. Bakushinsky, Iterative methods for nonlinear operator equations without regularity, new approach, Dokl. Russian Acad. Sci. 330 (1993) 282–284. [2] S.R. Arridge, Optical tomography in medical imaging: topical review, Inverse Problems 15 (1999) R41–R93. [3] S.R. Arridge, J.C. Hebden, Optical imaging in medicine: 2. Modelling and reconstruction, Phys. Med. Biol. 42 (1997) 841–853. [4] B. Blaschke, A. Neubauer, O. Scherzer, On convergence rates for the iteratively regularized Gauss–Newton method, IMA J. Num. Anal. 17 (1997) 421–436. [5] H. Bremmer, Random volume scattering, Radiat. Sci. J. Res. 680 (1964) 967–981. [6] B. Chance, R.R. Alfano (Eds.), Photon migration and imaging in random media, and tissues, vol. 1888, SPIE, 1993. [7] B. Chance, R.R. Alfano (Eds.), Optical Tomography, Photon migration, and spectroscopy of tissue and model media: theory, human studies, and instrumentations, Part 1 and 2, vol. 2389, SPIE, 1995. [8] B. Chance, C.E. Cooper, D.T. Delpy, E.O.R. Reynolds, Near-infrared spectroscopy and imaging of living systems, Phil. Trans. R. Soc. Lond. B 352 (1997) 643–648. [9] R. Chandrasekhar, Radiation Transfer, Oxford, Clarendon, 1950. [10] National Research Council, Mathematics and Physics of Emerging Biomedical Imaging, National Academy Press, Washington DC, 1996.
170
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
[11] D.T. Delphy, M. Cope, Quantification in tissue near-infrared spectroscopy, Phil. Trans. R. Soc. Lond. B 352 (1997) 649–659. [12] J.E. Dennis, R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equation, Prentice-Hall, Englewood Cliffs, NJ, 1983. [13] M. Hanke, J. Nagy, C. Vogel, Quasi-Newton approach to nonnegative image restorations, Linear Algebra Appl. 316 (2000) 223–236. [14] J. Hebden, S.R. Arridge, D.T. Delphy, Optical imaging in medicine: 1. Experimental techniques, Phys. Med. Biol. 42 (1997) 825–840. [15] A. Hielscher, R. Alcouffe, R. Barbour, Comparison of finite-difference transport and diffusion calculations for photon migration in homogenous and heterogenous tissues, Phys. Med. Biol. 42 (1998) 1285–1302. [16] T. Hohage, Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and inverse scattering problem, Inverse Prob. 13 (1997) 1279– 1299. [17] A. Ishimaru, Single Scattering and Transport, Theory (Wave Propagation and Scattering in Random Media I), Academic, New York, 1978. [18] M. Klibanov, T.R. Lucas, Numerical solution of a parabolic inverse problem in optical tomography using experimental data, SIAM J. Appl. Math. 59 (5) (1999) 1763–1789. [19] K. Levenberg, A method for the solution of certain problems in least squares, Quart. Appl. Math. 2 (2) (1944) 164–168. [20] H.W. Lewis, Multiple scattering in an infinite medium, Phys. Rev. 78 (1950) 526–529. [21] D. Marquardt, An algorithm for least, squares estimation of nonlinear parameters, SIAM J. Appl. Math. 11 (2) (1963) 431–441. [22] S.G. Nash, A. Sofer, Linear and Nonlinear Programming, McGraw-Hill, New York, 1996. [23] F. Natterer, F. W€ ubbeling, Mathematical methods in image reconstruction, SIAM Monographs on Mathematical Modeling and Computation, 2001.
1D inverse problem in diffusion based optical tomography using iteratively regularized Gauss–Newton algorithm q T. Khan
a,*
, A. Smirnova
b
a
b
Department of Mathematical Sciences, O-201 Martin Hall, P.O. Box 340975, Clemson, SC 29634-0975, USA Department of Mathematics and Statistics, Georgia State University, Atlanta, GA 30303, USA
Abstract In this paper, we investigate an one-dimensional inverse problem in diffusion based optical tomography using iteratively regularized Gauss–Newton (IRGN) algorithm for ill-posed nonlinear problems. We devise a stable reconstruction algorithm for the inverse problem using iterative regularization with Armijo–Goldstein–Wolf (AGW) type line search strategy. We demonstrate the efficacy of the IRGN combined with AGW by reconstructing the scattering parameter relevant to the inverse problem in optical tomography. Ó 2003 Elsevier Inc. All rights reserved. Keywords: Inverse problems; Nonlinear ill-posed; Iteratively regularized Gauss–Newton; Biomedical imaging; Reconstruction algorithms
1. Introduction In the last decade research in the area of biomedical optics has flourished. In particular there has been considerable new development in biomedical imaging
q
This work is supported by NSF grant (award number DMS-0207050). Corresponding author. E-mail addresses: [email protected] (T. Khan), [email protected] (A. Smirnova).
*
0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.019
150
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
using optical tomography [2,8,10,18]. Optical tomography is a way to probe highly scattering media using low-energy visible or near infra-red light (NIR) and then to reconstruct images of these media. Light in the near-infrared range (wavelength from 700 to 1200 nm) penetrates tissue and interacts with it. The predominant effects are absorption and scattering [6,11,14]. We assume: given a set of measurements of transmitted light between pairs of points on the surface of an object, there exists a unique distribution of internal scatters and absorbers which would yield that set. The formation of an image for the optical properties of the tissue from a series of boundary measurements is the inverse problem of optical absorption and scattering tomography (OAST) or optical tomography (OT) for short. The widely accepted photon transport model is the radiative transfer equation (RTE). The transport equation is an integro-differential equation for the radiance and has spatially dependent diffusion and absorption parameters as coefficients which are a priori unknown. Hence the problem is to infer from the measurements of some function of the photon density on the boundary, the coefficients of absorption and diffusion in the tissue. A low order diffusion approximation to the RTE has been derived and studied in the last several years. This is an approximation to the RTE by a parabolic differential equation in the time domain and by an elliptic differential equation in the frequency domain [2]. The diffusion approximation to the transport equation has been widely used to calculate photon migration in biological tissues [15]. The existing computational methods for the inverse problem for photon migration in biological tissues are almost exclusively based on the diffusion approximation [7]. It is well known that the diffusion based inverse problem in optical tomography is exponentially ill-posed or unstable [2,23]. In order to understand the effect of iteratively regularized Gauss–Newton (IRGN) algorithm on the stability of the inverse problem, we investigated the one-dimensional inverse problem in diffusion based optical tomography. The outline of the paper is as follows. In Section 2, we discuss the transport model and its diffusion approximation. In Section 3, we describe system approximation using finite element method (FEM). In Section 4, we describe the iterative Gauss–Newton algorithm for the inverse problem. In Section 5, we discuss our results and work in progress.
2. Photon transport model In optical imaging, low-energy visible light is used to illuminate the biological tissue. The illumination of the tissue can be modelled as a photon transport phenomenon. The process is described by the most widely applied equation in optical imaging, the radiative transfer or transport equation (RTE)
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
151
[9,17]. Let X Rn , n P 2 be bounded with C 2 ––smooth boundary oX, mðxÞ denote the outward unit normal to oX at the point x 2 oX, and C is defined as C :¼ ðx; s; tÞ 2 oX S n1 ½0; T ; mðxÞ s > 0 : S n1 ½0; T Þ and f 2 L1 ðX Let a, b 2 L1 ðXÞ with a > b > 0, u 2 H 1 ðX S ½0; T Þ. Then the RTE is given by n1
1 ou ðx; s; tÞ þ s ruðx; s; tÞ þ aðxÞuðx; s; tÞ c ot Z bðxÞ Hðs s0 Þuðx; s0 ; tÞ ds0 ¼ f ðx; s; tÞ
ð1Þ
S n1
together with the initial and boundary conditions uðx; s; 0Þ ¼ 0 in X S n1 ;
ð2Þ
uðx; s; tÞ ¼ 0
ð3Þ
on C ;
where uðx; s; tÞ describes the density function of particles (photons) which travel in X at time t though the point x in the direction s 2 S n1 , unit sphere in Rn . The parameter a is the total cross section or attenuation, and b is the scattering cross section. The difference l :¼ a b has the physical meaning of an absorption cross section. The parameters a, b and l are assumed to be real, nonnegative and bounded functions of the position. The parameters a and b are the sought-for tissue parameters and c is the velocity of light. The function H is the scattering phase function characterizing the intensity of a wave incident in direction s0 scattered in the direction s. It is assumed to be a real, nonnegative function and is normalized to one. The inverse problem is to recover a and b from measurements of some given functionals of the outgoing density uj jCþ at the boundary oC for ms different set of source distribution fj , j ¼ 1; . . . ; ms . An example of a measured quantity in optical tomography is Z 1 mðxÞ suðx; s; tÞ ds; x 2 oX; t P 0; ð4Þ gðx; tÞ ¼ 4p S 2 where g is known for many incoming fluxes fj . Simpler deterministic models can be derived from RTE for a constant refractive index by expanding the density u and source f in spherical harmonics and retaining a limited number of terms [2,3,5,20]. Due to the prevalence of scattering, the flux is essentially isotropic a small distance away from the sources; i.e., it depends only linearly on s. Thus we may describe the process adequately by the first two moments
152
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
Z 1 uðx; s; tÞ ds; 4p S 2 Z 1 suðx; s; tÞ ds u1 ðx; tÞ ¼ 4p S 2
u0 ðx; tÞ ¼
of u. We assume that u depends only linearly on s, i.e., uðx; s; tÞ ¼ c0 u0 ðx; tÞ þ c1 s u1 ðx; tÞ;
ð5Þ
where the constants c0 ¼ 1 and c1 ¼ 3 are easily obtained by computing the is the mean scattering cosine moments of order 0 and 1. If we let H Z ¼ 1 H s0 sHðs s0 Þ ds0 4p S 2 which does not depend on s and if we introduce the reduced scattering cross and section b0 ¼ ð1 HÞ D¼
1 3ðl þ b0 Þ
which is the diffusion coefficient. We get diffusion approximation P0 (DA): 1 ou0 ðx; tÞ r DðxÞru0 ðx; tÞ þ la ðxÞu0 ðx; tÞ ¼ f0 ðx; tÞ c ot
ð6Þ
together with initial condition u0 ðx; 0Þ ¼ 0
ð7Þ
in X;
and the boundary condition u0 ðx; tÞ þ 2DðxÞ
ou0 ðx; tÞ ¼ 0; om
x 2 oX:
ð8Þ
Similarly the equation for the measurement reads gðx; tÞ ¼ mðxÞ u1 ðx; tÞ ¼ DðxÞmðxÞ ru0 ðx; tÞ ¼ DðxÞ
ou0 : om
ð9Þ
Frequency domain diffusion approximation can easily be obtained by Fourier transforming the time-domain equation. The frequency domain analogue of Eq. (6) is given by ix r DðxÞru0 ðx; xÞ þ la ðxÞ þ ð10Þ u0 ðx; xÞ ¼ f0 ðx; xÞ; c where we made use of the relation o ix: ot
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
153
The frequency domain DE is elliptic where the time-domain DE is parabolic, an important distinction for numerical solutions. Furthermore, the diffusion approximation to the radiative transfer model can be written in the time independent (dc) case as in [2] r DðxÞru0 ðxÞ þ la ðxÞu0 ðxÞ ¼ f0 ðxÞ:
ð11Þ
The associated boundary condition is u0 ðxÞ þ 2DðxÞ
ou0 ðxÞ ¼ 0; om
x 2 oX:
ð12Þ
If we let X be the domain under consideration with surface oX, the weak forward problem corresponding to Eq. (11) is to find uðxÞ 2 H 1 ðXÞ such that for all vðxÞ 2 H 1 ðXÞ, the following variational equation is satisfied Z Z Z Z 1 vu ds ¼ rv Dru dX þ vla u dX þ vf dX: ð13Þ X X oX 2 X Now, we can define the forward problem as: given sources fj in X and q in Q, a vector of model parameters, for example the coefficient of diffusion D and the T coefficient of absorption la (i.e. q ¼ ðD; la Þ ) that belongs to a parameter set Q, find the data u on oX and the inverse problem as: given data z on oX find q. We can recast the forward problem in an abstract setting as the following parameter dependent equation: r DðxÞruðx; qÞ þ la ðxÞuðx; qÞ ¼ f ðxÞ;
ð14Þ
where x is in X, uðqÞ is in an appropriate abstract space H , and f represents a source or a forcing distribution. In general, measurement of uðqÞ may not be possible, only some observable part CuðqÞ of the actual state uðqÞ may be measured. In this abstract setting, the objective of the inverse or parameter estimation problem is to choose a parameter q in Q, that minimizes an error e Q criterion or cost functional J ðuðqÞ; CuðqÞ; qÞ over all possible q in Q subject to uðqÞ satisfying the diffusion approximation. A typical observation operator is m m ou 1 uðxi ; q; f Þ Cf uðqÞ ¼ D ðxi ; q; f Þ ¼ ; ð15Þ om 2 i¼1 i¼1 where xi is in oX, m is the number of measurements, and the second equality comes from the boundary condition (12). For example, in one dimension with X ¼ ½0; L , we have m 2 f1; 2g and x1 ¼ 0, x2 ¼ L. A typical cost functional Jk is given as ms X m
2 1X
fj f
2 ð16Þ Jk ðqÞ ¼
Ci uðqÞ zi j þ kkq ^qk ; 2 j¼1 i¼1
154
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 f
where zi j is the measured data at the boundary for a given source fj , k is the regularization parameter, and ^ q is an a priori choice, for example could be just the nominal values for the parameter. Now composing uðqÞ and CuðqÞ we obtain the parameter-to-output mapping: T ½q ¼ Cu. This is the nonlinear mapping of diffusion based optical tomography in abstract setting. 2.1. A simple example In what follows, we illustrate the theoretical concepts with a simple example in one dimension. Let X ¼ ½0; L , the diffusion approximation with constant background is the Strum–Louiville equation: r2 u þ q2 u ¼
f ; D
ð17Þ
where q2 ¼ la =D is constant, with the Rubin boundary condition: uð0Þ 2Druð0Þ ¼ 0;
ð18Þ
uðLÞ þ 2DruðLÞ ¼ 0:
The inverse problem is to estimate the scalar q from the data z measured at x ¼ 0 or x ¼ L. For example, the GreenÕs function (solution for a delta distribution source f ðxÞ ¼ dðx gÞ where xs ¼ g is the location of the source) for the problem (17) and (18) can be solved analytically and the solution for x < g is: uðx; g; qÞ ¼
ðeqg ceqg Þðeqx beqx Þ ; 2qDðb cÞ
ð19Þ
where c¼
e2Lq ; b
b¼
1 2Dq : 1 þ 2Dq
If we measure Dru m at x ¼ 0, then the inverse problem is to estimate q from the parameter to output map which is given by T ½q ¼ Cuðx; g; qÞ ¼ D
ou ðeqg ceqg Þð1 þ bÞ ð0; g; qÞ ¼ ; om 2ðb cÞ
ð20Þ
which is a nonlinear function of the parameter q as discussed. In Fig. 1, we plot Jk ðqÞ for the 1D diffusion approximation with Rubin boundary conditions for a homogeneous background medium with la ¼ 0:012 mm1 and D ¼ 0:33 mm. This is the simulation offfi the 1D diffusion approximation on the interval pffiffiffiffiffiffiffiffiffiffi ð0; 43:0Þ with q ¼ la =D. We computed cost functional Jk ðqÞ ¼ jTq z1 j2 þ kkqk2 for q0 ¼ 0:1907 mm1 (corresponding to z1 ¼ Tq0 , x1 ¼ 0 la ¼ 0:012 mm1 , and D ¼ 0:33 mm) over a range of q starting from 0:14 to 0:4. The solid
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
7
x 10
155
–6
cost functional Jλ(q)
6 5 4 3 2 1 0 0.15
0.2
0.25 0.3 parameter value q
0.35
0.4
Fig. 1. Parameter estimation and regularization.
curve represents J without regularization (k ¼ 0) and the broken curve represents Jk with regularization parameter k ¼ 105 . From Fig. 1, it is clear that without regularization ðk ¼ 0Þ the function J is rather insensitive to parameter ffi pffiffiffiffiffiffiffiffiffiffi q ¼ la =D (i.e. the numerical method starting with an overestimate of the true parameter is bound to fail). But with regularization ðk ¼ 105 Þ, Jk is more convex. We note here that the regularization has changed the problem so that we are solving for a minimum qk that is no longer the same as the problem without regularization, mainly q0 . This simple example for a constant background illustrates the complexity of a nonlinear ill-posed inverse problem. In the rest of the paper we investigate further the one-dimensional inverse problem in general with spatially varying parameters D and la .
3. System approximation The general analytic method for solving ODEs or PDEs containing a delta distribution source is the GreenÕs function method. However for complex geometries, the analytic solution is intractable. Therefore one requires numerical solutions. The FEM is more versatile than other methods including the finite difference method because of its ease in complex geometries and modelling boundary effects. The FEM is a variational method used to approximate the solution by a family of finite dimensional basis functions. Then the forward problem is reduced to one of linear algebra and one computes the approximate solution using FEM codes.
156
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
If we let X ¼ ½0; L be the 1D domain under consideration with surface oX ¼ f0; Lg, the weak formulation of the problem corresponding to Eq. (11) is to find uðxÞ 2 H 1 ðXÞ such that for all vðxÞ 2 H 1 ðXÞ, the following variational equation is satisfied Z Z 1 1 rv Dru dX þ la vu dX þ vð0Þuð0Þ þ vðLÞuðLÞ 2 2 X X Z ¼ vf dX: ð21Þ X
The FEM is derived by projecting the weak form of (21) onto a finite dimensional function space. The finite dimensional function space is the set of continuous and twice differentiable, piecewise cubic polynomials on subintervals of the domain X. We approximate u with the functions uNp ðxÞ ¼
N p þ1 X
cfi Bi ðxÞ;
ð22Þ
i¼1
where Bi ðxÞ ¼ B
x x i
h
with partitions x1 < x0 ¼ 0 < x1 < < xNp 1 < xNp ¼ L < xNp þ1 ; and stepsize h ¼ xi xi1
for 0 6 i < Np þ 1;
and BðxÞ is the cubic B-spline. By using the test functions v ¼ Bj on the weak form of the PDE in (21), we obtain the Galerkin approximation for all 1 6 j 6 Np þ 1 Z Z 1 1 rBj DruNp dX þ la Bj uNp dX þ Bj ð0ÞuNp ð0Þ þ Bj ðLÞuNp ðLÞ 2 2 X X Z ¼ Bj f dX: ð23Þ X
From the above we obtain a system of equations for cfi where the coefficients of the system equations are given by ij þ Mij þ Pij ; Aij ¼ D where the diffusion inner product Z ij ¼ D rBi DrBj dX; X
ð24Þ
ð25Þ
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
the absorption inner product Z Bi la Bj dX; Mij ¼
157
ð26Þ
X
and the pointwise inner product is 1 ð27Þ Pij ¼ ðBi ð0ÞBj ð0Þ þ Bi ðLÞBj ðLÞÞ: 2 and M using the Gaussian quadrature. We computed the integrals for D Therefore we obtain the system ACf ¼ F;
ð28Þ
where the forcing vector is Z Bi f dX; Fi ¼ X
and Cf is the ðNp þ 3Þ-vector containing the unknown coefficients cfi : Cf ¼ ðcf1 ; cf0 ; . . . ; cfNp þ1 ÞT :
ð29Þ
This linear system then can be solved using any one of the number of techniques. 3.1. Forward simulation To test our forward solver, we first computed the approximate solutions in two cases: (1) constant values of D and la , using f ¼ ex (see Fig. 2, and compare to solution in Fig. 3 where f is shown in Fig. 4 and absolute error is shown in Fig. 5) 3 computed solution 2.5 2 1.5 1 0.5 0 0
10
20
30
40
50
Fig. 2. Constant values of D and la , f ¼ expðxÞ.
158
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 3
exact solution
2.5 2 1.5 1 0.5 0
0
10
20
30
40
50
Fig. 3. Constant values of D and la , f ¼ expðxÞ.
1 forcing distribution
0.8 0.6 0.4 0.2 0
0
10
20
30
40
50
Fig. 4. Constant values of D and la , f ¼ expðxÞ.
1
x 10
–3
0.5
absolute error
0
–0.5
1
0
10
20
30
40
50
Fig. 5. Constant values of D and la , f ¼ expðxÞ.
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
159
Computed solution for approx delta distrib source 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
50 jxxs j2 22
A ffi Fig. 6. Constant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
.
jxx j2 A s f ¼ pffiffiffiffiffiffiffiffiffi e 22 ; 2p2
ð30Þ
with A ¼ 1 and ¼ 0:1 where xs ¼ 42:5 is the source location (see Fig. 8). Here this particular f is an approximation to the dðx xs Þ source for which the GreenÕs function solution is available (see Fig. 6 for approximate solution and compare it to analytic solution in Fig. 7. The absolute error is given in Fig. 9) Fig. 13 illustrates numerical simulation with nonconstant D and la (see Figs. 10 and 11), and f as in (30) with A ¼ 1, ¼ 0:1, and xs ¼ 42:5 (Fig. 12).
4. Inversion using iterative regularized Gauss–Newton algorithm In the last section, we outlined the numerical method for the forward problem. In this section, we will describe the computational method for the
Exact solution for delta distributed source 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
50 jxxs j2 22
A ffi Fig. 7. Constant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
.
160
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 Approximated delta distributed source 3.5 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
Fig. 8. Constant values of D and la , f ¼
50 jxxs j2 22
.
jxxs j2 22
.
A ffi pffiffiffiffiffiffi e 2p2
Absolute error 0.05
0
0.05
0.1
0.15
0
10
20
30
40
50
A ffi Fig. 9. Constant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
The coefficient of absorption mua –3
12
x 10
11 10 9 8 7 6
0
10
20
30
40
50 jxxs j2 22
A ffi Fig. 10. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
.
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 The coefficient of diffusion D 0.6
0.5
0.4
0.3 0
10
20
30
40
50 jxxs j2 22
.
jxxs j2 22
.
jxxs j2 22
.
A ffi Fig. 11. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
Approximated delta distributed source f 3.5 3 2.5 2 1.5 1 0.5 0
0
10
20
30
40
50
A ffi Fig. 12. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
Computed solution for variable D and mu a 2.5 2 1.5 1 0.5 0 0
10
20
30
40
50
A ffi Fig. 13. Nonconstant values of D and la , f ¼ pffiffiffiffiffiffi e 2p2
161
162
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
inverse problem. We approximate the infinite dimensional parameters D and la with the functions DðxÞ ¼
N X
^ k ðxÞ; dk B
ð31Þ
k¼1
la ðxÞ ¼
N X
^ k ðxÞ; ek B
ð32Þ
k¼1
^ k are the basis functions may be on a different set of grid points than where B the grid points for the basis functions Bk . The model parameter vector q is approximated by T
qN ¼ ðd1 ; . . . ; dN ; e1 ; . . . ; eN Þ ;
ð33Þ
where qN is to be identified from the minimization of the cost functional Jk ms X m
2 1X
fj Np N f
2 Jk qN ¼ ð34Þ
Ci u q zi j þ kkqN ^qN k ; 2 j¼1 i¼1 where T ½qN ¼ CuNp ðqN Þ and k is the regularization parameter. 4.1. Frechet derivative Any numerical method for solving a nonlinear inverse problem requires computation of the Frechet derivative and it is adjoint via linearization. T ½q
associates with q ¼ ðD; la ÞT say the detector response of the form T ½q ¼ CuðqÞ ¼ D
ou ; om
ð35Þ
where u is the solution to Eq. (11) with the forcing distribution f . The inverse problem is to determine D; la from a set of f
f
Ci j uðqÞ ¼ zi j
ð36Þ
detector outputs for a set of sources fj , j ¼ 1; . . . ; ms . Therefore we derive the linearization of T ½q . Let us assume that approximations D0 and l0 for D and la is available such that D ¼ D0 þ h1 , la ¼ l0 þ h2 with h1 and h2 small, and for simplicity, D ¼ D0 near oX. Let v0 be the solution of r D0 ðxÞrv0 ðxÞ þ l0 ðxÞv0 ðxÞ ¼ f0 ðxÞ
ð37Þ
with boundary condition v0 ðxÞ þ 2D0 ðxÞ
ov0 ðxÞ ¼ 0; om
x 2 oX:
ð38Þ
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
163
Let u ¼ v0 þ v where u solves Eqs. (11) and (12) and subtracting Eqs. (37) and (38) from Eqs. (11) and (12), we obtain ignoring second-order terms in h1 , h2 , and v r D0 ðxÞrvðxÞ þ l0 ðxÞvðxÞ ¼ r h1 ðxÞrv0 ðxÞ h2 ðxÞv0 ðxÞ
ð39Þ
with boundary condition vðxÞ þ 2DðxÞ
ov ðxÞ ¼ 0; om
x 2 oX:
The linearization T 0 ½q is given by ov h 0 T ½q 1 ¼ D ; h2 om
ð40Þ
ð41Þ
where v is the solution to Eq. (39). The adjoint of this linearized equation (39) can be derived using the GreenÕs formula Z ð r Drv þ la vÞ w vðr Drw þ la w Þ dX ð42Þ X ! Z ow ov v w D ¼ ds ð43Þ om om oX and the adjoint of T 0 ½q is h rv0 rw T 0 ½q 1 r ¼ ; h2 v0 w
ð44Þ
where w is a solution of the adjoint equation r DðxÞrw ðxÞ þ la ðxÞw ðxÞ ¼ 0
ð45Þ
with the boundary condition w ðxÞ þ 2DðxÞ
ow ðxÞ ¼ r ; om
x 2 oX:
ð46Þ
4.2. Discrete sensitivity relations Now we will derive the discrete sensitivity relations which will be used for the Jacobian calculation for our numerical scheme below. If we plug in the expression for DðxÞ and lðxÞ Eqs. (31) and (32) into the equations for the system matrices D and M Eqs. (25) and (26) respectively, we get ij ¼ D
N X k¼1
k;ij ; dk D
ð47Þ
164
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
where Z
k;ij ¼ D
^ k rBi rBj dX; B
ð48Þ
ek Mk;ij ;
ð49Þ
^ k Bi Bj dX: B
ð50Þ
X
and similarly Mij ¼
N X k¼1
where Mk;ij ¼
Z X
Now if we differentiate Eq. (28) with respect to the parameters qk and if we assume F is independent of qk we get A
oCf oA þ Cf ¼ 0; oqk oqk
ð51Þ
and if we solve for the partial derivatives with respect to the parameters qk oCf þ M þ P Þ1 Dk Cf ; ¼ ðD oqk
1 6 k 6 N;
ð52Þ
oCf þ M þ P Þ1 Mk Cf ; ¼ ðD oqk
N þ 1 6 k 6 2N ;
ð53Þ
þ M þ P, D k and Mk are the ðNp þ 3Þ ðNp þ 3Þ matrices with where A ¼ D k;ij and Mk;ij respectively. elements D 4.3. Numerical scheme Now we describe our numerical scheme for the inverse problem. The inverse problem in Eq. (34) can be casted as a finite dimensional nonlinear minimization problem mainly ms X mtot m
2 1 1X 1X
fj Np N f
2 min J qN ¼ Ci u q zi j ¼ kR qN k ¼ ri2 qN
N 2N 2 j¼1 i¼1 2 2 i¼1 q 2R T
with zero or nonzero residual, where, RðqN Þ ¼ ðr1 ðqN Þ; r2 ðqN Þ; . . . ; rmtot ðqN ÞÞ f f with mtot ¼ m ms , rn ðqN Þ ¼ Ci j uNp ðqN Þ zi j with n ¼ i þ ðj 1Þm, and the N ðqN Þ ¼ K ðqN ÞK ðqN Þ, where K ðqN Þ norm kRk is called the residual at q . Let H N is the Jacobian matrix of Rðq Þ, orn ðqN Þ ð54Þ K qN ¼ ; n ¼ 1; . . . ; mtot ; k ¼ 1; . . . ; N : oqk
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
The gradient is given by rJ qN ¼ K qN R qN ;
165
ð55Þ
and the Hessian is given by qN þ A qN ; ð56Þ r2 J qN ¼ H P m tot where AðqN Þ ¼ i¼1 ri r2 ri ðqN Þ. Thus, if the Jacobian matrix K ðqN Þ is avail of the Hessian r2 J ðqN Þ without requiring any able, we know the first part H second-order information. For example, if we use an observation operator as in Eq. (15) we get Np orn ðqN Þ o ou 1 ouNp N ðxi ; q ; fj Þ ¼ ¼ D ðxi ; qN ; fj Þ: ð57Þ oqk oqk 2 oqk om If we now plug the expression for uNp from Eq. (22) into Eq. (57) we get Np þ1 f orn ðqN Þ 1 X ocl j T oCfj ; ¼ Bl ðxi Þ ¼ B oqk 2 l¼1 oqk oqk
ð58Þ
where ¼ 1 B1 ðxi Þ; B0 ðxi Þ; . . . ; BN ðxi Þ; BN þ1 ðxi Þ T ; B p p 2
ð59Þ
and Cfj is given by Eq. (29). Therefore using the sensitivity equations (52) and (53) for the expression oCfj =oqk in Eq. (58) we can compute the Jacobian matrix KqN for our inverse calculation. One can also use the adjoint formulation to derive this Jacobian matrix which can be beneficial in two and three dimensions. However in one dimension speed of the forward sensitivity relations is adequate to compute the Jacobian of the inverse problem. Given the Jacobian of the inverse problem, we can define the IRGN method as follows:
qNkþ1 ¼ qNk ck ½K qNk K qNk þ kk I 1 K qNk R qNk þ kk ðqNk ^qN Þ ; where kk is a regularizing sequence that satisfies the conditions: kk > 0;
16
kk 6 r; kkþ1
lim kk ¼ 0:
k!1
The reader may consult [1,4,16] for the convergence analysis of IRGN under different types of conditions on the initial guess q0 . The best convergence rate was achieved when we used the iteratively regularized analog of a line search procedure suggested in [22] for the well-posed case. Namely, a backtracking strategy with ck ¼ 1; 1=2; 1=4, until the following two requirements were simultaneously fulfilled: Jkk ðqkþ1 Þ 6 Jkk ðqk Þ þ dchrJkk ðqk Þ; sk i
ð60Þ
166
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
which is the Armijo–Goldstein strategy and krJkk ðqkþ1 Þk 6 qkrJkk ðqk Þk
ð61Þ
which is the Wolfe type strategy. Here sk is the solution to ðK ðqk ÞKðqk Þ þ kk IÞsk ¼ ðK ðqk ÞRðqk Þ þ kk qk Þ. Our code uses the parameters d ¼ 0:0001 and q ¼ 0:99 as in [13]. 4.4. Simulation results The optimization of the functional Jk was performed using the IRGN iterative method described above. After 20 iterations using equally spaced grid points with Np ¼ 40 and N ¼ 20 with k0 ¼ 1:0E 6 with the choice of cubic Bsplines satisfying Rubin boundary condition for the solution and Neumann boundary condition for the parameter, we reconstructed the D profile for a set of 10 sources. We modelled the sources as in Eq. (30) with xs ¼ 0:1, 0.45, 0.8, 1.15, 1.5, 41.0, 41.3750, 41.75, 41.1250, 42.5. The reconstructed profile is depicted in Fig. 14. The line corresponds to true D, the circles correspond to the optimal profile for the space of N -dimensional cubic B-spline basis functions for D, and the stars correspond to the constructed D using IRGN. Since the problem is exponentially ill-posed the proper choice of the initial value of kk was extremely important for the numerical simulations. In Figs. 15– 18 the dependence of the accuracy of the reconstructed solutions on the choice of k0 is illustrated. For a posteriori stopping of IRGN iterations the generalized discrepancy principle (see [4]) was implemented, i.e. the iterations were stopped at the first index ^k, for which the residual kRðq^Nk Þk is less than or equal to sr, s > 1: 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5
exaact D ****** computed D oooo optimal D initial D
0.4 0.3 0.2
0
10
20
30
Fig. 14. lk ¼ 106 /k.
40
50
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
167
1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 15. lk ¼ 103 /k.
1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 16. lk ¼ 104 /k.
1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 17. lk ¼ 105 /k.
kRðq^Nk Þk2 6 sr < kRðqNk Þk2 ;
0 6 k < ^k;
where r is the error of the simulated data in L2 -norm.
ð62Þ
168
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 1.2 1
D
0.8 0.6 exact D computed D initial D
0.4 0.2 0
10
20
30
40
50
Fig. 18. lk ¼ 106 /k.
Table 1 Effect of regularization parameter k0
IRGN iteration
Inner iterations
Relative error
Residual
1.0E)1 1.0E)2 1.0E)3 1.0E)4 1.0E)5 1.0E)6 1.0E)7
10 10 10 10 4 6 1
56 56 56 36 11 28 1
0.12677 0.11652 0.09710 0.08179 0.07483 0.07024 div
7.5060E)4 4.2413E)4 2.3323E)4 1.10069E)4 3.5921E)5 4.5446E)5 div
In Table 1 one can see the summary of measured accuracy and work done for values of k0 in the range from 1.0E)4 to 1.0E)7. For every k0 IRGN approximations were calculated until either (62) was satisfied with sr ¼ 1:0E)4 or the number of iterations k became greater than 10. One can notice that the required accuracy of the residual was attained for two values of k0 :1.0E)5 and 1.0E)6. For k0 ¼ 1:0E)4 the residual is very close to 1.0E)4 but still greater. For other regularizing sequences kk , for example kk ¼ k0 =k n , 0:5 6 n, and kk ¼ k0 = expðkÞ the optimal values of the parameters can also be found, however when kk ¼ k0 =k the accuracy of reconstruction is higher.
5. Conclusion In this paper, we investigated an one-dimensional inverse problem in diffusion based optical tomography. We devised a stable reconstruction algorithm for the inverse problem using IRGN with Armijo–Goldstein–Wolf (AGW) type line search. The algorithm worked very well for the ill-posed inverse problem, for example in comparison to the MATLAB Levenberg–Marquardt trust region nonlinear least square routine LSQNONLIN [12,19,21] with
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170 1.2 1.1
1.2 1.1
1
1
0.9
0.9
0.8
0.8
0.7 0.6
0.7 0.6
0.5
0.5
exact D
****** computed D oooo optimal D
0.4 0.3
(a)
169
initial D
0.3
0.2 0
10
20
30
40
exact D
****** computed D oooo optimal D
0.4
50
(b)
0.2
initial D
0
10
20
30
40
50
Fig. 19. D reconstructed using (a) IRGN with AGW algorithm, (b) Levenberg–Marquardt trust region algorithm.
Np ¼ 40 and N ¼ 20, the reconstruction D profile is depicted in Fig. 19(b). The line corresponds to true D, circles correspond to the optimal D, and stars correspond to the constructed D. The comparison of Fig. 19(a) and (b) demonstrates that the Levenberg–Marquardt trust region did not converge to the optimal solution where IRGN/AGW did converge accurately. Therefore, we have demonstrated the efficacy of our algorithm with accurate reconstruction of the scattering parameter relevant to optical tomography. We are currently investigating theoretical convergence and implementation of the proposed IRGN–AGW algorithm in 2D and 3D cases with clinical data for early detection of breast cancer at the Biomedical Optics Laboratory at Clemson University.
References [1] A.B. Bakushinsky, Iterative methods for nonlinear operator equations without regularity, new approach, Dokl. Russian Acad. Sci. 330 (1993) 282–284. [2] S.R. Arridge, Optical tomography in medical imaging: topical review, Inverse Problems 15 (1999) R41–R93. [3] S.R. Arridge, J.C. Hebden, Optical imaging in medicine: 2. Modelling and reconstruction, Phys. Med. Biol. 42 (1997) 841–853. [4] B. Blaschke, A. Neubauer, O. Scherzer, On convergence rates for the iteratively regularized Gauss–Newton method, IMA J. Num. Anal. 17 (1997) 421–436. [5] H. Bremmer, Random volume scattering, Radiat. Sci. J. Res. 680 (1964) 967–981. [6] B. Chance, R.R. Alfano (Eds.), Photon migration and imaging in random media, and tissues, vol. 1888, SPIE, 1993. [7] B. Chance, R.R. Alfano (Eds.), Optical Tomography, Photon migration, and spectroscopy of tissue and model media: theory, human studies, and instrumentations, Part 1 and 2, vol. 2389, SPIE, 1995. [8] B. Chance, C.E. Cooper, D.T. Delpy, E.O.R. Reynolds, Near-infrared spectroscopy and imaging of living systems, Phil. Trans. R. Soc. Lond. B 352 (1997) 643–648. [9] R. Chandrasekhar, Radiation Transfer, Oxford, Clarendon, 1950. [10] National Research Council, Mathematics and Physics of Emerging Biomedical Imaging, National Academy Press, Washington DC, 1996.
170
T. Khan, A. Smirnova / Appl. Math. Comput. 161 (2005) 149–170
[11] D.T. Delphy, M. Cope, Quantification in tissue near-infrared spectroscopy, Phil. Trans. R. Soc. Lond. B 352 (1997) 649–659. [12] J.E. Dennis, R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equation, Prentice-Hall, Englewood Cliffs, NJ, 1983. [13] M. Hanke, J. Nagy, C. Vogel, Quasi-Newton approach to nonnegative image restorations, Linear Algebra Appl. 316 (2000) 223–236. [14] J. Hebden, S.R. Arridge, D.T. Delphy, Optical imaging in medicine: 1. Experimental techniques, Phys. Med. Biol. 42 (1997) 825–840. [15] A. Hielscher, R. Alcouffe, R. Barbour, Comparison of finite-difference transport and diffusion calculations for photon migration in homogenous and heterogenous tissues, Phys. Med. Biol. 42 (1998) 1285–1302. [16] T. Hohage, Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and inverse scattering problem, Inverse Prob. 13 (1997) 1279– 1299. [17] A. Ishimaru, Single Scattering and Transport, Theory (Wave Propagation and Scattering in Random Media I), Academic, New York, 1978. [18] M. Klibanov, T.R. Lucas, Numerical solution of a parabolic inverse problem in optical tomography using experimental data, SIAM J. Appl. Math. 59 (5) (1999) 1763–1789. [19] K. Levenberg, A method for the solution of certain problems in least squares, Quart. Appl. Math. 2 (2) (1944) 164–168. [20] H.W. Lewis, Multiple scattering in an infinite medium, Phys. Rev. 78 (1950) 526–529. [21] D. Marquardt, An algorithm for least, squares estimation of nonlinear parameters, SIAM J. Appl. Math. 11 (2) (1963) 431–441. [22] S.G. Nash, A. Sofer, Linear and Nonlinear Programming, McGraw-Hill, New York, 1996. [23] F. Natterer, F. W€ ubbeling, Mathematical methods in image reconstruction, SIAM Monographs on Mathematical Modeling and Computation, 2001.