Extraction of tissue optical property and blood flow from speckle contrast diffuse correlation tomography (scDCT) measurements.

Measurement of blood flow in tissue provides vital information for the diagnosis and therapeutic monitoring of various vascular diseases. A noncontact, camera-based, near-infrared speckle contrast diffuse correlation tomography (scDCT) technique has been recently developed for 3D imaging of blood flow index (αDB) distributions in deep tissues up to a centimeter. A limitation with the continuous-wave scDCT measurement of blood flow is the assumption of constant and homogenous tissue absorption coefficient (μ a ). The present study took the advantage of rapid, high-density, noncontact scDCT measurements of both light intensities and diffuse speckle contrast at multiple source-detector distances and developed two-step fitting algorithms for extracting both μ a and αDB. The new algorithms were tested in tissue-simulating phantoms with known optical properties and human forearms. Measurement results were compared against established near-infrared spectroscopy (NIRS) and diffuse correlation spectroscopy (DCS) techniques. The accuracies of our new fitting algorithms with scDCT measurements in phantoms (up to 16% errors) and forearms (up to 23% errors) are comparable to relevant study results (up to 25% errors). Knowledge of μ a not only improved the accuracy in calculating αDB but also provided the potential for quantifying tissue blood oxygenation via spectral measurements. A multiple-wavelength scDCT system with new algorithms is currently developing to fit multi-wavelength and multi-distance data for 3D imaging of both blood flow and oxygenation distributions in deep tissues.


