1. Introduction

In the last decade coherent Raman scattering (CRS) microscopy has emerged as important label-free and chemically sensitive imaging technology for biological systems [2–4]. Coherent Raman scattering probes intrinsic vibrational resonances to provide chemical contrast, which alleviates the need of perturbing the natural environment of biomolecules by the addition of labels such as fluorophores used for fluorescence imaging.

The two most common CRS imaging techniques are coherent anti-Stokes Raman scattering (CARS) [2, 3, 5] and stimulated Raman scattering (SRS) [6–9]. In both methods, the coherent vibration of the biomolecules is driven by the interference between two electromagnetic fields called pump and Stokes fields. The pump field is then Raman scattered by the driven vibrations. In CARS, the anti-Stokes part of the scattered field is acquired which is free from excitation background, while in SRS the scattering is measured by the loss of the pump beam or the gain of the Stokes beam. The nonlinear character of CARS and SRS enables sub-micron three dimensional resolution and the coherent superposition of the scattering of all molecules in the focal volume yields a significantly increased signal compared to spontaneous Raman scattering, allowing for fast imaging.

In both techniques, the measurement of a vibrational spectrum over a significant number of spectral points is required to provide a good chemical specificity in the presence of several chemical compounds. Such hyperspectral vibrational microscopy can be achieved by sequential imaging at different vibrational frequencies [10, 11]. For CARS intensity data, measuring the vibrational spectrum over a sufficiently large spectral range is actually required to extract quantitative chemical information. This is due to the fact that the quantity which is linear in the chemical composition is the complex susceptibility, while the CARS intensity is proportional to the absolute square of the susceptibility, having a modified spectrum compared to spontaneous Raman [4]. Several methods are available to retrieve the complex susceptibility spectra from CARS intensity spectra [1, 12, 13]. These methods determine the phase of the susceptibility from the CARS intensity, which is a spectrally non-local operation, and consequently the resulting errors decrease with increasing CARS spectral range [1]. Such a retrieval method is not required for SRS which is proportional to the imaginary part of the susceptibility. However, the SRS signal can contain significant background due to Kerr-lensing, thermal lensing, or induced absorption [14]. Furthermore, we have recently shown that the real part of the susceptibility measured in CARS provides an internal calibration standard, and can also be used to distinguish chemical substances which do not show vibrational resonances in the measured spectral range [1].

To date, the fastest SRS imaging was shown in [15], acquiring images of 500 × 500 spatial pixel at a rate of 30 frame per seconds (fps) with frame-by-frame frequency tunability. Considering the 91 frequencies acquired in the range 2800–3100 cm⁻¹ and 10 repetitions needed for sufficient signal-to-noise ratio, this corresponds to 0.03 hyperspectral frames per second (hfps). CARS imaging can reach similar speeds [2, 5]. Significantly faster rates are required for real-time imaging with quantitative chemical information, and for high-throughput high content screening, where a large number of similar structures need to be measured to provide statistically relevant results.

Here we demonstrate a method to reduce the number of measured frequencies in the hyper-spectral image, and consequently to increase the hfps rate, without losing significant spectral information. Such sparse sampling methods reconstructing signals from far fewer samples than Nyquist sampling have recently attracted considerable attention [16], and applications cover many areas, such as geological remote sensing and medical imaging. Our work is developing and applying a sparse sampling method for coherent Raman hyperspectral imaging.

2. Experimental

CARS hyperspectral images have been acquired on a home-built multi-modal laser-scanning microscope based on an inverted Nikon Ti-U. A detailed description of the set-up can be found in [10]. A broadband (660–970 nm) laser beam from a 5 fs Ti:Sa laser is split into pump and Stokes beams for CARS excitation with wavelength ranges of 660–730 nm and 730–900 nm, respectively. Hyperspectral imaging is achieved by spectral focussing [17–19]. The pump and Stokes pulses are equally linearly chirped, resulting in a constant instantaneous frequency difference (IFD) within the pulse duration. Excitation of different vibration energies in the range 1200–3800 cm⁻¹ with a resolution of 10 cm⁻¹ is achieved by controlling the delay between pump and Stokes. The data discussed in this paper were taken over a (2400–3700) cm⁻¹ range with a 60× 1.27 NA water immersion objective (Nikon CFI Plan Apochromat IR λS series) and a 1.4 NA oil immersion condenser for signal collection in forward direction, with a resulting spatial resolution for the CARS susceptibility (full-width at half-maximum of the point-spread function amplitude) of 0.6 (2) μm in the lateral (axial) direction. The CARS light is discriminated by a pair of Semrock band-pass filters FF01-562/40 and then detected by a Hamamatsu H7422-40 photomultiplier. CARS hyperspectral images of a human osteosarcoma U2OS cell were acquired for IFDs in the range (2400–3700) cm⁻¹ with 5 cm⁻¹ steps. An example image at 2805 cm⁻¹ is shown in Fig. 1.

 

