OSA's Digital Library

Biomedical Optics Express

Biomedical Optics Express

  • Editor: Joseph A. Izatt
  • Vol. 1, Iss. 1 — Aug. 2, 2010
  • pp: 209–222
« Show journal navigation

Corrections to linear methods for diffuse optical tomography using approximation error modelling

Tanja Tarvainen, Ville Kolehmainen, Jari P. Kaipio, and Simon R. Arridge  »View Author Affiliations


Biomedical Optics Express, Vol. 1, Issue 1, pp. 209-222 (2010)
http://dx.doi.org/10.1364/BOE.1.000209


View Full Text Article

Acrobat PDF (947 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Linear reconstruction methods in diffuse optical tomography have been found to produce reasonable good images in cases in which the variation in optical properties within the medium is relatively small and a reference measurement with known background optical properties is available. In this paper we examine the correction of errors when using a first order Born approximation with an infinite space Green’s function model as the basis for linear reconstruction in diffuse optical tomography, when real data is generated on a finite domain with possibly unknown background optical properties. We consider the relationship between conventional reference measurement correction and approximation error modelling in reconstruction. It is shown that, using the approximation error modelling, linear reconstruction method can be used to produce good quality images also in situations in which the background optical properties are not known and a reference is not available.

© 2010 Optical Society of America

1. Introduction

In diffuse optical tomography (DOT), the goal is to produce quantitative images of the optical properties of an unknown medium based on near-infrared light transport measurements made on the surface of the object [1

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

,4

4. S.R. Arridge and J.C. Schotland, “Optical tomography: forward and inverse problems,” Inv. Probl. 25, 123010, (2009).

,7

7. A. Gibson, J.C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical tomography,” Phys. Med. Biol. 50, R1–R43, (2005).

]. Since the photons travel randomly inside the medium, image reconstruction is a highly ill-posed problem and needs to be approached within the framework of inverse problems. In these methods, light propagation is modelled using the radiative transport equation (RTE) or its approximations [1

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

,5

5. G. Bal, “Inverse transport theory and applications,” Inv. Probl. 25, 053001, (2009).

,13

13. A. D. Klose and A. H. Hielscher, “Optical tomography with the equation of radiative transfer,” International Journal of Numerical Methods for Heat & Fluid Flow 18, 443–464, (2008).

]. In turbid media, the most commonly applied model is the diffusion approximation (DA) to the RTE. The accurate modelling of light transport and usage of sophisticated image reconstruction methods have enabled production of good quality optical images which strongly supports the view that DOT is one of the most promising new medical imaging modalities. However, DOT still needs further development, and, in particular, accurate and fast image reconstruction methods that could be used simultaneously with measurements and in online monitoring need to be developed.

The existing approaches in DOT can be classified into those based on inverse scattering theory (IST) and analytical (Green’s function) solutions to the light propagation model [21

21. J. C. Schotland and V. Markel, “Inverse scattering with diffusing waves,” J. Opt. Soc. Am. A 18, 2767–2777, (2001).

], and those based on data fitting using optimisation techniques and numerical solutions to the underlying partial differential equations [1

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

], also known as the parameter identification method (PIM). Although these two methods are thought of as quite different, we will show in this paper that both can be put into a common framework when we consider a Bayesian approach to analysing the underlying model errors.

In both the IST and PIM frameworks, either a linear or a non-linear approach to the inverse problem can be developed, and the central computational requirement for the forward problem is to construct a suitable model for the sources and detectors which can be used to match the experimental data taken on the turbid medium’s surface. However, it is generally true that the difficulty of producing highly accurate computational models can be significant, and is counterintuitive in regard to the inherently poor resolution of DOT. In practice, the models used are only approximate and have to be ‘corrected’ (or to use a more principled term ‘calibrated’) to actually work in practice. Usually, the correction is achieved using reference measurements.

Linearisation of the inverse problem is explicit in the IST approach and in iterative Newton type algorithms in the PIM approach. In the former technique, the linearisation is built from analytic Green’s functions which solve the DA in predetermined domain such as an infinite space [2

2. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37, 1531–1560, (1992).

, 20

20. 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).

]. These linear methods are limited by the fact that they cannot predict large changes in the optical parameters and that they assume the background optical properties to be known [6

6. D. A. Boas, “A fundamental limitation of linearized algorithms for diffuse optical tomography,” Opt. Express 1, 404–413, (1997).

]. However, such methods have been found to produce good quality images in situations in which these limitations are fulfilled and a reference measurement from known background medium is available [15

15. S.D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16, 5048–5060, (2008).

]. The linearisation error in an another diffuse imaging modality, electrical impedance tomography, has been addressed in [23

23. N. Polydorides, “Linearization error in electical impedance tomography,” Progress In Electromagnetics Research, PIER 93, 323–337, (2009).

].

In this paper, we use the approximation error modelling theory to explain the connection between the IST and PIM approaches. The approximation error approach for the treatment of modelling errors was introduced in [11

11. J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York, 2005.

, 12

12. J. Kaipio and E. Somersalo, “Statistical inverse problems: Discretization, model reduction and inverse crimes,” J. Comput. Appl. Math. 198, 493–504, (2007).

]. The method has been applied in DOT in compensating modelling errors due to discretisation, domain mismodelling and in compensating modelling errors between the RTE and the DA [3

3. S.R. Arridge, J.P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inv. Probl. 22, 175–195, (2006).

, 8

8. J. Heino, E. Somersalo, and J.P. Kaipio, “Compensation for geometric mismodelling by anisotropies in optical tomography,” Opt. Express 13, 296–308, (2005).

, 14

14. V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S.R. Arridge, and J.P. Kaipio, “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A 26, 2257–2268, (2009).

, 24

24. T. Tarvainen, V. Kolehmainen, A. Pulkkinen, M. Vauhkonen, M. Schweiger, S.R. Arridge, and J.P. Kaipio, “Approximation error approach for compensating for modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography,” Inv. Probl. 26, 015005, (2010).

]. Outside optical imaging, the approximation error approach has been utilised in electrical impedance tomography [11

11. J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York, 2005.

, 16–19

16. A. Lehikoinen, S. Finsterle, A. Voutilainen, L.M. Heikkinen, M. Vauhkonen, and J.P. Kaipio, “Approximation errors and truncation of computational domains with application to geophysical tomography,” Inverse Problems and Imaging 1, 371–389, (2007).

] and it has been extended to time-dependent inverse problems [9

9. J.M.J. Huttunen and J.P. Kaipio, “Approximation errors in nonstationary inverse problems,” Inverse Problems and Imaging 1, 77–93, (2007).

, 10

10. J.M.J. Huttunen and J.P. Kaipio, “Model reduction in state identification problems with an application to determination of thermal parameters,” Applied Numerical Mathematics 59, 877–890, (2009).

]. In this paper, we develop an approach in which the approximation error approach is used to correct the linear reconstruction method in DOT. In this approach, a reference measurement is not needed and the background optical properties do not need to be known but a model for the uncertainty of the background coefficients is needed.

The rest of the paper is organised as follows. In Section 2, linear reconstruction methods in DOT are reviewed. In Section 3, the approximation error approach is described and its relation to reference measurement approach is explained. In Section 4, simulation results are shown and in Section 5 conclusions are given.

2. Deterministic model based reconstruction and correction methods

2.1. Definitions

We are considering the fitting of a model

y=A(x)+e,A:𝕏m
(1)

where y ∈ ℝm is discrete data, x ∈ 𝕏 is a function in some space 𝕏 representing the absorption and scattering images, and e denotes the noise which is usually modelled to be Gaussian distributed e ~ 𝒩(0,Γe) with zero mean and covariance Γe ∈ ℝm×m. In the ‘true’ physical problem we assume 𝕏 = L(Ω) but in a practical implementation the solution and forward mapping are represented in discrete vector spaces

xxhn,AAh
(2)

where Ah : ℝn → ℝm is a linear or non-linear operator.

Throughout this paper, within both the IST and PIM frameworks we assume that the light propagation model is the diffusion approximation

(·D+μa+iωc)Φ=0,inΩ
(3)
Φ+2ζDΦν=J,onΩ
(4)

which is denoted as

𝓛(x)Φ=J.
(5)

In Eqs. (3)–(5), D = (3(μa + μs))-1 is the diffusion coefficient, μs is the reduced scattering coefficient, μa is the absorption coefficient, Φ is the photon density, J- is prescribed input source function, i is the imaginary unit, c is the speed of light in the medium and ω is the angular modulation frequency of the input signal. Furthermore, in this paper we use the following notations

background optical parameters:x0 = (μa,0,μs,0)T
perturbed optical parameters:x=(μa,μs)T=x0 + δx
estimated optical parameters:x̂
measured data :gobs
reference data :gref
modelled data :y
modelled reference data :y0
discretized (incorrect) background optical parameters:xh,0
discretized (incorrect) optical parameters:xh = xh,0 + δxh
exact model for the forward operator in domain Ω:A
approximate model in (incorrect) domain Ω̃:Ah
augmented model :Ãh
Jacobian:J
modelling error:ε

2.2. Linear reconstruction model within inverse scattering theory

In the IST approach we represent the operator A in Eq. (1) as the measurement of the Green’s function

A(x)=𝓜𝓖(x)J
(6)

where 𝓖(x) is the Green’s operator, parameterised by x, that solves for the internal field Φ given the source function and the light propagation model (5) and ℳ is the measurement of Φ at the boundary of the domain Ω. In practice, measurements are needed for multiple input sources functions {J1-, J2-,…} where the functions {Ji-} form a basis for general functions on the boundary Ω, leading to a vector form for the forward model

A(x)=(𝓜1𝓖(x)J1𝓜2𝓖(x)J2𝓜S𝓖(x)JS).
(7)

Under this notation, given a reference state x0, a series solution for a perturbed state x0 + δx is given by

A(x0+δx)=𝓜[𝓘𝓖(x0)𝓥(δx)]1𝓖(x0)
=A(x0)+𝓜[𝓖(x0)𝓥(δx)+(𝓖(x0)𝓥(δx))2+]𝓖(x0),
(8)

where 𝓘 is an identity operator and 𝒱(δx) is a potential operator determined by the perturbation. Eq. (8) is a general Born series for diffuse photon density waves.

Linear inversion methods in DOT based on the series Eq. (8) consider only the first term in the series and make the approximation

A(x0+δx)=A(x0)+𝖩(x0)δx,
(9)

where the Jacobian 𝖩(x0) is the discrete representation of the Fréchet derivative of the nonlinear mapping A(x), and seek to invert the linear expression

gobsy0=𝖩(x0)δx,
(10)

where gobs are measurements and y0 = A(x0) is assumed to be an accurate prediction of the measurements that would correspond to the reference state x0.

It is natural to consider both the accuracy of the prediction of y0 at measurement position rm due to source located rs as well as the accuracy of representation of the linearisation 𝖩(x0) which both depend on the accuracy of the Green’s functions

y0|rm=𝓜ΩG(rm,rs;x0)J(rs)drs
(11)
(yy0)|rm=𝓜ΩG(rm,r;x0)𝓥(δx)G(r,rs;x0)J(rs)drs.
(12)

In many approaches, G is taken to be an analytic expression such as the infinite domain Green’s function for the diffusion equation

G(r1,r2)=eσr1r2r1r2
(13)

where σ=(μa+iω/cD)1/2 and |r1 - r2| is the distance between points r1 and r2. Expression (13) is the one employed in this paper, although Green’s functions in other geometries can be defined [2

2. S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37, 1531–1560, (1992).

], including general domains using the Kirchoff Approximation [22

22. J. Ripoll, V. Ntziachristos, and M. Nieto-Vesperinas, “The Kirchhoff approximation for diffusive waves,” Phys. Rev. E 64, 1–8, (2001).

].

2.3. Linear reconstruction model within the parameter identification method

In the parameter identification method, we adopt a non-linear data fitting procedure, such as the regularized output least square method

x̂=argminx{12gobsA(x)Γe2+Ψ(x)}
(14)

where x̂ denotes the estimated optical parameters and Ψ(x) is a regularizing penalty functional. Taking only the linear approximation Eq. (9) leads to

x̂=x0+argminδx{12gobsA(x0)𝖩(x0)δxΓe2+Ψ(δx)}.
(15)

2.4. Model correction using a reference measurement

Within both approaches, the model used is not necessarily accurate. Let A(x0) represent an ‘exact’ model for the forward operator in domain Ω and state x0 and let Ah(xh,0) represent an approximate model in the incorrect domain Ω̃ with incorrect optical parameters xh,0. Consider an additive correction ε0 so that the two models agree at the reference state x0

y˜h=Ah(xh)+e+ε0.
(16)

In the literature of optical tomography, it is common practice to use a measurement set gref from background x0 and then to use the augmented model

A˜h(xh)=Ah(xh)+grefAh(xh,0).
(17)

We can now thus define the Reference measurement correction

ε0=grefAh(xh,0)gref=Ah(xh,0)+ε0.
(18)

Using the reference measurement correction in Eq. (9) leads to the corrected linear inversion method

gobsgref=𝖩(xh,0)δxh
(19)

where δxh = (xh - xh,0). Thus, the minimisation is

x̂=argminxh{12gobsAh(xh)gref+Ah(xh,0)Γe2+Ψ(xh)}
(20)

which in the case of the linear approximation can be written as

x̂=xh,0+argminδxh{12gobsgref𝖩(xh,0)δxhΓe2+Ψ(δxh)}.
(21)

3. Bayesian reconstruction and approximation error modelling

Because of the computational limitations the discrete mapping may contain modelling errors. The general principle behind approximation error modelling is to examine the difference

ε(xh)=A(x)Ah(xh)
(22)

between an ‘exact’ model A(x) and its computational approximation Ah(xh). We thus define the ‘augmented model’

A˜h(xh)=Ah(xh)+ε(xh).
(23)

The simplest way to improve the reference measurement correction is to use the average ε̅ between many solutions of the approximative model and measurements (or accurate forward solutions) from known targets as a reference. In this case, first order statistics of the error model (EM-1) is used and the linear minimization problem is of the form

x̂=xh,0+argminδxh{12gobsAh(xh,0)𝖩(xh,0)δxhε¯Γe2+Ψ(δxh)}
(24)

where Ah(xh,0) is the solution of the linear forward model such as the Green’s function and 𝖩(xh,0) is the Jacobian calculated with optical properties xh,0.

In the approximation error approach [11

11. J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York, 2005.

,12

12. J. Kaipio and E. Somersalo, “Statistical inverse problems: Discretization, model reduction and inverse crimes,” J. Comput. Appl. Math. 198, 493–504, (2007).

], second order statistics of the error model (EM-2) is considered. Thus, both the mean and the covariance of the errors between the ‘exact’ model and approximative model are taken into account. This leads to the minimisation

x̂=argminxh{12gobsAh(xh)ε¯Γe+ε2+Ψ(xh)}
(25)

where {ε̅ε} are the mean and covariance of the distribution of ε(xh) and Γe+ε = Γe + Γε. In the linear case, this minimisation problem is written as

x̂=xh,0+argminδxh{12gobsAh(xh,0)𝖩(xh,0)δxhε¯Γe+ε2+Ψ(δxh)}.
(26)

In order to determine estimates of {ε̅ε}, samples are drawn from an appropriate prior distribution model of the unknowns and used to determine samples from both the ‘exact’ and ‘approximate’ data distributions. In [3

3. S.R. Arridge, J.P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inv. Probl. 22, 175–195, (2006).

] this approach was used to examine the model error between fine and coarse meshes, whereas in this paper it is used to examine the error between a finite element (FE) solution of the DA and Green’s function model. The computation of the approximation error statistics is explained in more detail in Section 4.1.

Remark 1 (Approximation error vs. reference data)

By comparison of Eq. (25) with Eq. (20) we see that the reference measurement correction makes the assumptions

εgrefAh(xh,0),Γε𝖨.
(27)

Remark 2 (Reference method and the Rytov approximation)

Another commonly used approach is to consider the logarithm of the forward model in Eq. (6)

Alog(x)=𝓜log𝒢(x)+e
(28)

with linearisation

Alog(x0+δx)=Alog(x0)+𝖲𝖩(x0)δx+e
(29)

with 𝖲=diag[1y0]. Introducing the reference measurement correction we find

logA˜h(xh)=logAh(xh,0)+logA(x)logA(x0)+e
(30)
=log[y0grefA(x)]+e
(31)

This is a multiplicative zeroth order correction. Linearising both these models we arrive at a scaling for the Jacobians

diag[1y0]𝖩(xh,0)diag[1gref]𝖩(x0).
(32)

This is a multiplicative correction but with a diagonal scaling matrix, rather than the covariance which is found in the approximation error model using second order statistics.

4. Results

4.1. Model implementation

4.1.1. The target distribution

We tested the approach with simulations using a cylinder domain with radius r = 35mm and height h = 110mm. The sources and detectors were placed on two rings located 6mm above and below the central xy-plane of the cylinder. Both rings contained 16 sources and 16 detectors. The target consisted of homogeneous background with two small cubic inclusions. The edge lengths of the inclusions were d = 12.2mm and the heights were h = 16.4mm. The inclusions were located in such a way that the central xy-plane of one inclusion was located at the level of the upper source-detector ring and the central xy-plane of the other inclusion was located at the level of the lower source-detector ring. The simulation domain is illustrated in Fig. 1.

We tested a situation in which the background optical properties were known and situations in which the background optical properties were mismodelled and the true homogeneous values were 10%, 20%, 30% and 40% larger or smaller than the expected values. In all of the simulations the modulation frequency of the input signal was ω = 100MHz and the refractive index of the medium was 1.56. The absorption and reduced scattering coefficients of the background medium and the inclusions are given in Table 1 for the case 1 in which the true background optical properties were known and the case 2 in which the background optical properties were -30% off from the linearisation point xh,0 = (μa,0,μs,0)T = (0.01,1)T.

Fig. 1. Simulation domain with two inclusions. Locations of sources and detectors are marked with black dots.

Table 1. The absorption coefficient μa(mm-1) and the reduced scattering coefficient μs(mm-1) of the background medium and the inclusions.

table-icon
View This Table

The data was generated using the finite element method (FEM) with the diffusion approximation as light transport model. The data types considered were logarithm of amplitude and phase shift of the measured light. Gaussian distributed noise e with standard deviation of 0.5% of corresponding amplitude and phase was added to the simulated data.

4.1.2. Computational forward models

For the ‘exact’ model A(x) we take a FE-approximation of the DA with a dense mesh. Let {ℕ,𝕋,𝕌} represent a set of nodes ℕ = {Nk;k = 1…nh}, elements 𝕋 = {τj;j = 1…ne}, and shape functions 𝕌 = {uk(r);k = 1…nh} and 𝖦 = 𝖪-1 where 𝖪 the nh × nh system matrix. We used a tetrahedra mesh with nh = 148276 nodes. For the representation of the optical parameters xh we used 8748 cubic voxels both for absorption and scattering.

The ‘approximate’ model Ah(xh) is taken to be

Ah(xh)=G(xh,0)+𝖩(xh,0)δxh
(33)

where G(xh,0) is the infinite domain Green’s function, Eq. (13). The Green’s function Jacobian 𝖩 = [𝖩(μa)𝖩(μs)] is constructed by sampling the forward and adjoint fields on the same mesh node points {Nk} :

𝖩(μa)=(G*(rk,rd,i)·G(rk,rs,j))
(34)

and

𝖩(D)=(G*(rk,rd,i)·G(rk,rs,j))
(35)

where G* is the adjoint Green’s function and utilising the chain rule to obtain the derivatives for the scattering. In all of the approaches, the Jacobian was calculated using the homogeneous optical properties μa,0 = 0.01 mm-1 and μs,0 = 1mm-1.

4.1.3. Approximation error statistics

The approximation error mean and covariance were calculated similarly as in [14

14. V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S.R. Arridge, and J.P. Kaipio, “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A 26, 2257–2268, (2009).

]. Thus, we used a proper smoothness prior model as the prior distribution π(x) and draw r = 200 samples x(ℓ)h from it. The correlation length for both absorption and scattering was set as 14 mm which means that this is (roughly) the prior estimate of the spatial size of the inhomogeneities in the target domain. The prior means for absorption and scattering were set as μ̄a = 0.01 mm-1 and μ̄s = 1mm-1, and the marginal variances for the inhomogeneity and background part were set such as 2 s.t.d. limits for the contrast of the inhomogeneities and background corresponded to 1 % of the mean values for both absorption and scattering. The approximation error was computed from the samples as

ε=A(x())Ah(xh())
(36)

where the accurate model was the FE-solution of the DA and the reduced model was the solution of the Green’s function, Eq. (33). Furthermore, the mean and the covariance of the approximation error were computed as

ε¯=1r=1rεh()
(37)
Γε=1r1=1rε()ε()TεεT.
(38)

We note that the variance of the contrast of the inhomogeneities in the prior model was chosen unrealistically small since the assumption behind the linear model needed to be fulfilled when creating the approximation error statistics. This however, as shown later in Sec. 4.2, does not limit the applicability of the approach to reconstruct larger variations in optical parameters. For more information about the Gaussian smoothness prior model, see e.g. [3

3. S.R. Arridge, J.P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inv. Probl. 22, 175–195, (2006).

,14

14. V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S.R. Arridge, and J.P. Kaipio, “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A 26, 2257–2268, (2009).

].

Furthermore, we note that the setting up the error model is computationally intensive task and the time required is roughly the number of samples r multiplied by the time for the forward solutions of the accurate and approximative models [14

14. V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S.R. Arridge, and J.P. Kaipio, “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A 26, 2257–2268, (2009).

]. The error model, however, needs to be estimated only once for a fixed measurement setup and it can be done off-line.

4.2. Reconstructions

Reconstructions were calculated using three different approaches: the reference measurement correction (RM), the approximation error using first order statistics (EM-1), and the approximation error using second order statistics (EM-2).

Fig. 2. Horizontal slices at the heights z = 59mm (two columns from the left) and z = 50mm (two columns from the right) from 3D reconstructions of absorption (left) and scattering (right) distributions with known background. The images from top to bottom are: target distribution (top row), and reconstructions using RM (second row), EM-1 (third row) and EM-2 (fourth row).

In all of the reconstructions, the penalty functional Ψ was the proper smoothness prior model with otherwise same parameters as in creating the approximation error statistics except that the marginal variances for the inhomogeneity were set such as 2 s.t.d. limits for the contrast of the inhomogeneities corresponded to 40 % of the mean values for absorption and scattering and the marginal variances for the background part were set such as 2 s.t.d. limits for the contrast of the background corresponded to 20% of the mean values for absorption and scattering which is more realistic choice when considering absorption and scattering variations typical in DOT.

The reconstructions using the RM correction were obtained as the solution of the minimisation problem (21) where the gobs was the simulated measurement data, gref was the reference data which is the FE-solution of the DA with optical properties of the linearisation point xh,0 which were for the absorption and scattering μa = 0.01mm-1 and μs = 1mm-1 throughout the domain. Further, 𝖩(xh,0) was the Jacobian which was calculated as described in Sec. 4.1, Eqs. (34) and (35),using optical parameters of the linearisation point xh,0 (μa = 0.01mm-1 and μs = 1mm-1).

The reconstructions using the EM-1 were obtained by solving the minimisation problem (24) and the reconstructions using the EM-2 were obtained by solving the minimisation problem (26). In the approaches, the approximative model Ah(xh) was the infinite domain Green’s function, Eq. (33) and the Jacobian 𝖩(xh,0) was calculated similarly as in the RM. The mean and the covariance of the approximation error statistics were constructed as described in Section 4.1.

4.2.1. Reconstructions when the background properties are known

As the first case we investigated the situation in which the background optical properties were known, thus xh,0 = x0, see case 1 of Table 1. Two horizontal and two vertical slices of the 3D reconstructions are shown in Figs. 2 and 3. The locations of the slices were chosen such that they cut through the inclusions.

Fig. 3. Vertical slices at the depths x = 21mm (two columns from the left) and x = 45mm (two columns from the right) from 3D reconstructions of absorption (left) and scattering (right) distributions with known background. The images from top to bottom are: target distribution (top row), and reconstructions using RM (second row), EM-1 (third row), and EM-2 (fourth row).

As it can be seen, RM, EM-1 and the EM-2 give similar results. In all of the reconstructions, the inclusions can be well distinguished and they look similar to the target distribution. Furthermore, in all of the reconstructions, the estimated absorption and scattering parameters are close to the original values. There exists some crosstalk between the absorption and scatter images. This, however, is typical for DOT reconstructions and utilising more sophisticated approximation error and prior models does not seem to reduce the amount of the cross talk in linear reconstruction as was shown to happen with the DA based reconstructions in [3

3. S.R. Arridge, J.P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inv. Probl. 22, 175–195, (2006).

].

Thus, the results show that the linear reconstruction method with a reference measurement correction can reconstruct the absorption and scattering distributions well if the background optical properties are known. On the other hand, the approximation error modelling, in which statistical properties of the modelling error are utilised instead of a reference measurement, can be used to obtain as good quality reconstructions as the reference approach.

Fig. 4. Horizontal slices at the heights z = 59mm (two columns from the left) and z = 50mm (two columns from the right) from 3D reconstructions of absorption (left) and scattering (right) distributions with mismodelled background. The images from top to bottom are: target distribution (top row), and reconstructions using RM (second row), EM-1 (third row), and EM-2 (fourth row).

4.2.2. Reconstructions with mismodelled background

As the second case we investigated a situation in which the background optical properties were mismodelled. The measurement data was simulated using optical properties where the true values of the homogeneous background x0 were 10%, 20%, 30% and 40% larger or smaller than the linearisation point (μa = 0.01mm-1 and μs = 1mm-1). The absorption and scattering values of the case 2, in which the background optical properties are 30% smaller than the linearisation point, are given in Table 1. The reference measurement gref, Jacobian 𝖩(xh,0) and the approximation error statistics were calculated using the expected optical properties xh as described in Sec. 4.1.

Thus, the results show that when the background optical properties are mismodelled, the RM and the EM-1 which both utilise only the mean of the statistical information of the error model, do not produce good quality reconstructions. On the other hand, EM-2, which takes into account correlation of the modelling errors, is capable of estimating the optical parameters with good accuracy even if the background values are 30% off from the expected ones.

Fig. 5. Vertical slices at the depths x = 21mm (two columns from the left) and x = 45mm (two columns from the right) from 3D reconstructions of absorption (left) and scattering (right) distributions with mismodelled background. The images from top to bottom are: target distribution (top row), and reconstructions using RM (second row), EM-1 (third row), and EM-2 (fourth row).

4.2.3. Comparison and discussion

To compare the capability of the different approaches to estimate the values of optical parameters, the relative errors between the target and estimated optical properties, x and x̂, respectively, were calculated as

êx=(xx̂)2x2·100
(39)

for different levels of deviation in optical parameters between the true homogeneous background and the linearisation point. The relative errors of different models against the difference in background value are shown in Fig. 6.

The figure shows that when the background optical properties are known, all three approaches give similar reconstruction errors. When the error in background value compared to the linearisation point is increased, the errors of the RM and EM-1 estimates increase clearly. The errors of the EM-2 show also some increase but significantly less than the other approaches. Thus, the EM-2 method gives clearly more accurate estimates than the two other approaches and can compensate the errors caused by using mismodelled background optical properties in linear reconstruction.

Fig. 6. Relative errors between the target optical properties and estimated optical properties against the difference in background value calculated with the RM (∗), EM-1 (∘) and EM-2 (+).

5. Conclusions

