Optical imaging of phantoms from real data by an approximately globally convergent inverse algorithm

Abstract

A numerical method for an inverse problem for an elliptic equation with the running source at multiple positions is presented. This algorithm does not rely on a good first guess for the solution. The so-called ‘approximate global convergence’ property of this method is shown here. The performance of the algorithm is verified on real data for Diffusion Optical Tomography. Direct applications are in near-infrared laser imaging technology for stroke detection in brains of small animals.

Keywords:

• approximate global convergence property
• inverse problem
• diffusion optical tomography
• real data

AMS Subject Classifications::

• 65B21
• 65D10
• 65F10

1. Introduction

We consider a Coefficient Inverse Problem (CIP) for a partial differential equation (PDE) – the diffusion model with the unknown potential. The boundary data for this CIP, which model measurements, are originated by a point source running along a part of a straight line. This PDE governs light propagation in a diffusive medium, such as, e.g. biological tissue, smog. Thus, our CIP is one of problems of Diffusion Optical Tomography (DOT). We are interested in applications of DOT to the detection of stroke in small animals using measurements of near infrared light originated by lasers. Hence, the above point source is the light source in our case. The motivation of imaging of small animals comes from the idea that it might be a model case for the future stroke detection in humans, at least in brains of neonatals, via DOT. We apply our numerical method to a set of real data for a phantom medium modelling the mouse's brain. Although this algorithm was developed in earlier publications 1–4 of this group, its experimental verification is new and is the main contribution of this article.

As to our numerical method, we introduce a new concept of the ‘approximate global convergence’ property. In the previous publication 1 of this group about this method, the approximate global convergence property was established in the continuous case. Compared with 1, the main new analytical result here is that, we establish this property for the more realistic discrete case. In the convergence analysis of 1 the Schauder theorem 5 was applied for C^2+α – solutions of certain elliptic equations arising in our method. Now, however, since we consider the discrete case, we use the Lax–Milgram theorem. Here and below C^k+α are Hölder spaces 5, where k ≥ 0 is an integer and α ∈ (0, 1).

CIPs are both nonlinear and ill-posed. These two factors cause very serious challenges in their numerical treatments. Indeed, corresponding least squares Tikhonov regularization functionals usually suffer from multiple local minima and ravines. As a result, conventional numerical methods for CIPs are locally convergent ones, see, e.g. 6 and references cited there. To have a guaranteed convergence to a true solution, a locally convergent algorithm should start from a point which is located in a small neighbourhood of this solution. However, it is a rare case in applications when such a point is known. The main reason why our method avoids local minima is that it uses the structure of the underlying PDE operator and does not use a least squares functional.

Therefore, it is important to develop such numerical methods for CIPs, which would have a rigorous guarantee of providing a point in a small neighbourhood of that exact solution without any a priori knowledge of this neighbourhood. It is well known, however, that the goal of the development of such algorithms is a substantially challenging one. Hence, it is unlikely that such numerical methods would not rely on some approximations. This is the reason why the notion of the approximate global convergence property was introduced in the recent work 7; also see Section 1.1.2 of the book 8 and Subsection 4.1 below. The verification of our approximately globally convergent numerical method on computationally simulated data was done in 1–4. In this article, we make the next step: we verify this method on real data. Regardless on some approximations we have here, the main point is that our numerical method does not rely on any a priori knowledge of a small neighbourhood of the exact solution.

In 7,9–13 a similar numerical method for CIPs for a hyperbolic PDE was developed. These results were summarized in the book 8. In particular, it was demonstrated numerically in Test 5 of 9, Section 8.4 of 12, on page 185 of 8 and in Section 5.8.4 of 8 that the approximately globally convergent method of 7,9–13 outperforms such locally convergent algorithms as quasi-Newton method and gradient method. The main difference of the technique of these publications with the one of the current paper is in the truncation of a certain integral. In 7,9–13, the integral was truncated at a high value of the parameter s > 0 of the Laplace transform of the solution of the underlying hyperbolic PDE. Indeed, in the hyperbolic case the truncated residual of that integral, which we call the ‘tail function’, is , as , i.e. this residual is automatically small in this case. Being different from the hyperbolic case, in the elliptic PDE with the source running along a straight line, represents the distance between the source and the origin. Because of this, the tail function is not automatically small in our case. Thus, a special effort to ensure this smallness was undertaken in 1, and is repeated in the current publication. It was shown numerically in 1 that this special effort indeed improves the quality of the reconstruction, compare Figure 2(b) and (c) in 1.

Our experimental data were collected at the boundary of a 2-D cross-section of the 3-D domain of interest (Section 6). Therefore, we should work only with the 2-D mathematical model. The reason is that limitations of our device do not allow us to arrange the data collection, which is necessary for the use of the 3-D mathematical model. Hence, we have imaged this 2-D cross-section only and have ignored the dependence on the third variable. Also, our theory requires the use of many sources placed along a straight line. However, the configuration of our device allows us to use only three sources on this line. In addition to experimental limitations, if using the 3-D model, we would face a serious challenge of approximating the tail function. More precisely, we do not know yet how to produce 3-D analogs of formulas (47), (48) (Section 5.2) for the tail function. Likewise, while the current procedure of approximating tail functions is well tested computationally in the 2-D case 1, it is yet unclear, how a 3-D analog of this procedure would work numerically. As to the convergence analysis in 3-D, it is very similar to the 2-D case. Regardless on these discrepancies, we believe that the accuracy of our results (Table 8.2) provides a good confirmation of the validity of our assumptions as well as manifests a well known fact that there are always some mismatches between the analysis and computational studies of numerical methods, especially in the cases of real data.

