OSA's Digital Library

Optics Express

Optics Express

  • Editor: J. H. Eberly
  • Vol. 2, Iss. 6 — Mar. 16, 1998
  • pp: 213–226
« Show journal navigation

A gradient-based optimisation scheme for optical tomography

Simon R. Arridge and Martin Schweiger  »View Author Affiliations


Optics Express, Vol. 2, Issue 6, pp. 213-226 (1998)
http://dx.doi.org/10.1364/OE.2.000213


View Full Text Article

Acrobat PDF (545 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Optical tomography schemes using non-linear optimisation are usually based on a Newton-like method involving the construction and inversion of a large Jacobian matrix. Although such matrices can be efficiently constructed using a reciprocity principle, their inversion is still computationally difficult. In this paper we demonstrate a simple means to obtain the gradient of the objective function directly, leading to straightforward application of gradient-based optimisation methods.

© Optical Society of America

1. Introduction

Optical tomography is potentially a powerful tool to noninvasively obtain spatially resolved data of the optical parameters of tissue, from which physiologically relevant information such as local oxygenation can be calculated. Primary applications of this new imaging modality are the monitoring of cerebral blood and tissue oxygenation of newborn and preterm infants [1

1. A. D. Edwards, J. S. Wyatt, C. E. Richardson, D. T. Delpy, M. Cope, and E. O. R. Reynolds, “Cotside measurement of cerebral blood flow in ill newborn infants by near infrared spec-troscopy,” Lancet 2, 770–771 (1988). [CrossRef] [PubMed]

,2

2. J. S. Wyatt, M. Cope, D. T. Delpy, C. E. Richardson, A. D. Edwards, S. C. Wray, and E. O. R. Reynolds, “Quantitation of cerebral blood volume in newborn infants by near infrared spectroscopy,” J. Appl. Physiol. 68, 1086–1091 (1990). [PubMed]

] to prevent death or permanent brain damage caused by asphyxiation during birth, functional mapping of brain activation during physical or mental exercise [3

3. M. Tamura, “Multichannel near-infrared optical imaging of human brain activity,” in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC1996) Vol. 2, pp. 8–10.

], and imaging of the breast to detect tumours [4

4. J. C. Hebden, R. A. Kruger, and K. S. Wong, “Time resolved imaging through a highly scattering medium,” Appl. Opt. 30, 788–794 (1991). [CrossRef] [PubMed]

]. Recent developments can be found in several review articles in the recent Royal Society meeting [5

5. Near-infrared spectroscopy and imaging of living systems, special issue of Philos. Trans. R. Soc. London Ser. B, Vol. 352 (1997).

], and other journals [6

6. J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine: I. Experimental techniques,” Phys. Med. Biol. 42, 825–840 (1997). [CrossRef] [PubMed]

,7

7. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42, 841–853 (1997). [CrossRef] [PubMed]

].

The principal problem that arises in optical tomography is the dominance of scattering, which causes light to propagate diffusely in tissue and thus prohibits the application of direct reconstruction methods using the Radon transform. Various reconstruction methods, including ad hoc backprojection [8

8. S. B. Colak, G. W. Hooft, D. G. Papaioannou, and M. B. van der Mark, “3D backprojection tomography for medical optical imaging,” in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC1996) Vol. 2, pp. 294–298.

,9

9. S. A. Walker, S. Fantini, and E. Gratton, “Back-projection reconstructions of cylindrical inho-mogeneities from frequency domain optical measurements in turbid media,” in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC1996) Vol. 2, pp. 137–141.

] and semi-analytic [10

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

] schemes have been proposed, but increasingly attention is turning to iterative, optimisation based reconstruction methods [11–13

11. S. R. Arridge, M. Schweiger, and D. T. Delpy, “Iterative reconstructionof near-infrared absorption images,” in Inverse Problems in Scattering and Imaging, M. A. Fiddy, ed., Proc. SPIE 1767, 372–383 (1992). [CrossRef]

].

In this approach, the problem is regarded as the optimisation of an objective function representing the sum-squared difference of the data to the model, plus additional regularisation terms representing prior knowledge. The most prevalent approach to date has been a Newton-like method such as Levenberg-Marquardt which involves repeated construction and inversion of the problem Jacobian. This method rapidly becomes intractable as the size of the problem domain increases, for example when a fully 3D problem is considered. An alternative is to consider gradient-based algorithms, such as conjugate gradients (CG). In this paper we present a method to derive the gradient efficiently using an adjoint source scheme.

2. Problem definition

Let Ω be the domain under consideration, with surface ∂Ω. We consider S source positions ξ⃗jΩ (j = 1…S) and Mj measurement positions ξ⃗j,i Ω for each source j (i = 1…Mj ), resulting in a total number of measurements M TOT = Σj=1S Mj , with M UNIQM TOT distinct measurement positions. The forward problem is non-linear and is represented by:

F=𝛲[μ,κ]
(1)

