1. Introduction

In recent years there has been an increased interest in using diffuse optical tomography (OT) for imaging the human breast and brain [1–5]. This novel medical imaging technique attempts to reconstruct the spatial distribution of optical properties (absorption and transport scattering coefficients, µ_a and µ’_s) within the body from the measurements of transmitted and reflected near-infrared light. Three kinds of techniques are employed to obtain the information used to reconstruct the image: time resolved method [1,3,5], frequency domain method [2,4], and continuous intensity method (DC case) [6,7]. In time resolved case, an ultra-short pulsed laser is used as the light source and the temporal distribution of the light emerging from the tissue surface is detected using either a synchronous scan streak camera with picosecond time resolution or a time-correlated single photon counting (TCSPC) system. In frequency domain method, a radio frequency modulated light source is used to illuminate an object and the phase and modulation depth of the transmitted light are measured, usually using a heterodyne detection method. Both the time resolved and frequency domain schemes exhibit high temporal resolutions (in the picosecond range for the time resolved and nanosecond for the frequency domain) [3]. However, tomographic measurements made with the time resolved technology require data-acquisition time on a time frame of minutes to derive data having acceptable low noise levels [8]. In addition, they have the disadvantages of high cost and complexity in the equipment. In DC case, a point laser source with constant intensity is used to continuously illuminate the object and the amplitude of the outgoing light on the boundary is measured. The continuous intensity system has the advantage of low cost and high dynamic range, as well as a relative high signal to noise ratio (SNR) [8]. The drawbacks include the difficulty in simultaneous unique recovery of diffusion and absorption coefficients [9] (though a normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography has been developed by Y. Pei. [10]) and inability to probe objects embedded deep inside the medium [3,11].

Various reconstruction methods for OT have been proposed, which include the perturbation based [7,12,13] and the gradient based methods [1,2,4–6,14]. The advantage of the gradient-based iterative image reconstruction scheme as compared with the perturbation methods has been discussed in detail by A. H. Hielscher [1]. The gradient based method contains two important components: the efficient calculation of the gradient of the objective function and the optimization algorithm used for minimizing the objective function. The adjoint source [5] and adjoint differentiation [14] methods, both called adjoint models, are widely employed in the gradient calculation. The efficiency of different adjoint models have been analyzed by A. D. Klose [15] extensively. The Conjugate Gradient (CG) method [1,5] can be adopted for the optimization of the objective function. Arridge [5] have made a comparison between the CG method, Steepest Descent (SD) method, and the Algebraic Reconstruction Technique (ART) for a time resolved case, which shows that the CG method is preferable. The gradient-based method has been successfully applied to both the time resolved scheme and frequency domain scheme [1,5]. For the DC case, A. J. Davies et [14] demonstrated, using simple simulated data, that the CG method was time consuming and sometimes did not converge.

In this report we describe a spatial location weighted, gradient-based optimization scheme for DC imaging, and we demonstrate the ability to efficiently recover the depth information using a one-dimensional line search with different weighting factors. To illustrate the performance of this method, we simulate a circular object with one or three absorbing objects embedded in a homogeneous background. The reconstructions results clearly demonstrate that this method can be applied to recover the absorbing objects in a highly scattering medium for the DC case. The organization of this paper is as follows. Section 2 gives a review of the theory of forward and reverse problems in optical tomography. The gradient calculation using the adjoint model is also discussed. In Section 3, our spatial location weighted optimization method is presented, followed by a demonstration of several cases of reconstructing absorption and diffusion images using weighted and unweighted optimization in Section 4. Finally, we discuss the acceleration of convergence and improvement of localization associated with our spatial location weighted optimization scheme.

2. Theory review

2.1 Forward and inverse problems in OT

a) Forward model

Light propagation in a tissue medium can be described by the well-known diffusion equation [16]. For a domain Ω having a boundary ∂Ω, this is represented by the following expression:

where Φ(r) is the isotropic photon density at position r, r_s is the position of an isotropic DC source, and D(r) and µ_a(r) are the position-dependent diffusion and absorption coefficients, respectively. The boundary measurement Γ(ξ) at ξ∈∂Ω is related to Φ(r) by:

where n^⃑ is the outer normal vector of ∂Ω at ξ. The collimated source incident at ζ∈∂Ω is commonly represented by a diffuse point source located at a transport mean free path l_tr beneath the surfaces [17].

b) The inverse formulation

Let us consider an experimental setting that include S point light sources located at ξ_j∈∂Ω (j = 1, …, S), and M_j measurement positions ς_j,i∈∂Ω (i = 1, …, M_j) for each source j, and define the following objective function:

where Γ_j,i represents the photon intensity measured at position ς _j,i with the incident light source located at ξ_j. The subscript c denotes the values calculated by the forward simulation and the subscript me denotes the experimental (or in this case synthetically generated) values. The optical inverse problem can be solved by minimizing the objective function.

2.2 The gradient calculation

To solve the optimization problem using a conjugate gradient method, the gradient of the objective function needs to be calculated. Arridge et al [5] have developed an adjoint source scheme to calculate the gradient. They represent the gradient vector z⃑ as the following form:

where the Jacobian J of the forward problem is a c×N_TOT matrix, where $M_{\mathit{TOT}}=\sum _{j=1}^s{M}_j$ is the total measurement number and N_TOT is the number of the coefficients to be reconstructed. b⃑ is the residual data error for the measurement. Since the use of matrix J explicitly to derive the gradient ∇E is prone to be too computationally demanding (requires 2M_TOT forward computations), they developed a much more efficient way to compute ∇E directly using two forward computations for source vector and weighted measurement vector (requires only 2S forward computation) [18]. However, the computation of J is still of interest for sensitivity analysis of the imaging problem. The essential calculation of the J is based on the establishment of the PMDF (Photon Measurement Density Function, as defined in [11]). In [11], the Jacobian calculation utilizing the reciprocity principle in the FEM (Finite Element Method) framework has been described in detail. It has also been shown that the PMDF has maxima close to the source and detector positions and falls off toward the interior of the object for a homogeneous background, especially for DC case. Since the residual data error b⃑ in Eq. (4) is independent of the points being reconstructed, the gradient ∇E, which is used for parameters updating, will favor the points near to the surface. When a conjugate gradient method is applied, the convergence would be very slow for objects embedded in deeper layers.

A more clear interpretation to the problem just mentioned can be found when we take the reverse differentiation scheme (also called adjoint differentiation), which is used by A. H. Hielscher [1] and A. J. Davies et al [14], to derive the gradient. In [14], the derivative of the objective function Eq. (3) with respect to the optical properties D (can either be absorption or diffusion coefficient) is given by:

where K_el and K*_el are local stiffness matrix and adjoint matrix separately [₁₄]. m = 1,…N_TOT is the node number. The first sum is carried out for all elements and the second sum for all nodes belong to each element. From the calculation procedure of K* (for detail see [14]), we can know that it has large value on nodes near to the source and to the surface for a homogeneous initial guess. The second term in the sum is only related to the area of element [13]. So the same problem as the adjoint source scheme will be encountered when a CG method is applied. In [14], A. J. Davies demonstrated that the CG method would terminate after 9000 iterations before convergence had occurred (in their second example). In addition, the CG method fails to converge for some initial guesses.

3. The spatial location weighted reconstruction scheme

The problem discussed in the previous section implies that the CG method alone does not work satisfactorily for DC image reconstruction. In order to improve the performance of image reconstruction using the CG method, we have developed a spatial location weighted gradient-based optimization scheme which employs a three-step process. First, the gradient vector of the objective function is calculated using the reverse differentiation scheme, as described in A. J. Davies’s report [14]. Second, a steepest descent method is employed for the initial iteration. Instead of using the descent direction in Euclidean-measurement, we constructed a weighting matrix A and make the descend direction in the A-measurement. This gives more weight to the interior points in minimizing the objective function. In addition, the new algorithm also serves to accelerate convergence. In the third step we return to the CG method to get the final results.

3.1 Steepest descent directions: A-norm of z⃑

