1. Introduction

Diffuse optical tomography (DOT) is finding widespread application in the diagnosis and characterization of breast cancer [1–5], monitoring brain function [6–14], and small animal imaging for the study of disease detection, progression and treatment [15–18]. Specifically in DOT of human brain function, an imaging cap that uses an array of continuous-wave, near-infrared (NIR) light sources and detectors are placed on the scalp surface, allowing measurement of light level variance due to underlying functional changes [19]. A minimum of two sets of measurements are typically required to allow dynamic imaging, i.e. an activation set and a reference set, where the differences between these sets are used to recover spatial maps of hemodynamic changes within the brain, which relies on a forward model for light propagation in tissue. The image reconstruction problem in this case is single-step and linear, but remains ill-conditioned and ill-posed [20,21].

The two primary chromophore species of interest in cerebral hemodynamics studies are oxyhemoglobin (HbO₂) and deoxyhemoglobin (HbR). Given their difference in the near-infrared absorption spectra (Fig. 1 ), measurements at two wavelengths are typically sufficient to recover the change in concentration for both chromophores [22]. Most existing optical neuroimaging systems are equipped with two wavelengths of light sources, including the CW5 (690 and 830 nm) [14,23], the DYNOT (760 and 830 nm) [12,14,24], and the high-density-DOT imaging system at Washington University School of Medicine (750 and 850 nm) [9,10,13,25]. For this reason our analyses in this work will be focused on the dual-wavelength setup using 750 and 850 nm.

 

Fig. 1 Absorption spectra for oxyhemoglobin (HbO₂) and deoxyhemoglobin (HbR) in the near-infrared spectrum. The molar extinction coefficients used in this work are indicated by the four red-dots in this figure, which are 0.1193 mm⁻¹ and 0.3236 mm⁻¹ for HbO₂ and HbR at 750 nm, and 0.2436 mm⁻¹ and 0.1592 mm⁻¹ for HbO₂ and HbR at 850 nm.

Download Full Size | PDF

The spectrally constrained image reconstruction method was first documented by Corlu [26] in non-linear DOT, wherein absolute concentrations of chromophores are solved iteratively with a minimization function [21]. When compared with non-spectral methods, whereby images of absorption and scattering coefficients are first recovered and then un-mixed to provide chromophore concentrations, spectrally constrained reconstruction techniques are found to suppress artifacts in water and scattering power images, and also reduce crosstalk between chromophores and scatter parameters in breast imaging [27–29]. Li [30] further demonstrated that when considering dynamic (also known as temporal or differential) imaging, linear DOT algorithms that incorporate spectral priors directly in the image reconstruction, correlate the absorption coefficients across multiple wavelengths, thereby providing a solution that is more stable, albeit with an increase in crosstalk between chromophores. Such crosstalk is partially due to the suboptimal regularization and update of the Jacobian matrix (which relates a small change in measured data to changes in optical parameters) as pointed out by Eames [31], who proposed a Jacobian normalization technique in non-linear DOT that allows more uniform regularization and parameter update in image reconstruction, hence reducing imaging crosstalk. However in DOT of brain function where the imaging problem is linear and measurements are mostly reflectance rather than transmittance, such normalization techniques would inappropriately bias the measurement sensitivities towards deeper regions within the brain and cause significant imaging error.

In this work we address and resolve this crosstalk issue by proposing a novel regularization technique based on the singular value decomposition (SVD). While previous DOT literature has focused on using the SVD as a tool to understand the relationship between DOT measurements (in terms of number and location) and image quality [32] and to guide optimization of detector placement [17,33], the approach we take is to integrate the SVD directly into the image reconstruction algorithm. Previous applications in the use of SVD-based truncation methods for Jacobian inversion are well documented [34], however, our approach differs from previous work since we are not using the SVD approach to achieve a pseudo-inverse of the spectral Jacobian, but instead using it to regularize the wavelength dependant Jacobian, and then calculate the pseudo-inverse of its spectral counterpart.

