2.1.

Truncation

Raman spectrometers are frequently outfitted with high-pass filters (e.g., interference optical filters) to remove Rayleigh-scattered light that is typically orders of magnitude more intense than inelastic scattering. The filter cutoff must be as sharp and as close as possible to the excitation wavelength to limit excitation light bleed-through while maintaining sensitivity to Raman shifts $<400$ to $500{\mathrm{cm}}^{-1}$. At the same time, the acquisition spectrometer window may extend past the filter cutoff point, toward shorter wavelengths. This results in a spectral shift region at the beginning of every measured signal where the high-pass filter transition and the spectrometer window overlap [Fig. 2(a)], camera pixels 0–50). This region is truncated and removed from the measured accumulations and background.

2.2.

Cosmic Ray Removal

Cosmic rays randomly hit spectrometer camera pixels during the acquisition of Raman signals, resulting in the appearance of sharp artifacts. Short of reducing exposure time, nothing can be done to mitigate the presence of cosmic rays during acquisitions. Worse, the typically long acquisition times required for biological samples (seconds, up to minutes in some cases) makes the presence of cosmic ray artifacts likely. Fortunately, they are relatively easy to remove using one of the following methods.

The first approach, implemented in ORPL’s crfilter_single() function, relies on the localized nature of the cosmic ray artifact. Since their spectral span is usually limited to 3–5 camera pixels, it is possible to (1) use the numerical derivative of the spectrum, (2) identify cosmic ray artifacts using an adaptive threshold, and (3) remove the artifact from the original spectrum using interpolation. This method should be limited to cases where it is impossible to acquire several spectra during an acquisition (i.e., a single accumulation, $N=1$), or for removal of cosmic rays from a background signal. This is because it can be difficult to tune the algorithm parameters to remove every cosmic ray while keeping the signal of interest intact.

The second approach, implemented in ORPL’s crfilter_multi() function, relies on the random nature of cosmic rays. Because it is exceedingly unlikely that two spectra or accumulations exhibit a cosmic ray artifact at the same wavelength (same camera pixel), it is possible to (1) compare intensity across spectra (from different accumulations) wavelength per wavelength to compute the disparity, (2) identify outliers as artifacts, and (3) replace the region neighboring the artifact using interpolation. Related spectra can be the different accumulations from a single site using a point-probe system, or spectra acquired over a small sample region using an imaging system such as a Raman microscope.

In general, the second method should be prioritized whenever possible, as it tends to identify cosmic rays more reliably while minimally affecting the filtered signals and is easier to tune. However, it can only be used effectively on multiple signals or accumulations, which are not always available.

2.3.

Background Removal and Combining Accumulations

Backgrounds are measured before an acquisition to account for contributions to the signals that are not related to the sample, such as ambient light. They should be measured with the same experimental parameters (e.g., exposure time and instrument position whenever possible) as the main acquisitions, but with the excitation source turned off. It is also important to make sure no cosmic ray artifacts are present on a background spectrum to avoid introducing a downward spike in the signal. This can be done using the crfilter_single() on the background before removal when a single background signal is measured, or with crfilter_multi() otherwise. After this, the background is directly subtracted from the spectrum, after normalizing for exposure time. When multiple accumulations are measured for the same acquisition site, they are combined into a unique spectrum by computing their arithmetic mean, leading to improved signal-to-noise ratio (SNR). If that is the case, the background should be subtracted after the accumulation average to avoid unnecessary computational steps.

2.4.

$Y$-axis Calibration

Instrument response correction is a necessary step for any spectroscopic measurement. Typical methods involve measuring using a spectrum from a reference calibration light source of known emission profile to recover the instrument response function (IRF). However, this approach can be challenging to integrate in most Raman spectroscopy acquisition workflows due to the requirement of additional equipment and the difficulty of positioning the calibration lamp.⁴⁰ Instead, an alternate method for calibrating Raman instruments (e.g., correct for filter, detector etaloning, and quantum efficiency effects⁴¹) is based on measurement on a standard reference material (SRM) calibrated and manufactured by NIST (SRM-2241 for 785-nm excitation).³⁵ Because this method does not involve the use of additional instruments (irradiance source), it is better suited for use in clinical environments where time and space are limited. Additionally, reference materials can be encased in a custom 3D-printed enclosure to facilitate positioning of probes or instruments and improve systems calibration repeatability. With this method, the IRF is computed from a Raman spectrum measured on the SRM using experimental parameters (laser power and exposure time) that maximize the sensor’s dynamic range. Then, the measured spectrum is processed using the steps previously described: truncation, cosmic ray removal, background removal. After those steps, the instrument’s IRF is computed from the known theoretical fluorescence response (${\mathrm{SRM}}_{\text{theoretical}}$), given as polynomial coefficients by the manufacturer and the measured signal (${\mathrm{SRM}}_{\text{measured}}$) as

