The Second-Order Born Approximation in Diffuse Optical Tomography

Abstract

Diffuse optical tomography is used to find the optical parameters of a turbid medium with infrared red light. The problem is mathematically formulated as a nonlinear problem to find the solution for the diffusion operator mapping the optical coefficients to the photon density distribution on the boundary of the region of interest, which is also represented by the Born expansion with respect to the unperturbed photon densities and perturbed optical coefficients. We suggest a new method of finding the solution by using the second-order Born approximation of the operator. The error analysis for the suggested method based on the second-order Born approximation is presented and compared with the conventional linearized method based on the first-order Born approximation. The suggested method has better convergence order than the linearized method, and this is verified in the numerical implementation.

1. Introduction

Diffuse optical tomography involves the reconstruction of the spatially varying optical properties of a turbid medium. It is usually formulated as inverse problem with respect to the forward problem describing photon propagation in the tissue for given optical coefficients [1].

The forward model is described by the photon diffusion equation with the Robin boundary condition. In the frequency domain, it is given by −∇⋅(𝜅∇Φ)+𝜇_𝑎+𝑖𝜔𝑐Φ=𝑞inΩ,Φ+2𝑎𝜈⋅(𝜅∇Φ)=0on𝜕Ω,(1.1) where Ω is a Lipschitz domain in ℝ^𝑛, 𝑛=2,3,…, 𝜕Ω is its boundary, 𝜈 is the unit outward normal vector on the boundary, Φ is the photon density, 𝑞 is a source term, 𝑎 is a refraction parameter, and 𝜇_𝑎, 𝜇^_𝑠, and 𝜅=1/3(𝜇_𝑎+𝜇^_𝑠) are the absorption, reduced scattering, and diffusion coefficients, respectively. Assume that 𝑎 is a constant and 𝜅, 𝜇_𝑎, 𝜇^_𝑠 are scalar functions satisfying0<𝐿≤𝜅,𝜇_𝑎,𝜇^_𝑠,𝑎≤𝑈(1.2) for positive constants 𝐿 and 𝑈. The unique determination of the optical coefficients is studied in electrical impedance tomography problem [2–5] and some elliptic problem [6], which is applicable to diffuse optical tomography problem also. Let us denote 𝑥=(𝜇_𝑎,𝜅) and Φ=Φ(𝑥) to emphasize the dependence of Φ on the optical coefficient 𝑥.

Assuming we know some a priori information 𝑥₀ about the structural optical coefficients 𝑥 and the perturbation of the optical coefficients 𝛿𝑥=𝑥−𝑥₀, the diffuse optical tomography problem is to find the perturbation of the optical coefficients 𝛿𝑥 from the difference Φ(𝑥+𝛿𝑥)−Φ(𝑥) between the perturbed and unperturbed photon density distribution on the boundary 𝜕Ω. The relation between 𝛿𝑥 and Φ(𝑥+𝛿𝑥)−Φ(𝑥) is given by the following Born expansion [7, 8]: Φ(𝑥+𝛿𝑥)−Φ(𝑥)=ℛ¹(𝑥,𝛿𝑥)+ℛ²𝑥,(𝛿𝑥)²+⋯,(1.3) where ℛ⁰(𝑥)=Φ(𝑥),ℛ^𝑖(𝛿𝑥)^𝑖=ℛ𝛿𝑥,ℛ^𝑖−1(𝛿𝑥)^𝑖−1,𝑓,𝑖=1,2,…,ℛ(𝛿𝑥,𝑓)=ℛ_𝜇𝑎(𝛿𝑥,𝑓)+ℛ_𝜅(𝛿𝑥,𝑓),ℛ_𝜇𝑎(𝛿𝑥,𝑓)=_Ω𝛿𝜇_𝑎(𝜂)𝑅(⋅,𝜂)𝑓(𝜂)𝑑𝜂,ℛ_𝜅(𝛿𝑥,𝑓)=_Ω𝛿𝜅(𝜂)∇𝑅(⋅,𝜂)⋅∇𝑓(𝜂)𝑑𝜂,(1.4) and 𝑅(⋅,𝜂) is the Robin function for a source at 𝜂, which is the solution of (1.1) for the optical coefficient 𝑥 when 𝑞 is the Dirac delta function. By definition of (1.4), the operator ℛ and ℛ¹ are different in the following sense: ℛ¹(𝛿𝑥)=ℛ𝛿𝑥,ℛ⁰=ℛ(𝛿𝑥,Φ).(1.5) Let the perturbation of the coefficients be 𝛿𝑥^† when we neglect second-order terms and higher in the Born expansion (1.3). We can then formulate the linearized diffuse optical tomography problem to find 𝛿𝑥^† from the following equation, which is the first-order Born approximation: ℛ¹𝛿𝑥^†=Φ(𝑥+𝛿𝑥)−Φ(𝑥).(1.6) This linearized diffuse optical tomography problem is simple to implement and widely used [9, 10].

In this paper, a new method, which is more accurate than the linearized method, will be suggested (1.6), which is based on the second-order Born approximation. And the method is faster than the full nonlinear method [11]. Let the solution of the proposed method in this paper be 𝛿𝑥^𝐵, and let 𝛿𝑥 be sufficiently small. Then, the error for the linearized solution 𝛿𝑥^† and the proposed solution 𝛿𝑥^𝐵 is given by‖‖𝛿𝑥^†−𝛿𝑥‖‖_𝒜≤𝐶^†‖𝛿𝑥‖²_𝒜,(1.7a)‖‖𝛿𝑥^𝐵−𝛿𝑥‖‖_𝒜≤𝐶^𝐵‖𝛿𝑥‖³_𝒜,(1.7b) where 𝒜=𝐿^∞(Ω)×𝐿^∞(Ω) and 𝐶^† and 𝐶^𝐵 are constants which are independent of 𝛿𝑥. Hence, the error of the proposed solution 𝑥^𝐵 in (1.7b) is of the order 𝑂(‖𝛿𝑥‖³_𝒜), which is higher than the order of the error of the linearized solution 𝑥^†, 𝑂(‖𝛿𝑥‖²_𝒜).

The detailed statement with proof will be proved in Section 2. Numerical algorithm involving the detailed computation of the second-order term is given in Section 3. Numerical implementation of the proposed method and the linearized method is given in Section 4, and the conclusion of the paper is given in Section 5.

2. Error Analysis

Instead of solving linearized solution 𝛿𝑥^† in (1.6), we suggest the second order solution 𝛿𝑥^𝐵 satisfying ℛ¹𝛿𝑥^𝐵=(Φ(𝑥+𝛿𝑥)−Φ(𝑥))−ℛ²𝛿𝑥^†²,(2.1) or equivalently, ℛ¹𝛿𝑥^𝐵−𝛿𝑥^†=−ℛ²𝛿𝑥^†².(2.2) In this section, we analyze the error for the linearized solution 𝛿𝑥^† and the suggested solution 𝛿𝑥^𝐵.

Let ℬ=𝐻¹(Ω); then, the operator ℛ and ℛ^𝑖, 𝑖=1,2,…, are considered to be the operators from 𝒜×ℬ→ℬ and 𝒜^𝑖(=𝒜×⋯×𝒜^{𝑖times})→ℬ, respectively, by the definition given in (1.4). For the detailed explanation about the definitions of higher-order Fréchet derivative in diffuse optical tomography and its relation to the Born expansion, see [7].

Proposition 2.1. Let Φ be the solution of (1.1) for the given optical coefficients 𝜇_𝑎, 𝜅, source 𝑞, and modulating frequency 𝜔. Then one gets the following relation between the operators between ℛ and ℛ^𝑖, 𝑖=1,2,…: ‖‖ℛ^𝑖‖‖_{𝒜𝑖→𝐵}≤‖ℛ‖^𝑖_{𝒜×ℬ→𝐵}‖Φ‖_ℬ(2.3) for 𝑖=1,2,….

Proof. By the induction argument on 𝑖=1,2,… and using (1.5), we get the following inequality: ‖‖ℛ¹‖‖_𝒜→ℬ≤‖ℛ‖_{𝒜×ℬ→ℬ}‖Φ‖_ℬ,(2.4) which is (2.3) for 𝑖=1. Suppose that (2.3) holds for 𝑖=1,2,…,𝐼−1. Then we obtain ‖‖ℛ^𝐼(𝛿𝑥)^𝐼‖‖_ℬ≤‖ℛ‖_{𝒜×ℬ→ℬ}‖𝛿𝑥‖_𝒜‖‖ℛ^𝐼−1(𝛿𝑥)^𝐼−1‖‖_ℬ≤‖ℛ‖_{𝒜×ℬ→ℬ}‖𝛿𝑥‖_𝒜‖‖ℛ^𝐼−1‖‖_𝒜𝐼−1‖𝛿𝑥‖^𝐼−1_𝒜≤‖ℛ‖^𝐼_{𝒜×ℬ→ℬ}‖𝛿𝑥‖^𝐼_𝒜‖Φ‖_ℬ.(2.5) Using (2.5) and the definition of the operator norm ‖⋅‖_{𝒜𝐼→ℬ}, we obtain (2.3) for 𝑖=𝐼. Therefore, by the induction argument, we have proved (2.3) for 𝑖=1,2,….

By [7], ‖ℛ‖_{𝒜×ℬ→ℬ} is bounded, and thus ‖ℛ^𝑖‖_{𝒜𝑖→ℬ}, 𝑖=1,2,…, are also bounded by Proposition 2.1. Let us assume that there exists a bounded operator (ℛ¹)^† from ℬ to 𝒜 such that (ℛ¹)^†(ℛ¹)=𝑖𝑑_𝒜. (ℛ¹)^† is usually called the left inverse of ℛ¹. Let us denote ‖𝛿𝑥‖∶=‖𝛿𝑥‖_𝒜,‖Φ(𝑥)‖∶=‖Φ(𝑥)‖_ℬ,‖ℛ‖∶=‖ℛ‖_{𝒜×ℬ→ℬ},‖‖ℛ^𝑖‖‖∶=‖‖ℛ^𝑖‖‖_{𝒜𝑖→ℬ},𝑖=1,2,…,‖‖ℛ¹^†‖‖∶=‖‖ℛ¹^†‖‖_ℬ→𝒜,(2.6) for brevity.

Using Proposition 2.1 and the assumption on the left inverse, the main theorem of this paper is given as follows.

Theorem 2.2. Assume that there exists (ℛ¹)† such that (𝑅¹)^†ℛ¹=𝑖𝑑 and ‖(ℛ¹)^†‖ is bounded, and let‖𝛿𝑥‖≤12‖ℛ‖.(2.7a) Then, ‖‖𝛿𝑥^†−𝛿𝑥‖‖≤𝐶^†‖𝛿𝑥‖²,(2.7b)‖‖𝛿𝑥^𝐵−𝛿𝑥‖‖≤𝐶^𝐵‖𝛿𝑥‖³,(2.7c) where 𝐶^†∶=2‖‖ℛ¹^†‖‖‖ℛ‖²‖Φ‖,𝐶^𝐵∶=𝐶^†³4‖ℛ‖+𝐶^†²+‖ℛ‖𝐶^†=𝐶^†‖ℛ‖𝐶^†2‖ℛ‖+1².(2.8)

Proof. By (1.3) and (1.6), we obtain ℛ¹𝛿𝑥^†−𝛿𝑥=ℛ²(𝛿𝑥)²+ℛ³(𝛿𝑥)³+⋯.(2.9) Therefore we arrive at (2.7b) by the following inequality: ‖‖𝛿𝑥^†−𝛿𝑥‖‖≤‖‖ℛ¹^†‖‖(‖ℛ‖‖𝛿𝑥‖)²‖Φ‖1−(‖ℛ‖‖𝛿𝑥‖)≤𝐶^†‖𝛿𝑥‖².(2.10) From (2.7a) and (2.7b), we obtain the following upper bound of ‖𝛿𝑥^†‖: ‖‖𝛿𝑥^†‖‖≤1+𝐶^†‖𝛿𝑥‖‖𝛿𝑥‖≤1+‖‖ℛ¹⁺‖‖‖ℛ‖‖Φ‖‖𝛿𝑥‖.(2.11) Using (2.2) and (2.9), we obtain ℛ¹𝛿𝑥−𝛿𝑥^𝐵=ℛ²(𝛿𝑥)²−ℛ²𝛿𝑥^†²+ℛ³(𝛿𝑥)³+ℛ⁴(𝛿𝑥)⁴+⋯.(2.12) The second-order term on the righthand side of (2.12) is analyzed as follows: ℛ²(𝛿𝑥)²−ℛ²𝛿𝑥^†²=ℛ(𝛿𝑥,ℛ(𝛿𝑥,Φ))−ℛ𝛿𝑥^†,ℛ𝛿𝑥^†,Φ=ℛ𝛿𝑥,ℛ𝛿𝑥−𝛿𝑥^†,Φ+ℛ𝛿𝑥−𝛿𝑥^†,ℛ𝛿𝑥^†,Φ.(2.13) From (2.12), we obtain ‖‖𝛿𝑥−𝛿𝑥^𝐵‖‖≤‖‖ℛ¹^†‖‖‖‖ℛ²(𝛿𝑥)²−ℛ²𝛿𝑥^†²‖‖+‖‖ℛ³(𝛿𝑥)³+ℛ⁴(𝛿𝑥)⁴+⋯‖‖.(2.14) By using (2.3), (2.10), (2.11), (2.13), and the definition of 𝐶^†, (2.7c) is achieved from (2.14) as follows: ‖‖𝛿𝑥−𝛿𝑥^𝐵‖‖≤‖‖ℛ¹^†‖‖‖Φ‖‖ℛ‖²‖‖𝛿𝑥−𝛿𝑥^†‖‖‖𝛿𝑥‖+‖‖𝛿𝑥^†‖‖+(‖ℛ‖‖𝛿𝑥‖)³1−‖ℛ‖‖𝛿𝑥‖≤‖‖ℛ¹^†‖‖‖Φ‖‖ℛ‖²‖𝛿𝑥‖³𝐶^†2+‖‖ℛ¹^†‖‖‖ℛ‖‖Φ‖+2‖ℛ‖≤𝐶^†‖𝛿𝑥‖³𝐶^†1+𝐶^†4‖ℛ‖+‖ℛ‖≤𝐶^𝐵‖𝛿𝑥‖³.(2.15)

3. Numerical Algorithm

Assume that we can measure the photon density distribution Φ(𝑥+𝛿𝑥) and Φ(𝑥) on the entire boundary 𝜕Ω. That is to say, we have infinite detectors and one source. Then, the numerical algorithm is given as follows.

The detailed computation of the integral operators ℛ¹ and ℛ², which is introduced in (1.5), is as follows: ℛ¹(𝛿𝑥)=ℛ_𝜇𝑎𝛿𝜇_𝑎,Φ+ℛ_𝜅(𝛿𝜅,Φ),(3.1a)ℛ²(𝛿𝑥)=ℛ_𝜇𝑎𝛿𝜇_𝑎,ℛ_𝜇𝑎𝛿𝜇_𝑎,Φ+ℛ_𝜇𝑎𝛿𝜇_𝑎,ℛ_𝜅(𝛿𝜅,Φ),+ℛ_𝜅𝛿𝜅,ℛ_𝜇𝑎𝛿𝜇_𝑎,Φ+ℛ_𝜅𝛿𝜅,ℛ_𝜅(𝛿𝜅,Φ).(3.1b)

3.1. Discretization

Algorithm 1 is based on one source and infinite detectors. However, for practical reasons, we need to discretize Algorithm 1 to obtain the numerical algorithm for finite sources and finite detectors for finite frequencies. The following notations will be used for the discretization: (i)𝑁_𝑑 detector positions: 𝑟_𝑖𝑑 for 𝑖_𝑑=1,2,…,𝑁_𝑑, (ii)𝑁_𝑠 source functions: 𝑞_𝑖𝑠=𝛿_𝑖𝑠(Dirac delta function) for 𝑖_𝑠=1,2,⋯,𝑁_𝑠, (iii)𝑁_𝜔 frequencies: 𝜔_𝑖𝜔 for 𝑖_𝜔=1,2,…,𝑁_𝜔,(iv)𝑁_𝑒 elements: 𝑇_𝑖𝑒 for 𝑖_𝑒=1,2,…,𝑁_𝑒,(v)𝑁_𝑛 nodes: 𝑡_𝑖𝑛 for 𝑖_𝑛=1,2,…,𝑁_𝑛,(vi) the measurement index: 𝑗=(𝑖_𝜔−1)𝑁_𝑠𝑁_𝑑+(𝑖_𝑠−1)𝑁_𝑑+𝑖_𝑑, (vii) the optical coefficient index: 𝑘=(𝑖_𝜇𝜅−1)𝑁_𝑒+𝑖_𝑒, where 𝑖_𝜇𝜅 is 1 (the absorption coefficient) or 2 (the diffusion coefficient).

(I) Compute the solution Φ(𝑥) and the Robin function 𝑅(𝑥) and its first and second derivatives. (II) Find 𝛿𝑥† by solving ℛ1(𝛿𝑥†)=Φ(𝑥+𝛿𝑥)−Φ(𝑥) as in (1.6). (III) Find 𝛿𝑥Δ=𝛿𝑥𝐵−𝛿𝑥† by solving ℛ1(𝛿𝑥Δ)=−ℛ2(𝛿𝑥†) as in (2.2). (IV) Compute 𝛿𝑥𝐵 by adding 𝛿𝑥† and 𝛿𝑥Δ.

Algorithm 1 

Numerical algorithm (continuous version).

If we use piecewise linear or bilinear finite element method, the finite element solution is represented by𝑢_ℎ(𝑥)=^𝑁𝑛_𝑖𝑛=1𝑢_ℎ𝑖_𝑛𝜙_𝑖𝑛(𝑥),(3.2) where 𝜙_𝑖𝑛 is the piecewise linear or the bilinear function which is 1 on the 𝑖_𝑛th node and 0 on all the other nodes. Assume 𝜇_𝑎 and 𝜅 are piecewise constant function, which is constant for each 𝑁_𝑒 finite elements. Therefore, in diffuse optical tomography inverse problem, we have 𝑁_𝜔𝑁_𝑠𝑁_𝑑 measurement information contents and 2𝑁_𝑒 variables to find.

We should discretize ℛ¹ and ℛ² to obtain a discretized version of Algorithm 1. Let the Jacobian and Hessian matrices, which is the discretization of integral operators ℛ¹ and ℛ², be 𝐽 and 𝐻. The relation between higher order derivatives for the diffusion operator and higher order terms of Born expansions including ℛ¹ and ℛ² is analyzed in [7].

Firstly, let us discretize 𝛿𝑥, Φ, and the Robin function 𝑅 as follows: 𝛿𝑥≈⎛⎜⎝^𝑁𝑒_𝑖_𝑒=1𝛿𝜇_𝑖𝑒𝜒_𝑇𝑖𝑒,^𝑁𝑒_𝑖𝑒=1𝛿𝜅_𝑖𝑒𝜒_𝑇𝑖𝑒⎞⎟⎠,(3.3a)Φ^{𝑖𝜔,𝑖𝑠}≈^𝑁𝑛_𝑖𝑛=1Φ^{𝑖𝜔,𝑖𝑠}_𝑖𝑛𝜙_𝑖𝑛,(3.3b)𝑅^𝑖𝜔⋅,𝑟_𝑖𝑠≈^𝑁𝑛_𝑖𝑛=1𝑅^{𝑖𝜔,𝑖𝑠}_𝑖𝑛𝜙_𝑖𝑛.(3.3c) Since we chose the source function 𝑞_𝑠 as the Dirac delta function at the 𝑖_𝑠th source point, Φ^{𝑖𝜔,𝑖𝑠}=𝑅^𝑖𝜔(⋅,𝑟_𝑖𝑠). However, we will discriminate these two functions in this paper, since they are different for general source function 𝑞 which is different from the Dirac delta function. We will use 𝛿𝜇 instead of 𝛿𝜇_𝑎 for notational convenience.

Let the vector 𝛾₀ which corresponds to the discretization of 𝛿𝑥 in (3.3a) be defined as 𝛾₀=𝛿𝜇₁,𝛿𝜇₂,…,𝛿𝜇_𝑁𝑒,𝛿𝜅₁,𝛿𝜅₂,…,𝛿𝜅_𝑁𝑒.(3.4) By the adjoint method [12], 𝑅^𝑖𝜔(𝑟_𝑖𝑑,⋅)=(𝑅^𝑖𝜔(⋅,𝑟_𝑖𝑑))^∗, where ∗ denotes complex conjugate. Likewise for (3.3a), let 𝛾, 𝛾^†, 𝛾^Δ, and 𝛾^𝐵 be the discretization of 𝛿𝑥, 𝛿𝑥^†, 𝛿𝑥^Δ, and 𝛿𝑥^𝐵, respectively.

For a function 𝑓 and a measurable set 𝑇, let us denote 𝑓∈𝑇 if the intersection of the support of 𝑓 and 𝑇 is not empty. The discretization of the linearized solution 𝛾^† is attained by solving the following equation: 𝐽𝛾^†=𝑏,(3.5) where 𝐽(𝑗,𝑘)=_{𝜙𝑖𝑛1∈𝑇𝑖𝑒}_{𝜙𝑖𝑛2∈𝑇𝑖𝑒}𝑅^{𝑖𝜔,𝑖𝑑}_𝑖𝑛1^∗𝐸_𝑖𝑒𝑖_𝑛1,𝑖_𝑛2Φ^{𝑖𝜔,𝑖𝑠}_𝑖𝑛2when𝑖_𝜇𝜅=1,𝐽(𝑗,𝑘)=_{𝜙𝑖𝑛1∈𝑇𝑖𝑒}_{𝜙𝑖𝑛2∈𝑇𝑖𝑒}−3𝑅^{𝑖𝜔,𝑖𝑑}_𝑖𝑛1^∗𝜅²_𝑖𝑒𝐹_𝑖𝑒𝑖_𝑛1,𝑖_𝑛2Φ^{𝑖𝜔,𝑖𝑠}_𝑖𝑛2when𝑖_𝜇𝜅=2,𝑏(𝑗)=Φ^{𝑖𝜔,𝑖𝑠}(𝑥+𝛿𝑥)𝑟_𝑖𝑑−Φ^{𝑖𝜔,𝑖𝑠}(𝑥)𝑟_𝑖𝑑,𝐸_𝑖𝑒𝑖_𝑛1,𝑖_𝑛2=_𝑇𝑖𝑒𝜙_𝑖𝑛1(𝜉)𝜙_𝑖𝑛2(𝜉)𝑑𝜉,𝐹_𝑖𝑒𝑖_𝑛1,𝑖_𝑛2=_𝑇𝑖𝑒∇𝜙_𝑖𝑛1(𝜉)⋅∇𝜙_𝑖𝑛2(𝜉)𝑑𝜉.(3.6)

