Direct Sampling Method for Diffusive Optical Tomography
aa r X i v : . [ m a t h - ph ] O c t Direct Sampling Method for Diffusive Optical Tomography
Yat Tin Chow ∗ Kazufumi Ito † Keji Liu ∗ Jun Zou ‡ Abstract
In this work, we are concerned with the diffusive optical tomography (DOT) problem in the casewhen only one or two pairs of Cauchy data is available. We propose a simple and efficient directsampling method (DSM) to locate inhomogeneities inside a homogeneous background and solve theDOT problem in both full and limited aperture cases. This new method is easy to implement andless expensive computationally. Numerical experiments demonstrate its effectiveness and robustnessagainst noise in the data. This provides a new promising numerical strategy for the DOT problem.
Mathematics Subject Classification (MSC2000) : 35J67, 35R30, 65N21, 78A70, 78M25.
Keywords : direct sampling method, diffusive optical tomography.
Diffusive optical tomography (DOT) is a popular non-invasive imaging technique that measures theoptical properties of a medium and creates images which show the distribution of absorption coefficientinside the body. In medical applications, such tomography is usually done in the near-infrared (NIR)spectral window. Chromophores in NIR window such as oxygenated and deoxygenated hemoglobin,water and lipid, are abundant in our body tissue, and a weighted sum of their contributions gives adifferent absorption coefficient [6]. This provides us with detailed information about the concentrationof various substances inside the tissue, thus revealing any pathological situation. In particular, it isknown that cancerous tumors absorb light more than the surrounding tissues, therefore a very importantmedical application of DOT is the early diagnosis of breast cancer. With the development of cost-efficient, compact and portable commercial instruments, we are now able to obtain accurate absorptionmeasurements for bedside monitoring with a relatively low cost. This propels the development of fastand accurate numerical algorithms to reconstruct the internal distribution of absorption coefficient so asto retrieve tissue properties in an accurate manner. Medical applications of DOT include breast cancerimaging [7, 9, 11, 12, 13, 18, 22, 26, 27, 32, 33, 42], brain functional imaging [10], stroke detection [3, 8, 17],muscle functional studies [19, 39], photodynamic therapy [37, 38], and radiation therapy monitoring [36].More information about medical application can be found in [6, 16].Over the past few decades, much effort has been made to develop efficient numerical algorithms tosolve the DOT problem. Both the static and time-resolved models have been extensively studied. For thetime-resolved model, the problem is also known as photon migration, and it has been widely consideredin literature, for example [5, 23, 24]. A popular approach to solve the DOT problem is to formulate aperturbation model, which comes from the linearization of the inverse problem [31]. Another standardapproach reduces the problem to an integral equation of the first kind, the kernel of which being generatedby the Green’s functions; see for instance [1, 2, 5, 29]. It is, however, well known that such an integralequation is ill-posed. Therefore, the least-squares minimization with regularization is a usual technique ∗ Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong. ([email protected],[email protected]). † Department of Mathematics and Center for Research in Scientific Computation, North Carolina State University,Raleigh, North Carolina ([email protected]). ‡ Department of Mathematics, Chinese University of Hong Kong, Shatin, N.T., Hong Kong. The work of this author wassubstantially supported by Hong Kong RGC grants (projects 405513 and 404611). ([email protected]).
1o handle the problem. Born type iterative method was introduced in [40], based on the minimizationof such a functional by regularized conjugate gradient method. Nonetheless, a major drawback of theseapproaches to image small inclusions is the necessity to use a fine mesh, resulting in a time-consumingprocess to invert highly ill-conditioned dense matrices. Various other approaches have been investigated,e.g., the fast Fourier transform methods [35], the multigrid and Bayesian methods [28, 41]. Effort hasalso been made to the three-dimensional reconstructions in optical tomography [4, 34], as well as how todramatically reduce the high computational complexity in handling large data set in solving the inverseproblem [25]. Other alternative directions have also been developed, including the conversion of theoriginal inverse problem to a well-posed system of elliptic partial differential equations [20, 21], or to aVolterra-type integral differential equation [30].In this work, we develop a novel direct sampling method (DSM) to solve the DOT problem. AssumeΩ ⊂ R d ( d = 2 ,
3) is an open bounded connected domain with a piecewise C boundary, representingthe absorption medium. Suppose that µ ∈ L ∞ (Ω) is a non-negative function representing the absorptioncoefficient in Ω, and µ is the absorption coefficient of the homogeneous background medium. We shalloften write the support of µ − µ as D , which represents the inhomogeneous inclusions sitting inside Ω.Then the potential u ∈ H (Ω) which represents the photon density field in the DOT model satisfies thefollowing Neumann problem: − ∆ u + µu = 0 in Ω ,∂u∂ν = g on ∂ Ω , (1.1)where g ∈ H − ( ∂ Ω) represents the surface flux along the boundary of Ω.Assume that we have only a partial measurement of the surface potential u | Γ at a relatively opensubset of the boundary Γ ⊂ ∂ Ω , while the surface potential of ∂ Ω \ Γ is unknown to us. Here andthroughout the paper, we shall often write the surface potential over Γ as f , i.e., f := u | Γ . Our inversemedium problem is then formulated as follows: given a single or a small number (e.g., 2 to 3) of pairs( g, f ) of the Cauchy data, we recover the locations and the number of inhomogeneities D inside thehomogeneous background absorption medium.We first introduce a family of probing functions { η x } x ∈ Ω based on the monopole potential. Using theseprobing functions, we can then define an index function I ( x ) which gives a clear image of the location of D inside Ω. From this index function, we can determine the number and locations of inclusions effectively.We then move on to give an alternative characterization of the index function, which help enhance theimage quality as well as reduce computational complexity. The new method is fast, reliable and cheap,so it can serve as an effective initialization for any iterative refinement algorithms, such as the Born typeiteraiton [40] and other more refined optimisation type methods [28, 41]. The new DSM is very differentfrom the previous DSMs [14, 15] developed for inverse scattering problems in the sense that the probingfunctions which we introduce here may be regarded as the dual of the Green’s function for the originalforward problem (1.1) under an appropriate choice of Sobolev scale. Therefore our index function is inthe dual form with a different Sobelov scaling instead of the original L inner product form adopted inthe existing DSMs [14, 15].The paper is organized as follows. In section 2, we introduce the general philosophy of the DSM aswell as defining the probing and index functions that we shall use in our DSM method for solving theDOT problem. Then we give a representation of the index function in section 3. Explicit expressions ofthe probing function in some special domains are given in section 4, therefore giving us a clear pictureon the performance of our newly proposed method in these special domains. In section 5, we introducea technique for improving the imaging as well as the computational complexity of the index function.Numerical results are presented in section 6, and some concluding remarks are given in section 7.2 The direct sampling method
In this subsection, we give a brief introduction of the general philosophy of the DSM. Direct samplingtype methods were studied in [14, 15] for inverse scattering problems, based on the well-known fact thatthe scattered field can be approximated by a finite sum of the fundamental solutions centered at theinhomogeneous scatterers. With the help of a family of probing functions which is nearly orthogonal insome inner product space, an index function can be defined in such a way that it attains a large valueinside the inhomogeneous scatterers but a small value outside. This index function has been shown tobe a very effective method to reconstruct clustered scatterers in two- and three-dimensional scatteringmedia with a limited number of incident plane waves. In this work, we shall develop a novel DSM forsolving a very different inverse problem, the DOT problem, based on some important probing functions,which may be regarded as the dual of the Green’s function for the original forward problem (1.1) underan appropriate choice of Sobolev scale. In this sense, the new index function here is introduced in a dualform with a proper Sobolev scaling, instead of the standard L inner product used in the previous DSMs[14, 15].In what follows, we derive a general framework of the DSM for the DOT problem. We first considerthe case when the absorption coefficient of homogeneous background is non-zero ( µ = 0). We shall oftenwrite by u the potential satisfying the diffusion equation with homogeneous coefficient µ , i.e., − ∆ u + µ u = 0 in Ω ; ∂u ∂ν = g on ∂ Ω (2.1)and by f the surface potential of u along Γ, i.e. f := u | Γ . We will also denote by G x the Green’sfunction of the homogeneous diffusion equation: − ∆ G x + µ G x = δ x in Ω ; ∂G x ∂ν = 0 on ∂ Ω . (2.2)Now if u is the solution to (1.1), it follows from the definition of u that − ∆( u − u ) + µ ( u − u ) + ( µ − µ ) u = 0 in Ω ; ∂ ( u − u ) ∂ν = 0 on ∂ Ω . (2.3)Then by the Green’s representation, we obtain the well-known Lippmann-Schwinger representation forthe scattered potential f − f of the DOT problem:( f − f )( ξ ) = Z Ω G y ( ξ ) c ( y ) dy = Z D G y ( ξ ) c ( y ) dy ∀ ξ ∈ Γ , (2.4)where function c is defined by c ( y ) := − ( µ ( y ) − µ ) u ( y ) , y ∈ Ω . (2.5)Using (2.4) and some numerical quadrature rule, we can approximate the scattered potential f − f bya finite sum of fundamental solutions in the following form,( f − f )( ξ ) ≈ X k a k G x k ( ξ ) ∀ ξ ∈ Γ , (2.6)where { x k } are some quadrature points located inside D and { a k } are some coefficients.With the above approximation of the scattered field, we shall construct some index functions whichmay help locate the inhomogeneities of the medium. In the subsequent discussion, we will often write ∇ Γ , ∆ Γ and dσ Γ as the surface gradient, the surface Laplace operator and the volumetric element on the3urface Γ. Consider the following duality product h · , · i s which is well-defined for the space H s (Γ) × L (Γ)for some s ∈ R as follows: h φ, ψ i s := Z Γ ( − ∆ Γ ) s φ ψ dσ Γ for all φ ∈ H s (Γ) , ψ ∈ L (Γ) , (2.7)then the space L (Γ) can be considered as a subspace of the dual space of H s (Γ) in the following sense: h · , ψ i s ∈ (cid:0) H s (Γ) (cid:1) ∗ ∀ ψ ∈ L (Γ) . (2.8)Assume that | · | Y is an algebraic function of semi-norms in H s (Γ), and we can select a set of probingfunctions { η x } x ∈ Ω ⊂ H s (Γ) such that they are nearly orthogonal to the family { G y | Γ } y ∈ Ω with respectto h· , ·i s and | · | Y in the following sense, namely for any y ∈ Ω, the function x K ( x, y ) := h η x , G y i s | η x | Y , x ∈ Ω (2.9)attains maximum at x = y and decays when x moves away from y . Under this assumption, we are nowready to define the following index function I ( x ) := h η x , f − f i s | η x | Y , x ∈ Ω . (2.10)Directly substituting (2.6) into (2.10), we arrive at I ( x ) = h η x , f − f i s | η x | Y ≈ X k a k h η x , G x k i s | η x | Y = X k a k K ( x, x k ) , (2.11)where K is the function defined as in (2.9).For the case when the absorption coefficient of the homogeneous background vanishes, namely µ = 0,a similar representation of the index function I can be obtained following the same argument as above.In this case, the Green’s function G x of the homogeneous diffusion equation is of the form: − ∆ G x = δ x in Ω ; ∂G x ∂ν = 1 | ∂ Ω | on ∂ Ω ; Z ∂ Ω G x dσ = 0 , (2.12)whereas the incident potential u satisfies − ∆ u = 0 in Ω ; ∂u ∂ν = g on ∂ Ω ; Z ∂ Ω u dσ = 0 . (2.13)Similarly to the previous case, we can see that the difference u − u satisfies (2.3) with µ = 0. Hence weobtain the following Lippmann-Schwinger representation of f − f by the Green’s representation: ( f − f )( ξ ) = Z Ω G y ( ξ ) c ( y ) dy − | ∂ Ω | Z ∂ Ω ( f − f ) dσ = Z D G y ( ξ ) c ( y ) dy − | ∂ Ω | Z ∂ Ω f dσ, ξ ∈ Γ , (2.14) where the function c is defined as in (2.5). Therefore, as in the previous case, we obtain the followingapproximation of the scattered potential f − f :( f − f )( ξ ) ≈ X k a k G x k ( ξ ) − | ∂ Ω | Z ∂ Ω f dσ , ξ ∈ Γ , (2.15)for some quadrature points { x k } located inside D and some coefficients { a k } . Suppose that the familyof probing functions { η x } x ∈ Ω ⊂ H s (Γ) have the following property: h η x , i s = 0 ∀ x ∈ Ω , (2.16)4here the notation stands for the constant function ( ξ ) = 1 for all ξ ∈ Γ. Substituting this expressioninto (2.10), we obtain the following for the index function I ( x ): I ( x ) = h η x , f − f i s | η x | Y ≈ X k a k h η x , G x k i s | η x | Y − R ∂ Ω f dσ | ∂ Ω | h η x , i s | η x | Y = X k a k K ( x, x k ) , (2.17)where K is the function defined as in (2.9).From representations (2.11) and (2.17) of the index function I ( x ), we can see that the magnitude of I ( x ) is relatively large inside D , while it is relatively small outside. Therefore, if the magnitude of theindex function I ( x ) is relatively large at a point x ∈ Ω , it is most likely that the point x lies inside D .On the contrary, if the magnitude of I ( x ) is relatively small at x , it is then very likely that the point is atthe homogeneous background. Therefore the index function provides us with an estimate of the locationof D , and thus the number and locations of inclusions. In this subsection, we introduce a family of probing functions and the index function which we willuse in the DSM for the DOT problem. For this purpose, we first define, for any given x ∈ Ω, a function w x satisfying the following system − ∆ w x + µ w x = δ x in Ω ; w x = 0 on Γ ; ∂w x ∂ν = 0 on ∂ Ω \ Γ . (2.18)For convenience, we may sometimes write the function w x ( · ) as w ( x, · ) . Now, given any x ∈ R d , let Φ x be the fundamental solution to the homogeneous diffusion equation, i.e., − ∆Φ x ( y ) + µ Φ x ( y ) = δ x ( y ) , y ∈ R d (2.19)subject to the decay condition: Φ x ( y ) → | y | → ∞ . (2.20)Subtracting (2.18) from (2.19), we readily obtain a representation of the function w x based on thefundamental solution Φ x : w x = Φ x − ψ x , (2.21)where ψ x is a function satisfying the following system: − ∆ ψ x + µ ψ x = 0 in Ω ; ψ x = Φ x on Γ ; ∂ψ x ∂ν = ∂ Φ x ∂ν on ∂ Ω \ Γ . (2.22)With these notions, we are now ready to define the family of probing functions for our purpose. For afixed point x ∈ Ω, we define the probing function η x as the surface flux of w x over Γ, i.e., η x ( ξ ) := ∂w x ∂ν ( ξ ) , ξ ∈ Γ . (2.23)The probing function can be treated as the dual of the Green’s function for the original forward problem(1.1) along the boundary Γ under a proper choice of Sobolev scaling. Some mathematical justificationand detailed analyses for this choice of probing functions will be provided in section 4 for two specialdomains.Now, we shall define an index function that will be used to solve the DOT problem. To be morespecific, we confine ourselves to the case with s = 1, while similar definitions can be done for the index s = 1. Then the duality product h· , ·i for the space H (Γ) × L (Γ) can be explicitly given as the following: h φ, ψ i := − Z Γ ∆ Γ φ ψ dσ Γ , for all φ ∈ H (Γ) , ψ ∈ L (Γ) . (2.24)5s a remark, if Γ is a closed surface and function ψ ∈ H (Γ), then the above duality product actuallyequals to the well-known H (Γ) semi-inner product, i.e., h φ, ψ i = ( φ, ψ ) H (Γ) , where( φ, ψ ) H (Γ) := Z Γ ∇ Γ φ · ∇ Γ ψdσ Γ , for all φ, ψ ∈ H (Γ) . (2.25)One choice of the algebraic function |·| Y in (2.9) of semi-norms can be | · | / H (Γ) | · | / H (Γ) , with | · | H (Γ) := | · | L (Γ) and | · | H (Γ) := ( · , · ) / H (Γ) . Substituting these expressions into (2.9) and (2.10), we have thefollowing expressions of the kernel K and the index function I : K ( x, y ) := h η x , G y i | η x | H (Γ) | η x | H (Γ) ∀ x, y ∈ Ω , (2.26) I ( x ) := h η x , f − f i | η x | H (Γ) | η x | H (Γ) ∀ x ∈ Ω . (2.27)We will show in section 4 that, for two special domains Ω, namely the circular domain and the rectangulardomain, this index function has the desired property as stated in section 2.1: its magnitude is relativelylarge inside D but relatively small outside. Using this property, the inhomogeneities of the absorptionmedium Ω can be clearly located. Unlike the most conventional methods for the DOT problem, which usethe regularized least-squares methods or iterative methods, this method is computationally inexpensive,since it involves only the evaluation of a duality product h · , · i , which is an integral over the measure-ment surface Γ. Moreover, the noise which enters into the boundary data is averaged and smoothedby integration, and no matrix inversion is required. In section 5, some technique will be introduced tofurther improve the quality of imaging as well as sufficiently reducing the computational cost by avoidingthe calculation of the family of probing functions. Therefore this method is expected to be efficient com-putationally and robust against the noise in the data, which are demonstrated by numerical examples insection 6. In this section we try to give a representation of the h · , · i s duality product between the probing func-tions η x and the measurement boundary data f − f , and therefore the index function. This representationprovides a more detailed understanding of the behaviour of the index function.We first consider the case with s = 0, then the duality product h · , · i is the standard L (Γ) innerproduct. Given a point x ∈ Ω , from (2.3), (2.18), (2.23) and the Green’s identity, we have the followingrepresentation: ( η x , f − f ) = Z Γ η x ( u − u ) dσ Γ = Z Γ ( u − u ) ∂w x ∂ν dσ Γ = Z ∂ Ω (cid:18) ( u − u ) ∂w x ∂ν − w x ∂ ( u − u ) ∂ν (cid:19) dσ = Z Ω (( u − u )∆ w x − w x ∆( u − u )) dy = − ( u − u )( x ) − Z D ( µ − µ ) uw x dy ≈ − ( u − u )( x ) − X k λ k ( µ ( y k ) − µ ) u ( y k ) w x ( y k ) , (3.1)6or some quadrature points y k sitting inside in the inhomogeneous support D and some quadratureweights λ k of the respective quadrature rule. Therefore we can find that if x is near the point y k forsome k , then the magnitude of the duality product h η x , f − f i is relatively large, otherwise it would berelatively small.Next, we consider the case with s = 1 and the duality product h· , ·i as defined in (2.24). We considerthe rectangular domain Ω := (0 , h ) × ( − L, L ) for some h, L ∈ R and Γ = { , h } × ( − L, L ) ⊂ ∂ Ω. Thesubsequent derivations apply naturally to d -dimensional rectangular domains.From (2.3), we can see that the function ∂ x x ( u − u ) satisfies the following equation: − ∆ (cid:0) ∂ x x ( u − u ) (cid:1) + µ ∂ x x ( u − u ) + ∂ x x [( µ − µ ) u ] = 0 in Ω . (3.2)Now assuming u − u ∈ H (Ω), then using that ∇ Γ = ∂ x and ∆ Γ = ∂ x x on Γ, we have for any x ∈ Ωby the Green’s identity on Γ the following duality product of the probing function η x and the scatteredpotential f − f : h η x , f − f i = − Z Γ ∆ Γ η x ( u − u ) dσ Γ = Z Γ ∇ Γ ( u − u ) · ∇ Γ η x dσ Γ − ( u − u )(0 , L ) ∂η x ∂x (0 , L ) + ( u − u )(0 , − L ) ∂η x ∂x (0 , − L ) − ( u − u )( h, L ) ∂η x ∂x ( h, L ) + ( u − u )( h, − L ) ∂η x ∂x ( h, − L ) . (3.3)However, from (2.18) and the fact that (0 ,
1) is normal to the surface ∂ Ω \ Γ, we have that ∂η x ∂x (0 , L ) = ∂ w x ∂x ∂x (0 , L ) = 0 . (3.4)Similarly we can derive ∂η x ∂x (0 , L ) = ∂η x ∂x (0 , − L ) = ∂η x ∂x ( h, L ) = ∂η x ∂x ( h, − L ) = 0 . (3.5)Therefore we have h η x , f − f i = Z Γ ∇ Γ ( u − u ) · ∇ Γ η x dσ Γ = Z Γ ∂ x ( u − u ) ∂ x η x dσ Γ = − Z Γ ∂ x x ( u − u ) η x dσ Γ + η x (0 , L ) ∂ ( u − u ) ∂x (0 , L ) − η x (0 , − L ) ∂ ( u − u ) ∂x (0 , − L )+ η x ( h, L ) ∂ ( u − u ) ∂x ( h, L ) − η x ( h, − L ) ∂ ( u − u ) ∂x ( h, − L ) . (3.6)However, again, from the fact that (0 ,
1) is normal to the surface ∂ Ω \ Γ, we get from (1.1) and (2.1) that ∂ ( u − u ) ∂x (0 , L ) = ∂u∂x (0 , L ) − ∂u ∂x (0 , L ) = g (0 , L ) − g (0 , L ) = 0 . (3.7)Similarly we can derive ∂ ( u − u ) ∂x (0 , L ) = ∂ ( u − u ) ∂x (0 , − L ) = ∂ ( u − u ) ∂x ( h, L ) = ∂ ( u − u ) ∂x ( h, − L ) = 0 . (3.8)7sing (3.6), (2.18), (2.23) and the Green’s identity, we obtain h η x , f − f i = − Z Γ ∂ x x ( u − u ) η x dσ Γ = − Z Γ ∂ x x ( u − u ) ∂w x ∂ν dσ Γ = − Z Γ (cid:18) ∂ x x ( u − u ) ∂w x ∂ν − w x ∂∂ν (cid:0) ∂ x x ( u − u ) (cid:1)(cid:19) dσ Γ = − Z ∂ Ω (cid:18) ∂ x x ( u − u ) ∂w x ∂ν − w x ∂∂ν (cid:0) ∂ x x ( u − u ) (cid:1)(cid:19) dσ = − Z Ω (cid:0) ∂ x x ( u − u )∆ w x − w x ∆ (cid:0) ∂ x x ( u − u ) (cid:1)(cid:1) dy + ϕ ( x )= ∂ x x ( u − u )( x ) + Z D ∂ x x [( µ − µ ) u ] w x dy (3.9) ≈ ∂ x x ( u − u )( x ) + X k λ k ∂ x x [( µ − µ ) u ] ( y k ) w x ( y k ) (3.10)for some quadrature points y k sitting inside in the inhomogeneous support D and some quadratureweights λ k of the respective quadrature rule. We can now observe that h η x , f − f i is relatively large if x is near y k for some k . Following the framework introduced in section 2.1, in order to justify the DSM for a given domainΩ, we should show that the family of probing functions defined as in (2.23) satisfies two properties: thefirst one is (2.16), which will then validate (2.17); the second one is that they are nearly orthogonal tothe family { G y | Γ } y ∈ Ω with respect to h· , ·i s and | · | Y as stated in Subsection 2.1. It is not easy todeduce these properties for a general domain Ω. In this section, we will calculate explicitly the probingfunction defined in (2.23) and the kernel as in (2.26) for two special domains Ω, or two most frequentlyused sampling domains: a circular domain and a rectangular domain, which help us establish the twonecessary properties to justify the DSM. Let Ω be the the circular domain B (0), i.e., the open ball of radius 1 centered at origin in R .Consider the boundary measurement on Γ = ∂B (0) = S , and the case with µ = 0 first. Then for x ∈ Ω, we can see from (2.18) that the function w x is actually the Green’s function for the unit disk withthe Dirichlet boundary condition: w ( x, y ) = 12 π (cid:18) log | x − y | − log (cid:12)(cid:12)(cid:12)(cid:12) x | x | − | x | y (cid:12)(cid:12)(cid:12)(cid:12) (cid:19) , y ∈ B (0) . (4.1)Therefore the probing function η x defined in (2.23) is the Poisson kernel with the explicit form: η x ( y ) = 1 − | x | π | x − y | , y ∈ S . (4.2)Next, we shall calculate the kernel (2.26) for the better understanding of the index function I definedas in (2.27). For the notational sake, we shall write the Fourier coefficient of a function φ ∈ L ( S ) as[ F ( φ )]( n ) := 12 π Z π φ ( θ ) e − inθ dθ , n ∈ Z . (4.3)8ow for any φ, ψ ∈ H ( S ), their duality product h · , · i is equal to the H semi-inner product and canbe expressed in the Fourier coefficients of φ and ψ as12 π h φ, ψ i = X n ∈ Z n [ F ( φ )]( n )[ F ( ψ )]( n ) , (4.4)From the fact that both the functions G z and η x are in H ( S ), we can obtain an explicit expression ofthe kernel K ( x, z ) for x, z ∈ B (0) by calculating the Fourier expansions of G z and η x . For this purpose,we consider for a given function g the following Neumann problem: − ∆ v = 0 in B (0) ; ∂v∂ν = g on S ; Z S v dσ = 0 . (4.5)On one hand, from the Green’s representation formula, function v can be expressed by v ( z ) = Z S g ( y ) G z ( y ) dσ ( y ) − π Z S v ( y ) dσ ( y ) = Z S g ( y ) G z ( y ) dσ ( y ) , y ∈ S , (4.6)where G z is defined as in (2.12). On the other hand, by a separation of variables, we can explicitlycalculate the Fourier expansion of function v as follows: v ( z ) = 12 π X n ∈ Z \{ } r | n | z e inθ z | n | Z π g ( θ y ) e − inθ y dθ y = 12 π Z S g ( y ) X n ∈ Z \{ } r | n | z | n | e in ( θ z − θ y ) dy , (4.7)where z = ( r z , θ z ) is in the polar coordinate. Then comparing (4.7) with (4.6), we obtain the Fourierexpansion of G z as G z ( y ) = 12 π X n ∈ Z \{ } r | n | z | n | e in ( θ z − θ y ) = 12 π X n ∈ Z \{ } r | n | z e − inθ z | n | e inθ y , y ∈ S . (4.8)With a similar technique, we try to obtain the Fourier expansions of w x . For a given function h , considerthe following Dirichlet problem: − ∆ v = 0 in B (0) ; v = h on S . Then, by the definition of w x in (2.23) and the Green’s representation formula, we can see that thefunction v can be expressed as v ( x ) = Z S h ( y ) ∂w x ∂ν ( y ) dσ ( y ) , x ∈ S . (4.9)On the other hand, we can obtain the following Fourier expansion of function v by a separation ofvariables: v ( x ) = 12 π X n ∈ Z r | n | x e inθ x Z π g ( θ y ) e − inθ y dθ y = 12 π Z S g ( y ) X n ∈ Z r | n | x e in ( θ x − θ y ) ! dy , x ∈ S , (4.10)where x = ( r x , θ x ) is in the polar coordinate. Hence, we can compare (4.10) with (4.9) to obtain thefollowing Fourier expansion for the probing function defined as in (2.23): η x ( y ) = ∂w x ∂ν ( y ) = 12 π X n ∈ Z r | n | x e in ( θ x − θ y ) = 12 π X n ∈ Z r | n | x e − inθ x e inθ y , y ∈ S . (4.11)From (4.8) and (4.11), we get the Fourier coefficients of G z and η x :[ F ( G z )]( n ) = 12 π r | n | z e − inθ z | n | ∀ n ∈ Z \{ } ; [ F ( G z )](0) = 0 , (4.12)9nd [ F ( η x )]( n ) = 12 π r | n | x e − inθ x ∀ n ∈ Z . (4.13)Substituting the above two expressions into (4.4), we readily obtain for all x, z ∈ B (0), h η x , G z i = 2 π X n ∈ Z n [ F ( η x )]( n )[ F ( G z )]( n )= 12 π X n ∈ Z \{ } | n | r | n | z r | n | x e in ( θ x − θ z ) = 1 π Re (cid:16) ∞ X n =1 nr nz r nx e in ( θ x − θ z ) (cid:17) = 1 π Re (cid:16) r z r x e i ( θ x − θ z ) ∞ X n =1 nr n − z r n − x e i ( n − θ x − θ z ) (cid:17) , (4.14)which gives h η x , G z i = 1 π Re (cid:16) r z r x e i ( θ x − θ z ) (cid:0) − r z r x e i ( θ x − θ z ) (cid:1) (cid:17) = 1 π r z r x cos( θ x − θ z )(1 + r z r x ) − r z r x (1 − r z r x cos( θ x − θ z ) + r z r x ) . (4.15)An interesting point to observe is that the duality product h η x , G z i is actually symmetric about x and z . Moreover, the L -norm of η x can be obtained from the Plancherel identity as | η x | H ( S ) = 2 π X n ∈ Z (cid:12)(cid:12)(cid:12) [ F ( η x )]( n ) (cid:12)(cid:12)(cid:12) = 12 π X n ∈ Z r | n | x = 12 π (cid:18) r x − r x + 1 (cid:19) = 12 π (cid:18) r x − r x (cid:19) , (4.16)whereas the H semi-norm of η x can be calculated explicitly by (4.4) as | η x | H ( S ) = 2 π X n ∈ Z n (cid:12)(cid:12)(cid:12) [ F ( η x )]( n ) (cid:12)(cid:12)(cid:12) = 12 π X n ∈ Z n r | n | x = 1 π ∞ X n =0 (cid:0) ( n + 2)( n + 1) − n + 1) + 1 (cid:1) r nx = 1 π − − r x ) + (1 − r x ) (1 − r x ) = 1 π r x ( r x + 1)(1 − r x ) . (4.17) −1 −0.5 0 0.5 1−1−0.8−0.6−0.4−0.200.20.40.60.81 00.10.20.30.40.50.60.70.8 Figure 1:
Kernel K ( · , z ) for a circular domain with z = ( − . , .
32) marked as a blue star.
Substituting (4.15)-(4.17) into (2.26), we obtain the explicit expression of kernel K for all x, z ∈ B (0): K ( x, z ) = h η x , G z i | η x | H (Γ) | η x | H (Γ) = C (cid:18) − r x r x (cid:19) (cid:18) (1 − r x ) r x ( r x + 1) (cid:19) (cid:18) r z r x cos( θ x − θ z )(1 + r z r x ) − r z r x (1 − r z r x cos( θ x − θ z ) + r z r x ) (cid:19) . (4.18)From this expression, we can observe that the kernel K ( x, z ) has a relatively large value when x ≈ z .This validates the DSM for the cirular domain B (0). We remark that we may not be able to obtain10ull boundary measurements in practice. Instead, we usually encounter the case where we only havemeasurements along the upper half of the boundary, i.e., when Γ = { x ∈ S ; x > } . Figure 1 showsthe values of distribution of kernel K ( x, z ) with this Γ and z = ( − . , .
32) for x ∈ B (0). This figureindeed shows the distribution kernel K ( x, z ) has a relatively large value when x ≈ z , and justifies theDSM for the circular domain with partial measurements.Now we consider the case with µ = 0. For a given h , consider the following Dirichlet problem: − ∆ v + µ v = 0 in B (0) ; v = h on S . Then, from the definition of w x as stated in (2.23) and the Green’s representation formula, we can seethat the function v can be expressed as v ( x ) = Z S h ( y ) ∂w x ∂ν ( y ) dσ ( y ) , x ∈ S . (4.19)However, by a separation of variables, we can obtain the following expansion of v for x ∈ S , v ( x ) = 12 π X n ∈ Z J n ( i √ µ r x ) J n ( i √ µ ) e inθ x Z π g ( θ y ) e − inθ y dθ y = 12 π Z S g ( y ) X n ∈ Z J n ( i √ µ r x ) J n ( i √ µ ) e in ( θ x − θ y ) ! dy , (4.20)where J n are the Bessel functions of first kind of order n and x = ( r x , θ x ) is in its polar coordinate.Comparing (4.19) with (4.20) we obtain the Fourier expansion for the probing function in (2.23): η x ( y ) = ∂w x ∂ν ( y ) = 12 π X n ∈ Z J n ( i √ µ r x ) J n ( i √ µ ) e in ( θ x − θ y ) , y ∈ S . (4.21)From the fact that J m ( t ) (cid:30) p π | m | (cid:18) et | m | (cid:19) | m | → m → ∞ (4.22)for all t >
0, we have J m ( i √ µ r x ) J m ( i √ µ ) (cid:30) r | m | x → m → ∞ . (4.23)Hence, we observe that the coefficient in the series decays exponentially. Therefore we may approximatethe probing function η x effectively by truncating the series (4.21). This gives us an efficient method tocalculate the probing function for µ = 0, which is very useful for practical purpose.Next, we shall calculate the kernel K ( x, z ) for x, z ∈ B (0) in (2.26). For this, we first compute theFourier expansions of G z defined as in (2.2) and η x . For a given g , consider the Neumann problem: − ∆ v + µ v = 0 in B (0) ; ∂v∂ν = g on S . (4.24)Again, from the Green’s representation formula, the function v can be expressed in the following form v ( z ) = Z S g ( y ) G z ( y ) dσ ( y ) , y ∈ S . (4.25)On the other hand, we can calculate the Fourier expansion of function v by a separation of variables: v ( z ) = Z π g ( θ y ) − J ( i √ µ r z ) iπ √ µ J ( i √ µ ) + X n ∈ Z \{ } iπ √ µ J | n | ( i √ µ r z ) J | n |− ( i √ µ ) − J | n | +1 ( i √ µ ) e in ( θ z − θ y ) dθ y , (4.26) z = ( r z , θ z ) is in the polar coordinate. Comparing (4.26) with (4.25), we obtain G z ( y ) = − J ( i √ µ r z ) iπ √ µ J ( i √ µ ) + X n ∈ Z \{ } iπ √ µ J | n | ( i √ µ r z ) J | n |− ( i √ µ ) − J | n | +1 ( i √ µ ) e in ( θ z − θ y ) , y ∈ S . (4.27)Then it follows from (4.21) and (4.27) that[ F ( G z )]( n ) = 1 iπ √ µ J | n | ( i √ µ r z ) J | n |− ( i √ µ ) − J | n | +1 ( i √ µ ) e − inθ z ∀ n ∈ Z \{ } , (4.28)[ F ( G z )](0) = − J ( i √ µ r z ) iπ √ µ J ( i √ µ ) , (4.29)[ F ( η x )]( n ) = 12 π J n ( i √ µ r x ) J n ( i √ µ ) e − inθ x ∀ n ∈ Z . (4.30)Substituting the above expressions into (4.4), we have the following duality product for all x, z ∈ B (0), h η x , G z i = 2 π X n ∈ Z n [ F ( η x )]( n )[ F ( G z )]( n )= 1 iπ √ µ X n ∈ Z \{ } n J n ( i √ µ r z ) J n ( i √ µ r x ) (cid:0) J n − ( i √ µ ) − J n +1 ( i √ µ ) (cid:1) J n ( i √ µ ) e in ( θ x − θ z ) , (4.31)whereas the L norm and H semi-norm of η x can be obtained similarly as above. In this subsection, we turn our attention to the case when Ω is a rectangular domain (0 , h ) × ( − L, L )in R for some h, L ∈ R and calculate the probing functions introduced in (2.23) and the kernel in (2.26).We consider the case of partial measurements along Γ = { , h } × ( − L, L ) ⊂ ∂ Ω with µ = 0. We knowthe fundamental solution Φ x for x = ( x , x ) ∈ Ω, satisfying (2.19)-(2.20) with µ = 0, is given byΦ x ( y ) = − π log | y − x | , y ∈ Ω . (4.32)Then we obtain from (2.18) and the reflection principle that the function w x is of the following form w ( x, y ) = − π X k =1 , X I ∈ Z log | y − x ( k ) I | + 12 π X k =3 , X I ∈ Z log | y − x ( k ) I | , y ∈ Ω , (4.33)where x ( k ) I is defined, for k = 1 , , , I = ( i , i ) ∈ Z , as follows: x (1) I = ( x + 2 h i , x + 4 L i ) ; x (2) I = ( x + 2 h i , L − x + 4 L i ) ; x (3) I = ( − x + 2 h i , x + 4 L i ) ; x (4) I = ( − x + 2 h i , L − x + 4 L i ) . Therefore, writing x ( k ) I = ( x ( k ) I, , x ( k ) I, ) for k = 1 , , , I = ( i , i ) ∈ Z , we have the following seriesexpansion for the probing function defined as in (2.23) for all y ∈ { h } × ( − L, L ) ⊂ Γ: η x ( y ) = ∂w x ∂y ( y ) = − π X k =1 , X I ∈ Z h − x ( k ) I, | ( h, y ) − x ( k ) I | − X k =3 , X I ∈ Z h − x ( k ) I, | ( h, y ) − x ( k ) I | , (4.34)12hereas for y ∈ { } × ( − L, L ) ⊂ Γ, with a similar argument, we have η x ( y ) = − ∂w x ∂y ( y ) = − π X k =1 , X I ∈ Z x ( k ) I, | (0 , y ) − x ( k ) I | − X k =3 , X I ∈ Z x ( k ) I, | (0 , y ) − x ( k ) I | . (4.35)This gives an explicit expression of the probing function η x for all x ∈ Ω. But the above series represen-tation of η x is tedious to work with if we intend to calculate the kernel K defined in (2.26). However, wecan substitute (4.32) into (2.21) to obtain the following representation of w x for x ∈ Ω: w ( x, y ) = − π log | x − y | − ψ ( x, y ) , x, y ∈ Ω , where ψ ( x, y ) := ψ x ( y ) is defined as in (2.22) and is a smooth function with respect to x, y ∈ Ω. Thenwe can derive the following expression of the probing function for x ∈ (0 , h ) × ( − L, L ): η x ( y ) = ∂w x ∂y ( y ) = − π h − x | ( h, y ) − x | − ∂ψ∂y ( x, y ) , y ∈ { h } × ( − L, L ) , (4.36) η x ( y ) = − ∂w x ∂y ( y ) = − π x | (0 , y ) − x | + ∂ψ∂y ( x, y ) , y ∈ { } × ( − L, L ) . (4.37)On the other hand, subtracting (2.12) from (2.19) and substituting (4.32) into it, we can obtain thefollowing representation of G z for z ∈ Ω based on the fundamental solution Φ z : G z ( y ) = Φ z ( y ) − χ ( z, y ) = − π log | z − y | − χ ( z, y ) , (4.38)where χ ( z, · ) is a function satisfying the following system: − ∆ χ ( z, · ) = 0 in Ω ; ∂χ ( z, · ) ∂ν = ∂ Φ z ∂ν − | Γ | on ∂ Ω ; Z Γ χ ( z, · ) dσ = Z Γ Φ z dσ , (4.39)and is a smooth function with respect to z, y ∈ Ω. From (4.36)-(4.38) and the Green’s identity, we cansee that the h · , · i duality product of η x and G z can be expressed in the following form for all x, z ∈ Ω: h η x , G z i = (I) + (II) (4.40)where the terms (I) and (II) as as follows: (I) = Z L − L ∂∂y (cid:18) − π h − x | ( h, y ) − x | − ∂ψ∂y (cid:16) x, ( h, y ) (cid:17)(cid:19) ∂∂y (cid:18) − π log | ( h, y ) − z | − χ (cid:16) z, ( h, y ) (cid:17)(cid:19) dy + Z L − L ∂∂y (cid:18) − π x | (0 , y ) − x | + ∂ψ∂y (cid:16) x, (0 , y ) (cid:17)(cid:19) ∂∂y (cid:18) − π log | (0 , y ) − z | − χ (cid:16) z, (0 , y ) (cid:17)(cid:19) dy = Z L − L (cid:18) π ( h − x )( y − x ) | ( h, y ) − x | − ∂ ψ∂y y (cid:16) x, ( h, y ) (cid:17)(cid:19) (cid:18) − π y − z | ( h, y ) − z | − ∂χ∂y (cid:16) z, ( h, y ) (cid:17)(cid:19) dy + Z L − L (cid:18) − π x ( y − x ) | (0 , y ) − x | + ∂ ψ∂y y (cid:16) x, (0 , y ) (cid:17)(cid:19) (cid:18) − π y − z | (0 , y ) − z | − ∂χ∂y (cid:16) z, (0 , y ) (cid:17)(cid:19) dy = C Z L − L (cid:18) ( h − x )( y − x )( y − z ) | ( h, y ) − x | | ( h, y ) − z | + x ( y − x )( y − z ) | (0 , y ) − x | | (0 , y ) − x | + ε ( x, z, y ) (cid:19) dy , (4.41)where C is a constant and ε ( x, z, y ) is a function whose order of magnitude is smaller than the firsttwo terms in the integrand; whereas with the help of (3.5), the second term can be calculated directly as (II) = (cid:20) ∂η x ∂x ( h, y ) G z ( h, y ) (cid:21) Ly = − L − (cid:20) ∂η x ∂x (0 , y ) G z (0 , y ) (cid:21) Ly = − L = 0 . (4.42)13sing the same technique, the L norm and H semi-norm of η x can be expressed as follows: | η x | H (Γ) = C Z L − L ( h − x ) | ( h, y ) − x | + x | (0 , y ) − x | + ε ( x, y ) dy , (4.43)while | η x | H (Γ) = C Z L − L ( h − x ) ( y − x ) | ( h, y ) − x | + x ( y − x ) | (0 , y ) − x | + ε ( x, y ) dy , (4.44)where, for i = 2 , ε i ( x, y ) are some functions with its order of magnitude being smaller than the firsttwo terms of the integrand, and C i are some constants. Substituting (4.40)-(4.44) into (2.26), we can seethat the kernel K ( x, z ) is relatively large when x ≈ z , while it is relatively small outside. Figure 2 showsthe values of the kernel K ( x, z ) with z = (0 . , − . x ∈ (0 , h ) × ( − L, L ). This supports us touse the index function defined in (2.27) for locating inhomogeneities in the rectangular domain for theDOT problem. −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 0.20.40.60.81
Figure 2:
Kernel K ( · , z ) on a rectangular domain with z = (0 . , − . In this section we shall present an alternative characterization of the index function defined in (2.27)for the following two special cases, or two most frequently used sampling domains: Ω is an arbitraryopen connected domain with Γ being a closed surface on ∂ Ω; and Ω = (0 , h ) × ( − L, L ) is a rectangulardomain with Γ = { , h } × ( − L, L ). With these characterizations, we will introduce an additional scalingfunction which will greatly improve the imaging of the index function as well as essentially reducing thecomputational effort for the index function.We first define an auxiliary function which will be very useful to represent the index function I . Let ϕ be a function which satisfies the following homogeneous diffusion system: − ∆ ϕ + µ ϕ = 0 in Ω ; ϕ = ∆ Γ ( f − f ) on Γ; ∂ϕ∂ν = 0 on ∂ Ω \ Γ . (5.1)Then from the definition of the duality product h · , · i and the Green’s identity on Γ, we have h η x , f − f i = − Z Γ η x ∆ Γ ( f − f ) dσ Γ − Z ∂ Γ ( f − f ) ∂η x ∂n dσ ∂ Γ + Z ∂ Γ η x ∂∂n ( f − f ) dσ ∂ Γ (5.2)for any x ∈ Ω, where dσ Γ and dσ ∂ Γ represent the respective volume and surface elements on Γ and ∂ Γwhile the vector n is the unit vector field normal to ∂ Γ.When Ω is an arbitrary open connected domain but Γ is a closed surface, we know ∂ Γ = ∅ . Then itfollows from (2.18), (2.23), (5.1)-(5.2) and the Green’s formula that h η x , f − f i = − Z Γ η x ∆ Γ ( f − f ) dσ Γ = − Z Γ ϕ ∂w x ∂ν dσ Γ = Z ∂ Ω (cid:18) w x ∂ϕ∂ν − ϕ ∂w x ∂ν (cid:19) dσ = Z Ω ( w x ∆ ϕ − ϕ ∆ w x ) dy = ϕ ( x ) . (5.3)14n the other hand, when Ω is a rectangular domain, Ω = (0 , h ) × ( − L, L ), and Γ = { , h } × ( − L, L ), wehave from (5.1)-(5.2) that h η x , f − f i = − Z Γ η x ∆ Γ ( f − f ) dσ Γ − ( u − u )(0 , L ) ∂η x ∂x (0 , L ) + ( u − u )(0 , − L ) ∂η x ∂x (0 , − L ) − ( u − u )( h, L ) ∂η x ∂x ( h, L ) + ( u − u )( h, − L ) ∂η x ∂x ( h, − L )+ η x (0 , L ) ∂ ( u − u ) ∂x (0 , L ) − η x (0 , − L ) ∂ ( u − u ) ∂x (0 , − L )+ η x ( h, L ) ∂ ( u − u ) ∂x ( h, L ) − η x ( h, − L ) ∂ ( u − u ) ∂x ( h, − L ) . (5.4)However, substituting (3.5) and (3.8) into the above expression, we can directly see that all the boundaryterms in the expression are zero, hence we arrive at the following from (2.18), (2.23) and the Green’sformula that h η x , f − f i = − Z Γ η x ∆ Γ ( f − f ) dσ Γ = − Z Γ ∂w x ∂ν ϕ dσ Γ = Z ∂ Ω (cid:18) w x ∂ϕ∂ν − ϕ ∂w x ∂ν (cid:19) dσ = ϕ ( x ) . (5.5)Therefore, we can see from both (5.3) and (5.5) that, in both cases we consider, we have the followingrepresentation of the duality product h · , · i of η x and f − f : h η x , f − f i = ϕ ( x ) . (5.6)Substituting this expression into (2.27), we derive the following formula for the index function I ( x ) = h η x , f − f i | η x | H (Γ) | η x | H (Γ) = ϕ ( x ) | η x | H (Γ) | η x | H (Γ) , x ∈ Ω . (5.7)For this expression, we suggest an appropriate scaling function S ( x ) = 1 || ϕ || L ∞ (Ω) + | ϕ ( x ) | , x ∈ Ω . (5.8)One may note from the calculations similar to (4.43)-(4.44) that the L (Γ) norm and H (Γ) semi-normof η x can be actually approximated by those of the fundamental solution Φ x defined as in (2.19)-(2.20),namely | η x | H (Γ) ≈ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Φ x ∂ν (cid:12)(cid:12)(cid:12)(cid:12) H (Γ) , | η x | L (Γ) ≈ (cid:12)(cid:12)(cid:12)(cid:12) ∂ Φ x ∂ν (cid:12)(cid:12)(cid:12)(cid:12) L (Γ) , ∀ x ∈ Ω . (5.9)With these approximations, we propose the following modified version of the index function in (5.7): e I ( x ) := 1 (cid:12)(cid:12) ∂ Φ x ∂ν (cid:12)(cid:12) H (Γ) (cid:12)(cid:12) ∂ Φ x ∂ν (cid:12)(cid:12) H (Γ) ϕ ( x ) || ϕ || L ∞ (Ω) + | ϕ ( x ) | , x ∈ Ω . (5.10)The numerical experiments in the next section will confirm that this modified index function improvesthe images in locating inhomogeneities in the DOT problem.In what follows, we would like to give a comparison between the original index function I and themodified index function e I . For this purpose, given z ∈ Ω, let µ z be an absorption coefficient of thefollowing form µ z ( x ) = µ + µ χ B ε ( z ) ( x ) , x ∈ Ω , (5.11)15here ε ∈ R is a very small number, µ ∈ R is a non-zero constant, B ε ( z ) is an open ball of radius ε andcenter z , and χ B ε ( z ) is its characteristic function. Then from (2.6) and (2.15), we can see that, if µ = 0,the scattered potential f along Γ satisfies( f − f )( ξ ) ≈ a z G z ( ξ ) , ξ ∈ ∂ Ω , (5.12)for some a z ∈ R ; whereas if µ = 0, we have( f − f )( ξ ) ≈ a z G z ( ξ ) − | ∂ Ω | Z ∂ Ω f dσ , ξ ∈ ∂ Ω . (5.13)In either case, we can choose a value of µ such that a z = 1. Then, from (2.11) and (2.17), we can seethat I ( x ) ≈ a z K ( x, z ) = K ( x, z ) , x ∈ Ω , (5.14)where K is the function defined as in (2.9). Now, given z ∈ Ω, with this choice of µ z , we can also calculatethe modified index function e I ( x ) in (5.10) from the the scattered potential f along Γ. We denote thisindex function e I ( x ) as e K ( · , z ), i.e., e K ( x, z ) := e I ( x ) , x ∈ Ω . (5.15)We can now compare the behaviors of the two functions K ( · , z ) and e K ( · , z ) to get a better understandingof the two index functions and compare their performances. Figure 3 shows the values of the function e K ( x, z ) with z = (0 . , − . x ∈ (0 , h ) × ( − L, L ). −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 0.20.40.60.81 Figure 3:
The function e K ( · , z ) on a rectangular domain with z = (0 . , − . Comparing Figure 3 with Figure 2, we observe the maximum point of the function K ( · , z ) in Figure 2is shifted a bit upward and away from the point z , whereas the maximum point of e K ( · , z ) in Figure3 matches perfectly with the point z . This phenomenon is of crucial importance, and accounts for thebetter performance and higher accuracy of this modified index in (5.10) to locate inhomogeneities for theDOT problem.Moreover, we would like to emphasize that the evaluation of the modified index function e I ( x ) definedas in (5.10) is computationally much cheaper than (2.27). The function ϕ in the expression (5.10) can befound by solving a forward elliptic system (5.1) once a priori for which fast solvers with nearly optimalcomplexity are available for general domains, such as multigrid methods and domain decomposition meth-ods, while the semi-norms of the fundamental solution Φ x can be efficiently approximated by numericalquadrature rules since the closed form of Φ x are explicitly known.Therefore this modified index function b I shall serve as a very efficient and robust method to locateinhomogeneities for the DOT problem. In this section, we shall present several numerical examples to illustrate the effectiveness of the newlyproposed DSM method for the reconstruction of inhomogeneous inclusions in the DOT problem.We consider the sampling region Ω = ( − , × (0 , . µ = 0 in Ω. In each of the following examples, there are some inhomogeneous inclusions16laced inside Ω, with their absorption coefficients always set to be µ = 50. Our choice of Γ is ( − , ×{ , . } . In order to collect our observed data of the forward problem, we solve (1.1) by second ordercentral finite difference method with a fine mesh of size 0 .
004 and the boundary flux g = 1 on ( − , ×{ . } , g = − − , × { } and g = 0 on ∂ Ω \ Γ. The boundary flux g is chosen as above so as tocreate a uniform flow from the bottom to the top of the domain. We have observed from our numericalexperiments that this choice of boundary flux is very effective in obtaining reasonable Cauchy data forour concerned imaging. The scattered potential f s := f − f is then measured along Γ. We would like toemphasize that, in each of our following examples, we only collect the scattered potential from a singlechoice of boundary flux, therefore only one pair of the Cauchy data is available for our reconstruction.So the resulting inverse problem is a severely ill-posed problem. In order to test the robustness of ourreconstruction algorithm, we introduce some random noise in the scattered potential as follows: f δs ( x ) = f s ( x ) + εδ max x | f s ( x ) | , (6.1)where δ is uniformly distributed between − ε corresponds to the noise level in the data, whichis always set to be ε = 5% in all our examples.From the noisy observed data f δs , we then use the DSM method to solve the DOT problem bycalculating the index function I introduced in (2.27) and the modified index function e I in (5.10). Fromour numerical experiments, we observe that the square of modified index e I gives us sharper images of theinclusions, hence providing a more accurate estimate of the support of the inhomogeneities D . Therefore,in each of the following examples, we provide the reconstructed images from all the three indices I , e I and e I . For the better illustrative purpose and comparison, we normalize each of the aforementionedindices with a normalizing constant such that the maximum of each index is 1. The mesh size in ourreconstruction process is chosen to be 0 . −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.4 01020304050 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 Figure 4:
From top to bottom: exact medium in Example 1; index function I ; modified index function e I ; squareof modified index function e I . Example 1 . This example tests a medium with two circular inclusions of radius 0 . − . , .
25) and (0 . , . I , the modified index function e I andthe square of e I . We can see from the figures that the recovered scatterers are well separated, and thelocations of the scatterers are accurately reconstructed. We may also observe that the artifacts due to thepresence of noise are effectively removed by the modified index function. Moreover, the squared modifiedindex function provides the sharpest image of the inclusions and gives most accurate locations of theinclusions as well as being the most robust to noise. Example 2 . In this example, our medium consists of 4 circular inclusions of radius 0 .
065 withtheir corresponding positions: ( − . , . − . , . . , .
1) and (0 . , . I , the modified index e I and its square e I are shown in Figure 5(second to bottom). We observe that both the locations and sizes of the reconstructed scatterers agreewell with those of the true ones. We can also see that the modified index function e I gives a betterestimate of the sizes and locations of the inclusions as well as removes all the artifacts from the indexfunction I . In addition, the square of modified index provides the best image, because it not only givesmore focused and accurate shapes and locations of the inclusions, but also significantly reduce the heavyshadows from the image. Considering the severe ill-posedness of the problem and the challenging caseof 4 scatterers staying very close to each other in the vertical direction and close to the boundaries, ourreconstruction seem to be quite satisfactory with only one pair of the Cauchy data. −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.4 01020304050 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 Figure 5:
From top to bottom: exact medium in Example 2; index function I ; modified index function e I ; squareof modified index function e I . Example 3 . A medium with 4 circular inclusions of radius 0 .
065 centered at ( − . , . − . , . , .
3) and (0 . , .
1) is investigated in this example; see Figure 6 (top). Figure 6 (second to bottom)displays the three reconstructed images from the index I , the modified index e I and its square e I respec-tively. Although we observe some minor shifting of the recovered scatterers in the images, they are stillwell separated and located with reasonable accuracy, and the artifacts for the index function I are allremoved by the modified index function e I . Furthermore, the squared modified index function e I givesthe best reconstruction by providing the most accurate estimate of locations and shapes of the scatterersas well as removing part of the heavy shadows from e I .18 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.4 01020304050 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 Figure 6:
From top to bottom: exact medium in Example 3; index function I ; modified index function e I ; squareof modified index function e I . Example 4 . This example investigates a rectangular inclusions of width 0 . . , . I , e I and e I . We can observe from the reconstructed images that the reconstructed inclusion ismildly shortened and shifted downwards, and that the modified indicator function e I removes the artifactsfrom the index function I . In addition, the squared index e I significantly removes the shadows from e I ,hence providing the most reliable reconstruction of the inclusion. Example 5 . In this last example, two rectangular inclusions of width 0 .
05 and length 0 . − . , . − . , . I , e I and e I are respectively shown in Figure 8 (second to bottom). Due to the ill-posed nature of theDOT, the reconstructed inclusions are not free from shadows and oscillations. However, the overallprofile stands out clearly, and both the locations and the shapes of the inclusions are well reconstructed.Furthermore, the shadows are greatly reduced by the squared modified index function e I .In all of the above examples, we may observe that the reconstructions seem to be rather satisfactory,considering the severe ill-posedness of the inverse problem, and the facts that the 5% random noise isadded in the observed data and only one pair of the Cauchy data is available. Moreover, the square e I ofthe modified index function provides the most accurate locations for the inclusions and is the most robustagainst noise. Considering the good accuracy of reconstruction, computational efficiency and robustnessto noise, the new DSM proves to be an effective new numerical tool for the reconstruction of the DOTproblem. 19 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.4 01020304050 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 Figure 7:
From top to bottom: exact medium in Example 4; index function I ; modified index function e I ; squareof modified index function e I . −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.4 01020304050 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.80.10.20.3 00.20.40.60.81 Figure 8:
From top to bottom: exact medium in Example 5; index function I ; modified index function e I ; squareof modified index function e I . Concluding remarks
In this work we have presented a novel direct sampling method to locate inhomogeneities inside ahomogeneous background for the diffusive optical tomography problem, applicable even to one scatteredpotential data and in both full and limited aperture cases. It involves computing only a well-posed ellipticsystem involving the scattered potential as its boundary condition, therefore it can be viewed as a directmethod. This method is very easy to implement and computationally inexpensive since it does not requireany matrix operations solving ill-posed integral equations or any optimisation process as most existingalgorithms do. In particular, numerical experiments have demonstrated the effectiveness and robustnessof the proposed method. This new method provides an efficient numerical strategy and a new promisingdirection for solving the DOT problem.
References [1] R. Aarbour, H. Graber, J. Chang, S. Barbour, S. Koo and R. Aroson,
MRI-guided optical tomography ,IEEE Comp. Sci. Eng. 2, pp. 63-77 (1995).[2] S. Arridge and M. Schweiger,
Sensitivity to prior knowledge in optical tomographic reconstruction ,Proc. SPIE, 2389, pp. 378-388 (1995).[3] D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang,
Imaging the body with diffuse optical tomography , IEEE Sig. Proc. Mag. 18, pp. 57-75 (2001).[4] W. Cai, B.B. Das, F. Liu; F. Zeng; M. Lax, R. R. Alfano,
Three-dimensional image reconstructionin highly scattering turbid media , Proc. SPIE Vol. 2979, Optical Tomography and Spectroscopy ofTissue: Theory, Instrumentation, Model, and Human Studies II, 241, pp. 240-247 (1997).[5] W. Cai, D. Das, F. Lui, M. Zevallos, M. Lax and R. Alfano,
Time-resolved optical diffusion tomo-graphic image reconstruction in highly scattering turbid media , Proc. Natl. Acad. Sci. USA, 93(24),pp. 13561-13564 (1996).[6] R. Choe, Diffuse optical tomography and spectroscopy of breast cancer and fetal brain, Ph.D. Diss.University of Pennsylvania, (2005).[7] J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B.Chance, and A. G. Yodh,
Three-dimensional diffuse optical tomography in the parallelplane trans-mission geometry: Evaluation of a hybrid frequency domain/continuous wave clinical system forbreast imaging , Med. Phys. 30, pp. 235-247 (2003).[8] J. P. Culver, T. Durduran, T. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh,
Diffuse opticaltomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia , J. Cereb.Blood Flow Metab. 23, pp. 911-924 (2003).[9] H. Dehghani, B.W. Pogue, S. D. Jiang, B. Brooksby, and K. D. Paulsen,
Three-dimensional opticaltomography: resolution in small-object imaging , Appl. Opt. 42, pp. 3117-3128 (2003).[10] T. Durduran,
Non-invasive measurements of tissue hemodynamics with hybrid diffuse optical meth-ods , Ph.D. Diss. University of Pennsylvania (2004).[11] M. A. Franceschini, K. T. Moesta, S. Fantini, G. Gaida, E. Gratton, H. Jess,W.W. Mantulin, M.Seeber, P. M. Schlag, and M. Kaschke,
Frequency-domain techniques enhance optical mammography:Initial clinical results , Proc. Natl. Acad. Sci. U. S. A. 94, pp. 6468-6473 (1997).[12] D. Grosenick, H. Wabnitz, K. T. Moesta, J. Mucke, M. Moller, C. Stroszczynski, J. Stossel, B.Wassermann, P. M. Schlag, and H. Rinneberg,
Concentration and oxygen saturation of haemoglobinof 50 breast tumours determined by time-domain optical mammography , Phys. Med. Biol. 49(7), pp.1165-1181 (2004). 2113] X. J. Gu, Q. Z. Zhang, M. Bartlett, L. Schutz, L. L. Fajardo, and H. B. Jiang,
Differentiation ofcysts from solid tumors in the breast with diffuse optical tomography , Acad. Radiol. 11, pp. 53-60(2004).[14] K. Ito, B. Jin and J. Zou,
A direct sampling method to an inverse medium scattering problem , Inv.Prob., 28(2012), 025003 (11pp).[15] K. Ito, B. Jin and J. Zou,
A direct sampling method for inverse electromagnetic medium scattering ,Inv. Prob., 29 (2013) 095018 (19pp).[16] A.H. Hielschera, A.Y. Bluestonea, G.S. Abdoulaeva, A.D. Klosea, J. Laskera, M. Stewartb, U. Netzcand J. Beuthanc
Near-infrared diffuse optical tomography , Disease Markers, 18(5), pp. 313-337(2002).[17] D. M. Hueber, M. A. Franceschini, H. Y. Ma, Q. Zhang, J. R. Ballesteros, S. Fantini, D. Wallace, V.Ntziachristos, and B. Chance,
Non-invasive and quantitative near-infrared haemoglobin spectrometryin the piglet brain during hypoxic stress, using a frequency domain multi-distance instrument , Phys.Med. Biol. 46, pp. 41-62 (2001).[18] H. Jiang, Y. Xu, N. Ifitimia, J. Eggert, K. Klove, L. Baron, and L. Fajardo,
Three dimensionaloptical tomographic imaging of breast in a human subject , IEEE Trans. Med. Imaging 20, pp. 1334-1340 (2001).[19] J. R. Kappa, H. D. Berkowitz, R. Seestedt, and B. Chance,
Evaluation of calf muscle oxygen content.using NIRS in patients with peripheral vascular disease , Medical intelligence unit: vascular surgery2000, pp. 10-16 (1991).[20] M.V. Klibanov, T. R. Lucas, and R. M. Frank.,
A fast and accurate imaging algorithm in opti-cal/diffusion tomography , Inv. Prob., 13(1997), pp. 1341-1361.[21] M.V. Klibanov, T. R. Lucas, and R. M. Frank and M Robert,
New imaging algorithm in diffusiontomography , Optical Tomography and Spectroscopy of Tissue: Theory, Instrumentation, Model, andHuman Studies II, Britton Chance, Proc. SPIE Vol. 2979, pp. 272-283 (1997).[22] A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M.Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas,
Tomographic optical breast imaging guidedby three-dimensional mammography , Appl. Opt. 42, pp. 5181-5190 (2003).[23] Q. Luo, S. Nioka and B. Chance,
Functional near-infrared imager , Optical Tomography and Spec-troscopy of Tissue: Theory, Instrumentation, Model, and Human Studies II, Proc. SPIE 2979, pp.84-93 (1997).[24] J.S. Maier, A.E. Cerussi, S. Fantini, and E. Gratton,
Experimental recovery of absorption, scattering,and fluorescence parameters in highly-scattering media from a single frequency measurement , Proc.SPIE 2979, pp. 6-13 (1997).[25] V.A. Markel, Z-M Wang, J.C. Schotland,
Optical Diffusion Tomography With Large Data Sets ,Photonic Applications in Biosensing and Imaging, SPIE 5969, pp. 59691B-59691B-11 (2005).[26] V. Ntziachristos and B. Chance.,
Probing physiology and molecular function using optical imaging:applications to breast cancer , Breast Cancer Res. 3, pp. 41-46 (2001).[27] V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance,
Concurrent MRI and diffuse opticaltomography of breast after indocyanine green enhancement , Proc. Natl. Acad. Sci. U. S. A. 97, pp.2767-2772 (2000). 2228] S. Oh, A.B. Milstein, C.A. Bouman and K.J. Webb,
Multigrid inversion algorithms with applica-tions to optical diffusion tomography , Proc. of 36th Asilomar Conference on Signals, Systems, andComputers, pp. 901-905 (2002).[29] M. O’Leary, D. Boas, B. Chance and A. Yodh,
Images of inohomogeneous turbid media using diffusephoton density waves , Optical Society of America Proc. ‘Advances In Optical Imaging And PhotonMigration’, 21, pp. 106-115 (1994).[30] N. Pantong, J. Su, H. Shan, M.V. Klibanov and H. Liu,
Globally accelerated reconstruction algorithmfor diffusion tomography with continuous-wave source in an arbitrary convex shape domain , J Opt.Soc. Am. A Opt. Image Sci. Vis. 26(3), pp.456-472 (2009).[31] Y. Pei, Y. Yao, F-B Lin and R.L. Barbour,
Frequency domain fluorescence optical imaging using afinite element method , Optical Tomography and Spectroscopy of Tissue: Theory, Instrumentation,Model, and Human Studies II, Proc. SPIE 2979, pp. 164-175 (1997).[32] A. Pifferi, P. Taroni, A. Torricelli, F. Messina, R. Cubeddu, and G. Danesini,
Fourwavelength time-resolved optical mammography in the 680-980-nm range
Opt. Lett. 28, pp. 1138-1140 (2003).[33] B. W. Pogue, S. P. Poplack, T. O. McBride, W. A. Wells, K. S. Osterman, U. L. Osterberg, andK. D. Paulsen,
Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: Pilotresults in the breast , Radiology 218, pp. 261-266 (2001).[34] S. Takahashi, D. Imai, Y. Tanikawa and Y. Yamada,
Fundamental 3D FEM analysis of light propaga-tion in head model toward 3D optical tomography , Optical Tomography and Spectroscopy of Tissue:Theory, Instrumentation, Model, and Human Studies II, Proc. SPIE 2979, pp. 250-260 (1997).[35] T.L. Troy, J.S. Reynolds and E.M. Sevick-Muraca,
Photon migration imaging using multipixel mea-surements , Optical Tomography and Spectroscopy of Tissue: Theory, Instrumentation, Model, andHuman Studies II, Proc. SPIE 2979, pp. 111-121 (1997).[36] U. Sunar, J. Zhang, J. Du, T. Durduran, C. Zhou, G. Yu, A. Kilger, H. Quon, R. Lustig, L. Loevner,S. Nioka, A.G. Yodh, and B. Chance,
Clinical responses of head and neck tumors to radiation therapyby NIR spectroscopy , Biomedical Topical Meetings (The Optical Society of America, Washington,DC) FB7., OSA, (2004).[37] H.Wang, M. E. Putt, M. J. Emanuele, D. B. Shin, E. Glatstein, A. G. Yodh, and T. M. Busch.
Treatment-induced changes in tumor oxygenation predict photodynamic therapy outcome , CancerRes. 64, pp. 7553-7561 (2004).[38] B. C.Wilson, M. S. Patterson, and L. Lilge.,
Implicit and explicit dosimetry in photodynamic therapy:a new paradigm , Lasers Med. Sci. 12(3), pp. 182-199 (1997).[39] U. Wolf, M. Wolf, J. H. Choi, L. A. Paunescu, L. P. Safonova, A. Michalos, and E. Gratton,
Mappingof hemodynamics on the human calf with near infrared spectroscopy and the influence of the adiposetissue thickness , Adv. Exp. Med. Biol. 510, pp. 225-230 (2003).[40] Y. Yao, Y. Pei ; Y. Wang and R.L. Barbour,
Born type iterative method for imaging of heterogeneousscattering media and its application to simulated breast tissue , Optical Tomography and Spectroscopyof Tissue: Theory, Instrumentation, Model, and Human Studies II, Proc. SPIE, 2979, pp. 232-240(1997).[41] J.C. Ye, C.A. Bouman, K.J. Webb and R. P. Millane,
Nonlinear multigrid algorithms for Bayesianoptical diffusion tomography , IEEE Transections on Image Processing, 10(5), pp. 909-922 (2001).[42] Q.I. Zhu, M. M. Huang, N. G. Chen, K. Zarfos, B. Jagjivan, M. Kane, P. Hedge, and S. H. Kurtzman,