Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012, Article ID 637209, 15 pages
Research Article

The Second-Order Born Approximation in Diffuse Optical Tomography

Department of Mathematics, Dongguk University, Seoul 100715, Republic of Korea

Received 21 October 2011; Accepted 8 December 2011

Academic Editor: Chang-Hwan Im

Copyright © 2012 Kiwoon Kwon. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 [25] 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 𝑥0 about the structural optical coefficients 𝑥 and the perturbation of the optical coefficients 𝛿𝑥=𝑥𝑥0, 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(𝑥,𝛿𝑥)+2𝑥,(𝛿𝑥)2+,(1.3) where 0(𝑥)=Φ(𝑥),𝑖(𝛿𝑥)𝑖=𝛿𝑥,𝑖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 1 are different in the following sense: 1(𝛿𝑥)=𝛿𝑥,0=(𝛿𝑥,Φ).(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𝛿𝑥=Φ(𝑥+𝛿𝑥)Φ(𝑥).(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𝛿𝑥𝛿𝑥𝒜𝐶𝛿𝑥2𝒜,(1.7a)𝛿𝑥𝐵𝛿𝑥𝒜𝐶𝐵𝛿𝑥3𝒜,(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 𝑂(𝛿𝑥3𝒜), which is higher than the order of the error of the linearized solution 𝑥, 𝑂(𝛿𝑥2𝒜).

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 1𝛿𝑥𝐵=(Φ(𝑥+𝛿𝑥)Φ(𝑥))2𝛿𝑥2,(2.1) or equivalently, 1𝛿𝑥𝐵𝛿𝑥=2𝛿𝑥2.(2.2) In this section, we analyze the error for the linearized solution 𝛿𝑥 and the suggested solution 𝛿𝑥𝐵.

Let =𝐻1(Ω); 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: 1𝒜𝒜×Φ,(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 (1) from to 𝒜 such that (1)(1)=𝑖𝑑𝒜. (1) is usually called the left inverse of 1. Let us denote 𝛿𝑥=𝛿𝑥𝒜,Φ(𝑥)=Φ(𝑥),=𝒜×,𝑖=𝑖𝒜𝑖,𝑖=1,2,,1=1𝒜,(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 (1) such that (𝑅1)1=𝑖𝑑 and (1) is bounded, and let𝛿𝑥12.(2.7a) Then, 𝛿𝑥𝛿𝑥𝐶𝛿𝑥2,(2.7b)𝛿𝑥𝐵𝛿𝑥𝐶𝐵𝛿𝑥3,(2.7c) where 𝐶=212Φ,𝐶𝐵=𝐶34+𝐶2+𝐶=𝐶𝐶2+12.(2.8)

Proof. By (1.3) and (1.6), we obtain 1𝛿𝑥𝛿𝑥=2(𝛿𝑥)2+3(𝛿𝑥)3+.(2.9) Therefore we arrive at (2.7b) by the following inequality: 𝛿𝑥𝛿𝑥1(𝛿𝑥)2Φ1(𝛿𝑥)𝐶𝛿𝑥2.(2.10) From (2.7a) and (2.7b), we obtain the following upper bound of 𝛿𝑥: 𝛿𝑥1+𝐶𝛿𝑥𝛿𝑥1+1+Φ𝛿𝑥.(2.11) Using (2.2) and (2.9), we obtain 1𝛿𝑥𝛿𝑥𝐵=2(𝛿𝑥)22𝛿𝑥2+3(𝛿𝑥)3+4(𝛿𝑥)4+.(2.12) The second-order term on the righthand side of (2.12) is analyzed as follows: 2(𝛿𝑥)22𝛿𝑥2=(𝛿𝑥,(𝛿𝑥,Φ))𝛿𝑥,𝛿𝑥,Φ=𝛿𝑥,𝛿𝑥𝛿𝑥,Φ+𝛿𝑥𝛿𝑥,𝛿𝑥,Φ.(2.13) From (2.12), we obtain 𝛿𝑥𝛿𝑥𝐵12(𝛿𝑥)22𝛿𝑥2+3(𝛿𝑥)3+4(𝛿𝑥)4+.(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𝛿𝑥𝛿𝑥𝛿𝑥+𝛿𝑥+(𝛿𝑥)31𝛿𝑥1Φ2𝛿𝑥3𝐶2+1Φ+2𝐶𝛿𝑥3𝐶1+𝐶4+𝐶𝐵𝛿𝑥3.(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 1 and 2, which is introduced in (1.5), is as follows: 1(𝛿𝑥)=𝜇𝑎𝛿𝜇𝑎,Φ+𝜅(𝛿𝜅,Φ),(3.1a)2(𝛿𝑥)=𝜇𝑎𝛿𝜇𝑎,𝜇𝑎𝛿𝜇𝑎,Φ+𝜇𝑎𝛿𝜇𝑎,𝜅(𝛿𝜅,Φ),+𝜅𝛿𝜅,𝜇𝑎𝛿𝜇𝑎,Φ+𝜅𝛿𝜅,𝜅(𝛿𝜅,Φ).(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).

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 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 𝐽 and 𝐻. 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 𝛿𝑥, Φ, 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 𝛾0 which corresponds to the discretization of 𝛿𝑥 in (3.3a) be defined as 𝛾0=𝛿𝜇1,𝛿𝜇2,,𝛿𝜇𝑁𝑒,𝛿𝜅1,𝛿𝜅2,,𝛿𝜅𝑁𝑒.(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𝜅2𝑖𝑒𝐹𝑖𝑒𝑖𝑛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:

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 𝑅0(𝜉,𝜂) defined by 𝑅0(𝜉,𝜂)=1(𝑝2)𝑔𝑝𝜅(𝜂)||𝜉𝜂||2𝑝𝑝3,1𝜔2𝜅(𝜂)log2𝑆||𝜉𝜂||𝑝=2,(3.10) where 𝑔𝑝 is the hypersurface area of the unit sphere in 𝑝, 𝑝=2,3, and 𝑆=sup𝜉,𝜂Ω|𝜉𝜂|. Some important relations between 𝑅 and 𝑅0 are found in [13].

Second, when 𝑖𝑒1𝑖𝑒2, the Robin function 𝑅 and 𝜙𝑖𝑛 are approximated by constant values 𝑅0(𝑐(𝑖𝑒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=𝑅0𝑐𝑖𝑒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)𝑇log2𝑆||𝜉𝜂||𝑑𝜉𝑑𝜂14𝜋1+log4𝑆2𝜋||𝑇||||𝑇||2,𝑝=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 𝑇0=𝑇𝐵𝜉, 𝑇+=𝑇𝐵𝜉, and 𝑇=𝐵𝜉𝑇. Noting that |𝑇+|=|𝑇|, we obtain 𝑇||𝜉𝜂||𝑚𝑑𝜂=𝑇0||𝜉𝜂||𝑚𝑑𝜂+𝑇+||𝜉𝜂||𝑚𝑑𝜂𝑇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𝜋2𝜅𝑐𝑖𝑒11+log4𝑆2𝜋||𝑇𝑖𝑒1||||𝑇𝑖𝑒1||2𝑝=2.(3.16)

4. Numerical Implementation

In the numerical implementation, the following parameters are used: (i)Ω=[0,6]×[0,6](cm2),(ii)𝑁𝑑=16,(iii)𝑁𝑠=16,(iv)𝑁𝜔=1,(v)𝑁𝑥=𝑁𝑦=16,(vi)𝑁𝑒=𝑁𝑥𝑁𝑦,(vii)𝑁𝑛=(𝑁𝑥+1)(𝑁𝑦+1),(viii)𝜇𝑎=0.05+(0.20.05)𝜒𝐷(cm1),(ix)𝜇𝑠=8(cm1),(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−1) compared to the background absorption coefficient (0.05 cm−1). Four cases of the obstacle 𝐷 are considered in Figures 1, 2, 3, and 4. The reconstruction of the absorption coefficient 𝜇𝑎=0.05+(0.20.05)𝜒𝐷 (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 𝜇𝑎 and source/detector locations are plotted. The initial guess (𝜇𝑎0 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 (𝜇𝑎 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.


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


  1. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Problems, vol. 25, pp. 1–59, 2009. View at Google Scholar · View at Zentralblatt MATH
  2. V. Isakov, “On uniqueness in the inverse transmission scattering problem,” Communications in Partial Differential Equations, vol. 15, no. 11, pp. 1565–1587, 1990. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  3. K. Kwon, “Identification of anisotropic anomalous region in inverse problems,” Inverse Problems, vol. 20, no. 4, pp. 1117–1136, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  4. K. Kwon and D. Sheen, “Anisotropic inverse conductivity and scattering problems,” Inverse Problems, vol. 18, no. 3, pp. 745–756, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  5. J. Sylvester and G. Uhlmann, “A global uniqueness theorem for an inverse boundary value problem,” Annals of Mathematics, vol. 125, no. 1, pp. 153–169, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  6. H. Kang, K. Kwon, and K. Yun, “Recovery of an inhomogeneity in an elliptic equation,” Inverse Problems, vol. 17, no. 1, pp. 25–44, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  7. K. Kwon and B. Yazıcı, “Born expansion and Fréchet derivatives in nonlinear diffuse optical tomography,” Computers & Mathematics with Applications, vol. 59, no. 11, pp. 3377–3397, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  8. V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. I. Fourier-Laplace inversion formulas,” Journal of the Optical Society of America A, vol. 18, no. 6, pp. 1336–1347, 2001. View at Publisher · View at Google Scholar
  9. G. Murat, K. Kwon, B. Yazici, E. Giladi, and X. Intes, “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: I,” Inverse Problems, vol. 23, no. 3, pp. 1115–1133, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. G. Murat, K. Kwon, B. Yazici, E. Giladi, and X. Intes, “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging. II,” Inverse Problems, vol. 23, no. 3, pp. 1135–1160, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  11. K. Kwon, B. Yazıcı, and M. Guven, “Two-level domain decomposition methods for diffuse optical tomography,” Inverse Problems, vol. 22, no. 5, pp. 1533–1559, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems, vol. 15, no. 2, pp. R41–R93, 1999. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  13. C. Miranda, Partial Differential Equations of Elliptic Type, Springer, New York, NY, USA, 1970.