Qualitative and quantitative evaluations are conducted among non-spectral method and spectral methods (using conventional and SVD-based regularization respectively) first on a two dimensional (2D) circular model and then on a three dimensional (3D) subject specific head geometry. The 2D model provides an initial proof-of-concept analysis whereas the 3D model has more degrees of freedom, is highly ill-determined and more realistic to in vivo diffuse optical tomography of human brain studies.

2. Theory

2.1 Non-spectral method

The non-spectrally constrained image reconstruction problem consists of two steps. The goal of the first step is the recovery of a temporal change of absorption $\Delta {\mu }_{\lambda }$ at wavelength λ for a selection of wavelengths. This requires two sets of boundary data ${\Phi }_{\lambda ,ref}$ and ${\Phi }_{\lambda }$ at each wavelength λ which are acquired before and after a change of absorption, $\Delta {\mu }_{\lambda }$. Given that $\Delta {\mu }_{\lambda }$ is relatively small, the forward problem can be linearized by the Rytov approximation (whereby the log of the intensity measurements are considered) [20]

where $J_{\lambda }$ is the wavelength dependent Jacobian matrix for wavelength λ, with a size of number of measurements by number of nodes. Since $J_{\lambda }$is non-invertible, the Moore-Penrose generalized inverse is applied

where ${\alpha }_{\lambda }{}^2$ is the Tikhonov regularization parameter [20,35]. Assuming a dual-wavelength system and two chromophores of interest, we have

where $M_s$ is the spectral composition matrix (containing spectral prior information) of the following form

where ${\epsilon }_{c,\lambda }$ is the molar extinction coefficient of the absorbing chromophore c at wavelength λ (Fig. 1). Hence the second step inverse problem is a spectral decomposition process

This step provides the mapping from images of change in wavelength dependent absorptions to images of change in chromophore concentrations.

2.2 Spectral method using conventional regularization

Instead of reconstructing changes in absorption at each wavelength (Eq. (2)) and then spectrally decomposing these into changes of chromophore concentrations (Eq. (5)), the spectral prior information $M_s$ can be directly incorporated into the forward problem by constructing a spectral Jacobian matrix $J_s$

By doing so the forward problem can be expressed as

and the spectral image reconstruction using conventional Tikhonov regularization becomes

where ${\alpha }_s{}^2$ is the spectral Tikhonov regularization parameter. Here it is worth noting that the weight of regularization is directly associated with image resolution [33]. Previous works have found that ${\alpha }_s{}^2$ (in Eq. (8)) must be smaller than ${\alpha }_{\lambda }{}^2$ (in Eq. (2)) in order to achieve equivalent image resolution between non-spectral and spectral methods [30]. The spectral method allows a direct mapping between changes in measurements and changes in chromophore concentrations, bypassing the absorption-based transition step in the non-spectral method. Furthermore, measurements at multiple wavelengths are utilized simultaneously and absorption coefficients are correlated across these wavelengths within the image reconstruction. However the regularization technique utilized here regularizes not only the underlying wavelength dependent Jacobian matrices $J_{\lambda }$ but also the spectral prior information $M_s$ incorporated within the spectral Jacobian matrix $J_s$, which results in numerical error in the form of crosstalk on the reconstructed images, as demonstrated later in Section 3.1.

2.3 Spectral method using SVD based regularization

Instead of constructing the spectral Jacobian matrix first and performing the regularization second, we propose an alternative regularization technique that reverses the order of these two operations. This is achieved by first applying the singular value decomposition (SVD) on each wavelength dependent Jacobian matrix, which yields a triplet of matrices

where U and V are orthonormal matrices containing the singular vectors of $J_{\lambda }$, and S is a diagonal matrix that contains the i number of singular values of $J_{\lambda }$, i.e., ${\sigma }_{\lambda ,i}$. Since applying Tikhonov regularization of weight ${\alpha }^2$ changes the singular values of $J_{\lambda }$ from ${\sigma }_{\lambda ,i}$ to $\sqrt{{\sigma }_{\lambda ,i}{}^2+{\alpha }^2}$(Fig. 2 ), we can recompose the regularized $J_{\lambda }$, denoted ${\widehat{J}}_{\lambda }$, as

 

