Open Access
Add to BibSonomyAbstract
Diffuse optical tomography (DOT) allows tomographic (3D), non-invasive reconstructions of tissue optical properties for biomedical applications. Severe under-sampling is a common problem in DOT which leads to image artifacts. A large number of measurements is needed in order to minimize these artifacts. In this work, we introduce a compressed sensing (CS) framework for DOT which enables improved reconstructions with under-sampled data. The CS framework uses a sparsifying basis, ℓ1-regularization and random sampling to reduce the number of measurements that are needed to achieve a certain accuracy. We demonstrate the utility of the CS framework using numerical simulations. The CS results show improved DOT results in comparison to “traditional” linear reconstruction methods based on singular-value decomposition (SVD) with ℓ2-regularization and with regular and random sampling. Furthermore, CS is shown to be more robust against the reduction of measurements in comparison to the other methods. Potential benefits and shortcomings of the CS approach in the context of DOT are discussed.
© 2010 Optical Society of America
1. Introduction
Diffuse optical tomography (DOT) has received significant interest due to its potential to non-invasively image tissue function. For example, significant advances in optical mammography were reported where benign and malignant lesions could be identified and their changes over time in response to cancer therapy could be monitored [1, 2]. In DOT, near-infrared (NIR) light is delivered and detected to the interrogated tissue by sources and detectors placed over the tissue surface. These measurements are then utilized in image reconstruction algorithms to obtain two or three-dimenstional (3D) maps of the tissue optical properties [3,4,5,2]. However, since the general problem of image reconstruction in turbid media is ill-posed and often under-determined, complex image reconstruction algorithms are needed [3]. An overview of various inverse approaches to the optical tomography problem can be found in references [5] and [6].
One such, common, approach to the linearized inverse problem is to use a singular value decomposition (SVD) based algorithm with the Tikhonov regularization (ℓ2-norm) [7, 4, 8, 9]. Several improvements such as a spatially-dependent regularization [10], have been proposed to overcome its shortcomings. However, it is still desirable to improve the fidelity of the reconstructions despite under-sampling.
Compressed sensing (CS) or compressive sampling [11,12] has recently emerged as an attractive method for signal recovery from a relatively small number of random linear combinations of the signal samples -far fewer than dictated by the Shannon-Nyquist theorem [13, 14]. The key concept behind CS is to exploit the prior knowledge that a signal is sparse or compressible in a known transform domain (i.e., a large number of the unknown coefficients are close to or equal to zero) and look for the sparsest solution [15]. It turns out that it is better to use the ℓ1-norm instead of the ℓ2-norm for the penalty term since it has been proven that minimal the ℓ1-norm signifies the sparsest solution [15]. Since the introduction of CS, a diverse spectrum of signals were shown to be suitable for use with the CS-framework [16, 17, 18, 19].
For medical image reconstruction, CS has been used in magnetic resonance imaging [20], interior tomography [21], photo-acoustic tomography [22, 23, 24], computed tomography [25], terahertz imaging [26] and in holography [27]. A formulation of DOT reconstruction problem as sparse representation was previously presented [28] which satisfies only one of the CS steps. A similar formulation with sparsity regularization was also discussed and an improvement in spatial resolution was demonstrated, again taking a step towards a CS framework [29]. Furthermore, a related approach, ℓ1-penalized sparse constructions, has been utilized for fluorescent DOT [30] and bioluminescence tomography [31]. While these studies have demonstrated the potential of the CS-framework for general tomographic approaches, the CS-framework has not yet been demonstrated for DOT of endogenous contrast.
In this work, we describe an implementation of the CS-framework in the context of DOT of endogenous contrast [32]. We propose CS as an efficient methodology for reconstructing 3-D images of the optical properties of tissue from a relatively small number of measurements. The quality of the reconstructed images is compared with “traditional” SVD reconstructions using several image quality metrics. Furthermore, in order to investigate the effect of “random” sampling in the “traditional” SVD approach, both regularly under-sampled and randomly under-sampled cases are tested. By starting with a relatively large number of measurements and working our way down, we show that as the number of measurements is reduced, the CS based algorithm exhibits higher fidelity results compared to both randomly and regularly under-sampled SVD based reconstructions.
The paper is organized as follows: The DOT linear inverse problem and the implementation of the compressed sensing (CS) theory, are presented in Section 2. The simulation of the DOT reconstruction problem using both CS and SVD is described in Section 3.1. Several image quality metrics, used to evaluate the quality of the reconstructed images, are also presented. The reconstruction performance of the proposed methodology is assessed in Section 4 and discussed in Section 5.
2. Theory
2.1. The Forward and Inverse Problem in DOT
An inhomogeneous photon diffusion equation (IPDE) is used to model the propagation of photons through turbid media, such as biological tissue. In order to make the problem of reconstructing a tomographic image of the unknown optical properties more tractable, the IPDE is often linearized. The linearization involves expanding the IPDE where we seek to construct an image of small “perturbations” from a homogeneous “background” or “baseline” medium. In this formalism, the linearized DOT methods aim at reconstructing a three-dimensional (3-D) image of the “perturbed” optical properties [5]. This linearization is often achieved by using either the Born or the Rytov approximation (RA) [33].
In this work, we focus on the Rytov approximation, since it has been shown to approximate larger perturbations compared to the Born approximation [34]. In the presence of inhomogeneities, the IPDE is solved [33, 34] by approximating the “measured” photon fluence, at position r, as Φ(r) = Φ0(r)exp(Φsc(r)). Here, Φ0(r) denotes the incident wave and Φsc(r) the scattered wave due to an inhomogeneity. Φ(r) and Φ0(r) are calculated by employing either analytical or numerical solutions of the IPDE (forward problem).
The objective is to estimate Φsc(r) based on measured data (at position r) and associate it to the spatial distribution of the optical properties of the interrogated volume, in our case, to the relative absorption coefficient, Δμa(rV) = μa,I – μa,B (rV), where μa,I, μa,B denote the (absolute) absorption coefficients inside the inhomogeneity and in the background region (surrounding medium), respectively. rV represents a position in the 3-D space, where the reconstruction method is performed. For notational simplicity, the position vectors (rV or r, etc.) will be omitted hereafter, unless needed for clarity.
If the medium is further discretized into volume elements (voxels, each of volume V), then we obtain the following matrix formulation for the DOT image reconstruction problem:
The problem of the reconstruction of the relative absorption coefficient Δμa (optical image) is often referred to, as the linear inverse problem, where J is the Jacobian matrix or the weight matrix. Following this RA formulation, the Jacobian can be written as [34] where v is the speed of light in the medium, D0 is the background diffusion coefficient and G denotes a Green’s function solution. The indices are i = 1, .., ND, j = 1, .., NS, k = 1, ..., N and l = 1, .., NS × ND, where NS, ND, N, denote the number of sources, detectors, and voxels, respectively. Finally, rj, ri, rk refer, correspondingly, to the positions of the sources, detectors, and voxels. To construct the above Jacobian, appropriate Green’s function solutions must be provided, based on the corresponding boundary conditions. The Jacobian relates the performed measurements with the optical image that needs to be reconstructed. Each row (lth) of the Jacobian corresponds to a given source (jth) and detector (ith) pair. Each column (kth) corresponds to a given voxel.The above linear inverse problem is ill-conditioned and often regularization is required to stabilize the solution for the reconstructed image Δμa, by formulating it as a minimization problem, as follows:
where λ denotes the regularization parameter and ||.||2 is the ℓ2-norm. Singular value decomposition (SVD) approaches have been frequently employed to solve the above problem by inverting the Jacobian [8,35]. The regularization parameter is commonly optimized using an L–curve analysis where the trade off between error in the minimization (residual) and the image noise is optimized [36].2.2. Compressed Sensing
The fundamental concept behind compressed sensing (CS), provided that specific theoretical conditions hold [11, 12] (see below), is that the unknown signal can be recovered with good accuracy even when sampled at a rate significantly below the Shannon-Nyquist sampling rate [13, 14]. The Shannon-Nyquist theorem has been widely considered, so far, as a lower limit to the number of samples that are required for a signal to be accurately reconstructed. However, recent developments in the signal recovery theory led to a change in the acquisition paradigms in various research fields [17]. This resulted in the drastic reduction of the required number of measurements [37, 25].
To further explain the CS framework, we consider an image vector x ∈ Rn (of n unknowns) and a measurement vector y ∈ Rm (of m measurements). We assume that they are linearly related by a measurement matrix Φ, such that, y = Φx. The CS theory asserts that the signal (image) x can be efficiently recovered provided that the following conditions hold:
- The image vector x is k – sparse or can be sparsified using a known orthogonal transformation T (e.g. the discrete Fourier transform), such that x = Tx̄ and x̄ is k – sparse. k – sparsity implies that x̄ can be represented by a linear combination of only k basis vectors (the remaining n – k transformation coefficients are zero). The locations of the k non-zero elements are not known a priori, as the image vector x itself is unknown. The fulfilment of the condition m ≥ k log(n) guarantees that the signal can be exactly recovered with high probability.
- The product matrix ΦT is known as the CS-matrix and fulfills the incoherence property [12], i.e. the measurement matrix Φ must be incoherent with respect to the orthogonal transformation T. Furthermore, linear combinations of the columns of the CS-matrix should act like noise and be linearly independent. It is known that ideally distributed random matrices are maximally incoherent with respect to any basis.
- The recovery of the unknown signal (image) is performed using a nonlinear reconstruction scheme that enhances sparsity. A fundamental such technique involves the use of a ℓ1-regularized, nonlinear, constrained minimization for the unknown sparse image x̄. This can be formulated as follows [17, 38]: where ||x̄||1 denotes the ℓ1-norm of x̄, i.e. ||x̄||1 = Σi|x̄i|.
We will now discuss how the CS theory could be applied to the DOT reconstruction problem.
2.3. Application of CS to the DOT Reconstruction Problem
The implementation of the CS theory within the framework of DOT, involves the construction of a linear convex optimization problem, i.e.:
where Δμa is a vector of unknown absorption coefficients at each voxel. It also corresponds to the 3D image of the spatial distribution of the relative absorption coefficient, that needs to be recovered. Furthermore, Δμa is assumed to have a sparse representation in a given basis T, such that Δμa = TΔμ̄a. If this is true, then based on Equation (4), Δμa can be recovered by solving a minimization problem of the form [17, 38]The above formulation denotes an unconstrained minimization problem, in which, the regularization parameter Λ should be selected (optimized) such that the error between the residual and the image is minimized. As described in the following Section 3.3, this is calculated in a heuristic manner.
As described in the previous section, certain conditions must be satisfied to apply the theory of CS to the DOT reconstruction problem. First, we investigate the sparsity requirement. For the linearized DOT problem, sparsity can, generally, be pre-assumed since the linearization itself assumes that the perturbation from a homogeneous background or reference medium is relatively small in volume and contrast [34, 33]. However, this may not be true in biological media unless a very good reference medium is available. Therefore, we have used the Fourier basis as a transformation basis in order to reinforce sparsity. Alternative bases have been also tested, such as the Haar/wavelets basis, but no significant differences were observed.
The second property is the assumption of incoherence between the transformation matrix T and J. A coherence function, μ′, is formed as the maximum value of the inner products between columns selected from the aforementioned matrices [17, 37], i.e.
where i and j denote columns within the corresponding matrices. A small value of μ′ implies incoherence and has been used in applications of CS to magnetic resonance imaging (MRI) [20] and photoacoustic tomography (PAT) [24, 23]. A more generalized measure of coherence, not commonly utilized, is known as the cumulative coherence function (CCF) and is formulated as: where m denotes the order of the CCF and Ω is the set of all possible subsets of m columns selected from the Jacobian, whose size is (NS · ND) × N (see Equation (2)). Furthermore, k ranges over a specific subset of Ω at a given order and j defines columns of the transformation matrix for all possible subsets [39]. It has been shown, that a slowly increasing CCF versus order m signifies incoherence [39]. It should be noted, that the coherence function μ′ defined in Equation (7), represents an order-1 CCF (m = 1).The third requirement is the application of a nonlinear optimization scheme based on ℓ1-regularization to recover the signal (Δμa). The convex minimization problem described in Equation (6), is known as sparse ℓ1-regularization in the context of CS [11, 12], since the ℓ1-norm is employed. Various other schemes have been also proposed to solve this minimization problem [17,20,40] but ℓ1-norm minimization is accepted to be the most applicable one since it was shown that the minimum ℓ1-norm solution is also the sparsest solution [15].
3. Methods
3.1. Numerical Simulations of Forward Data
A spherical inhomogeneity of radius 0.5 cm was considered at a depth of 1.5 cm from the plane of source-detectors. The forward data was generated using the analytical solution of a sphere embedded in an infinite medium in the frequency (w=70 MHz) domain [41]. The center of the sphere was at (x, y, z) = (−1.0, 1.0, 1.5). A 5 × 5 rectangular array of sources (S) and detectors (D) was placed at z = 0, resulting in 625 pairs of sources and detectors (S/D). This particular arrangement of sources and detectors (basic set-up) is then used as the basis for our measurement reduction schemes (see below). The absorption and reduced scattering coefficient of the background medium were taken to be μa,B = 0.02 cm−1 and μ′s = 8 cm−1, respectively. A wide range of values was assumed for the absorption coefficient μa,I within the spherical inhomogeneity, i.e., μa,I = {0.02, 0.06, 0.08, 0.14, 0.18} cm−1. A sketch illustrating the geometry is shown in Fig. 1.
Fig. 1 Illustrating the geometry of the simulated optical domain (infinite-medium). For the field-of-view, 3-D volume of size 4 × 4 × 3 cm3 and voxel dimensions 0.4 × 0.4 × 0.6 cm3 is considered. A spherical inhomogeneity is embedded in the medium at z = −1.5 cm. A 5 × 5 array of sources and detectors is placed on a plane at z = 0 cm.
For image reconstruction, a 3-D volume (Lx × Ly × Lz) of size 4 × 4 × 3 cm3 and voxel dimensions 0.4 × 0.4 × 0.6 cm3 (500 voxels) was simulated as the optical domain (or field-of-view). A Jacobian of size 500 × 625 was constructed according to Equation (2) using the Green’s functions for an infinite medium.
The images of the investigated volume were reconstructed based on the implementation of the CS theory, as described in Section 2.3. The full set, as well as, a reduced number of measurements was considered in order to assess the performance of the proposed technique in response to a reduction of the number of measurements. The quality of the reconstructed images was compared with that of the corresponding SVD reconstructions.
As described previously, CS relies on random sampling of signals. For this reason, a number of S/D pairs was randomly removed every time, by ignoring random rows from the Jacobian. The percentage of the removed number of measurements that were investigated are ∼ 0, 15, 25, 40, 50, 70, 80, 85, 90, 92, 95, 99% out of the original 625 S/D pairs. In a similar manner, using random sampling, the conventional SVD approach was also tested. This method will be simply referred to as “SVD”.
Since random removals of S/D pairs in SVD might imply a CS-like improvement in its performance, a “regular”-sampling method was also employed for the SVD reconstructions. This typically involves a regular grid of S/D pairs located on the measurement plane over a given area. The spacing between the S/D pairs on the grid is each time modified (increased) in a regular manner, in order to obtain a reduced set of measurement data. This method of regular sampling in SVD will be referred to, hereafter, as “SVD-reg”. Based on this method, regular grids were constructed using 625, 300, 150, 144, 72, 50, 36, 24, 12 and 4 S/D pairs, corresponding to the following percentages of removed S/D pairs: ∼ 0, 52, 66, 77, 89, 92, 95, 97, 99, 99.4%, respectively. It should be noted, that SVD-reg samples the S/D plane in a regular manner, but not the 3-D volume necessarily. The latter could be achieved by modifying the number of neighbours used, but such elaboration is out of the scope of this work [42].
Figure 2 shows the positions of the S/D pairs on the z = 0 plane for (a) the entire set of S/D pairs (0% removal), (b) the remaining S/D pairs after ∼ 90% “random” sampling (removal) for CS and SVD, and (c) the remaining S/D pairs after ∼ 90% “regular” sampling for SVD-reg is shown in Fig. 2(c). Note that the sample shown in (b) is a special case where all S/D pairs from the remaining sources and detectors were utilized. We further illustrate the random removal case for a 3 source (S1, S2, S3) and 2 detector (D1, D2) experiment. We start with all possible pairs, i.e. S1D1, S1D2, S2D1, S2D2, S3D1, S3D2. If we remove 50% of these pairs, we may end up with S1D2, S2D2, S3D1 or S1D1, S3D1, S3D2 or other possible combinations in a random fashion. Each row of the Jacobian corresponds to one such pair.
Fig. 2 (a) Full set of S/D pairs (0% removal), (b) remaining S/D pairs after ∼ 90% “random” sampling (removal) for CS and SVD, (c) remaining S/D pairs after ∼ 90% “regular” sampling for SVD-reg.
3.2. Performance Evaluation of the Reconstructed Images
To evaluate the performance of the reconstructed images several image quality metrics were used. Two such metrics are the observed contrast (Co) and the contrast-to-noise ratio (CNR) [43, 44], which have been modified to fit the requirements of the 3-D DOT reconstruction problem. The observed contrast (Co) is defined as
where Δμ̂a,I, Δμ̂a,B denote the spatially averaged reconstructed relative absorption inside the spherical inhomogeneity and in a similar-sized 3-D region in the background (surrounding medium), respectively. The background region is taken to be as far as possible away from the sphere (along the diagonal on the (x, y)-plane), so that the spreading of the reconstructed spherical inhomogeneity does not interfere with it.The contrast-to-noise ratio (CNR), which is a measure of the detectability of an inhomogeneity given background “noise” or image artifacts is defined as [43, 44]:
where σI and σB denote the standard deviations of the reconstructed relative absorption inside and outside the inhomogeneity (background), respectively, in a similar manner as described previously. In general, the CNR represents the combined effects of the reconstruction errors and the presence of measurement noise. It should be noted that in this work, no measurement noise has been added to the forward data, thus, the calculated CNR can be considered as an upper limit of what can be achieved in practice.We also calculate the normalized root mean square error (nRMSE) compared to the expected/simulated value in the 3-D space, as [23]
where Δμa,true(rV) is the true value of the relative absorption coefficient assumed at position rV, N denotes the total number of voxels and , the maximum and minimum values, correspondingly, of the reconstructed relative absorption within the investigated volume.Finally, the localization error (LE) was calculated, as the absolute 3-D distance between the position rmax where the maximum reconstructed relative absorption ( ) occurs and the true position rtrue,I of the center of the simulated inhomogeneity, i.e.
Overall, these metrics allow us to evaluate the performance of the three approaches we test from multiple angles highlighting differences in the optimization of the regularization, the presence of image artifacts and the quantification capabilities of each approach.
3.3. Implementation of the CS Algorithm for DOT
The re-formulation of the CS framework for DOT was presented in Section 2.3. Both the forward and inverse parts were developed in GNU R [45, 46] using standard libraries. In order to solve Equation 6, we have implemented an ℓ1 penalized conjugate gradient algorithm as described in reference [20]. SVD-based reconstructions were also carried out using standard libraries in [45, 46]. Samples of the codes and a detailed write-up of the functions are available on request [46].
As already mentioned, a heuristic approach was considered to optimize the conjugate gradient algorithm and the regularization parameter Λ (see Equation 6). We have monitored the objective function (the argument in Equation 6 over a large number of conjugate gradient (CG) iterations until it has reached a plateau. The Λ value, Λ = 5 × 10−10, that led to the minimal objective function over all the iterations was used in the simulations. An L-curve analysis was also applied (valid for ℓ2-regularization) [47], but it was not finally adopted since it did not perform efficiently (a uniform distribution of Λ can lead to an uneven sampling of the L-curve). Semi-automated schemes for calculating the error have been proposed in the literature as an alternative but these suggestions did not behave well for our particular problem [23]. A more complex parameter selection procedure may be necessary for the optimization of the ℓ1-regularization [47] which is beyond the scope of this work. We note that the formalism developed in reference [47] suggest that it is possible to move beyond a heuristic approach in the future.
4. Results
We have developed a procedure for linearized DOT image reconstruction within the CS-framework in Section 2.3. In order to ensure that we satisfy the CS requirements described in Section 2.2, the first step that we sought was to transform the unknown DOT image (Δμa) into a sparse one (Δμ̄a). This was achieved by using the Fourier basis as a transformation basis (T) where the resulting product matrix of the Jacobian and T is the “CS matrix”. We have argued that an important property of this matrix is that it should be incoherent with respect to the orthogonal transformation T with linearly independent columns. To test this condition, we have proposed to use the cumulative coherence function (CCF), defined in Equation (8). Due to the high computational complexity of the above formula, the CCF was calculated for an order of up to mmax = 16. Likewise, the number of columns of the Jacobian J was taken to be N = 16 (see Equation (2)). The normalized CCF is shown in Fig. 3. The slowly increasing trend, which is clearly identified, validates the incoherence property of the CS matrix JT [39]. More elaborated techniques for the construction of an efficient CS matrix have been proposed [48], however, they are beyond the scope of this paper. Other conditions such as ℓ1-regularization and random sampling were also implemented as described in Section 2.3.
Fig. 3 Cumulative coherence function (CCF) between the Jacobian (J) and the discrete Fourier transform (T) for several scaled orders m/N (m = 2 – 15, N = 16). The CCF was scaled with respect to order m = 1, i.e. (CCF(m/N) – CCF(1/N))/CCF(m/N).
We now demonstrate the performance of the CS algorithm and compare it to those using SVD and SVD-reg by using the image quality metrics, as described in Equations (9)–(12). Fig. 4 shows a set of sample images that were obtained for a case of ∼50% removed S/D pairs. They are representative of the other sets that were tested and repeated runs have confirmed that the results for CS are independent of the specific S/D pairs that were removed. The “forward” or expected images are shown in the top row for μa,I = 0.08 cm−1. Five layers (slices) of 0.6 cm steps were considered along the z-axis (the S/D plane was taken to be at z = 0 cm). The DOT reconstructed images using SVD-reg, SVD, and CS, are shown in the second-to-forth rows. It can be observed, that in the CS images, the simulated inhomogeneity was accurately reconstructed at z = −1.5 cm. However, in the SVD and SVD-reg reconstructed images, the inhomogeneity is reconstructed at z = −0.9 cm, i.e., a shift towards the S/D plane, is exhibited. Such a shift is known to occur in regularized and linearized DOT reconstructions and can be attributable to the highly increased sensitivity of the Jacobian near the optodes and several regularization schemes, including spatially-variant regularization [10], have been proposed to reduce this effect. Other iterative, linear inversion algorithms may produce different systematic artifacts [34, 3, 4, 5]. The ability of the proposed CS implementation to recover the simulated inhomogeneity on the correct z-layer might be due to the effects of random sampling and ℓ1-regularization that are involved in the CS-framework which relax the sampling requirements.
Fig. 4 From top to bottom: “Forward” imaging data, SVD-reg, SVD, and CS reconstructed images for a simulated inhomogeneity of μa,I = 0.08 cm−1 and a case of ∼ 50% reduced S/D pairs. Five layers along the z-axis were considered, ranging from z = −0.3 cm (S/D plane) to z = −2.7 cm (from left to right). The inhomogeneity was taken to be localized at z = −1.5 cm (layer 3). Please note the use of differences in colorbars in each row.
To evaluate the performance of the reconstructed images, the observed contrast (Co, Equation (9)) was initially investigated, for μa,I = 0.08 cm−1. The results are shown in Fig. 5(a) for the investigated range of removed S/D pairs, for all three methods (CS, SVD, SVD-reg). The contrast was normalized according to that calculated by the full set of measurements (0% removal) for each case. A slower decrease of the contrast, indicating a higher robustness, is observed for the proposed CS method, compared to both SVD and SVD-reg, as the number of measurements is reduced. Furthermore, the superiority of SVD with respect to SVD-reg should be pointed out, which can be attributed to the effects of “random sampling”. As illustrated in Fig. 5(b), if we compare two cases, i.e., 0% and 92% removals, for a range of the inhomogeneity μa,I values, we see that CS with 92% of pairs removed performs nearly as well as both SVD and SVD-reg without any removals for the whole range.
Fig. 5 (a) Normalized observed contrast versus the investigated range of removed S/D pairs for CS, SVD and SVD-reg. The contrast was normalized according to that for 0% removal and the absorption coefficient of the inhomogeneity was μa,I = 0.08 cm−1. (b) Observed contrast versus the true absorption coefficient of the inhomogeneity μa,I, for 0% and 92% S/D removal.
The normalized contrast-to-noise ratio (CNR) according to Equation (10), is shown in Fig. 6 for μa,I = 0.08 cm−1. In this case, all three methods exhibit similar behaviour with the number of removed S/D pairs. As we will discuss further in the following section, this may be due to differences in the regularization parameters that were employed and their differing effects on the contrast and noise.
Fig. 6 Normalized contrast-to-noise ratio (CNR) versus the range of removed S/D pairs. The CNR was normalized according to that for 0% removal. The absorption coefficient of the inhomogeneity was assumed to be μa,I = 0.08 cm−1.
Figure 7(a) shows the normalized root mean square error (nRMSE, Equation (11)) versus the investigated range of removed S/D pairs, for μa,I = 0.08 cm−1. An improvement of approximately ∼ 10% in the calculated nRMSE can be clearly seen for the CS method compared to SVD and SVD-reg. This can be mainly attributable to the shift of the reconstructed inhomogeneity towards the S/D plane that is observed in SVD methods, as described in the beginning of this section (see also Fig. 4). The localization shift is more evident in Fig. 7(b), where the localization error (LE, Equation (12)) between the center of the simulated inhomogeneity and that of the reconstructed inhomogeneity, is depicted for all three methods. The LE as calculated for the CS method, remains zero for a wide range of S/D removal cases (0 – 80%), indicating an exact reconstruction of the center of the simulated spherical inhomogeneity. For the same range of removed measurements, the LE for SVD and SVD-reg, was calculated to be ∼ 0.6 cm. However, it should be noted, that if the z-axis shifting of the center of the inhomogeneity is ignored, then all three methods result in similar localization errors.
Fig. 7 (a) Normalized root mean square error (nRM SE) and (b) localization error (LE) for several percent values of removed S/D pairs. The absorption coefficient of the inhomogeneity was assumed μa,I = 0.08 cm−1.
5. Discussion
An implementation of the compressed sensing (CS) theory was presented for diffuse optical tomography (DOT) reconstruction and was compared to a “traditional” linear image reconstruction algorithm. In “traditional” approaches, in order to minimize the image artifacts due to undersampling, a large number of measurements is required. CS method, on the other hand, randomly samples the S/D plane (see Fig. 2) and by applying an appropriate sparsifying basis, we have shown that the CS method can provide more accurate 3-D images of the relative absorption coefficient from a significantly reduced number of measurements. The quality of the reconstructed images was compared to that of SVD reconstructions using several image quality metrics. Both randomly and regularly under-sampled SVD were investigated.
The proposed CS method was shown to exhibit superior performance than both SVD and SVD-reg, according to the calculated linear signal strength, observed contrast, contrast-to-noise ratio, normalized root mean square error and localization error. The CS reconstructions were demonstrated to be more robust than the corresponding SVD reconstructions, as the number of measurements was reduced. Furthermore, the location of the simulated spherical inhomogeneity was accurately reconstructed at z = −1.5 cm using CS, as opposed to the shifted localization towards z = −0.9 cm, that commonly occurs in SVD methods. The favorable effects of random sampling involved in the CS theory, were also demonstrated by the superior performance of SVD with respect to SVD-reg.
An ℓ1-penalized minimization process is utilized in the CS method which leads to increased computational requirements compared to of other linearized DOT algorithms. Significant research has been directed to develop efficient inverse reconstruction algorithms that could solve the ℓ1-regularization and optimization problem at a lower computational cost [38, 47, 49].
As discussed in the previous section, the SVD and SVD-reg reconstructions were performed based on an optimal regularization parameter λ, obtained by L-curve analysis. The value of λ is known to impose a tradeoff between the goodness-of-fit in the inverse problem and the associated smoothness [36]. The adopted λ in our simulations resulted in smoother image reconstructions for SVD and SVD-reg, compared to CS, which inevitably, decreased the calculated observed contrast and increased the contrast-to-noise ratio. It should be noted, that different λ values were investigated, but did not alter significantly the qualitative comparison between the three techniques under consideration. On the other side, the corresponding optimal regularization parameter Λ involved in the CS inverse algorithm, was selected based only on a heuristic process. An automated and optimized approach for the selection of Λ is required, so that the full capabilities of the proposed method can be identified (see Section 2.3).
In the described CS implementation, the discrete Fourier transform (DFT) was applied as a simple, yet appropriate sparsifying basis function. However, the selection of the optimal basis to sparsify the associated signal is an important factor in the theory of CS [48, 50]. In cases where the optical image cannot be efficiently approximated by a sparse representation (e.g. in optical mammography [2]), more complex basis functions might be required. Studying the optimal sparse expansion of the investigated signal will be the topic of future research.
Finally, as discussed earlier in this paper, measurement noise was not considered in our simulations. In practice, measurement noise would be present, thereby degrading the performance metrics. Therefore, the results of our analysis can be regarded as an upper bound of what can be achieved in absolute terms. Note, that the performance of CS in case of simulated noise was studied for dynamic MRI [51] and it was shown that CS continues to perform better despite descreased signal-to-noise ratio and reduction of number of measurements. In DOT experiments, noise presents itself in many different forms and there are other confounding factors such as the modeling of the boundaries. This would require re-optimization of all three algorithms that we have compared in this paper. There are different approaches to CS that could be utilized to deal with this point which beyond the goals of this paper. The robustness of the CS method based on real, experimental measurement data will be investigated in a subsequent paper. Furthermore, application of the proposed CS framework to more complex, nonlinear DOT models will be examined.
The results of our studies presented in the previous section, indicate a promising perspective for the future application of the CS theory to the DOT problem. Further studies are required to identify the full capabilities of this methodology within the DOT context.
Acknowledgments
This work was partially funded by Fundaciò Cellex Barcelona, Marie Curie IRG ( FP7-PEOPLE-2009-RG: 249223 RPTAMON) and Instituto de Salud Carlos III (FIS). The authors wish to acknowledge the useful discussions with Udo Weigel and Peyman Zirak of the Institute of Photonic Sciences (ICFO), with Saurav Pathak and Regine Choe of the University of Pennsylvania and with Cornelius Weber of the University of Hamburg.
References and links
1. D. R. Leff, O. J. Warren, L. C. Enfield, A. Gibson, T. Athanasiou, D. K. Patten, J. Hebden, G. Z. Yang, and A. Darzi, “Diffuse optical imaging of the healthy and diseased breast: a systematic review.” Breast Cancer Res Treat 108, 9–22 (2008). [CrossRef]
2. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Reports on Progress in Physics 73, 076701 (2010). [CrossRef]
3. S. R. Arridge, “Optical tomography in medical imaging,” Inverse problems 15, R41–R93 (1999). [CrossRef]
4. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol 50, 1–43 (2005). [CrossRef]
5. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Problems 25, 123010 (2009). [CrossRef]
6. R. J. Gaudette, D. H. Brooks, C. A. DiMarzio, M. E. Kilmer, E. L. Miller, T. Gaudette, and D. A. Boas, “A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient,” Phys Med Biol 45, 1051–70 (2000). [CrossRef] [PubMed]
7. S. R. Arridge and M. Schweiger, “A gradient-based optimisation scheme for optical tomography,” Opt. Express 2, 213–226 (1998). [CrossRef] [PubMed]
8. J. P. Culver, V. Ntziachristos, M. J. Holboke, and A. G. Yodh, “Optimization of optode arrangements for diffuse optical tomography: A singular-value analysis,” Optics Letters 26, 701–703 (2001). [CrossRef]
9. J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: Evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Medical Physics 30, 235 (2003). [CrossRef] [PubMed]
10. B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Applied optics 38 (1999). [CrossRef]
11. E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Communications on Pure and Applied Mathematics 59, 1207 (2006). [CrossRef]
12. D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory52, 1289–1306 (2006). [CrossRef]
13. C. E. Shannon, “Communication in the presence of noise,” Proceedings of the IRE 37, 10–21 (1949). [CrossRef]
14. H. Nyquist, “Certain topics in telegraph transmission theory,” Transactions of the American Institute of Electrical Engineers p. 617 (1928). [CrossRef]
15. D. L. Donoho, “For most large underdetermined systems of linear equations the minimal l1-norm solution is also the sparsest solution,” Communications on Pure and Applied Mathematics 59, 797–829 (2006). [CrossRef]
16. D. Needell and R. Vershynin, “Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit,” Foundations of computational mathematics 9, 317–334 (2009). [CrossRef]
17. R. Baraniuk, “Compressive sensing,” Lecture notes in IEEE Signal Processing magazine 24, 118–120 (2007). [CrossRef]
18. J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Transactions on Information Theory 53, 4655 (2007). [CrossRef]
19. D. Needell and J. A. Tropp, “Cosamp: Iterative signal recovery from incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis 26, 301–321 (2009). [CrossRef]
20. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse mri: The application of compressed sensing for rapid mr imaging,” Magnetic Resonance in Medicine 58, 1182–1195 (2007). [CrossRef] [PubMed]
21. H. Yu and G. Wang, “Compressed sensing based interior tomography,” Physics in Medicine and Biology 54, 2791–2805 (2009). [CrossRef] [PubMed]
22. J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography.” IEEE transactions on medical imaging 28, 585–594 (2008). [CrossRef]
23. Z. Guo, C. Li, L. Song, and L. V. Wang, “Compressed sensing in photoacoustic tomography in vivo,” Journal of Biomedical Optics 15, 021311 (2010). [CrossRef] [PubMed]
24. D. Liang, H. F. Zhang, and L. Ying, “Compressed-sensing photoacoustic imaging based on random optical illumination,” International Journal of Functional Informatics and Personalised Medicine 2, 394–406 (2009). [CrossRef]
25. G. H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (piccs): a method to accurately reconstruct dynamic ct images from highly undersampled projection data sets,” Medical physics 35, 660 (2008). [CrossRef] [PubMed]
26. Z. Xu and Y. L. Edmund, “Image reconstruction using spectroscopic and hyperspectral information for compressive terahertz imaging,” J. Opt. Soc. Am. A 27, 1638–1646 (2010). [CrossRef]
27. D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17, 13040–13049 (2009). [CrossRef] [PubMed]
28. J. C. Ye, S. Y. Lee, and Y. Bresler, “Exact reconstruction formula for diffuse optical tomography using simultaneous sparse representation,” in “Biomedical Imaging: From Nano to Macro, 2008. ISBI 2008. 5th IEEE International Symposium on,” (2008), pp. 1621–1624.
29. N. Cao, A. Nehorai, and M. Jacob, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Optics Express 15, 13695–13707 (2007). [CrossRef] [PubMed]
30. P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Applied Optics 46, 1679–1685 (2007). [CrossRef] [PubMed]
31. H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 1: l1 regularization,” Opt. Express 18, 1854–1871 (2010). [CrossRef] [PubMed]
32. M. Süzen, A. Giannoula, P. Zirak, N. Oliverio, U. M. Weigel, P. Farzam, and T. Durduran, “Sparse image reconstruction in diffuse optical tomography: An application of compressed sensing,” in “OSA Biomedical Topicals,” (Miami, FL, USA, 2010).
33. A. C. Kak and M. Slaney, “Principles of computerized tomographic imaging,” New York (1999).
34. M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. thesis, University of Pennsylvania (1996).
35. G. H. Golub and C. Reinsch, “Singular value decomposition and least squares solutions,” Numerische Mathematik 14, 403–420 (1970). [CrossRef]
36. P. C. Hansen, “Analysis of discrete ill-posed problems by means of the l-curve,” SIAM review 34, 561–580 (1992). [CrossRef]
37. E. Candès and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Problems 23, 969–985 (2007). [CrossRef]
38. S. Kim, K. Koh, M. Lustig, S. Boyd, and D. Gorinevsky, “An interior-point method for large-scale l1-regularized least squares. selected topics in signal processing,” IEEE Journal of 1, 606–617 (2007).
39. J. A. Tropp, “Greed is good: Algorithmic results for sparse approximation,” IEEE Transactions on Information Theory50, 2231–2242 (2004). [CrossRef]
40. M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing mri,” IEEE Signal Processing Magazine 25, 72–82 (2008). [CrossRef]
41. D. A. Boas, M. A. O’Leary, B. Chance, and A. G. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications,” Proceedings of the National Academy of Sciences of the United States of America 91, 4887 (1994). [CrossRef] [PubMed]
42. H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Applied optics 48, 137–143 (2009). [CrossRef]
43. H. Ponnekanti, J. Ophir, and Y. Huang, “Fundamental mechanical limitations on the visualization of elasticity contrast in elastography,” Ultrasound in Medicine and Biology 21, 553–543 (1995). [CrossRef]
44. X. M. Song, B. W. Pogue, S. D. Jiang, M. M. Doyley, H. Dehghani, T. D. Tosteson, and K. D. Paulsen, “Automated region detection based on the contrast-to-noise ratio in near-infrared tomography,” Applied optics 43, 1053–1062 (2004). [CrossRef] [PubMed]
45. R. D. C. Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria (2010). ISBN 3-900051-07-0.
46. M. Süzen and T. Durduran, “Basic dot,” GNU R Simulation Software for Diffuse Optical Tomography and Compressive Sampling (2009,2010).
47. E. van den Berg and M. P. Friedlander, “Probing the pareto frontier for basis pursuit solutions,” SIAM Journal on Scientific Computing 31, 890–912 (2008). [CrossRef]
48. R. A. DeVore, “Deterministic constructions of compressed sensing matrices,” Journal of Complexity 23, 918–925 (2007). [CrossRef]
49. E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted l1 minimization,” Journal of Fourier Analysis and Applications 14, 877–905 (2008). [CrossRef]
50. J. M. Duarte-Carvajalino and G. Sapiro, “Learning to sense sparse signals: Simultaneous sensing matrix and sparsifying dictionary optimization,” IEEE Trans. Image Processing18, 1395–408 (2009). [CrossRef]
51. U. Gamper, P. Boesiger, and S. Kozerke, “Compressed sensing in dynamic mri,” Magnetic Resonance in Medicine 59, 365–373 (2008). [CrossRef] [PubMed]