In this paper we examined the correction of errors when using a first order Born approximation with an infinite space Greens function model as the basis for linear reconstruction in DOT, when real data is generated on a finite domain with possibly unknown background optical properties. We considered the relationship between conventional reference measurement correction and the approximation error modelling in DOT reconstruction. The approaches were tested with simulations in cases in which the background optical properties were known and mismodelled. The results show that when the background optical properties are known, the conventional reference correction method and the approximation error approaches using first and second order statistics give similar results, and thus all of the approaches can be used in DOT reconstruction with a linear model. However, when the background optical properties are mismodelled, the conventional reference method and the approximation error approach using first order statistics fail to give good estimates of the target optical properties. In this case, the approximation error approach using second order statistics is capable of providing feasible reconstructions of the optical properties even when the background values are 30% smaller than the expected values. Thus, the approximation error approach using second order statistics can compensate the errors caused by inaccurately known background optical properties in linear DOT reconstruction.

Acknowledgments

This work was supported by the Academy of Finland (projects 122499, 119270, 140731 and 213476, Finnish Center of Excellence in Inverse Problems Research) and by EPSRC grant EP/E034950/1.

References and links

1.

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

2.

S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37, 1531–1560, (1992).

3.

S.R. Arridge, J.P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, “Approximation errors and model reduction with an application in optical diffusion tomography,” Inv. Probl. 22, 175–195, (2006).

4.

S.R. Arridge and J.C. Schotland, “Optical tomography: forward and inverse problems,” Inv. Probl. 25, 123010, (2009).

5.

G. Bal, “Inverse transport theory and applications,” Inv. Probl. 25, 053001, (2009).

6.

D. A. Boas, “A fundamental limitation of linearized algorithms for diffuse optical tomography,” Opt. Express 1, 404–413, (1997).

7.

A. Gibson, J.C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical tomography,” Phys. Med. Biol. 50, R1–R43, (2005).

8.

J. Heino, E. Somersalo, and J.P. Kaipio, “Compensation for geometric mismodelling by anisotropies in optical tomography,” Opt. Express 13, 296–308, (2005).

9.

J.M.J. Huttunen and J.P. Kaipio, “Approximation errors in nonstationary inverse problems,” Inverse Problems and Imaging 1, 77–93, (2007).

10.

J.M.J. Huttunen and J.P. Kaipio, “Model reduction in state identification problems with an application to determination of thermal parameters,” Applied Numerical Mathematics 59, 877–890, (2009).

11.

J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York, 2005.

12.

J. Kaipio and E. Somersalo, “Statistical inverse problems: Discretization, model reduction and inverse crimes,” J. Comput. Appl. Math. 198, 493–504, (2007).

13.

A. D. Klose and A. H. Hielscher, “Optical tomography with the equation of radiative transfer,” International Journal of Numerical Methods for Heat & Fluid Flow 18, 443–464, (2008).

14.

V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S.R. Arridge, and J.P. Kaipio, “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A 26, 2257–2268, (2009).

15.

S.D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, “Imaging complex structures with diffuse light,” Opt. Express 16, 5048–5060, (2008).

16.

A. Lehikoinen, S. Finsterle, A. Voutilainen, L.M. Heikkinen, M. Vauhkonen, and J.P. Kaipio, “Approximation errors and truncation of computational domains with application to geophysical tomography,” Inverse Problems and Imaging 1, 371–389, (2007).

17.

A. Nissinen, L.M. Heikkinen, and J.P Kaipio, “The Bayesian approximation error approach for electrical impedance tomography - experimental results,” Meas. Sci. Technol. 19, 015501, (2008).

18.

A. Nissinen, L.M. Heikkinen, V. Kolehmainen, and J.P Kaipio, “Compensation of errors due to discretization, domain truncation and unknown contact impedances in electrical impedance tomography,” Meas. Sci. Technol. 20, 015504, (2009).

19.

A. Nissinen, V. Kolehmainen, and J.P Kaipio, “Compensation of modelling errors due to unknown domain boundary in electrical impedance tomography,” IEEE Trans. Med. Imag, Submitted.

20.

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

21.

J. C. Schotland and V. Markel, “Inverse scattering with diffusing waves,” J. Opt. Soc. Am. A 18, 2767–2777, (2001).

22.

J. Ripoll, V. Ntziachristos, and M. Nieto-Vesperinas, “The Kirchhoff approximation for diffusive waves,” Phys. Rev. E 64, 1–8, (2001).

23.

N. Polydorides, “Linearization error in electical impedance tomography,” Progress In Electromagnetics Research, PIER 93, 323–337, (2009).

24.

T. Tarvainen, V. Kolehmainen, A. Pulkkinen, M. Vauhkonen, M. Schweiger, S.R. Arridge, and J.P. Kaipio, “Approximation error approach for compensating for modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography,” Inv. Probl. 26, 015005, (2010).

OCIS Codes
(100.3190) Image processing : Inverse problems
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.6960) Medical optics and biotechnology : Tomography
(290.7050) Scattering : Turbid media

ToC Category:
Image Reconstruction and Inverse Problems

