Compensation for geometric mismodelling by anisotropies in optical tomography

We propose an approach for the estimation of the optical absorption coefficient in medical optical tomography in the presence of geometric mismodelling. We focus on cases in which the boundaries of the measurement domain or the optode positions are not accurately known. In general, geometric distortion of the domain produces anisotropic changes for the material parameters in the model. Hence, geometric mismodelling in an isotropic case may correspond to an anisotropic model. We seek to approximate the errors due to geometric mismodelling as extraneous additive noise and to pose a simple anisotropic model for the diffusion coefficient. We show that while geometric mismodelling may deteriorate the estimates of the absorption coefficient significantly, using the proposed model enables the recovery of the main features. © 2005 Optical Society of America OCIS codes: (170.3010) Image reconstruction techniques, (170.3880) Medical and biological imaging, (170.6960) Tomography References and links 1. S. R. Arridge, ”Optical tomography in medical imaging,” Inv. Probl. 15, R41-R93 (1999). 2. I. Nissilä, T. Noponen, J. Heino, T. Kajava, and T. Katila ”Diffuse optical imaging,” in Advances in Electromagnetic Fields in Living Systems 4, J. Lin, ed. (Kluwer, New York, to be published) 3. H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, ”Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results,” Appl. Opt. 42, 135-145 (2003). 4. A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, ”Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42, 5181-5190 (2003). 5. B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Welss, K. S. Osterman, U. L. Ostenberg, and K. D. Paulsen, ”Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: Pilot results in the breast,” Radiology 218, 261-266 (2001). 6. M. A. Franceschini, K. T. Moesta, S. Fantini, G. Gaida, E. Gratton, H. Jess, W. W. Mantulin, M. Seeber, P. M. Schlag, and M. Kaschke, ”Frequency-domain techniques enhance optical mammography: initial clinical results,” Proc. Natl. Acad. Sci. USA 94, 6468-6473 (1997). (C) 2005 OSA 10 January 2005 / Vol. 13, No. 1 / OPTICS EXPRESS 296 #5530 $15.00 US Received 20 October 2004; revised 28 December 2004; accepted 2 January 2005 7. B. Chance, E. Anday, S. Nioka, S. Zhou, L. Hong, K. Worden, C. Li, T. Murray, Y. Ovetsky, D. Pidikiti, and R. Thomas, ”A novel method for fast imaging of brain function, non-invasively, with light,” Opt. Exp. 2, 411-423 (1998). 8. M. A. Franceschini, V. Toronov, M. E. Filiaci, E. Gratton, and S. Fantini, ”On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Exp. 6, 49-57 (2000). 9. A. Y. Bluestone, G. Abdoulaev, C. H. Schmitz, R. Barbour, and A. H. Hielscher, ”Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Exp. 9, 272-286 (2001). 10. J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, ”Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol.47, 4155-4166 (2002). 11. J. C. Hebden, A. Gibson, T. Austin, R. Yusof, N. Everdell, D. T. Delpy, S. R. Arridge, J. H. Meek, and J. S. Wyatt, ”Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys. Med. Biol. 49, 1117-1130 (2004). 12. J. J. Stott, J. P. Culver, S. R. Arridge, and D. A. Boas, ”Optode positional calibration in diffuse optical tomography,” Appl. Opt.42, 3154-3162 (2003). 13. D. A. Boas, T. Gaudette, and S. R. Arridge, ”Simultaneous imaging and optode calibration with diffuse optical tomography,” Opt. Exp. 8, 263-270 (2001). 14. T. Tarvainen, V. Kolehmainen, M. Vauhkonen, A. Vanne, J. P. Kaipio, A. P. Gibson, S. R. Arridge, and M. Schweiger, ”Computational calibration method for optical tomography,” Applied Optics, 2005, in press. 15. J. P. Kaipio and E. Somersalo, Computational and Statistical Methods for Inverse Problems (Applied Mathematical Sciences 160, Springer-Verlag, New York, ISBN 0-387-22073-9, 2004). 16. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, ”A finite element approach for modeling photon transport in tissue,” Med. Phys. 20, 299-309 (1993). 17. M. Schweiger, S. R. Arridge, M. Hiraoka, 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). 18. J. Heino and E. Somersalo, ”Estimation of optical absorption in anisotropic background,” Inv. Probl. 18, 559-73 (2002). 19. J. Sylvester, ”An anisotropic inverse boundary value problem,” Comm. Pure Appl. Math. 43, 201-232 (1990). 20. V. Kolehmainen, Department of Applied Physics, University of Kuopio, P.O.B 1627, FIN-70211 Kuopio, Finland, M. Lassas, and P. Ola are preparing a manuscript to be called ”Inverse conductivity problem with an imperfectly known boundary.” 21. S. Järvenp ̈ aä, Implementation of PML Absorbing Boundary Condition for Solving Maxwell’s Equations with Whitney Elements (PhD Thesis, University of Helsinki, 2001). 22. J. Heino and E. Somersalo, ”A modelling error approach for the estimation of optical absorption in the presence of anisotropies,” Phys. Med. Biol. 49, 4785-4798 (2004). 23. Åke Björck,Numerical Methods for Least Squares Problems (SIAM, Philadelphia, 1996). 24. J. E. Dennis and R. B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations (SIAM, Philadelphia, 1996).