The instrument response is then corrected for by dividing the Raman spectrum acquired by the IRF vector.

2.5.

Baseline Removal

2.5.1.

Standard algorithms

Baseline removal is the most difficult yet most critical step in the processing of a Raman spectrum. Because the probability of Raman scattering is orders of magnitude smaller compared to elastic scattering and fluorescence (in biological materials), even small artifacts introduced during this step can have a disastrous effect and completely overshadow the Raman signal. A common baseline removal method consists in using a polynomial fit (typically of order 5 or 6) coupled with a peak rejection rule (e.g., iModPoly algorithm⁴²). The peak rejection rule is used to exclude regions of the signal that feature Raman peaks. Ideally, the entire Raman signal is excluded, leaving only a smooth baseline to be fitted with the polynomial function. This method is simple and can be effective in cases where Raman peaks are easily detected with simple peak-finding algorithms and if the baseline can be well-modeled by a polynomial function. However, it is not often the case for biological samples. Although intrinsic fluorescence can be adequately modeled as a polynomial function, the same cannot be said for absorption and scattering within the visible-IR range. Furthermore, biological Raman peaks are weak and complex, which makes them far harder to automatically detect compared to pure chemical compounds or inorganic materials. The result is that polynomial methods for baseline removal are difficult to tune and often introduce important artifacts when used for spectra measured on biological samples [Fig. 3(a)]. A method developed by Perez-Pueyo based on morphological processing (referred to as MorphBR moving forward) has shown to be better suited than polynomial fitting for biological applications.⁴³ However, baselines removed using this algorithm feature a jagged staircase effect and a “hill,” that is as wide as the filtering window, near the beginning of the spectrum [Fig. 3(b)].

Fig. 3

Demonstration of problems with common baseline removal algorithms. The spectrum used in this example was measured with a Renishaw Raman microscope on saliva samples.¹³ (a) iModPoly (polynomial fit) and (b) MorphBR (morphological baseline removal).

2.5.2.

BubbleFill algorithm

Inspired by the approach of Perez-Pueyo (MorphBR), we developed a new algorithm based on morphological processing that results in a smoother baseline fit, introduces fewer artifacts and for which fitting parameters can be tuned to different levels across the $x$-axis. We named this new algorithm BubbleFill because the fitting process uses circular bubbles to fill the region underneath the spectrum. A detailed flowchart of the algorithm is shown in Fig. 4 and an illustration of the bubble growth loop is presented in Fig. 5. First, the overall slope of the spectrum is removed using a linear fit and the result is scaled for the $x$-axis and $y$-axis to span the same range (square aspect ratio). The baseline estimate is initialized at 0 over the entire $x$-axis. Then, circular bubbles are iteratively grown under the scaled spectrum, starting with a bubble of diameter equal to the spectrum’s length aligned on the center of the $x$-axis. The bubble pops when it reaches the spectrum, and two new bubbles start to grow on each side of the contacting point. A bubble’s diameter and alignment depend on the $x$-axis region where it is grown. After every bubble has popped, the new baseline estimate is updated as the maximum value between the bubble and the current baseline. This process repeats until every new bubble has reached the critical minimum diameter specified as a tuning parameter. Finally, the baseline estimate is smoothed using a Savitzky–Golay filter and scaled back to the original spectrum $y$-axis range. The final Raman signal is then obtained by the subtraction of the baseline estimate to the original spectrum. The only tuning parameter corresponds to the smallest allowed bubble diameter. The smaller they are allowed to grow, the more aggressive is the baseline fitting process and vice versa. Additionally, it is possible to specify different bubble diameters across the $x$-axis, effectively achieving a different degree of sensitivity over different regions of the input signal.

Fig. 4

Flowchart of the BubbleFill baseline removal algorithm. The algorithm inputs are: (a) A raman signal ($S_0$) represented as a vector of length $N$; (b) the minimum width of the bubbles allowed to grow (${\varphi }_{\min }$); and (c) the order of the polynomial fit used to remove the global slope of the signal ($f_{\text{order}}$). First, the global slope of the signal ($S_0$) is removed using a polynomial fit and the signal intensity is normalized to obtain a square aspect ratio ($\min (S_0)=0$, $\max (S_0)=1$). Then, bubbles of increasingly smaller size are grown underneath the signal and the baseline fit ($B$) is iteratively updated. Once the bubble growth loop is completed, the square aspect ratio normalization is reversed and the baseline fit is smoothed using a Savitzky–Golay filter. The output Raman signal is computed as the subtraction of the baseline fit to the input signal.

Fig. 5

First six iterations of the bubble growth and baseline update iterative process that is at the core of the BubbleFill algorithm. Only the first six iterations are displayed, resulting in a small, obtuse angle toward the right side of the peak. Additional iteration reveals that a small bubble is eventually grown in that region and significantly improve the baseline fit. The Savitzky–Golay filter applied to the baseline fit after the bubble growth process further smooths possible remaining sharp angles.

2.6.

$X$-axis Calibration

For the same reason that the $y$-axis calibration of a Raman instrument should ideally not be made using a calibration light source, the $x$-axis Raman shift needs to be computed from a reference sample’s spectrum. Reference samples used for this purpose should have a Raman to baseline ratios (RBR) that is the ratio between the intensity of the tallest Raman peak to the maximum of the baseline signal, of at least 0.2. Additionally, it is preferable to use samples featuring narrow Raman peaks that can be easily identified with an automatic peak-finding tool, such as SciPy’s find_peaks() function. Acetaminophen and Nylon are two reference samples well suited for this purpose. They both feature narrow peaks that are uniformly spread over the 0 to $2000{\mathrm{cm}}^{-1}$ region, have an RBR ratio above 0.5, are shelf stable while not requiring tedious maintenance and can easily be brought into sterile environments.

To calibrate a Raman spectrometer’s $x$-axis, a spectrum is measured from a chosen reference sample, then the spectrum is processed using the steps described in Secs. 2.1–2.5. The camera pixel location index of the most prominent Raman peaks are identified using a peak-finding tool [e.g., SciPy’s find_peaks() function]. Finally, a polynomial fit of order 2 or 3 is used to create a map conversion between camera pixel index and Raman shifts. An important note, when calibrating the $x$-axis for machine learning applications, is that it is critical to interpolate every spectrum over a common $x$-axis. Otherwise, there is no guarantee that all spectra are expressed in the same vector space or, in other words, that the intensity vectors across spectra correspond to the same wavelength/energy.

2.7.

Smoothing and Normalization

Smoothing and normalization are highly dependent on the end use of a Raman spectrum. When spectra are intended for visual analysis, they are smoothed using an average moving filter and normalized so that the resulting spectrum’s minimum is 0 and maximum, or intensity at a chosen Raman shift, is 1. This type of normalization helps with visual comparison and assessment of spectra acquired on samples that cover a wide range of absorption and scattering. When spectra are used in machine learning applications, however, they are not smoothed and are normalized using the standard normal variate (SNV) method. The SNV normalized ($s^{\star }$) transformation of a signal (s) is given as

2.8.

Spectrum Quality Factor

The definition of a general spectra quality factor is necessary for comparison and ranking of signals from different datasets. However, common metrics such as the SNR or the RBR tend to be highly dependent on the instrument (spectrometer slit and resolution and detector efficiency) and acquisition parameters (excitation power, and exposure time). Because of this, Raman signals measured from different sample types or using different systems tend to have vastly different SNR and RBR. Often, a top-quality spectrum from one dataset would have an SNR and RBR smaller than every spectrum from another dataset. It is therefore impracticable to use either metric for general spectrum quality comparison across multiple datasets or instruments. Instead, the quality of the Raman spectra presented in this work was measured using the average signed squared intensity (ASSI) as defined as

3.1.

Comparison of Baseline Removal Algorithms: Synthetic Tool

Synthetic spectra ($S_i$) were modeled as a combination of three signal components: Raman ($R_i$), baseline ($B_i$) and noise ($N_i$). Generation of benchmark signals in this manner allows a fine control on the desired SBR and SNR. Additionally, outputs given by a baseline removal algorithm can be compared to the original signal components with which the input spectrum was generated to compute overall fitting error. However, the Raman and baseline components used for the generation of synthetic spectrum need to be independent of the baseline removal algorithm to be tested to limit possible biases. The formula used to generate synthetic spectra is

The Raman component used for generation of benchmark spectra were experimentally measured on acetaminophen, Nylon and PDMS samples using a point probe system⁴⁴ and were processed using the aforementioned workflow. The experimental Raman spectra were hand-fitted as a superposition of Gaussian curves [Fig. 6(a)] to ensure only clean Raman peaks remained, and to limit eventual biases introduced by the baseline removal algorithm. For the baseline component, spectra were measured on aluminium and nigrosin. As neither of the signals featured noticeable Raman peaks, it was concluded that their spectral responses consisted essentially of pure fluorescence (baseline) signals. These experimental baselines were smoothed using a 50-pixel wide average moving filter to further remove residual traces of noise resulting in the Baseline components ($B_i$) used for the generation of benchmark spectra. Together, the two baseline components and three Raman components [Fig. 6(b)] enable the creation of a wide variety of possible synthetic spectrum suited for use as quantitative benchmark tests for baseline removal algorithm validation and optimization [Fig. 6(c)]. Finally, the noise component ($N_i$) is generated following a normal distribution with an average of 0 and a specified standard deviation.

Fig. 6

Baseline algorithm benchmark with synthetic spectrum generated from the aluminium baseline and acetaminophen Raman with a RBR of 0.05. (a) Manual fitting of a signal composed of Gaussian peaks over a experimental spectrum. (b) Synthetic benchmark spectra can be created by combining a Raman signal component (nylon, PDMS or tylenol) and a baseline signal component (nigrosin or aluminium). (c) The Raman to Baseline Ration can be adjusted to be representative of different tissue types or applications.

The metric chosen to evaluate the performance of the different baseline removal algorithms implemented in ORPL is the normalized mean squared error (nMSE). The nMSE between a Raman signal computed from a baseline removal algorithm ($R_i^c$) and the Raman target ($R_i^t$) used in the generation of a benchmark signal is defined as

3.2.

Baseline Removal Benchmark without Noise

A benchmark spectrum is generated using the acetaminophen Raman and aluminium baseline components with an SBR of 0.05. The spectrum’s baseline is removed with the BubbleFill, MorphBR and iModPoly algorithms and their respective Raman outputs are compared to the target. The results presented are the optimal fit of each algorithm (smallest nMSE) obtained by sweeping the possible range of value for their respective tuning parameter. This method is used to compare baseline removal algorithms in a ‘best-case scenario’ leveled playing field.

This first example highlights the shortcomings of polynomial fitting methods when removing baselines that feature “localized bumps,” such as the ones seen in many biological samples (Fig. 7). Even though not all biological sample feature a baseline as difficult to remove as aluminium, localized bumps can still be introduced in the measured signal when acquisitions are made over an aluminium substrate. This is relevant and problematic because aluminium slides are becoming an increasingly popular and cheaper alternative to calcium fluoride (${\mathrm{CaF}}_2$) slides for Raman microscopy.^11,41,45 The Raman spectrum given by the BubbleFill algorithm closely matches the target except for the two smaller peaks near $1500{\mathrm{cm}}^{-1}$. In most instances, using a different tuning (e.g., larger bubbles) would result in fewer small peaks being removed at the cost of a worse baseline fit. This is the general trade-off when tuning any baseline removal algorithm. However, with BubbleFill, it is possible to use a different tuning for different regions of the spectrum. For instance, the bubbles grown in the 1400 to $1600{\mathrm{cm}}^{-1}$ region could be much larger than for the rest of the spectrum. This would result in a similar baseline fit, but the small peaks near $1500{\mathrm{cm}}^{-1}$ would be preserved. This feature is, to our knowledge, unique to BubbleFill and could lead to significant improvements to baseline removal in some applications. Furthermore, additional testing with different combinations of baseline and Raman components and for different SBRs indicate that BubbleFill outperforms other tested algorithms in numerous instances (Fig. 8). The largest difference in performance was observed for spectra generated with the aluminium baseline and acetaminophen or nylon Raman components. In these cases, BubbleFill has an nMSE nearly two orders of magnitude smaller than the second-best performing algorithm. It fell behind only for spectra generated with nigrosin baseline and nylon or PDMS Raman components where it still remained the best option when the RBR approached 0.

Fig. 7

Baseline algorithm benchmark with synthetic spectrum generated from the aluminium baseline and acetaminophen Raman. Signal-to-noise : 0.05. Noise standard deviation : 0.

Fig. 8

nMSE as a function of the RBR for the three algorithms tested. Each graph is for spectra generated with a different combination of baseline and Raman signal components.

The different algorithms tested in this work along with the results presented in this section are summarized in Table 1.

Table 1

Baseline algorithm comparison summary. Execution time was measured on a benchmark spectrum generated with the Nylon Raman, aluminium baseline, Raman/baseline ratio of 0.5 and noise standard deviation of 0.1.

3.3.

Baseline Removal Benchmark with Noise

Baseline removal algorithms have been tested on spectra with noise to confirm that the results shown thus far can translate to real-world applications. However, the method of comparing algorithms’ performance via the mean squared error becomes impractical when noise is added to the input signal. The problem is that the MSE between an algorithm’s computed Raman and the target is overwhelmingly correlated with the added noise itself. This makes the computed MSE seemingly identical between the three algorithms in every scenario tested, even when differences remained visually noticeable in the respective outputs. Instead of comparing nMSE across the different algorithms on noisy spectra, two examples have been chosen to illustrate tendencies also observed in the experimental datasets presented in the following section. In both examples, the input signal was generated using baseline and Raman components as described before, but with the addition of noise. The noise signal added follows a normal distribution with average 0 and standard deviation of 0.01, i.e., 1% of the baseline’s maximum.

Figure 9 shows the baseline removal benchmark on a spectrum generated from nigrosin and nylon with Raman-to-noise ratio of 0.15. This example was chosen because it is a best-case scenario for the iModPoly algorithm—the nigrosin baseline is smooth, and the nylon peaks are easily identifiable and located toward the center of the spectrum. Yet, both BubbleFill and MorphBR performed similarly to iModPoly.

Fig. 9

Baseline algorithm benchmark with synthetic spectrum generated from the nigrosin baseline and nylon Raman. RBR : 0.15. Noise standard deviation : 0.01.

Figure 10 shows the baseline removal benchmark on a spectrum generated from aluminium and PDMS. Similarly to the tests performed on signals without noise, iModPoly failed to correctly remove the aluminium baseline while BubbleFill and MorpthBR managed a near-perfect recovery of the Raman target except for a small discrepancy (shared by every algorithm) for the first peak near $500{\mathrm{cm}}^{-1}$. In general, none of the algorithm’s performance was significantly affected by the addition of noise to the input spectrum.

Fig. 10

Baseline algorithm benchmark with synthetic spectrum generated from the aluminum baseline and PDMS Raman. RBR : 0.15. Noise standard deviation : 0.01.

3.4.

Testing ORPL with Different Experimental Datasets

We tested the capability of ORPL on four datasets acquired with three different Raman platforms to cover a wide range of sample types, instruments, and sample preparation methods. The first dataset (1719 spectra) consists of in vivo brain tissue spectra measured with the hand-held probe from Reveal Surgical and provides a good reference for spectra acquired in a surgical workflow. The second dataset (524 spectra) consists of in vivo and ex vivo prostate tissue spectra measured with a custom lab-built system using a commercial EmVision LLC handheld probe. The third (8670 spectra) and fourth datasets (7774 spectra) consist of paraffin-fixed prostate tissue slices and dried saliva samples respectively and were measured with a commercial Renishaw Raman microscope as 3D Raman maps. Combined, the four datasets amount to a total of 18,687 individual acquisitions (not counting repeated accumulations) that cover a wide range of signal-to-background and SNR as confirmed by the spectral quality of each signal measured using the ASSI metric defined earlier (Fig. 11).

Fig. 11

Spectral quality histogram of the experimental datasets (brain, prostate Point Probe, prostate Wide-Field and saliva) measured with the ASSI quality metric.

The raw data (first row of Fig. 12) shows that the spectra of each dataset covered nearly the entire dynamic range of the acquisition instruments. This is frequent with biological samples and is explained by the large variability of fluorescence signal strength, absorption, and scattering. Nevertheless, Raman pre-processing using the workflow presented in this work and the BubbleFill algorithm for the removal of intrinsic fluorescence (step 3 of Fig. 12) resulted in averaged spectra with small deviations. Many common Raman active bands can be identified and are common across all datasets, including vibrational modes typically found in proteins and lipids.

Fig. 12

Raw data: measured raw spectra from instruments with accumulations combined. Step 1: Truncation (for datasets 1 and 2) and cosmic ray removal. Step 2: Background removal (for datasets 1 and 2) and calibration of $x$ and $y$ axis. Step 3: Baseline removal with BubbleFill. Step 4: SNV normalization. Average Spectra: average spectra computed from step 4 results with +- standard deviation zone represented as shadow.

Finally, spectra of all datasets were clustered in groups of high, average, and low quality based on the ASSI metric (Fig. 13). As the quality increased, the standard deviation became smaller, converging toward the average spectrum. This behavior was observed for all datasets. Furthermore, some specific Raman bands (peak at $1300{\mathrm{cm}}^{-1}$ for the prostate wide field dataset) gained in intensity while artifacts (400 to $600{\mathrm{cm}}^{-1}$ for the prostate point probe dataset) disappeared. These results indicate that the ASSI metric has the potential to be used in machine learning applications to discard low-quality spectra based on a threshold or as a general signal quality metric during acquisition to facilitate troubleshooting of instrument and software.

Fig. 13

Average pre-processed spectra of each experimental tested datasets clustered in high (top 20% ASSI), average (middle 20% ASSI) and low (bottom 20% ASSI) quality.