OSA's Digital Library

Optics Express

Optics Express

  • Editor: J. H. Eberly
  • Vol. 4, Iss. 10 — May. 10, 1999
  • pp: 353–371
« Show journal navigation

Truncated Newton’s optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation

R. Roy and E. M. Sevick-Muraca  »View Author Affiliations


Optics Express, Vol. 4, Issue 10, pp. 353-371 (1999)
http://dx.doi.org/10.1364/OE.4.000353


View Full Text Article

Acrobat PDF (214 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

The development of non-invasive, biomedical optical imaging from time-dependent measurements of near-infrared (NIR) light propagation in tissues depends upon two crucial advances: (i) the instrumental tools to enable photon “time-of-flight” measurement within rapid and clinically realistic times, and (ii) the computational tools enabling the reconstruction of interior tissue optical property maps from exterior measurements of photon “time-of-flight” or photon migration. In this contribution, the image reconstruction algorithm is formulated as an optimization problem in which an interior map of tissue optical properties of absorption and fluorescence lifetime is reconstructed from synthetically generated exterior measurements of frequency-domain photon migration (FDPM). The inverse solution is accomplished using a truncated Newton’s method with trust region to match synthetic fluorescence FDPM measurements with that predicted by the finite element prediction. The computational overhead and error associated with computing the gradient numerically is minimized upon using modified techniques of reverse automatic differentiation.

© Optical Society of America

1. Introduction

Conventional imaging modalities such as magnetic resonance imaging (MRI) and x-ray computer-aided tomography (x-ray CT) provide high-resolution medical imaging enabled by the direct geometric correlation between incident and detected radiation. Yet the high cost of operating hospital MRI facilities, and the inability to detect important diseases through x-ray CT imaging suggests an opportunity for the development of near-infrared (NIR) biomedical optical imaging. NIR biomedical optical imaging, or optical tomography, depends upon the low absorbance, yet high scattering of non-mutagenic near-infrared light in the “therapeutic wavelength window” (600–1000 nm) enabling it to safely propagate through several centimeters of tissues. While the propagation of low-energy NIR light can occur safely, the high scattering properties of tissues renders the direct geometric correlation between incident and detected irradiation invalid. The forward imaging problem, (i.e., prediction of the propagation of NIR light through tissues when a map of the tissue optical properties is known) can be described in terms of the diffusion approximation to radiative transfer. The inverse imaging problem (i.e., prediction of the interior optical properties from measurements of light propagation made at the exterior tissue-air interface positions) requires solution of a series of equations that are non-linear in the optical properties to be estimated. In this contribution, we describe a novel algorithm to estimate a solution to the nonlinear, inverse imaging problem for fluorescence frequency-domain photon migration (FDPM). In the following sections of this introduction, we briefly (i) present the background of frequency-domain photon migration (FDPM) imaging, (ii) compare prior work and our current approach towards solution of its inverse imaging problem, and (iii) introduce the concept of fluorescent contrast agents which can accelerate the convergence of the inverse imaging solution.

1.1 Frequency-domain photon migration: forward and inverse problems

Frequency-domain photon migration, depends upon launching intensity modulated (30-200 MHz) light at the interface of a highly scattering medium (such as tissue), and detecting the intensity-modulated wave that successfully propagates to the detector located a distance ρ away from the incident source. Depending upon the spatial distribution of interior absorption and scattering optical properties, the detected light is both phase-delayed by θ, and amplitude attenuated by factor M when compared to the incident light (see Figure 1). The solution to the forward imaging problem involves predicting the phase-delay, θ, and amplitude demodulation, M as a function of angular modulation frequency, ω, and position, r, along the tissue surface given a known spatial distribution of absorption and scattering optical properties. The absorption coefficient, μax , representing the inverse mean absorption path and the isotropic scattering coefficient, μ′sx, equivalent to the inverse isotopic scattering length, govern the propagation of light through scattering media, such as tissues. Since near-infrared light is multiply scattered in tissues, the numerical solution to the diffusion approximation of the radiative transfer equation provides solution to the forward biomedical optical imaging problem. For a propagating intensity modulated wave of light, the optical diffusion equation is written as

·[Dx(r)Φx(r,ω)]+[iωc+μaxi(r)+μaxf(r)]Φx(r,ω)=0onΩ
(1)

where Φx represents the complex number describing the scalar flux of photons at position r; ω is the modulation frequency; c is the speed of light within the medium; μaxi is the absorption due to chromophores and μaxf is the absorption due to fluorophores (μax = μaxi + μaxf). Dx is an optical diffusion coefficient equivalent to 1/[3(μax + μ′sx)]. From numerical solution of Eqn (1) with a Robin boundary condition, the complex scalar flux at the surface, Φx = (=Me iθ) and the phase delay, θ, and amplitude attenuation, M, can be determined. It is important to note that the diffusion equation applies when μax ≪ μ′sx and Dx is approximated by 1/[3(μ′sx)].

Figure 1 Schematic of fluorescence photon migration.

The successful implementation of biomedical optical imaging involves the effective solution of the inverse imaging problem, i.e., determining the interior optical property map of absorption and scattering given measurements of θ and M along the tissue surface. To date, approaches have focused upon the linearization of the problem using first order Born or Rytov approximations [1

1. D. A. Boas, M. A. O’Leary, B. Chance, and A.G. Yodh, “Scattering of diffuse photon density waves by spherical heterogeneities within turbid media: analytic solutions and applications,” Proc. Natl. Acad. Sci. , 91, 4887’91 (1994). [CrossRef] [PubMed]

, 2

2. M. A. O’Leary, D. A. Boas, B. Chance, and A.G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusion photon tomogrpahy,” Opt. Lett. 20, 426’428 (1995). [CrossRef]

]. Non-linear optimization employing iterative Born [3

3. R. L. Barbour, H. Graber, Y. Wang, J. Chang, and R. Aronson, “Perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data,” in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, J. Beuthan, E. Gratton, M. Kashke, B. Masters, S. Svanberg, and P. van der Zee, eds. (SPIE Press, Bellingham, WA., 1993), pp 87’120.

, 4

4. Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R.L. Barbour, “Frequency-domain optical imaging of absorption and scattering by a Born iterative method,” J. Opt. Soc. Am. A 14, 325’342 (1997). [CrossRef]

], or distorted Born methods [3

3. R. L. Barbour, H. Graber, Y. Wang, J. Chang, and R. Aronson, “Perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data,” in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, J. Beuthan, E. Gratton, M. Kashke, B. Masters, S. Svanberg, and P. van der Zee, eds. (SPIE Press, Bellingham, WA., 1993), pp 87’120.

,5

5. W. C. Chew and Y. M. Wang, “Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method,” IEEE Trans. On Medical Imaging. 9, 218’225 (1995). [CrossRef]

] require accurate information regarding the normal “background” optical properties which the tissue “heterogeneity” resides. Since the normal tissue background optical properties can be expected to vary greatly, this a priori information is not realistically feasible for biomedical optical imaging.

While the difficulties associated with the nonlinear and large scale of inverse imaging problem can be addressed with our inversion strategy, additional significant challenges exist when attempting to reconstruct absorption and scattering properties that are physiologically feasible. In a study involving actual optical property measurements of normal and diseased tissues, Troy et al. [15

15. T. L. Troy, D. L. Page, and E. M. Sevick-Muraca, “Optical properties of normal and diseased breast tissues: prognosis for optical mammography,” J. Biomed Opt. 1, 342’355 (1996). [CrossRef] [PubMed]

] found that the optical properties of absorption and isotropic scattering of normal and diseased breast tissues, for example, may not provide sufficient contrast for successful image reconstruction. Indeed, the small influence of “heterogeneity” absorption and scattering on FDPM measurements measured at the boundary can often be within or close to the measurement error of 1 degree in phase and 1% in amplitude attenuation.

1.2 Fluorescence frequency-domain photon migration: forward and inverse problems

In order to over come the difficulties associated with insufficient contrast and costly nonlinear optimization, we and other investigators have developed fluorescence FDPM imaging approaches and inversion strategies [8

8. D.Y. Paithankar, A. U. Chen, B.W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media,” Appl. Opt. 36, 2260’2272 (1997). [CrossRef] [PubMed]

, 16

16. E. M. Sevick-Muraca, G. Lopez, T. L. Troy, J. S. Reynolds, and C. L. Hutchinson, “Fluorescence and absorption contrast mechanisms for biomedical optical imaging using frequency-domain techniques,” Photochem. and Photobio. 6655’64 (1997). [CrossRef]

, 17

17. M. A. O’Leary, D. A. Boas, B. Chance, and A.G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158’160 (1996). [CrossRef]

, 18

18. J. Chang, H. L. Graber, and R.L. Barbour, “Luminescence optical tomography of dense scattering media,” J. Opt. Soc. Am. A. 14, 288’299 (1997). [CrossRef]

]. Fluorescence FDPM physically depends upon the administration of a fluorescent contrast agent, or another agent or gene vector, which results in the expression of a fluorescent protein that emits in the near-infrared wavelength regime. When activated by the absorption of an intensity modulated wave, Φ x (r, ω), an intensity modulated fluorescent wave, Φ m (r, ω), is generated at the same modulation frequency, but is amplitude attenuated and phase-delayed relative to the activating excitation wave. The phase-delay, θm, and amplitude attenuation, Mm, of the generated fluorescent wave can be as high as 90 degrees and a factor of 10’s –100’s of the incident excitation wave owing to the fluorescence properties (lifetime and quantum efficiency) of the dye. Specifically, as the biochemical environment of the fluorophore changes, its lifetime, τ (i.e., the lifetime of the activated fluorophore, or the mean time between absorption of excitation and release of fluorescent light) also alters, providing discrimination of diseased tissues possibly on the basis of biochemical environment. The detected intensity modulated fluorescence wave that has been generated within and which has propagated to the tissue boundary possesses a phase-delay, θm, and amplitude attenuation, Mm, that is exquisitely more sensitive to embedded “heterogeneities” than possible with absorption FDPM imaging described above [16

16. E. M. Sevick-Muraca, G. Lopez, T. L. Troy, J. S. Reynolds, and C. L. Hutchinson, “Fluorescence and absorption contrast mechanisms for biomedical optical imaging using frequency-domain techniques,” Photochem. and Photobio. 6655’64 (1997). [CrossRef]

].

Predictions of fluorescence FDPM measurements of phase-delay, θm, and amplitude attenuation, Mm are achieved through the solution of the complex fluorescence fluence, Φ m (r, ω), at the boundary from the fluorescence optical diffusion equation:

·[Dm(r)Φm(r,ω)]+[iωc+μam(r)]Φm(r,ω)=ϕμaxm11iωτΦx(r,ω)onΩ
(2)

where Dm is the optical diffusion equation at the emission wavelength (=1/3⌊μam+μ′sm⌋)≈1/3⌊μ′sm⌋; μam and μ′sm are the absorption and isotropic scattering coefficients at the fluorescent wavelength; and the right hand term describes the generation of fluorescence within the medium. The term ϕ represents the quantum efficiency of the fluorescence process and the absorption owing to fluorophore is represented by the coefficient μaxm. Here we consider first order relaxation decays only. Note that the source term requires coupling with the solution of excitation fluence described by Eqn (1), where by the absorption coefficient at the excitation light, μax is now provided by the naturally occurring, non-fluorescent chromophores, μaxi, and the fluorescent contrast agent, μaxm. The numerical solution for the excitation fluence distribution Eqn (1) with the Robin boundary condition enables prediction of Φ m from Eqn (2) at the medium boundary and determination of θm and Mm. It is important to note that the diffusion equation applies when μam<<μ′sm.

With the added unknowns at the emission wavelength, the disadvantages using Newton-Raphson optimization for fluorescence FDPM imaging may become more severe. Paithankar et al., [8

8. D.Y. Paithankar, A. U. Chen, B.W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media,” Appl. Opt. 36, 2260’2272 (1997). [CrossRef] [PubMed]

] and Troy and Sevick-Muraca [19

19. T. L. Troy and E. M. Sevick-Muraca, “Fluorescence lifetime imaging and spectroscopy in random media,” in Applied Fluorescence in Chemistry, Biology, and Medicine, Rettig, Strehmel, Shrader, and Seifert, eds., Springer Verlag, pp. 3’36 (1999). [CrossRef]

] employed Newton-Raphson, multi-grid finite difference reconstructions of μaxm τ, and ϕμaxm, given correct priors on other optical properties. Recently, Jiang [20

20. H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element based algorithm and simulations,” Appl. Opt. 37, 5337’5343 (1998). [CrossRef]

] performed similar work, employing dual meshing finite element methods instead of finite difference methods. Both Paithankar et al. [8

8. D.Y. Paithankar, A. U. Chen, B.W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media,” Appl. Opt. 36, 2260’2272 (1997). [CrossRef] [PubMed]

] and Jiang [20

20. H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element based algorithm and simulations,” Appl. Opt. 37, 5337’5343 (1998). [CrossRef]

] resorted to removal of “spurious” optical property values or “filtering” using empirically chosen parameters to achieve the “correct” image distorted through inefficient numerical computation of the Jacobian. O’Leary et al., [17

17. M. A. O’Leary, D. A. Boas, B. Chance, and A.G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158’160 (1996). [CrossRef]

] and Chang, et al., [18

18. J. Chang, H. L. Graber, and R.L. Barbour, “Luminescence optical tomography of dense scattering media,” J. Opt. Soc. Am. A. 14, 288’299 (1997). [CrossRef]

] employed Born and iterative Born which again not only required known “background” optical properties, but which also was unable to handle fluorescence in the “background.” Later Chang et al., [21

21. Chang, H. L. Graber, and R.L. Barbour, “Improved reconstruction algorithm for luminescence optical tomography when background luminophore is present,” Appl Opt. 37 , 3547’3552 (1998). [CrossRef]

] demonstrated the ability to reconstruct fluorescence lifetime under conditions of “imperfect” uptake given the spatial distribution of absorption owing to fluorophore; and Lee and Sevick-Muraca [22

22. J. Lee and E. M. Sevick-Muraca, “Lifetime and absorption imaging with fluorescence FDPM,” Time-resolved fluorescence spectroscopy and imaging in tissues, E. M. Sevick-Muraca (ed.)., Proc. Soc. Photo-Opt. Instrum. Eng., 3600: (to be published), (1999).

] employed an iterative Born-type solution of simultaneous absorption and lifetime from fluorescence FDPM synthetic data.

1.3 Fluorescence FDPM imaging: current work

In this contribution, we also focus upon fluorescence FDPM imaging, owing to the significant influence of interior fluorescent heterogeneities on measured signals at the medium boundary and adapt our inversion strategy for the non-linear, large scale optimization required for biomedical optical imaging of absorption and fluorescence lifetime. In the following section, the numerical methodology for solving the FDPM forward imaging problem and the simulator used for generation of synthetic data sets for testing our inversion algorithm are presented. Section 3 describes our formulation of the inverse imaging problem as an optimization problem using truncated Newton’s methods with trust region and the reverse automatic differentiation in order to reduce the computational overhead associated with optimization. The performance of our algorithm using synthetic data sets is presented in Section 4 and the conclusions, prognosis for FDPM imaging and ongoing work are summarized in Section 5.

2.0 Model prediction of FDPM

2.1 Boundary conditions for FDPM forward simulator

Prediction of FDPM measurements on surface dΩ resulting from diffusion of intensity-modulated waves within volume Ω, is accomplished from solution of Eqn (1) and (2) using the Robin boundary condition. We employ the Robin boundary condition of the form given by Ishimaru[23

23. A. Ishimaru, Wave propagation and scattering in random media,(Academic Press, New York, 1978).

, 24

24. M. Schweiger, S. R. Arridge, M. Hiraka, and D. T. Delpy, “The finite-element method for the propagation of light in scattering media- boundary and source conditions,” Med. Phys. 22, 1779’1792 (1995) [CrossRef] [PubMed]

]:

Φxrω2γDx(r)Φxrωn+Sδ(rrs)=0ondΩ
(3)

Where S is the complex flux density associated with a point source of excitation light located at rs on dΩ; n is the unit vector normal to the surface, and γ accounts for refractive index mismatch at the boundary which through Snell’s reflection, results in internal reflection at the surface.

Groenhuis et al. [25

25. R. A. J. Groenhuis, H. A. Ferwerda, and J. J. Ten Bosch, “Scattering and absorption of turbid material determined from reflection measurements,” Appl. Opt. 22, 2456’2462 (1983). [CrossRef] [PubMed]

] have developed an empirical approach to determine the value of the parameter γ, viz.

γ=(1+rd)(1rd)
(4)

where

rd=1.44nrel2+0.72nrel1+0.668+0.063nrel
(5)

nrel is the relative index of refraction of the scattering medium with respect to the surrounding medium.

2.2 Galerkin finite element formulation for solution of Φx

We employ the Robin boundary condition in our numerical solution of Φx,m employing the Galerkin finite element model with linear triangular elements. In our formulation, we assume μax,m<<μsx,m, so that Dx,m ≈1/3μsx,m. The formulation of the Galerkin finite element solution of Φx begins by multiplying Eqn (1) by a weighting function, wj, and integrating over the domain of interest:

Ω[·(DxΦx)+(iωc+μaxi+μaxf)Φx]wjdΩ=0j=1,2,,N
(6)

where wj {wj : j = 1,⋯,N} is a set of linearly independent weighting functions. N is the number of unknows.

Using Green’s theorem, equation (6) becomes

Ω[Dx(Φx)·(wj)+(iωc+μaxi+μaxf)Φxwj]dΩΓDxwjΦxndΓ=0
(7)

where the second term is a line integral along the external boundary and wj is assumed to be continuous over the whole domain. Now introducing equation (3) in the line integral

ΓDxwjΦxndΓ=1γΓ(Φx+S)wjdΓ=1γΓΦxwjdΓ+1γΓSwjdΓ
(8)

Upon combing Eqns (8) and (7), once can write

Ω[Dx(Φx)(wj)+(iωc+μaxi+μaxf)Φxwj]dΩ1γΓΦxwjdΓ=1γΓSwjdΓ
(9)

Suppose that the solution domain Ω is divided into M 1 triangular elements. Now we can obtain the element equation from Eqn (9) as follows:

el=1M[Ωel[Dxel(Φxel)(wj)+(c+μaxi+μaxfel)Φxelwj]+1γΓelΦxelwj]=
elM[1γΓelSwj]
(10)

Assuming that Φxel , Dxel and μelaxf vary linearly within each triangular element el and μaxi is constant, we may then write Φxel , Dxel and μelaxf for each element as:

Φxel=j=13Lj(Φx)j
(11)

where the Li are the natural co-ordinates for the triangle [26

26. O. C. Zienkiewcz and R. L. Taylor,The finite element methods in engineering science, (McGraw-Hill, New York, 1989).

]. According to the Bubnov-Galerkin method, the weighting functions are chosen to be the same as the approximation functions used to represent Φxel , that is,

wj=Ljforj=1,2,,N

The equation (10) can be written as:

el=1M[Ωel[Dxel(ΦxelxLjx+ΦxelyLjy)+(c+μaxi+μaxfel)Lj]Φxel+1γΓelΦxelLj]=
el=1M[1γΓelSLj]
(12)

After integration of Eqn (12) and rearrangement, the resulting element stiffness equations become:

el=1M[K1el+K2el+K3el]Φxel=elMrel
(13)

Assembly of these element stiffness matrices Kiel and local vector r el yields the system equations KΦ̄ x =b [26

26. O. C. Zienkiewcz and R. L. Taylor,The finite element methods in engineering science, (McGraw-Hill, New York, 1989).

].

Since S denotes the complex point source applied to the surface dΩ. at position rs and S is expressed as:

S=Sδ(rrS)
(14)

where rs is the position vector of the point source. By definition, we have

ΓelSδ(rrs)dΓ={SrsΓel0otherwise
(15)

where el is the element of the finite element mesh. From Eqn(12) we have

rel=1γΓelLj(rrs)=1γLj(rs)S

If the position rs of the point source coincides with the position of node j, when Lj(rs) = 1, i.e. we have

(rel)j=1γS

2.3 Galerkin formulation of finite element solution Φm

The forward finite element solution of the emission diffusion equation is formulated similarly to the excitation diffusion equation with the exception that (i) the complex excitation fluence must be first computed, and that (ii) the formulation of fluorescence lifetime poses some additional considerations. As in the formulation of the solution for Φx, parameters of μaxm, τ, and Dm are expected to vary linearly within the triangular elements. In order to facilitate the formulation of the Galerkin finite element method, the RHS term of Eqn (2) is expressed in terms of a binomial expansion assuming ωτ<1 and higher order terms are neglected. Hence, the RHS source term of Eqn (2) becomes Sm (r,ω)=ϕμaxm[1 + iω(r)τ]Φ x (r, ω). This assumption is appropriate, since the phase-delay, θm, of the generated fluorescent intensity modulated wave reaches a maximum of 90 degrees with respect the local excitation intensity modulated wave and becomes insensitive to fluorescence lifetime when ωτ>1. The formulation of the emission finite element solution for Φ m is identical to that presented in section 2.2 for Φx with the exception that the parameters Φ m , Dm, μaxm, and τ are assumed to vary linearly within each element. The finite element formulation of emission diffusion is given below:

el=1M[Ωel[Dmel(ΦmelxLjx+ΦmelyLjy)+(iωc+μam)ΦmelLj]dΩ+1γΓelΦmelLj]=
el=1el[ΩelΦxelμaxmelϕ(1+ωτel)LjdΩ]
(16)

It is important to emphasize that the absorption owing to fluorophore is equivalent in our case to μaxm and to μaxf. We differentiate between the two parameters when reconstructing absorption from synthetic excitation, μaxf, and absorption from emission data μax → m.

2.4 Solution of the system of finite element equations

For both excitation and emission finite element formulations, the local stiffness matrices, Kjel, are formed and then combined into a complex, global stiffness matrix, K, which depends upon the absorption and scattering properties at each wavelength. For solution of excitation fluence, Φ̄x, the vector b reflects the complex source of excitation light from boundary point sources at positions rs, while for solution of the emission fluence, Φ̄m, vector b represents the source generated from within volume Ω.

To solve for Φ̄x,m , the complex set of sparse, linear equations generated by the finite element method,

KΦ̄x,m=b
(17)

are efficiently solved using a subroutines ZGBTRF and ZGBTRS (LAPACK subroutines). ZGBTRF computes a LU factorization of a complex M 2 by N band matrix using partial pivoting with row interchanges. (M 2 = 2*KL + KU + 1, KL and KU are the numbers of subdiagonals and superdiagonals within the band of K, respectively). ZGBTRS solves a system of linear equations with a general band matrix K using the LU factorization computed by ZGBTRF. The main advantage of this method is that once the decomposition is done, the solution vector can be obtained for any number of right hand side vectors. The accuracy of the solution for Φ̄x,m depends crucially on the mesh size, h, with the accuracy improving with h 2.

3.0 Inverse solution for fluorescence FDPM

The variables in the inverse problems are the absorption coefficients, μaxf (for excitation), μaxm (for emission), and the lifetime τ at each nodal point of a finite element model. If there are N nodal points then the number of measurements must be greater than N. The number of measurements made at NB boundary nodes for Ns sources, requires that

(NB1)NSN
(18)

Due to the reciprocality of sources and detectors, i.e., when a detector position is swapped with a source position equivalent measurements are registered, then we require:

(NB1)NS2N
(19)

If l =1,…,Ns denotes the different distribution of source s and j =1,…, NB denotes the boundary nodes at which measurements are made, the error function for optimization of absorption and lifetime imaging may be taken as

Ex,m(μaxf)=12l=1Nsj=1NBjl((((Φx,m)l)j)c(((Φx,m)l)j)me(((Φx,m)l)j)me)((((Φx,m*)l)j)c(((Φx,m*)l)j)me(((Φx,m*)l)j)me)
(20)

The subscript c denotes the values calculated by the forward simulator problem and the subscript me denotes the experimental (or in this case synthetically generated) fluence values. The superscript * denotes the complex conjugate of the complex function Φx,m. The subscript l and j denote fluence values that arise between source l and detector j.

It is seen from the equation (20) that the calculation of error function E involves Ns solutions of the direct problem. The global stiffness matrix K x,m for excitation and emission remains constant with only the vector b changing with each excitation source. Since the global stiffness matrix K x depends on μ̄axf and b is independent of μ̄axf, LU decomposition of K x for solution of Φx need only be performed once in each iteration of the optimization problem.

An efficient way of calculating Ex,m for excitation and emission equations consists of the following steps:

  • Given μaxf calculate the local stiffness matrices K el (μ̄axf)
  • K=elKeland set F = 0
  • Decompose K into the LU decomposition LU = K

    For each source l

  • calculate b(s)
  • Solve, LUΦ̄x,m =b gives (Φ̄x,m)
  • Ex,m=12jdΩ((Φx,m)c(Φx,m)me(Φx,m)me)j((Φx,m*)c(Φx,m*)me(Φx,m*)me)j
  • F = F + Ex,m

As most of the efficient optimization algorithms, the truncated Newton method requires the function value Ex,m and the gradient of the error function, ∇Ex,m , in order to efficiently update parameter estimates. The gradients of the error function in Eqn (20) represent the direction in which the error function is increasing/decreasing most rapidly. Since absorption coefficients μaxf μaxm and τ are the unknown parameters, the gradients are obtained by taking partial derivatives of Eqn (20) with respect to μaxf, μax→m and τ.

Ex=Ex(μaxf)=RejdΩ(((Φx*)c(Φx*)me(Φx)me(Φx*)me),(Φx)c(μaxf))
(21)
Em=Em(τ)=RejdΩ(((Φm*)c(Φm*)me(Φm)me(Φm*)me),(Φm)c(τ))
(22)
Em=Em(μaxm)=RejdΩ(((Φm*)c(Φm*)me(Φm)me(Φm*)me),(Φm)c(μaxm))
(23)

〈.,.〉 denotes the inner product of two complex vectors and Re(.) denotes the real part of a complex. number. The gradient of the real valued error function with respect to real valued μ̄axf μax→m , and τ remains real valued. At minimum these gradients should be zero, so the real part of the gradients are zero. Therefore, we consider only the real part of the gradients. In the following, the optimization strategy using the truncated Newton method with the gradients computed by reverse automatic differentiation is described.

3.1 Truncated Newton’s method as an optimization strategy for FDPM imaging

The truncated Newton method with trust region is used for non-linear optimization problem. The reader interested in background literature is referred to the references [27

27. L. C. W. Dixon and R. C. Price, “Numerical experience with the truncated Newton method for unconstrained optimization,” JOTA , 56, 245’255 (1988). [CrossRef]

, 28

28. R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, England,(1996).

].

Newton’s method is based upon approximating the function Ex,m (μ̄ak +d) (μa denotes μaxf {for excitation} or μaxm {for emission} in this section) at the kth iteration by the quadratic model:

Ex,m(μ¯ak+d)=Ex,m(μ¯ak)+Q(d)
Q(d)=gkTd+12dTGkd
(24)

Gkd=gk
(25)

A sequence of approximate solutions of the above equation is generated by the conjugate gradient method until a vector d i (see Step 3) is obtained for which the following condition on the relative residual is satisfied as suggested by Dembo and Steihaug [29

29. R. S. Dembo and T. Steihaug, “Truncated Newton algorithms for large-scale unconstrained optimization,” Math Programming 26, 190’212 (1983). [CrossRef]

];

rigkmin(1k,gk)

G(x)d=1σ[g(x+σd)g(x)]
(26)

where σ=machine precisiond. This avoids the calculation and storage of the Hessian but requires an additional gradient evaluation for each minor iteration. The pseudo codes for truncated Newton optimization is as follows:

  • Initialization

    An initial guess (μ̄a0 is made of the solution: The radius R 1 = 0.01 of the trust region is defined around (μ̄a0), E 0 = E(μ̄a0), g 0 = g(μ̄a0) are computed.

  • Stopping Criteria

    If ∥g k∥ ≤ ε exit, ε is specified by the user

  • Computing Newton Direction

    Apply conjugate gradient method to find an approximate solution d k of the Newton equation G k d k = -g k such that ∥d k∥ ≤ Rk and dkT g k ≤ -ε1d k∥∥g k∥, typically ε1 = 0.001 [30

    30. R. C. Price, Sparse matrix optimization using automatic differentiation, Ph. D. thesis, University of Hertfordshire, U. K., (1987).

    ]. The conjugate gradient algorithm ensures that at the j-th iteration d j+1 =d j + αj p j remains within the trust region (αj steplength and pj is the conjugate direction)

  • Line Search

    Condition 2:

    E(μ̄a + αd) - E(μ̄a) ≤ ε2αd T g where ε2 is a constant, such that 0.5 > ε2 0 (typically ε2 = 0.1).

    Condition 3

    E(μ̄a + 2αd) > E(μ̄a + ε3αd T g) where ε3 is a constant, such that 0.5 > ε3 0

  • New Approximation

    Compute μ̄ak+1 = μ̄ak + λk d k, E x,mk+1, g k+1

  • Update the trust region

    Once λk has been found the radius of the trust region R k+1 is altered as follows:

    R1=0.01
    Rk+1=2Rkifλk1.0
    Rk+1=13Rkifλk<1.0

A similar formula applies for lifetime τ and absorption coefficient μax→m.

3.2 Reverse automatic differentiation for computation of ∇Ex,m

The reverse automatic differentiation (RAD) method is used to calculate the gradient g = g(μ̅a or τ̅ or μ̅axm). The reader interested in background literature is referred to References [28

28. R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, England,(1996).

, 33

33. B. Christianson, A. J. Davies, L. C. W. Dixon, R. Roy, and P. van der Zee, “Giving reverse differentiation a helping hand,” Opt. Meth. And Software 8, 53’67 (1997). [CrossRef]

, and 34

34. A. J. Davies, B. Christianson, L. C. W. Dixon, R. Roy, and P. van der Zee, ”Reverse differentiation and the inverse diffusion problem,” Adv. In Eng. Software 28, 217’221 (1997). [CrossRef]

]. A brief description of the reverse differentiation and analytic implementation of this method is discussed in the following section.

Reverse automatic differentiation is discussed in terms of the Wengert [35

35. R. E. Wengert, “A simple automatic derivative evaluation program,” Comm. A. C. M. , 7, 463’464 (1964).

] list. This list decomposes complicated function of many variables into a sequence of simpler functions of one or two variables. Functions of two variables are called the “binary function”. Binary functions are either addition, or subtraction or multiplication. Functions of one variable are either reciprocation, raise of power, exponential, logarithm, trigonometric, or hyperbolic. These functions are defined as unary functions.

If f(x) is a function of the n variables x 1,…,xn then a set of new variables x n+1,…x P are introduced, where P is the number of arithmetic operations involved in calculating the function f(x). A Wengert list for the calculation for the calculation of f(x) consists of a list of these unary and binary function processed in a given order;

Givenxi,i=1,,nFori=n+1,,PthenifFiisbinaryxi=Fi(xj,xk),j,k<iandifFiisunaryxi=F(xj),.j<if(x)=xn+P

Reverse differentiation can be explained as follows. Given the Wengert list for the calculation of f a set of adjoint variables = ∂f /∂xi is introduced and the order of the operations reversed. The method is based on the function of a function rule of calculus that implies

fxi=jfxjFjxij>i

The automatic method involves a pass in the reverse direction through the list after the function has been evaluated in a forward sweep. Hence the method is:

Givenxi,i=1,,Psetx̂i=0,i=1,.,n+P1andx̂n+P=1fori=n+P,,n+1
then, ifFiisbinary,x̂j=x̂j+x̂iFixji>jandx̂k=x̂k+x̂iFixki>jelse ifFiisunary,x̂j=x̂j+x̂iFixji>jderivativesgi=x̂ii=1,,n

Griewank [14

14. A. Griewank, “On automatic differentiation,” edited M. Iri and K. Tanaka, Mathematical programming: Recent developments and application, (Kluwer Academic Publishers, 1989) pp 83’108.

] showed that it is always possible to calculate any gradient vector of a function in less than 3 times the computational cost of the function by the reverse automatic differentiation method irrespective of the dimension of the problem. The main disadvantage of the automatic code is that it requires large storage and significant data accessing overheads [28

28. R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, England,(1996).

, 33

33. B. Christianson, A. J. Davies, L. C. W. Dixon, R. Roy, and P. van der Zee, “Giving reverse differentiation a helping hand,” Opt. Meth. And Software 8, 53’67 (1997). [CrossRef]

, and 34

34. A. J. Davies, B. Christianson, L. C. W. Dixon, R. Roy, and P. van der Zee, ”Reverse differentiation and the inverse diffusion problem,” Adv. In Eng. Software 28, 217’221 (1997). [CrossRef]

]. In order to overcome these difficulties, the reverse differentiation method is extended and has been performed analytically. This does not require large storage and at the same time gradients are calculated in less than three times the operation count of the function. The gradients of Ex,m with respect to parameter μaxf are calculated using the reverse differentiation method as follows.

  • Fx,m=0,K̂=Ex,mK=0,b̂=Ex,mb=0
  • Calculate element matrix K el (μ̄a)
  • Assemble global matrix K
  • Decompose K = LU

    For each light source s:

  • Calculate b(s)
  • Solve LU Φ̄ x,m = b(s) to give (Φ̄x,m)c
  • Ex,m=j12((Φx,m)c(Φx,m)me(Φx,m)me)j((Φx,m*)c(Φx,m*)me(Φx,m*)me)j
  • Fx,m = Fx,m + Ex,m
  • Φ¯̂x,m=Ex,mΦx,m={((Φx,m*)c(Φx,m*)me(Φx,m)me(Φx,m*)me)jfor alljdΩ0for alljdΩ
  • Solve LUv¯x,m=Φ¯̂x,m
  • Update = - v̄x,mTΦ̄x,m, = + x,m

    (μ̂axf)p=Ex,m(μaxf)p=Ex,mKK(μaxf)p+Ex,mbb(μaxf)p

  • =eli,j(K̂el)i,j(Ki,jel(μaxf)p)+eljb̂jbj(μaxf)p

The above gradients of error function Ex,m are calculated using the global stiffness matrices, K, and the matrix b assembled from the excitation and emission diffusion equations as described in Section 2.2 and 2.3. (The reverse differentiation method described above uses the fact that if Φ x,m is the solution of the equation KΦ̅ x,m = b then = vT x,m Φ̂ x,m and b = x,m where x,m is the solution of the equation Kv¯x,m=Φ¯̂x,m, see Appendix A). With ∇Ex,m and Ex,m computed, the Newton directions can be computed as described in Section 3.1. Similar formulas hold for τ̂ and μ̂ax→m

4 Conclusions

We have formulated the solution to the absorption and fluorescence FDPM optical imaging problem as a non-linear optimization scheme. We develop the truncated Newton method for minimizing the sum of errors between measured and synthetic data generated using a simplistic finite element model with a minimal number of sources and detectors. Since the truncated Newton method requires only the product of the Hessian with a direction at any point, it is computationally efficient for large optimization problems. We couple the truncated Newton method with a finite element solver which assumes a polynomial series in absorption and lifetime in order to accurately compute Φx,m and employ principles of reverse differentiation in order to efficiently compute the gradient of the error function, ∇Ex,m . Using synthetic data with simulated measurement noise, we demonstrate in Part II the optimization strategy on a series of studies involving a minimum number of sources and detectors.

Appendix A

Calculation of adjoint variables

In this appendix we derive the expressions for the adjoint variables K̂=Ex,mk and b̂=Ex,mb.

  1. Adjoint variable

    K̂=Ex,mK=Ex,mΦ¯x,mΦ¯x,mK=Φ¯̂x,mΦ¯x,mK=Φ¯̂x,mΦx,mμ¯axfμ¯axfK
    (A1)

    (since K is a function of μaxf, we differentiate K with respect of μaxf ). Now we use the global equation KΦ̅x,m = b to find μ¯axfK as follows:

    KΦ̄x,m=b
    Kμ̄axfΦ¯x,m+KΦ¯xfμ¯axf=0

    Since the source term b is constant in excitation equation.

    Kμ¯axf=-KΦ¯x,mμ¯xfΦ¯x,m
    μ¯xfK=K1Φ¯x,mΦ¯x,mμ¯ax,m
    (A2)

    On substituting Eqn (A2) in Eqn (A1) we have

    K̂=Φ¯̂x,mK1Φx,m=v¯TΦ¯x,m
    (A3)

    where

    v¯=Φ¯̂K1
    Kv¯=Φ¯̂x,m
    (A4)

    We can solve the Eqn(A4) and obtain . Finally is obtained using Eqn(A3).

  2. Adjoint variable

    b̂=Ex,mb=Ex,mΦ¯x,mΦ¯x,mb=Φ¯̂x,mΦ¯x,mb=Φ¯̂x,mΦ¯x,mτ¯τ¯b
    (A5)

    (since b is a function of τ̅, we differentiate b with respect of τ̅). Now we use the global equation KΦ̅ x,m=b to find τ¯b as follows:

    KΦ¯x,m=b
    KΦ¯x,mτ¯=bτ¯
    τ¯b=K11Φ¯x,mτ¯
    (A6)

    Substitute Eqn(A6) into Eqn(A5) we have

    Kb̂=Φ¯̂x,m
    (A7)

    Eqn(B7) is the same as Eqn (B4). Hence = v̅

Mamplitude modulation
rposition
cspeed of light
S(r s )source at location rs and modulated frequency ω.
Doptical diffusion coefficient = 1/3 [μa + μ′s]
nnormal
rs source position
rd parameter describing γ
nrel relative refractive functions
Wj independent weighting functions
Lj co-ordinates for triangular elements
Kiel element stiffness matrices
r el local vector in FEM stiffness equations
blocal vector in FEM stiffness equation containing source terms
Zvariance
G(0,1)random number with Gaussian distribution of zero mean and unit variance
Nnumber of nodal points
Ns number of sources
NB number of boundary nodes
ρdistance between source and detector pair
θphase-delay
ωangular modulation frequency
μa absorption coefficient
μs isotropic scattering coefficient
Φcomplex fluence
ϕquantum efficiency
τfluorescence lifetime
Ωvolume
γparameter to account for refractive index mismatch
Γsurface
εstopping criteria for truncated Newton’s method
εerror function for all sources
Ferror function for individual sources
g k gradient vector at iteration k = ∇Ek
dsearch direction
GHessian matrix

References

1.

D. A. Boas, M. A. O’Leary, B. Chance, and A.G. Yodh, “Scattering of diffuse photon density waves by spherical heterogeneities within turbid media: analytic solutions and applications,” Proc. Natl. Acad. Sci. , 91, 4887’91 (1994). [CrossRef] [PubMed]

2.

M. A. O’Leary, D. A. Boas, B. Chance, and A.G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusion photon tomogrpahy,” Opt. Lett. 20, 426’428 (1995). [CrossRef]

3.

R. L. Barbour, H. Graber, Y. Wang, J. Chang, and R. Aronson, “Perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data,” in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, J. Beuthan, E. Gratton, M. Kashke, B. Masters, S. Svanberg, and P. van der Zee, eds. (SPIE Press, Bellingham, WA., 1993), pp 87’120.

4.

Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R.L. Barbour, “Frequency-domain optical imaging of absorption and scattering by a Born iterative method,” J. Opt. Soc. Am. A 14, 325’342 (1997). [CrossRef]

5.

W. C. Chew and Y. M. Wang, “Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method,” IEEE Trans. On Medical Imaging. 9, 218’225 (1995). [CrossRef]

6.

K. D. Paulsen and H. Jiang, “Spatially varying optical property reconstruction using a finite element diffusion equation approximation,” Med. Phys. 22, 691’701 (1995). [CrossRef] [PubMed]

7.

H. Jiang, K. D. Paulsen, U. L. Osterberg, B.W. Pogue, and M. S. Patterson, “Optical image reconstruction using frequency-domain data simulations and experiments,” J. Opt. Soc. Am. A. 13, 253’266 (1996). [CrossRef]

8.

D.Y. Paithankar, A. U. Chen, B.W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media,” Appl. Opt. 36, 2260’2272 (1997). [CrossRef] [PubMed]

9.

M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the finite-element method for the forward and inverse models in optical tomography,” J. Math. Img. Vision 3, 263’283 (1993). [CrossRef]

10.

J. J. McKeown ,“On algorithms for sums of squares problems,” Towards global optimization, edited by L. C. W. Dixon and G. P. Szeeg, (North-Holland Amsterdam, Holland, 1975).

11.

M. T. Vespucci, “An efficient code for the minimization of highly nonlinear and large residual least squares functions,” Optimization 18, 825’855 (1987). [CrossRef]

12.

R. R. Meyer, “Theoretical and computational aspects of nonlinear regression,” Nonlinear Programming, eds. J. B. Rosen, O. L. Mangasarian, and K. Ritter, (Academic Press, New York, 1970).

13.

L. B. Rall, Automatic differentiation: Techniques and application, Lecture notes in computer science (Springer Verlag, 1981) p. 120. [CrossRef]

14.

A. Griewank, “On automatic differentiation,” edited M. Iri and K. Tanaka, Mathematical programming: Recent developments and application, (Kluwer Academic Publishers, 1989) pp 83’108.

15.

T. L. Troy, D. L. Page, and E. M. Sevick-Muraca, “Optical properties of normal and diseased breast tissues: prognosis for optical mammography,” J. Biomed Opt. 1, 342’355 (1996). [CrossRef] [PubMed]

16.

E. M. Sevick-Muraca, G. Lopez, T. L. Troy, J. S. Reynolds, and C. L. Hutchinson, “Fluorescence and absorption contrast mechanisms for biomedical optical imaging using frequency-domain techniques,” Photochem. and Photobio. 6655’64 (1997). [CrossRef]

17.

M. A. O’Leary, D. A. Boas, B. Chance, and A.G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158’160 (1996). [CrossRef]

18.

J. Chang, H. L. Graber, and R.L. Barbour, “Luminescence optical tomography of dense scattering media,” J. Opt. Soc. Am. A. 14, 288’299 (1997). [CrossRef]

19.

T. L. Troy and E. M. Sevick-Muraca, “Fluorescence lifetime imaging and spectroscopy in random media,” in Applied Fluorescence in Chemistry, Biology, and Medicine, Rettig, Strehmel, Shrader, and Seifert, eds., Springer Verlag, pp. 3’36 (1999). [CrossRef]

20.

H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element based algorithm and simulations,” Appl. Opt. 37, 5337’5343 (1998). [CrossRef]

21.

Chang, H. L. Graber, and R.L. Barbour, “Improved reconstruction algorithm for luminescence optical tomography when background luminophore is present,” Appl Opt. 37 , 3547’3552 (1998). [CrossRef]

22.

J. Lee and E. M. Sevick-Muraca, “Lifetime and absorption imaging with fluorescence FDPM,” Time-resolved fluorescence spectroscopy and imaging in tissues, E. M. Sevick-Muraca (ed.)., Proc. Soc. Photo-Opt. Instrum. Eng., 3600: (to be published), (1999).

23.

A. Ishimaru, Wave propagation and scattering in random media,(Academic Press, New York, 1978).

24.

M. Schweiger, S. R. Arridge, M. Hiraka, and D. T. Delpy, “The finite-element method for the propagation of light in scattering media- boundary and source conditions,” Med. Phys. 22, 1779’1792 (1995) [CrossRef] [PubMed]

25.

R. A. J. Groenhuis, H. A. Ferwerda, and J. J. Ten Bosch, “Scattering and absorption of turbid material determined from reflection measurements,” Appl. Opt. 22, 2456’2462 (1983). [CrossRef] [PubMed]

26.

O. C. Zienkiewcz and R. L. Taylor,The finite element methods in engineering science, (McGraw-Hill, New York, 1989).

27.

L. C. W. Dixon and R. C. Price, “Numerical experience with the truncated Newton method for unconstrained optimization,” JOTA , 56, 245’255 (1988). [CrossRef]

28.

R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, England,(1996).

29.

R. S. Dembo and T. Steihaug, “Truncated Newton algorithms for large-scale unconstrained optimization,” Math Programming 26, 190’212 (1983). [CrossRef]

30.

R. C. Price, Sparse matrix optimization using automatic differentiation, Ph. D. thesis, University of Hertfordshire, U. K., (1987).

31.

L. Armijo “Minimization of functions having Lipschitz continuous first partial derivatives,” Pacific J. Mathematics 16, 1’3 (1966).

32.

P. Wolfe, “Convergence condition for ascent method,” SIAM Rev. , 11226’253 (1969). [CrossRef]

33.

B. Christianson, A. J. Davies, L. C. W. Dixon, R. Roy, and P. van der Zee, “Giving reverse differentiation a helping hand,” Opt. Meth. And Software 8, 53’67 (1997). [CrossRef]

34.

A. J. Davies, B. Christianson, L. C. W. Dixon, R. Roy, and P. van der Zee, ”Reverse differentiation and the inverse diffusion problem,” Adv. In Eng. Software 28, 217’221 (1997). [CrossRef]

35.

R. E. Wengert, “A simple automatic derivative evaluation program,” Comm. A. C. M. , 7, 463’464 (1964).

OCIS Codes
(100.3190) Image processing : Inverse problems
(170.0170) Medical optics and biotechnology : Medical optics and biotechnology
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(260.2510) Physical optics : Fluorescence

ToC Category:
Research Papers

History
Original Manuscript: March 11, 1999
Published: May 10, 1999

Citation
R. Roy and E. Sevick-Muraca, "Truncated Newton's optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation," Opt. Express 4, 353-371 (1999)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-4-10-353


Sort:  Journal  |  Reset  

References

  1. D. A. Boas, M. A. OLeary, B. Chance and A.G. Yodh, "Scattering of diffuse photon density waves by spherical heterogeneities within turbid media: analytic solutions and applications," Proc. Natl. Acad. Sci. 91, 4887-91 (1994). [CrossRef] [PubMed]
  2. M. A. OLeary, D. A. Boas, B. Chance and A.G. Yodh, "Experimental images of heterogeneous turbid media by frequency-domain diffusion photon tomogrpahy," Opt. Lett. 20, 426-428 (1995). [CrossRef]
  3. R. L. Barbour, H. Graber, Y. Wang J. Chang and R. Aronson, "Perturbation approach for optical diffusion tomography using continuous-wave and time-resolved data," in Medical Optical Tomography: Functional Imaging and Monitoring, G. Muller, B. Chance, R. Alfano, J. Beuthan, E. Gratton, M. Kashke, B. Masters, S. Svanberg and P. van der Zee, eds. (SPIE Press, Bellingham, WA., 1993), pp 87-120.
  4. Y. Yao, Y. Wang, Y. Pei, W. Zhu and R.L. Barbour, "Frequency-domain optical imaging of absorption and scattering by a Born iterative method," J. Opt. Soc. Am. A 14, 325-342 (1997). [CrossRef]
  5. W. C. Chew and Y. M. Wang, "Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method," IEEE Trans. On Medical Imaging. 9, 218-225 (1995). [CrossRef]
  6. D. Paulsen and H. Jiang, "Spatially varying optical property reconstruction using a finite element diffusion equation approximation," Med. Phys. 22, 691-701 (1995). [CrossRef] [PubMed]
  7. H. Jiang, K. D. Paulsen, U. L. Osterberg, B.W. Pogue and M. S. Patterson, "Optical image reconstruction using frequency-domain data simulations and experiments," J. Opt. Soc. Am. A. 13, 253-266 (1996). [CrossRef]
  8. D.Y. Paithankar, A. U. Chen, B.W. Pogue, M. S. Patterson and E. M. Sevick-Muraca, "Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media," Appl. Opt. 36, 2260-2272 (1997). [CrossRef] [PubMed]
  9. M. Schweiger, S. R. Arridge. D. T. Delpy, " Application of the finite-element method for the forward and inverse models in optical tomography," J. Math. Img. Vision 3, 263-283 (1993). [CrossRef]
  10. J. J. McKeown,"On algorithms for sums of squares problems," Towards global optimization, edited by Dixon, L. C. W. and Szeeg, G. P., (North-Holland Amsterdam, Holland, 1975).
  11. M. T. Vespucci, "An efficient code for the minimization of highly nonlinear and large residual least squares functions," Optimization 18, 825-855 (1987). [CrossRef]
  12. R. R. Meyer, "Theoretical and computational aspects of nonlinear regression," Nonlinear Programming, eds. Rosen, J. B., Mangasarian, O. L. and Ritter, K. (Academic Press, New York, 1970).
  13. L. B. Rall, Automatic differentiation: Techniques and application, Lecture notes in computer science (Springer Verlag, 1981) p. 120. [CrossRef]
  14. A. Griewank, "On automatic differentiation," edited Iri, M. and Tanaka, K., Mathematical programming: Recent developments and application, (Kluwer Academic Publishers, 1989) pp 83-108.
  15. T. L. Troy, D. L. Page and E. M. Sevick-Muraca, "Optical properties of normal and diseased breast tissues: prognosis for optical mammography," J. Biomed Opt. 1, 342-355 (1996). [CrossRef] [PubMed]
  16. E. M. Sevick-Muraca, G. Lopez,, T. L. Troy, J. S. Reynolds and C. L. Hutchinson, "Fluorescence and absorption contrast mechanisms for biomedical optical imaging using frequency-domain techniques," Photochem. Photobiol. 66 55-64 (1997). [CrossRef]
  17. M. A. OLeary, D. A. Boas, B. Chance and A. G. Yodh, "Fluorescence lifetime imaging in turbid media," Opt. Lett. 21, 158-160 (1996). [CrossRef]
  18. J. Chang, H. L. Graber and R.L. Barbour, "Luminescence optical tomography of dense scattering media," J. Opt. Soc. Am. A. 14, 288-299 (1997). [CrossRef]
  19. T. L. Troy and E. M. Sevick-Muraca, "Fluorescence lifetime imaging and spectroscopy in random media," in Applied Fluorescence in Chemistry, Biology, and Medicine, Rettig, Strehmel, Shrader, Seifert, eds., Springer Verlag, pp. 3-36 (1999). [CrossRef]
  20. H. Jiang, "Frequency-domain fluorescent diffusion tomography: a finite-element based algorithm and simulations," Appl. Opt. 37, 5337-5343 (1998). [CrossRef]
  21. Chang, H. L. Graber and R.L. Barbour, "Improved reconstruction algorithm for luminescence optical tomography when background luminophore is present," Appl. Opt. 37, 3547-3552 (1998). [CrossRef]
  22. J. Lee and E. M. Sevick-Muraca, "Lifetime and absorption imaging with fluorescence FDPM," Time-resolved fluorescence spectroscopy and imaging in tissues, E. M. Sevick-Muraca (ed.)., Proc. Soc. Photo-Opt. Instrum. Eng., 3600: (to be published), (1999).
  23. A. Ishimaru, Wave propagation and scattering in random media,( Academic Press, New York, 1978).
  24. M. Schweiger, S. R. Arridge, M. Hiraka and D. T. Delpy, "The finite-element method for the propagation of light in scattering media- boundary and source conditions," Med. Phys. 22, 1779-1792 (1995). [CrossRef] [PubMed]
  25. R. A. J. Groenhuis, H. A. Ferwerda and J. J. Ten Bosch, "Scattering and absorption of turbid material determined from reflection measurements," Appl. Opt. 22, 2456-2462 (1983). [CrossRef] [PubMed]
  26. O. C. Zienkiewcz and R. L. Taylor, The finite element methods in engineering science, (McGraw-Hill, New York, 1989).
  27. L. C. W. Dixon and R. C. Price, "Numerical experience with the truncated Newton method for unconstrained optimization," JOTA 56, 245-255 (1988). [CrossRef]
  28. R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, England,(1996).
  29. R. S. Dembo and T. Steihaug, "Truncated Newton algorithms for large-scale unconstrained optimization," Math Programming 26, 190-212 (1983). [CrossRef]
  30. R. C. Price, Sparse matrix optimization using automatic differentiation, Ph. D. thesis, University of Hertfordshire, U. K., (1987).
  31. L. Armijo "Minimization of functions having Lipschitz continuous first partial derivatives," Pacific J. Mathematics 16, 1-3 (1966).
  32. P. Wolfe, "Convergence condition for ascent method," SIAM Rev., 11 226-253 (1969). [CrossRef]
  33. B., Christianson, A. J., Davies, L. C. W. Dixon, R. Roy and P. van der Zee, "Giving reverse differentiation a helping hand," Opt. Meth. And Software 8, 53-67 (1997). [CrossRef]
  34. A. J., Davies, B. Christianson,, L. C. W. Dixon,, R. Roy and P. van der Zee, "Reverse differentiation and the inverse diffusion problem," Adv. In Eng. Software 28, 217-221 (1997). [CrossRef]
  35. R. E. Wengert, "A simple automatic derivative evaluation program," Comm. A. C. M. 7, 463-464 (1964).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Figure 1
 

Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited