# The Second-Order Born Approximation in Diffuse Optical Tomography

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

---

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


---

## Body

## 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(1.1)-∇⋅(κ∇Φ)+(μa+iωc)Φ=qinΩ,Φ+2aν⋅(κ∇Φ)=0on⁡∂Ω,
where Ω is a Lipschitz domain in ℝn, n=2,3,…, ∂Ω is its boundary, ν is the unit outward normal vector on the boundary, Φ is the photon density, q is a source term, a is a refraction parameter, and μa, μs′, and κ=1/3(μa+μs′) are the absorption, reduced scattering, and diffusion coefficients, respectively. Assume that a is a constant and κ, μa, μs′ are scalar functions satisfying(1.2)0<L≤κ,μa,μs′,a≤U
for positive constants L and U. 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 x=(μa,κ) and Φ=Φ(x) to emphasize the dependence of Φ on the optical coefficient x.Assuming we know some a priori informationx0 about the structural optical coefficients x and the perturbation of the optical coefficients δx=x-x0, the diffuse optical tomography problem is to find the perturbation of the optical coefficients δx from the difference Φ(x+δx)-Φ(x) between the perturbed and unperturbed photon density distribution on the boundary ∂Ω. The relation between δx and Φ(x+δx)-Φ(x) is given by the following Born expansion [7, 8]: (1.3)Φ(x+δx)-Φ(x)=R1(x,δx)+R2(x,(δx)2)+⋯,
where (1.4)R0(x)=Φ(x),Ri((δx)i)=R(δx,Ri-1((δx)i-1,f)),i=1,2,…,R(δx,f)=Rμa(δx,f)+Rκ(δx,f),Rμa(δx,f)=∫Ωδμa(η)R(⋅,η)f(η)dη,Rκ(δx,f)=∫Ωδκ(η)∇R(⋅,η)⋅∇f(η)dη,
and R(·,η) is the Robin function for a source at η, which is the solution of (1.1) for the optical coefficient x when q is the Dirac delta function. By definition of (1.4), the operator ℛ and ℛ1 are different in the following sense: (1.5)R1(δx)=R(δx,R0)=R(δx,Φ).
Let the perturbation of the coefficients be δx† 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 δx† from the following equation, which is the first-order Born approximation: (1.6)R1(δx†)=Φ(x+δx)-Φ(x).
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 δxB, and let δx be sufficiently small. Then, the error for the linearized solution δx† and the proposed solution δxB is given by(1.7a)‖δx†-δx‖A≤C†‖δx‖A2,(1.7b)‖δxB-δx‖A≤CB‖δx‖A3, where 𝒜=L∞(Ω)×L∞(Ω) and C† and CB are constants which are independent of δx. Hence, the error of the proposed solution xB in (1.7b) is of the order O(∥δx∥𝒜3), which is higher than the order of the error of the linearized solution x†, O(∥δx∥𝒜2).The detailed statement with proof will be proved in Section2. 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δx† in (1.6), we suggest the second order solution δxB satisfying (2.1)R1(δxB)=(Φ(x+δx)-Φ(x))-R2(δx†)2,
or equivalently, (2.2)R1(δxB-δx†)=-R2(δx†)2.
In this section, we analyze the error for the linearized solution δx† and the suggested solution δxB.Letℬ=H1(Ω); then, the operator ℛ and ℛi, i=1,2,…, are considered to be the operators from 𝒜×ℬ→ℬ and 𝒜i(=𝒜×⋯×𝒜︷itimes)→ℬ, 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 μa, κ, source q, and modulating frequency ω. Then one gets the following relation between the operators between ℛ and ℛi, i=1,2,…:
(2.3)‖Ri‖Ai→B≤‖R‖A×B→Bi‖Φ‖B
for i=1,2,….Proof.
By the induction argument oni=1,2,… and using (1.5), we get the following inequality:
(2.4)‖R1‖A→B≤‖R‖A×B→B‖Φ‖B,
which is (2.3) for i=1. Suppose that (2.3) holds for i=1,2,…,I-1. Then we obtain
(2.5)‖RI(δx)I‖B≤‖R‖A×B→B‖δx‖A‖RI-1(δx)I-1‖B≤‖R‖A×B→B‖δx‖A‖RI-1‖AI-1‖δx‖AI-1≤‖R‖A×B→BI‖δx‖AI‖Φ‖B.
Using (2.5) and the definition of the operator norm ∥·∥𝒜I→ℬ, we obtain (2.3) for i=I. Therefore, by the induction argument, we have proved (2.3) for i=1,2,….By [7], ∥ℛ∥𝒜×ℬ→ℬ is bounded, and thus ∥ℛi∥𝒜i→ℬ, i=1,2,…, are also bounded by Proposition 2.1. Let us assume that there exists a bounded operator (ℛ1)† from ℬ to 𝒜 such that (ℛ1)†(ℛ1)=id𝒜. (ℛ1)† is usually called the left inverse of ℛ1. Let us denote (2.6)‖δx‖:=‖δx‖A,‖Φ(x)‖:=‖Φ(x)‖B,‖R‖:=‖R‖A×B→B,‖Ri‖:=‖Ri‖Ai→B,i=1,2,…,‖(R1)†‖:=‖(R1)†‖B→A,
for brevity.Using Proposition2.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(ℛ1)† such that (R1)†ℛ1=id and ∥(ℛ1)†∥ is bounded, and let(2.7a)‖δx‖≤12‖R‖.
Then,
(2.7b)‖δx†-δx‖≤C†‖δx‖2,(2.7c)‖δxB-δx‖≤CB‖δx‖3, where
(2.8)C†:=2‖(R1)†‖‖R‖2‖Φ‖,CB:=(C†)34‖R‖+(C†)2+‖R‖C†=C†‖R‖(C†2‖R‖+1)2.Proof.
By (1.3) and (1.6), we obtain
(2.9)R1(δx†-δx)=R2(δx)2+R3(δx)3+⋯.
Therefore we arrive at (2.7b) by the following inequality:
(2.10)‖δx†-δx‖≤‖(R1)†‖(‖R‖‖δx‖)2‖Φ‖1-(‖R‖‖δx‖)≤C†‖δx‖2.
From (2.7a) and (2.7b), we obtain the following upper bound of ∥δx†∥:
(2.11)‖δx†‖≤(1+C†‖δx‖)‖δx‖≤(1+‖(R1)+‖‖R‖‖Φ‖)‖δx‖.
Using (2.2) and (2.9), we obtain
(2.12)R1(δx-δxB)=R2(δx)2-R2(δx†)2+R3(δx)3+R4(δx)4+⋯.
The second-order term on the righthand side of (2.12) is analyzed as follows:
(2.13)R2(δx)2-R2(δx†)2=R(δx,R(δx,Φ))-R(δx†,R(δx†,Φ))=R(δx,R(δx-δx†,Φ))+R(δx-δx†,R(δx†,Φ)).
From (2.12), we obtain
(2.14)‖δx-δxB‖≤‖(R1)†‖[‖R2(δx)2-R2(δx†)2‖+‖R3(δx)3+R4(δx)4+⋯‖].
By using (2.3), (2.10), (2.11), (2.13), and the definition of C†, (2.7c) is achieved from (2.14) as follows:
(2.15)‖δx-δxB‖≤‖(R1)†‖‖Φ‖[‖R‖2‖δx-δx†‖(‖δx‖+‖δx†‖)+(‖R‖‖δx‖)31-‖R‖‖δx‖]≤‖(R1)†‖‖Φ‖‖R‖2‖δx‖3[C†(2+‖(R1)†‖‖R‖‖Φ‖)+2‖R‖]≤C†‖δx‖3[C†(1+C†4‖R‖)+‖R‖]≤CB‖δx‖3.


## 3. Numerical Algorithm

Assume that we can measure the photon density distributionΦ(x+δx) and Φ(x) 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ℛ1 and ℛ2, which is introduced in (1.5), is as follows: (3.1a)R1(δx)=Rμa(δμa,Φ)+Rκ(δκ,Φ),(3.1b)R2(δx)=Rμa(δμa,Rμa(δμa,Φ))+Rμa(δμa,Rκ(δκ,Φ)),+Rκ(δκ,Rμa(δμa,Φ))+Rκ(δκ,Rκ(δκ,Φ)).


### 3.1. Discretization

Algorithm1 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)
Nd detector positions: rid for id=1,2,…,Nd,(ii)
Ns source functions: qis=δis(Dirac delta function) for is=1,2,⋯,Ns,(iii)
Nω frequencies: ωiω for iω=1,2,…,Nω,(iv)
Ne elements: Tie for ie=1,2,…,Ne,(v)
Nn nodes: tin for in=1,2,…,Nn,(vi)
the measurement index:j=(iω-1)NsNd+(is-1)Nd+id,(vii)
the optical coefficient index:k=(iμκ-1)Ne+ie, where iμκ is 1 (the absorption coefficient) or 2 (the diffusion coefficient).Algorithm 1:Numerical algorithm (continuous version).
(I) Compute the solutionΦ(x) and the Robin function R(x) and its first and second derivatives.(II) Findδx† by solving ℛ1(δx†)=Φ(x+δx)-Φ(x) as in (1.6).(III) FindδxΔ=δxB-δx† by solving ℛ1(δxΔ)=-ℛ2(δx†) as in (2.2).(IV) ComputeδxB by adding δx† and δxΔ.If we use piecewise linear or bilinear finite element method, the finite element solution is represented by(3.2)uh(x)=∑in=1Nnuh(in)ϕin(x),
where ϕin is the piecewise linear or the bilinear function which is 1 on the inth node and 0 on all the other nodes. Assume μa and κ are piecewise constant function, which is constant for each Ne finite elements. Therefore, in diffuse optical tomography inverse problem, we have NωNsNd measurement information contents and 2Ne variables to find.We should discretizeℛ1 and ℛ2 to obtain a discretized version of Algorithm 1. Let the Jacobian and Hessian matrices, which is the discretization of integral operators ℛ1 and ℛ2, be J and H. The relation between higher order derivatives for the diffusion operator and higher order terms of Born expansions including ℛ1 and ℛ2 is analyzed in [7].Firstly, let us discretizeδx, Φ, and the Robin function R as follows: (3.3a)δx≈(∑ie=1NeδμieχTie,∑ie=1NeδκieχTie),(3.3b)Φiω,is≈∑in=1NnΦiniω,isϕin,(3.3c)Riω(⋅,ris)≈∑in=1NnRiniω,isϕin. Since we chose the source function qs as the Dirac delta function at the isth source point, Φiω,is=Riω(·,ris). However, we will discriminate these two functions in this paper, since they are different for general source function q which is different from the Dirac delta function. We will use δμ instead of δμa for notational convenience.Let the vectorγ0 which corresponds to the discretization of δx in (3.3a) be defined as (3.4)γ0=(δμ1,δμ2,…,δμNe,δκ1,δκ2,…,δκNe).
By the adjoint method [12], Riω(rid,·)=(Riω(·,rid))*, where * denotes complex conjugate. Likewise for (3.3a), let γ, γ†, γΔ, and γB be the discretization of δx, δx†, δxΔ, and δxB, respectively.For a functionf and a measurable set T, let us denote f∈T if the intersection of the support of f and T is not empty. The discretization of the linearized solution γ† is attained by solving the following equation: (3.5)Jγ†=b,
where (3.6)J(j,k)=∑ϕin1∈Tie∑ϕin2∈Tie(Rin1iω,id)*Eie(in1,in2)Φin2iω,iswheniμκ=1,J(j,k)=∑ϕin1∈Tie∑ϕin2∈Tie-3(Rin1iω,id)*κie2Fie(in1,in2)Φin2iω,iswheniμκ=2,b(j)=Φiω,is(x+δx)(rid)-Φiω,is(x)(rid),Eie(in1,in2)=∫Tieϕin1(ξ)ϕin2(ξ)dξ,Fie(in1,in2)=∫Tie∇ϕin1(ξ)⋅∇ϕin2(ξ)dξ.The discretized solutionγΔ is obtained by solving the following equation: (3.7)JγΔ=(γ†)tHγ†,
where (3.8)H(j,ie1,ie2)=∑ϕin1∈Tie1∑ϕin2∈Tie2(Rin1iω,id)*(Hμμ+Hμκ+Hκμ+Hκκ)(ie1,ie2;in1,in2)Φin2iω,is,
where Hμμ, Hμκ, Hκμ, and Hκκ are the discretization of corresponding terms in (3.1b) such that
(3.9)Hμμ(ie1,ie2;in1,in2)=∫Tie1∫Tie2ϕin1(ξ)Riω(ξ,η)ϕin2(η)dξdη,Hμκ(ie1,ie2;in1,in2)=∫Tie1∫Tie2ϕin1(ξ)∇ηRiω(ξ,η)⋅∇ηϕin2(η)dξdη,Hκμ(ie1,ie2;in1,in2)=∫Tie1∫Tie2∇ξϕin1(ξ)⋅∇ξ(Riω(ξ,η))ϕin2(η)dξdη,Hκκ(ie1,ie2;in1,in2)=∫Tie1∫Tie2∇ξϕin1(ξ)⋅[∇ξ∇ηRiω(ξ,η)]∇ηϕin2(η)dξdη.
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:Algorithm 2:Numerical algorithm (discretized version).
(I) Compute the solutionΦ(γ0)iniω,is and the Robin function R(γ0)id,iniω for iω=1,⋯,Nω,is=1,…,Ns,in=1,…,Nn 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γB by adding γ† and γΔ.


### 3.2. Approximation of Hessian

In this subsection we approximateHμμ, by assuming κ and μs′ are constant in Ω. The approximation is progressed in three ways.First, we approximate the Robin functionR(ξ,η) when (ξ,η)∈Ω∖∂Ω by its leading term R0(ξ,η) defined by (3.10)R0(ξ,η)={1(p-2)gpκ(η)|ξ-η|2-pp≥3,1ω2κ(η)log⁡(2S|ξ-η|)p=2,
where gp is the hypersurface area of the unit sphere in ℝp, p=2,3,… and S=supξ,η∈Ω|ξ-η|. Some important relations between R and R0 are found in [13].Second, whenie1≠ie2, the Robin function R and ϕin are approximated by constant values R0(c(ie1),c(ie2)) and ϕin(c(ie)) in Tie, respectively, where c(ie) of the center of the element Tie. That is to say, when ie1≠ie2, (3.9) is approximated as follows: (3.11)Hμμ(ie1,ie2;in1,in2)=R0(c(ie1),c(ie2))∫Tie1ϕin1(ξ)dξ∫Tie2ϕin2(η)dη.Third, whenie1=ie2, we use the following lemma.Lemma 3.1.
Let the measurable setT be contained in ℝp, p=2,3,…, and 0<m<p; then, the following inequality holds for T: (3.12a)∬T|ξ-η|-mdξdη≤p1-m/pp-mgpm/p|T|2-m/p,p≥2,(3.12b)∬Tlog⁡⁡(2S|ξ-η|)dξdη≤14π(1+log⁡(4S2π|T|))|T|2,p=2, where |T| is the volume of T.Proof.
If a ball with a radiusr has the same volume as T, we have
(3.13)r=(|T|pgp)1/p
for the space dimensions p=2,3,…. Let the ball of radius r with center ξ∈T be Bξ. Let T0=T∩Bξ, T+=T∖Bξ, and T-=Bξ∖T. Noting that |T+|=|T-|, we obtain
(3.14)∫T|ξ-η|-mdη=∫T0|ξ-η|-mdη+∫T+|ξ-η|-mdη≤∫T0|ξ-η|-mdη+∫T-|ξ-η|-mdη=∫Bξ|ξ-η|-mdη≤∫0rρp-m-1gpdρ=gpp-mrp-m≤gpp-m(|T|pgp)1-m/p
for all ξ∈T. Therefore,
(3.15)∬T|ξ-η|-mdηdξ≤gp|T|p-m(|T|pgp)1-m/p=p1-m/pp-mgpm/p|T|2-m/p.
Equation (3.12b) is derived in the same manner.Therefore, whenie1=ie2, (3.9) is approximated using the inequality in Lemma 3.1 as follows: (3.16)Hμμ(ie1,ie1;in1,in2)≈ϕin1(cie1)ϕin2(cie1)⋅{p2/p|Tie1|1+2/p2(p-2)gp2/pκ(c(ie1))p≥3,18π2κ(c(ie1))(1+log⁡(4S2π|Tie1|))|Tie1|2p=2.


## 3.1. Discretization

Algorithm1 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)
Nd detector positions: rid for id=1,2,…,Nd,(ii)
Ns source functions: qis=δis(Dirac delta function) for is=1,2,⋯,Ns,(iii)
Nω frequencies: ωiω for iω=1,2,…,Nω,(iv)
Ne elements: Tie for ie=1,2,…,Ne,(v)
Nn nodes: tin for in=1,2,…,Nn,(vi)
the measurement index:j=(iω-1)NsNd+(is-1)Nd+id,(vii)
the optical coefficient index:k=(iμκ-1)Ne+ie, where iμκ is 1 (the absorption coefficient) or 2 (the diffusion coefficient).Algorithm 1:Numerical algorithm (continuous version).
(I) Compute the solutionΦ(x) and the Robin function R(x) and its first and second derivatives.(II) Findδx† by solving ℛ1(δx†)=Φ(x+δx)-Φ(x) as in (1.6).(III) FindδxΔ=δxB-δx† by solving ℛ1(δxΔ)=-ℛ2(δx†) as in (2.2).(IV) ComputeδxB by adding δx† and δxΔ.If we use piecewise linear or bilinear finite element method, the finite element solution is represented by(3.2)uh(x)=∑in=1Nnuh(in)ϕin(x),
where ϕin is the piecewise linear or the bilinear function which is 1 on the inth node and 0 on all the other nodes. Assume μa and κ are piecewise constant function, which is constant for each Ne finite elements. Therefore, in diffuse optical tomography inverse problem, we have NωNsNd measurement information contents and 2Ne variables to find.We should discretizeℛ1 and ℛ2 to obtain a discretized version of Algorithm 1. Let the Jacobian and Hessian matrices, which is the discretization of integral operators ℛ1 and ℛ2, be J and H. The relation between higher order derivatives for the diffusion operator and higher order terms of Born expansions including ℛ1 and ℛ2 is analyzed in [7].Firstly, let us discretizeδx, Φ, and the Robin function R as follows: (3.3a)δx≈(∑ie=1NeδμieχTie,∑ie=1NeδκieχTie),(3.3b)Φiω,is≈∑in=1NnΦiniω,isϕin,(3.3c)Riω(⋅,ris)≈∑in=1NnRiniω,isϕin. Since we chose the source function qs as the Dirac delta function at the isth source point, Φiω,is=Riω(·,ris). However, we will discriminate these two functions in this paper, since they are different for general source function q which is different from the Dirac delta function. We will use δμ instead of δμa for notational convenience.Let the vectorγ0 which corresponds to the discretization of δx in (3.3a) be defined as (3.4)γ0=(δμ1,δμ2,…,δμNe,δκ1,δκ2,…,δκNe).
By the adjoint method [12], Riω(rid,·)=(Riω(·,rid))*, where * denotes complex conjugate. Likewise for (3.3a), let γ, γ†, γΔ, and γB be the discretization of δx, δx†, δxΔ, and δxB, respectively.For a functionf and a measurable set T, let us denote f∈T if the intersection of the support of f and T is not empty. The discretization of the linearized solution γ† is attained by solving the following equation: (3.5)Jγ†=b,
where (3.6)J(j,k)=∑ϕin1∈Tie∑ϕin2∈Tie(Rin1iω,id)*Eie(in1,in2)Φin2iω,iswheniμκ=1,J(j,k)=∑ϕin1∈Tie∑ϕin2∈Tie-3(Rin1iω,id)*κie2Fie(in1,in2)Φin2iω,iswheniμκ=2,b(j)=Φiω,is(x+δx)(rid)-Φiω,is(x)(rid),Eie(in1,in2)=∫Tieϕin1(ξ)ϕin2(ξ)dξ,Fie(in1,in2)=∫Tie∇ϕin1(ξ)⋅∇ϕin2(ξ)dξ.The discretized solutionγΔ is obtained by solving the following equation: (3.7)JγΔ=(γ†)tHγ†,
where (3.8)H(j,ie1,ie2)=∑ϕin1∈Tie1∑ϕin2∈Tie2(Rin1iω,id)*(Hμμ+Hμκ+Hκμ+Hκκ)(ie1,ie2;in1,in2)Φin2iω,is,
where Hμμ, Hμκ, Hκμ, and Hκκ are the discretization of corresponding terms in (3.1b) such that
(3.9)Hμμ(ie1,ie2;in1,in2)=∫Tie1∫Tie2ϕin1(ξ)Riω(ξ,η)ϕin2(η)dξdη,Hμκ(ie1,ie2;in1,in2)=∫Tie1∫Tie2ϕin1(ξ)∇ηRiω(ξ,η)⋅∇ηϕin2(η)dξdη,Hκμ(ie1,ie2;in1,in2)=∫Tie1∫Tie2∇ξϕin1(ξ)⋅∇ξ(Riω(ξ,η))ϕin2(η)dξdη,Hκκ(ie1,ie2;in1,in2)=∫Tie1∫Tie2∇ξϕin1(ξ)⋅[∇ξ∇ηRiω(ξ,η)]∇ηϕin2(η)dξdη.
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:Algorithm 2:Numerical algorithm (discretized version).
(I) Compute the solutionΦ(γ0)iniω,is and the Robin function R(γ0)id,iniω for iω=1,⋯,Nω,is=1,…,Ns,in=1,…,Nn 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γB by adding γ† and γΔ.


## 3.2. Approximation of Hessian

In this subsection we approximateHμμ, by assuming κ and μs′ are constant in Ω. The approximation is progressed in three ways.First, we approximate the Robin functionR(ξ,η) when (ξ,η)∈Ω∖∂Ω by its leading term R0(ξ,η) defined by (3.10)R0(ξ,η)={1(p-2)gpκ(η)|ξ-η|2-pp≥3,1ω2κ(η)log⁡(2S|ξ-η|)p=2,
where gp is the hypersurface area of the unit sphere in ℝp, p=2,3,… and S=supξ,η∈Ω|ξ-η|. Some important relations between R and R0 are found in [13].Second, whenie1≠ie2, the Robin function R and ϕin are approximated by constant values R0(c(ie1),c(ie2)) and ϕin(c(ie)) in Tie, respectively, where c(ie) of the center of the element Tie. That is to say, when ie1≠ie2, (3.9) is approximated as follows: (3.11)Hμμ(ie1,ie2;in1,in2)=R0(c(ie1),c(ie2))∫Tie1ϕin1(ξ)dξ∫Tie2ϕin2(η)dη.Third, whenie1=ie2, we use the following lemma.Lemma 3.1.
Let the measurable setT be contained in ℝp, p=2,3,…, and 0<m<p; then, the following inequality holds for T: (3.12a)∬T|ξ-η|-mdξdη≤p1-m/pp-mgpm/p|T|2-m/p,p≥2,(3.12b)∬Tlog⁡⁡(2S|ξ-η|)dξdη≤14π(1+log⁡(4S2π|T|))|T|2,p=2, where |T| is the volume of T.Proof.
If a ball with a radiusr has the same volume as T, we have
(3.13)r=(|T|pgp)1/p
for the space dimensions p=2,3,…. Let the ball of radius r with center ξ∈T be Bξ. Let T0=T∩Bξ, T+=T∖Bξ, and T-=Bξ∖T. Noting that |T+|=|T-|, we obtain
(3.14)∫T|ξ-η|-mdη=∫T0|ξ-η|-mdη+∫T+|ξ-η|-mdη≤∫T0|ξ-η|-mdη+∫T-|ξ-η|-mdη=∫Bξ|ξ-η|-mdη≤∫0rρp-m-1gpdρ=gpp-mrp-m≤gpp-m(|T|pgp)1-m/p
for all ξ∈T. Therefore,
(3.15)∬T|ξ-η|-mdηdξ≤gp|T|p-m(|T|pgp)1-m/p=p1-m/pp-mgpm/p|T|2-m/p.
Equation (3.12b) is derived in the same manner.Therefore, whenie1=ie2, (3.9) is approximated using the inequality in Lemma 3.1 as follows: (3.16)Hμμ(ie1,ie1;in1,in2)≈ϕin1(cie1)ϕin2(cie1)⋅{p2/p|Tie1|1+2/p2(p-2)gp2/pκ(c(ie1))p≥3,18π2κ(c(ie1))(1+log⁡(4S2π|Tie1|))|Tie1|2p=2.


## 4. Numerical Implementation

In the numerical implementation, the following parameters are used:(i)
Ω=[0,6]×[0,6](cm2),(ii)
Nd=16,(iii)
Ns=16,(iv)
Nω=1,(v)
Nx=Ny=16,(vi)
Ne=Nx*Ny,(vii)
Nn=(Nx+1)*(Ny+1),(viii)
μa=0.05+(0.2-0.05)χD(cm-1),(ix)
μs′=8(cm-1),(x)
κ=1/3*(μa+μs′)=1/3*(0.05+8),(xi)
ω=2π*300 MHz,(xii)
a=1,(xiii)
Jindex=Nx*Ny*0.4.Since the diffusion coefficientκ is constant, the right-hand side b is a Ns*Nd column vector, Jacobian J is a (Ns*Nd)×Ne matrix, the Hessian H is a Ne×(Ns*Nd)×Ne third-order tensor, and (γ†)tHγ† is Ns*Nd column vector in (3.5) and (3.7). H=Hμμ is approximated by (3.11) and (3.16).In the above setting, we reconstruct the obstacleD which has different absorption coefficient (0.2 cm−1) compared to the background absorption coefficient (0.05 cm−1). Four cases of the obstacle D are considered in Figures 1, 2, 3, and 4. The reconstruction of the absorption coefficient μa=0.05+(0.2-0.05)χD (cm−1) 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 μa and source/detector locations are plotted. The initial guess (μa0 or γ0) 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 (μa† or γ†) and the second approximation (μaB or γB) are plotted, respectively.Figure 1
Jindex=Nx*Ny*0.4, Jalpha=6.4920e-009, 10% noise, sources (*), and detectors (o).Figure 2
Jindex=Nx*Ny*0.4, Jalpha=6.5711e-009, 10% noise, sources (*), and detectors (o).Figure 3
Jindex=Nx*Ny*0.3, Jalpha=1.8227e-7, 10% noise, sources (*) and detectors (o).Figure 4
Jindex=Nx*Ny*0.4, Jalpha=6.5029e-009, 10% noise, sources (*), and detectors (o).In all four cases, 10% noise is added. Truncated singular value decomposition(SVD) is used.Jindex 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 Jindexth 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 Figures1 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.



---

*Source: 637209-2012-03-04.xml*