Reconstructions for Continuous-Wave Diffuse Optical Tomography by a Globally Convergent Method

In this paper, a novel reconstruction method is presented for Near Infrared (NIR) 2-D imaging to recover optical absorption coefficients from laboratory phantom data. The main body of this work validates a new generation of highly efficient reconstruction algorithms called “Globally Convergent Method” (GCM) based upon actual measurements taken from brain-shape phantoms. It has been demonstrated in earlier studies using computer-simulated data that this type of reconstructions is stable for imaging complex distributions of optical absorption. The results in this paper demonstrate the excellent capability of GCM in working with experimental data measured from optical phantoms mimicking a rat brain with stroke.


Introduction
The studies using Near Infrared light (NIR) for biomedical imaging have become quite extensive in the past 15 -20 years.Earlier applications of NIR imaging started with breast imaging [1][2][3], since NIR has better penetration depths in soft tissues than visible light.Applications of NIR technologies have made tremendous progress, including the detection of brain injury/trauma [4], the determination of cerebrovascular hemodynamics and oxygenation [5][6][7][8], and functional brain imaging in response to a variety of neurological activations [9][10][11].Recently, the technique was also successfully used for measuring the efficacy of photodynamic therapy for prostate cancer [12].Among several common NIR imaging mechanisms Frequency Domain (FD) imagers were developed and used in patients [13], and Time Domain (TD) methods were used for brain studies in [14,15].As a low-cost alternative to the TD and FD imaging systems, Continuous-Wave (CW) NIR imaging systems have become more commonly used in recent years.
To spatially quantify brain hemodynamic activities resulting from functional neuronal signals, it is desirable to extract distributions of light absorption from light intensity measurements through mathematical models.Since these optical properties are described by coefficients in the light diffusion model [16,17], one needs to solve an inverse problem of the corresponding partial differential equation, the diffusion equation, to obtain diffuse optical tomography (DOT).
The inverse reconstruction for DOT has been mathematically challenging.The problem is nonlinear as well as ill-posed.As reviewed by Hielscher in [18], the majority of the inverse reconstruction algorithms used for DOT has utilized a perturbation approach involving the inversion of large Jacobian matrices, for examples [19][20][21].The successful reconstruction of the unknown coefficients depends highly on the initial guess being close to the true solution since the least-square residues might have multiple local minima [22].
Our current focus of theoretical research is on numerical methods that have no restrictions on the initial guess, namely, the development of a Globally Convergent Method (GCM) as initiated by Klibanov and Timonov [22], [23].Our earlier generation of a GCM is called a "convexification method," as introduced in [23], based on modified residual error estimates.When a proper Carleman weight function is added to the error terms, those residues depend on the unknown variable in a convex manner.Our current approach (2 nd generation of a GCM) is called a continuation method or homotopy method, as defined in [24].The homotopy connects the sought system with a similar system that is easier to solve.In this approach, our inverse reconstruction is a continuation of reconstructions from a DOT problem where light sources are far away yielding a "tail function" [25,26].This new generation method has achieved satisfactory results in reconstructions of simulated CW data [27].A similar GCM was developed for the case of the reconstruction of a coefficient in a wave-like time dependent partial differential equation by Beilina and Klibanov [28,29].This work proceeds to show the validity of this method based on physical measurements using laboratory optical phantoms with the shape of a rat head.In [30], we presented the preliminary results where we proved the convergence of the method and showed a few examples using phantom data.The focus there was to show the convergence of the method is mathematically guaranteed by the "Approximate Globally Convergent Theorem" and it is feasible to reconstruct in phantom experiments.The current paper is a more in-depth and comprehensive study of the reconstruction method based on phantom data.It contains the capability and limitations of the proposed method, error analysis in using the GCM in practical experiments.The study shows that there are two major areas that GCM clearly is superior to the conventional DOT methods.1. GCM can reconstruct the absolute value of absorption coefficient vs. conventional method can only reconstruct the relative change of absorption coefficient.2. GCM can tolerate discrepancy and fluctuation of reduced scattering coefficient, and still give quality reconstructions of absorption distribution.In practice, it gives the justification of use DOT in stroke and other applications where the scattering does fluctuate.This is the main contribution of the method to the DOT research community.
We present the model in Section 2. Section 3 presents the experimental setup for the optical phantom study and its data acquisition scheme in 3-D geometry to be processed for 2-D tomography.In Section 4 we present results of reconstructed absorption coefficients.We summarize by conclusions and discussions in Section 5.

The Optical Diffusion Model
An optical light source can be used in three forms in DOT: 1) the light is amplitude-modulated at a radio frequency (RF) in Frequency-Domain NIR imaging, 2) the light has a short pulse with a pulse width of a few picoseconds used in Time-Domain NIR imaging, and 3) the light has a Continuous-Wave (CW) form and the amplitude is independent of time.In this paper we only discuss the CW-based problem.The time-independent diffusion equation and boundary conditions of the CW case can be found in [16]: x s can be chosen at several points by varying s .For 2-dimenional inverse problem of Equation (1), the source is modeled by a Dirac function in the domain 0 Ω in the same 2-demensional plane.The forward problem domain 0 Ω should in theory be the entire 2-dimensional plane, but in calculation it is a truncated rec- tangle [29].We call 0 Ω the background domain throughout this paper.Equation (1b) is commonly known as Robin boundary condition [16].
In real applications, the domain of interest is a bounded and irregularly shaped domain contained in the background domain 0 Ω .We call it the physical domain and write it as Ω , e.g. the cross section of a rat brain.The distribution of the absorption parameter ( , ) a x y µ is unknown in the physical domain Ω .It is closely related to oxygenated (HbO) and deoxygenated (HbR) hemoglobin concentrations [15], which is an important index for imaging in functional brain research.The locations of the light sources are outside of the physical domain Ω .When the forward problem ( 1) is restricted into the physical domain Ω we have where ( ) ( ) , , , , x y s x y ϕ ∈ ∂Ω is the light intensities which can be measured on the boundary of domain Ω .The right hand side of (2a) is equal to 0 since the light sources are located outside of the physical domain on the boundary of Ω and ( , ) a x y is a new unknown coefficient, defined as: where ∆ represents the Laplacian operator.For the rest of this paper the focus is on the reconstruction of the parameter distribution ( , ) a x y .This parameter ( , ) a x y reflects both the light absorption and the reduced scat- tering coefficients.When the scattering coefficient ' ' ( , ) with a unit of cm −2 .The Influence of inhomogeneity and discrepancy of scattering coefficient will be discussed in the last section.

If ( , )
x, y s Φ  is known in the whole physical domain Ω we can compute the coefficient ( , ) a x y in Equation (3) in a straightforward way by computing the second order derivative of ( , ) x, y s Φ  . A better and more stable solution [25], which involves only the first order derivative of ( , ) x, y s a x, y x, y s dxdy where η is taken as any quadratic finite element function.
However, we only know ( , ) x, y s Φ  on the boundary of Ω in real applications.So the idea of the reconstruction for the coefficient ( , ) a x y is to construct an approximation of ( , ) x, y s Φ  in the whole physical domain Ω from the boundary measurement data.In order to achieve this goal, a sequence of transformations of Equation ( 3) is performed.The transformed equation can be solved in a regular domain, e.g. a rectangle domain which contains the physical domain Ω [31].Once we have an approximation of ( , ) x, y s Φ  , we use Equation (5) to compute the coefficient ( , ) a x y .