Introduction
Optical tomography [1,2] is a relatively recent method for non-invasive medical imaging with potential applications, for example, in the detection of breast cancer [3][4][5][6], imaging of the hemodynamic changes in the brain [7][8][9], and detection and localisation of problems in blood perfusion and oxygenation in the head of premature infants [10,11].In optical tomography, generally the purpose is to reconstruct the spatial distribution of the optical properties inside the body based on the properties of the measured light on the surface of the body.The measurements are done by guiding light onto the boundary of the tissue, commonly using optical fibers, and observing the part of the light traversed through the tissue at the measurement locations around the body.
A common approximation when solving the inverse problem in optical tomography is to assume that the body is isotropic and that the boundary shape and the optode locations are known.However, in reality the boundary shape and the optode locations are only known up to certain accuracy.The discrepancy between the model and the reality is likely to cause artefacts and degrade the quality of the obtained estimates.The situation is especially difficult when difference imaging, that is, estimating the changes in the optical properties between two measurement sets with fixed geometry, is not suitable.In order to reduce the effects of the discrepancies in modelling, the inaccuracies can be taken into consideration with different approaches.One approach to is apply calibration during the reconstruction.In [12], an algorithm for calibration for the optode locations was presented, and calibration algorithms for the optode gains (due to various effects such as laser power, fiber-tissue coupling and superficial skin properties) have been discussed in [13,14].Another approach is to build a statistical model for the discrepancies and compute how these are shown in the measurements.It is then possible to redefine the measurement error model in a way that it takes these model induced errors into account [15].
In this paper, we propose an approach which is related to the general framework that was treated in [15].We start by considering a mapping that explains certain types of discrepancies between the model and the reality, and subsequently modify the original observation model.The paper is organized as follows.In section 2, we investigate what kind of effects the geometric discrepancies induce on the usual computational model that is used for inversion.It is known that geometric discrepancies may transform an isotropic medium into an anisotropic one.We illustrate this effect with a numerical example.In section 3, we consider the inverse problem of the estimation of the absorption coefficient with geometric mismodelling and explain how the measurement model can be modified to take the geometric mismodelling into account.In section 4 we conclude the paper by discussing the extensions of the proposed methods to the more general Bayesian framework of [15].

The basic equations
Let Ω ∈ R n , n = 2orn = 3, be a bounded domain with a connected boundary.In this work, we model the light propagation in strongly scattering media by the diffusion equation, which is a commonly used model for light propagation in optical tomography [1,[16][17][18].For an intensity modulated source the diffusion equation can be written as where u is the photon density, κ ∈ R n×n is the diffusion coefficient, µ a the absorption coef- ficient and k = ω/c, where ω is the angular modulation frequency of the source term q and c is the speed of light.In the general form, κ ∈ R n×n presents anisotropic light propagation.By anisotropic, we mean that light propagation depends on the absolute direction in tissue, as opposed to the isotropic case, where light propagation does not depend on the absolute direction.The isotropic case is obtained by setting in Cartesian coordinates κ(i, i)=κ ii = κ 0 , κ ij = 0, i = j, 1 ≤ i, j ≤ n.Commonly, the isotropic case is represented simply using a scalar diffusion coefficient κ 0 .
On the boundary of the body ∂ Ω, we assume the Robin boundary condition, Here, the source is represented either as an equivalent interior source ( f = 0, q = 0) or as a diffuse boundary source (q = 0, f = 0).The constant γ = γ n depends on the dimension n and has values γ 2 = 1/π, γ 3 = 1/4.The inverse problem is to determine the optical parameters κ and µ a , or at least one of them, from the measurements on the boundary of the body.The boundary data consists of measured outward flux u out = −n•κ∇u at optode locations x j ∈ ∂ Ω, 1 ≤ j ≤ N on the boundary.Using the Robin boundary condition, the outward flux is obtained simply as u out (x j )=2γ u(x j ), 1 ≤ j ≤ N. Generally, one or more measurement types may be derived from the data for the estimation of the optical parameters.In this work, we use the natural logarithm of the amplitude ln A = ln|u out | and the phase angle ϕ = argu out of the complex boundary flux.
We study here the following problem: Assume that either the shape of Ω or the sensor/source locations are not well known.However, we have a model Ω with sensors and sources at fixed positions, and we use this model for inversion.The question is, what kind of errors the possible inaccuracies of the model bring and how should they be taken into account when calculating the estimates of the optical parameters.
Assume that Ω and Ω are diffeomorphic, that is, there is a one-to-one mapping Φ : Ω → Ω such that Φ as well as its inverse Φ −1 : Ω → Ω are differentiable.Furthermore, we assume that Φ maps the boundary ∂ Ω to ∂ Ω.Let us denote Furthermore, let u = u(x) denote the physical solution of the diffusion equation in Ω.B y physical solution, we mean the actual solution in the real word, i.e., e.g., in the human body or in a phantom.Now, we may write Hence, in the model domain Ω , there exists a corresponding function that we denote by u and that operates on the model coordinates x ∈ Ω .A schematic representation of the mappings is depicted in Fig. 1.Thus, u (x )=u(x).Then u satisfies also a diffusion equation.To see this, let us write the diffusion equation for u in the weak form.If ψ is a smooth test function, multiplying both sides of the Eq. ( 1) by ψ and integrating by parts, we arrive at the equation where the notation µ = µ a − ik was introduced.Taking into account the boundary condition (Eq.( 2)) we have further 2γ We make a change of variables in this integral.Let us denote A straightforward application of the chain rule shows that and in the same relation applies between ψ and ψ , where ψ (x )=ψ(x).Here, D denotes the differential, i.e., DΦ or DΦ −1 is a matrix, whose element ( j, k) is given by DΦ Therefore, Eq. ( 6) is transformed into where κ , µ , q , f and γ are given by Above, φ −1 in Eqs. ( 14) and ( 15) denotes the restriction of Φ −1 on the boundary ∂ Ω. Comparing the weak forms (Eqs.( 6) and ( 10)) we deduce that transforming the original domain Ω by a diffeomorphism Φ −1 , Eq. ( 1) along with the boundary condition (Eq.( 2)) remain valid if we modify the diffusion tensor, the absorption coefficient, the source and the boundary condition according to the Eqs.( 11)-( 15), respectively.There are two important issues worth noticing here.First, even if κ is a scalar, κ in general is not.In other words, modelling an isotropic body in a different geometry induces anisotropy into the model.Second, µ has an imaginary part that may not be constant.Physically, this corresponds to non-constant propagation speed.

The natural question thus arises: When solving the inverse problem based on a model that differs geometrically from the true body, should we allow anisotropies and non-constant propagation speeds even when the true object is isotropic and having constant propagation speed?
In this article, we study this question in the light of numerical simulations.The emphasis here is on anisotropies and we shall not discuss non-constant propagation speeds.
It is also worth noting, that given any two domains Ω and Ω , the mapping Φ is generally not uniquely defined.Hence, choosing any Φ : Ω → Ω that fullfills our requirements is merely one possibility, as well as the transformed quantities defined in the Eqs.( 11)-( 15) represent merely one possible set of quantities that give us the same solution in Ω as the original quantities in Ω.We remind that in general, inverse boundary value problems for anisotropic media suffer from non-uniqueness, one source of non-uniqueness being diffeomorphic transformations that leave the boundary intact (see [19]).In fact, the cited article shows that for the electrical impedance tomography (EIT) problem in two dimensions, this is the only source of non-uniqueness.In general, the characterization of the non-uniqueness of anisotropic problems is an open problem.In [20], the authors consider the related anisotropic EIT problem where they look for the minimally anisotropic solution that is uniquely defined.An extension of this approach to optical tomography has not been considered yet.

Numerical example
Let us illustrate the generation of anisotropies by a numerical example.The numerical calculations in this work are carried out in the two-dimensional (2D) space.Assume that the real domain Ω has the resemblance of a circle of radius 4 cm.However, assume that we do not know the boundary exactly and use a 4 cm radius circle Ω as the domain in the computational model.Let Ω and Ω be related by the mappings Φ : Ω → Ω and Φ −1 : Ω → Ω .
For this example, let the actual (unknown) deformation in spherical coordinates be Φ : (r , θ ) → (r, θ ): which defines the real domain of this numerical example.Figure 2 illustrates both the employed model and the actual domains and the corresponding optode positions.
Let the diffusion coefficient κ in the real domain Ω be isotropic.To compute the corresponding anisotropic κ (Eq. ( 11) in the modelled domain, we need to calculate the Jacobian DΦ −1 : where we use (x 1 , x 2 ) to denote the Cartesian coordinates and the indices j and k to refer to the Cartesian coordinates.Now let us denote (y 1 , y 2 )=(r, θ ), (y 1 , y 2 )=(r , θ ), and use the indices x 2 [cm] Fig. 2. Transformation defined in Eq. ( 16).The actual domain Ω.The data used in the numerical inversion example is generated using this model with constant isotropic κ = κ 0 = 0.05 cm −1 and background absorption µ a,bg = 0.1cm −1 .In the inclusions, µ a = 0.2cm −1 .The simulated data set was collected by placing the source at each optode location at turn, and using the rest of the optodes for measurement, resulting in 240 measurements of amplitude and phase each.
( i and in the polar coordinates.We then write The matrices corresponding to are readily obtained using the mappings from spherical to Cartesian and Cartesian to spherical coordinates.The matrix is obtained from Eq. ( 16) using D( Φ−1 )=(D Φ) −1 .Hence, we do not need to calculate the inverse mapping Φ−1 explicitly.
We may now use Eq. ( 11) to evaluate κ .Note that κ is evaluated in Ω , that is, in the (x 1 , x 2 )- coordinates.However, κ and the matrix ∂ y i ∂ x k are obtained in (x 1 , x 2 )-coordinates, so to evaluate Eq. ( 11) we need to expand κ( In practise the numerical simulations in this work are conducted using the finite element method (FEM).In the FEM framework, we define the diffusion tensor in a piecewise constant basis, that is, κ | e j = κ j = const, where Ω = ∪ L j=1 e j is divided into triangular elements e j .The evaluation of κ is done numerically in each element e j using the center of the element for the coordinate values.In Fig. 3, the diffusion tensor κ generated by Eq. ( 16) is illustrated in some elements.The real diffusion tensor κ was a constant scalar κ = κ 0 = 0.05 cm −1 .It is clear that the anisotropy generated by the associated geometric mismodelling is severe.Thus, using an isotropic model with the modelled geometry Ω as such can lead to significant artefacts to the estimated scalar absorption coefficient.16).The tensor distribution is illustrated by ellipses.The axes of an ellipse correspond to the directions of the diffusion tensor eigenvectors.The diffusion tensor is written as κ = UΛU T , where U =[ᾱ 1 ᾱ2 ] contains the eigenvectors ᾱi , i = 1, 2, and Λ = diag(λ 1 , λ 2 ) the eigenvalues.The axis corresponding to ᾱ1 and the larger eigenvalue λ 1 is given a constant length, and the second axis is scaled by λ 2 /λ 1 .The colour reflects the value of the larger eigenvalue λ 1 .

Estimation of the absorption coefficient
Consider the inverse problem in the presence of discrepancies in the boundary shape and optode positions between the model and the reality.More precisely, we assume that the diffusion and absorption coefficients in the domain of interest Ω are unknown.The main interest is in the spatial distribution of the optical absorption coefficient.The boundary data consists of frequency domain data (ln A, ϕ), where A and ϕ are the measured amplitudes and phase shifts, respectively, and the model with circular domain Ω is used in the inversion.
The observation model used in this work is where y is a boundary data vector consisting of measured ln A and ϕ, G(µ a , κ) is the FEM model of the noiseless observations and ν is additive measurement noise.The data is created using the isotropic diffusion model in the geometry that is depicted in Fig. 2(b).The sources are approximated by virtual interior point sources under the surface so that we have f = 0 and q(x)=qδ (x s ), where x s is the virtual source position and δ is the Dirac delta.For data generation, we use a mesh consisting of 7014 elements and for the reconstructions we use a different mesh of 2993 elements.Both meshes we created using a bubble mesh generator [21].
The noise level was assumed to be ∼ 1 % for the amplitude, and ∼ 1 • (constant level) for the phase.Artificial noise corresponding to these levels was added to the computed signal.
In the following we consider two approaches for the estimation of optical absorption.The first approach is the conventional approach that is based on isotropic diffusion model and assuming that the employed geometry is correct.In the second approach we first pose a very crude statistical model for the (background) anisotropy with specified uncertainty for the associated parameters.This model is then used to generate the model for the extraneous measurement errors.

Conventional isotropic reconstruction
The traditional reconstruction is performed here using an isotropic model and truncated Levenberg-Marquardt (LM) iteration.The (truncated) LM iteration is described in more detail in the Appendix under Optimization.In the following, we give a few points related to the implementation in this work.The optical parameters µ a and κ 0 are discretized to have a constant values in each element of the FEM mesh.Thus, in the LM iteration, we solve for x =[µ a ; κ 0 ], where µ a and κ 0 are vectors of length L, the number of elements in the FEM mesh.We assume that we only have independent Gaussian additive measurement noise, so the noise covariance is a diagonal matrix with the standard deviations of the logarithm of the amplitude and the phase given by σ lnA = 0.01 and σ ϕ = 0.0175, respectively.
The initial values for the iteration of µ a and κ 0 were the spatial constants µ a = 0.05 cm −1 and κ 0 = 0.05 cm −1 = κ 0,bg .In practise, the minimisation is performed in two stages, where first the scalar background values for µ a and κ 0 were searched for, after which the iteration was carried out in pixel basis.
Figure 4(a) displays the reconstruction of the absorption coefficient.The positions of the perturbations are clearly shifted, and the lower perturbation is not very clear.Figure 4(b) displays the same estimate mapped into the actual domain.Although in the reconstruction κ 0 was given its true value to start with, the estimate of κ 0 contains large artefacts, which compensate for the differences between the model and reality.Indeed, if κ 0 is not included in the reconstruction but fixed to its true value, the estimate of µ a is much poorer.Since the estimate of κ 0 is of no interest here, the estimate is not shown.

Reconstruction using modelling error approach
Inspired by the fact that the differences between the model and reality may be interpreted as the generation of anisotropies, we attempt the reconstruction of the absorption coefficient by including anisotropies into the estimation scheme.It is clear, however, that our knowledge on the potential anisotropies generated is very limited and possibly of qualitative nature rather than quantitative.Therefore, using any fixed guess for the structure of the anisotropy in the estimation is likely to fail.To overcome this problem, we use a method in which the uncertainty in the anisotropy is first modelled, a modified measurement error model is then constructed, and finally the truncated LM iteration employed.The computational scheme is described in more detail in the Appendix under Model uncertainties.Other examples and more details may also be found in [22].For the general approach dealing with the construction of the so-called enhanced measurement error models, see [15].With the notations used in the Appendix, we have the quantity of primary interest x = µ a with the nuisance parameter z representing anisotropy.To implement the method, we write the diffusion tensor in the form κ = Udiag(λ 1 , λ 2 )U T , where U contains information on the princi- pal direction of the anisotropy θ ; U =[cos θ − sin θ ; sin θ cos θ ], and λ 1 and λ 2 on the strength of the diffusion.Hence, we have z =[λ 1 ; λ 2 ; θ ] ∈ R 3L×1 , where using the same discretization as above, λ 1 , λ 2 and θ are vectors containing the values of the respective parameters in each element of the FEM mesh.All parameters are assumed to be mutually independent, and the covariance matrix Γ z ∈ R 3L×3L is assumed to be of the form Γ z = diag(σ 2 z (i)), i = 1 ...3L, where σ z =[σ λ 1 ; σ λ 2 ; σ θ ] ∈ R 3L×1 .In order to choose the numerical values for z * =[λ 1 * ; λ 2 * ; θ * ] and σ z , for simplicity, we first assume that the values are spatially constant.Then, for each para- meter, we consider a (large) plausible interval in which the parameter values may vary, and set the midpoints and variances so that a corresponding Gaussian distribution covers well the interval.The numerical values are listed in Table 1.For the measurement noise, we have again σ lnA = 0.01 and σ ϕ = 0.0175.
For the angle θ , it might be possible to make an intelligent guess based on the geometry.Light propagates better in the preferential direction of the anisotropy, which can be pictured to correspond to a smaller distance.Hence, e.g., in the example of Fig. 2, one might deduce that choosing θ * = π/2 (vertical direction) is likely to give better results.Here, we have included two cases of estimation.In the first one, we have used θ * = π/2, corresponding to an intelligent guess, and in the second one, θ * = 0, which can be taken as a "worst case" scenario.The initial value of µ a used in both reconstructions was µ a = 0.05 cm −1 .In practise, the minimisation is again performed in two-stages, where firstly, the scalar background value for µ a is searched for, and secondly, the iteration is changed into pixel basis.The estimates of the absorption coefficient with the modelling error approach are depicted in Fig. 5.

Conclusions
In this paper, we have proposed an approach for the estimation of optical absorption coefficient when information on the boundaries of the measurement domain and the optode positions is inaccurate.By posing a simple model for the background anisotropy and computing an approximation for the enhanced error model including the modelling error as noise, we are able to improve the quality of the estimated distribution of µ a .This paper should be taken as a proof (C) 2005 OSA 10 January 2005 / Vol. 13, No. 1 / OPTICS EXPRESS 305 of concept, as the modelling and the numerical experiments were conducted in 2D, whereas reconstructions with real data will generally require a full 3D model to be successful.Also, including anisotropies into the model presents one possibility to help the reconstruction when the measurement geometry is not accurately known.Calibration algorithms for optode positions present another possibility.In this work, we have used Gaussian priors for the anisotropy parameters for the sake of computational convenience.However, a uniform distribution for the angle θ could be more appropriate when the prior information available on the geometry does not suggest any preferential direction.Here we have used only an approximation of the second moment of the modelling error.Based on our understanding of a related problem of electrical impedance tomography studied in [15], it is possible that the mean of the modelling error plays a significant role.The computation of the mean would necessitate the use of stochastic simulation which is out of the scope of this paper.Equivalently, using non-Gaussian priors would also require stochastic simulation, for example by using Monte Carlo Markov Chain (MCMC) methods, which is computationally much heavier.
The extensions of the proposed approach employing more general models for the induced anisotropy all require stochastic simulation methods for the computation of the estimate of the absorption coefficient.However, we could then pose a realistic probabilistic model on the boundary and thus use a more complex model for the background anisotropy.Furthermode, this would allow for a more precise model for the enhanced error model.
where D z G(x c , z * ) is the differential, or Jacobian matrix, of the mapping G with respect to the variable z.Note that above we ignore the associated expectation of the modelling error E{ε(x, z)} whose computation would require the use of stochastic simulation methods such as Markov chain Monte Carlo methods, which are out of scope of this paper.If Γ z is the a priori covariance of the parameter z expressing the degree of uncertainty of the value z * , cov(z)=E δ zδ z T = Γ z , we find that the covariance of the modelling error is, within the linearized approximation, Assuming statistical independence between the noise n and the parameters x and z, we then get an enhanced additive error model where Γ ν is the additive noise covariance.Hence, we have reduced the problem to a standard additive noise model.For brevity, we write G(x, z * )=G(x) in the sequel.In general, the above formula can also be used as a starting point for any iterative estimation scheme with regularization.

Optimization
Consider first an observation model with additive error, where y is a vector that we measure, x is a vector representing the quantity of primary interest and e is the additive noise.If we assume that e is Gaussian, zero mean and has the covariance matrix Γ, the likelihood density of the measurement y with given x is The maximum likelihood estimator (ML) of x is the value that maximizes the likelihood of the observed outcome y.Evidently, such an estimator is found by minimizing the functional The Levenberg-Marquardt (LM) method is a classical iterative method for computing the ML estimator.Assume again that x c is the current estimate of x.We write a local linearization of G around (C) 2005 OSA 10 January 2005 / Vol. 13, No. 1 / OPTICS EXPRESS 297
(C) 2005 OSA 10 January 2005 / Vol. 13, No. 1 / OPTICS EXPRESS 300 Fig. 2. Transformation defined in Eq. (16).(a) The domain Ω used in the inversion.The 16 optodes are located at equal distances on the boundary.The figure also depicts two inclusions inside the domain.(b)The actual domain Ω.The data used in the numerical inversion example is generated using this model with constant isotropic κ = κ 0 = 0.05 cm −1 and background absorption µ a,bg = 0.1cm −1 .In the inclusions, µ a = 0.2cm −1 .The simulated data set was collected by placing the source at each optode location at turn, and using the rest of the optodes for measurement, resulting in 240 measurements of amplitude and phase each.

Fig. 3 .
Fig.3.Diffusion tensor generated by the transformation defined in Eq. (16).The tensor distribution is illustrated by ellipses.The axes of an ellipse correspond to the directions of the diffusion tensor eigenvectors.The diffusion tensor is written as κ = UΛU T , where U =[ᾱ 1 ᾱ2 ] contains the eigenvectors ᾱi , i = 1, 2, and Λ = diag(λ 1 , λ 2 ) the eigenvalues.The axis corresponding to ᾱ1 and the larger eigenvalue λ 1 is given a constant length, and the second axis is scaled by λ 2 /λ 1 .The colour reflects the value of the larger eigenvalue λ 1 .

1 ]Fig. 4 .
Fig. 4. (a) Estimate of the absorption coefficient using the isotropic LM iteration in domain Ω used for inversion, and (b) the estimate in (a) mapped onto the actual domain Ω.

1 ]Fig. 5 .
Fig. 5. Estimates for the optical absorption using the modelling error approach combined with LM iteration.(a) Case 1 with θ * = π/2 in the domain Ω , and (b) the estimate in (a) mapped into the real domain Ω.(c) Case 2 with θ * = 0i nΩ , and (d) the estimate in (c) mapped into Ω.

2 Γ − 1 subject
x c , G(x)=G(x c + δ x) ≈ G(x c )+D G(x c )δ x. (29)Assuming that the local linearization is reliable within a ball of radius r > 0 around the current estimate, the Levenberg-Marquardt step comprises solving the approximate constrained problem,Minimize y − (G(x c )+D G(x c )δ x) to δ x ≤r.(30)This is a least squares inequality constrained (LSQI) problem, and by using the Lagrange multipliers, it leads to the problem of minimizing the functionaly − (G(x c )+D G(x c )δ x) 2 Γ −1 + λ δ x 2 .(31) (C) 2005 OSA 10 January 2005 / Vol. 13, No. 1 / OPTICS EXPRESS 307

Table 1 .
The parameter values used in the modelling error approach.