Real-time and high-throughput Raman signal extraction and processing in CARS hyperspectral imaging

We present a new collection of processing techniques, collectively"factorized Kramers--Kronig and error correction"(fKK-EC), for (a) Raman signal extraction, (b) denoising, and (c) phase- and scale-error correction in coherent anti-Stokes Raman scattering (CARS) hyperspectral imaging and spectroscopy. These new methods are orders-of-magnitude faster than conventional methods and are capable of real-time performance, owing to the unique core concept: performing all processing on a small basis vector set and using matrix/vector multiplication afterwards for direct and fast transformation of the entire dataset. Experimentally, we demonstrate that a 703026 spectra image of chicken cartilage can be processed in 70 s (approximately 0.1 ms / spectrum), which is approximately 70 times faster than with the conventional workflow (approximately 7.0 ms / spectrum). Additionally, we discuss how this method may be used for machine learning (ML) by re-using the transformed basis vector sets with new data. Using this ML paradigm, the same tissue image was processed (post-training) in approximately 33 s, which is a speed-up of approximately 150 times when compared with the conventional workflow.


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][2][3][4][5][6]. With an imaging speed as fast as 50 000 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][9][10][11]. The NRB, however, behaves as a stable homodyne amplifier for the Raman-generated signal; thus, reducing the NRB also reduces the Raman 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; thus, it may amplify weak signals above the noise floor. Furthermore, there is 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 generate robust, quantitative CARS Raman spectral data 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 datathough this is readily testable as will be discussed.

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 C ARS , may be described as [1]: where E p , E S , and E pr are the frequency-domain (ω) pump, Stokes, and probe fields, respectively; χ (3) is the third-order nonlinear susceptibility, which is a summation of resonant ( χ r ; Raman vibrational) and nonresonant ( χ nr ; electronic) components, andC st is the system response function that incorporates such properties as laser source profiles, optical filter transmission profiles, and detector response. In the right-hand part of Eq. 1, the tilde above χ (3) is used to indicate that the nonlinear susceptibility is convolved with the probe [17]; though, in the remainder of this manuscript, it will not be explicitly used. '⋆' and ' * ' are the cross-correlation and convolution operators, respectively. The overarching goal of phase retrieval methods is to extract Im{ χ (3) (ω)} from I C ARS (ω), which is the equivalent material property probed by traditional Raman spectroscopy [21]. If C st (ω) and I N RB (ω) 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 C ARS (ω) χ (3) (ω)/ χ nr (ω). Using the KK formalism and assuming, for the moment, that the I N RB (ω) of the sample itself is measurable [17]: whereĤ is the Hilbert transform. To approximate the NRB, one uses a surrogate/reference material with nonlinear susceptibility χ r e f (ω), which leads to a CARS signal, I r e f (ω). One can model this relationship between the actual NRB and the surrogate as I r e f (ω) = Ξξ(ω)I N RB (ω), where ξ(ω) is a frequency-dependent function and Ξ is a constant. Both are real valued. It should be noted that these terms encompass differences in both the material properties as well as any optical system response changes (e.g., related toC st ). Applying this new scenario to Eq. 2: From this equation one will notice that the use of a reference material has led to a multiplicative amplitude error and an additive phase error [17], which are themselves related by a KK relation. Thus, baseline detrending of Im{K(ω)} is not appropriate. There is a solution: PEC. Under the assumption of a slowly-varying ξ(ω), one may find the phase error using detrending methods, such as asymmetric least squares (ALS) [22,23], and remove it and the associated amplitude error (within a constant Ξ): A err (ω) = 1/ Ξξ(ω) = exp −Ĥ {φ err (ω)} /Ξ. PEC does not account for and remove Ξ as the Hilbert transform of a constant is 0 (i.e.,Ĥ {ln Ξξ(ω)} =Ĥ {ln ξ(ω)}). Finding the constant Ξ is the role of SEC [17]. This may be calculated from the real part of K(ω) after PEC: where ' · · · ' indicates the mean over the frequency. Due to computational distortion of the numerical Hilbert transform, one usually does not simply use the mean but rather a trendline [17]. In summary, using the KK relation, PEC, and SEC, one can calculate K C ARS (ω) from K(ω) as: Thus, without directly measuring the NRB, one can find the ratio χ r (ω)/ χ nr (ω) at every pixel because every pixel is self-referenced to its own nonresonant component. The full conventional workflow is illustrated in Fig. 1(a). This ratio is maintained even in the presence of absorption and scatter as both the Raman and NRB components are equally affected; though, the SNR deteriorates.

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 Figure 1 The first step in this process is factorization of the input data. In this work, SVD decomposes a matrix A into three matrices as A = USV H . The H-superscript indicates the Hermitian transpose; U and V are unitary matrices whose columns are the left-and right singular vectors, respectively; and S is a diagonal matrix whose entries are known as singular values (we will denote the vector containing just the singular values as 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 U is composed of spatial basis vectors while V, spectral basis vectors. Further, 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 U, S, and 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 U, S, or V were altered for denoising; though, all derivations remain valid.
For the fKK, one would conceptually apply the KK relation to the spectral basis vectors in 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 ln I

[m]
C ARS (ω)/I r e f (ω) = a m (ω), where the m-superscript denotes the m th spectrum, which leads to the a m -vector. For the following derivation, we assume that we have M spectra, and the ω vector is N-frequency increments long. Thus, A may be written as: Assuming a "reduced" SVD implementation, we have more spectra than the length of the frequency vector ω (i.e., M > N); thus, U ∈ R M, N , s ∈ R N ; diag(s) = S ∈ R N, N , and V ∈ R N, N . As the Uand S-elements act as constant weighting terms to the right singular vectors (V's columns) and the Hilbert transform is a linear operator [26], this is equivalent to only applying the transform to the right singular vectors: The total fKK process without PEC or SEC may be described 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 (ω) to A(ω), 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 an image of all the same spectrum, the standard deviation (σ A (ω)) of the previously defined A may be approximated as [28]: where · · · indicates the mean spectrum, α is a Poisson noise multiplier, and σ g is the standard deviation of the additive white Gaussian noise. We have assumed that the I r e f (ω) used is effectively noiseless as the reference spectra is often an averaged and/or denoised version of repeated measurements of a surrogate material. Applying this scaling term, the fKK would be re-written as: I r e f (ω) for all m (10) where A m,: is the m th spectrum (row) in A.
In the remainder of this manuscript, we will include f (ω) 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 (ω) 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 A or V T ).

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: where D is a detrending operator, and Φ PEC is a basis set describing phase error. We do not want to detrend every spectrum as described in the proceeding equation and the orthonormal V singular vectors are not readily usable for baseline detrending as they often have positive and negative values with no clear baseline. Rather we will take the approach of sub-sampling U (to form U ss ), calculate Φ err , and regress to approximate Φ PEC . This dramatically reduces the computational burden compared to using the full U. Our current practice, inspired by vertex component analysis (VCA) [29], is to sub-sample U by keeping the rows of U that have the highest and lowest values for each column:, and optionally a sub-sample between. For the maximum and minimum: U ss = U q,: where the ':' indicates all row or column entries, and q = q min ∪ q max indicates the union row indices. q can also contain a sub-sample between the max-and min-values for each column U. From this: where λ is a non-negative scalar regularization weight and I is an identity matrix. The left-hand statement in Eq. 17 is an ordinary least-squares regression using a [pseudo]-inverse. In practice, however, this result is unstable owing to significant multicollinearity in the singular vectors. These collinearities cause erroneously large Φ PEC entries, especially those corresponding to the smallest singular values. One solution to this problem is ridge regression (also known as Tikhonov regularization) as shown on the right side of Eq. 17. The action of the combined fKK and fPEC without fSEC can be described as: noting that the amplitude and phase terms are still related by a Hilbert transform.

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: It should be noted that the mean of the first expression in the previous equation can be solved analytically, for example, using partial fraction decomposition, assuming that χ r (ω) = m A m /(Ω m − ω − iΓ m ), A m , Ω m , and Γ m are real and positive-valued, and χ nr is constant and positive, real-valued.
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: where Ξ ∈ R M is a vector of constants.
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 V T SEC into the previous expression: A solution for this matrix is the subtraction of the mean of the PEC-corrected right singular vector: Thus, if the mean of each corrected right singular vector is zero, the mean of the magnitude will also be zero. As we previously mentioned, due to numerical errors in the Hilbert transform, rather than a strict mean, we use a trendline function, which was previously implemented as a large-window, small-order Savitzky-Golay filter [17]. Thus: where M is a trendline (or mean function).

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 K C ARS (akin to Eq. 5 for the conventional implementation). Applying Eqs. 8, 18, and 22: again noting that U, S, and V may be reduced in size from the original SVD for the purposes of denoising.

The 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 (ω), U, S, V, Φ PEC , and V T SEC . One assumes that upcoming images will comprise the same 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 U new (which describes the spatial mixtures of SV T ). From Eq. (7) for the new data, A new that incorporates f (ω) as well, and solving for U new applying ridge regression: Now, one can simply apply the U new to Eq. 23. 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 ∼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.

Broadband CARS (BCARS) imaging platform and software
Images were collected on an in-house-developed BCARS microscope that is described in detail elsewhere [1]. The picosecond probe laser and femtosecond supercontinuum were 13 mW and 7.1 mW on-sample. The CCD integration time was set to 3.5 ms, which corresponds to 5 ms per pixel owing to data transfer time, stage movement, and CCD cleaning time.
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.

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 µm to 40 µm, as measured by BCARS imaging ("XZ" images).

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.

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.

Simulation
We simulated a noiseless 3-chemical mixture with the concentration map shown in Figure 2 (a) and a ternary plot of concentrations shown in 2 (b). Chemical 1, 2, and 3 are displayed in red, green, and blue, respectively. The base dataset is 74 pixels x 246 pixels (18 204 total spectra). To analyze the fKK-EC performance versus number of spectra, this dataset is sidescaled by a factor of 0.5, 1, 2, 3, and 4; for a total of 4 551, 18 204, 72 816, 163 836, and 291 264 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 χ nr (ω)'s were quadratic polynomials with randomly generated non-negative coefficients, and χ r e f (ω) 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 2 500 cm −1 sampled 810 times; though, Raman peaks could only be assigned between 500 cm −1 to 1700 cm −1 . The stimulation profilẽ C st (ω) in Eq. 1 was set to a constant for simplicity. 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 A × max(M, N) × ǫ), where M and N are the row and column dimensions of the SVD-input matrix A, and ǫ 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: ≤100 µs for SVD and selecting basis vectors, ≤140 µs for the KK, ≤3.2 ms for the PEC, and ≤140 µs for the SEC; for a total of ≤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 ≥40 for all but the smallest dataset. For the 291 264 spectra simulation, for example, the total time was <25 seconds for all 3 replicate simulations (86 µs / spectrum). The most significant difference is the time to perform phase retrieval, with the conventional KK requiring ≈ 40 s and the new fKK ≈ 4.3 ms -an over 9000× improvement. The fPEC was over 1250× faster than the PEC, and the fSEC was over 3150× faster than the SEC. For the factorized workflow, the reconstruction step only added 3.3 s. Fig. 2(f) 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 C ARS in Eq. 5 or Eq. 23) and the known Im{ χ/ χ nr } at each pixel. The mean RSS, RSS , is shown in Fig. 2(e). For reference, the RSS if K C ARS (ω) = 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. Figs. 2(d)-(f) 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 ( RSS <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). Fig. 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 ≈150× speedup, which was <30 µs per spectrum for all dataset sizes. Thus, this could be performed in real-time as the data is acquired. Fig. 2(e) shows that the machine learning implementation provides equivalent RSS to the non-ML fKK-EC method. Fig. 2(f) shows the computational fraction of each step. Finally, Fig. 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 ≈3% overlap (per side) with neighboring images. The stitched image is 846 pixels by 831 pixels (703 026 pixels total). Fig. 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 716 cm −1 , subtracting a linear interpolated baseline between (691 to 738) cm −1 . Tentatively, we assign this peak to the nucleotide adenine [30]. 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 [31]) peak relative to the trough at 900 cm −1 . Lipids were highlighted using the intensity at 2837 cm −1 (CH 2 -symmetric stretch [32]). 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 ≈70× faster using the fKK-EC methods (average of 3 repeats ± 1 standard deviation: conventional method ≈ 4973 s ± 26 s total [≈ 7.0 ms / spectrum]; fKK-EC ≈ 70 s ± 3.0 s total [≈ 99 µ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 portions 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 78 114 spectra. 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× faster than the conventional workflow, requiring <50 µs / spectrum to process the entire image. Though these images could have been analyzed in real-time, they were processed after acquisition.

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 χ 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 [33], 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. Additionally, we are looking into methods to create a universal basis vector set that could be re-used without training on the current sample. We are also exploring active learning machine 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 highthroughput 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.