Fig. 1 CARS intensity I_C image of a human osteosarcoma U2OS cell acquired at an IFD of 2805 cm⁻¹. Logarithmic gray scale as shown. The pixel size is 0.1 × 0.1 μm², and the pixel dwell time was 10 μs. The pump and Stokes beam power at the sample were 53 mW and 31 mW, respectively. Spectra at the positions indicated in the Fig. are shown in Fig. 3. The scale bar represents 5 μm. The full dataset of $\sqrt{I_{\text{C}}}$ is shown in the panel a) of the hyperspectral movie Media 1.

Download Full Size | PDF

3. Results and discussion

The method proceeds as follows. First a hyperspectral data set A with a large spectral range and high spectral resolution is taken over a small number of spatial points. Singular value decomposition is used to determine the important spectra above the noise in A. Then a number of important spectral points equal to the number of important spectra is determined to yield minimum deviation over A in the spectra reconstructed from the important spectral points. Subsequently, hyperspectral data sets over large number of spatial points taken only at the important spectral points can be reconstructed.

3.1. Determination of important spectra and spectral points

A hyperspectral CARS intensity (I_C) image A with a spectral step size at the Nyquist-limit of the spectral resolution is acquired using either a coarse spatial sampling or a small region of the sample containing a representative ensemble of the chemical components present. The subsequent processing is done on the square root of I_C, which has a whitened shot noise [1]. The $\sqrt{I_{\text{C}}}$ data are formatted as a (S × P′) matrix with the elements d_s,p, where S is the number of spectral points and P′ is the number of spatial points measured. The matrix d is then factorized using a compact singular value decomposition (SVD) into d = uσv^T, where the (S×S) matrix u is the rotation from the new spectral basis given by the singular spectra, to the old spectral basis given by the measured IFDs. The (S × S) matrix σ is diagonal containing the singular values, and the (S × P′) matrix v^T contains the spatial images in the basis of the singular vectors. As example we use as hyperspectral data A the data shown in Fig. 1 having S = 261, but with a coarse spatial step of 1 μm, containing only 1% of the spatial points of the frame, yielding P′ = 951.

We determine the number S_max of singular values not dominated by noise using the algorithm reported in [1]. It consists of fitting the singular values σ_j,j (see Fig. 2(a) with a linear dependence on the index j for j ≥ S/2, and determine the largest index j = S_max for which σ_j,j is larger than $\sqrt{2}$ times the fit. For the example in Fig. 2(a) we find S_max = 8. Components with j > S_max are attributed to noise and we discard them by defining a σ̃ in which the diagonal values are set to zero for j > S_max and define the noise-filtered data d̃ = uσ̃v^T.

 

Fig. 2 a) Diagonal elements of the singular value matrix σ of dataset A, normalized to the largest value. The solid line is a linear fit to σ_{j, j} for j ≥ S/2. The dashed line is the linear fit multiplied by η. b) Error ε in the spectral reconstruction of the data as a function of the number S′ of spectral points used, for the different methods discussed in the text as labeled. The inset shows the evolution of ε during the random walk for S′ = 8, where i indicates the step number. The green surrounding lines bars indicate the standard deviation of ε over 10 independent repetitions of the walk.

Download Full Size | PDF

In order to reduce the number of spectral points to be measured from S to S′ ≪ S, we use only the most important S′_max ≤ S_max singular spectra, and determine their weights in the reconstruction by measurements at S′ ≥ S′_max spectral points. S′_max and the specific spectral points to be measured are determined from d̃. In detail, we create a (S′ × P′) matrix d̂ by selecting S′ out of the S spectral points of d̃, which are given by a vector s = (s₁,.., s_S′) of ascending spectral point indices, i.e. d̂_n,p = d̃_sn,p with n = 1..S′. An approximation of the matrix d̃ can be obtained by determining the weights of the S′_max most important singular spectra using d̂. To do this, we determine a S′ × S′ rotation matrix û by taking the first S′ singular spectra of u at the spectral points identified by s, and using Gram-Schmidt orthogonalization of the resulting matrix in sequence of increasing singular spectrum index, i.e. decreasing importance. We then calculate a reduced matrix of coefficients â such as d̂ = ûâ. Having the coefficients of the singular spectra, we can reconstruct the data over all spectral points using the (S × S′) sub-matrix ũ of u containing the first S′ singular spectra. Similar to the SVD noise filtering, we restrict the reconstruction to S′_max ≤ S′ important spectra by setting in ũ the singular spectra above S′_max to zero. The noise filtered reconstructed data is then given by d^rec = ũû⁻¹d̂