History
Original Manuscript: June 7, 2010
Revised Manuscript: July 9, 2010
Manuscript Accepted: July 9, 2010
Published: July 16, 2010

Virtual Issues
Optical Imaging and Spectroscopy (2010) Biomedical Optics Express

Citation
Tanja Tarvainen, Ville Kolehmainen, Jari P. Kaipio, and Simon R. Arridge, "Corrections to linear methods for diffuse optical tomography using approximation error modelling," Biomed. Opt. Express 1, 209-222 (2010)
http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-1-1-209


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. S. R. Arridge, "Optical tomography in medical imaging," Inv. Probl. 15, R41-R93 (1999).
  2. S. R. Arridge, M. Cope, and D. T. Delpy, "The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis," Phys. Med. Biol. 37, 1531-1560 (1992).
  3. S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, T. Tarvainen, and M. Vauhkonen, "Approximation errors and model reduction with an application in optical diffusion tomography," Inv. Probl. 22, 175-195 (2006).
  4. S. R. Arridge and J. C. Schotland, "Optical tomography: forward and inverse problems," Inv. Probl. 25, 123010 (2009).
  5. G. Bal, "Inverse transport theory and applications," Inv. Probl. 25, 053001 (2009).
  6. D. A. Boas, "A fundamental limitation of linearized algorithms for diffuse optical tomography," Opt. Express 1, 404-413 (1997).
  7. A. Gibson, J. C. Hebden, and S. R. Arridge, "Recent advances in diffuse optical tomography," Phys. Med. Biol. 50, R1-R43 (2005).
  8. J. Heino, E. Somersalo, and J. P. Kaipio, "Compensation for geometric mismodelling by anisotropies in optical tomography," Opt. Express 13, 296-308 (2005).
  9. J. M. J. Huttunen and J. P. Kaipio, "Approximation errors in nonstationary inverse problems," Inv. Probl. Imaging 1, 77-93 (2007).
  10. J. M. J. Huttunen and J. P. Kaipio, "Model reduction in state identification problems with an application to determination of thermal parameters," Appl. Numer. Math. 59, 877-890 (2009).
  11. J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, (Springer, New York, 2005).
  12. J. Kaipio and E. Somersalo, "Statistical inverse problems: Discretization, model reduction and inverse crimes," J. Comput. Appl. Math. 198, 493-504 (2007).
  13. A. D. Klose and A. H. Hielscher, "Optical tomography with the equation of radiative transfer," Int. J. Numer. Meth. Heat Fluid Flow 18, 443-464 (2008).
  14. V. Kolehmainen, M. Schweiger, I. Nissilä, T. Tarvainen, S. R. Arridge, and J. P. Kaipio, "Approximation errors and model reduction in three-dimensional diffuse optical tomography," J. Opt. Soc. Am. A 26, 2257-2268 (2009).
  15. S. D. Konecky, G. Y. Panasyuk, K. Lee, V. Markel, A. G. Yodh, and J. C. Schotland, "Imaging complex structures with diffuse light," Opt. Express 16, 5048-5060 (2008).
  16. A. Lehikoinen, S. Finsterle, A. Voutilainen, L. M. Heikkinen, M. Vauhkonen, and J. P. Kaipio, "Approximation errors and truncation of computational domains with application to geophysical tomography," Inv. Probl. Imaging 1, 371-389 (2007).
  17. A. Nissinen, L. M. Heikkinen, and J. P Kaipio, "The Bayesian approximation error approach for electrical impedance tomography - experimental results," Meas. Sci. Technol. 19, 015501 (2008).
  18. A. Nissinen, L. M. Heikkinen, V. Kolehmainen, and J. P Kaipio, "Compensation of errors due to discretization, domain truncation and unknown contact impedances in electrical impedance tomography," Meas. Sci. Technol. 20, 015504 (2009).
  19. A. Nissinen, V. Kolehmainen, and J.P Kaipio, "Compensation of modelling errors due to unknown domain boundary in electrical impedance tomography," IEEE Trans. Med. Imag, Submitted.
  20. 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).
  21. J. C. Schotland and V. Markel, "Inverse scattering with diffusing waves," J. Opt. Soc. Am. A 18, 2767-2777 (2001).
  22. J. Ripoll, V. Ntziachristos, and M. Nieto-Vesperinas, "The Kirchhoff approximation for diffusive waves," Phys. Rev. E 64, 1-8 (2001).
  23. N. Polydorides, "Linearization error in electical impedance tomography," Prog. Electromagn. Res., PIER 93, 323-337 (2009).
  24. T. Tarvainen, V. Kolehmainen, A. Pulkkinen, M. Vauhkonen, M. Schweiger, S. R. Arridge, and J. P. Kaipio, "Approximation error approach for compensating for modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography," Inv. Probl. 26, 015005 (2010).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited