1. Introduction

Diffuse optical tomography (DOT) is a non-invasive medical imaging technique, which uses low power light in the near infrared (NIR) range to detect hemodynamic changes inside the brain. In particular, the absorption spectra of oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (Hb) in the NIR range differ distinctly, resulting in the detectability of both forms of hemoglobin in the targeted cortical regions using two or more wavelengths. In general, DOT provides a noninvasive, portable, and cost-effective imaging solution with a higher temporal resolution. However, the spatial resolution and depth penetration of DOT are limited compared to functional magnetic resonance imaging (fMRI). Recently, DOT has gained considerable attention and become an important tool for functional brain imaging [1–4] to be used in various biomedical applications [5,6].

Similar to many other medical imaging modalities, reconstructing a high quality DOT image from a limited number of measurements is a great challenge. A variety of solutions have been proposed and explored to improve the spatial resolution of DOT by (1) optimally selecting optode arrangements [7–9], (2) applying different regularization terms in image reconstruction algorithms [9–11], and (3) reducing the measurement noise after noise information is calculated and used to weight/balance the measurements from different channels [12,13]. However, acquisition of the noise information usually requires a large amount of measurement samples, which might be difficult to collect in many DOT applications.

In this paper, we wish to introduce a novel noise quantification method, which can be utilized with noise normalization and sparsity regularization to greatly improve DOT spatial resolution. This will be demonstrated in following sections by both computer simulation and laboratory phantom results. Specifically, we have derived a noise quantification expression by associating the Rytov linear approximation of the forward diffusion model with noise analysis. The relative noise variance related to the Rytov approximation can be computed via an equation consisting of the magnitude of original measurements and the variance of measurement noise. This novel noise quantification method is practical to most of DOT applications, where the Rytov approximation can be applied. Furthermore, this noise quantification method is implemented with noise normalization and sparsity regularization-based reconstruction to form an effective method for DOT image reconstruction. The Noise normalization reduces distortion of reconstructed objects and the sparsity regularization based reconstruction improves spatial resolution of the image. Hence when combining the two, the quality of the reconstructed images is improved greatly compared to that of using the conventional method.

The rest of the paper is organized as follows: In section 2, we will review image reconstruction algorithms for DOT. In section 3, we derive the novel method to quantify noise variance that can be used to enhance DOT image reconstruction. In Section 4, we will present both computational and experimental DOT results to illustrate the effectiveness of our new method, followed by final discussion and conclusion.

2. Review of image reconstruction algorithms for DOT

A medical imaging problem can be generally divided into two sub-problems: the forward problem, which deals with the physical model of the imaging system, and the inverse problem, which converts the measurements into an image to show hidden objects or underlying information [4].

When light propagates in biological media, scattering and absorption are two fundamental interactions between light and tissue. Given that the reduced scattering coefficient is much greater than the absorption coefficient for near infrared light, this process can be modeled as the diffusion approximation (DA) of the Boltzmann transport equation [14]. The DA of the frequency domain system is given by:

In functional brain imaging, the changes of light scattering induced by neural activation are relatively small and therefore can be ignored without causing significant errors [3,15]. In the field of diffuse optical tomography or near infrared spectroscopy, two common forms are used to solve inverse problem [3]. The Born approach writes U(r) = U₀(r) + U_sc(r), while the Rytov approach writes U(r) = U₀(r) exp[U_sc(r)], where U₀(r) and U(r) are the spatial part of the photon fluence rate in the homogeneous background and in homogeneous medium, while U_sc(r) is the correction to U₀ as a result of the heterogeneities. The linearized Rytov approach, i.e., U(r) = U₀(r) × exp[U_sc(r)], leads to a matrix equation, which can be written as U_sc(rd, rs) = ln[U(rd, rs)/U₀(rd, rs)]. Since most of functional brain imaging experiments measure relative changes of optical signals with respect to a baseline reference, the measured signals can be well expressed by or matched with U(r_d, r_s) and U₀(r_d, r_s) for further data/image processing, where r_s and r_d represent a light source and detector position. Therefore, Rytov linear approximation is a better fit for the DOT brain image studies. The solution of Eq. (1) can be expressed by the linear Rytov approximation [16].

On the other hand, the modified Beer-Lambert law [17] can be written as I = I_s10^-µaL and thus OD = log(I_s/I) = μ_aL, where μ_a is the absorption coefficient along the optical pathlength, L, from a light source at position r_s to a light detector at position r_d, I_s and I are the incident and detected light intensity. If a relative absorption change, Δμ_a, occurs within the medium along the optical pathlength from the light source to the detector between two measurement states, Φ₀ and Φ, the Rytov approximation is related to the modified Beer-Lambert law by the following equation [15]:

Equation (2) can be further discretized into the following linear equation:

2.1 Regularization algorithms used for DOT

Due to the limited numbers of sources and detectors used in DOT, the number of measurements obtained is much fewer than the number of pixels to be reconstructed. As a result, Eq. (3) is highly ill-posed; its solution is not unique [4,19]. Therefore, regularization to its solutions is often required. The conventional way to regularize the solutions is to use ||x||₂, which is the l₂ norm of x, as defined by ||x||₂ = $\sqrt{\Sigma i{\left|x_i\right|}^2}$. This approach is also known as the Tikhonov regularization method [20], where the image can be reconstructed from the following equation:

However, using the l₂ norm as the regularization term for DOT image reconstruction over-penalizes the pixels that have near-zero coefficients, producing many pixels with unwanted small coefficients. This approach tends to blur reconstructed images and results in poor spatial resolution. As a result, a reconstructed object in DOT often appears to have low spatial resolution and unclear boundaries compared to the actual imaged object(s) and to those obtained by other medical imaging modalities (e.g., CT or MRI). Hence, one of the main challenges in DOT is to improve the spatial resolution of reconstructed images.

In general, many images in DOT problems are usually sparse, namely, most of the coefficients in matrix x are zero or close to zero. For example, in the application of functional brain imaging [15], a reconstructed DOT image reflects relative changes in light absorption within imaged areas between two different brain functional states. For better visualization, the measured region interrogated by the optical probe array is often chosen to be reasonably larger than the activation area. Thus, a large portion within the reconstructed image will have little changes in absorption, or to be close to zero with respect to the activation area. Given this condition, researchers have utilized sparsity regularization or l₁ norm regularization to reconstruct DOT images for improving its respective spatial resolution [9,11,21,22]. Let us briefly inspect the specificity and/or suitability of each approach for DOT image reconstruction.

Based on the brief review given above, it is clear that the l₁ norm of image coefficients serves more effectively as the regularization term for sparse images than the l₂ norm [24]. Hence, in our development, we first utilized the sparsity regularization (SR) for DOT by replacing the l₂norm of x in Eq. (4) with an l₁ norm term. The image reconstruction task now becomes to solve the following equation [25]:

2.2 Measurement noise removal for DOT

Measurement noise always exists, shields the true response signals, and leads to wrong conclusions when DOT is utilized to image brain activities stimulated by given tasks. In current practice, instrumental noise is usually inspected after data acquisition of functional DOT. The channels showing significant signal oscillations or with a low signal to noise ratio (SNR) are removed prior to removal of physiological noises and image reconstruction process [30,31], where SNR is calculated as 10 × [log_{10 $(\frac{signal power}{noise variance})$}. However, the noise level within the remaining channels may still be substantial and vary among them. Methods to design a weighting or balancing matrix for the Tikhonov method were explored and discussed [12,13]. In ref [13], the so called generalized least-squares (GLS) method modified the Tikhonov method in Eq. (4) as follows:

3. Novel development on quantification of noise variance

In practice, the measurement noise is always unavoidable. As in many of previous reports, the noise within b given in Eq. (3) is usually treated as white Gaussian noise (WGN) [9,11,32]. However, this might not be true in actual experimental measurements when Rytov linear approximation is applied. This is because WGN should be added to or included in the raw measurement data, not to the logarithm of a ratio of two raw measurements. For example, Fig. 1 shows noise variances seen at different channels with variable S-D separations of a laboratory phantom experiment [7], where 100 data samples were collected from all channels with and without the object within the liquid phantom. Then the natural logarithm between the two sets of 100 samples at two measurement states were calculated to obtain 100 relative measurements b [e.g., b(1) = ln(Φ₀(1)/Φ(1)), b(2) = ln(Φ₀(2)/Φ(2)),…]. Then the sample variance of these 100 data points of each channel were calculated and shown in Fig. 1. From this figure, it clearly shows the variation of noise variances among different S-D separations; namely, a larger variance is obtained as the S-D separation increases, particularly when the separation is larger than 4 cm. Since the noise variance is not identical among different channels, the relative noise is no longer WGN.

 

Fig. 1 Noise variance of different S-D pairs from a phantom experiment.

Download Full Size | PDF

To better analyze noise variances of DOT measurements, in this section, we first derive a quantitative expression for noise variance, and then utilize it as a weight-matrix to normalize measurement signals, followed by completion of noise-suppression, sparsity-regularized DOT image reconstruction. Our novel development reported in this paper lies on two specific aspects of DOT images: (i) varying variances of noise and (ii) sparse reconstructed objects. In brief, our approach integrates (1) quantification and (2) normalization of noise variance with (3) sparsity regularization, as noted QNNVSR.

3.1 Noise variances of DOT signals

Following the concept given in Eq. (8), we modify Eq. (7) by including a weight matrix D, which normalizes noise variances of the fidelity term [13,33], as given below:

The noise covariance matrix might be easy to obtain from both simulations and phantom experiments with given objects to be imaged. In computer simulations, the statistical information is known in advance. Even in a laboratory phantom experiment, each measurement of individual S-D pairs can be acquired with relatively stable values; hence, the noise variance of each S-D pair or channel can be found with conventional estimation method from a number of measurement samples. However, such a covariance matrix might be difficult to obtain in biological applications where the measurements are dynamic or do not have enough acquisition samples available for estimation of the covariance information, particularly for in vivo imaging of brain activities and functions. Hence an alternative approach to quantify the noise variance is needed.

In functional brain imaging by DOT, researchers often measure or determine changes in light intensity or photon density at different time instances. It follows that with the linear Rytov approximation, the noise-free, relative signal changes measured by DOT [as seen in Eq. (3)] can be expressed by

Mathematically, this equation implies that the relative noise vector, w, is affected not only by noises at both time instants, w₁ and w₂, but also by actual signals from imaged objects or activities at both temporal states, Φ₀ and Φ. Consequently, the assumption that the noise variances (i.e., elements of w) of different channels at different S-D separations have a distribution of WGN [9,11,32] is not accurate or appropriate in this case.

As shown in ref [34], incorporating the noise variance information into image reconstruction is a meaningful step to achieve high-quality DOT. Before starting to derive a novel noise quantification method in the next sub-section, we state the following assumptions that simplify the noise quantification process when the covariance matrix of Eq. (12) is approximated:

The first assumption is acceptable if the noise or signal fluctuation of each light source is much smaller than the respective light intensity and each detector takes independent measurement from one another at independent locations. The second assumption is also valid since the amplitude of instrument noise is usually less than that of the measured signal; for example, the system noise is usually assumed to be a few percent of the measured signals [5,36].

3.2 Quantification of noise variance for DOT measurements

Now, based on Eq. (12) and assumption (A1), the noise added to the i-th S-D channel is

Next, let a = w₁(i)/ϕ₀(i) and b = w₂(i)/ϕ(i), where w₁(i) and w₂(i) are random variables with a zero-mean Gaussian distribution and variances of σ₁² and σ₂², respectively. Then, because ϕ₀(i) and ϕ(i) are to be constants, a and b are also zero-mean Gaussian random variables with variances of σ_a² = σ₁² /ϕ₀²(i) and σ_b² = σ₂² /ϕ²(i), respectively.

Therefore, Eq. (13) becomes

Based on Assumption (A2), both a and b are Gaussian random variables with small variance much less than 1. Next, we let

Applying Taylor expansion to Eq. (14) yields:

Following the same steps, we have

Since the variance of x is much less than one, Eq. (15) can be approximated as:

Hence, the variance of the relative noise of one S-D channel, i, can be quantified as shown in Eq. (17), consisting of original noise variances, σ₁² and σ₂², and the amplitude of raw measurements, ϕ₀(i) and ϕ(i), at two different time instances. Note that even if the noise terms, w₁ and w₂, are WGN, due to the ϕ₀(i) and ϕ(i) terms, the noise variances of Eq. (12) are not necessarily the same at different channels, as derived and expressed in Eq. (17).

3.3 Normalization of noise variance with sparsity regularization for DOT reconstruction

Note that in Eq. (12), the measurement noise w₁ and w₂ are the noise from the same S-D channel at two different time instances. If further assuming that the noise does not change largely over time, namely, σ₁² = σ₂² = σ_w², Eq. (17) can be modified as

It is noteworthy that estimated variance of relative noise at each channel can be considered as its noise power to signal power ratio of the respective channel since σ_a² = σ₁² /ϕ₀²(i) and σ_b² = σ₂² /ϕ²(i). Hence, its reciprocal or the weighting coefficient given in matrix D is closely associated with or equivalent to the SNR of each channel. As the weighting matrix D is applied in Eq. (9), measurement channels with a higher SNR will gain higher weights compared to those with a low SNR. Equivalently, matrix D will suppress the signals from the measurement channels with a lower SNR, thus reducing the impact of noise. It should also be noted that for the applications where the original noise level cannot be estimated, one can approximately assume that noise variances among different channels are nearly identical (i.e., σ₁ = σ₂ = a small constant) for all channels. Then, measured amplitudes can be used to estimate relative noise variances and to build the weighting matrix.

In Section 3 above, we derived and explained our newly developed algorithm on quantification and normalization of noise variance with sparsity regularization, QNNVSR, for DOT image reconstruction. Specifically, in order to perform QNNVSR, the following steps need to be followed:

4. Algorithm validation by numerical simulations and laboratory experiments

4.1 Procedures and setup for computer simulations

In the computer simulation experiment, a CW DOT system was considered. The geometry of sources and detectors is illustrated in Fig. 2, where 25 bifurcated optodes were placed as a 5 × 5 square grid over the top surface of a semi-infinite medium and centered at the origin. The distance between every two closest optodes was 1.4 cm. By using the first to sixth nearest S-D pairs, 188 channels could be obtained, based on the laboratory phantom measurements (see Section 4.3). White Gaussian noise with a variance of 10⁻⁴, which is approximately equal to the noise variance found in phantom experiments (see Fig. 1), was added to both background measurements, Φ₀, and inhomogeneous measurements, Φ.

 

Fig. 2 Probe array setup for both simulation and phantom experiments. The unit in both x and y directions is cm.

Download Full Size | PDF

The absorption and reduced scattering coefficients of the homogeneous background medium were chosen in a reasonable range [37] as μ_a = 0.06 cm⁻¹ and μ_s’ = 8.2 cm⁻¹, respectively. The forward model was computed analytically as in refs [18,38]. Since the main focus of this study is to improve the spatial resolution for functional brain imaging, only relative changes in absorption within activated brain areas are considered. Thus, the difference in absorption, Δμ_a, between the imaged brain area and background is what we need to reconstruct for DOT images. An “L” shape absorber with Δμ_a = 0.1 cm⁻¹ was placed in the center of optode array as shown in Fig. 3, at a depth of 1.5 cm.

 

Fig. 3 Original absorbing object used for simulation experiment; the x and y dimension is in cm. Color bar represents Δμ_a in cm⁻¹.

Download Full Size | PDF

In order to demonstrate the advantage of the QNNVSR method, we compared the reconstruction results with the conventional Tikhonov regularization method as given in Eq. (4), with SR in Eq. (7), and with the GLS method in Eq. (8). The regularization parameters in Tikhonov regularization methods were selected with the L-curve approach [19], as it is effective to provide the optimal solution and thus commonly used with Tikhonov regularization methods. For all the other methods, the regularization parameters were selected from a number of trials that yielded the least distortion of the imaged object's shape. Moreover, methods besides Tikhonov regularization were all performed by the SALSA algorithm [28], which is a solver for convex optimization problems with the same stopping criteria.

Two quantitative metrics were used to evaluate the reconstruction results: the structural similarity index metric (SSIM) [39] and the contrast-to-noise ratio (CNR) [40].

The SSIM is given by

On the other hand, by defining the area of targeted object as region of interest (ROI) and the rest area as background (BKG), the CNR is calculated by:

4.2 Results from computer simulations

4.2.1 Validation of noise quantification

In order to verify the noise quantification method just derived in Section 3.2, two sets of data from two spherical objects were simulated and collected: Φ₀ readings were related to the base line measurements from all S-D channels in the absence of the absorbers. Corresponding Φ values were acquired from the measurements in the presence of the objects inside the medium.

Both Φ₀ and Φ were simulated 100 sampling times, separately. The mean and variance of Φ₀ and Φ over 100 samples were computed and used to calculate the relative noise variance using Eq. (17). These calculated noise variances were compared with the actual ones directly estimated from 100 simulated relative measurement samples, i.e. b, as shown in Fig. 4. It shows clearly that the approximated noise variances quantified by Eq. (17) are consistent with the actual ones. The mean square error between the estimated and actual noise variance was calculated and equal to $1.76 × 10⁻⁹. Note that since the noise term here corresponds to relative changes in light absorption in a log scale [see Eq. (2)], it is derived from a fraction of two raw measurements. Thus, it does not have a unit.

 

Fig. 4 Comparison of noise variance from all simulated channels. Calculated noise variances by Eq. (17) are represented by y-axis and the actual noise variances computed from simulated 100 measurements are shown in x-axis.

Download Full Size | PDF

4.2.2 Comparison of reconstructed images using different reconstruction algorithms

Several reconstructed images from computational experiments using four different reconstruction algorithms are shown in Figs. 5(a) to 5(d). The result by Tikhonov regularization is shown in Fig. 5(a) with blurry and diffuse image boundary; the shape of constructed image is completely distorted from the “L” shape [see Fig. (3)] into an elliptical object. After performing normalization of noise variance with Eq. (17), the GLS method improves the shape of the reconstructed object [see Fig. 5(b)] compared to the Tikhonov method. However, the recovered object is still fuzzy and has much background noise. It is clear that the SR approach shows a less blurred image as shown in Fig. 5(c), but the reconstructed object tends to have a “C” shape, still somewhat deviating from the “L” shape. Finally, the reconstruction with the QNNVSR method gives rise to a much improved reconstructed image that recovers the original shape almost completely as shown in Fig. 5(d). Note that since each reconstruction algorithm has added different weighting factors or matrix, the absolute values of absorption changes reconstructed by different algorithms could not be compared consistently. Thus, all the reconstructed images are normalized to their respective maximum for fair comparison.

 

Fig. 5 Comparison of normalized reconstruction results from the simulation experiment with: (a) Tikhonov method, (b) GLS, (c) SR, and (d) QNNVSR.

Download Full Size | PDF

The SSIM [see Eq. (20)] and CNR [Eq. (21)] values of the reconstructed object by all four methods were computed and listed in Table 1. Notice that QNNVSR gives rise to the highest SSIM value among all the methods tested, meaning that quantification and normalization of noise variance with sparsity regularization can best enhance image reconstruction in DOT, because it is able to recover the object shape very close to the actual one. Moreover, a high CNR indicates that most of non-zero pixels in the reconstructed image are consistently located within the expected ROI. Even though the SR method yielded a clean background, its CNR is the lowest. This is partially due to the distortion of the reconstructed object's shape; some of the non-zero pixels are fallen outside the expected ROI, so they are considered as background noise in the CNR calculation. In sharp contrast, our QNNVSR method recovers the object in both location and shape, almost perfectly having a clean background. Therefore, its CNR is the highest among all four methods.

Table 1. SSIM and CNR values of reconstructed images from numerical experiments with Tikhonov, GLS, SR, and NQNSR

View Table | View all tables in this article

4.3 Setup for laboratory experiments

In addition to computational experiments, laboratory phantom experiments were also performed to test our newly derived reconstruction method with real experimental noise.

In the phantom experiment, the experimental setup and procedures were similar to those used in [7,38]: A CW-based DOT imaging system (DYNOT, NIRx, New York) was used to obtain the measurements from a container of 15 × 10 × 10 cm³ filled with 1% Intralipid solution. The optodes were placed on the very top surface of the intralipid solution such that the tips of optodes were just touching the liquid phantom surface. In this way, no air gap existed between the tips and liquid to minimize the refractive index (RI) mismatch. This setup provided us with a very similar boundary condition to that in a light-tissue interaction situation. The wavelength of incident light was 830 nm and the absorption and reduced scattering coefficients of the homogeneous background medium were measured (Model 96208, ISS Inc., Champaign, Illinois) as μ_a = 0.12 cm⁻¹ and μ_s’ = 8.8cm⁻¹, respectively. The geometry of sources and detectors was the same as that used for the computer simulation as illustrated in Fig. 2, where 25 bifurcated optodes were placed as a 5 × 5 square grid over the top surface of the Intralipid medium and centered at the origin. By considering all the possible combination of different optodes, 300 measurement channels could be obtained and the SNRs of all these channels were calculated as shown in Fig. 6.

 

Fig. 6 Calculated SNRs among different channels in the phantom experiment (“+” means outliers). X-axis represents S-D separations, while discrete boxes from Y axis mark the 1st, 2nd, 3rd, … and nth nearest S-D pairs.

Download Full Size | PDF

It is seen that SNRs of the first (1.4-cm separation) to the sixth (4.2-cm separation) nearest S-D pairs have a median value larger than 45 dB, relatively higher than the rest of the separations. Therefore, to avoid channels with a significantly low SNR level, only first to sixth nearest S-D pairs were selected for data processing (i.e., channels with S-D separations equal to or less than 4.2 cm). As a result, 188 channels were used for image analysis and reconstruction. Two sets of objects were used separately: (1) Two spherical absorbers with identical sizes (diameter = 1 cm) were placed along the y axis (see Fig. 2), 1.5 cm below the measurement surface. The two absorbers used were two plastic black beads, which therefore should be treated as “pure” absorbers. The center-to-center distance between the two absorbers was 1.5 cm. (2) A black “L” shaped object, as shown in Fig. 7, was placed along the y axis (see Fig. 2) and at 1.5 cm below the measurement surface. The total grid area on the background of Fig. 7 is 25 mm × 25 mm. The L-shaped absorber was made from a vinyl eraser. Similar to the first phantom experiment, the absorber was also black and thus should be treated as a “pure” absorber.

 

Fig. 7 Absorber used in the second phantom experiment, where the size of the background grid has a total area of 25 mm × 25 mm, with a grid unit of 5 mm in each direction.

Download Full Size | PDF

Moreover, strong light scattering of biological tissue causes the detection sensitivity of DOT to attenuate exponentially with increased depth, resulting in depth localization in DOT to be poor in general. Thus, it becomes unrealistic to accurately quantify the local absorption perturbation without knowing its actual depth. A variety of efforts have been made by different research groups to improve the accuracy of in-depth localization. One common method that has been used in functional brain imaging was to combine DOT with prior anatomical information from MRI. It is well known that most of optical signals measured by functional DOT come from cortical regions of the human brain below the measurement surface [15]. Therefore, a mathematical constraint to force the DOT images to be reconstructed on the cortical layers is added as a spatial prior based on structural and functional MRI [41] to reduce DOT partial volume error and to improve quantitative accuracy of reconstructed DOT images. This means that in DOT reconstruction, a 2D image is reconstructed at a pre-selected depth (between 1 and 1.5 cm). In this way, a 3D reconstruction problem is reduced to a 2D slice reconstruction problem. This is how we obtain 2D reconstructed images from our 3D experimental space. Since the thickness of extra-cranial or extra-cerebral layers is relatively constant across adult human heads (between 1 and 1.5 cm), such a cortical depth constraint can be applied in practical and/or clinical applications without causing too large errors. This depth constraint needs to be adjusted for pediatric and adolescent populations based on corresponding MRI images.

Figure 8 shows reconstructed images of the first phantom experiment with two spherical objects by using (a) Tikhonov regularization, as given by Eq. (4), where the regularization parameter was selected by the L-curve method [18,19] with λ = 10⁻⁴, (b) GLS with λ = 3 × 10⁻⁶, (c) SR, as given in Eq. (5), λ = 3 × 10⁻⁴, and (d) QNNVSR with λ = 3 × 10⁻⁵. Despite slight distortion of one absorber, SR successfully reconstructs two objects along the vertical (i.e., y) axis, while Tikhonov and GLS methods fail to recover two distinct absorbers inside the medium. This demonstrates that sparsity regularized reconstruction indeed improves the spatial resolution of DOT images. However, due to the measurement noise, one of the two objects in the image tends to be smaller than its actual size. When QNNVSR is applied, the two imaged objects are retrieved with a comparable size to the actual one with minor distortion in shape as seen in Fig. 8(d).

 

Fig. 8 Normalized reconstructed images of the first phantom experiment (two spherical objects) obtained by using (a) Tikhonov regularization, (b) GLS, (c) SR, and (d) QNNVSR.

Download Full Size | PDF

In order to emphasize and compare the spatial contrast, cross sections of non-normalized image amplitudes recovered with Tikhonov, GLS, SR, and QNNVSR at x = 0 (i.e., along y-axis of Figs. 2 or 8) of the image plane are compiled and plotted in Fig. 9. As mentioned in Section 4.2.2, since each reconstruction algorithm has added different weighting factors or matrix, the absolute values of absorption changes reconstructed by different algorithms could not be compared consistently. The major inspection here is to see which algorithm can well separate the two imaged objects. Figure 9 illustrates clearly that two targeted objects are well separated with SR and QNNVSR, while Tikhonov and GLS could not resolve two objects at all with much blurry cut-off edges. In this one-dimensional plot along this specific cross section, we cannot observe much clear distinction between the reconstruction quality by SR and QNNVSR.

 

Fig. 9 Recovered amplitudes of absorption changes along y-axis of image plane (see Fig. 2 or Fig. 8) with Tikhonov, GLS, SR, and QNNVSR reconstruction algorithm at x = 0.

Download Full Size | PDF

The reconstruction results of the second laboratory experiment with the L-shaped absorber are obtained and compared in Fig. 10. The regularization parameters used for Tikhonov, GLS, SR, and QNNVSR method are λ = 10⁻⁵, 8 × 10⁻⁷, 2 × 10⁻⁵, and 8 × 10⁻⁶, respectively. Though all four methods are able to recover an absorber in the center of the images, the quality of the reconstructed background and object vary across these four cases. Even though S-D pairs with larger separations (i.e., S-D separation >4.2 cm or after the sixth nearest S-D pairs) were removed to minimize significant noise, the reconstructed object in Fig. 10(a) by Tikhonov regularization is surrounded by many background artifacts with a distorted shape of the object. In Fig. 10(b), a blurred “L” shaped object can be roughly identified, but the background is still filled with artifacts. The result in Fig. 10(c) reconstructed with SR, on the other hand, shows a clearer background and an improved object shape. Finally, the reconstruction result by our QNNVSR method, as shown in Fig. 10(d), demonstrates that not only the shape of the reconstructed object is improved (shown to be similar to the actual one given in Fig. 7), but also the background noise is further suppressed compared to the one using SR alone.

 

Fig. 10 Normalized images of the second phantom experiment with an L-shape object (see Fig. 7) reconstructed by (a) Tikhonov, (b) GLS, (c)SR, and (d) QNNVSR. The data used for reconstruction were taken from the first to sixth nearest S-D pairs (up to 4.2 cm), having a total number of 188 channels.

Download Full Size | PDF

Furthermore, corresponding SSIM and CNR values from all four reconstructed images are calculated and listed in Table 2. The top four rows correspond to the cases with S-D separation up to 4.2 cm, namely with 188 S-D channels from the first to the sixth nearest S-D pairs. This table illustrates clearly that QNNVSR offers the best value of SSIM and CNR among all the reconstruction algorithms tested.

Table 2. SSIM and CNR values of reconstructed images from phantom experiments with Tikhonov, GLS, SR, and NQNSR

View Table | View all tables in this article

5. Discussion and future work

In this paper, we have derived and demonstrated a novel method for quantification and normalization of noise variance and then combined it with sparsity regularized reconstruction to enhance spatial resolution and image contrast of DOT images. According to our step-by-step derivation, the noise variance can be approximately estimated from relatively limited resources available, provided a short-period of time is given for baseline data acquisition [see Eqs. (17) or (18)]. Furthermore, a weight matrix D can be generated and used for normalization of noise variance [see Eq. (19)], which enhances signals with a large SNR and suppresses signals with a low SNR. Our novel work rests on the following findings: a combination between normalization of measurement data by estimated noise variance and sparsity regularization [i.e., combining Eqs. (9) and (19)] enhances the spatial resolution and the shape fidelity considerably for DOT images. Specifically, after we define and quantify the weight matrix D, we can solve the sparsity regularized DOT reconstruction using SALSA algorithm [28]. Throughout this paper, we have demonstrated our points by reconstructed images in computer simulations and laboratory phantom experiments.

This study has a few of limitations. First, we do acknowledge that a major weakness of the experimental study or setup was the use of “dark” or “pure” absorbers, which would not rigorously satisfy the mathematical requirements given by Rytov approximation, as given by Eq. (2). Such mismatch between the mathematical requirement and experimental conditions may lead to potentially inaccurate results and conclusions. To overcome this problem, further laboratory experiments using absorbers with much lower absorption coefficients are needed to confirm QNNVSR. Second, in practical DOT brain imaging applications, light absorption changes due to brain activations are usually less than 5% of the baseline signals [7,42], which is much less than the absorption change used in this laboratory experiment. Hence, in actual functional brain study, the contrast will be much lower than the phantom experiment we have utilized. It implies that we may be able to use only the 1st and 2nd nearest S-D separations for sensitive DOT measurements, depending on the actual optode setup and probe geometry. Third, since the overall development and derivation is based on CW setting with ω = 0 in Eq. (1), our method is not able to separate light scattering from light absorption, and thus it works mainly for functional brain imaging, namely, to measure only changes in HbO and Hb concentrations, assuming constant light scattering. Last, the current results were obtained using the trial-and-error approach to select the regularization parameter; we have not yet discovered a systematic method to optimally determine the regularization parameter, which will be the next step of this research project.

While our development focuses on functional brain imaging, where the images are often sparse in the spatial domain for absolution changes, the methodology can be valid and modified accordingly for broader applications. For example, if the image to be reconstructed is not sparse in the spatial domain, one can map or transform the image into another domain where the image will become sparser (e.g., for some transform S, S(x) will be much sparser than the original x). Such techniques have been utilized and demonstrated in other compressive sensing applications. Further development in this direction can be our future work that will broaden the usefulness of QNNVSR for other DOT applications.