Under the usual assumptions of a system corrupted by multivariate Gaussian noise, and a maximum-likelihood approach to the solution, we may define the inverse problem as the optimisation of an objective function

Ψ=12j=1Si=1Mj(yj,i𝛲j,i[μ,κ]σj,i)2
(2)

which we write in vector form as

Ψ=12(yF)TR2(yF)=12bTb
(3)

where yj,i is the data for the ith measurement from source j with standard deviation σj,i and bj,i = σj,i1 (yj,i - Fj,i ) is the residual data for this measurement. We use subscripted vectors y⃗j ,F⃗j ,σ⃗j to represent all the relevant quantities for a single source j.

Thus Pj [μ, κ] is the projection operator for source j, F⃗j is the projection data obtained by sampling Pj at the discrete measurement positions {ξj,1,ξj,2ξj,M} and b⃗j = Rj1 (y⃗j - F⃗j ). Bold vectors y⃗, P⃗, F⃗, b⃗ represent the combined quantities from all sources. R is the data-space correlation matrix, which we here take to be

R=diagR1R2RS=diag(σ1,1,σ1,2,σj,i,σS,Ms)
(4)

3. Transport model

We assume here that the propagation model is the diffusion approximation to the radiative transfer equation which in the frequency domain is:

κ(r)Φ̂rω+μ(r)Φ̂(r,ω)+ιωcΦ̂rω=q̂0rt,
(5)

and in the time domain :

κ(r)Φrt+μ(r)Φ(r,t)+1cΦ(r,t)t=q0rt,
(6)

where Φ is the isotropic photon density, q 0 is an isotropic source distribution and c is the speed of light in the medium. For a detailed discussion of the physical and mathematical models employed in this field, refer to the recent review articles by Hebden and Arridge [6

6. J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine: I. Experimental techniques,” Phys. Med. Biol. 42, 825–840 (1997). [CrossRef] [PubMed]

,7

7. S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42, 841–853 (1997). [CrossRef] [PubMed]

]. The model is characterised by the two spatially varying functions μ(r⃗) (absorption) 1 and κ(r⃗) (diffusion), which gives rise to the dual search-space nature of the optimisation problem defined in Eq. 2.

Γ(ξ)=(ξ)n̂Φ(ξ),
(7)

where n̂ is the outer normal of Ω at ξ. We use the Robin-type boundary condition

Φ(ξ)+καn̂Φ(ξ)=0,
(8)

where α is a term to incorporate boundary reflections as a result of a refractive index mismatch at Ω [16

16. J. D. Moulton, Diffusion modelling of picosecond laser pulse propagation in turbid media, M. Eng. thesis, McMaster University, Hamilton, Ontario (1990).

,18

18. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

]. A collimated source incident at ζ ∈ Ω is commonly represented by a diffuse point source q 0(r⃗) = δ(r⃗-r⃗s ) where r⃗s is located at a depth of one scattering length below the surface.

4. Finite element approach

Although the above approach is general, and can be applied to analytical forward models, the most successful approaches use a numerical approach to solve the Forward problem. Here we assume it is an FEM approach.

The domain Ω is divided into P elements, joined at D vertex nodes. The solution Φ is approximated by the piecewise linear function Φh(r⃗) = ΣiD Φi ui (r⃗) ∈ ʊh, where ʊh is a finite dimensional subspace spanned by basis functions ui (r⃗), i = 1…D chosen to have limited support. The problem of solving for Φh becomes one of sparse matrix inversion for which standard methods such as Cholesky decomposition or conjugate gradient solvers are readily available. The advantage of the FEM approach is its versatility which makes it applicable to complex geometries and highly inhomogeneous parameter distributions.

As developed in [17

17. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modelling photon transport in tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

,18

18. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

], the diffusion equation in the FEM framework is expressed, in the frequency domain as :

(K(κ)+C(μ)+αA+ιωB)Φ(ω)=Q(ω)
(9)

and in the time-domain as :

(K(κ)+C(μ)+αA)Φ(t)+BΦ(t)t=Q(t)
(10)

where the system matrices K,C, A and B have entries given by :

Kij=Ωκ(r)ui(r)uj(r)dnr
(11)
Cij=Ωμ(r)ui(r)uj(r)dnr
(12)
Bij=1cΩui(r)uj(r)dnr
(13)
Aij=Ωui(r)uj(r)d(Ω)
(14)

To relate the FEM approach to the the forward model, we define :

Fj=𝚳[Φj]
(15)

where Φ⃗j is the solution to Eq. 9 or 10 for the jth source, and M : L2(Ω) → L2(Ω) is a measurement operator.

5. Measurement types

The simplest measurement operator to consider is the boundary operator B which implements Eq. 7. In the discrete case this is a Mj × D matrix. Since B is linear it can be applied to either temporal data Φ(t) or frequency domain data Φ(ω). We define the integrated intensity measurement ε to represent both ε[Φ(ω)] = B[Φ(ω = 0)] (i.e. the solution to the DC steady problem) or ε[Φ(t)] = B[∫ Φ(t)dt] (the sum over time of all arriving photons).

In our approach we have advocated the use of measurement operators that are self-normalised

𝚳̅=τε
(16)

where τ is the integral of Γ(t) weighted by a temporal filter. Examples that have been suggested are :

time-gatedintensity:𝚳E¯(T)[Φ(t)]=1ε[Φ(t)]0TB[Φ(t)]dt,
(17)
n-thtemporalmoment:𝚳tn[Φ(t)]=1ε[Φ(t)]0tnB[Φ(t)]dt,
(18)
n-thcentralmoment:𝚳cn[Φ(t)]=1ε[Φ(t)]0(tt)nB[Φ(t)]dt,
(19)
normalisedLaplacetransform:𝚳L̅(s)[Φ(t)]=1ε[Φ(t)]0estB[Φ(t)]dt.
(20)

Additional features could be chosen, such as the logarithmic slope of the temporal decay of Γ, the peak intensity, etc. However, we aim to chose features that satisfy two conditions: robustness of the experimental measurements with respect to systematic errors such as fluctuations of the source power, detector sensitivity, or fibre coupling losses, and secondly the ability of the forward model to generate the corresponding data efficiently. Both conditions are satisfied for measurement types (18), (19) and (20), since these are normalised and therefore not dependent on absolute intensity measurements, and they can be calculated directly by our forward model without the need for explicitly generating the temporal profile of Γ(t) [22

22. S. R. Arridge and M. Schweiger, “Direct calculation of the moments of the distribution of photon time of flight in tissues with a finite-element method,” Appl. Opt. 34, 2683–2687 (1995). [CrossRef] [PubMed]

,23

23. M. Schweiger and S. R. Arridge, “Direct calculation of the Laplace transform of the distribution of photon time of flight in tissue with a finite-element method,” Appl. Opt. 36, 9042–9049 (1997). [CrossRef]

].

6. Reconstruction methods

To solve the optimisation problem, we consider the gradient of the objective function. Letting x represent either μ or κ the kth component of the gradient is given by :

Ψxk=j=1Si=1Mj(yj,iPj,i[μ,κ]σj,i2)(Pj,i[μ,κ]xk)
(21)

leading to the vector form

z=j=1SP'jTRj1bj=P'TR1b
(22)
=j=1SJjTbj=JTb
(23)

The Jacobian 𝖩 of the forward problem is a M TOT × N TOT matrix where M TOT = Σ Mj and N TOT is the number of coefficients of μ and κ expressed in some basis. In this paper we take the basis to be the element values, although other choices are possible. 𝖩 has the structure

b1b2...bS=J1,(μ)J1,(κ)J2,(μ)J2,(κ)......JS,(μ)JS,(κ)ΔμΔκ
(24)

where we define 𝖩j as the Mj × N TOT matrix that is the Jacobian for the objective function corresponding to source j and 𝖩j,(μ),𝖩j,(κ) as the sub-matrices corresponding to the two variables μ, κ.

The Jacobian may be found in a row-by-row fashion if we establish the notation that the vector corresponding to the i th row of 𝖩j is the Photon Measurement Density Function PM→DF(j,i)T, as defined in [19

19. S. R. Arridge, “Photon measurement density functions. Part 1: Analytic forms,” Appl. Opt. 34, 7395–7409 (1995). [CrossRef] [PubMed]

,20

20. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite element calculations,” Appl. Opt. 34, 8026–8037 (1995). [CrossRef] [PubMed]

] which we represent :

Jji=PMDF(j,i)T=[PMDF(j,i),(μ)TPMDF(j,i),(κ)T]
(25)

Algorithm 1 Nonlinear conjugate-gradient method

Set residual p⃗ (0) = d⃗ (0)

Define termination criterion ε

Set iteration counter n = 0

repeat

Find γ (n) that minimises Ψ(x⃗ (n) + γ (n) d⃗ (n))

β(n+1)=max(p(n+1)T(p(n+1)p(n))p(n)Tp(n),0)

d⃗ (n+1) = p⃗ (n+1) + β (n+1) d⃗ (n)

n = n + 1

untilz⃗(x⃗ (n))∥ < ε

7. Gradient calculation

If the gradient is found by the explicit form of Eq. 23 then no advantage is gained. The use of an adjoint scheme to calculate z⃗ has been discussed previously for the time-domain diffusion model using a Finite Difference Scheme [26

26. S. S. Saquib, K. M. Hanson, and G. S. Cunningham, “Model-based image reconstruction from time-resolved diffusion data,” in Medical Imaging: Image Processing, K. M. Hanson, ed., Proc SPIE 3034, 369–380 (1997). [CrossRef]

], the time-domain transport model using a Finite Difference Scheme [27

27. O. Dorn, Das inverse Transportproblem in der Lasertomographie, Ph. D. thesis, University of Münster, 1997.

], and the steady-state DC case using a Finite Element scheme [28

28. R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, 1997.

]. A summary of these approaches and their relation to other reconstruction methods is given in [29

29. S. R. Arridge and M. Schweiger, “A general framework for iterative reconstruction algorithms in optical tomography, using a finite element method,” in Computational Radiology and Imaging: Therapy and DiagnosisC. Borgers and F. Natterer, eds., IMA Volumes in Mathematics and its Applications (Springer1998, in press).

]. Here we give a derivation for the normalised measurement operators that we use in our application.