Transformations
The first transformation is to let ln ( , ) u x, y s = Φ  .Then Equation (3a) becomes the nonlinear elliptic equation ( , , ) ( , , ) ( , , ) ( , ) 0. u x y s u x y s u x y s a x y The second transformation is to let 2 / v u s = and ( , , ) ( , , ) , where s is the second coordi-nate of the light source ( , ) m x s .We obtain that ( ) s s s s q q x y d s q x y d s q x y d T s T q T s s The term ( , ) T x y is the so-called "tail function" and is defined to be ( , ) ( , , ) T x y v x y s = .This is an integral-differential equation without the unknown coefficient ( , ) a x y .If we are able to approximate well both functions T and q , and their derivatives, then we can approximate well the target coefficient ( , ) a x y via Equation (5).We use the asymptotic relation to approximate ( , ) T x y and then solve a discrete version of Equa- tion (7) to solve the problem.As discussed in the introduction, the mathematical underpinning of transforming the inverse problem into Equation ( 7) is the mathematical method of Homotopy, i.e., using a free parameter (in our case, the light source location s), we create a continuous passage from one difficult problem at one end of the parameter to another "easier" problem at the other end.Once we have solved the relatively easier problem, we change the free parameter continuously and continue to solve the perturbed problem until we reach the other end.The relatively easier mathematical problem here is that the light source moves to infinity, where the light intensity (photon flux) at the unknown medium has an asymptotical formula to be discussed below.This is the additional information required to solve Equation (7).For more details about the algorithm, we refer to [32].

The Phantom Experiments
Our DOT imaging system is set for surface data collection and 3-D reconstruction through 2-D tomography of parallel cross-sections.The envisioned scenario is that of a rat brain which is covered with an optical mask that is filled with a "matching material", a gelatin tissue phantom that has similar optical properties to those of the rat's skull/skin.The purpose of the experiment is to determine the unknown location and severity of a stroke in the rat brain phantom.Two optical phantoms are specifically designed and made for testing GCM-based imaging reconstructions, as a feasibility study for GCM applications in rat stroke models.The 3-D geometry of the phantoms is depicted in Figure 2.Both phantoms are shaped as a hemisphere with a diameter of 13 mm on top of a cube with dimensions 30 mm × 30 mm × 30 mm.The top hemisphere, the meshed shape in Figure 2(a), mimics a rat head, and the cube of the phantom emulates an optical mask filled with the "matching material".Spherical hollows of 2.5 -3.0 mm diameter are within the phantoms.Because the said cross section is obstructed by the top surface of the hemisphere in the CCD's field of view, we can only collect the light intensity data at the boundary of a 2-D cross section of the presumed "animal head" as depicted in Figure 2(b).
In our experiments the optical coefficients of the phantoms are 1/ (3 ) 0.037  for the inclusions by adjusting the ink-Intralipid mix.An inclusion with pure black ink has a theoretical value of a µ that tends towards infinity.These optical parameters are standard values used for our rat brain phantom ex- periments.

Reconstruction Results
The reconstruction is done without assuming any prior knowledge of the total number, location or optical parameter values of the inclusions.We collected six groups of "phantom data" with inclusions and an additional group of measurement data for the phantom without any inclusion as the background reference.The data are the light intensity at the boundary of the physical domain Ω collected by a CCD camera.They served as Dirichlet boundary condition while solving the inverse models.Within the six groups, the first three have one inclusion each of which has different absorption contrasts; see the upper panel of incl bk a a = and 4, respectively.Figure 3 shows the reconstructed images of all six groups of data.The recovered positions in the first three groups are within the circle of the actual position of the single inclusion, the dashed circle.The recovered positions in the last three groups are also approximate to the actual locations of the inclusions; the two dashed circles, except that the distance between the two reconstructed inclusions is closer than the actual distance.However, we observe from the lower panel of Figure 3 that the distance between two reconstructed inclusions becomes larger when the contrast of the inclusions increases (from left to right in the figure).Especially, in the figure of two inclusions with contrast 4:1, the lower right panel, the two inclusions can be better resolved from each other.
Listed in Table 1 are the maximal contrast values and relative errors for all six cases.It shows that we have imaged the contrast values of the inclusions in the same scale of magnitude.These contrast values are important as they can be further used to calculate HbO and HbR levels.In Table 2, we show the localization errors of our reconstructions.The absolute localization error is defined as the distance between the actual center actual C and the center location recon C of a reconstructed inclusion where ( , ) a x y attains a maximal contrast value.The relative localization error is defined as where R is radius of Ω .We have observed that the localization errors tend to increase with the number of inclusions.

Conclusions and Discussions
We have applied the GCM algorithm of [32] for a Coefficient Inverse Problem to DOT phantom measurements.These data mimic imaging of a blood clot in a rat brain.Reconstructions are performed without any prior knowledge of the inclusions.In our 6 cases, we found that reconstructed images of inclusions from NIR experimental data are quite accurate, including both locations of inclusions and inclusions-to-background contrasts.The main contribution of this paper is the validation of GCM by acquiring the phantom data from phantom boundary measurements, and reconstructing without assuming any prior knowledge of the total numbers, locations and the optical parameter values of inclusions.
Our GCM compares favorably with other conventional locally convergent methods.We have tested the same data sets with a conventional DOT reconstruction method which is a Tikhnov regularization based method with   depth compensation [33,34] and the regularization parameter is 0.001.The conventional DOT method was performed using a publicly available software HomER (http://www.nmr.mgh.harvard.edu/PMI/resources/homer/home.htm)[35].The computational times for both methods are similar, about 1 -2 minutes [29].
In Figure 4, we show the reconstructed images by the conventional DOT reconstruction method for the oneinclusion cases.The locally convergent method can identify the changes in the absorption coefficient well; however, note that the conventional DOT algorithm cannot reconstruct the absolute value of ( , ) a x y due to the limitations of the method.The images in Figure 4  in the range [1,10] which is in the same range as actual change.For the two-inclusion cases, the conventional DOT did not differentiate the two inclusions.
As stated in Section 2, we assume the scattering coefficient ' ( , ) s μ x y is known and uniform in our experiments.We also assume that ' s μ keeps unchanged during the background measurements and measurements with inclusions.The reconstructed contrast / incl bk a a reflects the change in absorption coefficient ( , ) a μ x y .How- ever, in application for small animal measurements the scattering coefficient is neither homogeneous nor unchanged with time.In order to quantify the errors in recovered contrasts in ( , ) a μ x y with respect to discrepan- cies in the assumed background scattering coefficient, numerical simulations have been done to provide the comparison.Commonly, 5% additive noise is included in the simulated light intensities of the forward problem.Light sources are located at (x,y) = (−2,0), (2,0),(0,−2), (0,2), (0.2,2), (0.4,2) in the background domain.The first 4 sources are used for obtaining the initial tail function.The last 3 sources are used for solving the integral-differential equation in Equation (7).The detector is supposed to be a CCD camera taking in the Dirichlet boundary value, the simulated light intensity data with noise on the boundary of the physical domain are the measurement data ( , , ) x y s ϕ in Equation (2b).The background domain for forward problem is [ 3,3]    a μ x y are pretty accurate (less than 3.0%) for all the test cases.
We then studied the case where the scattering coefficient changes after the background measurements were done.Its purpose is to test the capability of the method in dealing with animal data where scattering coefficient fluctuate with time during experiments.The changes of scattering coefficients occur in locations where inclusions are at.The error pattern is the same as the pattern shown in the middle of the upper panel in Figure 5.The reconstructions were also successful: the locations of the inclusions are precise.However the relative peak errors change with the magnitude of fluctuations in scattering coefficients.The recovered contrasts in ( , ) a μ x y are shown in Figure 6.It shows that in this case, the errors in the recovered contrasts depend on the errors in the scattering coefficient linearly.The increase of reconstruction errors is due to the sensitivity of the solution to the scattering coefficient variance.
In order to investigate the effect of cross talk for the locations of the reconstructed inclusions under the influence of random errors in scattering coefficient, a comparative study of GCM method and conventional method has been performed.Two circular inclusions centered at (−0.25,0) and (0.5,0) with radius 0.3 are considered.The maximum value of a µ in the inclusions is 0.3.The upper panel of Figure 7 shows the reconstructed loca- tions by using the conventional method and the lower panel shows the results by our GCM method.The accurate locations are marked out by dotted circles.It should be noted that in order to give a better understanding for the results of the conventional method, the results are not cut by a threshold.For GCM method, it requires a cut at the stage of computing the tail function.We can see that at the 15% error level in scattering coefficient, the conventional method failed to reconstruct the two inclusions, while our GCM is able to give the locations of the two inclusions despite of some artifacts in between the two inclusions.
In summary, we have validated our GCM reconstruction algorithm based on real data acquired by a NIR CCD Camera using optical phantoms, after testing successfully the GCM method on computer simulated data in previous papers [27,31,32].This study has confirmed the ability of GCM to reconstruct improved DOT images using laboratory rat-head phantoms.

Figure 1 (
a) is a photo of the measurement setup.The center of the picture displays the optical phantom which contains the hidden inclusions that are not visible in the photo.The four thin tube-like probes are laser fibers that provide the light sources.A diode laser (Coherent Inc. wavelength at 808 nm) is multiplexed to serve as the source by a multiplexer (Avantes Inc.Multiplex Channels 1 × 16), which delivers and controls light at four locations through the four laser fibers.The fiber on the right-hand side of Figure 1(a) can be moved to other positions by adjusting the mechanical arm that controls the lateral distance.A Charge-Coupled Device (CCD) camera for light intensity measurements is mounted directly above the setup.A schematic drawing of the experimental setup is shown in Figure 2(b).

1 Figure 1 .
Figure 1.(a) A photograph of the experimental setup (b) A 3-D schematic drawing of experimental setup for phantom study.

Figure 2 .
Figure 2. (a) A schematic drawing for the 3-D geometry of an optical phantom.(b) The data acquisition process of on the boundary of the physical domain (meshed area) for Equation (2b), from light intensity measurements of 3-D phantom surface.

Figure 5 .
The last three groups have two inclusions; see the lower panel of Figure3, also with three different contrasts.Let of the inclusion, and bk a be the background value.Then the contrasts in our experiments are / 2,3

Figure 3 .
Figure 3. Reconstructed a(x,y) for the six groups of data.Upper panel: one inclusion cases.Lower panel: two inclusions cases.The contrasts are arranged in the order 2:1, 3:1 and 4:1 from left to right.The actual locations of inclusions are marked out in the dashed circles.
depict ( , ) a x y ∆ , the change of the coefficient in convention-al DOT from the background in the range [−5 × 10 −4 , 20 × 10 −4 ], while the images in Figure 3 by GCM depict ( , ) a x y ∆ [ 3,3]  − × − and the mesh is a 120 by 120 rectangular grid.The interested physical domain for inverse problem is[ 1,1] [ 1,1] − × − and the mesh for inverse problem is a 40 by 40 rectangular grid.Three different error patterns in scattering coefficient have been studied, as shown in the upper panel of Figure 5.The amount of errors

Figure 4 .
Figure 4. Reconstructed results (changes of a(x,y) from background) for the first three groups of data by the conventional DOT algorithm.The contrasts are arranged in the order 2:1, 3:1 and 4:1.Actual positions of inclusions are marked out by dashed circles.

Figure 5 .
Figure 5. Upper panel: three typical error patterns in the background scattering coefficient are shown.Lower panel: the recovered contrasts (lines with square symbols) and actual contrasts in ( , ) a μ x y are shown in each figure for tests with contrasts 2, 3 and 4. In each case, 31 simulations are done with respect to the errors in the scattering coefficient that range from −15.0% to 15.0%.

Figure 6 .Figure 7 .
Figure 6.The recovered contrasts (lines with square symbols) and actual contrasts in ( , ) a μ x y in the cases of that the scattering coefficient changes after the background measurements