We apply a steepest descent method in our new algorithm for the first iteration. Instead of adopting a descent direction d⃑ in Euclidean measurement, as indicated by the following equation:

we seek a descent direction in A-measurement. This is accomplished by constructing a diagonal positive definite matrix A of N_TOT×N_TOT whose elements are in direct proportion to the distance from the points to the surface.

Whereas a variety of weighting methods could be adopted, we have chosen one that has the form shown in Eq. (7). For a circular object, L_i can be the radial distance between the point i and the surface. β is a weighting factor.

Using this new method, the minimizing problem is defined as:

The solution can be expressed as:

And then, a one-dimensional line minimization along the current descent direction is performed (In our algorithm, the minimum along the search direction is approximated using a Golden Section search procedure and stopped after a predefined termination criterion is satisfied). We will find, from Eq. (9), that the descent direction will give more weight to the interior points when an appropriate β is selected.

3.2 The choice of the weighting factor β

In order to find the optimal weighting factor, one-dimensional line search along different directions corresponding to different values of β is performed to minimize the objective function. Various techniques can be applied to perform such minimization [19]. In this report, we use a very simple method: set initial β = 0 and increse β by 1 each time until the optimum β is found. When β = 0, it is the same as the conventional steepest descent method. The larger the β is, the more weight the interior points have. The computational burden of the objective function for a specific β value is comparable with that of one iteration in a conjugate gradient method. After the minimum of the objective function is reached with a specific β, the initial guess of the coefficients are updated. Then we return to the conjugate-gradient method to get the final results. The following algorithm briefly describes the process of our new method.

4. Simulation results

To investigate the efficiency and the accuracy of Algorithm 1, we used finite element method (FEM) to simulate a circular object having a radius of 25 mm and one or three absorbing objects embedded in a homogeneous background, each with a radius of 2.5 mm, as shown in Table 1. The finite-element grid used comprised 958 elements for one absorbing object (see Fig. 1), and 4880 for three absorbing objects. The optical properties of the background medium used for the initial guess were µ_a = 0.025 mm^-1 and µ’ _s = 2.0 mm^-1 (D = 0.16461 mm). All reconstructions started from a homogeneous background.

Table 1. Optical parameters of the circular test object

View Table | View all tables in this article

The forward data used in the reconstruction of these simulations were calculated for 36 source and 36 detector positions, equally spaced along the circumference every 10 degree, using a FEM light transport model in continuous intensity mode [20]. The FEM forward engine was based on the FEMLAB2.1 (COMSOL Inc., Burlington, MA 01803, USA, a FEM solution in MATLAB (The MathWorks, Inc.)), which provided a powerful interactive environment for modeling and solving scientific and engineering problems based on partial differential equations (PDEs). A model of Helmholtz’s equation, which has the form:

was adopted to solve Eq. (1). The coefficient c and a in Eq. (10) can be determined according to Eq. (1). ƒ is a point source with a unit amplitude. All the codes were programmed in MATLAB6.0.

The objective of these simulations is to demonstrate that the spatial location weighted optimization technique is more efficient, accurate, and computationally less expensive than the unweighted optimization method. Furthermore, this method accelerates convergence.

Figure 1 shows the reconstruction results with our new method described above. The left column of Fig. 1 shows the reconstruction results using the conventional CG method after several iterations. The middle column is the results using the spatial location weighted steepest descent method with different weighting factors β. Row one of the right column is the geometry and mesh for the simulation. Row two of the right column is the original target image for the embedded objects. The objective functions, E, of the reconstructed images with different weighting factors are shown in row three of the right column. We find that the minimum objective function occurs at β = 8, which is the optimum weighting factor. The computation time for an iteration of the spatial location weighted steepest descent method is approximately equal to that of the CG method. The one dimensional line search of the weighting factor takes no more than 9 steps (we take 11 for illustration) to find the optimum β.

 

Fig. 1. Reconstructed images after several iterations using conjugate gradient method (col. 1) and weighted steepest descent method with different weighting factors (col. 2), geometry and mesh of the simulation, target image, and the objective function E as a function of weighting factor (col. 3)

Download Full Size | PDF

