1. Introduction

Diffuse optical tomography (DOT) is a non-invasive functional imaging modality whose goal is to estimate the optical properties of human tissue. It provides useful physiological information about blood volume and oxygenation and has applications in optical mammography [1, 2] and functional brain imaging [3, 4].

The forward problem in DOT describes the photon propagation in tissue and the inverse problem involves estimating the absorption and scattering coefficients of tissue from light measurements on the surface. The inverse problem is ill posed and highly underdetermined, and thus regularization is typically required through which certain prior assumptions on the solution (e.g., an assumption on its smoothness or a bound on its norm) can be imposed. The most common choice is the Tikhonov-type regularization [5, 6, 7], where the least-square residual is regularized using the L ₂ norm of the unknown parameters. This adjustment aims to reduce high-frequency noise in the reconstructed images; however, it tends to produce an over-smooth solution and performs best when the real solution assumes a Gaussian distribution. Different methods have been developed to overcome these drawbacks. Paulsen and Jiang proposed a non-linear regularization method that minimizes the total-variation norm [8]. The anisotropic diffusion regularization was applied in [9, 10], which provides better performance in preserving geometrical structures such as edges or corners and can readily incorporate prior structural information from other imaging modalities. Pogue et al. employed spatially variant regularization parameters to compensate for the spatial dependence of the contrast and resolution in the reconstruction [6], and several shape-based approaches were proposed in [11, 12] with their performances analyzed using the Cramér-Rao bound [13].

In this paper, we propose an alternative method to improve the spatial resolution of the reconstructed images in DOT. Our method is based on the observation that abnormalities in the domain (e.g., tumors in the breast or activations in the brain) are typically spatially concentrated and thus sparse. Therefore, it would be helpful to incorporate such information in solving the inverse problem. In recent years, the topic of sparse signal representation and estimation has developed in a variety of applications, including image reconstruction and restoration [14], feature selection and classification in machine learning [15, 16], radar signal processing and imaging [17, 18], and biomedical imaging [19, 20, 21, 22, 23, 24]. It was also commonly known as LASSO (Least Absolute Shrinkage and Selection Operator) in the field of regression [25]. In [15], the authors used a Laplacian prior to enforce sparsity and interpreted it using a hierarchical Bayesian method. In [19], a re-weighted minimum norm algorithm called FOCal Underdetermined System Solver (FOCUSS) was proposed, whereby a low resolution estimate of the sparse signal was first obtained and then pruned to a sparse signal representation. Different level-set algorithms were developed to incorporate the sparse nature of the activations in DOT and consider both the support and values of the activated regions [23, 24].

We propose imposing the sparseness of the optical abnormality using the L ₁ norm regularization and solving it using the expectation-maximization (EM) algorithm. The EM algorithm was originally formulated to obtain the maximum likelihood estimation with missing data. In our case, we reformulate the measurement model in such a way that we first consider the real absorption perturbation as a “clean” perturbation corrupted by white Gaussian noise, regard it as a type of missing data, and then estimate the pure perturbation iteratively using the EM algorithm. In particular, we apply a soft-thresholding approach in the maximization step to simplify the computation. We validate our method using 3D simulated examples where the scattering coefficients are assumed to be spatially constant and known and we recover only the absorption coefficient. We compare our approach with standard methods, including the Tikhonov regularization and simultaneous iterative reconstruction technique (SIRT) [26, 27], as well as a level-set algorithm [23]. The level-set method represents the support of the abnormality using a level-set scheme and makes the inverse problem well-posed by exploiting the spatially concentrated nature of the abnormality. We show that our method provides better resolution than the standard Tikhonov-type regularization and is efficient in estimating two closely-spaced abnormalities.

This paper is organized as follows: In Section 2 we briefly review the forward and measurement models in DOT. In Section 3, we propose the image reconstruction method using sparse regularization and provide the solution using the EM algorithm. Numerical examples are given in Section 4, where we present our reconstruction results and compare the performance with other methods. Finally, we offer conclusions in Section 5.

2. Forward and measurement models

Photon propagation in human tissue is mathematically described by the Boltzmann transport equation. It can be simplified under certain assumptions using the diffusion approximation (DA) equation, which can be solved using the finite element method (FEM) [28, 29], finite difference method (FDM) [30], and Born or Rytov approximations [26]. In this paper, we employ the Rytov solution to DA, where the scattered field is assumed to be slowly varying in space. Furthermore, we assume that the scattering coefficient μ _s is spatially constant and known, and we focus solely on the reconstruction of the absorption coefficient μ _a.

Denote N _s the number of light sources, N _d the number of detectors, and N _m = N _s × N _d the number of measurements. Assuming that the domain is divided into N voxels with a constant absorption coefficient μ _ai in the i ^th voxel, the Rytov solution to DA is expressed as [26]

where the superscript “^c” denotes complex values, ϕ ^c ∈ C^N_m the Rytov phase, A ^c ∈ C^{N_m × N} the weighting matrix, and μ ∈ R^N the change of μ _a in each voxel compared with the background. More specifically, we can write μ = [δμ _a1, δμ _a2,…, δμ _aN]^T, where δμ _ai = μ _ai - μ _a0 with μ _a0 denoting the homogeneous background absorption coefficient. Note that Eq. (1) can be used for different simple geometries such as infinite, semi-infinite, or slab, by modifying the elements in A ^c using the method of image sources [31] and extrapolated zero boundary conditions [32, 33]. It can also be easily extended to consider multiple frequencies.

Separating the real and imaginary parts of ϕ ^c and A ^c and letting ϕ = [Rϕ ^cT ,Iϕ ^cT]^T and A = [RA ^cT, IA ^cT]^T, the complex Eq. (1) can be rewritten as a real equation:

with ϕ ∈ R ^N_tot , A ∈ R^N_tot×N, and N _tot = 2 × N _m. Assuming additive Gaussian noise e, the measurement model is given as

where y ∈ R^N_tot denotes the measurement vector. In this paper, we assume e to be zero-mean and spatially uncorrelated with a covariance matrix σ ² I. Namely,

where N (½,Σ) denotes a Gaussian distribution with mean ½ and covariance matrix ∑.

3. Inverse problem and expectation-maximization algorithm

In this section, we first present the inverse problem and interpret it from a regularization perspective. We then introduce the concept of sparsity regularization and formulate the inverse problem using the L ₁ norm of the unknown vector μ. We propose solving the optimization problem using the EM algorithm and give the solution in the end.

3.1. Inverse problem from a regularization perspective

The inverse problem in DOT involves estimating μ from the measurements y. The problem is ill posed and highly underdetermined since N is typically much larger than N _tot. To solve this problem, regularization is usually required, where additional information is introduced about the solution to μ (e.g., an assumption on its smoothness or a bound on its norm). Using regularization, the inverse solution is formulated as

where the first term in brackets denotes the least-square error, g(μ) the regularization function, and γ the regularization parameter controlling the tradeoff between the noise residual and the prior. In the Tikhonov regularization,

which assumes that the solution μ is smooth over the domain. The corresponding solution is given as

where I denotes the identity matrix and A ^# the pseudo-inverse of A. It can be solved either directly or iteratively using, for example, the conjugate gradient (CG) method.

Equation (5) can also be interpreted from a Bayesian point of view. If we assume that the measurements y are corrupted by Gaussian noise and μ is normally distributed with a zero mean, Eq. (5) is equivalently expressed as

where p(y|μ) denotes the likelihood function with ln p(y|μ) ∝ -||y-Aμ||² and p(μ) the prior distribution with ln p(μ) ∝ -||μ||² ₂. Therefore, the solution using the Tikhonov regularization is equivalent to the Bayesian solution assuming a Gaussian prior on μ.

3.2. Introduction of the sparsity regularization

Note that in the image reconstruction of DOT, the distribution of μ (i.e., the change of the absorption coefficient) is not necessarily Gaussian in 3D. In particular, it appears in a sparse format where most of its elements are zeros (corresponding to the unchanged background), and only a few of them are nonzero (corresponding to the abnormalities). Based on this fact, we propose incorporating this sparseness prior into the image regularization formulation to improve the reconstructed resolution.

The ideal measure of the sparseness of a vector x is its L ₀ norm: ||x||₀, which is the number of non-zero entries in x. However, the minimization/maximization involving ||x||₀ is NP-hard (Nondeterministic Polynomial-time hard) and can be solved only using a combinatorial optimization approach. Instead, it is usual to use the L ₁ norm ||x||₁ as an approximated measure, where

It has been shown that if the solution μ is sparse enough, the solution to the L ₁ norm minimization is equivalent to minimizing the L ₀ norm [34].

Under the measurement model (3) and the sparsity assumption on μ, our inverse problem is expressed as

where ζ ² is a parameter on the noise level. It can be equivalently expressed in an unconstrained form as

Equation (11) is a convex problem but it is non-differentiable. It can be solved using linear or quadratic programming, but in this paper we propose solving it using the expectation-maximization (EM) algorithm. By decomposing measurement model (3), the EM method takes advantage of certain structures and results in computationally simple solutions.

3.3. Equivalent model and EM algorithm

We propose solving the optimization problem in (11) using the EM algorithm. The EM algorithm is a general method to obtain the maximum penalized log-likelihood estimator (MPLE) by introducing missing data and maximizing the complete penalized log-likelihood. In our case, the MPLE is expressed as

which provides the same solution as (11). To apply EM, we first reformulate the measurement model (3) so that

where α is a positive constant. This reformulation is equivalent to expressing

where e ₁ and e ₂ are independent, such that

Note that in order for the matrix σ ² I - α ² AA^T to be positive definite (so that it is a valid covariance matrix), we must have

where β ₁ denotes the largest eigenvalue of AA^T.

The decomposition in (14) introduces a hidden variable x, which is the noisy version of the true absorption perturbation μ. By regarding x as the missing data, we can estimate the real μ using the EM algorithm. The EM algorithm produce a sequence of estimates μ̂^(k),k = 1,2,… by alternating two steps (see below) until some stopping criterion is met.

3.4. Some comments

A few comments deserve mention here regarding the aforementioned image reconstruction algorithm.

3.4.1. On the model (14)

By introducing the hidden variable x, we decompose the inverse problem (the mapping of y → μ) into two parts: The first part can be interpreted from the typical image reconstruction point of view, which removes the effects of measurement noise (i.e., the mapping of y → x); the second part can be regarded as a denoising procedure (the mapping of x → μ). This formulation can help improve the resolution of the reconstructed images, as we will show in Section 4.

3.4.2. On the convergence of the EM algorithm

Based on the information matrices, it has been shown in [36] that each iteration of EM is guaranteed to increase the penalized log-likelihood function, namely,

However, it does not mean that the sequence will converge to the maximum likelihood estimator of μ. If the distribution is multimodal, Wu showed in [37] that the EM algorithm is guaranteed to converge to a stationary point (i.e., the local maximum or saddle point) provided that the Q function and the penalty term are continuous in both μ and μ̂. The convergence performance also depends highly on the choice of the starting value μ ⁽¹⁾. To escape from the stationary point, we can use several different random initial points and also incorporate the prior information regarding the abnormality distribution.

3.4.3. On the soft-threshold method

In the M-step, we employed the soft-threshold method to implement the maximization. This approach is commonly used in the wavelet denoising scenario for image processing [15, 38, 39, 40], where the goal is to solve the problem of argmin_x||y-Ax||² ₂ + γ||x||₁ with A being orthogonal and x sparse. Certainly, other algorithms, such as the subgradient method and the interior-point method can also be used to solve (23), but we choose the thresholding algorithm for the sake of simplicity. From (20) and (24), we can see that our method has a computational complexity of O(N ²) for the update in the E-step and O(N) in the M-step.

3.4.4. On the choice of the parameters

We need to consider the following issues when determining the model parameters: (1) the convergence rate of the EM algorithm and (2) the effect of the sparsity regularization. The first factor is affected by α. According to the theory of the convergence rate of EM [36], α should be made as large as possible for a faster convergence. In our case, since we need to satisfy the condition α ² ≤ σ ²/β ₁ (see Subsection 3.3) for the validity of the EM model (14), we should choose α to be as close to σ/√β ₁ as possible.

The effect of the sparsity regularization is controlled by α, γ, and σ as shown in Eqs. (13) and (22). Once σ is given and a selected based on (18), it is obvious that the higher the γ the sparser the reconstructed image. According to the first-order optimality condition in the convex optimization, Kim et al. showed in [41] that for the general problem of argmin_x ||Ax - y||² ₂ + γ||x||₁, an upper bound on the useful range of γ is given as (γ)_max = ||2A^Ty||_∞, where ||x||_∞ = max_i{x_i}. For σ that is higher than (γ)_max, the estimated x̂ would have all its elements be zero. After applying this result to Eqs. (13) and (22), we derive that σ, γ, and α should satisfy the following equations:

Note that Eqs. (26) and (27) provide only an upper bound on the model parameters. Their main function is to avoid selecting a γ that is too large for obtaining a meaningful result. The right-hand side of Eq. (27) can be easily computed once the model is determined. Regarding the right-hand side of Eq. (26), although x denotes the noisy version of the true perturbation μ, its maximum can be estimated given certain prior information. For example, according to [26], the perturbation in the abnormality’s absorption value ranges from 0.02 cm^-1 to 0.3 cm^-1 for most biological tissue, depending on the tissue type. Based on this information, we would be able to obtain a reasonable estimate of the bound on α ² γ. In practice, we first determine the range of α, γ, and σ based on Eqs. (18), (26), and (27) and then select their values experimentally based on the signal-to-noise ratio to obtain the optimal reconstruction results.

4. Numerical examples

In this section, we provide numerical examples to illustrate the reconstruction results using our sparsity regularization and compare it with other methods. We considered a 3D transmission geometry as shown in Fig. 1, where the origin is at the center of the bottom surface (z = 0 cm) and the sides are of length 8 cm, 8 cm, and 6 cm along the x, y, and z directions, respectively. We placed 25 sources on the bottom surface at 1.75 cm intervals and 25 detectors on the top surface (z = 6 cm) at 1.5 cm intervals. We chose a source modulation frequency of 200 MHz and wavelength λ = 750 nm.

We chose the background optical parameters as μ _a0 = 0.05 cm^-1, μ′_s0 = 9.5 cm^-1, and the speed of light c = 22 cm/ns [42]. We divided the domain in to small voxels with size 4 × 4 × 5 mm³ and set up the forward model (2) using the method of image sources with the extrapolated zero boundary condition. We added Gaussian noise with standard deviation σ = 0.01 and performed image reconstruction using four methods: (a) our ZL ₁-EM approach, (b) the level-set algorithm proposed in [23], (c) the Tikhonov regularization, and (d) the simultaneous iterative reconstruction technique (SIRT). In the following, we first demonstrate the reconstruction results of four methods assuming one absorbing abnormality and then discuss their performances in separating two closely spaced abnormalities.

 

Fig. 1. Simulation setup. The domain is composed of an 8 × 8 × 6 cm³ cube, where twenty-five sources (circles) are placed on the bottom surface (z = 0 cm) and twenty-five detectors (black dots) on the top (z = 6 cm). The range along the x and y axes is [-4,4] cm.

Download Full Size | PDF

4.1. Image reconstruction results

We assumed that there is one spherical absorption abnormality centered at [- 1.5, 1.25, 2.9] cm with radius R = 1 cm and absorption perturbation δμ _a = 0.2 cm^-1. The abnormality covers 54 voxels out of a total of 4,800 voxels in the domain. The original μ _a distribution is shown in Fig. 2, where each small image shows the cross-section layer at different z values at 0.5 cm intervals.

For our L ₁-EM approach, we selected α = 3.5 × 10^-4, γ= 10⁴, and stopped the EM iteration when |μ̂^(k+1) - μ̂^(k)| ≤ 10^-3. Under the current setup, the solution converges after 400 iterations, taking about 45 seconds using Matlab 7 on a PC with Pentium 4 2.6 GHz CPU and 1 G of RAM. For the level-set algorithm, we used κ = 0.4, λ = 0, μ = 2 × 10^-4, and stopped the iteration when the relative change e(k) is less than 10^-8; see [23] for the notation and implementation. We selected γ = 0.05 for the Tikhonov regularization (see Eq. (7)) and terminated the SIRT after 800 iterations.

The results are shown in Fig. 3. We can see that all of the four methods can recover the center location of the abnormality correctly. However, clear differences in performance can also be observed. More specifically,

From the above example, we can see that our method can provide image reconstruction result with fairly good resolution. Combined with other methods using Tikhonov regularization, it can be useful in determining the range of the activation with more accuracy.

4.2. Performance analysis

In this subsection, we study the performance of our method for separating two closely spaced absorbing abnormalities. We assumed two spherical absorptive abnormalities with radius R = 0.75 cm and δμ _a = 0.2 cm^-1. We considered the following two cases: In the first case, the centers of these two spheres are at [-1.5, 0.8, 2] cm and [1.5, -0.8, 4] cm, with a distance of 4 cm between the two centers and 2.5 cm between the two closest points on the spheres. In the second case, the centers are moved to [-1, 0.5, 2.5] cm and [1, -0.5, 3.5] cm, and the distance becomes 2.45 cm between centers and 1 cm between the closest points. The original δμ _a distributions are shown in Fig. 4.

The reconstruction results are shown in Fig. 5 and Fig. 6, respectively. We can see that for the first case, all four methods can separate the two abnormalities with different levels of accuracy. In a similar fashion as the results shown in Fig. 3, the L ₁-EM method and level-set algorithm reconstruct the problem with the best resolution, and the Tikhonov regularization and SIRT method show wider side lobes than the true distribution. When the two abnormalities get closer, the difference in the reconstruction performance becomes more obvious. We can tell that the SIRT method fails to differentiate the two abnormalities: There is only one big sphere showing around the midpoint along the two abnormalities. Both the L ₁-EM and level-set algorithm can separate the abnormalities well. However, the level-set algorithm underestimates the size in this case. This might be due to the small size of the activations, leading to a poor approximation of curve evolution; the discretization of the problem on a finer grid may yield better estimates. The result using the Tikhonov regularization again shows bigger side lobes.

 

Fig. 3. Reconstructed δμ _a distributions in cm^-1, assuming only one absorbing abnormality as shown in Fig. 2. (a) Using the sparse regularization with EM algorithm; (b) using the level-set algorithm; (c) using the Tikhonov regularization; (d) using SIRT. Small images show the cross-section layers at different z values at 0.5 cm intervals.

Download Full Size | PDF

 

Fig. 4. Original δμ _a distributions in cm^-1, assuming two spherical absorptive abnormalities.

Download Full Size | PDF

 

Fig. 5. Reconstructed δμ _a distributions in cm^-1, assuming two spherical absorptive abnormalities as shown in Fig. 4(a). (a) Using the sparse regularization with EM algorithm; (b) using the level-set algorithm; (c) using the Tikhonov regularization; (d) using SIRT.

Download Full Size | PDF

5. Conclusions

We proposed an image reconstruction method for diffuse optical tomography by introducing the sparsity regularization. We formulated the inverse problem by regularizing the L ₁ norm of the unknown absorption coefficients and solved the optimization problem using the expectation-maximization (EM) method with a soft-threshold approach. We compared our method with other image reconstruction approaches and showed its efficiency in improving the reconstruction resolution.

In the current work, we applied our method to a transmission geometry, using the forward model based on the Rytov approximation to the diffuse approximation equation. This framework is useful for breast tumor detection. Our future work includes applying the proposed methods to the FEM forward model with a spherical domain shape and testing it using the real optical measurements from brain imaging. In order to apply our method to a nonlinear forward solver like FEM, a linearized model using the Taylor series needs to be obtained first for the applicability of the decomposition (14). More details on model linearization and computation of the Jacobian matrix for the FEM model can be found in [28, 29].

 

Fig. 6. Reconstructed δμ _a distributions in cm^-1, assuming two spherical absorptive abnormalities as shown in Fig. 4(b). (a) Using the sparse regularization with EM algorithm; (b) using the level-set algorithm; (c) using the Tikhonov regularization; (d) using SIRT.

Download Full Size | PDF