The discretized solution 𝛾^Δ is obtained by solving the following equation: 𝐽𝛾^Δ=𝛾^†^𝑡𝐻𝛾^†,(3.7) where 𝐻𝑗,𝑖_𝑒1,𝑖_𝑒2=_{𝜙𝑖𝑛1∈𝑇𝑖𝑒1}_{𝜙𝑖𝑛2∈𝑇𝑖𝑒2}𝑅^{𝑖𝜔,𝑖𝑑}_𝑖𝑛1^∗𝐻_𝜇𝜇+𝐻_𝜇𝜅+𝐻_𝜅𝜇+𝐻_𝜅𝜅𝑖_𝑒1,𝑖_𝑒2;𝑖_𝑛1,𝑖_𝑛2Φ^{𝑖𝜔,𝑖𝑠}_𝑖𝑛2,(3.8) where 𝐻_𝜇𝜇, 𝐻_𝜇𝜅, 𝐻_𝜅𝜇, and 𝐻_𝜅𝜅 are the discretization of corresponding terms in (3.1b) such that 𝐻_𝜇𝜇𝑖_𝑒1,𝑖_𝑒2;𝑖_𝑛1,𝑖_𝑛2=_{𝑇𝑖𝑒1}_{𝑇𝑖𝑒2}𝜙_𝑖𝑛1(𝜉)𝑅^𝑖𝜔(𝜉,𝜂)𝜙_𝑖𝑛2(𝜂)𝑑𝜉𝑑𝜂,𝐻_𝜇𝜅𝑖_𝑒1,𝑖_𝑒2;𝑖_𝑛1,𝑖_𝑛2=_{𝑇𝑖𝑒1}_{𝑇𝑖𝑒2}𝜙_𝑖𝑛1(𝜉)∇_𝜂𝑅^𝑖𝜔(𝜉,𝜂)⋅∇_𝜂𝜙_𝑖𝑛2(𝜂)𝑑𝜉𝑑𝜂,𝐻_𝜅𝜇𝑖_𝑒1,𝑖_𝑒2;𝑖_𝑛1,𝑖_𝑛2=∫_{𝑇𝑖𝑒1}∫_{𝑇𝑖𝑒2}∇_𝜉𝜙_𝑖𝑛1(𝜉)⋅∇_𝜉𝑅^𝑖𝜔(𝜉,𝜂)𝜙_𝑖𝑛2(𝜂)𝑑𝜉𝑑𝜂,𝐻_𝜅𝜅𝑖_𝑒1,𝑖_𝑒2;𝑖_𝑛1,𝑖_𝑛2=_{𝑇𝑖𝑒1}_{𝑇𝑖𝑒2}∇_𝜉𝜙_𝑖𝑛1(𝜉)⋅∇_𝜉∇_𝜂𝑅^𝑖𝜔(𝜉,𝜂)∇_𝜂𝜙_𝑖𝑛2(𝜂)𝑑𝜉𝑑𝜂.(3.9) Even though the Hessian is not discretized, we obtain the following discretized numerical algorithm (Algorithm 2), expecting the Hessian is simply discretized and approximated in the next subsection:

(I) Compute the solution Φ(𝛾0)𝑖𝜔,𝑖𝑠𝑖𝑛 and the Robin function 𝑅(𝛾0)𝑖𝜔𝑖𝑑,𝑖𝑛 for 𝑖𝜔=1,⋯,𝑁𝜔,𝑖𝑠=1,…,𝑁𝑠,𝑖𝑛=1,…,𝑁𝑛 as in (3.3b) and (3.3c), respectively. (II) Find 𝛾† by solving the equation (3.5). (III) Find 𝛾Δ by solving the equation (3.7). (IV) Compute 𝛾𝐵 by adding 𝛾† and 𝛾Δ.

Algorithm 2 

Numerical algorithm (discretized version).

3.2. Approximation of Hessian

In this subsection we approximate 𝐻_𝜇𝜇, by assuming 𝜅 and 𝜇^_𝑠 are constant in Ω. The approximation is progressed in three ways.

First, we approximate the Robin function 𝑅(𝜉,𝜂) when (𝜉,𝜂)∈Ω⧵𝜕Ω by its leading term 𝑅₀(𝜉,𝜂) defined by 𝑅₀(𝜉,𝜂)=⎧⎪⎪⎨⎪⎪⎩1(𝑝−2)𝑔_𝑝𝜅(𝜂)||𝜉−𝜂||^2−𝑝𝑝≥3,1𝜔₂𝜅(𝜂)log2𝑆||𝜉−𝜂||𝑝=2,(3.10) where 𝑔_𝑝 is the hypersurface area of the unit sphere in ℝ^𝑝, 𝑝=2,3,… and 𝑆=sup_{𝜉,𝜂∈Ω}|𝜉−𝜂|. Some important relations between 𝑅 and 𝑅₀ are found in [13].

Second, when 𝑖_𝑒1≠𝑖_𝑒2, the Robin function 𝑅 and 𝜙_𝑖𝑛 are approximated by constant values 𝑅₀(𝑐(𝑖_𝑒1),𝑐(𝑖_𝑒2)) and 𝜙_𝑖𝑛(𝑐(𝑖_𝑒)) in 𝑇_𝑖𝑒, respectively, where 𝑐(𝑖_𝑒) of the center of the element 𝑇_𝑖𝑒. That is to say, when 𝑖_𝑒1≠𝑖_𝑒2, (3.9) is approximated as follows: 𝐻_𝜇𝜇𝑖_𝑒1,𝑖_𝑒2;𝑖_𝑛1,𝑖_𝑛2=𝑅₀𝑐𝑖_𝑒1,𝑐𝑖_𝑒2_{𝑇𝑖𝑒1}𝜙_𝑖𝑛1(𝜉)𝑑𝜉_{𝑇𝑖𝑒2}𝜙_𝑖𝑛2(𝜂)𝑑𝜂.(3.11)

Third, when 𝑖_𝑒1=𝑖_𝑒2, we use the following lemma.

Lemma 3.1. Let the measurable set 𝑇 be contained in ℝ^𝑝, 𝑝=2,3,…, and 0<𝑚<𝑝; then, the following inequality holds for 𝑇: _𝑇||𝜉−𝜂||^−𝑚𝑑𝜉𝑑𝜂≤𝑝^{1−𝑚/𝑝}𝑝−𝑚𝑔^𝑚/𝑝_𝑝||𝑇||^{2−𝑚/𝑝},𝑝≥2,(3.12a)_𝑇log2𝑆||𝜉−𝜂||𝑑𝜉𝑑𝜂≤14𝜋1+log4𝑆²𝜋||𝑇||||𝑇||²,𝑝=2,(3.12b) where |𝑇| is the volume of 𝑇.

Proof. If a ball with a radius 𝑟 has the same volume as 𝑇, we have 𝑟=|𝑇|𝑝𝑔_𝑝^1/𝑝(3.13) for the space dimensions 𝑝=2,3,…. Let the ball of radius 𝑟 with center 𝜉∈𝑇 be 𝐵_𝜉. Let 𝑇₀=𝑇∩𝐵_𝜉, 𝑇⁺=𝑇⧵𝐵_𝜉, and 𝑇⁻=𝐵_𝜉⧵𝑇. Noting that |𝑇⁺|=|𝑇⁻|, we obtain _𝑇||𝜉−𝜂||^−𝑚𝑑𝜂=_𝑇0||𝜉−𝜂||^−𝑚𝑑𝜂+_𝑇+||𝜉−𝜂||^−𝑚𝑑𝜂≤_𝑇0||𝜉−𝜂||^−𝑚𝑑𝜂+_𝑇−||𝜉−𝜂||^−𝑚𝑑𝜂=_𝐵𝜉||𝜉−𝜂||^−𝑚𝑑𝜂≤^𝑟₀𝜌^{𝑝−𝑚−1}𝑔_𝑝𝑑𝜌=𝑔_𝑝𝑝−𝑚𝑟^𝑝−𝑚≤𝑔_𝑝𝑝−𝑚||𝑇||𝑝𝑔_𝑝^{1−𝑚/𝑝}(3.14) for all 𝜉∈𝑇. Therefore, _𝑇||𝜉−𝜂||^−𝑚𝑑𝜂𝑑𝜉≤𝑔_𝑝||𝑇||𝑝−𝑚||𝑇||𝑝𝑔_𝑝^{1−𝑚/𝑝}=𝑝^{1−𝑚/𝑝}𝑝−𝑚𝑔^𝑚/𝑝_𝑝||𝑇||^{2−𝑚/𝑝}.(3.15) Equation (3.12b) is derived in the same manner.

Therefore, when 𝑖_𝑒1=𝑖_𝑒2, (3.9) is approximated using the inequality in Lemma 3.1 as follows: 𝐻_𝜇𝜇𝑖_𝑒1,𝑖_𝑒1;𝑖_𝑛1,𝑖_𝑛2≈𝜙_𝑖𝑛1𝑐_𝑖𝑒1𝜙_𝑖𝑛2𝑐_𝑖𝑒1⋅⎧⎪⎪⎨⎪⎪⎩𝑝^2/𝑝||𝑇_𝑖𝑒1||^1+2/𝑝2(𝑝−2)𝑔^2/𝑝_𝑝𝜅𝑐𝑖_𝑒1𝑝≥3,18𝜋²𝜅𝑐𝑖_𝑒1⎛⎜⎜⎝1+log⎛⎜⎜⎝4𝑆²𝜋||𝑇_𝑖𝑒1||⎞⎟⎟⎠⎞⎟⎟⎠||𝑇_𝑖𝑒1||²𝑝=2.(3.16)

4. Numerical Implementation

In the numerical implementation, the following parameters are used: (i)Ω=[0,6]×[0,6](cm²),(ii)𝑁_𝑑=16,(iii)𝑁_𝑠=16,(iv)𝑁_𝜔=1,(v)𝑁𝑥=𝑁𝑦=16,(vi)𝑁_𝑒=𝑁𝑥∗𝑁𝑦,(vii)𝑁_𝑛=(𝑁𝑥+1)∗(𝑁𝑦+1),(viii)𝜇_𝑎=0.05+(0.2−0.05)𝜒_𝐷(cm⁻¹),(ix)𝜇^_𝑠=8(cm⁻¹),(x)𝜅=1/3∗(𝜇_𝑎+𝜇^_𝑠)=1/3∗(0.05+8),(xi)𝜔=2𝜋∗300 MHz,(xii)𝑎=1,(xiii)𝐽𝑖𝑛𝑑𝑒𝑥=𝑁𝑥∗𝑁𝑦∗0.4.

Since the diffusion coefficient 𝜅 is constant, the right-hand side 𝑏 is a 𝑁_𝑠∗𝑁_𝑑 column vector, Jacobian 𝐽 is a (𝑁_𝑠∗𝑁_𝑑)×𝑁_𝑒 matrix, the Hessian 𝐻 is a 𝑁_𝑒×(𝑁𝑠∗𝑁𝑑)×𝑁_𝑒 third-order tensor, and (𝛾^†)^𝑡𝐻𝛾^† is 𝑁𝑠∗𝑁𝑑 column vector in (3.5) and (3.7). 𝐻=𝐻_𝜇𝜇 is approximated by (3.11) and (3.16).

In the above setting, we reconstruct the obstacle 𝐷 which has different absorption coefficient (0.2 cm⁻¹) compared to the background absorption coefficient (0.05 cm⁻¹). Four cases of the obstacle 𝐷 are considered in Figures 1, 2, 3, and 4. The reconstruction of the absorption coefficient 𝜇_𝑎=0.05+(0.2−0.05)𝜒_𝐷 (cm⁻¹) is implemented using two algorithms. One is the suggested Algorithm 2 based on the second-order Born approximation. The other is linearized method based on the first-order Born approximation, which is equivalent to the step I and II in Algorithm 2. We denoted these two methods in the figures: the 2nd order approximation and the 1st-order approximation, respectively. On the upper-left part of the figures, original 𝜇_𝑎 and source/detector locations are plotted. The initial guess (𝜇_𝑎0 or 𝛾₀) for the absorption coefficient is plotted on the upper-right part of the figures. In the lowerleft and lowerright part of each figure, reconstructed absorption coefficients by the first approximation (𝜇^†_𝑎 or 𝛾^†) and the second approximation (𝜇^𝐵_𝑎 or 𝛾^𝐵) are plotted, respectively.

Figure 1 

𝐽𝑖𝑛𝑑𝑒𝑥=𝑁𝑥∗𝑁𝑦∗0.4, 𝐽𝑎𝑙𝑝ℎ𝑎=6.4920𝑒−009, 10% noise, sources (*), and detectors (o).

Figure 2 

𝐽𝑖𝑛𝑑𝑒𝑥=𝑁𝑥∗𝑁𝑦∗0.4, 𝐽𝑎𝑙𝑝ℎ𝑎=6.5711𝑒−009, 10% noise, sources (*), and detectors (o).

Figure 3 

𝐽𝑖𝑛𝑑𝑒𝑥=𝑁𝑥∗𝑁𝑦∗0.3, 𝐽𝑎𝑙𝑝ℎ𝑎=1.8227𝑒−7, 10% noise, sources (*) and detectors (o).

Figure 4 

𝐽𝑖𝑛𝑑𝑒𝑥=𝑁𝑥∗𝑁𝑦∗0.4, 𝐽𝑎𝑙𝑝ℎ𝑎=6.5029𝑒−009, 10% noise, sources (*), and detectors (o).

In all four cases, 10% noise is added. Truncated singular value decomposition(SVD) is used. 𝐽𝑖𝑛𝑑𝑒𝑥 is the number of largest singular values used in the truncated SVD method. We used the Tikhonov regularization parameter Jalpha as the value of the 𝐽𝑖𝑛𝑑𝑒𝑥th largest singular values.

As is shown in the figure, the discrimination between background and the obstacle is clearer in the second-order approximation than the first-order approximation. The reconstructed image resolution depends on the distance from the boundary of the tissue, which is verified by comparing Figures 1 and 2 with Figures 3 and 4. And the resolution also depends on the size of obstacle, which is verified by comparing Figures 1 and 3 with Figures 2 and 4. Due to the diffusion property of near infrared light, the reconstructed image is much blurred especially in Figure 3. The sensitivity to the noise made some kind of irregular checkerboard pattern near the boundary (Figures 1, 3, and 4).

5. Conclusions

We derived a new numerical method based on the second-order Born approximation. The method is a method of order 3, which is more accurate than the well-known linearized method based on the first-order Born approximation. The error analysis for the method is proved, and the computation of the second-order term is explained using some approximation and integral inequalities. The comparison between the suggested and the linearized method is implemented for four different kinds of absorption coefficients. In the implementation, the suggested method shows more discrimination between the optical obstacle and the background than the linearized method. If more accurate numerical quadrature with more efficient approximation of the Robin function is used, the efficiency of the present method will be elaborated. The simultaneous reconstruction of the absorption and the reduced scattering coefficients based on the proper approximation on the second derivatives of the Robin function would be an interesting topic.

Acknowledgment

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2010-0004047).