## Abstract

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

## 1. 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–6], imaging of the hemodynamic changes in the brain [7–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 de-grade 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].

## 2. Generation of anisotropies

#### 2.1. The basic equations

Let Ω∊ℝ
^{n}
, *n*=2 or *n*=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–18]. For an intensity modulated source the diffusion equation can be written as

where *u* is the photon density, *κ*∊ℝ^{n×n} is the diffusion coefficient, *µ*
_{a} the absorption coefficient 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, *κ*∊ℝ^{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 φ=arg*u*
_{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 Ω. By 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

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Φ
_{jk}
=*∂*Φ
_{j}
/*∂x*′
_{k}
=*∂x*_{j}
/*∂x*′
_{k}
or D${\mathrm{\Phi}}_{\mathit{\text{jk}}}^{-1}$
=*∂*${\mathrm{\Phi}}_{j}^{-1}$
/*∂x*_{k}
=*∂x*′
_{j}
/*∂x*_{k}
, respectively. Furthermore, we have

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.

#### 2.2. 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 $\tilde{\Phi}$
: (*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

*i* and *ℓ* in the polar coordinates. We then write

The matrices corresponding to $\frac{\partial {x\prime}_{j}}{\partial {y\prime}_{\ell}}$ and $\frac{\partial {y}_{i}}{\partial {x}_{k}}$ are readily obtained using the mappings from spherical to Cartesian and Cartesian to spherical coordinates. The matrix $\frac{\partial {y\prime}_{\ell}}{\partial {y}_{i}}=D\left({\tilde{\Phi}}^{-1}\right)$ is obtained from Eq. (16) using D($\tilde{\Phi}$
^{-1})=(D$\tilde{\Phi}$
)^{-1}. Hence, we do not need to calculate the inverse mapping $\tilde{\Phi}$
^{-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 $\frac{\partial {y}_{i}}{\partial {x}_{k}}$ are obtained in (*x*
_{1},*x*
_{2})-coordinates, so to evaluate Eq. (11) we need to expand *κ*(*x*
_{1},*x*
_{2})=*κ*(*x*
_{1}(*x*′_{1},*x*′_{2}),*x*
_{2}(*x*′_{1},*x*′_{2})) and likewise for $\frac{\partial {y}_{i}}{\partial {x}_{k}}$.

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 Ω′=${\cup}_{j=1}^{L}$
*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.

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

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

#### 3.2. 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 *κ*=*U*diag(*λ*
_{1},*λ*
_{2})*U*
^{T}, where*U* contains information on the principal 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};*θ*]∊ℝ^{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}
∊ℝ^{3L×3L} is assumed to be of the form Γ
_{z}
=diag(${\mathrm{\sigma}}_{z}^{2}$
(*i*)), *i*=1…3*L*, where ${\sigma}_{z}=[{\sigma}_{{\lambda}_{1}};{\sigma}_{{\lambda}_{2}};{\sigma}_{\theta}]\in {\mathbb{R}}^{3L\times 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 parameter, 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.

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

## Appendix

In this section we summarize the computational aspects of the article, the model uncertainties and optimization methods employed in the examples.

*Model uncertainties*

Consider a model for a noisy measurement of a quantity with additive noise,

Here, *ν* is a vector representing the noise, the vector *x* is the quantity of primary interest and *z* is a vector of poorly known parameters of secondary interest. The mapping *G*, linear or non-linear, is assumed to be known.

Assume that, e.g., for saving computational efforts, we try to estimate *x* only from the observed values of *y*, while *z* is considered as a nuisance parameter. We would like to fix the parameter *z* to a value *z*∗ that, to our knowledge, represents a typical or plausible value, but, however, the uncertainty about this value is considerable. We write

The term *ε* (*x*, *z*) represents a modelling error that, in general, cannot be neglected. The difficulty is that it depends on unknowns *x* and *z*. If we model *x* and *z* as random variables as is customary in statistical inversion theory [15], we can in principle estimate its probability distribution. In general, this is a computationally challenging task.

A strongly simplified but useful approximation for the probability density is obtained by linearization and Gaussian approximation. Assume that *x*
_{c} is the current approximation for the value of *x*. We write

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*∗,

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 *x*
_{c},

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,

This is a least squares inequality constrained (LSQI) problem, and by using the Lagrange multipliers, it leads to the problem of minimizing the functional

The Lagrange parameter *λ*≥0 is solved from the inequality constraint, see [23]. The LM updated value is then *x*
_{+}=*x*
_{c}+*δx*, where *δx* is the minimimizer of the above quadratic functional,

If *x*
_{0} is the initial value of the iteration, we obtain a sequence of estimates, *x*
_{1},x2,…. There are automatic ways of selecting *λ* when the trust region parameter *r* is given, see [24]. The ML estimator, from the point of view of ill-posed inverse problems, is useless since small errors in the observations produce huge oscillations in the ML estimator. There are several ways to overcome these difficulties. The classical methods include various regularization techniques. One of these methods is the truncation of the iteration after some step. This is one of the most popular approaches in optical tomography and what we also adopt here.

When using well-defined probabilistic model for the associated uncertainties, it is possible to turn to statistical inversion methods that are based on the Bayesian framework of statistics. The Bayesian approach is to introduce a prior probability density, *π*
_{pr}(*x*) and consider the posterior density that is, according to the Bayes’ formula,

$$\propto \mathrm{exp}(-\frac{1}{2}{\parallel y-G\left(x\right)\parallel}_{{\Gamma}^{-1}}^{2}-V\left(x\right)),$$

where

The *maximum a posteriori estimator* (MAP) is then the minimizer of the expression

In the language of classical approach, the term *V* can be seen as a regularizing penalty term. If properly chosen, it prevents the MAP estimate from being corrupted by the measurement noise and modelling errors.

## Acknowledgements

This work was supported by the Academy of Finland and the Emil Aaltonen foundation.

## References and links

**1. **S. R. Arridge, “Optical tomography in medical imaging,” Inv. Probl. **15**, R41–R93 (1999). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [PubMed]

**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). [CrossRef] [PubMed]

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

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

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

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

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

**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. [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**18. **J. Heino and E. Somersalo, “Estimation of optical absorption in anisotropic background,” Inv. Probl. **18**, 559–73 (2002). [CrossRef]

**19. **J. Sylvester, “An anisotropic inverse boundary value problem,” Comm. Pure Appl. Math. **43**, 201–232 (1990). [CrossRef]

**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ää, *Implementation of PML Absorbing Boundary Condition for Solving Maxwell’s Equations with Whitney Elements* (PhD Thesis, University of Helsinki, 2001). [PubMed]

**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) . [CrossRef] [PubMed]

**23. **Åke Björck, *Numerical Methods for Least Squares Problems* (SIAM, Philadelphia, 1996). [CrossRef]

**24. **J. E. Dennis and R. B. Schnabel, *Numerical Methods for Unconstrained Optimization and Nonlinear Equations* (SIAM, Philadelphia, 1996). [CrossRef]