To determine S′_max for the reconstruction, we have to consider that S_max was determined using the noise in a dataset with S spectral components. Measuring only S′ spectral points, the noise is expected to be a factor of $\sqrt{S/S^{\prime }}$ larger. Furthermore, we note that we can choose the spectral points s to minimize the influence of noise in the reconstruction. We accordingly determine s by minimizing the error of the reconstructed hyperspectral image ε = ||d̃ − d^rec||_F/||d̃||_F, where ||.||_F indicates the Frobenius norm, and S′_max = S_max is used to determine d^rec. For a typical number of S in the order of 100, and S′ in the order of 10, the number of possible combinations in s is S!/(S′!(S−S′)!) ∼ 10¹³, making it computationally unfeasible to find the exact minimum by calculating the error for each combination. We have therefore employed the following minimization algorithm. We start by choosing equidistant spectral points in s, which corresponds to an under-sampling of the spectral resolution, and we choose the s with minimum distance and centered in the spectral range. The minimization is then done by a random walk, choosing new spectral points s′ with s′_i = s_i + δ_i using a random value δ_i having a uniform distribution in the range s_i−1 − s_i < 2δ_i < s_i+1 − s_i, i.e. covering half the distance to the adjacent spectral points, such that the ascending order is conserved. The new spectral points s′ are kept for subsequent steps only if the error improved. We limit the number of steps in the random walk, to 200 in the results shown here. We also investigated other minimization algorithms leading generally to similar or inferior results in terms of convergency and computational complexity. The error ε of the reconstruction for equidistant spectral points, ε^eq, and for random walk optimized points, ε^rw, is shown in Fig. 2(b). The error reduces with increasing number of measured spectral points S′ approximately as $1/\sqrt{S^{\prime }}$.

The evolution of ε during the random walk minimization is shown in the inset of Fig. 2(b) for S′ = S_max = 8. Most of the improvement is realized within 10³ steps. The standard deviation of ε over 10 independent walks is generally below 0.2ε, so that we can conclude that the minimization algorithm is sufficiently stable. Interestingly, the error for equidistant spectral points is typically only a factor of two larger than after random walk minimization, showing that the exact choice of the spectral points is in general not important for the retrieval.

The number of important components S′_max for the reconstruction is then determined taking into account the modification of the noise by multiplying the linear fit of the noise components in σ (see solid line in Fig. 2) by the factor $\eta =\left({\epsilon }^{\text{rw}}/{\epsilon }^{\text{eq}}\right)/\sqrt{S/S^{\prime }}$, and is shown as dashed line in Fig. 2(a).

Spectral components which are not represented in the SVD basis, i.e. which are within the noise of the small spatial dataset, are not retrieved by the sparse sampling method.

3.2. Reconstruction of the spectra of a sparsely sampled hyperspectral image

Having determined s, a hyperspectral CARS intensity dataset B is acquired at the S′ spectral points of s over a large spatial region with a number of spatial points P ≫ P′. We used here the image shown in Fig. 1, for which P = 95120, and we chose S′ = S_max = 8, for which we find S′_max = 4 (see Fig. 2). We format the data in B as a S′ × P matrix D̂ of $\sqrt{I_{\text{C}}}$ and reconstruct the data over all spectral points using D^rec = ũû⁻¹D̂.

To evaluate the retrieval, we compare in Fig. 3(a) the spectra of I_C at the positions marked in Fig. 1 as obtained from the reconstruction algorithm (dashed lines) and directly measured and denoised (solid lines). The reconstruction reproduces the measured spectra up to small deviations, which are more important for the position #1 corresponding to a lipid droplet in the cell. We have verified that using S_max instead of S′_max components in the reconstruction results in a significantly increased noise in D^rec, as illustrated in Fig. 3(a) (gray line). The full datasets of the denoised and reconstructed $\sqrt{I_{\text{C}}}$ are shown in panel b) and c) of the hyperspectral movie Media 1, respectively.

 