Introduction
Blood flow (BF) facilitates the delivery of nutrients and oxygen to tissues and removal of metabolic byproducts from tissues. Continuous monitoring of BF in tissue provides information to medical professionals for the diagnosis and therapeutic monitoring of various diseases associated with tissue ischemia, such as peripheral arterial disease, cerebral vascular disease, and various types of wounds [1,2]. Tissue oxygenation level, on the other hand, reflects the balance between oxygen supply (via BF) and oxygen consumption, and thus is an important biomarker for tissue hypoxia associated with tissue injury [3,4]. Furthermore, simultaneous monitoring of BF and tissue oxygenation allows for the estimation of tissue oxygen consumption, another crucial biomarker for tissue pathophysiological conditions [3][4][5].
Use of a remote wide field illumination and 2D detection by a charge-coupled device (CCD) or complementary metal-oxide-semiconductor (CMOS) camera achieves fast, high-resolution, 2D mapping of tissue hemodynamics in a noncontact manner. For example, conventional laser speckle contrast imaging (LSCI) [16,17,37] and optical intrinsic signal imaging (OISI) [38] enable noncontact 2D mapping of blood flow and oxygenation distributions, respectively. However, both LSCI and OISI are mainly sensitive to the tissue surface with limited penetration depths (< 1 mm), therefore they are inadequate for noninvasive imaging of deep tissues hemodynamics.
A limitation with continuous-wave measurements of BF is the assumption of constant and homogenous µ a (λ) and µ ′ s (λ) (reduced scattering coefficient) from the literature [48][49][50]. Although this assumption is widely accepted in DCS/DCT, inter-subject heterogeneity and potential dynamic variations in µ a (λ) and µ ′ s (λ) may cause errors in calculating BF variations [51][52][53][54]. It has been shown that fitting DCS data of auto-correlation functions (g 1 ) at multiple S-D distances and multiple wavelengths enables simultaneous extractions of µ a (λ), µ ′ s (λ) and blood flow index (BFI) for a homogeneous tissue volume [51]. Because of the intrinsic relationship between diffuse speckle contrast and auto-correlation function, it is possible to extract BFI and optical properties from diffuse speckle contrast measurements at multiple S-D distances and/or at multiple exposure times. For example, Valdes et al. performed speckle contrast optical spectroscopy measurements at multiple S-D distances (5 to 20 mm with an increment of 2.5 mm) to quantitatively fit for BFIs in tissue-simulating phantoms and human forearms [53]. Similarly, Dragojevic et al. conducted speckle contrast optical tomographic reconstructions with multiple S-D distances (up to 6.5 mm) for quantification of BFIs in rodent's brains [55,56]. More recently, Liu et al. explored extracting µ a (λ), µ ′ s (λ), coherence factor (ß), and BFI from diffuse speckle contrast analysis at multiple S-D distances (2.5 to 10.5 mm with an incremental step of 1 mm) and with multiple exposure times (0.1 to 1 ms with an increment of 0.1 ms) [52]. The spectroscopic measurements (not imaging) were made in a semi-contact manner with a source fiber in contact with the tissue phantom. The fitting accuracy depended highly on data quality, which required long-time data acquisition and averaging: 1000 frames at each exposure time. Their system was tested in tissue phantoms with known optical properties and data acquisition time was tens of minutes for each step of multi-distance and multi-exposure measurements [52].
To overcome the limitations of existing systems, we took the advantage of rapid, high-density, noncontact scDCT measurements of both light intensities and diffuse speckle contrasts at multiple S-D distances and developed new algorithms for extracting both µ a (λ) and BFI. The noncontact hardware of scDCT makes data acquisition easy and fast with a high sampling density. scDCT usually scans a point light via an electronically-controlled galvo mirror (switching time: 1 ms) to multiple source positions and takes one or more images by the sCMOS camera at each source position with a single exposure time of a few milliseconds per image [40][41][42][43][44][45][46][47]. The fast and high-density sampling by the sCMOS camera sensor array (e.g., 2048 × 2048 pixels) enables rapid tomographic imaging of tissue hemodynamic distributions with improved signal-to-noise ratio (SNR) (via spatial data averaging) [39][40][41][42][43][44][45][46]. The optimized scDCT acquires one scanning over multiple source positions (e.g., 9 × 9) within ∼1.5 minutes with numerous S-D pairs at larger distances up to 20 mm, which benefits larger penetration depths. In the in-vivo test, scDCT is used to collect the data from the ROI on human forearm. (c) An illustrative raw intensity image taken from the selected ROI. The point source locates at the center of ROI. (d) A zoom-in view of a small area around the source. A pixel-window at the S-D distance of r 0 is defined in a cluster with adjacent 50 pixels marked with the same color (blue or red). Note that pixel windows in all directions around the source r 0 are used for data analysis.
The focus of the present study was to demonstrate the feasibility and accuracy of our new two-step fitting algorithms for the extraction of both µ a and αD B from scDCT measurements. Therefore, only one source located at the center of the ROI was used in this study for proof of the concept. Data acquisition time at one source position was less than 10 seconds. The new fitting algorithms were tested in tissue-simulating phantoms with known optical properties and human forearms. Measurement results were compared against established NIRS/DCS techniques. Knowledge of µ a (λ) not only improved the accuracy in calculating BFI but also provides the potential for quantifying tissue blood oxygenation via spectral measurements [36,57,58]. As mentioned earlier, oxygenation information can be determined by measuring µ a (λ) at multiple wavelengths using NIRS principles.
A sCMOS camera (Orca Flash 4.0LT, Hamamatsu, Japan) was used for data collection at 25 frames per second. A zoom lens (Zoom 7000, Navitar, NY) was mounted to the camera enabling easy adjustment of ROI and F number. The camera detected diffuse light from the tissue surface in a reflectance mode with a typical exposure time of 2 ms. A pair of polarizers (LPNIRE050-B and LPNIRE200-B, Thorlabs, NJ) were added across the source and detection paths to reduce direct reflection from sources. A long-pass filter (> 800 nm, 84-762, Edmund Optics, NJ) was added in front of the zoom lens to minimize the influence of ambient light.
The speckle contrast is defined as: Here, σ s and ⟨I⟩ are the spatial standard deviation and mean value of light intensities over a pixel window of 50 pixels, respectively ( Fig. 1(d)). r is the distance between the centers of the pixel-window and source. The diffuse speckle contrast K s (r) can be related to the temporal auto-correlation function g 1 (r, τ) via the following equation [39]: Here, T is the exposure time and τ is the correlation lag time. By substituting g 1 (r, τ) with its analytical solution in semi-infinite geometry, a nonlinear relationship between K s (r) and the BFI (αD B ) was derived in our previous publication [39]. K s (r) can be then written as a function of αD B , together with other parameters: Here, λ is the wavelength of the source, β is the coherence factor, and k 0 is the wavenumber. αD B is a combined term representing BFI, where α is a unitless ratio of dynamic scatterers to total scatterers and D B (unit: cm 2 /s) is the effective diffusion coefficient of moving scatterers.

New algorithms for extracting µ a and αD B
As K s (r) is a function of multiple parameters including αD B , µ a (λ) and µ ′ s (λ) Eq. (3), assumptions of optical properties (µ a (λ) and µ ′ s (λ)) from literature may lead to errors in αD B calculations.
Thus, a new approach with two fitting steps is proposed here to extract both µ a and αD B from scDCT measurements: (1) µ a (λ) is extracted by fitting the measured light intensities (I) at multiple S-D distances ( Fig. 2(a)) and (2) αD B is assessed by fitting the measured speckle contrasts (K s (r)) at multiple S-D distances ( Fig. 2(b)) with the input of known µ a (λ) obtained from the first step. The influence of µ ′ s (λ) assumption is discussed in Section 4. In the first step, the simplified solution of photon diffusion equation is given by [54,59,60]: Here, I(r) is the light intensity measured at the S-D separation r, and µ eff is the effective attenuation coefficient defined as: Based on Eq. (4), ln [I(r) r 2 ] has a linear relationship with r [54,59,60]. Thus, µ eff (λ) can be extracted by fitting the slope of logarithmic light intensities over multiple S-D separations. By assuming µ ′ s (λ) from literature or acquiring µ ′ s (λ) from other measurements, µ a (λ) can be calculated from µ eff (λ) (Fig. 2(a)) [48][49][50]. The second step is to fit αD B from the measured speckle contrasts at multiple S-D distances by minimizing the penalty term [52,53,61,62]: Here, K 2 m is the experimentally measured speckle contrasts and K 2 t is the theoretical speckle contrasts with an initial guess of αD B , calculated using Eq. (3). To eliminate the influence of linear scaling factor β in Eq. (3), both K 2 m and K 2 t are normalized to their mean values over all the selected r, respectively ( Fig. 2(b)). The normalization cancels the constant β from calculations [21,53].
Since there is only one wavelength (λ = 830 nm) involved in the experiments in this work, we write µ a (λ) and µ ′ s (λ) as µ a and µ ′ s from here on for simplicity.

Phantom tests with varied µ a and αD B
To examine performance of the new algorithms, tissue-like liquid phantom tests were first performed ( Fig. 1(a)). A glass tank was filled with a liquid solution that consisted of distilled water, India ink (Black India, MA) and Intralipid (Fresenius Kabi, Sweden) [40,63]. µ a and αD B were varied via ink titration and temperature variation of the liquid phantom, respectively. A FOV of 90 × 90 mm 2 was selected in the phantom center for scDCT measurements, with a single source position fixed in the center of the FOV (Fig. 1(c)). Camera exposure time was set as 2 ms and F# was 8. At each level of µ a and αD B , 50 frames were taken by scDCT within 2 seconds for averaging to improve SNR. For comparison purpose, µ a and αD B were also measured using an established hybrid instrument combining a commercial NIRS device (Imagent, ISS, IL) for µ a and a customized DCS device for αD B [57,58,63,64]. NIRS and DCS data were taken by a hybrid fiber-optic probe placing on the phantom surface immediately after scDCT measurement at each step.
Ink titration to change µ a . µ ′ s and αD B were set up and kept as 10.2 cm −1 and around 1.15×10 −8 cm 2 /s throughout titration steps. Indian ink was added step-by-step to create gradient increases in µ a (λ = 830 nm) from 0.04 cm −1 to 0.16 cm −1 with 0.03 cm −1 increment per step. At each titration step, the phantom was stirred thoroughly and then rested for 2 minutes to allow for stabilization before data collection by scDCT and NIRS/DCS. Temperature variation to change αD B . µ a and µ ′ s were set up and kept as 0.1 cm −1 and 10.2 cm −1 (λ = 830 nm), respectively, during temperature changes to create αD B variations. Brownian motions of Intralipid particles (αD B ) were varied via temperature changes of the phantom from ∼5°C (fridge temperature) to ∼25°C (room temperature). Data were collected at each increment of 5°C using scDCT and NIRS/DCS devices, sequentially. Phantom temperature was continuously monitored by a thermometer throughout the test.

In-vivo tests in human forearms
To validate new algorithms for in-vivo tests, resting baseline measurements were performed in five human forearms. This study was approved by the University of Kentucky Institutional Review Board. The subject's forearm was secured on a table with an arm support ( Fig. 1(b)). A FOV of 60 × 60 mm 2 was selected to adapt to the dimensions of forearm, while other settings were kept the same as in the phantom tests. In total, 250 frames were taken by noncontact scDCT within 10 seconds from the resting forearm. For comparison purpose, µ a and αD B were then measured by a contact hybrid NIRS/DCS probe placing on the same area.

Data analysis
Dark noise and shot noise were first removed from the raw scDCT images using the methods reported previously [39][40][41][42][43][44][45][46]. Spatial speckle contrasts K s (r) and light intensities at all directions around the source were calculated ( Fig. 1(d)). Different ranges of S-D distances were used for data analysis to adapt different sizes of FOV: 0.25 to 35 mm for phantoms and 0.25 to 20 mm for human forearms, with an increment of 0.5 mm. For each S-D distance r 0 , all pixels within the distance range of (r 0 -0.2) mm ≤ r ≤ (r 0 + 0.2) mm around the source were selected and segmented into pixel windows consisted of 50 adjacent pixels ( Fig. 1(d)). The speckle contrast K s (r) at each pixel window was calculated using Eq. (1). After excluding 50% outlier data (the highest and lowest 25%), light intensities and speckle contrasts at r 0 were then averaged to obtain mean values of I(r 0 ) and K s (r 0 ), respectively. Note that the outlier data removal is primarily used to reduce the influence of spatial in-vivo tissue heterogeneity. In addition to the spatial averaging, time-course data of I(r 0 ) and K s (r 0 ) were averaged over 50 frames to further improve SNRs.
Based on experimental results (see Section 3), effective S-D distances were determined for extracting µ a and αD B using the two-step fitting algorithms Eq. (4)-6). A constant µ ′ s was set up in all tissue phantoms based on the phantom recipe, while µ ′ s values in human forearms were measured by the commercial NIRS device (Imagent).
The results from scDCT and NIRS/DCS measurements were compared and percentage errors were calculated by taking NIRS/DCS measurements as standards. Linear regressions were conducted and Pearson correlation coefficients were calculated to evaluate correlations between the scDCT and NIRS/DCS measurements. Paired t-tests were used to test group differences between the two measurements. p < 0.05 was considered a statistically significant correlation between the two measurements.  b and d). For visualization and comparison, I(r) data were normalized to the source location (r = 0 mm) and K s 2 (r) data were normalized to the mean value of K s 2 (r) in the selected effective S-D ranges. Although variations existed, data from the phantom and human forearm had similar patterns.   c and d) to the models Eq. (4)-6), respectively. Based on our measurement results from all tissue phantoms and human forearms, effective S-D ranges used in this study were determined and listed in Table 1.   Figure 5 shows comparison results from tissue-simulating phantoms during ink titration (a to c) and temperature variation (d to f), measured by the scDCT and NIRS/DCS devices. During ink titration, significant correlations were observed between the scDCT and NIRS measurements of µ a changes (linear regression: R 2 > 0.99, p < 10 −4 ; a and c). Although αD B varied slightly during ink titration, no significant difference was observed between the scDCT and DCS measurements of αD B (paired t-test: p = 0.36; b). The small variation in αD B might result from the influence of small room temperature variation on Intralipid particle motions in the phantom during the test. When creating a large temperature variation in the phantom, significant correlations were observed between the scDCT and DCS measurements of αD B changes (linear regression: R 2 = 0.93, p < 0.01; e and f). No significant differences were observed between the scDCT and NIRS measurements of µ a (paired t-test: p = 0.62; d). Table 2 lists all measurement results from tissue phantoms by the scDCT and NIRS/DCS devices. scDCT measurement errors were less than 16% (assuming NIRS/DCS as the gold standard), except one large outlier error in αD B (45.42%) at the last step of ink titration. Figure 6 shows the results from five forearms, measured by the scDCT and NIRS/DCS devices. There was no significant difference between the mean values of µ a (paired t-test: p = 0.83; a), measured by the two devices. A significant correlation was observed between scDCT and NIRS measurements of µ a (linear regression: R 2 = 0.88, p = 0.018; b). Similarly, there was    Table 3 lists all measurement results from human forearms by the scDCT and NIRS/DCS devices. Note that individual αD B was calculated using individual forearm µ ′ s , measured by the NIRS device. scDCT measurement errors were less than 23%, assuming NIRS/DCS as the gold standard.

Discussion and conclusions
This study developed a new approach with two fitting steps to extract both µ a (λ) and BFI (αD B ) from innovative scDCT measurements of both light intensities and diffuse speckle contrasts at multiple S-D distances. The assessment of µ a from the measured light intensities improved the accuracy of αD B calculations from the measured speckle contrasts. Moreover, knowledge of µ a (λ) at multiple wavelengths would enable quantification of tissue blood oxygenation [36,57,58]. The noncontact camera-based hardware of scDCT makes data acquisition easy and fast with a high sampling density, thus allowing for tomographic imaging [39][40][41][42][43][44][45][46][47]. Importantly, the fully noncontact measurements improve the versatility of scDCT in clinical applications where contact measurements may result in an infection risk or hemodynamic disruption on soft, vulnerable tissues. The scDCT measurements with new fitting algorithms were examined in tissue-simulating phantoms with known optical properties and human forearms against standard NIRS/DCS measurements. We selected S-D distance ranges for light intensity and speckle contrast analyses based on satisfactory of data fitting in ln [I(r) r 2 ] and K s 2 (r) to our models Eq. (4-6, and Fig. 4), respectively. As a result, the effective S-D ranges for light intensity and speckle contrast analyses were slightly different (Table 1). In the future, we may explore using the overlapped S-D range for both light intensity and speckle contrast analyses.
Phantom test results measured by the scDCT and NIRS/DCS devices were correlated and agreed well with expected values of µ a and αD B based on the phantom recipe over a relatively wide range of variations, which are commonly observed in biological tissues. Most of scDCT measurement errors in µ a and αD B were less than 16% (assuming NIRS/DCS as the gold standard), except one large outlier error in αD B (45.42%). Compared to phantom test results, in vivo test results in human forearms demonstrated slightly larger measurement discrepancies between the NIRS/DCS and scDCT measurements (up to 23%). These discrepancies might be partly attributed to temporal variation of µ a and αD B since the two measurements were taken sequentially. Tissue spatial heterogeneity might also contribute to the discrepancy as the two measurements covered different regions with different sizes/areas. In addition, NIRS/DCS probe-tissue contact pressure may result in tissue hemodynamic variations. The accuracies of our new fitting algorithms with scDCT measurements of multiple parameters are comparable to relevant study results. Up to 25% errors in fitting µ a and αD B were observed during DCS measurements at multiple S-D distances and multiple wavelengths and diffuse speckle contrast measurements at multiple S-D distances and multiple exposure times in manipulated tissue phantoms [51,52].
We note that use of a longer acquisition time for each measurement would improve SNR via temporal data averaging. However, it also impacts the temporal resolution of scDCT, especially when using multiple source positions for tomographic measurements. Moreover, possible physiological variations in the measured tissue over a long acquisition time may bias data collection.
In addition, we have also explored simultaneous extraction of µ a , µ ′ s and αD B from scDCT measurements of light intensities and speckle contrasts at multiple S-D distances. The fitting accuracy was very sensitive to data quality and SNR (data are not shown). The crosstalk among µ a , µ ′ s and αD B was also observed, as founded in previous studies [52,54,61]. To compare fitting sensitivities, we generated two noise-free reference curves (blue color): ln [I(r) r 2 ] in Fig. 7(a) and K s 2 (r) in Fig. 7(b), with the given parameters: µ a = 0.1 cm −1 , µ ′ s = 8 cm −1 , αD B = 1 × 10 −8 cm 2 /s, r = 4 to 15 mm. We then varied µ a from 0.04 cm −1 to 0.16 cm −1 to generate more testing curves. Apparently, variations in ln [I(r) r 2 ] resulting from µ a changes were much larger than in K s 2 (r), indicating a better sensitivity of ln [I(r) r 2 ] to µ a changes ( Fig. 7(a) and Fig. 7(b)). By contrast, K s 2 (r) signals are intrinsically sensitive to BFI. Therefore, our two-step fitting approach generated more robust results in µ a (by fitting ln [I(r) r 2 ]) and BFI (by fitting K s 2 (r)) than simultaneous fitting algorithms.
We recognize limitations of this pilot study in exploring new fitting algorithms for data analysis. The small number of subjects (n = 5) impacts the study power, which will be addressed in future studies. The effective S-D distances were determined manually based on the optimal fitting of experimental data to the models ( Fig. 4 and Table 1). More experiments in a larger number of subjects are needed to generate a universally applicable method for automatic determination of effective S-D distances for data fittings. In this study, µ ′ s values were obtained from the phantom recipe or in vivo NIRS measurements. As mentioned earlier, simultaneous fitting of µ a , µ ′ s and αD B using the current scDCT data generated unstable results with crosstalk. In the future, we will explore adding more wavelengths into the scDCT system and fit multi-wavelength and multi-distance data to extract all parameters simultaneously [51]. Finally, we will also explore 3D imaging of both blood flow and oxygenation distributions by developing a multiple-wavelength scDCT system.

Disclosures. The authors declare no conflicts of interest.
Data availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.