{'fraction_non_alphanumeric': 0.03412759415833974, 'fraction_numerical': 0.015987701767870868, 'mean_word_length': 5.170500355703107, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Singular value decomposition based regularization prior to spectral mixing improves crosstalk in dynamic imaging using spectral diffuse optical tomography

The spectrally constrained diffuse optical tomography (DOT) method relies on incorporating spectral prior information directly into the image reconstruction algorithm, thereby correlating the underlying optical properties across multiple wavelengths. Although this method has been shown to provide a solution that is stable, the use of conventional Tikhonov-type regularization techniques can lead to additional crosstalk between parameters, particularly in linear, single-step dynamic imaging applications. This is due mainly to the suboptimal regularization of the spectral Jacobian matrix, which smoothes not only the image-data space, but also the spectral mapping space. In this work a novel regularization technique based on the singular value decomposition (SVD) is presented that preserves the spectral prior information while regularizing the Jacobian matrix, leading to dramatically reduced crosstalk between the recovered parameters. Using simulated data, images of changes in oxygenated and deoxygenated hemoglobin concentrations are reconstructed via the SVD-based approach and compared with images reconstructed by using non-spectral and conventional spectral methods. In a 2D, two wavelength example, it is shown that the proposed approach provides a 98% reduction in crosstalk between recovered parameters as compared with conventional spectral reconstruction algorithms, and 60% as compared with non-spectrally constrained algorithms. Using a subject specific multilayered model of the human head, a noiseless dynamic simulation of cortical activation is performed to further demonstrate such improvement in crosstalk. However, with the addition of realistic noise in the data, both non-spectral and proposed algorithms perform similarly, indicating that the use of spectrally constrained reconstruction algorithms in dynamic DOT may be limited by the contrast of the signal as well as the noise characteristics of the system.


Introduction
Diffuse optical tomography (DOT) is finding widespread application in the diagnosis and characterization of breast cancer [1][2][3][4][5], monitoring brain function [6][7][8][9][10][11][12][13][14], and small animal imaging for the study of disease detection, progression and treatment [15][16][17][18]. Specifically in DOT of human brain function, an imaging cap that uses an array of continuous-wave, nearinfrared (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 2 ) and deoxyhemoglobin (HbR). Given their difference in the nearinfrared 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 highdensity-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 dualwavelength setup using 750 and 850 nm.
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 unmixed 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][28][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 SVDbased 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.

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 and λ Φ at each wavelength λ which are acquired before and after a change of absorption, λ µ ∆ . Given that λ µ ∆ is relatively small, the forward problem can be linearized by the Rytov approximation (whereby the log of the intensity measurements are considered) [20] ,ref where J λ is the wavelength dependent Jacobian matrix for wavelength λ, with a size of number of measurements by number of nodes. Since J λ is non-invertible, the Moore-Penrose generalized inverse is applied ( ) where 2 λ α is the Tikhonov regularization parameter [20,35]. Assuming a dual-wavelength system and two chromophores of interest, we have where s M is the spectral composition matrix (containing spectral prior information) of the following form where , c λ ε 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.

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 s M can be directly incorporated into the forward problem by constructing a spectral Jacobian matrix s J By doing so the forward problem can be expressed as and the spectral image reconstruction using conventional Tikhonov regularization becomes ( ) where 2 s α 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 2 s α (in Eq. (8)) must be smaller than 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 nonspectral 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 λ but also the spectral prior information s M incorporated within the spectral Jacobian matrix s J , which results in numerical error in the form of crosstalk on the reconstructed images, as demonstrated later in Section 3.1.

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 λ , and S is a diagonal matrix that contains the i number of singular values of J λ , i.e., ,i λ σ . Since applying Tikhonov regularization of weight 2 We can then construct the regularized s J , denoted ˆs J , 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 s M is exactly preserved and therefore contains all information needed for the wavelength unmixing operation. In addition, 2 λ α in the non-spectral method is now directly used in the spectral method, avoiding the operation to find an equivalent 2 s α for comparative studies which can be problematic due to its subjective choice.

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 −1 , µ′ s = 0.74 mm −1 at 750 nm, and µ a = 0.019 mm −1 , µ′ s = 0.64 mm −1 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 2 = +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 ( s α in Eq. (8)) as compared to the non-spectral method ( λ α in Eq. (2)) and the proposed SVD approach ( λ α in Eq. (10)), which was in line with previous finding [30]. Specifically we chose  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 2 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.
For quantitative analysis, the crosstalk of ΔHbO 2 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 2 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 ~3 2 10 − × for 'Svd-Spec', representing a 60% reduction from ~3 5 10 − × for 'Non-Spec', and a 98% reduction from ~10 −1 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 3 5 10 − × , 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.

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 T2weighted (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. 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.
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 2 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 λ α = 10 −2 × the maximum singular value of J λ , 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.
Reconstructed chromophore concentration images of simulated brain activation using nonspectral ('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 2 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 5 4.9 10 − × for 'Svd-Spec', representing a 95% reduction from 3 1.1 10 − × 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 3 2 10 − × ( Fig. 9(d)) showing little difference between non-spectral and spectral reconstruction techniques.

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 preregularized 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 2 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 s M 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.

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.