Consider first the frequency domain problem. We have that

z=JTb
(26)
=j=1SJjTbj
(27)
=j=1Sj=1MjPMDF(j,i)Tbj,i
(28)

If we define the Adjoint problem

(K+C+αAιωB)Φm+(ω)=Qi+(ω)
(29)

where Qi+(ω) is the adjoint source [20

20. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite element calculations,” Appl. Opt. 34, 8026–8037 (1995). [CrossRef] [PubMed]

], then we have that PMD⃗F(j,i)T is a vector in two parts, whose kth components are given by σj,i1 Φi+ (ω) × Φ⃗j (ω) and σj,i1 ∇Φ⃗j (ω) × ∇Φ⃗j (ω) respectively. Here the notation a⃗ × b⃗ as defined in [20

20. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite element calculations,” Appl. Opt. 34, 8026–8037 (1995). [CrossRef] [PubMed]

] is taken to mean the vector of samples of the product of the fields a⃗, b⃗ and ∇a⃗ × ∇b⃗ is the vector of samples of the scalar product of the fields ∇a ∙ ∇b. We may now write

z(μ)=j=1Si=1Mjbj,iσj,iΦi+(ω)×Φj(ω)
(30)
z(κ)=j=1Si=1Mjbj,iσj,iΦi+(ω)×Φj(ω)
(31)

But Φm+ (ω) is the solution to Eq. 29, so we form a vector

νj+(ω)=i=1Mjbj,iσj,iQi+(ω)
(32)

and solve

(K+C+αAιωB)ηj+(ω)=νj+(ω)
(33)

leading to:

z(μ)=j=1Sηj+(ω)×Φj(ω)
(34)
z(κ)=j=1Sηj+(ω)×Φj(ω)
(35)

In the following we extend the notation of §2. by adding a superscript to indicate the measurement type. Thus PjΜ is the projection operator for measurement Fj,iΜ, and bj,iΜ, σj,iΜ are the forward data, residual data, and standard deviation of the error for the ith measurement from source j of type 𝚳. To find the gradient for the normalised measurement types, we use the result from[20

20. S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite element calculations,” Appl. Opt. 34, 8026–8037 (1995). [CrossRef] [PubMed]

] that the Jacobian for 𝚳¯ can be constructed from

Jji𝚳̅=1Fj,iε(Jji𝚳Fj,i𝚳̅Jjiε)
(36)

from which we can write

J𝚳̅Tb𝚳̅=jSiMjJji𝚳̅bj,i𝚳̅
(37)
=jSiMj1Fj,iε(Jji𝚳Fj,i𝚳̅Jjiε)bj,i𝚳̅
(38)

which leads to the gradient :

z(μ)𝚳̅=jSiMjbj,i𝚳̅σj,i𝚳̅Fj,iε(τ[Φi+×Φj]Fj,i𝚳̅Φi+×Φj)
(39)
z(κ)𝚳̅=jSiMjbj,i𝚳̅σj,i𝚳̅Fj,iε(τ[Φi+×Φj]Fj,i𝚳̅Φi+×Φj)
(40)

We therefore need to construct two adjoint sources :

νj𝚳̅+(0)=i=1Mjbj,i𝚳̅Fj,i𝚳̅σj,i𝚳¯Fj,iεQi+
(41)
νj𝚳̅+(1)=i=1Mjbj,i𝚳̅σj,i𝚳¯Fj,iεQi+
(42)

From here we proceed as before, finding two solutions η⃗j 𝚳¯+ (0), η⃗j 𝚳¯+ (1) and computing the gradient as :

z(μ)𝚳̅=jSτ[ηj𝚳̅+(1)×Φj]ηj𝚳¯+(0)×Φj
(43)
z(κ)𝚳̅=jSτ[ηj𝚳̅+(1)×Φj]ηj𝚳¯+(0)×Φj
(44)

This method is applicable for measurement types (18) and (20). The nth central moment given by (19) can be constructed as a linear combination of the first n temporal moments and so the gradient requires n adjoint sources, constructed in a straightforward manner.

8. Results

We present simultaneous reconstructions of μa and μs using the nonlinear conjugate gradient algorithm in conjunction with the adjoint source scheme for two different test cases: a circular object with three absorbing and/or scattering objects embedded in a homogeneous background, and the model of a slice of the head of a neonate, with a complex distribution of different tissue types and embedded lesions.