Fig. 2 Singular value spectra of the 750nm (blue line) and 850 nm (red line) wavelength Jacobian matrices of the 2D circular model as described in section 3.1, before (solid line) and after (dashed line) Tikhonov regularization.

Download Full Size | PDF

We can then construct the regularized $J_s$, denoted ${\widehat{J}}_s$, as

The inverse problem therefore becomes

The advantage here is that regularization is applied to the Jacobian matrix while (unlike all other previously proposed methods) the spectral prior information $M_s$ is exactly preserved and therefore contains all information needed for the wavelength unmixing operation. In addition, ${\alpha }_{\lambda }{}^2$ in the non-spectral method is now directly used in the spectral method, avoiding the operation to find an equivalent ${\alpha }_s{}^2$ for comparative studies which can be problematic due to its subjective choice.

3. Methods and results

3.1 Two-dimensional circular model

To compare the image performance on crosstalk among all three methods as described in sections 2.1-2.3, a set of proof-of-concept simulations were carried out on a 2D circular model. Using a well-defined numerical model has the benefit of knowing the exact location and magnitude of the targets, allowing accurate analysis of the reconstructed images. The simulated 2D model consists of a uniform circular mesh of radius 43 mm with 1785 nodes corresponding to 3418 linear triangular elements. Sixteen co-located source/detector fibers are modeled equidistant on the external boundary and were used for continuous-wave data collection, giving rise 240 (16 × 15 where the source fiber is not used for detection) differential measurements per wavelength (Fig. 3 ). The model was assumed to be homogeneous with µ_a = 0.017 mm⁻¹, µ′_s = 0.74 mm⁻¹ at 750 nm, and µ_a = 0.019 mm⁻¹, µ′_s = 0.64 mm⁻¹ at 850 nm. A target of 5 mm in radius was located at three different depths, i.e. 13 mm, 28 mm and 43 mm, corresponding to a modeled ΔHbO₂ = +0.05mM (Figs. 5(a) –5(c)), Target). Forward data was generated using NIRFAST [36], which is a modeling and image reconstruction toolbox based on the Diffusion Approximation. For regularization, we found that in order to achieve comparable image resolution and contrast between these three methods, a lower value of α was required for the conventional spectral method (${\alpha }_s$ in Eq. (8)) as compared to the non-spectral method (${\alpha }_{\lambda }^{}$ in Eq. (2)) and the proposed SVD approach (${\alpha }_{\lambda }^{}$ in Eq. (10)), which was in line with previous finding [30]. Specifically we chose ${\alpha }_{\lambda }={10}^{-2}\times$the maximum singular value of $J_{\lambda }$ and ${\alpha }_s=5\times {10}^{-3}\times$the maximum singular value of $J_s$ in this comparative study.

 

Fig. 3 Schematic view showing the placement of 16 co-located sources (red squares) and 16 detectors (blue cross) on the boundary of a 2D circular model. Note for each source excitation, the same fiber is not used as detector, giving rise to 240 differential measurements in total.

Download Full Size | PDF

 

Fig. 5 Same as Fig. 4, but with 0.2% added noise.

Download Full Size | PDF

Linear, single-step reconstructed images of hemodynamic changes using non-spectral (denoted ‘Non-Spec’), spectral with conventional regularization (‘Conv-Spec’), and spectral with SVD-based regularization (‘Svd-Spec’) are shown in Figs. 4 –6 with 0%, 0.2% and 0.5% added Gaussian random noise respectively. Within each figure, three scenarios (a-c) having the ΔHbO₂ target moving towards the center of the circular model are presented, emulating a gradual decrease in signal-to-noise ratio (SNR) in the measurement; that is as the target moves deeper within the domain, the magnitude of the detected differential signal decreases. It is evident from Fig. 4 that in the case of 0% noise a region of crosstalk can be easily identified in the recovered ΔHbR images along with some imaging artifacts due to measurement sensitivity distribution for all three methods, with ‘Svd-Spec’ demonstrating the least magnitude of crosstalk and imaging artifacts while ‘Conv-Spec’ produces the most. As noise increases to 0.2%, imaging artifacts due to noise appear although at smaller magnitudes as compared to the crosstalk. When noise reaches 0.5% the artifacts begin to dominate the ΔHbR concentration images for both the ‘Svd-Spec’ and ‘Non-Spec’. When the target is at the center of the model (43 mm from boundary, Fig. 6(c)), the recovered images from ‘Svd-Spec’ and ‘Non-Spec’ become indistinguishable. In comparison, crosstalk and imaging artifacts in ‘Conv-Spec’ still dominate the ΔHbR concentration images even at 0.5% noise level.

 

Fig. 4 Reconstructed images of ΔHbO₂ (upper row) and ΔHbR (lower row) concentration at three different locations (a-c) using methods ‘Svd-Spec’, ‘Non-Spec’ and ‘Conv-Spec’ with 0% noise.

Download Full Size | PDF

 

Fig. 6 Same as Fig. 4, but with 0.5% added noise.

Download Full Size | PDF

For quantitative analysis, the crosstalk of ΔHbO₂ into ΔHbR, defined as the mean value of the known region of crosstalk in the ΔHbR concentration image divided by the mean value of the known target region in the ΔHbO₂ concentration image, are calculated. In Fig. 7 the plotted crosstalk versus depth of target is shown for all cases presented in Figs. 4–6. Again it can be seen that crosstalk increases along with depth of target due to decreasing magnitude of the differential signal and/or SNR. On average, crosstalk is ~$2\times {10}^{-3}$ for ‘Svd-Spec’, representing a 60% reduction from ~$5\times {10}^{-3}$ for ‘Non-Spec’, and a 98% reduction from ~10⁻¹ for ‘Conv-Spec’. At 0.5% noise the imaging artifacts dominate over the crosstalk for both ‘Svd-Spec’ and ‘Non-Spec’, particularly when the target is located at the center of the model (Fig. 6(c)), where crosstalk for both methods are approximately $5\times {10}^{-3}$, indicating little difference between the two methods, however still much better as compared to conventional spectral reconstruction, ‘Conv-Spec’. These simulations were repeated with a ΔHbR target and similar results were obtained.

 

Fig. 7 Crosstalk versus depth of target: 13, 28, 43 mm, representing scenario (a)–(c) in Figs. 4–6, respectively.

Download Full Size | PDF

3.2 Three-dimensional head model

Next we extended our analysis to a more realistic, three-dimensional finite element head model (FEM) created using the Mimics software package [37]. The 1-mm-resolution model consisted of 1,087,223 nodes corresponding to 6,289,566 linear tetrahedral elements, with each node labeled by one of the five segmented head tissue types, namely scalp, skull, cerebrospinal (CSF), gray matter and white matter, as shown in Figs. 8(a) -8(b). These structural priors were obtained from T1-weighted MPRAGE (echo time (TE) = 3.13 ms, repetition time (TR) = 2400 ms, flip angle = 8°, 1 x 1 x 1 mm isotropic voxels) and T2-weighted (TE = 84 ms, flip angle = 120°, 1 x 1 x 4 mm voxels) scans of the same subject collected from a Siemens Trio (Erlagen, Germany) 3T scanner at Washington University School of Medicine. In-house automated algorithms, which involve a series of iterative thresholding, region growing, and masking technique for tissue specification, were applied to both MRI data sets. Tissue optical properties assigned to this heterogeneous head model were values used in previous in vivo study [25] at 750 nm and 850 nm (Table 1 ), which are the two wavelengths equipped in the current HD-DOT system at Washington University School of Medicine [38] and can be adapted to other and more wavelengths.

 

Fig. 8 (a) Posterior surface rendered view of the 3D FEM head model, (b) segmented tissues shown on an axial slice taken through the 3D model noted by the black solid line in (a), with scalp, skull, CSF, gray and white matter indicated from dark to light in grayscale, (c) posterior and (d) lateral schematic view showing the placement of the high-density imaging array over the visual cortex with 24 sources (red squares) and 28 detectors (blue circles), (e) first to third nearest neighbor measurement (1NN-3NN) definitions with separation of 13, 30, and 40 mm respectively.

Download Full Size | PDF

Table 1. Head tissue optical properties at 750 nm and 850 nm

View Table

A high-density (HD) imaging array which contains 24 sources and 28 detectors as described previously [9,38], was placed on the back of the head over the visual cortex, Figs. 8(c) and 8(d). Within this optode arrangement, nearest neighbor measurements can be defined based on source-detector separations, Fig. 8(e). In this study, 260 differential measurements per wavelength were used to image the hemodynamic change of human brain function, which include first, second and third nearest neighbor measurements, corresponding to a maximum source-detector separation of 40 mm.

Instead of presenting a ΔHbO₂ concentration target as in the 2D circular model, here we modeled a ΔHbR concentration target for the 3D head model. To generate the forward data, we simulated a focal activation of ΔHbR = −0.01mM and 1 cm in radius on the right hemisphere of the visual cortex (Figs. 9(a) and 9(b)) as the type of brain activation one would expect from a retinotopic mapping study [38]. Gaussian random noise of 0.1%, 0.14% and 1% in amplitude was added to first, second and third nearest neighbor measurements respectively to mimic our current in vivo performance [39]. Similar to the current in vivo data collection strategy (whereby data is collected at a sampling rate of ~10 Hz and then block averaged), ten sets of noise added data were generated and individually reconstructed, and averaged to produce the final image. Image reconstruction was performed with a ‘whole brain’ structural constraint that limits the recovered activation on the gray and white matter only [11,40]. The optimal values of regularization were ${\alpha }_{\lambda }$ = 10⁻² × the maximum singular value of $J_{\lambda }$, which were found to provide good imaging quality based on our previous human [38] and animal [41] DOT studies. The conventional spectral method was excluded in this part of the study due to earlier shown severe crosstalk between chromophores, which was also seen but not presented.

 

Fig. 9 (a) Posterior surface rendered view of the 3D FEM brain model, (b) a regional field of view (FOV) focused on the right hemisphere of visual cortex enclosed within the black window in (a), showing simulated chromophore target, (c) reconstructed images of ΔHbO₂ (upper row) and ΔHbR (lower row) concentration for ‘Svd-Spec’ and ‘Non-Spec’ method using noise free data, and (d) using noise added data.

Download Full Size | PDF

Reconstructed chromophore concentration images of simulated brain activation using non-spectral (‘Non-Spec’) and spectral method with SVD-based regularization (‘Svd-Spec’) with noise free and noise added data are displayed on the surface rendered FEM model of the brain in Fig. 9. Similar to our finding with the 2D model, a region of crosstalk of ΔHbR into ΔHbO₂ is identified when noiseless data are used for image reconstruction, with the magnitude of crosstalk from ‘Svd-Spec’ significantly smaller than ‘Non-Spec’ (Fig. 9(c)). The crosstalk (as defined in the 2D case) in Fig. 9(c) is $4.9\times {10}^{-5}$ for ‘Svd-Spec’, representing a 95% reduction from $1.1\times {10}^{-3}$ for ‘Non-Spec’. When measurement noise based on a current imaging system specification is added, the two methods produce images that are similar with crosstalk both at approximately $2\times {10}^{-3}$ (Fig. 9(d)) showing little difference between non-spectral and spectral reconstruction techniques.

4. Discussion

In this study we have presented a new regularization technique for linear image reconstruction in spectral diffuse optical tomography. This technique utilizes the singular value decomposition (SVD), allowing regularization on the spectral Jacobian matrix to be applied without altering the underlying spectral prior information. Specifically, this is achieved by regularizing the wavelength dependant Jacobian matrix prior to mixing with the spectral coefficients, a completely new method where the regularization does not affect the mapping of absorption changes onto chromophore concentrations, thereby providing a much more accurate spectral inversion technique. This method is different to previously published work, whereby the Jacobian is first spectrally mixed and then regularized and inverted.

Through a series of proof-of-concept analyses using a 2D circular model, we have shown that the SVD-based regularization technique followed by spectral mixing dramatically reduces the crosstalk between chromophores. In comparison to conventional spectral method, the crosstalk is reduced by approximately 98%. When compared with non-spectral method, this technique also demonstrates 60% crosstalk reduction with noise free data and demonstrated consistently robust performance against crosstalk as noise increases, before the image is dominated by imaging artifacts due to poor SNR (Fig. 6(c)). We extended our evaluation analysis on a 3D subject specific head model with simulated focal activation on the visual cortex, with the results showing 95% reduction in crosstalk when using the proposed algorithm as compared to non-spectral method.

In the case of spectral techniques, using conventional regularization (Eq. (8)) the regularization parameter is applied to a spectral Jacobian (one that relates a small change in measurements to a change of either oxyhemoglobin or deoxyhemoglobin concentrations) and therefore any smoothing applied will also have an effect on the spectral mixing. However, using the proposed SVD-based spectral method (Eq. (12)) the spectral Jacobian is a pre-regularized matrix (which only relates a small change in measurements to change in absorption) and the spectral mixing is unaffected by this regularization, therefore implying that the latter method should provide better crosstalk in parameter recovery. In the case of noise-added data, it is possible to expand Eq. (12) by adding noise to the final term on the right-hand-side and then demonstrate that at a certain noise level, the improvement in crosstalk are gradually minimized as the effect of noise on the (un-regularized) spectral matrix becomes dominant.

Validation of the SVD-based method as presented in this paper has been limited at a specific wavelength pair, i.e. 750 and 850 nm. In this case the crosstalk between ΔHbR and ΔHbO₂ is systematically small, owing to the large and opposite difference between their molar extinction coefficients at 750 and 850 nm (Fig. 1). In the concept of wavelength optimization [26], this means that the condition number of their spectral composition matrix $M_s$ is relatively low as compared to other wavelength pairs. Nevertheless since the crosstalk improvement by our SVD-based method relies on no assumption on either specific spectral range or specific combinations of the molar extinction coefficients, we expect the improvement to be consistent when using other pairs of wavelengths. A full validation of this could be carried out over an extensive range of wavelength pairs as an extension of this work.

When compared with the linear, one-step DOT as presented here, non-linear DOT reconstruction has the added benefit of iteratively updating the recovered chromophore concentrations simultaneously with a minimization function until the solution finally converges to a predefined threshold. The improvement at each iteration is therefore accumulative and the final recovered imaged are shown to provide better contrast and robustness to noise. However, in the case of dynamic (or temporal or differential) imaging, as shown here, appropriate regularization is shown to be critical to minimize crosstalk while maintaining the benefits of spectral constraints. We expect that in the presence of larger magnitude differential data and lower noise levels, this proposed reconstruction algorithm to be superior to non-spectral techniques. Finally, we expect that in line with other published works in spectrally constrained non-linear imaging techniques, the addition of more wavelengths will provide better contrast and minimize crosstalk further.

5. Conclusion

While incorporating spectral prior information in spectral image reconstruction correlates the absorption coefficients across multiple wavelengths, providing a solution that is stable, the use of conventional regularization techniques can result in additional crosstalk between chromophores. A novel regularization technique that regularizes the Jacobian using the singular value decomposition (SVD) prior to spectral mixing has been presented in an attempt to reduce this crosstalk effect. Simulations have shown that using SVD-based regularization can dramatically reduce the crosstalk presented in images recovered using conventional regularization technique in linear, single-step spectral DOT. Although our analysis has shown that given the current HD-DOT imaging system specification, such improvement in crosstalk over the non-spectral method may not be substantial, it is evident that future instruments with higher SNR measurements would only yield better image quality for the spectral method over the non-spectral method. The use of SVD in matrix regularization as described in this study is also potentially applicable to other one-step linear imaging problems, offering an alternative approach to the conventional Tikhonov regularization and should also play an important part in non-linear spectral imaging techniques in DOT.