1. Introduction

Though long promised, coherent anti-Stokes Raman scattering (CARS) spectroscopic microscopy (microspectroscopy) has only recently demonstrated broadband hyperspectral biological imaging at acquisition rates far in excess of what traditional Raman microspectroscopy can provide [1–6]. With an imaging speed as fast as 50000 spectra per second [7], a new fundamental challenge arises: high throughput extraction of Raman vibrational information from the raw CARS spectra.

CARS spectra are quintessentially a coherent mixture of photons generated through vibrationally resonant (Raman) and nonresonant (electronic) processes. The electronic contribution is typically referred to as the "nonresonant background" (NRB) and is the root cause of CARS spectral distortion. Thus, a significant effort was made in the early years of CARS microscopy development to reduce the NRB via optical means [8–11]. The NRB, however, behaves as a stable homodyne amplifier for the Raman-generated signal; thus, reducing the NRB also reduces the Raman components of the signal. So important is the NRB's role in signal amplification [12], that without it CARS may show little to no benefit over spontaneous Raman spectroscopy for biological imaging [13].

Unlike additive fluorescent background signals in Raman spectroscopy, the NRB is coherent with the co-generated Raman-resonant CARS components. Beneficially, the NRB may amplify weak signals above the noise floor. Due to this coherent mixing that induces distortions in the spectral shapes, however, the NRB cannot be simply subtracted from the CARS spectra. There is, though, a fixed phase relationship between the Raman- and NRB-components. This inherent property led to the realization that computational methods could be used to extract the Raman portion of the CARS spectra using so-called "phase-retrieval methods": the Kramers–Kronig relation (KK) [14] or the maximum entropy method [15]. These early works assumed that the NRB was either known a priori or the NRB of a surrogate material (e.g. coverslip glass, water, salt [16]) was appropriate. Later, it was demonstrated that using surrogate materials for NRB approximations led to amplitude and phase errors that were linked analytically [17]. These errors could be corrected using "phase-error" correction (PEC) and "scale-error" correction (SEC) methods [17], which also reveals the relationship between the actual NRB and the surrogate. Importantly, this relationship demonstrated that CARS is unique among imaging techniques: it is inherently self-referencing. The spectral ratio of the Raman component to the actual NRB is an inherent property of a molecular system; thus, this ratio is maintained even in the case of sample scatter or absorption – just the signal-to-noise ratio (SNR) is affected. This enables one-to-one comparison of spectra between samples and even different CARS architectures (with different laser systems and wavelengths) [17]. Other coherent Raman methods, most notably stimulated Raman scattering (SRS) microscopy/spectroscopy [18], do not co-generate an NRB and do not have this internal referencing ability. Thus, SRS spectra are undistorted and useful for chemical identification, but the spectral amplitudes are not necessarily directly comparable with other results, potentially challenging quantitative analysis.

To extract robust, quantitative Raman component spectral data from the CARS spectra and to support the rapidly increasing data rates and volumes, we have developed a series of new methods collectively referred to as "factorized Kramers–Kronig and error correction" (fKK-EC). The new, unique principle of fKK-EC is that raw CARS spectral data can be factorized/decomposed into a small set of basis vectors on which the necessary processing steps will actually be performed. In this work, we use singular value decomposition (SVD) for its robust, accurate decomposition of matrices, although it is possible to use others as well. Previously, SVD has been used for denoising [1,17,19], but the remainder of operations were performed on the individual spectra. Additionally, matrix factorization, such as non-negative matrix factorization (NMF) / multivariate curve resolution (MCR) have been applied to post-processed data for analysis [19,20].

The fKK-EC is composed of three parts that will be described theoretically in more detail below: phase retrieval via a factorized KK relation (fKK), factorized PEC (fPEC), and factorized SEC (fSEC). These three parts operate on the basis vectors; thus, the image data is not reconstituted between each step. This limited operation on a small number of basis vectors is economical in terms of speed and memory usage without losing the spectral information, compared with the previous methods. Furthermore, basis vector sets can be re-used on new data; thus, the fKK-EC method can be used like a machine learning method, ML:fKK-EC, for short. In this paradigm, the full fKK-EC is performed (“trained") on a portion of data (e.g., the first image), and subsequent images are able to be processed (in full) via matrix multiplication. This factorized method enables new data to be processed on-the-fly in real-time during acquisition: denoised, phase-retrieved, and phase- and scale-error corrected. Like all ML methods, this process does require that the training data reflect what will be contained in upcoming data – though this is readily testable as will be discussed.

2. Theory

2.1 Background: conventional post-processing for a single CARS spectrum

CARS is a third-order nonlinear scattering phenomenon in which two photons ("pump" and "Stokes") excite a Raman vibrational mode from which a third photon ("probe") inelastically scatters [2]. Furthermore, this process does not happen in isolation and other nonlinear processes, such as degenerate four-wave mixing, may occur, leading to the generation of a so-called nonresonant background (NRB). So ubiquitous is the NRB that theoretical treatments of the CARS mechanism automatically incorporate the NRB, and the term "CARS signal" implies a coincident NRB. Thus, in this manner the CARS signal, $I_{\mathrm {CARS}}$, may be described as [1]:

The overarching goal of phase retrieval methods is to extract Im$\{\chi ^{(3)}(\omega )\}$ from $I_{\mathrm {CARS}}(\omega )$, which is the equivalent material property probed by traditional Raman spectroscopy [21]. If $\tilde {C}_{\mathrm {st}}(\omega )$ and $I_{\mathrm {NRB}}(\omega )$ are quantitatively measurable/known, this goal would be achievable [14]. However, this has not thus far been demonstrated. A more capable solution that also leads to the aforementioned self-referencing of CARS, is to calculate $K_{\mathrm {CARS}}(\omega ) \triangleq \chi ^{(3)}(\omega )/\chi _{\mathrm {nr}}(\omega )$. Using the KK formalism and assuming, for the moment, that the $I_{\mathrm {NRB}}(\omega )$ of the sample itself is measurable [17]:

In summary, using the KK relation, PEC, and SEC, one can calculate $K_{\mathrm {CARS}}(\omega )$ from $K(\omega )$ as:

 

Fig. 1. (a) Workflow for the conventional denoising, KK, PEC, and SEC, where $m$ is the total number of spatial pixels (“flattened") and $n$ is the number of frequency channels. (b) Workflow for the fKK-EC where the processing steps are performed on basis vectors rather than the underlying spectra. (c) Workflow for the ML:fKK-EC in which only the training data is processed via the fKK-EC and regression is used to transform new data.

Download Full Size | PDF

2.2 SVD factorization, denoising, and fKK

The proposed fKK-EC enables high-throughput and real-time Raman signal extraction from spectroscopic CARS data via factorization, which dramatically reduces the number of vectors for which each processing step is applied. For example, rather than independently applying to one million spectra in a one-megapixel image, the processing may be applied to 100 basis vectors. A flow chart that describes the fKK-EC workflow is shown in Fig. 1(b).

The first step in this process is factorization of the input data. In this work, SVD decomposes a matrix $\mathbf {A}$ into three matrices as $\mathbf {A} = \mathbf {U}\mathbf {S}\mathbf {V}^H$. The $H$-superscript indicates the Hermitian transpose; $\mathbf {U}$ and $\mathbf {V}$ are unitary matrices whose columns are the left- and right singular vectors, respectively; and $\mathbf {S}$ is a diagonal matrix whose entries are known as singular values (we will denote the vector containing just the singular values as $\mathbf {s}$). In this work, we explicitly assume that the dataset is oriented so that row-number ($m$) corresponds to spatial components (see Fig. 1(a)) and column-number ($n$) to frequency. Thus $\mathbf {U}$ is composed of spatial basis vectors while $\mathbf {V}$, spectral basis vectors. Further, $\mathbf {A}$ is real; thus, the Hermitian transpose ($H$) is a transpose ($T$) as will be indicated in the remaining derivations. The SVD [1,17,19,24,25] is widely used for denoising by removing noise-dominant singular vectors that [ideally] only contribute to noise. This is accomplished by either setting singular values corresponding to noise-related singular vectors to 0, or equivalently creating new $\mathbf {U}$, $\mathbf {S}$, and $\mathbf {V}$ matrices that exclude the non-desired singular values and vectors, which leads to reduced data volumes. We have implemented the latter in our simulations and experiments. Note that in the remaining derivations we do not explicitly denote whether $\mathbf {U}$, $\mathbf {S}$, or $\mathbf {V}$ were altered for denoising; though, all derivations remain valid. Also note, that in the conventional workflow, SVD is an optional denoising method (see Fig. 1(a)), whereas in the fKK-EC workflow, SVD is the necessary factorization step (See Fig. 1(b)), which can also provide denoising.

For the fKK, one would conceptually apply the KK relation to the spectral basis vectors in $\mathbf {V}$. However, this is not appropriate because of the log-function in the KK (Eq. (2)) and the orthonormal nature of SVD singular vectors (positive- and negative-values). Rather we apply the SVD to $\frac {1}{2}\ln {[I^{[m]}_{\mathrm {CARS}}(\omega )/I_{\mathrm {ref}}(\omega )]} = \mathbf {a}_m(\omega )$, where the $m$-superscript denotes the m$^{th}$ spectrum, which leads to the $\mathbf {a}_m$-vector. For the following derivation, we assume that we have $M$ spectra, and the $\omega$ vector is $N$-frequency increments long. Thus, $\mathbf {A}$ may be written as:

As an addendum to this derivation, we will discuss considerations under the case of mixed Poisson-Gaussian noise (heteroscedastic noise generally). In previous work [17], denoising was improved via the use of an Anscombe transformation prior to SVD. As Poisson noise is not additive [27], SVD is often impaired in separating signal and noise. The Anscombe transform aims to convert a signal with mixed noise into a signal with unit variance. Though advantageous, this nonlinear transform is not compatible with the current fKK derivation. Thus, to improve denoising, there are 2 options: (1) denoise before the fKK using the Anscombe transformation and SVD (then reconstruction), or (2) apply a scaling term $f(\omega )$ to $A(\omega )$, which is the same for each spectrum. In simulations and experiments below, we apply the latter. The scaling term we selected was inspired by the purpose of the Anscombe transformation: normalizing variance. Suppose we have a homogeneous sample and have recorded numerous spectra containing both additive white gaussian noise and shot (Poisson) noise, the standard deviation ($\sigma _{\mathbf {A}}(\omega )$) of the previously defined $\mathbf {A}$ may be approximated as [28]:

In the remainder of this manuscript, we will include $f(\omega )$ in the derivations; though, this factor can be set to one in the case of pure additive white Gaussian noise. Mathematical notation note: we are explicitly writing $f(\omega )$ to emphasize that it is a single spectrum, and when it divides a matrix, it is applied along the spectral axis (e.g., each row of $\mathbf {A}$ or $\mathbf {V}^T$).

2.3 Factorized PEC (fPEC)

PEC is the process of finding the phase error caused by using a surrogate reference material as an approximation for the sample NRB. In the factorized context:

The action of the combined fKK and fPEC without fSEC can be described as:

2.4 Factorized SEC (fSEC)

In the conventional form of the SEC, the PEC-corrected spectra are divided by the mean of the real part as described in Eq. (5). An alternative and equivalent approach is to calculate the mean of the natural log of the magnitude of the PEC-corrected spectra:

The left-hand expression in Eq. (19) for the dataset is equivalent to the magnitude of the term inside the exponential function in Eq. (18) as:

For the fSEC, we want to avoid calculating the mean for each spectrum and to operate on the PEC-corrected right singular vector. Thus, we will incorporate an fSEC correction matrix $\mathbf {V}^T_{\mathrm {SEC}}$ into the previous expression:

2.5 Reconstruction and the full fKK-EC

Using the previous descriptions of the fKK, fPEC, and fSEC, we can assemble the full fKK-EC workflow and reconstruct an approximate $\mathbf {K}_{\mathrm {CARS}}$ (akin to Eq. (5) for the conventional implementation). Applying Eqs. (8), 18, and 22:

2.6 Machine learning (ML) paradigm ML:fKK-EC

As previously described, the fKK-EC methods enable high-throughput analysis at significantly higher rates than the coventional workflow. Another significant benefit of the fKK-EC methods is that they can be trained as a machine learning (ML) model, i.e., the fKK, fPEC, and fSEC are fully applied to a sub-set of data, and new data is simply projected onto the derived basis vectors (as schematically described in Fig. 1(c)). That is to say that new data can be transformed into denoised-Raman-retrieved (fKK, fPEC, fSEC) without explicitly applying these methods, but rather with simple matrix multiplication. We will call this workflow "ML:fKK-EC".

Hypothetically, we are going to collect many images of a sample. We will apply the full fKK-EC method to the first (or first few) images (i.e., "training"). This provides us with: $f(\omega )$, $\mathbf {U}$, $\mathbf {S}$, $\mathbf {V}$, $\mathbf {\Phi }_{\mathrm {PEC}}$, and $\mathbf {V}^T_{\mathrm {SEC}}$. One assumes that upcoming images will comprise the same or similar chemical content but in differing concentrations and mixture profiles. In the ML:fKK-EC method, we will not re-derive the SVD, but rather regress a new left singular vector matrix $\mathbf {U}_{\mathrm {new}}$ (which describes the spatial mixtures of $\mathbf {S}\mathbf {V}^T$). From Eq. (7) for the new data, $\mathbf {A}_{\mathrm {new}}$ that incorporates $f(\omega )$ as well, and solving for $\mathbf {U}_{\mathrm {new}}$ applying ridge regression:

Of importance is the selection of the training data. As in all “supervised" machine learning methods, the training set has to include a certain variety and depth of training data [30]. In the ML:fKK-EC, the training set needs to contain enough variety in chemical components so that the spectral basis set ($\mathbf {V}$) can reconstruct new spectra in future images. Various strategies for this could exist: coarsely image a large area and use that for training, apply expert knowledge to select a few small images for training, use the first image for training and retrain later — it is sample and application specific. Fortunately, it is facile to test whether the trained basis set is adequate for new data by comparing $\mathbf {A}_{\mathrm {new}}$ as measured and as reconstructed applying Eq. (25) to Eq. (24). The precise comparison could be performed programmatically such as through calculating the residual sum-of-squares (RSS) or by manual inspection. If the original basis set is inadequate, one could simply append the unsupported spectrum (or spectra) to the original training set and recalculate the SVD, i.e., it does not necessitate retraining on an additional huge swath of data. Other strategies and on-going work in this context will be provided in the Discussion section.

The ML:fKK-EC method, as will be demonstrated in simulation and experiment below, is extremely fast. Firstly, the time-consumption of the individual steps is limited to a training dataset that is much smaller than the full dataset. Secondly, new data does not need to be subjected to the fKK, fPEC, or fSEC, but rather is converted through a series of matrix multiplications: solving for Eq. (25) and applying to Eq. (23), where all the other matrices were calculated during training. For example, on the broadband CARS (BCARS) system used to collect data for this paper, spectra require $\approx$ 5 ms to record, but applying the ML:fKK-EC to a new spectrum requires 10's of microseconds; thus, it can be applied to new data as it is acquired, as opposed to after all data is acquired. This advancement in CARS microscopy affords many new opportunities not previously available, such as on-the-fly evaluation of imaging quality and rapid identification of regions-of-interest and chemical constituents.

3. Materials and methods

3.1 Broadband CARS (BCARS) imaging platform and software

Images were collected on an in-house-developed BCARS microscope that is described in detail elsewhere [1]. In brief, the picosecond probe laser ($\approx$ 3.2 ps, 771 nm, 40 MHz repetition rate; Toptica, FemtoFiber pro NIR) and femtosecond supercontinuum ($\approx$ 16 fs, 850 nm to 1350 nm, 40 MHz repetition rate; Toptica, FemtoFiber pro UCP) were 13 mW and 7.1 mW on-sample, respectively, after the water-immersion microscope objective (Olympus, UPlanSApo 60XW IR). The anti-Stokes photons were collected (Olympus, LUCPlanFL N 40X), filtered via short-pass filters (Semrock) and focused onto a charge-coupled device (CCD)-equipped (Andor, DU970N-FI) spectrograph (Acton, SpectroPro2300i). The CCD integration time was set to 3.5 ms, which corresponds to $\approx$ 5 ms per pixel owing to data transfer time, stage movement, and CCD cleaning time. The sample is mounted onto a two-stage system: one for high-speed, high-resolution XYZ-imaging (Physik Instrumente, P-545) and a long-travel XY-stage (Physik Instrumente, M-545.2PO). The maximum size of a single XY-image is 200 $\mu$m x 200 $\mu$m (limited by the high resolution stage); thus, movement with the long-travel stage and stitching is necessary for larger images.

The BCARS system was controlled by custom LabView software written in-house. Data files were processed in Python using NumPy, SciPy, scikit-learn, and the open-source CRIkit2 software package for Python (https://github.com/CCampJr/CRIkit2). Processing was performed on a Dell Precision 7730 laptop with a 6-core i7-8850H processor at 2.6 GHz and 64 GB of RAM.

3.2 Chicken tissue preparation

Chicken legs were procured from a local grocer. Hyaline cartilage tissue was harvested from the knee joint above the tibia using a scalpel. The resected tissue varied in thickness from approximately 20 $\mu$m to 40 $\mu$m, as measured by BCARS imaging ("XZ" images).

3.3 Simulation software

The simulations were written in Python and performed from within a Jupyter Notebook. The NumPy, SciPy, scikit-learn, Pandas, Seaborn, and CRIKit2 software packages for Python were used for processing and visualization. Simulation software will be furnished upon request and will be available in a forthcoming open-source software package for Python. The simulations were performed on the same laptop as the image processing described above.

4. Results

Below we present simulations and experiments to demonstrate the enhanced performance (throughput) of fKK-EC and the comparability of its results with the conventional workflow. Additionally, within the experimental results, we demonstrate the application and results from the ML:fKK-EC.

4.1 Simulation

We simulated a noiseless 3-chemical mixture with the concentration map shown in Fig. 2(a) and a ternary plot of concentrations shown in Fig. 2(b). Chemical 1, 2, and 3 are displayed in red, green, and blue, respectively. The base dataset is 74 pixels x 246 pixels (18204 total spectra). To analyze the fKK-EC performance versus number of spectra, this dataset is side-scaled by a factor of 0.5, 1, 2, 3, and 4; for a total of 4551, 18204, 72816, 163836, and 291264 spectra, respectively. Synthetic Raman spectra were generated using a summation of complex Lorentzian functions with number of peaks, amplitude, central frequency, and width being selected stochastically. Further, the real-valued $\chi _{\mathrm {nr}}(\omega )$’s were quadratic polynomials with randomly generated non-negative coefficients, and $\chi _{\mathrm {ref}}(\omega )$ from a linear polynomial. This approach was not chosen because of its physical realism, but rather to challenge the method — especially the detrending algorithm. The random number generator seed was fixed across experiments so that the same random spectra were generated. The simulated CARS spectra (and NRB) are shown in Fig. 2(c). The chemical spectra contain 22, 25, and 10 peaks, respectively. The spectral range of simulation was $-$500 cm$^{-1}$ to 2500 cm$^{-1}$ sampled 810 times; though, Raman peaks could only be assigned between 500 cm$^{-1}$ to 1700 cm$^{-1}$. The stimulation profile $\tilde {C}_{\mathrm {st}}(\omega )$ in Eq. (1) was set to a constant for simplicity.

 

Fig. 2. (a) Concentration map of simulated dataset composed of three chemical constituents, colored as red (Chem 1), green (Chem 2), and blue (Chem 3). (b) Ternary plot showing the concentration of the simulation. (c) CARS (solid lines) and NRB (dashed) spectra of the three pure constituents. The spectra of a reference surrogate is shown in black. (d) Processing speed enhancement of fKK-EC and ML:fKK-EC with respect to conventional processing. Each trace shows the mean enhancement over three runs with the shading showing $\pm$1 standard deviation. (e) Fraction of computing time of each step. Note: only fKK-EC and ML:fKK-EC have a reconstruction (Reconst.) step. Also, only the ML:fKK-EC has a regression (Regress.) step. (f) Mean RSS showing the factorized methods show relatively similar, if not improved, RSS values from the conventional workflow. (g)–(i) Comparison of single-pixel spectra processed using the conventional methods and the fKK-EC. (j)–(l) Comparison of single-pixel spectra processed using the fKK-EC and the ML:fKK-EC.

Download Full Size | PDF

Figure 2(d) shows the speed enhancement of the factorized methods relative to the conventional workflow. For all methods, the number of kept singular values/vectors was determined by the singular values larger than ($\max {\mathbf {A}} \times \max (M,N) \times \epsilon$), where $M$ and $N$ are the row and column dimensions of the SVD-input matrix $\mathbf {A}$, and $\epsilon$ is the "machine epsilon" for the given data type. This is the same cutoff used to estimate rank in NumPy and MATLAB software. For comparison, the time per spectrum for the conventional workflow was approximately: $\leq$100 $\mu$s for SVD and selecting basis vectors, $\leq$140 $\mu$s for the KK, $\leq$3.2 ms for the PEC, and $\leq$140 $\mu$s for the SEC; for a total of $\leq$3.6 ms / spectrum. In each conventional-method simulation run, 6 basis vectors were kept per the previously described cutoff threshold. In all factorized-method simulation runs 35 to 50 singular vectors were kept, depending on the image size. For the fKK-EC, the enhancement was $\geq$40 for all but the smallest dataset. For the 291264 spectra simulation, for example, the total time was $<$25 seconds for all 3 replicate simulations (86 $\mu$s / spectrum). The most significant difference is the time to perform phase retrieval, with the conventional KK requiring $\approx$ 40 s and the new fKK $\approx$ 4.3 ms — an over 9000$\times$ improvement. The fPEC was over 1250$\times$ faster than the PEC, and the fSEC was over 3150$\times$ faster than the SEC. For the factorized workflow, the reconstruction step only added 3.3 s. Figure 2(e) gives a graphical representation of the fraction of total computational time for each method. Of course it should be noted that for the ML:fKK-EC, the training fraction will reduce as more non-training data is processed.

We also compared the spectra obtained by the fKK-EC method with that of the conventional method. To that end, we calculated the residual sum-of-squares (RSS) between the extracted Raman-to-NRB ratio spectra ($K_{\mathrm {CARS}}$ in Eq. (5) or Eq. (23)) and the known Im$\{\chi /\chi _{\mathrm {nr}}\}$ at each pixel. The mean RSS, $\langle$RSS$\rangle$, is shown in Fig. 2(f). For reference, the RSS if $K_{\mathrm {CARS}}(\omega )~=~0$ (“Null RSS") is also shown. One can see that the fKK-EC and conventional workflow return similar results, with the fKK-EC being slightly better (lower). Whether this is intrinsic or due to imperfect hyperparameter tunings for each processing step (e.g., ALS parameters) will be investigated in the future as the current goal was to demonstrate approximately equivalent results. Figures 2(g)–2(i) compare the spectra retrieved by the conventional method and the new fKK-EC (versus the ideal) at the pixels with the maximum concentration of each simulated chemical species. In each instance, the fKK-EC spectrum returns a result closer to the ideal than the conventional method. It was determined that all errors were due to the phase error-correcting steps: the ALS could closely but not perfectly retrieve the phase error. Under a separate simulation using constant-valued NRB’s, the ideal, conventional workflow, and fKK-EC all agreed ($\langle$RSS$\rangle$ $<$10$^{-14}$).

Next, we performed the same comparisons using the ML:fKK-EC implementation. The training portion of the dataset is identified in Fig. 2(a), and it is performed independently for each dataset size. Figure 2(d) shows the speed enhancement of ML:fKK-EC versus the conventional workflow, both including ("+train") and excluding ("$-$train") the time used for the training portion. Thus, for a trained ML:fKK-EC system, we calculate an $\approx$ 150$\times$ speedup, which was $<$30 $\mu$s per spectrum for all dataset sizes. Thus, this could be performed in real-time as the data is acquired. Figure 2(e) shows the computational fraction of each step, and Fig. 2(f) shows that the machine learning implementation provides equivalent RSS to the non-ML fKK-EC method. Finally, Figs. 2(j)–2(l) compare the retrieved spectra from the ML:fKK-EC and non-ML version: the results are indistinguishable.

4.2 Experimental: chicken cartilage tissue imaging

Next, we analyzed a stitched series of BCARS images (9) of hyaline cartilage excised from chicken knee tissue. The individual original images are 300 pixels x 300 pixels, with $\approx$ 3 % overlap (per side) with neighboring images. The stitched image is 846 pixels by 831 pixels (703026 pixels total). Figure 3(a) shows a pseudocolor image from the fKK-EC process, colorizing DNA, collagen, and lipids. The DNA was highlighted utilizing the peak at 720 cm$^{-1}$. To maximize contrast between DNA and other chemical components, we used the side of this peak 713 cm$^{-1}$, subtracting a linear interpolated baseline between (691 to 738) cm$^{-1}$. Tentatively, we assign this peak to the nucleotide adenine [31]. We did not see a strong peak at 785 cm$^{-1}$, which corresponds, in part, to phosphodiester stretch of the DNA backbone; thus, we hypothesize, that DNA-nucleases may have degraded the DNA as this is not fresh chicken tissue, but rather grocery store procured. The collagen was highlighted by 855 cm$^{-1}$ (proline ring C-C-stretch [32]) peak relative to the trough at 900 cm$^{-1}$. Lipids were highlighted using the intensity at 2837 cm$^{-1}$ (CH$_2$-symmetric stretch [33]).

 

Fig. 3. (a) Pseudocolor image derived from fKK-EC processed CARS image, highlighting DNA (yellow), collagen (cyan), and lipids (red). Inset shows zoomed-in region. Scale bar is 100 $\mu$m. (b) Single pixel spectra for locations identified by arrows in (a). (c) Comparison of spectrum processing time between conventional and fKK-EC workflow.

Download Full Size | PDF

Spectra retrieved using the conventional method and the fKK-EC are shown in Fig. 3(b) with the locations identified in Fig. 3(a). The spectra are qualitatively the same. Differences were identified as a result of the different response of the SVD to raw BCARS spectra versus that of the log-CARS-to-Reference dataset. Retrieving such similarly denoised and processed spectra was $\approx$ 70$\times$ faster using the fKK-EC methods (average of 3 repeats $\pm$ 1 standard deviation: conventional method $\approx$ 4973 s $\pm$ 26 s total [$\approx$ 7.0 ms / spectrum]; fKK-EC $\approx$ 70 s $\pm$ 3.0 s total [$\approx$ 99 $\mu$s / spectrum]). It should be noted that for the conventional processing, computer memory limitations precluded the processing of the entire image at once; thus, the speed was estimated by performing the KK, PEC, and SEC on 10000 spectra of the image and scaling up the time. The SVD/denoising was performed on the whole image. The fKK-EC and ML:fKK-EC were performed on the entire image.

Next we processed the same image using the ML:fKK-EC, using 1 of the 9 images as the training image (see Fig. 4(a)). The training image contained 78114 spectra. This image was selected for training as even under brightfield observation it appeared to contain collagen, lipid droplets, and nuclei: three components expected from our knowledge of hyaline cartilage. Other images would have served as well. Again the retrieved spectra, see Fig. 4(b), show qualitatively similar results to the conventional workflow with slight noise and baseline differences. Excluding the training time ($<$ 10 s), this method was approximately 150$\times$ faster than the conventional workflow, requiring $\approx$ 47 $\mu$s / spectrum to process the entire image. Though these images could have been analyzed in real-time, they were processed after acquisition.

 

Fig. 4. Pseudocolor image derived from ML:fKK-EC processed CARS image, highlighting DNA (yellow), collagen (cyan), and lipids (red). The dashed-white box indicates the sub-image used for training. Scale bar is 100 $\mu$m. Arrows identify single-pixels used for spectral comparison in (b) between conventional and ML:fKK-EC workflow, which shows close agreement. (b) Single pixel spectra for locations identified by arrows in (a). (c) Comparison of spectrum processing time between conventional and ML:fKK-EC workflow.

Download Full Size | PDF

5. Discussion and conclusion

Traditionally, the acquisition of CARS spectra was slow, requiring at least tens of milliseconds per spectrum, and most CARS hyperspectral imagery was for a small data size (up to 256 pixels x 256 pixels). Therefore, the speed of individual spectrum-based processing methods was sufficient for the old type of CARS hyperspectral imaging. However, now that the advanced CARS imaging can collect much larger images at a much faster speed, new hyperspectral image processing methods are needed. An additional complication, owing to the inherent distortion of raw CARS spectra, is that the quality and results of an imaging experiment cannot be ascertained until after processing. This, of course, has been a big incentive to use alternative modalities, such as SRS. But as previously described, those alternative modalities do not have the self-referencing ability of CARS, which may be a boon for quantitative analysis. Thus, the aim of this work is the development of high throughput, robust self-referenced Raman signal extraction from CARS spectra with real-time capability.

Though this work demonstrates that the factorization approaches are supremely efficient and capable of being used in a machine learning paradigm, there are still many improvements possible and areas of inquiry for these methods. From a physics/chemistry perspective, we are actively modeling and investigating the nature of the NRB and differences between NRBs of different materials. Further, we are examining the degree to which the real-valued $\chi _{\mathrm {nr}}$ assumption is valid in light of multiphoton resonances often found in biomolecules. This information would not only improve quantitative analysis, but as related to this work, it could enable the creation of optimal detrending functions for PEC and SEC (whether factorized version or not).

There are also many computational lines of inquiry. For example, we are exploring random sampling ("randomized") SVD as a factorization method [34], which can approximate the SVD over large datasets orders-of-magnitude faster than traditional SVD. This development could enable real-time processing during all acquisitions (via the ML:fKK-EC) by initially training with few spectra and retraining when it is calculated that the current basis vectors do not adequately support new data. To bolster this, we are examing methods to update the basis vector set without re-running SVD. Additionally, we are exploring alternative factorization approaches that could enable the ML:fKK-EC to more broadly process chemical components that were not contained within the training data. Similarly, we are working to develop a universal basis set that could be re-used without training on the current sample. We are also exploring "active learning" methods to take advantage of real-time processing that could identify and explore regions of interest during an acquisition.

In conclusion, this work presents the development of a series of new methods for extracting the self-referenced Raman signatures from raw CARS spectra. These new methods, in aggregate, are orders-of-magnitude faster than the conventional implementations and are amenable to high-throughput and even real-time processing with appropriate training data. This advancement facilitates on-the-fly visualization and analysis and would further support such opportunities as in vivo imaging and ad hoc selection of regions-of-interest.