The forward data used in the reconstruction of these test cases were simulated by the FEM light transport model, using highly resolved meshes to ensure good quality of the data. Coarser meshes were used for the inverse solution to speed up computation times. For all following reconstructions we used a combination of skew and Laplace transform data. We have argued previously that this yields the best results among the data types listed above [31

31. M. Schweiger and S. R. Arridge, “Optimal data types in optical tomography,” in Information Processing in Medical Imaging, Lect. Notes Comput. Sci. , 1230, 71–84 (1997). [CrossRef]

].

8.1 Circular test object

The test object was a circle of radius 25 mm and three embedded objects, each with a radius of 1.5 mm. The optical parameters of the background medium and the objects are listed in Table 1. The two images in the left column of Fig. 1 show the locations of the perturbations.

Forward data were calculated for 32 source and 32 detector positions, equally spaced along the circumference. Noise was simulated with a noise model presented previously [30

30. S. R. Arridge, M. Hiraoka, and M. Schweiger, “Statistical basis for the determination of optical pathlength in tissue,” Phys. Med. Biol. 40, 1539–1558 (1995). [CrossRef] [PubMed]

][31

31. M. Schweiger and S. R. Arridge, “Optimal data types in optical tomography,” in Information Processing in Medical Imaging, Lect. Notes Comput. Sci. , 1230, 71–84 (1997). [CrossRef]

], assuming that 104 photons were collected for each measurement. The meshes for the forward and inverse solvers contained 7200 and 2808 triangular elements, respectively. Figure 1 shows the target images, the reconstruction results obtained with the CG solver, and, for comparison, reconstructions obtained with a steepest descent and an ART (algebraic reconstruction technique) solver. All reconstructions started from the homogeneous background, and were terminated after 100 iterations.The ART solver used was Block-ART with random source ordering as described in [29

29. S. R. Arridge and M. Schweiger, “A general framework for iterative reconstruction algorithms in optical tomography, using a finite element method,” in Computational Radiology and Imaging: Therapy and DiagnosisC. Borgers and F. Natterer, eds., IMA Volumes in Mathematics and its Applications (Springer1998, in press).

].

All three algorithms manage to distinguish between absorption and scatter features with relatively little cross-talk, with the exception of the ART algorithm which shows a strong artefact in the absorption image in the position of the scattering feature. None of the methods recovers the absolute values of the perturbations very well, but the CG and ART methods obtain higher peak values than the steepest descent method after a given number of iterations.

Table 1. Optical parameters of circular test object.

table-icon
View This Table
| View All Tables
Figure 1. Target image (col. 1), reconstruction after 100 iterations with conjugate gradient method (col. 2), with steepest descent method (col. 3) and with ART method (col. 4). In all cases, the top image is μa , and the bottom image is μs . The animations linked to the figure show the first 50 reconstruction iterations, and the final 50 in steps of 10. [Media 1] [Media 2] [Media 3] [Media 4] [Media 5] [Media 6]

The L2 norms of the data errors, b2, for the three algorithms are shown in Fig. 2, as a function of the iteration number (left graph) and as a function of the reconstruction runtime (right graph). We find that ART performs better than the CG method at the beginning of the reconstruction, but levels out very quickly. After 100 iterations both reconstructions arrive approximately at the same error norm. With respect to runtime, the CG solver performs better than ART, because a single ART iteration is computationally more expensive than a CG iteration.

Figure 2. L2 data norms as a function of iteration number, for different algorithms, as a function of iteration number (left) and of runtime (right).

Figure 3 shows the L2 norms of the reconstructed images, (x⃗ - x⃗ target)2, for both absorption and scatter as a function of the iteration number. The results are similar to the data errors. Again we find that ART performs well during the first few iterations but quickly levels out, in particular for μs . For higher iterations the CG algorithm becomes superior. The steepest descent method performs poorly for μa and μs .

Figure 3. L2 solution norms as a function of iteration number, for different algorithms. Left: absorption, right: scatter.

8.2 Neonatal head model

To investigate a more complex case we use a mesh that simulates the parameter distribution in a cross-section of a neonatal head. The outline and internal boundaries between tissue types were derived from an MRI image. The sagittal diameter is 80 mm. We distinguish four tissue types and assigned literature values for their optical properties (Table 2). In addition we placed absorbing and scattering perturbations in random locations to simulate lesions or similar high-contrast features.

Skew and Laplace transform forward data for 32 sources and 32 detectors, equidistantly spaced at the boundary, were generated with a high-resolution mesh (41500 elements). For the reconstruction a coarser mesh of 10500 elements was used. The reconstruction was started from a homogeneous distribution, corresponding to the skin parameters (μa = 0.022, μs = 1). This is likely to be a worst-case scenario, since in practice it may be possible to assume more prior knowledge.

Table 2. Optical parameters of neonatal head model.

table-icon
View This Table
| View All Tables
Figure 4. Target image (left col.), reconstruction after 100 iterations with conjugate gradient method (right col.). Top row is μa , bottom row is μs . The animations linked to the figure show the first 50 reconstruction iterations, and the final 50 in steps of 10. [Media 7] [Media 8]

9. Conclusions

