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

Authors: Kiwoon Kwon
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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:(1a)-∇·(κ∇Φ)+(μa+iωc)Φ=qinΩ,(1b)Φ+2aν·(κ∇Φ)=0on∂Ω,where ν is an outer unit normal vector, Φ is a photon density distribution, κ=1/3(μa+μs′) is an diffusion coefficient, μa is an absorption coefficient, μs′ is a reduced scattering coefficient, and a is a reflection coefficient.If∂κ/∂ν=0 on ∂Ω, by setting Ψ=κΦ and k=Δκ/κ+μa/κ+i(ω/cκ) with Im(k)≥0, we have(2a)-ΔΨ+k2Ψ=qκinΩ,(2b)Ψ+2aν·(κ∇Ψ)=0on∂Ω.Ifk is constant and q(·)/κ=δ(·,rs) for some source point rs, we have the following solution of (2a):
(3)Ψ(r)=R(r,rs)=eik|r-rs|4π|r-rs|.DOT is to find the optical coefficientsμa and/or μs′ from the measurement information Φi,j which is the value of the solution of (1a) and (1b) at ri∈∂Ω when q(r)=δ(r,rj),rj∈∂Ω. The ri and rj are usually called source and detector point, respectively.In Section2, 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 (μa,μs) in (1a) and (1b) is equivalent to the determination of k in (2a) and (2b) when ω≠0.Whenμa,κ,∇κ have upper and lower bound and q is contained in H-1(Ω) or a Dirac delta function, (1a) and (1b) have a unique solution Φ∈H1(Ω) and (2a) and (2b) have a unique solution Ψ∈H1(Ω) [17, 18].Boundary value problem (1a) and (1b) with q(r)=δ(r,rs) is equivalent to boundary value problem with (1a) for q=0 and nonzero Robin boundary condition replacing (1b). This argument can be proved using the function H in [19]. Therefore, DOT is redescribed as to find the optical coefficients from the Robin-to-Dirichlet map defined as a map from H-1/2(∂Ω) to H1/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 Im(k)>0 for ω≠0, k 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,μa=1+(m-1)χD,m≠1).
The cases for positive constantm can be understood as special cases of Case 2. The limiting cases m=∞ and m=0 could be considered as (1a) in Ω∖D¯ with the Robin boundary condition (1b) on ∂Ω and the boundary condition on ∂D as Φ=0 (sound-soft case) and ∂Φ/∂ν=0 (sound-hard case) on ∂D. The uniqueness for sound-soft and sound-hard obstacle D is considered in [20, 22].Case 2 (κ=1,μa=1+(m(x)-1)χD, where m(x)≠1 on ∂D).
This case is called “inverse transmission problem” and the uniqueness is solved in [20].Case 3 (κ=1,μa=μa(x)).
In [16], nonuniqueness of μa(x) 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 (κ=In+(K(x)-In)χD,μa=1+(m(x)-1)χD).
Here,K(x) is a unknown positive-definite matrix function such that K(x)≠In on ∂D and m(x) is a positive function such that m(x)≠1 on ∂D. The uniqueness of D is solved in [21, 24–28] and the nonuniqueness of K(x) is reported in [24, 29, 30]. Therefore, although the domain of nonhomogenity D can be uniquely determined by infinite measurements, the nonhomogeneous anisotropic diffusion coefficient K(x) 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 informationx0=(κ0,μa0) about the optical coefficients x=(κ,μa) to find. Then, DOT is redescribed as to find the perturbation δx=x-x0=(δκ,δμa).LetΦ and Φ0 be the solutions of (1a) and (1b) for x and x0, respectively. When q(r)=δ(r,rs), Φ and Φ0 are the Robin function R(·,rs) and R0(·,rs) for the optical coefficients x and x0, respectively. R(·,rs) is expanded as the Born expansion, an infinite series, with respect to R0(·,rs) and δx as follows:
(4)R(r,rs)=R0(r,rs)+ℛR0(·,rs)+ℛ2R0(·,rs)+⋯,
where
(5)(ℛΨ)(r)=ℛ(δx)Ψ(r)=(ℛ1Ψ)(r)+(ℛ2Ψ)(r),(ℛ1Ψ)(r)=ℛ1(δμa)Ψ(r)=-∫Ωδμa(r′)R(r,r′)Ψ(r′)dr′,(ℛ2Ψ)(r)=ℛ2(δκ)Ψ(r)=-∫Ωδκ(r′)∇R(r,r′)∇Ψ(r′)dr′.It is proved that thenth order term in the above Born expansion (4) is the nth order Frechet derivative divided by n! [18]. The bounds for the operators ℛ1 and ℛ2 are given in the Sobolev, Lebesque, and weighted Sobolev spaces norms [18]. The estimate for the Lebesque space norms for ℛ1 is also given in [32]. Therefore, there exist positive constants C1 and C2 such that(6a)∥ℛ1∥w1,p(Ω)→w1,p(Ω)≤C1∥δμa∥L∞(Ω),(6b)∥ℛ2∥w1,p(Ω)→w1,p(Ω)≤C2∥δκ∥L∞(Ω),for p≥1. The detailed estimate for the coefficients C1,C2 in special cases is given in [18, 32].If(7)C3=C1∥δμa∥∞+C2∥δκ∥∞<1
holds, then the error of the mth order Born approximation BmR0(·,rs) is given by
(8)∥R(·,rs)-BmR0(·,rs)∥W1,p≤C3m∥R(·,rs)∥W1,p(Ω)
using the estimates (6a) and (6b).


## 4. Numerical Methods

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


### 4.1. Linearized Methods

Linearized DOT to findδx 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δx 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δx 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 nth order approximation [55].


## 4.1. Linearized Methods

Linearized DOT to findδx 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δx 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δx 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 nth 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.



---

*Source: 824501-2013-06-03.xml*