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=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 =κ ₀, κ_ij =0, i≠j, 1≤i, j≤n. Commonly, the isotropic case is represented simply using a scalar diffusion coefficient κ ₀.

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 γ ₂=1/π, γ ₃=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 lnA=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 Ω. 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.

 

Fig. 1. Schematic representation of the mapping Φ : Ω'→Ω and its inverse Φ^-1 : Ω→Ω′.

Download Full Size | PDF

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 ₁,x ₂) to denote the Cartesian coordinates and the indices j and k to refer to the Cartesian coordinates. Now let us denote (y ₁,y ₂)=(r,θ), (y′₁,y′₂)=(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.

 

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.05 cm^-1 and background absorption µ _a,bg=0.1 cm^-1. In the inclusions, µ _a=0.2 cm^-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.

Download Full Size | PDF

We may now use Eq. (11) to evaluateκ′. Note that κ′ is evaluated in Ω′, that is, in the (x′₁,x′₂)-coordinates. However, κ and the matrix $\frac{\partial y_i}{\partial x_k}$ are obtained in (x ₁,x ₂)-coordinates, so to evaluate Eq. (11) we need to expand κ(x ₁,x ₂)=κ(x ₁(x′₁,x′₂),x ₂(x′₁,x′₂)) 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.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.

 

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=[ᾱ₁ ᾱ₂] contains the eigenvectors ᾱ _i , i=1, 2, and Λ=diag(λ ₁,λ ₂) the eigenvalues. The axis corresponding to ᾱ₁ and the larger eigenvalue λ ₁ is given a constant length, and the second axis is scaled by λ ₂/λ ₁. The colour reflects the value of the larger eigenvalue λ ₁.

Download Full Size | PDF

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 (lnA,φ), 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 lnA 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 κ ₀ are discretized to have a constant values in each element of the FEM mesh. Thus, in the LM iteration, we solve for x=[µ _a;κ ₀], where µ _a and κ ₀ 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 κ ₀ were the spatial constants µ _a=0.05 cm^-1 and κ ₀=0.05 cm^-1=κ _0,bg. In practise, the minimisation is performed in two stages, where first the scalar background values for µ _a and κ ₀ 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 κ ₀ was given its true value to start with, the estimate of κ ₀ contains large artefacts, which compensate for the differences between the model and reality. Indeed, if κ ₀ is not included in the reconstruction but fixed to its true value, the estimate of µ _a is much poorer. Since the estimate of κ ₀ is of no interest here, the estimate is not shown.

 

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

Download Full Size | PDF

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 κ=Udiag(λ ₁,λ ₂)U ^T, whereU contains information on the principal direction of the anisotropy θ;U=[cosθ-sinθ; sinθ cosθ], and λ ₁ and λ ₂ on the strength of the diffusion. Hence, we have z=[λ ₁;λ ₂;θ]∊ℝ^3L×1, where using the same discretization as above, λ ₁, λ ₂ 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…3L, where ${\sigma }_z=[{\sigma }_{{\lambda }_1};{\sigma }_{{\lambda }_2};{\sigma }_{\theta }]\in {ℝ}^{3L\times 1}$. In order to choose the numerical values for z∗=[λ ₁∗;λ ₂∗;θ∗] 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.

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

View Table

 

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 θ∗=0 in Ω′, and (d) the estimate in (c) mapped into Ω.

Download Full Size | PDF

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.