Often it is thought that iterative optimisation based approaches to image reconstruction poses an insurmountable computational burden, especially when the fully three-dimensional and time-dependent problem is addressed. In this paper we have presented a method for deriving the gradient of the objective function in a Finite Element Method based scheme for optical tomography. We presented only the conjugate gradient optimisation method. The use of higher order methods such as Truncated-Newton, is being investigated as well as the choice of appropriate regularisation schemes.

We have presented simultaneous reconstructions of absorption and scatter from simulated data for two cases: a simple homogeneous circle with embedded absorbing and scattering features, and a complex neonatal head model with an inhomogeneous background and several high-contrast lesions. We have presented a comparison between conjugate gradient (CG), steepest descent and ART implementations and shown that after the first iterations, where ART is favourable, CG performs equivalently or even better, if runtime is considered.

We have also demonstrated that the reconstruction of complex problems such as the neonatal head model in the second example, is fundamentally more difficult than the often used blobs in a homogeneous background, and we suggest that such complex cases must be considered in evaluating reconstruction algorithms that are to be employed in clinical imaging problems.

Acknowledgements

We would like to thank Ken Hanson of Los Alamos National Laboratory, Frank Natterer of Münster University, and Ranadhir Roy of University College London for many useful discussions concerning adjoint methods in optimisation. This work was supported by funding from Action Research, and the Wellcome Trust.

Footnotes

1In biomedical optics the absorption parameter is usually denoted μ a, with the scattering parameter denoted μ s′, and the diffusion parameter given by 13(μa+μ's). We will use μ to denote f a throughout, to avoid a proliferation of subscripts.

References

1.

A. D. Edwards, J. S. Wyatt, C. E. Richardson, D. T. Delpy, M. Cope, and E. O. R. Reynolds, “Cotside measurement of cerebral blood flow in ill newborn infants by near infrared spec-troscopy,” Lancet 2, 770–771 (1988). [CrossRef] [PubMed]

2.

J. S. Wyatt, M. Cope, D. T. Delpy, C. E. Richardson, A. D. Edwards, S. C. Wray, and E. O. R. Reynolds, “Quantitation of cerebral blood volume in newborn infants by near infrared spectroscopy,” J. Appl. Physiol. 68, 1086–1091 (1990). [PubMed]

3.

M. Tamura, “Multichannel near-infrared optical imaging of human brain activity,” in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC1996) Vol. 2, pp. 8–10.

4.

J. C. Hebden, R. A. Kruger, and K. S. Wong, “Time resolved imaging through a highly scattering medium,” Appl. Opt. 30, 788–794 (1991). [CrossRef] [PubMed]

5.

Near-infrared spectroscopy and imaging of living systems, special issue of Philos. Trans. R. Soc. London Ser. B, Vol. 352 (1997).

6.

J. C. Hebden, S. R. Arridge, and D. T. Delpy, “Optical imaging in medicine: I. Experimental techniques,” Phys. Med. Biol. 42, 825–840 (1997). [CrossRef] [PubMed]

7.

S. R. Arridge and J. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol. 42, 841–853 (1997). [CrossRef] [PubMed]

8.

S. B. Colak, G. W. Hooft, D. G. Papaioannou, and M. B. van der Mark, “3D backprojection tomography for medical optical imaging,” in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC1996) Vol. 2, pp. 294–298.

9.

S. A. Walker, S. Fantini, and E. Gratton, “Back-projection reconstructions of cylindrical inho-mogeneities from frequency domain optical measurements in turbid media,” in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC1996) Vol. 2, pp. 137–141.

10.

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

11.

S. R. Arridge, M. Schweiger, and D. T. Delpy, “Iterative reconstructionof near-infrared absorption images,” in Inverse Problems in Scattering and Imaging, M. A. Fiddy, ed., Proc. SPIE 1767, 372–383 (1992). [CrossRef]

12.

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 (1995). [CrossRef]

13.

B. W. Pogue, M. S. Patterson, H. Jiang, and K. D. Paulsen, “Initial assessment of a simple system for frequency domain diffuse optical tomography,” Phys. Med. Biol. 40, 1709–1729 (1995). [CrossRef] [PubMed]

14.

D. T. Delpy, M. Cope, P. van der Zee, S. R. Arridge, S. Wray, and J. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33, 1433–1442 (1988). [CrossRef] [PubMed]

15.

B. Chance, M. Maris, J. Sorge, and M. Z. Zhang, “A phase modulation system for dual wavelength difference spectroscopy of haemoglobin deoxygenation in tissue,” in Time-resolved Laser Spectroscopy in Biochemistry II, J. R. Lakowicz, ed., Proc. SPIE 1204, 481–491 (1990). [CrossRef]

16.

J. D. Moulton, Diffusion modelling of picosecond laser pulse propagation in turbid media, M. Eng. thesis, McMaster University, Hamilton, Ontario (1990).

17.

S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modelling photon transport in tissue,” Med. Phys. 20, 299–309 (1993). [CrossRef] [PubMed]

18.

