Compressed sensing time-resolved spectrometer for quantification of light absorbers in turbid media.

Time-resolved (TR) spectroscopy is well-suited to address the challenges of quantifying light absorbers in highly scattering media such as living tissue; however, current TR spectrometers are either based on expensive array detectors or rely on wavelength scanning. Here, we introduce a TR spectrometer architecture based on compressed sensing (CS) and time-correlated single-photon counting. Using both CS and basis scanning, we demonstrate that-in homogeneous and two-layer tissue-mimicking phantoms made of Intralipid and Indocyanine Green-the CS method agrees with or outperforms uncompressed approaches. Further, we illustrate the superior depth sensitivity of TR spectroscopy and highlight the potential of the device to quantify absorption changes in deeper (>1 cm) tissue layers.

The concentration of a light absorber (i.e., chromophore) in a sample is directly related to the sample's absorption coefficient; however, in highly scattering media the optical pathlength is much longer than the distance between the points of light's entry into and exit out of the medium. Thus, the classic approach of quantifying chromophore concentrations with the Beer-Lambert law cannot be readily applied. Furthermore, when multiple chromophores are present-as is typical in living tissues-multi-wavelength and even hyperspectral capabilities may be required for accurate chromophore quantification. Though it is possible to quantify chromophores using only a few wavelengths, acquiring measurements at dozens of wavelengths enables the use of spectral features to better constrain spectroscopic analysis and reduce crosstalk between various chromophores, while also enabling the use of approaches such as derivative spectroscopy [17][18][19][20].
Optical spectroscopy techniques used to measure the concentration of light absorbers in turbid media can, despite their diversity, be divided into three categories: continuous-wave (CW), frequency-domain (FD), and time-resolved (TR). Of these techniques, TR methods -where picosecond light pulses are injected into the sample-have the highest information content and are considered the "gold-standard" for chromophore quantification in turbid media. Further, since TR techniques can quantify the time-of-flight of photons, they allow the use of time-windowing to separate contributions from different layers of the medium [21][22][23]. Nevertheless, these advantages come at the cost of added complexity and expense, particularly for hyperspectral TR systems.
Hyperspectral TR devices are typically built using a pulsed supercontinuum laser as a light source and a streak camera [24][25][26][27][28] or a time-correlated single-photon counting (TCSPC) module [19,[29][30][31] for detection. In TCSPC measurements, the intensity of the light source is typically adjusted such that the probability of detecting one photon per light pulse is substantially less than one, and the probability of detecting more than one photon is negligible; as such, TCSPC operates under the assumption that only one photon is detected per signal period [32]. Sensitive detectors-such as photomultiplier tubes (PMT), single-photon avalanche photodiode (SPAD) arrays, or silicon photomultipliers (SiPM)-and high-precision timing electronics are then used to measure the arrival time of each detected photon in order to build a histogram of arrival times (i.e., a distribution of times-of-flight of photons; DTOF). Previous implementations of TCSPC-based hyperspectral TR-NIR spectroscopy were based on wavelength scanning [19,30]; however, for biomedical applications, the broad spectral features of the main chromophores of interest (e.g., hemoglobin, water, fat) suggest that the measurements are sparse in the spectral domain and might be acquired more efficiently using compressive sensing (CS) [33]. CS leverages the sparsity of a signal to measure it with fewer samples than what is required using traditional data acquisition methods based on the Shannon-Nyquist sampling theorem [34][35][36]. This drastically reduces both the acquisition time and the size of the dataset. Given these benefits, the combination of CS and diffuse optical imaging has been used in a wide variety of biomedical applications, as summarized by Angelo et al. [37]. Notably, Pian et al. recently presented a hyperspectral TR fluorescence lifetime imaging system based on CS and TCSPC [38]. In their report, Pian et al. utilized CS for compression in the spatial domain; however, the application of CS to the spectral domain for TR-NIR spectroscopy remains unexplored.
Leveraging recent advances in TCSPC and CS, this study presents the development and characterization of a hyperspectral TR spectrometer based on a CS-capable architecture. The spectral accuracy and resolution of the spectrometer were characterized using a calibrated spectral source. Next, the spectrometer's ability to accurately quantify the absorption coefficient of turbid media over a wide spectral range was demonstrated in tissue-mimicking phantoms made of Intralipid and Indocyanine Green (ICG). ICG is an optical tracer with distinct spectral features in the NIR window of 650 -950 nm and is often used as a light absorber in tissue-mimicking phantoms [39][40][41][42]. Further, ICG is FDA approved and has been used in a variety of clinical applications [43], including assessment of hepatic function [44,45], evaluation of skin-flap viability [46,47], and quantitative perfusion assessment (e.g., dynamic contrast-enhanced nearinfrared spectroscopy) [48][49][50][51][52][53][54][55]. For future reference, we report the extinction coefficient spectrum for ICG in Intralipid over a broad spectral range (710 -830 nm). Further, through the addition of ICG to the bottom compartment of a two-layer phantom, we highlight the TR spectrometer's ability to quantify absorption changes in deeper (>1 cm) tissue layers.

Spectrometer design
A schematic of the TR spectrometer is shown in Figure 1. The output of a pulsed supercontinuum laser (SuperK Extreme K94-120-02, NKT Photonics, Denmark) is transmitted through two spectral filters (FEL0700, FES0950, Thorlabs Inc., NJ) to limit the emission bandwidth to 700 -950 nm. The filtered beam is then attenuated by a neutral density filter (NDF; NDC-50C-4M-B, Thorlabs Inc.) and sent to a microscope objective (RMS10X, Thorlabs Inc.) that couples it into an emission probe (fiber bundle, D = 3.7 mm, NA = 0.55, Fiberoptics Technology Inc., CT). The emission fiber bundle is attached to a 3D-printed probe holder for positioning on a sample. A custom-made round-to-linear detection fiber bundle (D = 1.55 mm [round], 0.25 mm × 7.50 mm [linear], NA = 0.55, ScienceTech Inc., Canada), also fixed to the probe holder, receives diffusely reflected light from the sample and guides it to a concave holographic diffraction grating (f = 100 mm, 520 lines/mm, 16 nm/mm dispersion, ScienceTech Inc.) which resolves the spectral content of the reflected light. The spectrum produced by the diffraction grating is then sent to a digital micromirror device (DMD; DLi4130 .7" VIS XGA High-Speed Development Kit, active area W × H = 14.2 mm × 10.7 mm, Digital Light Innovations, TX). The DMD contains a 2D array of micro-mirrors which can be individually tilted to either reflect light towards ("ON" position) or away ("OFF" position) from the detector. The DMD is used to spatially encode the signal from the diffraction grating and direct it through a series of lenses (LB1309-B, LB1723-B, LA1134-B, LA1951-B, magnification M = 0.4, Thorlabs Inc.) that focus the light onto a PMT (H7422P-50, active area D = 5 mm, Hamamatsu Photonics, Japan) coupled to a TCSPC module (SPC-130; Becker & Hickl GmbH, Germany). To avoid accidental detector saturation between measurements, a manually-controlled shutter is installed in front of the PMT.  1. Schematic of the hyperspectral TR system used in the study. Light from a pulsed supercontinuum laser is bandlimited by two filters (F1, F2), attenuated by a neutral density filter (NDF), and coupled into an emission fiber using a microscope objective (L1). Detected light is then dispersed by a diffraction grating and re-directed by a DMD into a series of lenses (L2 -L5) which focus the light onto a PMT coupled to a TCSPC module. The resulting output is a 3D TR spectrum.
For all measurements, the active array of the DMD was positioned to reflect the 700 -915 nm wavelength range from the diffraction grating. However, due to the limited area of both the shutter aperture and the entrance window of the photosensor module containing the PMT, adequate signal was only obtained in the 710 -830 nm range.

Data acquisition
The signal reflected from the diffraction grating is a time-dependent 2D distribution of light intensity I(x, y, t) whose wavelength λ only varies in x (Fig. 1). During data acquisition, preselected binary patterns were sequentially displayed on the DMD so that only a section of the spectrum I(x, y, t) was reflected towards the detector. For each pattern displayed on the DMD, a DTOF of the detected photons was measured. As a result, the full dataset produced by a single acquisition could be represented by a 2D matrix Y whose rows contained the DTOFs measured for each of the displayed patterns (i.e., each row contained one DTOF and only one DTOF was acquired per pattern). The sequence of binary patterns was chosen such that I(x, y, t) was spatially sampled only in the x direction. As such, the signal was binned in the y dimension and the final spectra were only functions of x and t. The simplest example of such a sequence of patterns is a stack where only a single column of mirrors in the y dimension is "ON" for any given pattern. This column is then shifted one unit in the x dimension for each subsequent pattern. With such a sequence, the spectrum would be wavelength-scanned and no additional image processing would be needed. More generally, the sampling of I(x, y, t) by an N-sized stack Φ of 1 × N sampling patterns (i.e., a measurement matrix Φ N×N ) can be represented as where the signal is spatially integrated over the height of the DMD (∆y) since the patterns only sample the spectrum in x.
One method to recover the spectrumÎ(x, t) from the data in Y is by performing a basis scan, i.e., choosing the sampling patterns which make up Φ from a known basis (e.g., Hadamard, Fourier, wavelet) [56][57][58][59]. Using this approach,Î(x, t) can be recovered through a simple inversion: While a basis scan offers simple image reconstruction, it typically suffers from lengthy acquisition times since, at a minimum, the number of acquired measurements must be equal to the number of desired data points. However, since NIR spectra acquired from tissue primarily contain broad, low-frequency spectral features, we expect that these spectra could be sparsely represented in the spectral domain. Under such conditions, a compressed sensing (CS) approach could be implemented where, for M ≪ N, an M-sized stack Φ of 1 × N sampling patterns (i.e., a measurement matrix Φ M×N ) is derived from a chosen basis and used to sample the spectrum. I(x, t) could then be recovered by solving a minimization-optimization problem [37,57]. In this article, both the full basis scan and compressed sensing approaches were investigated. For the compressed sensing approach,Î(x, t) was recovered using a total variation minimizing algorithm (TVAL3) [60], which is a common linear solver for CS signal recovery [57,59,61]. For the methods described above, the quality of signal recovery depends on the choice of basis used to generate the sampling patterns in Φ. We chose Φ to contain Hadamard basis patterns ranked from lowest to highest spatial frequency since, when used to reconstruct spatially-varying DTOF data, this basis has been shown to match or outperform a variety of other choices [59]. However, the implementation of Hadamard patterns has several practical challenges. First, Hadamard patterns are composed of ±1 elements. Given that it is not physically possible to directly implement both these element values onto a DMD, the measurements were acquired using a differential Hadamard single-pixel sampling approach [58]. Specifically, one measurement was acquired with a basis pattern in which −1 elements were replaced with 0's, a second measurement was acquired using the pattern's inverse, and the difference of the two measurements was taken. While this approach doubles the total number of measurements, it is also known to have superior noise suppression properties [58,62]. Second, there exists a boundary diffraction effect at the edges of "ON" mirrors. As such, the light energy collected from two adjacent mirror lines is not necessarily equal to two independent measurements of those lines and may lead to spectral inaccuracies following reconstruction. This inaccuracy can be mitigated by splitting every single pattern into two even-only and odd-only column components, which ensures that no two adjacent lines in any pattern will be "ON" [63].
Thus, the data acquisition protocol used 4N measurements to obtain one uncompressed acquisition of an N-point spectrum. Since acquisition time was not a major limitation for our application, spectral accuracy and noise suppression were prioritized. However, many other implementations-including different choices of Φ-can be used to optimize the system's performance if needed. Note that for all the results reported herein, we acquired 256-point spectrums (i.e., N = 256), which was a good compromise between achieving adequate spectral resolution and the number of required measurements (4N = 4 × 256 = 1024 measurements).

System characterization
Prior to characterization, the TR spectrometer was spectrally calibrated using a neon lamp. Four peaks with known wavelengths (724.5, 753.58, 794.32, and 811.5 nm) were identified in the measured spectrum and their wavelength values were fixed. The wavelengths of the remaining 252 points in the spectrum were then determined by linearly interpolating between these known wavelengths. As a result, each I(x, t) value was mapped to a corresponding I(λ, t) value. To validate this calibration procedure, the spectrum of the same neon lamp was measured with a calibrated commercial CW spectrometer (QE65000, Ocean Optics, FL) for comparison.
For system characterization, the DTOF collection time was set to 1 s, the laser power was adjusted until a count rate of 800 kHz was obtained, and an instrument response function (IRF) was acquired by placing a thin piece of white paper between the emission and detection probes. The IRF was then measured again after 9 hours of continuous system operation to assess the system's stability. Each IRF was first acquired using a full basis, and then re-acquired using a compressed basis built from the first 96 basis patterns (i.e., M = 0.375N for N = 256).

Phantom experiments
To test the ability of the TR spectrometer to measure the absorption coefficient of turbid media over a wide spectral range, experiments were conducted in tissue-mimicking phantoms using Intralipid-20%(Fresenius Kabi AG., Germany) as a scattering agent and ICG as a light absorber. The liquid phantom was placed in an opaque PVC container consisting of a shallow, removable top layer (10 × 10 × 4 cm, L × W × H) which was separated by a thin, transparent plastic film from a deeper bottom layer (10 × 10 × 7 cm, L × W × H). The bottom layer of the phantom contained an injection port and a magnetic stir bar to allow the addition of ICG to the Intralipid solution during an experiment without needing to remove the phantom's top layer. Before each set of phantom measurements, a neon lamp spectrum was acquired to calibrate the spectrometer as described in Section 2.3. After each phantom experiment, a second neon lamp spectrum was acquired to confirm the spectral stability of the system over the experimental period. In addition, the system's IRF was also acquired before and after the experiment to assess temporal stability.
For the first set of measurements, the top layer of the phantom was removed, and the bottom compartment was filled with a 0.8% Intralipid solution created by diluting Intralipid-20% with distilled water. The emission and detection probes were then fixed to the 3D-printed probe holder ( Fig. 1) at a distance of 28 mm. The holder was attached to a Manfrotto articulated arm (Vitec Imaging Solutions, Italy) to position and maintain the probe tips on the surface of the Intralipid solution. After acquiring a baseline measurement, the phantom's absorption coefficient was modified by adding ICG and mixing the solution for 60 s using a magnetic stirrer. The solution was then allowed to settle for 60 s before another measurement was acquired. This process was repeated for a total of 4 ICG additions (Table 1). Similar to the IRF measurements in Section 2.3, each phantom measurement acquired with a full basis was repeated using a compressed basis for comparison.
For the second set of measurements, the top layer of the phantom was filled with 0.8% Intralipid solution to a height of 12 mm, placed over the bottom layer, and each step of the protocol above was repeated. Note that the probe tips were positioned on the surface of the Intralipid solution in the top layer. For clarity, we will hereafter refer to the first and second set of measurements as homogeneous and two-layer phantom measurements, respectively.

Data analysis
Prior to reconstructingÎ(λ, t), each measured DTOF was denoised using a previously reported signal denoising algorithm [64]. Two different approaches were then used to reconstructÎ(λ, t) from the full and compressed acquisitions (see Section 2.2). These reconstructions were compared using the structural similarity (SSIM) index where C 1 and C 2 are regularization constants, and µ x , µ y , σ x , σ y and σ xy are the local means, the standard deviations, and cross-covariance for reconstructions x, y [65].
For the homogeneous phantom, absorption (µ a ) and reduced scattering (µ ′ s ) coefficients were estimated using the diffusion approximation [66]. Given that an experimentally measured DTOF can be expressed as a convolution between the system's IRF and the temporal point-spread function (TPSF) of the phantom (i.e., DTOF = IRF * TPSF), the IRF at each wavelength was convolved with the model solution to produce a theoretical DTOF. Next, a nonlinear optimization routine built using the fminsearch [67] function in MATLAB 2020a (The MathWorks Inc., Natick, Massachusetts, 2020) was used to fit each of the measured DTOFs with a theoretical DTOF by finding the set of µ a and µ ′ s values that provided the optimal fit. Fitting was conducted in several rounds to reduce error in the final µ a values. First, the routine was used to fit for µ a , µ ′ s , and an amplitude term at all wavelengths. Second, the amplitude term was fixed to its mean value, and the data was re-fit for just µ a and µ ′ s . The wavelength-dependent µ ′ s estimated from all the DTOFs was then fit with the power law [68] where λ is the wavelength, a is the scattering amplitude, b is the scattering power, and 785 nm has been used as a reference wavelength. The fitted µ ′ s were then fixed for the last round of fitting, from which we obtained the final µ a values. To reduce the effect of noise on the data analysis, the fitting range was set to 50% and 5% of the peak value on the leading and falling edges, respectively [69]. Changes in absorption due to the addition of ICG to the phantom were calculated as where µ a,IL (λ) and µ a,j (λ) are the absorption coefficients for λ at baseline (i.e., the solution with Intralipid only) and after the jth ICG addition, respectively. The fitting approach described above is the most common method used for estimating µ a of homogeneous turbid media; however, when applied to a layered medium, the presence of a superficial top layer skews the sensitivity of the derived optical properties [70]. Therefore, to determine the µ a values of the two-layer phantom, a two-layered analytical solution was used [71]. First, the optical properties of the top layer were estimated by fitting the DTOFs measured at baseline (i.e., only Intralipid in both phantom layers) using the homogeneous fitting approach described above. Next, the DTOFs measured at each ICG addition were fit with a theoretical DTOF obtained by convolving the system's IRF with a theoretical TPSF for a two-layer medium. Each layer's optical properties were set to an initial guess and a constrained nonlinear optimization routine (fminsearchbnd, MATLAB, The Mathworks Inc.) was used to change model parameters until an optimal fit was reached. As with the homogeneous phantom data, all fitting was done from 50% to 5% of each DTOF's peak value on the leading and falling edges, respectively. For both layers, amplitude and µ ′ s were fixed to the properties determined by fitting the baseline data using the homogeneous fitting approach. The µ a values obtained from the homogeneous fitting were also used as the initial guess for µ a in the phantom's top layer and constrained to vary within ±10% of their initial value. For the bottom layer, the initial guesses for µ a were estimated from a linear regression on the tail of the logarithm of each DTOF and left unconstrained during the fitting routine. Changes in the absorption of the bottom layer due to the addition of ICG were then computed by subtracting the baseline µ a values from the µ a values obtained with the two-layer fitting approach.
For comparison, the experimental DTOFs at each wavelength were temporally integrated (to obtain a dataset equivalent to CW NIRS), and the modified Beer-Lambert law was applied to calculate the change in µ a between baseline and each ICG addition. Optical pathlength was estimated for each wavelength as the product of the speed of light in the medium-assuming a refractive index of n = 1.4-and the first normalized moment of each DTOF (i.e., the mean photon time-of-flight) [72]. Figure 2 displays two spectra from the same neon light source, measured by a spectrally calibrated CW spectrometer and the TR spectrometer. Since the spectral source is a low-pressure lamp, the linewidths of its emission spectrum are very narrow and can be neglected when estimating the spectral resolution; under this assumption, the measured spectrum will be limited by the spectrometer resolution. Therefore, the theoretical spectral resolutions of the TR and CW spectrometers were calculated to be 6.0 nm and 4.7 nm, respectively, which are in agreement with the FWHM values in Table 2. The mean difference in the positions of four other peaks that were not used for calibration was 0.68±0.35 nm. Figure 3(A) shows raw and denoised versions of a DTOF measured from a single DMD pattern during an IRF acquisition. The group of raw DTOFs measured for a full set of patterns (N = 256) were used to reconstruct a "raw" IRF whose 750 nm component is shown in Fig. 3(B) (as detailed in Section 2.2). Similarly, a "denoised" IRF, whose 750 nm component is shown in Fig. 3(B) was reconstructed using denoised DTOFs. The curves in Fig. 3(B) illustrate how denoising raw data affects the quality of a reconstructed IRF. Fig. 3(C) compares the denoised 750 nm IRF from Fig. 3(B) to IRFs reconstructed from compressed acquisitions where only DTOFs for the first 12.5% (M = 32) and 37.5% (M = 96) of patterns were acquired. The strong qualitative resemblance between all three IRFs in Fig. 3(C)-even though the IRFs from compressed acquisitions were reconstructed using substantially less data than the one from the full acquisition-highlights that the majority of the data necessary for IRF reconstruction can be measured using only a small subset of patterns. In other words, the data is inherently sparse. Notably, Fig. 3(D) depicts the relationship between measurement compression and information loss for this dataset by plotting the SSIM index between IRF acquisitions with various compression levels. Note that while Figs. 3(B) and 3(C) show reconstructions of the same IRF at a single wavelength, Fig. 3(D) presents results that were averaged across IRFs at all wavelengths. Figure 4(A) shows a temporally-integrated IRF-equivalent to a spectrum acquired with a traditional CW spectrometer-reconstructed using denoised DTOFs. While DTOFs were denoised prior to IRF reconstruction, the temporally integrated spectrum still contained high-frequency oscillations. To limit the effect of these oscillations on the data, the spectrum underwent additional spectral denoising [64]. Note the consistency of the denoised spectrum: only a minor spectral shift (less than 1 nm) was detected after 9 h of continuous system operation. Fig. 4(B) compares the spectrally denoised spectrum from Fig. 4(A) to spectra reconstructed from compressed IRF acquisitions which have not been spectrally denoised. The qualitative similarity between the three curves shows that spectra reconstructed from compressed acquisitions are well-suited for capturing broad spectral features while also inherently suppressing undesirable high-frequency oscillations. This effect was confirmed by plotting the SSIM index between temporally-integrated IRF spectra-for both the original and the spectrally denoised case-from acquisitions with various compression levels ( Fig. 4(C)). Figure 5 displays a plot of the system's full IRF in 3D, along with its projection onto the spectral dimension (equivalent to a spectrum that would be obtained with a typical CW spectrometer), and components of the IRF at a few wavelengths (equivalent to multi-wavelength TR measurements). Interestingly, Fig. 5(C) shows notable differences in the propagation of different wavelengths through the spectrometer, with shorter wavelengths being detected at slightly later times. This is likely due to the angle of the DMD relative to the PMT which resulted in shorter wavelengths having a longer overall pathlength through the instrument compared to longer wavelengths. The presence of these differences illustrates the importance of rigorous characterization of the IRF at each wavelength before performing further quantitative analysis.

Intensity (a.u.)
Commercial spectrometer TR spectrometer Fig. 2. Comparison of spectra from the same neon lamp measured by a calibrated commercial CW spectrometer and the TR spectrometer.   , its projection onto the spectral dimension (i.e., collapsed across its temporal dimension) in (B), and IRF components at a few wavelengths (as in a typical multi-wavelength TR system) in (C). Note that the spectrum in (B) is equivalent to a spectrum that would be obtained with a typical CW spectrometer.

Phantom experiments
The results of the homogeneous phantom experiments are presented in Fig. 6. Fig. 6(A) shows µ a obtained by fitting the baseline measurements (i.e., phantom filled with Intralipid only) and the measurements obtained after one ICG addition. Dashed lines indicate the µ a obtained from fitting the compressed measurements (37.5% compression rate; M = 96) for both phantoms and µ a values for water [73] are also shown for reference. Note the similarity between the shape of the absorption spectrum of the Intralipid phantom and the reference water spectrum, which is not surprising given that the phantom was composed of 99.2% water by volume. To compare full and compressed acquisitions, we estimated the scattering amplitude and power of the Intralipid phantom by fitting Eq. (4) to Intralipid solution µ ′ s spectra measured from the two acquisition methods. The resulting values were in good agreement (full: a = 0.82, b = 0.69; compressed: a = 0.82, b = 0.72) as were the absorption spectra obtained using the two acquisition methods (Fig. 6(A)). This agreement is also seen in Fig. 6(B) which shows the four ICG spectra for the compressed and full acquisitions; note that these have been isolated by subtracting baseline absorption from phantom absorption following ICG addition. To assess the consistency of the shape of the recovered spectra, the mean and standard deviations of ICG curves normalized to their mean value were calculated (Fig. 6(C)). Fig. 6(C) also includes reference spectra for ICG in water and blood plasma [73]. It should be noted that the presence of Intralipid appears to alter the shape of the ICG spectrum compared to its spectrum in water: similar to the spectrum of ICG in plasma, there is a shift of the spectrum's peak towards longer wavelengths (Fig. 6(C)). Finally, using the four ICG spectra acquired with a full basis along with the concentrations from Table 1, the molar extinction coefficient for ICG in 0.8% Intralipid was calculated (Fig. 6(D)); note that the provided values follow the convention for decadic absorbance. For future reference, the numerical values of the ICG extinction coefficient are provided in Data File 1 [74] within the supplementary material.  Table 1 for concentrations). (A) shows µ a spectra for a variety of phantom states and acquisition types; a water spectrum is also included for reference. (B) shows ∆µ a between various ICG additions and baseline IL absorption (i.e., µ a of ICG). The average shapes of spectra reconstructed from both full and compressed acquisitions (37.5% compression rate) are compared to spectra of ICG absorption in water and in blood plasma (C); note that, for ease of visualization, all curves in (C) are normalized by their mean value. (D) shows a plot of the molar extinction coefficient of ICG in 0.8% Intralipid calculated from the full acquisition curves in (B); see Data File 1 [74] for underlying values of the extinction coefficient. Figure 7 shows the change in µ a (i.e., ∆µ a ) in the bottom compartment of a two-layer phantom, calculated using the two approaches described in Section 2.5. Each subplot from Fig. 7(A) -D corresponds to the spectra obtained for one of the ICG additions in Table 1. The corresponding spectra acquired in the homogeneous phantom at the same ICG concentration are shown for reference. All spectra were reconstructed from full acquisition data (N = 256), except for the yellow dashed curves reconstructed from compressed acquisitions (37.5% compression rate; M = 96). In all cases, the ∆µ a estimated using a fit to a two-layer analytical model was much closer to the ground truth-represented by the bottom layer absorption coefficient-than the ∆µ a estimated using the Beer-Lambert law, which assumes a homogeneous medium. Overall, good agreement between the reconstructions from the compressed and full acquisitions was observed at all ICG concentrations. However, Figs. 7(A) and B, show that spectra reconstructed from compressed acquisitions tend to have higher absorption near their peaks compared to spectra reconstructed from full acquisitions. The ∆µ a spectra obtained by fitting the full-acquisition data using the two-layer analytical model are displayed in Fig. 7(E), and show increased absorption as a function of the added amount of ICG. Similar to Fig. 6(C), Fig. 7(F) shows the mean and standard deviation of ∆µ a spectra from full and compressed acquisitions normalized to their mean value; these are compared to the normalized spectrum from the homogeneous phantom experiments. Compared to the spectra from the homogeneous phantom, the spectra from the bottom of the two-layer phantom have a slightly higher standard deviation, particularly between 820 -830 nm. Spectra from the homogeneous phantom and the bottom layer of the two-layer phantom exhibit a similar shift towards a typical spectrum of ICG in plasma compared to a typical spectrum of ICG in water.  Table 1 for concentrations). Spectra were recovered by using a two-layer analytical model to fit results from full (Two-layer model) and compressed acquisitions (Two-layer model -compressed; 37.5% compression rate). In addition, spectra calculated using the Beer-Lambert law as well as spectra acquired on the homogeneous phantom (i.e., the ground truth) are shown. (E) shows only the spectra recovered from full acquisition data using the two-layer analytical model. The average shapes of the spectra reconstructed from both full and compressed acquisitions are shown along with the average shape of the spectra from the homogeneous phantom results (F); note that, for ease of visualization, all curves in (F) are normalized by their mean value.

Discussion
TR spectrometers are powerful analytical tools that can measure both the temporal and spectral responses of turbid media to optical probing. This rich spectral and temporal information content can be used to get more accurate estimates of light absorber concentrations in these media. Moreover, given the broad spectral features of the main NIR light absorbers in tissue, CS can be used to substantially reduce both measurement time and dataset size in TR-NIR spectroscopy with minimal loss of critical information. To this end, this work introduces a CS-capable TR spectrometer characterized using a calibrated spectral source and demonstrates its ability to measure chromophore concentrations in homogeneous and layered tissue-mimicking phantoms.
Spectral characterization of the TR spectrometer using a neon light source showed good qualitative (Fig. 2) and quantitative ( Table 2) agreement with a calibrated CW spectrometer. Nevertheless, some subtle spectral features were not resolved by our system. For example, the spectrum measured by the CW unit clearly shows two small peaks at 735 nm and 747 nm, and a higher one at 750 nm, which are not resolvable in the spectrum acquired by the TR spectrometer. These differences are likely due to a trade-off between spectral range and resolution inherent to the spectrometer's design: the finite spectral width reflected by the thinnest lines of our chosen DMD patterns imposes a limit on the system's spectral resolution. Given the broad spectral features of interest for typical biomedical applications (e.g., absorption features of hemoglobin, ICG, and cytochrome c oxidase), the impact of these subtle spectral differences is expected to be minimal. If required for other applications, the resolution could be improved by selecting a diffraction grating with higher dispersion; however, this may limit the spectral range measured in a single acquisition.
The results presented in Figs. 3(A) and B showed that denoising raw DTOFs prior to data reconstruction substantially reduced noise in reconstructed DTOFs. Due to the statistical nature of TCSPC, the quality of raw DTOFs was expected to improve as DTOF collection time was increased. In practice, we found that collection times in the 300 -1000 ms range provided a good trade-off between measurement quality and acquisition time. Even after temporal denoising, high-frequency oscillations could be seen in temporally integrated spectra; however, the overall shape of the spectra remained stable even after 9 h of continuous operation and, if needed, these oscillations could easily be removed using spectral denoising (Fig. 4(A)). Notably, spectral denoising had an important impact on the similarity between spectra reconstructed from full acquisitions and their compressed counterparts: following spectral denoising, both the qualitative (Fig. 4(B)) and quantitative similarity ( Fig. 4(C)) between a full acquisition and a compressed acquisition that used only 12.5% of the total raw measurements was negligible. In contrast, comparison of a non-denoised (i.e., raw) spectrum to its counterparts acquired with compressed acquisitions resulted in a noticeable SSIM index drop-off ( Fig. 4(C)). This clearly demonstrated that broad spectral features-which remain after spectral denoising-are highly conserved when applying a compressed sensing approach, which is one of the key advantages offered by the configuration in Fig. 1. Overall, a good qualitative and quantitative agreement between reconstructions from full and compressed datasets across both the temporal (Figs. 3(C) and D) and spectral dimensions (Figs. 4(B) and C) was found.
Consistency between spectra acquired using compressed sensing and basis scan approaches was also confirmed by the homogeneous phantom experiments (Fig. 6(A)). As expected, the absorption spectrum of the phantom filled with 0.8% Intralipid showed similarities to a reference water spectrum. Each addition of ICG increased the phantom's µ a across all measured wavelengths; similar µ a values (Fig. 6(B)) and spectral shape ( Fig. 6(C)) were obtained from both compressed and full reconstructions. Further, the standard deviation of the recovered ICG spectral shape was low, which highlighted the spectrometer's ability to recover a consistent spectral shape for a variety of chromophore concentrations. However, the differences between the ICG absorption spectra reconstructed from full and compressed acquisitions increased at higher concentrations ( Fig. 6(B)), which was likely due to lower measurement signal-to-noise ratio (SNR) at higher absorptions. As described in Section 2.5, we note that wavelength-dependent effects of the device's IRF on the measured spectra were taken into account through convolution during the fitting procedure and that the amplitude term used in the fitting to account for attenuation/amplification factors (e.g., NDF, fiber coupling) was not expected to have a noticeable wavelength dependence.
In general, the absorption spectrum of ICG is known to depend both on its concentration and its solvent [39,73]. Yuan et al. previously measured absorption coefficients of 0 to 2.0 µM ICG in 0.6% Intralipid at 660, 780, and 830 nm [40]. One of their key findings was that the absorption of ICG in Intralipid at 830 nm was higher than its absorption at 660 nm. This is different from the absorption of ICG in water and is consistent with shifting of the ICG spectral shape towards longer wavelengths as observed herein (Fig. 6(C)) and in previous studies which used higher concentrations of ICG [39]. Considering that similar spectral shifting has been noted in ICG-plasma solutions [73] and is generally attributed to ICG adsorption to plasma proteins, Du et al. previously inferred that this effect may be related to the binding of dye to Intralipid particles [39]. Overall, the results of the homogeneous phantom experiment demonstrate the TR spectrometer's ability to recover absolute µ a values in a turbid, homogeneous medium across a range of wavelengths using both basis scanning and compressed sensing approaches. For future reference, we report the molar extinction coefficient of ICG in 0.8% Intralipid (Fig. 6(D)) calculated from the 4 full acquisitions obtained at the concentrations in Table 1.
One of the most important advantages of TR spectroscopy is the ability to remove the contribution of superficial layers when investigating chromophore absorption in deeper layers. The two-layer phantom results presented in Figs. 7(A) -D demonstrate the advantage of the TR approach by contrasting relatively small bottom layer absorption changes obtained using the Beer-Lambert Law with the larger changes obtained using a fit to a two-layer TR analytical model. Utilizing the more sophisticated model allows recovery of µ a values much closer to the ground truth acquired directly from the bottom layer. Nevertheless, this method still incorrectly estimates µ a with an average error of 24±5% across all wavelengths and ICG additions. Interestingly, though spectra recovered from compressed acquisitions generally resemble those from full acquisitions, they more closely matched the ground truth µ a values after the first few ICG additions; this resulted in a slightly lower overall µ a recovery error (21±7%). Similar to the homogeneous phantom results, each ICG addition increased the phantom's µ a across all measured wavelengths (Fig. 7(D)) with a similar spectral shape for both the compressed and full reconstructions (Fig. 7(E)). Nevertheless, the ICG spectrum recovered for the final ICG addition had a slightly different shape than the previous spectra, as reflected in its lower than expected absorption in the 730 -770 nm range. Further investigation revealed that this discrepancy was due to crosstalk between the top and bottom layer absorptions recovered by the two-layer fitting approach. Specifically, a notably higher top layer absorption was recovered between 730 -770 nm for this ICG concentration, compared to the three others. Since no changes to the top layer were made throughout the experiment, this discrepancy is likely due to the lower SNR at this higher ICG absorption, which resulted in lower data quality and, subsequently, less consistent fitting results. As in the homogeneous phantom experiments, the recovered spectra were shifted towards longer wavelengths compared to a typical ICG absorption spectrum in water. However, in contrast to the shift seen in the homogeneous phantom, the red shift of the absorption spectrum of the bottom layer was smaller and-in the 710 -750 nm range-more closely matched a shifted ICG-water absorption spectrum.
The two-layer analysis presented in this work has previously been used to separate cerebral hemoglobin concentrations and oxygen saturations from superficial contributions in both simulations [76] and subjects [77]. One approach to estimate top layer properties is the use of an additional measurement at a shorter S-D distance [77]. However, we avoided this approach since the TR system used in the study was not designed for multi source-detector measurements, and changing probe distances after setting up the phantom measurements could lead to larger errors due to changes in coupling between the probes and the phantom. To further model the error in the recovered bottom layer µ a values, we conducted TR simulations in a two-layer slab using NIRFAST [78,79]. The slab's top layer (12 mm) and bottom layer (70 mm) optical properties were set to the values recovered from baseline and post-ICG addition in the homogeneous phantom, respectively. Next, the percent error in recovered µ a values was calculated for three scenarios of top layer thickness input into the analysis: correctly estimated, underestimated, and overestimated (Table 3). In addition, for the case when the top layer thickness input was correct, changes in the percent error due to the addition of random noise to the raw data (Table 3) was also calculated.  Table 3 illustrates that the quality of the data, as well as errors in assumptions inherent to the analysis itself, may have substantially contributed to the discrepancy in the recovered µ a values of Fig. 7. Nevertheless, the advantage of using TR data instead of typical CW data to recover absorption changes in a turbid layered medium remain clear.
Compressive TR spectroscopy exploits the signal sparsity typically present in spectra of chromophores which play key roles in biomedical applications. As shown by the similarity of spectra acquired using basis scanning and compressed sensing approaches ( Fig. 6 and Fig. 7), this technique can be used to reduce measurement time and dataset size with minimal loss of information. In addition, the use of the DMD to reflect diffracted light into the system's detector allows for a rapid and flexible way to selectively attenuate specific wavelengths by manipulating what fraction of mirrors are "ON" in the y dimension (Fig. 1); this can be desirable when attempting to resolve changes in weaker spectral lines surrounded by intense ones [80]. It is important to note that, when used with a coherent light source, the periodic spacing of the DMD's micromirrors causes the DMD to act as a diffractive element [81]. Since diffraction is wavelength-dependent, the diffraction orders produced by the spectral lines reflected from the DMD will no longer maintain the same relative spacing as in the incident spectrum. However, the spectrometer operates under a single-pixel paradigm, which means that the DMD is responsible for spatially sampling the spectrum and establishing the spectral-to-spatial mapping. Given that the aforementioned disruption to wavelength spacing occurs after this mapping is established it should not affect the system's spectral accuracy. Further, any influence this effect could have on the light intensity can be quantified by acquiring an IRF (Fig. 5). While the spectrometer presented herein only investigated a spectral range of approximately 120 nm, this limitation was mainly due to suboptimal coupling between the DMD and the photosensor module (see Section 2.1). If needed, the device could be modified to accommodate a wider spectral range by exchanging the detector for one with a larger active area, selecting a wider DMD for the spectrometer, and optimizing the coupling between these two components. Further, while the system had slightly coarser resolution than a CW spectrometer (Table 2), resolution can be modified using finer or coarser DMD patterns during acquisitions. Future work will focus on testing the system using a more comprehensive array of chromophores and configurations that enable greater spectral ranges and resolutions. In addition, the system will be tested using more sensitive detectors, which are likely to improve the system's SNR and enable further reduction of acquisition time.

Conclusion
This work presents the development of a compressive sensing TR spectroscopy system for accurate quantification of light absorbers in turbid media. Notably, the system's ability to quantify the absorption coefficient of ICG in both homogeneous and layered tissue-mimicking Intralipid phantoms was demonstrated. Further, the phantom results revealed that compressed sensing matched or outperformed the traditional basis scanning method (i.e., uncompressed measurements) in all cases. Using data from the phantom experiments, the extinction coefficient spectrum for ICG in 0.8% Intralipid over a broad spectral range (710 -830 nm) was determined. Future work will focus on exploring a wider array of system configurations to accommodate wider spectral ranges, improve the system's SNR, and further reduce acquisition times before testing the device for neuromonitoring.