Because of the above-mentioned substantial challenges, the topic of the development of non-locally convergent numerical methods for CIPs is currently in its infancy. As to those methods for CIPs for elliptic PDEs, we refer to, e.g. publications 14–19 and references cited there. Papers 14,15,18 are concerned with the scattering amplitudes and publications 16,17,19 are concerned with the Dirichlet-to-Neumann map (DtN). We are not using neither the scattering amplitude nor the DtN here.

In a general setting, an ill-posed problem is a problem of solving an equation F(x) = y with a compact operator F. Here x ∈ B₁, y ∈ B₂, where B₁ and B₂ are certain Banach spaces. Naturally, the operator F varies from one problem to another one. It is well known that the existence theorem for such an equation is very tough to prove, since the range of any compact operator is ‘too narrow’ in certain sense, see, e.g. Theorem 1.2 in 8. Therefore, by one of fundamental concepts of the theory of Ill-Posed Problems, one should assume the existence of an exact solution of such a problem for the case of an ‘ideal’ noiseless data 8,20,21. Although this solution is never known in practice and is never achieved in practice (because of the noise in the real data), the regularization theory says that one needs to construct a good approximation for it. We assume that our CIP has an exact solution for noiseless data and also assume that this solution is unique.

The rest of this article is arranged as follows. In Section 2, we pose both forward and inverse problems and study some properties of the solution of the forward problem. In Section 3, we present our numerical method. In Section 4, we conduct the convergence analysis. In Section 5, we discuss the numerical implementation of our method. In Section 6, we describe the experiment. Our data processing is described in Section 7. In Section 8, we present reconstruction results. We briefly summarize the results in Section 9.

2. Statement of the problem

2.1. The inverse problem

Below x = (x, z) ∈ ℝ² and Ω ⊂ ℝ² is a convex bounded domain. In our convergence analysis we assume that ∂Ω ∈ C³. In numerical studies ∂Ω is piecewise smooth. Consider the following elliptic equation in ℝ² with the solution vanishing at infinity,(1) (2)

Inverse problem

Let k = const. > 0. Suppose that in (1) the coefficient a(x) satisfies the following conditions(3) Let be a straight line and Γ ⊂ L be an unbounded and connected subset of L. Determine the function a(x) inside of the domain Ω, assuming that the constant k is given and also that the following function ϕ(x, x₀) is given(4)

The authors are unaware about a uniqueness theorem for this inverse problem. Nevertheless, because of the above application, it is reasonable to develop a numerical method. Thus, we assume that uniqueness for this problem holds. Our numerical studies of both the past 1–4 and the current publication indicate that a certain uniqueness theorem might be established.

We assume that sources {x₀} are located outside of the domain of interest Ω, because this is the case of our measurements and because we do not want to work with singularities in our numerical method. The CIP (1)–(4) has an application in imaging using light propagation in a diffuse medium, such as biological tissues. Since the modulated frequency equals zero in our case, then this is the so-called continuous-wave (CW) light. The coefficient where is the reduced scattering coefficient and μ_a(x) is the absorption coefficient of the medium 6,22. In the case of our particular interest in stroke detections in brains of small animals, the area of an early stroke can be modelled as a small sharp inclusion in an otherwise slowly fluctuating background. Therefore, our focus is on the reconstruction of both locations of sharp small inclusions and the values of the coefficient a(x) inside of them, rather than on the reconstruction of slow changing background functions.

2.2. Some properties of the solution of the forward problem (1), (2)

2.1.1. Existence and uniqueness

First, we state the existence and uniqueness of the solution of the forward problem (1), (2). For brevity we consider only the case since this is the case of our Inverse Problem. Let K_p(z), z ∈ ℝ, p ≥ 0 be the Macdonald function. It is well known 23 that for y ∈ ℝ(5)

Theorem 2.1

Let Ω ⊂ ℝ² be the above bounded domain. Assume that the coefficient a(x) satisfies conditions (3). Then for each source position there exists unique solution u(x, x₀) of the problem (1), (2) such that(6) where the function u₀ is the fundamental solution of Equation (1) with a(x) ≡ k², the function satisfies (2), and In addition,

A similar result for the 3-D case was proven in 8, see Theorem 2.7.2 in this reference. Hence, we leave out the proof of Theorem 2.1 for brevity.

2.2.2. The asymptotic behaviour at |x₀| → ∞

It follows from (5) and (6) that the asymptotic behaviour of the function u₀(x − x₀) isDenote b(x) = a(x) − k². Then by (3) b(x) = 0 for x ∈ ℝ²∖Ω. Let M₁ be a positive constant. DenoteAlso, let the function p_∞(x) satisfies conditions(7) and be the solution of the following problem(8) (9) The uniqueness and existence of the solution of the problem (7)–(9) are similar to these of Theorem 2.1.

Lemma 2.2

1 + p_∞ > 0.

Proof

Let Then(10) Consider a sufficiently large number R > 0 such that for x ∈ {|x| ≥ R}. Then the maximum principle 5 applied to equation (10) for x ∈ {|x| < R} shows that in {|x| < R}.▪

Lemma 2.3 1

Let the function b ∈ B(M₁). Then there exists a constant M₂(M₁, Ω) > 0 such that

3. The numerical method

Since we can put the origin on the straight line L, if necessary, then without any loss of the generality we can set s ≔ |x₀|, assuming that only the parameter s changes when the source x₀ runs along Γ ⊂ L. Denote u(x, s) ≔ u(x, x₀), x ∈ Ω, x₀ ∈ Γ. Since and the point source x₀ ∈ Γ, then . Since by Theorem 2.1, , then we can consider the function for x ∈ Ω. We obtain from (1) and (4)(11) (12) where ϕ₁ = ln ϕ and are two positive numbers, which should be chosen in numerical experiments.

3.1. The integral differential equation