Fig. 3 (a) Comparison of the measured (solid lines) and reconstructed (dashed line) CARS intensity spectra at the positions indicated in Fig. 1. The light gray spectrum is obtained using S′_max = S_max. The blue dots represent the S′ = 8 wavenumbers given by s used in the reconstruction of the hyperspectral image. A hyperspectral movie of the measured and reconstructued $\sqrt{I_{\text{C}}}$ is shown in panels b) and c) of Media 1, respectively.(b) Corresponding retrieved $\Im (\overline{\tilde{\chi }})$ using the PCKK method. A hyperspectral movie of the measured and retrieved $\Im (\overline{\tilde{\chi }})$ is shown in panels b and c) of Media 2, respectively.

Download Full Size | PDF

To determine the CARS susceptibility, we have measured the CARS intensity I_ref in a non-resonant medium (the coverslip glass) to calculate the CARS ratio Ī_C = I_C/I_ref, which is then used to retrieve the imaginary part of the complex susceptibility $\Im (\overline{\tilde{\chi }})$ using the phase-corrected Kramers-Kronig (PCKK) method described in [1]. In Fig. 3 we compare the resulting $\Im (\overline{\tilde{\chi }})$ at the positions marked in Fig. 1 as obtained from the reconstruction (dashed lines) and directly measured (solid lines). We see that the reconstructed $\Im (\overline{\tilde{\chi }})$ reproduces the directly measured spectra up to small deviations. The full datasets of the retrieved $\Im (\overline{\tilde{\chi }})$ using unfiltered, SVD denoised or reconstructed Ī_C are shown in panels a), b) and c) of the hyperspectral movie Media 2, respectively.

We reconstruct $\sqrt{I_{\text{C}}}$ because it has a white noise and spectrally local to the measured I_C. One could also reconstruct $\text{log}(\overline{\tilde{\chi }})$, since it is linearly related to log(Īχ_C) via PCKK, but it is spectrally non-local and not of white noise.

To investigate the error of the reconstruction in more detail, we consider the spatial distribution of the spectral error $E_{{\text{S}}_{\text{p}}}=\sqrt{P}\left({\sum }_{j=1}^S{\left({\mathbf{D}}_{j,p}-{\mathbf{D}}_{j,p}^{\text{rec}}\right)}^2\right)/{\Vert \mathbf{D}\Vert }_{\text{F}}$ where D is the noise filtered data taken at full spectral and spatial resolution. Figure 4 shows the spatial dependence of the spectral error for $\sqrt{I_{\text{C}}}$ and $\Im (\overline{\tilde{\chi }})$. The largest deviations can be observed at the lipid droplets.

 

Fig. 4 Spectral error E_S of the $\sqrt{I_{\text{C}}}$ and $\Im (\overline{\tilde{\chi }})$ due to reconstruction from sparse sample data a linear grayscales as shown. The scale bar represents 5 μm.

Download Full Size | PDF

The demonstrated sparse sampling method corresponds to an improvement in the acquisition speed by a factor (S_max/S + ρ)⁻¹ where ρ = P′/P is the fraction of spatial positions in dataset A compared to B, For the example shown here we have S_max = 8, S = 261, ρ = 0.01 resulting in a 25 fold increase in acquisition speed. In our analysis it is required that ρP ≫ S for a correct SVD-based noise filtering procedure. In typical situations where samples of a given chemical variety are imaged over a large set of replica, such as in high throughput screening, dataset A needs to be acquired only once, such that the speedup is effectively only limited by S/S′, which is 33 in the present case. In absolute terms it is typically sufficient to measure at 5–10 spectral points to reconstruct most of the chemical information obtainable in a full spectral scan.

4. Conclusion

We have developed and demonstrated a sparse sampling method to reduce the acquisition time in CARS hyperspectral imaging. The method employs an adapted factorization of a large hyperspectral data matrix (S × P) into two smaller matrices of (S × P′) and (S′ × P), with S′ ≪ S and P′ ≪ P. We demonstrated that the acquisition time for a human osteosarcoma U2OS cell was reduced by a factor of 25 without significant loss of information. This method, combined with state of the art CARS microscopes, is expected to enable real-time chemical imaging and high-throughput high-content label-free microscopy. The method is suited also for other coherent vibrational microscopy techniques such as stimulated Raman scattering, and in general for hyperspectral imaging techniques with sequential spectral acquisition.