M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element model for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. 22, 1779–1792 (1995). [CrossRef] [PubMed]

19.

S. R. Arridge, “Photon measurement density functions. Part 1: Analytic forms,” Appl. Opt. 34, 7395–7409 (1995). [CrossRef] [PubMed]

20.

S. R. Arridge and M. Schweiger, “Photon measurement density functions. Part 2: Finite element calculations,” Appl. Opt. 34, 8026–8037 (1995). [CrossRef] [PubMed]

21.

M. S. Bazaraa, H. D. Sherali, and C. M. Shetty, Nonlinear Programming: Theory and Algorithms, 2nd edition (Wiley, New York, 1993).

22.

S. R. Arridge and M. Schweiger, “Direct calculation of the moments of the distribution of photon time of flight in tissues with a finite-element method,” Appl. Opt. 34, 2683–2687 (1995). [CrossRef] [PubMed]

23.

M. Schweiger and S. R. Arridge, “Direct calculation of the Laplace transform of the distribution of photon time of flight in tissue with a finite-element method,” Appl. Opt. 36, 9042–9049 (1997). [CrossRef]

24.

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 Imag. Vision 3, 263–283 (1993). [CrossRef]

25.

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]

26.

S. S. Saquib, K. M. Hanson, and G. S. Cunningham, “Model-based image reconstruction from time-resolved diffusion data,” in Medical Imaging: Image Processing, K. M. Hanson, ed., Proc SPIE 3034, 369–380 (1997). [CrossRef]

27.

O. Dorn, Das inverse Transportproblem in der Lasertomographie, Ph. D. thesis, University of Münster, 1997.

28.

R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, University of Hertfordshire, 1997.

29.

S. R. Arridge and M. Schweiger, “A general framework for iterative reconstruction algorithms in optical tomography, using a finite element method,” in Computational Radiology and Imaging: Therapy and DiagnosisC. Borgers and F. Natterer, eds., IMA Volumes in Mathematics and its Applications (Springer1998, in press).

30.

S. R. Arridge, M. Hiraoka, and M. Schweiger, “Statistical basis for the determination of optical pathlength in tissue,” Phys. Med. Biol. 40, 1539–1558 (1995). [CrossRef] [PubMed]

31.

M. Schweiger and S. R. Arridge, “Optimal data types in optical tomography,” in Information Processing in Medical Imaging, Lect. Notes Comput. Sci. , 1230, 71–84 (1997). [CrossRef]

OCIS Codes
(100.3190) Image processing : Inverse problems
(110.6960) Imaging systems : Tomography
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.6920) Medical optics and biotechnology : Time-resolved imaging
(170.6960) Medical optics and biotechnology : Tomography

ToC Category:
Focus Issue: Tomographic image reconstruction

History
Original Manuscript: January 20, 1998
Published: March 16, 1998

Citation
Simon Arridge and Martin Schweiger, "A gradient-based optimisation scheme foroptical tomography," Opt. Express 2, 213-226 (1998)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-2-6-213


Sort:  Journal  |  Reset  

References

  1. A. D. Edwards, J. S. Wyatt, C. E. Richardson, D. T. Delpy, M. Cope and E. O. R. Reynolds, \Cotside measurement of cerebral blood ow in ill newborn infants by near infrared spec- troscopy," Lancet 2, 770-771 (1988). [CrossRef] [PubMed]
  2. S. Wyatt, M. Cope, D. T. Delpy, C. E. Richardson, A. D. Edwards, S. C. Wray and E. O. R. Reynolds, \Quantitation of cerebral blood volume in newborn infants by near infrared spectroscopy," J. Appl. Physiol. 68, 1086-1091 (1990). [PubMed]
  3. M. Tamura, \Multichannel near-infrared optical imaging of human brain activity," in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC 1996) Vol. 2, pp. 8-10.
  4. J. C. Hebden, R. A. Kruger and K. S. Wong, \Time resolved imaging through a highly scattering medium," Appl. Opt. 30, 788-794 (1991). [CrossRef] [PubMed]
  5. Near-infrared spectroscopy and imaging of living systems, special issue of Philos. Trans. R. Soc. London Ser. B, Vol. 352 (1997).
  6. J. C. Hebden, S. R. Arridge and D. T. Delpy, \Optical imaging in medicine: I. Experimental techniques," Phys. Med. Biol. 42, 825-840 (1997). [CrossRef] [PubMed]
  7. S. R. Arridge and J. C. Hebden, \Optical imaging in medicine: II. Modelling and reconstruction," Phys. Med. Biol. 42, 841-853 (1997). [CrossRef] [PubMed]
  8. S. B. Colak, G. W. Hooft, D. G. Papaioannou and M. B. van der Mark, \3D backprojection tomography for medical optical imaging," in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC 1996) Vol. 2, pp. 294-298.
  9. S. A. Walker, S. Fantini and E. Gratton, \Back-projection reconstructions of cylindrical inho- mogeneities from frequency domain optical measurements in turbid media," in OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration, R. R. Alfano and J. G. Fujimoto, eds. (Optical Society of America, Washington, DC 1996) Vol. 2, pp. 137-141.
  10. M. A. O Leary, D. A. Boas, B. Chance and A. G. Yodh, \Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography," Opt. Lett. 20, 426-428 (1995). [CrossRef]
  11. S. R. Arridge, M. Schweiger and D. T. Delpy, \Iterative reconstruction of near-infrared absorption images," in Inverse Problems in Scattering and Imaging, M. A. Fiddy, ed., Proc. SPIE 1767, 372-383 (1992). [CrossRef]
  12. 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 (1995). [CrossRef]
  13. B. W. Pogue, M. S. Patterson, H. Jiang and K. D. Paulsen, \Initial assessment of a simple system for frequency domain diffuse optical tomography," Phys. Med. Biol. 40, 1709-1729 (1995). [CrossRef] [PubMed]
  14. D. T. Delpy, M. Cope, P. van der Zee, S. R. Arridge, S. Wray and J. Wyatt, \Estimation of optical pathlength through tissue from direct time of ight measurement," Phys. Med. Biol. 33, 1433-1442 (1988). [CrossRef] [PubMed]
  15. B. Chance, M. Maris, J. Sorge and M. Z. Zhang, \A phase modulation system for dual wave- length difference spectroscopy of haemoglobin deoxygenation in tissue," in Time-resolved Laser Spectroscopy in Biochemistry II, J. R. Lakowicz, ed., Proc. SPIE 1204, 481-491 (1990). [CrossRef]
  16. J. D. Moulton, Diffusion modelling of picosecond laser pulse propagation in turbid media, M. Eng. thesis, McMaster University, Hamilton, Ontario (1990).
  17. S. R. Arridge, M. Schweiger, M. Hiraoka and D. T. Delpy, \A finite element approach for modelling photon transport in tissue," Med. Phys. 20, 299-309 (1993). [CrossRef] [PubMed]
  18. M. Schweiger, S. R. Arridge, M. Hiraoka and D. T. Delpy, \The finite element model for the propagation of light in scattering media: Boundary and source conditions," Med. Phys. 22, 1779-1792 (1995). [CrossRef] [PubMed]
  19. S. R. Arridge, \Photon measurement density functions. Part 1: Analytic forms," Appl. Opt. 34, 7395-7409 (1995). [CrossRef] [PubMed]
  20. S. R. Arridge and M. Schweiger, \Photon measurement density functions. Part 2: Finite element calculations," Appl. Opt. 34, 8026-8037 (1995). [CrossRef] [PubMed]
  21. M. S. Bazaraa, H. D. Sherali and C. M. Shetty, Nonlinear Programming: Theory and Algorithms, 2nd edition (Wiley, New York, 1993).
  22. S. R. Arridge and M. Schweiger, \Direct calculation of the moments of the distribution of photon time of ight in tissues with a finite-element method," Appl. Opt. 34, 2683-2687 (1995). [CrossRef] [PubMed]
  23. M. Schweiger and S. R. Arridge, \Direct calculation of the Laplace transform of the distribution of photon time of ight in tissue with a finite-element method," Appl. Opt. 36, 9042-9049 (1997). [CrossRef]
  24. 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 Imag. Vision 3, 263-283 (1993). [CrossRef]
  25. 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]
  26. S. S. Saquib, K. M. Hanson and G. S. Cunningham, \Model-based image reconstruction from time-resolved diffusion data," in Medical Imaging: Image Processing, K. M. Hanson, ed., Proc SPIE 3034, 369-380 (1997). [CrossRef]
  27. O. Dorn, Das inverse Transportproblem in der Lasertomographie, Ph. D. thesis, University of M"unster, 1997.
  28. R. Roy, Image reconstruction from light measurements on biological tissue, Ph. D. thesis, Uni- versity of Hertfordshire, 1997.
  29. S. R. Arridge and M. Schweiger, \A general framework for iterative reconstruction algorithms in optical tomography, using a finite element method," in Computational Radiology and Imaging: Therapy and Diagnosis C. Borgers and F. Natterer, eds., IMA Volumes in Mathematics and its Applications (Springer 1998, in press).
  30. S. R. Arridge, M. Hiraoka and M. Schweiger, \Statistical basis for the determination of optical pathlength in tissue," Phys. Med. Biol. 40, 1539-1558 (1995). [CrossRef] [PubMed]
  31. M. Schweiger and S. R. Arridge, \Optimal data types in optical tomography," in Information Processing in Medical Imaging, Lect. Notes Comput. Sci., 1230, 71-84 (1997). [CrossRef]

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.

Supplementary Material


» Media 1: MOV (415 KB)     
» Media 2: MOV (611 KB)     
» Media 3: MOV (279 KB)     
» Media 4: MOV (602 KB)     
» Media 5: MOV (744 KB)     
» Media 6: MOV (485 KB)     
» Media 7: MOV (880 KB)     
» Media 8: MOV (436 KB)     

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited