Uniqueness, Born Approximation, and Numerical Methods for Diffuse Optical Tomography

Abstract

Diffuse optical tomogrpahy (DOT) is to find optical coefficients of tissue using near infrared light. DOT as an inverse problem is described and the studies about unique determination of optical coefficients are summarized. If a priori information of the optical coefficient is known, DOT is reformulated to find a perturbation of the optical coefficients inverting the Born expansion which is an infinite series expansion with respect to the perturbation and the a priori information. Numerical methods for DOT are explained as methods inverting first- or second-order Born approximation or the Born expansion itself.

1. Introduction

DOT is to find optical coefficients of tissue using near infrared light. DOT is known to be of low cost, portable, nonionized, and nonmagnetized. And DOT has higher temporal resolution and more functional information than conventional structural medical imaging modalities such as magnetic resonance imaging (MRI) and computerized tomography (CT). For the comparison to other functional imaging modalities such as functional MRI (fMRI), photon emission tomography (PET), and electroencephalogram (EEG), see [1]. DOT is used in the area of breast imaging [2–4], functional neuroimaging [5, 6], brain computer interface (BCI) [7, 8], and the study about seizure [9, 10], newborn infants [11, 12], osteoarthritis [13], and rat brain [14, 15].

In this paper, DOT is explained as an inverse problem with respect to a forward problem formulated as an elliptic partial differential equation. Propagation of light in biological tissues is usually described by diffusion approximation equation in the frequency domain, the simplest but nontrivial approximation of the Boltzmann equation, as follows:−∇⋅(𝜅∇Φ)+𝜇_𝑎+𝑖𝜔𝑐()Φ=𝑞inΩ,Φ+2𝑎𝜈⋅𝜅∇Φ=0on𝜕Ω, (1a)(1b)where 𝜈 is an outer unit normal vector, Φ is a photon density distribution, 𝜅=1/3(𝜇_𝑎+𝜇^′_𝑠) is an diffusion coefficient, 𝜇_𝑎 is an absorption coefficient, 𝜇^′_𝑠 is a reduced scattering coefficient, and 𝑎 is a reflection coefficient.

If 𝜕𝜅/𝜕𝜈=0 on 𝜕Ω, by setting Ψ=√𝜅Φ and 𝑘=Δ√𝜅/√𝜅+𝜇_𝑎/𝜅+𝑖(𝜔/𝑐𝜅) with 𝐼𝑚(𝑘)≥0, we have−ΔΨ+𝑘²𝑞Ψ=√𝜅()inΩ,Ψ+2𝑎𝜈⋅𝜅∇Ψ=0on𝜕Ω. (2a)(2b)

If 𝑘 is constant and 𝑞(⋅)/√𝜅=𝛿(⋅,𝑟_𝑠) for some source point 𝑟_𝑠, we have the following solution of (2a): Ψ(𝑟)=𝑅𝑟,𝑟_𝑠=𝑒^{𝑖𝑘|𝑟−𝑟𝑠|}4𝜋𝑟−𝑟_𝑠. (3)

DOT is to find the optical coefficients 𝜇_𝑎 and/or 𝜇^′_𝑠 from the measurement information Φ_𝑖,𝑗 which is the value of the solution of (1a) and (1b) at 𝑟_𝑖∈𝜕Ω when 𝑞(𝑟)=𝛿(𝑟,𝑟_𝑗),𝑟_𝑗∈𝜕Ω. The 𝑟_𝑖 and 𝑟_𝑗 are usually called source and detector point, respectively.

In Section 2, the unique determination of the optical coefficients is discussed and many known results are summarized for the uniqueness questions. In Section 3, DOT is reformulated as to find perturbation of the optical coefficient inverting the Born expansion. The errors of the Born approximation in the Lebesgue and Sobolev norms are given. In Section 4, numerical methods of DOT are mainly described as the methods inverting the first-, second-, and higher-order Born approximation and Born expansion itself.

2. Uniqueness

The research about unique determination of the optical coefficients in DOT is rare except [16], but it is a very important issue for DOT as an inverse problem. The determinaiton of optical coefficients (𝜇_𝑎,𝜇_𝑠) in (1a) and (1b) is equivalent to the determination of 𝑘 in (2a) and (2b) when 𝜔≠0.

When 𝜇_𝑎,𝜅,∇𝜅 have upper and lower bound and 𝑞 is contained in 𝐻⁻¹(Ω) or a Dirac delta function, (1a) and (1b) have a unique solution Φ∈𝐻¹(Ω) and (2a) and (2b) have a unique solution Ψ∈𝐻¹(Ω) [17, 18].

Boundary value problem (1a) and (1b) with 𝑞(𝑟)=𝛿(𝑟,𝑟_𝑠) is equivalent to boundary value problem with (1a) for 𝑞=0 and nonzero Robin boundary condition replacing (1b). This argument can be proved using the function 𝐻 in [19]. Therefore, DOT is redescribed as to find the optical coefficients from the Robin-to-Dirichlet map defined as a map from 𝐻^−1/2(𝜕Ω) to 𝐻^1/2(𝜕Ω). Using unique solvability of (1a) with the Dirichlet or Neumann boundary condition replacing (1b), the Robin-to-Dirichlet map is equivalent to the Neumann-to-Dirichlet map and to the Dirichlet-to-Neumann map. Since 𝐼𝑚(𝑘)>0 for 𝜔≠0, 𝑘 is not a Dirichlet eigenvalue of (2a). Therefore, knowing the Dirichlet-to-Neumann map is also equivalent to knowing farfield map in inverse scattering problem [20, 21].

Using above results and many known results for inverse scattering problem, we summarized the uniqueness and nonuniqueness results as follows.

Case 1 (𝜅=1,𝜇_𝑎=1+(𝑚−1)𝜒_𝐷,𝑚≠1). The cases for positive constant 𝑚 can be understood as special cases of Case 2. The limiting cases 𝑚=∞ and 𝑚=0 could be considered as (1a) in Ω⧵𝐷 with the Robin boundary condition (1b) on 𝜕Ω and the boundary condition on 𝜕𝐷 as Φ=0 (sound-soft case) and 𝜕Φ/𝜕𝜈=0 (sound-hard case) on 𝜕𝐷. The uniqueness for sound-soft and sound-hard obstacle 𝐷 is considered in [20, 22].

Case 2 (𝜅=1,𝜇_𝑎=1+(𝑚(𝑥)−1)𝜒_𝐷, where 𝑚(𝑥)≠1 on 𝜕𝐷). This case is called “inverse transmission problem” and the uniqueness is solved in [20].

Case 3 (𝜅=1,𝜇_𝑎=𝜇_𝑎(𝑥)). In [16], nonuniqueness of 𝜇_𝑎(𝑥) is shown by assuming that refractive index is not determined or we use only continuous wave light source. However, if refractive index is known and 𝜔≠0, unique determination of the optical coefficient is possible. This nonuniqueness is proved by using [23].

Case 4 (𝜅=𝐼_𝑛+(𝐾(𝑥)−𝐼_𝑛)𝜒_𝐷,𝜇_𝑎=1+(𝑚(𝑥)−1)𝜒_𝐷). Here, 𝐾(𝑥) is a unknown positive-definite matrix function such that 𝐾(𝑥)≠𝐼_𝑛 on 𝜕𝐷 and 𝑚(𝑥) is a positive function such that 𝑚(𝑥)≠1 on 𝜕𝐷. The uniqueness of 𝐷 is solved in [21, 24–28] and the nonuniqueness of 𝐾(𝑥) is reported in [24, 29, 30]. Therefore, although the domain of nonhomogenity 𝐷 can be uniquely determined by infinite measurements, the nonhomogeneous anisotropic diffusion coefficient 𝐾(𝑥) cannot be determined uniquely. Similar results are known for the nonuniqueness and illusion for anisotropic nonhomogeneous electric conductivity [24, 31].

3. Born Approximation

Suppose that we know a priori information 𝑥⁰=(𝜅⁰,𝜇⁰_𝑎) about the optical coefficients 𝑥=(𝜅,𝜇_𝑎) to find. Then, DOT is redescribed as to find the perturbation 𝛿𝑥=𝑥−𝑥⁰=(𝛿𝜅,𝛿𝜇_𝑎).

Let Φ and Φ⁰ be the solutions of (1a) and (1b) for 𝑥 and 𝑥⁰, respectively. When 𝑞(𝑟)=𝛿(𝑟,𝑟_𝑠), Φ and Φ⁰ are the Robin function 𝑅(⋅,𝑟_𝑠) and 𝑅⁰(⋅,𝑟_𝑠) for the optical coefficients 𝑥 and 𝑥⁰, respectively. 𝑅(⋅,𝑟_𝑠) is expanded as the Born expansion, an infinite series, with respect to 𝑅⁰(⋅,𝑟_𝑠) and 𝛿𝑥 as follows: 𝑅𝑟,𝑟_𝑠=𝑅⁰𝑟,𝑟_𝑠+ℛ𝑅⁰⋅,𝑟_𝑠+ℛ²𝑅⁰⋅,𝑟_𝑠+⋯, (4) where (ℛΨ)(𝑟)=ℛ(𝛿𝑥)Ψ(𝑟)=ℛ₁Ψ(𝑟)+ℛ₂Ψ(𝑟),ℛ₁Ψ(𝑟)=ℛ₁𝛿𝜇_𝑎Ψ(𝑟)=−_Ω𝛿𝜇_𝑎𝑟^′𝑅𝑟,𝑟^′Ψ𝑟^′𝑑𝑟^′,ℛ₂Ψ(𝑟)=ℛ₂(𝛿𝜅)Ψ(𝑟)=−_Ω𝛿𝜅𝑟^′∇𝑅𝑟,𝑟^′∇Ψ𝑟^′𝑑𝑟^′. (5)

It is proved that the 𝑛th order term in the above Born expansion (4) is the 𝑛th order Frechet derivative divided by 𝑛! [18]. The bounds for the operators ℛ₁ and ℛ₂ are given in the Sobolev, Lebesque, and weighted Sobolev spaces norms [18]. The estimate for the Lebesque space norms for ℛ₁ is also given in [32]. Therefore, there exist positive constants 𝐶₁ and 𝐶₂ such thatℛ_{1𝑤1,𝑝(Ω)→𝑤1,𝑝(Ω)}≤𝐶₁𝛿𝜇_{𝑎𝐿∞(Ω)},ℛ_{2𝑤1,𝑝(Ω)→𝑤1,𝑝(Ω)}≤𝐶₂‖𝛿𝜅‖_𝐿^∞_(Ω), (6a)(6b)for 𝑝≥1. The detailed estimate for the coefficients 𝐶₁,𝐶₂ in special cases is given in [18, 32].

If 𝐶₃=𝐶₁𝛿𝜇_𝑎∞+𝐶₂‖𝛿𝜅‖_∞<1 (7) holds, then the error of the 𝑚th order Born approximation 𝐵^𝑚𝑅⁰(⋅,𝑟_𝑠) is given by 𝑅⋅,𝑟_𝑠−𝐵^𝑚𝑅⁰⋅,𝑟_{𝑠𝑊1,𝑝}≤𝐶^𝑚₃𝑅(⋅,𝑟_𝑠)_{𝑊1,𝑝(Ω)} (8) using the estimates (6a) and (6b).

4. Numerical Methods

Numerical methods for DOT are explained in terms of solving the Born approximation given in Section 3. Other numerical methods not categorized as inverting Born approximation are commented in the final section.

4.1. Linearized Methods

Linearized DOT to find 𝛿𝑥 by solving the first-order Born approximation is studied by many researchers in the initial stage of DOT. Although this method lacks exact recovery in most cases, it is used still frequently when faster real-time computing is needed and very good a priori information is given. In this method, the discretized problem is an algebraic equation and it is essential to use efficient matrix solver. In the Jacobian matrix, the number of rows is the number of measurements and the number of columns is the number of elements to be determined. In most cases, the number of elements is larger than the number of measurements in order to obtain higher resolution image. Therefore, the algebraic equation is usually an underdetermined system. Efficient matrix solvers including arithmetic reconstruction technique (ART), simultaneous arithmetic reconstruction technique (SART), simultaneous iterative reconstruction technique (SIRT), and the Krylov space methods are studied [33, 34]. Since the linearized problems also have ill-posed property of DOT, there are many researches about regularization methods including the Tikhonov regularization. The dependence on the discretization error due to ill posedness of the linearized method is reported in [35].

4.2. Nonlinear Methods

Nonlinear method is to find 𝛿𝑥 by solving the Born expansion [36]. This method is usually formulated as optimization problem and is solved by the Newton-type method including Levenberg-Marquardt method [37, 38]. A few softwares based on nonlinear method with finite element forward solver with corresponding references are as follows: (i)TOAST (Time-Resolved Optical Absorption and Scattering Tomography) [39], (ii)NIRFAST (Near InfraRed Florescence and Spectral Tomography) [40, 41], (iii)PMI (Photon Migration Imaging) Toolbox [42]. Nonlinear method needs heavy computation due to large iteration numbers. To reduce the heavy computation, there are studies about efficient numerical techniques such as multigrid, domain decomposition [43], and adaptive [44] method.

4.3. Inverse Born Approximation of Order Higher Than Two

This method is to find 𝛿𝑥 by solving second- and higher-order Born approximation. In fact, solving the Born approximation of order higher than two is implicit but can be approximated by explicit inverse Born approximation. Formal inverse of the Born approximation is called inverse Born approximation. The first-order inverse Born approximation corresponds to linearized DOT and inverse Born expansion itself corresponds to (nonlinear) DOT. Higher-order methods improve the order of convergence for lower-order methods. The error of the inverse Born approximation is given and analyzed in terms of the Lebesgue space norms when 𝛿𝜅=0 [45] and in terms of the Sobolev space norms for the second order [46].

4.4. Other Methods

The solution of (1a) and (1b) is usually solved by finite element method [36]. Another approach for the forward problem is to compute directly the Robin function using the Fourier-Laplace transform [47–50]. The disadvantage of this method is that it depends on the special geometry of the region of interest and the inverse Fourier-Laplace transfrom is known to be severely ill posed.

Equations (1a) and (1b) could be replaced by probabilistic approach in the Monte Carlo method [51]. The method takes much more time than finite element method and highly depends on random number generator. The comparison of finite element method and the Monte Carlo method is done by many papers including [52].

The diffusion approximation (1a) and (1b) is the first-order approximation of radiative transfer equation. There are studies about DOT based on radiative transfer equation [53, 54] and its 𝑛th order approximation [55].

5. Conclusion

Unique determination of DOT is surveyed. The study about nonuniqueness for anisotropic diffusion coefficients and for unknown refractive index is also surveyed. The perturbation of photon density with respect to the perturbed optical coefficient is expanded using the Born expansion and the error analysis removing higher-order terms is given. The numerical methods for DOT are described by inverting first-, second-, and higher-order Born approximation, and the Born expansion itself is reviewed.

Acknowledgment

This work is supported by Basic Science Research Program through the National Research foundation of Korea (NRF) funded by the Ministry of Education (2010-0004047).