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

© 1998 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,2] to prevent death or permanent brain damage caused by asphyxiation during birth, functional mapping of brain activation during physical or mental exercise [3], and imaging of the breast to detect tumours [4]. Recent developments can be found in several review articles in the recent Royal Society meeting [5], and other journals [6,7].

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,9] and semi-analytic [10] schemes have been proposed, but increasingly attention is turning to iterative, optimisation based reconstruction methods [11–13].

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 *M*_{j}
measurement positions *ξ⃗*_{j,i}
*∂*Ω for each source *j* (*i* = 1…*M*_{j}
), resulting in a total number of measurements *M*
_{TOT} = ${\mathrm{\Sigma}}_{j=1}^{S}$
*M*_{j}
, with *M*
_{UNIQ} ≤ *M*
_{TOT} distinct measurement positions. The forward problem is non-linear and is represented by:

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

which we write in vector form as

where *y*_{j,i}
is the data for the *i*^{th}
measurement from source *j* with standard deviation *σ*_{j,i}
and *b*_{j,i}
= ${\sigma}_{j\mathit{,}i}^{-1}$ (*y*_{j,i}
- *F*_{j,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 *P*_{j}
[*μ, κ*] is the *projection operator* for source *j*, *F⃗*_{j}
is the projection data obtained by sampling *P*_{j}
at the discrete measurement positions {*ξ*⃗_{j,1},*ξ*⃗_{j,2}…*ξ*⃗_{j,M}} and *b⃗*_{j}
= ${\mathrm{R}}_{j}^{-1}$ (*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

## 3. Transport model

Data acquisition for optical tomography is not limited to steady-state attenuation measurements. Time-of-flight systems using an ultra-short pulsed laser as a light source and time-resolved detectors allow the boundary measurement of the temporal intensity response function with a resolution of a few picoseconds [14]. Frequency domain systems use a radio frequency modulated light source and measure the phase shift and modulation amplitude of the transmitted light [15]. While the data acquisition equipment for the steady-state, temporal and frequency domain are quite different, the methods are related from the theoretical viewpoint: the time-of-flight method includes the steady-state case via the integral over time of the detected light, and is related to the frequency domain case via its Fourier transform.

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

and in the time domain :

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,7]. 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.

For both the frequency-domain and time-domain cases, the boundary measurement Γ(*ξ*) at *ξ* ∈ *∂*Ω is related to Φ(*r⃗*) by

where $\overrightarrow{\hat{n}}$ is the outer normal of *∂*Ω at *ξ*. We use the Robin-type boundary condition

where *α* is a term to incorporate boundary reflections as a result of a refractive index mismatch at *∂*Ω [16,18]. 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⃗*) = ${\mathrm{\Sigma}}_{i}^{D}$ Φ_{i}
*u*_{i}
(*r⃗*) ∈ ʊ^{h}, where ʊ^{h} is a finite dimensional subspace spanned by basis functions *u*_{i}
(*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,18], the diffusion equation in the FEM framework is expressed, in the frequency domain as :

and in the time-domain as :

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

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

where Φ⃗_{j} is the solution to Eq. 9 or 10 for the *j*^{th}
source, and M : L_{2}(Ω) → L_{2}(*∂*Ω) 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 *M*_{j}
× *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*

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

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,23].

## 6. Reconstruction methods

To solve the optimisation problem, we consider the *gradient* of the objective function. Letting *x* represent either *μ* or *κ* the *k*^{th}
component of the gradient is given by :

leading to the vector form

The Jacobian 𝖩 of the forward problem is a *M*
_{TOT} × *N*
_{TOT} matrix where *M*
_{TOT} = Σ *M*_{j}
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

where we define 𝖩_{j} as the *M*_{j}
× *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→D${\mathrm{F}}_{(j\mathit{,}i)}^{\mathrm{T}}$, as defined in [19,20] which we represent :

Previously reported reconstruction algorithms are of Newton-type, a prototypical example being the Levenberg-Marquardt algorithm [24,25]. The feature of such methods is that 𝖩 needs to be explicitly created, and repeatedly inverted. Whenever the number of parameters of the optimisation problem is large, explicit calculation and inversion of the linearised step is regarded as intractable. To overcome this difficulty we have previously reported the use of ART (Algebraic Reconstruction Technique) methods that operate in a row-by-row fashion using only the implicit representation of 𝖩 given by Eq. 25. However these methods still have a storage overhead that is quite costly [29]. Instead, gradient methods require only the computation of the gradient *z⃗* defined in Eq. 22, the proto-typical example being conjugate gradients. In this method a set of conjugate search directions is generated to find the minimum of the objective function. At each iteration step a one-dimensional line minimisation along the current search direction is performed. Conjugate gradient methods are well-established in nonlinear optimisation. See for example [21]. A brief outline is given in Algorithm 1.

**Algorithm 1** Nonlinear conjugate-gradient method

Set search direction *d⃗*
^{(0)} = -*z⃗*(*x⃗*
^{(0)})

Set residual *p⃗*
^{(0)} = *d⃗*
^{(0)}

Define termination criterion *ε*

Set iteration counter *n* = 0

repeat

Find *γ*
^{(n)} that minimises Ψ(*x⃗*
^{(n)} + *γ*
^{(n)}
*d⃗*
^{(n)})

${\beta}^{\left(n+1\right)}=max(\frac{{\overrightarrow{p}}^{{\left(n+1\right)}^{T}}\left({\overrightarrow{p}}^{(n+1)}-{\overrightarrow{p}}^{\left(n\right)}\right)}{{\overrightarrow{p}}^{{\left(n\right)}^{T}}{\overrightarrow{p}}^{\left(n\right)}},0)$

*d⃗*
^{(n+1)} = *p⃗*
^{(n+1)} + *β*
^{(n+1)}
*d⃗*
^{(n)}

*n* = *n* + 1

**until** ∥*z⃗*(*x⃗*
^{(n)})∥ < *ε*

Different choices for *β* exist. Above we used the Polak-Ribière method. Note that this requires restarting of the conjugate gradient method whenever *β* < 0 to guarantee convergence. A variety of line search methods can be used, either utilising gradient information, such as the secant method, or using only function evaluations such as the quadratic fit method. Often an exact line search is too computationally expensive due to the large number of function or derivative computations. Inexact line search methods such as Armijo’s Rule define the bounds for acceptable step lengths which guarantee convergence. There the line search is terminated when the step length is within the valid range.

## 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], the time-domain transport model using a Finite Difference Scheme [27], and the steady-state DC case using a Finite Element scheme [28]. A summary of these approaches and their relation to other reconstruction methods is given in [29]. 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

If we define the *Adjoint problem*

where ${\overrightarrow{Q}}_{i}^{+}$(*ω*) is the *adjoint source* [20], then we have that PMD⃗${\mathrm{F}}_{(j,i)}^{\mathrm{T}}$ is a vector in two parts, whose *k*^{th}
components are given by ${\sigma}_{j\mathit{,}i}^{-1}$ ${\overrightarrow{\Phi}}_{i}^{+}$ (*ω*) × Φ⃗_{j} (*ω*) and ${\sigma}_{j\mathit{,}i}^{-1}$ ∇Φ⃗^{j} (*ω*) × ∇Φ⃗_{j} (*ω*) respectively. Here the notation *a⃗* × *b⃗* as defined in [20] 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

But ${\overrightarrow{\Phi}}_{m}^{+}$ (*ω*) is the solution to Eq. 29, so we form a vector

and solve

leading to:

In the following we extend the notation of §**2**. by adding a superscript to indicate the measurement type. Thus ${P}_{j}^{\mathcal{{\rm M}}}$ is the projection operator for measurement ${F}_{j\mathit{,}i}^{\mathcal{{\rm M}}}$, and ${b}_{j\mathit{,}i}^{\mathcal{{\rm M}}}$, ${\sigma}_{j\mathit{,}i}^{\mathcal{{\rm M}}}$ are the forward data, residual data, and standard deviation of the error for the *i*^{th}
measurement from source *j* of type 𝚳. To find the gradient for the normalised measurement types, we use the result from[20] that the Jacobian for $\overline{\mathbf{\Mu}}$
can be constructed from

from which we can write

which leads to the gradient :

We therefore need to construct two adjoint sources :

From here we proceed as before, finding two solutions *η⃗*_{j}
^{$\overline{\mathbf{\Mu}}$+} (0), *η⃗*_{j}
^{$\overline{\mathbf{\Mu}}$+} (1) and computing the gradient as :

This method is applicable for measurement types (18) and (20). The *n*^{th}
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].

#### 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][31], assuming that 10^{4} 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].

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.

The L_{2} norms of the data errors, **b**⃗^{2}, 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 3 shows the L_{2} 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}
.

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

Figure 4 shows the target images and CG reconstructions after 100 iterations. The high contrast objects were generally well localised, in particular for *μ*_{s}
, while on the other hand the central *μ*_{a}
object is very diffuse, and the small *μ*_{a}
object at the lower right was not recovered. There is also some cross-talk from the *μ*_{s}
object at the upper left into the *μ*_{a}
image. The fine detail of the background tissues was not recovered. This example demonstrates the difficulties one might expect in the reconstruction of clinical data obtained from complex tissue structures, and it shows that simplistic models of blobs in a homogeneous background, as presented in the previous example, may not be very relevant in the assessment of reconstruction performance of clinical images.

## 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

^{1} | In biomedical optics the absorption parameter is usually denoted μ
_{a}, with the scattering parameter denoted μ
_{s}′, and the diffusion parameter given by $\frac{1}{3\left({\mu}_{a}+{\mu \text{'}}_{s}\right)}$. 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 Diagnosis*C. 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]