In the above simulation, the embedded object is placed at 10 mm beneath the surface. Two other simulations, with the embedded object located at 5 mm and 15 mm beneath the surface separately, have also been conducted. The results of the three different cases are summarized in Table 2.

Table 2. Summaries of the results for the three simulation cases

View Table | View all tables in this article

where O: Origin of the embedded object

P: coordinate of the reconstructed peak µ_a value

Dev: Deviation of P from O

WSD: Weighted steepest descent method

M: Mean square of relative errors

N: Number of iteration

From Table 2, we will find that the WSD method exhibits an extraordinary fast convergence rate. In all three cases, the Dev (The deviation of the location of the peak µ_a value from the center of the original location of the embedded object) and M of the WSD method are smaller than that of the CG method, indicating that the WSD method recovers more precise position information of the embedded object than the CG method does. None of the methods recovers the absolute values of the embedded objects’ coefficients very well, but the new method gives better results. Though the error functions of the WSD method are larger than that of the CG method, it does not mean that their final results are inferior. Both algorithms are more sensitive to the perturbation closer to the boundary, but the new method can recover more deep heterogeneity compared to the conventional CG method.

Figure 2 shows the results from a more complex case. Three embedded objects with different parameters (see Table 1) are placed in a homogeneous background. It is only the bottom left object that has perturbations in both coefficients, while the top and bottom right objects have perturbations in only µ_a and only µ’ _s, respectively. The two images in the left column of Fig. 2 show the original target images of the embedded objects. The middle column shows the results using the conventional CG method after 100 iterations. The right column is the results using the WSD method with the optimal weighting factor (β = 7). It is worth noting that both the conventional CG method and the WSD method present cross-talk in the recover images. The increase in an object’s absorption coefficient (or scattering coefficient) is recovered as a decrease in its scattering coefficient (or absorption coefficient). So the intensity measurements alone are insufficient to distinguish effects of absorption from scattering, as Arridge [9] have claimed (Y. Pei [10] has reported a normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography. This method was not used in this article).

 

Fig. 2. Target image (col. 1), reconstructed images after 100 iterations using CG method (col. 2), and WSD method with optimal weighting factor (col. 3). In all cases, the top image is µ_a, and the bottom image is µ’_s.

Download Full Size | PDF

From the simulations conducted above, we can find that the WSD method can give more accurate location of the embedded target. However, the proper determination of the optimal weighting factor is dependent on the target locations. To further validate the practical applicability of the algorithm, we conducted simulation using multiple targets that are located both near-boundary and deep inside (see Fig. 3). Figure 3(a) shows the original target images of the embedded objects. Figure 3(b) shows the result using conventional CG method after 100 iterations. Figure 3(c) shows the result using one step WSD method with optimal weighting factor (β = 5). We find that this weighting factor can efficiently localize the near-boundary target, but fail to recover the other two objects. Since the near-boundary object exercises more influence on the surface measurements, the calculated gradient vector has large values for the points around the near-boundary objects. To reduce this effect, we first take 30 iterations of conventional CG method, which can greatly decrease the gradient values for the points around the near-boundary object. And then we perform another WSD step to find the second optimal weighting factor, as shown in Fig. 3(d) (β = 8). However, both the conventional CG method and the WSD method failed to recover the deepest object.

 

Fig. 3. Simulation results for multiple targets that are located with different depth. (a) target image. (b) reconstruction result using conventional CG method after 100 iterations. (c) reconstruction result using one step of WSD method and (d) final result using another WSD step.

Download Full Size | PDF

5. Conclusions

In this paper we have presented a method for image reconstruction in DC diffusion-based optical tomography. By imposing a spatial location weighted steepest descent method, we can precisely locate the embedded objects. In addition, the computation burden can also be greatly reduced. We presented the reconstruction of both absorption and scattering coefficient in steady state. In fact, this method can also be applies to other cases such as reconstruction of absorption and diffusion coefficients in frequency domain or time dependent case.

We have made a comparison between the conventional CG method and our new method using the simulated data. The new method can effectively locate the perturbations in absorptions and diffusions. Future work in our laboratory will involve experimental studies for clinical imaging problems.