We now ‘eliminate’ the coefficient a(x) from Equation (11) via the differentiation with respect to s. However, to make sure that the resulting the so-called ‘tail function’ is small, we use the above mentioned (Introduction) special effort of the paper 1. Namely, we divide (11) by s^p, p > 0. In principle, any number p > 0 can be used. But since in computations we took p = 2, both here and in 1, then we use below only p = 2, for the sake of definiteness. Denote w(x, s) = . By (11), equation for the function w(x, s) is(13) Next, let q(x, s) ≔ ∂_sw(x, s) = ∂_s(s⁻² ln u(x, s)), for . Then,(14) where is the so-called ‘tail function’. The exact expression for this function is of course . However, since the function is unknown, we will use an approximation for the tail function, see Subsection 5.2 as well as 1. By the Tikhonov concept for ill-posed problems 8,20,21, one should have some a priori information about the solution of an ill-posed problem. Thus, we can assume the knowledge of a constant M₁ > 0 such that the function a(x) − k² = b(x) ∈ B(M₁). Hence, it follows from Lemma 2.3 that 1(15) where the number M₂ (M₁, Ω) is independent on Differentiate (13) with respect to s. Using (14), we obtain(16) Next, by (13) and (14)(17) Comparing (17) with (16), we obtain(18) In addition, (12) implies that(19) (20)

If we approximate both functions q and T exactly, together with their derivatives up to the second order, then we can also approximate well the target coefficient a(x) via backwards calculations (43). Therefore, the main questions now is: How to approximate well both functions q and T using (14)–(20)?

3.2. Layer stripping with respect to the source position

In this subsection, we present a layer stripping procedure with respect to s for approximating the function q, assuming that the function T is known. Usually the layer stripping procedure is applied with respect to a spatial variable. However, sometimes the presence of a differential operator with respect to this variable in the underlying PDE results in the instability, since computing derivatives is an unstable procedure. The reason why our layer stripping procedure is stable is that we do not have a differential operator with respect to s in our PDE.

We approximate the function q(x, s) as a piecewise constant function with respect to the source position s. We assume that there exists a partition of the interval with a sufficiently small step size h such that(21) Hence,Let ψ_n(x) be the average of the function ψ(x) over the interval (s_n, s_n−1). Then, we approximate the boundary condition (19) as a piecewise constant function with respect to s,(22) Using (21), integrate Equation (18) with respect to s ∈ [s_n, s_n−1). We obtain for n ≥ 1(23) where A_k,n, k = 1, … , 4 are certain coefficients, see 1 for exact formulas for them.

Let(24) Then one can prove that (25) Assuming that h is sufficiently small and using (25), we assume below that(26)

It is convenient for our convergence analysis to formulate the Dirichlet boundary value problem (22), (23) in the weak form. Denote . Assume that there exists functions Ψ_n such that(27) Consider the function Then, we obtain an obvious analogue of Equation (23) for the function p_n. Multiply both sides of the latter equation by an arbitrary function and integrate over Ω using integration by parts and (26). We obtain(28) Hence, the function q_n ∈ H¹ (Ω) is a weak solution of the problem (22), (23) if and only if the function satisfies the integral identity (28). The question about existence and uniqueness of the weak solution of the problem (28) is addressed in Theorem 4.2.

We now, describe our algorithm of sequential solutions of boundary value problems (22), (23) for n = 1, … , N, assuming that an approximation for the tail function is found (see Subsection 5.2 for the latter). Recall that q₀ = 0. Hence, we have:

Step n ∈ [1, N]. Suppose that functions q₁, … , q_n−1 are computed. On this step, we find the weak solution of the Dirichlet boundary problem (22), (23) for the function q_n via the FEM with triangular finite elements.

Step N + 1. After functions q₁, … , q_N are computed, the function a(x) is reconstructed using backwards calculations at n ≔ N, i.e. for the lowest value of s ≔ s_N = s as described in Subsection 5.1.

4. Convergence analysis

4.1. Approximate global convergence

The central question which was addressed in above cited publications 1–4 was of constructing such a numerical method for the CIP (1)–(4), which would simultaneously satisfy the following three conditions.

1.	This method should deliver a good approximation for the exact solution of this CIP without an a priori knowledge of a small neighbourhood of this solution.
2.	A theorem should be proven which would guarantee of obtaining such an approximation.
3.	A good numerical performance of this technique should be demonstrated on computationally simulated data and, optionally, on real data.

We use the word ‘optionally’, because it is usually hard and expensive to actually collect real data. It is enormously challenging to construct such a numerical method. Therefore, as it is often done in mathematical modelling, conceptually, our approach consists of the following six steps.

	Step 1. A reasonable approximate mathematical model is proposed. The accuracy of this model cannot be rigorously estimated.
	Step 2. A numerical method is developed, which works within the framework of this model.
	Step 3. A theorem is proven, which guarantees that, within the framework of this model, the numerical method of Step 2 indeed reaches a sufficiently small neighbourhood of the exact solution. Naturally, the smallness of this neighbourhood depends on the level of the error, both in the data and in some additional approximations is sufficiently small.
	Step 4. The numerical method of Step 2 is tested on computationally simulated data.
	Step 5. The numerical method of Step 2 is tested on real data.
	Step 6. Finally, if results of Steps 4 and 5 are good ones, then that approximate mathematical model is proclaimed as a valid one.

Thus, results of the current paper (Step 5) in combination with the previous ones of 1–4 lead to the positive conclusion of Step 6. Step 6 is logical, because its condition is that the resulting numerical method is proved to be effective. It is sufficient to achieve that small neighbourhood of the exact solution after a finite (rather than infinite) number of iterations. We refer to page 157 of the book 24 where it is stated that the number of iterations can be regarded as a regularization parameter sometimes for an ill-posed problem. Still, in our computations both for simulated and real data, classical convergence in the Cauchy sense was achieved. These consideration lead to the following definition of the approximate global convergence property.

Definition 4.1 (approximate global convergence) 7,8

Consider a nonlinear ill-posed problem P. Suppose that this problem has a unique solution x* ∈ B for the noiseless data y*, where B is a Banach space with the norm ‖ · ‖_B. We call x* ‘exact solution’ or ‘correct solution’. Suppose that a certain approximate mathematical model is proposed to solve the problem P numerically. Assume that, within the framework of the model this problem has unique exact solution Also, let one of the assumptions of the model be that Consider an iterative numerical method for solving the problem P. Suppose that this method produces a sequence of points where N ∈ [1, ∞). Let the number θ ∈ (0, 1). We call this numerical method approximately globally convergent of the level θ, or shortly globally convergent, if, within the framework of the approximate model a theorem is proven, which guarantees that, without any a priori knowledge of a sufficiently small neighbourhood of x*, there exists a number such that(29) Suppose that iterations are stopped at a certain number Then the point x_k is denoted as x_k ≔ x_glob and is called ‘the approximate solution resulting from this method’.

With reference to the notion of the approximate mathematical model , as well as to the above Steps 1–6, it is worthy to mention here that one of the keys to the successful numerical implementation 14 of the non-local reconstruction algorithm of 18 was the use of a certain approximate mathematical model. The same is true for the two dimensional analogue of the Gel'fand–Levitan–Krein equation being applied in 25 to an inverse hyperbolic problem. Thus, it seems that whenever one is trying to construct an efficient non-locally convergent numerical method for a truly complicated nonlinear ill-posed problem, it is close to the necessity to introduce some approximations which cannot be rigorously justified and then work within the resulting approximate mathematical model then. Analogously, although the Huygens–Fresnel theory of optics is not rigorously supported by the Maxwell's system (see Section 8.1 of the classical textbook 26), the ‘diffraction part’ of the entire modern optical industry is based on the Huygens–Fresnel optics.

4.2. Exact solution

By one of the statements of Section 1 as well as by Definition 4.1, we should assume the existence and uniqueness of an ‘ideal’ exact solution a*(x) of our inverse problem for an ‘ideal’ noiseless exact data ϕ*(x, x₀) in (4). Next, in accordance with the regularization theory, one should assume the presence of an error in the data of a small level γ > 0 and construct an approximate solution for this case.

Since the exact solution was defined in 1, we outline it only briefly here for the convenience of the reader. Let the function a*(x) satisfying conditions (3) be the exact solution of our inverse problem for the noiseless data ϕ*(x, x₀) in (4). We assume that a*(x) is unique. Let the function u*(x, s) be the same as the function u(x, x₀) in Theorem 2.1, but for the case a(x) ≔ a*(x). DenoteThe function q* satisfies the analogue of equation (18) with the boundary condition (19)(30) where by (20) ψ* (x, s) = ∂_s(s⁻² ln u* (x, s)) for . We call q*(x, s) the exact solution. It should be noted that our real data ϕ(x,s) are naturally given with a random noise. It is well known that the differentiation of the noisy data is an ill-posed problem 8,21. However, because of some limitations of our device, we work only with three values of the parameter s in our real data (Subsection 5.2). Hence, we simply calculate the s – derivative via the finite difference. Our numerical experience shows that this does not lead to a degradation of our reconstruction results.

4.3. Our approximate mathematical model

In this model, we basically assume that the upper bound for the tail function can become sufficiently small independently on In addition, since we calculate the tail function separately from the function q, we assume that the function is given (but not the exact tail function Although the smallness of T is supported by (15), the independence of that upper bound of this function from does not follow from (15). Still this is one of two assumptions of our approximate mathematical model. Following the concept of Steps 1–6 of Subsection 4.1, we now verify this model on real data.

More precisely, our approximate mathematical model consists of the following two assumptions

Assumptions

1.	We assume that the number is sufficiently large and fixed. Also, the tail function is given, and tail functions have the forms(31) (32) (33) where ξ ∈ (0, 1) is a small parameter. Furthermore, the parameter ξ is independent on
2.	As the unknown vector x ∈ B in Definition 4.1, we choose the vector function in (21), (23). In this case, we will have only one iteration in Definition 4.1, i.e. in (29).

It is for the sake of our convergence analysis that we choose here the vector function rather than the coefficient a (x) as our unknown function. Indeed, finding a(x) in the weak form (43), (44) would require some interpolation estimates of the differences between functions and their interpolations via the finite dimensional subspace G_m (see Subsection 4.4 for ). The latter has an ‘underwater rock’ in terms of the convergence analysis.

By Definition 4.1, we need to prove uniqueness of the function q*(x, s) for Assuming that the exact tail function satisfying (31), (33) is given and the interval is sufficiently small, this can be done easily, using the fact that q*(x, s) is the solution of an obvious analog of the problem (18), (19). This is because Volterra-like integrals are involved in (18).

4.4. Approximate global convergence theorem

Let γ be the sum of the following errors: the level of the error in the function ψ(x, s) (by (20) this function is generated by the measurement data), errors in our finite element approximations of q_n, the magnitude of the norm of the tail function as well as the length of the interval .

Consider a finite dimensional subspace with dim G_m = m. A realistic example of G_m is the subspace, generated by triangular finite elements, which are piecewise linear functions. We assume that f_x, f_z ∈ L_∞ (Ω), ∀f ∈ G_m. Since all norms in the subspace G_m are equivalent, then there exists a constant C_m = C_m(G_m) such that(34) Here and below and the same for the L₂(Ω) norm. Theorem 4.2 works with the weak formulation (28). In the course of the proof of this theorem, we need to estimate products ▿q_j▿p_n. While this was done for solutions in 1 using Schauder theorem, in the weak formulation, we need to assume that functions(35) We introduce functions which are analogues of functions q_n, Ψ_n for the case a ≔ a*. Let We assume that(36) (37) (38) (39) (40) where the number C* > 1 is given and σ > 0 is a small parameter characterizing the level of the error in the data ψ(x, s). The connection between σ and the error in ψ(x, s) is clear from comparison of (22), (27) with (40). Condition (38) is a direct analogue of (26). Each function is the weak solution of an analog of the problem (28).

Theorem 4.2

Let Ω ⊂ ℝ² be a convex bounded domain with the boundary ∂Ω ∈ C³. Assume that Assumptions 1 and 2 of Subsection 4.3 hold. In addition, assume that conditions (24), (26), (35)–(40) hold. Let . Consider the error parameter γ = β + h + σ + ξ, where ξ is the small parameter defined in (31)–(33), which is independent on (Assumption 1). Let C_m ≥ 1 be the constant defined in (34). Then there exists a constant B = B(Ω) > 0 such that if(41) then for each n ∈ [1, N] there exists unique solution p_n ∈ G_m of the problem (28) and for functions q_n = p_n + Ψ_n the following estimate holds(42) where Thus, (42) implies that our numerical method has the approximate global convergence property of the level θ (Definition 4.1).

Remarks 4.1

1.	We omit the proof of this theorem, because it is similar with the proof of a similar theorem in 1 for the Cα space. The only technical challenge in theorem 4.2 is to obtain analogous arguments in the finite dimensional space Gm.
2.	The smallness of the parameter β in (41) is a natural requirement since the original Equation (18) contains Volterra-like integrals in nonlinear terms. It is well known from the standard ODE course that the existence of a Volterra-like nonlinear integral equation of the second kind can be proven only on a small interval.

5. Numerical implementation

Since this implementation was described in detail in 1, we outline it only briefly here for the convenience of the reader. As to the functions q_n, we have sequentially calculated them via the FEM solving Dirichlet boundary value problems (22), (23). As it was mentioned in the end of Subsection 3.2, we have used standard triangular finite elements. Two important questions which are discussed in this section are about approximating the function a(x) and the tail function T(x).

5.1. Approximation of the function a(x)

We reconstruct the target coefficient a(x) via backwards calculations as follows. First, we reconstruct the function w(x, s_N) from q_n, n = 1, 2, … , N and tail T. We have u(x, s_N) = exp [s²w(x, s_N)]. Next, since in (1) the source we use Equation (1) in the weak form as(43) where the test function η_p(x), p ∈ [1, P] is a quadratic finite element of a computational mesh with the boundary condition η_p(x) |_∂Ω = 0. The number P is finite and depends on the mesh we choose. Equalities (43) lead to a linear algebraic system which we solve. Since this formulation is complex but standard, it is omitted here. Interested readers can see 27. Finally, we let(44)

5.2. Construction of the tail function

The above construction of functions q_n depends on the tail function T. In this subsection, we state briefly our procedure of approximating the tail function, see 1–4 for details. This procedure consists of two stages. First, we find a first guess for the tail using the asymptotic behaviour of the solution of the problem (1), (2) as |x₀| → ∞ (Lemma 2.3), as well as boundary measurements. On the second stage we refine the tail.

Figure 1 shows the boundary data collection scheme in our experiment. We have used six locations of the light sources (Figure 1a). Sources number 1, 2 and 3 are the ones which model the source x₀ running along the straight line L (4). The distance between these sources was six (6) millimeters. Each light source means drilling a small hole in the phantom to fix the source position. It was impossible to place more sources in a phantom of this size that mimics actual mouse head. So, since the data for the functions q_n are obtained via the differentiation with respect to the source position, we have used only two functions q₁, q₂. Fortunately, the latter was sufficient for our goal of imaging of abnormalities. Sources 1, 4, 5, 6 (Figure 1a) were used to approximate the tail function as described below. We construct that approximation in the domain Ω₁ depicted on Figure 1(b). This domain is (units in millimeters)(45)

Figure 1. The schematic diagrams of inverse problem domain and light source locations. (a) Illustrates a layout of the inverse problem setting in 2-D. The circular disk Ω corresponds to a horizontal cross-section of the hemisphere part of phantom (the supposed ‘mouse head’ in animal experiments), The computational domain Ω₁ is a rectangle containing Ω inside, six light sources are located outside of the computational domain. Because of limitations of our device, we use only three locations of the light source along one line (number 1, 2, 3) to model the source x₀ running along the straight line L. Light sources numbered 1, 4, 5, and 6 are used to construct an approximation for the tail function. (b) Depicts the computational domain Ω₁ = {(x, z), |x| < 5.83, |z| < 5.83} (unit:mm) and its rectangular meshes for tail functions and inverse calculations for the numerical method of this article. Actual mesh is much dense than these displayed. The diagram is not scaled to actual sizes.

5.2.1. The first stage of approximating the tail

Since this stage essentially relies on positions of sources 1, 4, 5, 6, it makes sense to describe this stage in detail here, although it was also described in 1. The light sources are placed as far as possible from the domain Ω₁, so that the asymptotic approximations (47), (48) of solutions are applicable. Let s^(j) be the 2-D vector characterizing the position of the light source number j. First, consider the light source number 1. Denote(46) Using (5), Theorem 2.1 and Lemma 2.3, we obtain for the function (47) Since by (45) we use (47) to approximate the unknown function p_∞(x, z₁) as(48)

Formula (48) gives the value of p_∞(x, z₁) only at z = z₁. Since Ω₁ is a square, we set the first guess for the tail as the one which is obtained from (48) by simply extending the values at z = z₁ to the entire domain of Ω₁,where the function S(x, z, s⁽¹⁾) is given in (46). Next, we compute the function and get a⁽¹⁾(x, z), (x, z) ∈ Ω via (43), (44).

For light sources 4–6, we repeat the above procedure to get a⁽⁴⁾(x, z), a⁽⁵⁾ (x, z) and a⁽⁶⁾ (x, z) respectively. Then we consider the average coefficient and set (see (44))(49) Next, we solve the forward problem (1), (2) for the light s⁽³⁾ with this coefficient a(x, z) again to get The final approximate tail function obtained on the first stage is(50)

5.2.2. The second stage for the tail

The second stage involves an iterative process that enhances the first approximation for the tail (50). We describe this stage only briefly here referring to 1 for details. In this case, we use only one source number 3, s⁽³⁾. Recall that . Let the function a(x, z) be the one calculated in (49). Denote a₁(x) ≔ a(x, z). Next, we solve the following boundary value problemThen, we iteratively solve the following boundary value problemsWe set u_m ≔ u_m−1 + w_m. Next, using (43) and (44) with u ≔ u_m, we find the function a_m+1 (x). We have computationally observed that this iterative process provides a convergent sequence {a_m (x)} in L₂(Ω₁). We stop iterations at m ≔ m₁, where m₁ is defined via(51) where ϵ > 0 is a small number of our choice, and the norm in L₂ (Ω₁) is understood in the discrete sense. Next, assuming that , we set for the tail(52) Then, using (52), we proceed with calculating of functions q_n as described above.

6. Real data

We now describe our experimental setup for collecting the optical tomography data from an optical phantom. This phantom is a man-made subject that has the same optical property as small animals. Such a phantom is a well-accepted standard to test reconstruction methods for real applications before animal experiments.

Figure 2 is a photograph of our measurement setup. The center of the picture is the phantom (rectangular box with a hemisphere on top surface) which, in particular, contains a hidden inclusion inside (not visible in this photo). The hemisphere mimics the mouse head in animal experiments of stroke studies, and hidden inclusions of ink-mix mimic blood clots. The four needles are laser fibers that provide light sources in our experiments. The fiber on the right-hand side of Figure 2 can be moved to three other positions which are source positions 1, 2, 3 on Figure 1(a). Three other fibers represent positions 4, 5, 6 of the source on Figure 1(a). A CCD camera is mounted directly above the setup (not shown in photo), and the camera focus is on the top surface of the phantom. CCD stands for Charge Coupled Device. CCD camera with its sensitivity and response range is commonly used for near infrared laser imaging of animals.

Figure 2. A photograph of the experimental setup. An optical phantom is connected with four laser fibers, the one on the right-hand side is movable to locations 1, 2 and 3 (Figure 1a). The gelatin made phantom has the shape of the rectangular box with a hemisphere on its top surface. The hemisphere mimics the head of a mouse with a mask exposing the crown part of the mouse. A CCD camera mounted above the phantom (not shown) provides light intensity measurements of the top surface of the phantom.

The purpose of this experimental setup is to study the feasibility of using our numerical method in an animal stroke model. In the case of a real animal, the top hemisphere (meshed shape in Figure 6.2a) is to be replaced by a mouse head and the rectangular block of phantom to be replaced by an optical mask filled with a matching’ fluid (a tissue-like solution with the optical properties similar to the animal skin/skull). The image reconstruction should provide the spatial distribution of the optical coefficient a(x) (directly related to the blood content) in a 2-D cross-section of the animal brain.

The goal to is to accomplish a noninvasive imaging means to monitor and investigate hemodynamic dysfunction during ischemic stroke in animal models. The current experimental setup works only for small animals, such as rats, but not for human brains due to the limited penetration depth of NIR light through a larger size and thicker bone structure of the skull. The methodology developed here, however, may be applicable to a neonate head, while further studies are needed to confirm such an expectation.

Our inverse reconstruction is performed in the 2-D plane depicted in Figure 3(b), with the optical parameter distribution of the medium inside the circle (at the center of Figure 3b) as an unknown coefficient in the photon diffusion model. Should we need a diagnosis of a different cross-section of mouse brain in animal experiments, we need a different optical mask filled with matching fluid, and repeat both the experiment and the reconstruction.

Figure 3. The schematic diagrams of our data acquisition process. (a) schematically depicts a 3-D phantom and its hidden inclusion (shown in meshed surface) 5 mm below the phantom. The light source is placed at various locations on the top surface of the rectangular block (Figure 1a) and light intensity measurements are taken on the top surface of phantom. (b) The measurement surface by a CCD camera. The data are collected from the surface of the hemisphere as well as from the un-shaded area of top rectangle, both area are lifted from the phantom in this drawing for a better illustrative purpose. The rectangular figure in Figure 3(b) also illustrates the middle a 2-D cross-section (meshed circle) of the presumed ‘animal head’, at the boundary of which light intensity data are collected for the reconstruction. The light intensity at the 2-D cross-section (except for the boundary) is obstructed by the top surface of the hemisphere. Light sources are also located in the same plane as this cross-section area. The 2-D inverse problem is solved in this cross-section by ignoring the dependence on the orthogonal coordinate.

The geometry of the phantom (shown in a vertical central cross-section) is depicted in Figure 4 with dimensions specified. The phantom is shaped by a hemisphere (diameter 13 mm) on the top of a cube of 30 mm × 30 mm × 30 mm. A spherical hollow of 5 mm diameter is located inside the phantom, with its top 5mm below the top surface (Figure 4). The location of this hollow is symmetric here but can be in any place when needed. One and two spherical hollows of 3 mm diameter inside the phantom are also used in experiments (see Section 8 for detail). We fill the hollow with liquids made of different kinds of ink/intralipid mix to model strokes by blood clots. Intralipid is a product for fat emulsion, which mimics the response of human or animal tissue to light at wavelengths in the red and infrared ranges. The phantom is made of gelatin mixed with the intralipid. The percentage of intralipid content is adjustable. Thus the phantom has the same optical parameters as the background medium of the target animal model. At the location of the inclusion (the hollow), we inject ink/intralipid mixed fluids whose optical absorption coefficients are 2 times, 3 times and 4 times higher than in the background. Also, we use the pure black ink to test our reconstruction method for the case of the infinite absorption. Different levels of the inclusion/background absorption ratio are used to validate our method for its ability for different blood clots. Light intensity measurements are taken directly above the semicircle in Figure 4. We note that Figure 1(a) shows a top view of the layout but Figure 4 shows a cross-section from side view. The disk region Ω in Figure 1(a) corresponds to the meshed area in Figure 3(b) as well as a horizontal cross-section of the 13 mm hemisphere. The reason we use such a geometry is that we need to test our capability to reconstruct inclusions at several different depth (moving up and down vertically as in Figure 4) as well as at different locations (moving horizontally in Figure 3b).

Figure 4. Dimensions of the phantom, shown in a vertical cross-section at the center of the phantom from side view. The circle in Figure 1(a) corresponds to the boundary of the meshed circle in Figure 3(b). This also corresponds to the boundary of the 13 mm diameter circle at a top view above the rectangular box (but not related to circle or semicircle in this graph). The 5 mm diameter hollow sphere is for the placement of inclusions, filled with liquids of different compositions to emulate strokes. The hidden inclusion is 5 mm below the top surface.

7. Processing real data

As it is typical when working with imaging from real data, we need to make several steps of data pre-processing before applying our inverse algorithm.

7.1. Computational domains

We need to use several computational domains. The computational domains and meshes used in our numerical calculation involve four domains. Figures 5(a)–(c) and 1(b) are the finite element meshes on each of these domains Ω, Ω₀, Ω₀∖Ω, Ω₁, respectively. In actual computations we use more refined meshes for each domain than those illustrated. Therefore, figures are not scaled to actual sizes.

Figure 5. (a) The domain (units: mm) with a triangular mesh. The domain represents a cross-section of our phantom, and the inclusion inside Ω is to be reconstructed. (b) The domain Ω₀ = {(x, z), |x| < 23.32, |z| < 23.32} (units: mm). The domain Ω₀ contains Ω and light sources. Meshed small square is the domain Ω₁ which is the same as on Figure 1(b). Computations of the inverse problem are performed in Ω₁. (c) Shows Ω₀∖Ω the domain used for data processing. To smooth out the measurement noise, Equation (1) with the Dirichlet boundary conditions u ∣_∂Ω= ϕ(x, x₀), is solved in Ω₀∖Ω. Then the smoothed data along ∂Ω₁ is used for the inverse problem.

Here, the disk Ω is the domain of interest that contains the inclusion, corresponding to the meshed area as shown in Figure 3(b). It reflects a cross-section of the hemisphere of phantom as shown in Figure 3(b). The light intensity data used in our computations are originated by the measurements at ∂Ω by taking pictures with CCD camera on the top of phantom as discussed, shown in Figure 3(b). The CCD camera collected data is acquired from the surface of the hemisphere and the non-shaded area of the top rectangle, and we extract only the needed data for ∂ Ω, i.e. the circle in Figure 3(b) by getting the data at these locations as boundary values. Here, .

On Figure 5 Ω₀ is a large domain, which can be interpreted as a truncated plane ℝ². Indeed, one cannot practically solve the forward problem (1), (2) in the infinite plane ℝ². Light sources are located in Ω₀∖Ω. The background simulation and calibration of background parameters are performed in Ω₀. To smooth out the noise in the data, Equation (1) is solved in Ω₀∖Ω for each source position, which is similar with 2. Because of (3), we use a(x) ≔ k² for x ∈ Ω₀∖Ω. In this procedure the Dirichlet boundary condition at ∂Ω is taken from the real data, and we use the zero Dirichlet boundary condition at ∂Ω₀. As a result, we obtain the smoothed Dirichlet boundary condition at ∂Ω₁ for each source location. We solve the inverse problem in Ω₁, and Ω ⊂ Ω₁.

In solving the above boundary value problems for elliptic equations, we have used the finite element mesh for domain Ω₀ with 30,404 elements and 85,449 nodes. The exterior forward problem has used 26,782 elements and 78,316 nodes for domain Ω₀∖Ω. Quadratic triangular elements and eight-node rectangular serendipity elements are used in the finite element code. A reconstruction typically takes about 2 min to complete in a Pentium 4 duo processor.

7.2. Data pre-processing and approximations of boundary conditions

As it was stated in Subsection 7.1, to smooth out our measurement data, we solve the Dirichlet boundary problem for Equation (1) with a(x) ≡ k² for x ∈ Ω₀∖Ω (Figure 5c) for each source location. Namely,(53) Here, is the background value and ϕ(x, s⁽ⁱ⁾) is the experimentally measured data for the light source s⁽ⁱ⁾. Let be the trace of the solution of the Dirichlet boundary value problem (53) at x ∈ ∂Ω₁. We solve the inverse problem in the square Ω₁ with the smoothed data We have shown in 3,21 that solving inverse problems in the original domain Ω and extended domain Ω₁ are equivalent mathematically. But numerically, the noisy component in is much smaller than in ϕ(x, s⁽ⁱ⁾). This smoothing effect takes place because the inverse problem is solved in Ω₁.

7.3. The forward problem and calibration

We now address the question on which value of k² one should use when working with the real data. The optical properties of the background of the phantom (without the hidden inclusion) are known theoretically from the concentration of intralipid in the mix. However, there is a discrepancy between the theoretical value and actual measurements. Before we solve the inverse problem to image hidden inclusions, we calibrate our model by adjusting the background value of k² to the real data measured for the reference medium, which is the phantom without inclusion, i.e. the hollow is filled with the same intralipid solution as that in the phantom itself.

First, we numerically solve the forward problem with the source position s⁽¹⁾ in the domain Ω₀ without any inclusion,(54) where is the background value and δ(x − s⁽¹⁾) is approximated by the Gaussian functionwhere ϵ = 0.01. Then we calibrate the parameter μ_a (fixing ) as well as the amplitude A > 0 of light source in our model (54) to match the measured light intensity for s⁽¹⁾ for the uniform background. We do not know the number A. Thus, we choose A in such a way that u_comp (x_max, s⁽¹⁾) ≈ u_meas (x_max, s⁽¹⁾). Here x_max is the brightest point, i.e. the far right point on ∂Ω, being closest to the light source. Also, u_comp (x_max, s⁽¹⁾) and u_meas (x_max, s⁽¹⁾) are computed and measured light intensities respectively. Next, we should approximate the constant k². To do this, we take another sampling point x_min with the minimal light intensity, which is the farthest left point on ∂Ω. Next we consider ratios R_comp (k²), R_meas, whereThese ratios are independent on the number A in (54). We choose k² such that R_comp (k²) ≈ R_meas. As a result, the calibrated value of k² was k² = 2.403. This computed value matches quite well the theoretical value of 2.4 of the intralipid solution we have used.

8. Results of the reconstruction

Let a_incl = a(x) be the value of a(x) inside the inclusion and a_b = k² = 2.403 be the value of the coefficient a(x) in the background, which was computed in subsection 7.3. Our ratios a_incl/a_b where(55) The value a_incl/a_b = ∞ means that the inclusion was filled with a black absorber, i.e. black ink. table 1 gives a detailed information of our real data with different positions, numbers, diameters and contrast of inclusions. ‘Y’ means that the experiment was done for this setting, and ‘–’ means that the experiment was not done. Note that in Group 3, we have imaged two different inclusions simultaneously.

Table 1. Real data.

Inclusion groups	Positions	Diameters	Ratio2	Ratio3	Ratio4	Ratio∞
Group 1		5 mm	Y	Y	Y	Y
Group 2		3 mm	Y	Y	Y	–
Group 3		3 mm	Y	Y	Y	–

As described in Sub-subsection 5.2.1, we construct the ‘asymptotic tail’ on the first stage using light sources 1, 4, 5 and 6 as an initial approximation, to be followed by other subroutines for further refinements. The image of the function a(x, z) in (49) is depicted on Figure 6, for a phantom where the theoretical value of the inclusion/background contrast (55) was a_incl/a_b = 3. Images from the asymptotic tails for other phantoms listed in Table 8.1 were similar.

Figure 6. The function a(x, z) in (49) is depicted. This function was obtained on the first stage procedure for the tail (Subsubsection 5.2.1) as an initial approximation. The initial reconstruction is obtained from a phantom where the theoretical value of the inclusion/background contrast (55) was a_incl/a_b = 3.

Figures 7–9 depict the reconstructed images from our real data. Figure 7 depicts the 3-D plots of reconstructed images of group 1 for contrast values of: 2:1, 3:1, 4:1 and ∞ : 1, from the left to the right respectively. Figures 8 and 9 show 2:1, 3:1, 4:1, from left to right respectively for groups 2,3. In all our examples ϵ = 10⁻⁵ in (51). The finally reconstructed results for contrasts are listed in Table 2. Note that our above reconstruction algorithm does not use any knowledge of neither the location of the inclusion, nor the contrast value a_incl/a_b.

Figure 7. Reconstructed inclusion contrast of data group 1, actual contrasts are 2:1, 3:1, 4:1 and ∞ : 1, from left to right, respectively. The transparent frames show the theoretical values of inclusion/background contrasts of actual inclusions, which are made of different ink-intralipid mix.

Figure 8. Reconstructed inclusion contrast of data group 2, actual contrasts are 2:1, 3:1, 4:1, from left to right, respectively. The transparent frames show the theoretical values of inclusion/background contrasts of actual inclusions, which are made of different ink-intralipid mix.

Figure 9. Reconstructed inclusion contrast of data group 3, actual contrasts are 2:1, 3:1, 4:1, from left to right, respectively. The transparent frames show the theoretical values of inclusion/background contrasts of actual inclusions, which are made of different ink-intralipid mix.

Table 2. Reconstructed values of the contrast within imaged inclusions and they relative errors, a_b = 2.403, compare with (8.1).

	The true contrast aincl/ab		Relative error
Group 1	2	2.11	0.056
	3	2.9	0.032
	4	4.22	0.057
	∞	6.69	Unknown
Group 2	2	2.33	0.167
	3	3.29	0.0986
	4	4.57	0.143
Group 3	2	2.29	0.148
	3	3.49	0.164
	4	4.58	0.147

9. Summary

We have worked with a set of real data, using the approximately globally convergent numerical method of 1 for a CIP for an elliptic equation. These data mimic imaging of clots in the heads of a mouse. We have introduced a new concept of the approximate global convergence property and have established this property for the discrete case of a finite number of finite elements, unlike the continuous case of 1. Figures 7–9 as well as Table 2 demonstrate that our are quite accurate, including both inclusion/background contrasts and locations of inclusions. Note that accurate values of contrasts are usually hard to reconstruct via locally convergent algorithms. On the other hand, these values are especially important for our target application, since they might be used for monitoring stroke treatments. Therefore, following the concept of Steps 1–6 of Subsection 4.1, we conclude that our approximate mathematical model of Subsection 4.3 is a valid one.

Acknowledgements

The work of all authors was supported by the National Institutes of Health grant number 1R21NS052850-01A1. In addition, the work of MVK was supported by the US Army Research Laboratory and US Army Research Office under the grant number W911NF-11-1-0399.