Open Access
Add to BibSonomyAbstract
We present a measurement and reconstruction method for laser-scanning microscopy based on compressed sensing, which enables significantly higher frame rates and reduced photobleaching. The image reconstruction accuracy is ensured by including a model of the physical imaging process into the compressed sensing reconstruction procedure. We demonstrate its applicability to unmodified commercial confocal fluorescence microscopy systems and for Raman imaging, showing a potential data reduction of 10–15 times, which directly leads to improvements in acquisition speed, or reduction of photobleaching, without significant loss of spatial resolution. Furthermore, the reconstruction model is also robust to noise, and effective for low-light applications. This method has promising applications for all imaging modalities based on laser-scanning acquisition, including fluorescence, Raman, and nonlinear microscopy.
© 2016 Optical Society of America
1. Introduction
Recent applications in microscopy are requiring higher recording capabilities and speed, as in the case of time-lapse imaging, three-dimensional measurements, or high-throughput experiments, while improvements in resolution often imply noisier signals due to the lower amount of photons available per pixel. These constraints are pushing the development of faster acquisition approaches.
In the past years, compressed sensing (CS) has demonstrated potential to increase signal acquisition speeds by enabling efficient data retrieval from significantly fewer coefficients than the original data points with very low to no data loss [1, 2], provided that the signal is represented in a suitable basis [3, 4]. It has enabled a reduction in the amount of required measurements in various fields, such as magnetic resonance imaging [5] or computerized tomography [6]. This approach has also been extensively used in optics and imaging [7, 8], where images are usually highly compressible signals due to their low frequency content. Several different implementations of microscopes for CS were developed [9], leading to numerous applications, such as imaging with a single-point detector [10], sparse excitation for wide-field fluorescence [11], hyperspectral imaging [12], or spectroscopy [13]. However, the applications mentioned above usually require developing specific hardware to enable the acquisition of the linear combinations needed in CS. On the other hand, if the imaging process itself is inherently sparse, direct use of a CS reconstruction approach fits well with the nature of the acquired data. Examples of this include holography, by employing the underlying 3D information in the recording [14] or by relying on the properties of the reconstruction model [15, 16], or in stochastic optical reconstruction microscopy [17], where each frame contains very sparse information.
However, laser-scanning microscopy (LSM), which covers widely used methods such as confocal and two-photon fluorescence, is a point-based measurement approach with sequential acquisition which precludes the use of an encoding device to structure the light and generate the signal combination as required by CS. Furthermore, several CS implementations create the masks at the detection level by discarding light, which is not suitable for low-light applications, such as Raman spectroscopic imaging where additionally minimizing the incident power on the sample is also a goal, although more efficient masks based on the spectral properties of the sample can also be employed [18]. Other attempts include modulation of the excitation pattern, such as by employing sparse beamlets to generate some level of confocality in the measurement [19], but the loss of excitation power precludes this solution to be used for applications needing high local power, such as in nonlinear microscopy.
Here we propose a CS approach which does not require any encoding device, and can be implemented directly in existing LSM systems without any modification. The reduced acquisition rate is based on recording only a subset of the image through the targeting of spatial locations by addressing the scanning system accordingly. The reconstruction is then performed through an optimization constrained by minimization of the total variation (TV) of the reconstructed image. This approach alone is not directly suitable for use in a CS framework, but by additionally considering the point-spread function (PSF) of the system, the conditioning of the inversion problem can be significantly improved, allowing an image reconstruction without significant loss of spatial resolution. A similar approach based on subset measurements has been previously proposed for electron microscopy [20], where the reconstruction is based on a priori information of the sample structure, as opposed to the modeling of the system imaging function employed here.
This approach presents the advantage of allowing faster imaging based on the reduced sampling, and furthermore retains the tightly focused excitation, enabling its use for nonlinear microscopy where strong local power is required, and does not involve any loss of light.
We first test the performance of this spatial approach with state-of-the-art CS measurement schemes in simulations, by comparing their respective reconstruction fidelity, resolution, and robustness to noise. We then demonstrate its application to confocal fluorescence, acquired with a commercial instrument, and assess its ability to accurately retrieve fine cellular structures under significantly reduced sampling. The method is also demonstrated in hyperspectral Raman imaging, which requires a three-dimensional inversion due to the presence of the additional spectral dimension. In this context it allows, under significantly reduced sampling conditions, the reconstruction of both a Raman spectroscopic image and full spectral information at each location of the sample. These measurements also enable the possibility of improving the signal-to-noise of the data through regularization of the reconstructed signals, and confirm the ability to employ these measurements for quantitative chemometrics.
2. Theory
2.1 Compressed sensing principle
The fundamental principle of compressed sensing is the capability of exactly representing a signal x of size N with a significantly reduced amount M of coefficients, also denoted as measurement b, provided that the coefficients are expressed in a basis enabling a sparse representation. In its simplest form, a standard CS problem can be described as
where A is the measurement matrix of size M × N, which contains the coefficients of the linear combinations of the values of x, defining the measurements b. Under the assumption that x is sparse to some degree, Eq. (1) can be solved as a linear problem such as through convex minimization under a given set of constraints. In the case of imaging, it is common to consider TV as an optimization constraint, which is known to have excellent regularization properties, while preserving the sharp features of the image.The exactness of the reconstruction is then essentially determined by the basis defined by the vectors A = [a0;…;aM], which must be as close to orthogonal as possible. As noted above, the generation of the measurement patterns in microscopy usually requires the insertion of a device modulating the excitation or detected light pattern by discarding light, as illustrated in Fig. 1(a). Popular patterns include binary masks [10, 11] such as Hadamard functions generated with digital mirror devices or more complicated patterns [4, 21] that require spatial light modulators.
Fig. 1 (a) Measurement principle for a state-of-the-art compressed sensing scheme, where the light is modulated by a programmable device, generating excitation pattern to create the measurement matrix A through sequential measurements. (b) Classical laser-scanning microscope, where the inclusion of the point-spread function into the reconstruction model spreads the local measurements.
As the CS reconstruction principle inherently relies on the linear combinations of several samples for generating the measurement values, it has essentially been used for applications with continuous dimensions, such as wide-field microscopy. This approach cannot be translated to other widely used imaging methods based on scanning and point detection, such as confocal or nonlinear microscopy, since the sequential detection of precise locations prevents the measurement of linear combinations of signals from different regions in the sample. Furthermore, the data reduction procedure would imply that some locations in the sample are simply not measured, which is not consistent with the reconstruction approach.
The problem expressed in Eq. (1) in the context of missing spatial locations has however also been addressed with TV regularization in the context of image restoration and zooming [22, 23]. This approach has been shown to minimize discontinuities and oscillations occurring with other classical approaches, along with good properties in regards to noise robustness [24]. These models usually include a smoothing operation applied to the original fully sampled signal [25].
2.2 Sub-sampling for LSM
Reconstructing data obtained from LSM with missing samples is similar to using TV image restoration, where the system point-spread function (PSF), which characterizes the physical limitations of the instrumental setup, can be employed as the mathematical smoothing function. This principle is illustrated in Fig. 1(b), where the spread of the signal caused by the PSF implies that each measured point leads to a linear combination of values in the vicinity, as shown in Fig. 1(b2). While a straightforward sub-sampling implies that only a few of the locations are non-zero, leading to a sparse measurement matrix, the inclusion of the smoothing operator (i.e. the PSF) in the model spreads the signal contributions, significantly increasing the amount of non-zero coefficients. The measurement matrix A can then be made more homogeneous by employing a Fourier representation of the vertical vectors instead of the original measurements, providing a full pseudo-random measurement matrix that can be employed for CS reconstruction.
One should note that this approach differs from the standard CS procedure, where a physical device generates the vectors of the matrix A, which is then employed for the reconstruction of the signal [Fig. 1(a)]. Instead, the measurements here consist of an under-sampled version of x, which is reconstructed from a measurement matrix A generated off-line based on an imaging model of the system. This procedure is illustrated in Fig. 1(b1), where the measurements of different locations in a laser-scanning setup are used to generate the measurements, which are directly detected without modification. This implies that there is no loss of power neither at illumination nor at detection level.
The inclusion of the imaging model in the reconstruction means that each measurement is the result of a linear combination of the values in the vicinity. With this inherent spread in the measurement, the inversion problem can be described as
where I(x,y) is the fully sampled image 2D intensity values, H(x,y) is the 2D PSF function, * is the convolution operator, and vec[] is the vectorization operator. While Eq. (2) represents the 2D image data in a form related to the physical variables, it is possible to rewrite the problem in the form of Eq. (1), as all of the operations are linear.The convolution in Eq. (2) implies that the imaging process is modeled by smoothing the signal leading to the measurements b, and that during inversion, it is required to apply a deconvolution on the measurements to retrieve the candidate signal x at each iteration of the reconstruction procedure. We employ here a truncated division in the frequency domain, which can be described as
where Î, Ĥ are the Fourier transforms of I and H respectively, F−1 is the inverse Fourier transform, ω = (ωx,ωy) is the spatial frequency, Ĩ is the deconvolved image, and ε is a small positive value, chosen to truncate the high frequencies and reduce the noise amplification known to occur with deconvolution performed through direct frequency division [26]. By comparing Eqs. (1) and (2), one can see that the result of the reconstruction should be I∗H, corresponding to the classical diffraction-limited image without missing samples.There are several deconvolution approaches which have been proposed in the literature, which are more refined and potentially provide better results than the deconvolution described in Eq. (3) [27]. One should also note that the approach of applying deconvolution in compressed sensing has also been recently discussed in theory [28], and applied in practice to wide-field microscopy techniques [29]. Regularization by minimization of TV is also often used in deconvolution procedures [30, 31]. The choice here of a simpler deconvolution procedure is motivated by its fast computation, as this step has to be performed at each iteration of the CS reconstruction process. Furthermore, it can be assumed that the noise amplification induced by the frequency division can be compensated by the regularization provided by TV. We noticed in practice that while the value of ε does not have a major influence on the results within a reasonable range, too small values can cause the algorithm to diverge due to the noise amplification.
The problem described in Eq. (2) is solved by employing a total variation augmented Lagrangian solver (TVAL3 [32, 33]) implemented as a MATLAB code, which was originally designed to solve CS problems by minimizing either the l1- or l2-norm of the TV for both 1D and 2D signals. We modified the code to implement the generation of measurement matrices corresponding to the model described in Eqs. (2) and (3), and to support 3D signals in the case of Raman hyperspectral data.
3. Results and discussion
3.1 Simulations
We use simulated images to compare our approach with standard measurement matrices employed in CS. We employ two types of simulated images (256 × 256 pixels); the first image [Fig. 2(a1)] contains both smooth and sharp features enabling the estimation of reconstruction accuracy, while the second image [Fig. 2(a2)] simulates beads which provide information about spatial resolution. The imaging process is simulated by the convolution with a Gaussian smoothing kernel (full width at half maximum (FWHM) of 5 pixels), with the addition of Poisson photon statistics, for a relatively small amount of photons in the first image (maximum intensity of 1000 photons), and a total of 104 photons per bead in the second image.
Fig. 2 Comparison of reconstruction of a simulated image with different sub-sampling approaches. (a1) Simulated image with both sharp and smooth features and (a2) simulated beads, with a FWHM of 5 pixels. (b) Differences between the reconstructions (25%) and the original image A1. (c1) RMSE of the different reconstruction methods compared to the original image (a1); lines connecting the median errors are shown for the sake of visibility. (c2) RMSE of the sub. PSF incl. method for different PSF sizes and sub-sampling ratios. (d) Reconstruction accuracy for the case of simulated beads, measured as (d1) RMSE compared to original image and (d2) FWHM of the reconstructed beads (N = 20).
The images are sub-sampled, either by randomly removing coefficients in the spectral domain after a Fourier or Hadamard transform (rates ranging from 5% to 50%), or by homogeneously under-sampling the image for spatial approaches (corresponding rates of 1/4, 1/9 or 1/16), either for direct spatial sub-sampling (Spatial sub.) or by including the imaging model in the reconstruction (Sub. PSF incl.). The reconstruction performance is then estimated by calculating the root-mean-square error (RMSE) compared to the noise-free signal, where 10 realizations are considered for spectral methods which rely on random sub-sampling. The RMSE here accounts both for reconstruction errors and the effect of regularization, which minimizes the error compared to the noise-free signal by reducing the noise.
The resolution is assessed by extracting 20 beads profiles, on which a Gaussian curve is fitted to estimate the FWHM of the reconstructed signal. The reconstruction parameters for each method, and in particular regularization parameters are adjusted to minimize the variation in FWHM resolution.
The differences between the reconstruction and the original (noise-free) image displayed in Fig. 2(b) show that, while spectral methods based on Fourier and Hadamard bases have a rather constant error pattern throughout the image, spatial approaches have large differences mostly concentrated around edges. The overall RMSE values are smaller for spatial approaches [Fig. 2(c1)], in particular at lower sub-sampling rates, indicating that spatial methods tend to reduce noise more strongly, which is also consistent with the difference pattern showing less variation in the background regions. Spatial approaches also have a faster convergence rate and require approximately half or less of the amount of iterations compared to spectral methods.
In the case of the image of simulated beads, the spatial approach leads to very large errors, showing that it cannot accurately reconstruct sharp objects. On the other hand, including the PSF in the spatial model leads to results comparable with spectral methods in terms of RMSE [Fig. 2(d1)].
Similarly, the estimation of the reconstruction resolution shows a large degradation of performance in the spatial-only case. The other methods, both spectral and spatial with PSF included, can recover the original FWHM resolution at 25% sampling [Fig. 2(d2)]. The resolution then degrades at lower sub-sampling rates, in particular for the spatial approach including the PSF, where it tends to provide larger FWHM with smaller variation, and a lower RMSE.
These results show as expected that the performance of CS is dependent on the sub-sampling rate, and significantly degrades when too few samples are used. However, we can also expect in the case where the imaging model is included that the performance also depends on the extent of the PSF. To estimate the relation between the sub-sampling rate and the PSF size, we calculated the RMSE of the reconstruction of the simulated image without noise for PSFs between 3 and 11 pixels FWHM and for various sampling rates [Fig. 2(c2)]. There are clearly two regimes in the reconstruction performance, first where the RMSE stays below a certain level at sampling rates finer than the PSF size, and then increases exponentially for rates approximately at or below the PSF size, though the transition between these two regimes does itself depend somewhat on the extent of the PSF. By considering that the PSF FWHM approximately corresponds to twice the resolution, this indicates that our approach can be employed down to sampling rates two times smaller than the Nyquist frequency without significant effects on the error.
3.2 Beads measurements and PSF estimation
As mentioned in Section 2.2, our approach relies on the knowledge of the PSF to perform an accurate reconstruction based on a sub-sampled version of the image. We estimate the PSF of our imaging system experimentally by employing 100 nm fluorescent beads (exc 540 nm / em 560 nm) as point sources. The beads are deposited on a glass substrate and observed with both the confocal and Raman systems presented below.
In both cases, the beads are measured at the maximum available resolution, by employing the maximum optical zoom (10x) on the confocal microscope (yielding 20.72 nm/px), or by using the largest possible camera resolution on the Raman system (pixel size Δx = 6.5 μm, yielding 104 nm/px). An image of the beads measured by the Raman system is then generated by averaging the 900–1000 cm−1 region, located at the center of the fluorescence emission.
To generate the experimental point spread functions, 20–30 beads are manually extracted from the measurements. The images are then automatically registered by employing the StackReg package [34] freely available for the ImageJ program. The PSF is then obtained by averaging all the bead images, re-sampling at the pixel size at which measurements are performed, and scaling to compensate for the emission wavelength used in the reconstruction algorithm. A typical field of view of the fluorescent beads measured with the confocal microscope is shown in Fig. 3(a).
Fig. 3 (a) Typical field of view of deposited fluorescent beads. (b) Beads distribution in area and intensity, showing a cluster with normal distribution at the bottom-left, including ∼74% of the beads.
To ensure that the measurements correspond to single beads, images were analyzed to extract the mean intensity and area of each bead by employing the Analyze Particles function of ImageJ. It is possible to identify a normal distribution both in the mean intensity (μ = 1010, σ = 120) and area (μ = 22.7, σ = 3.5) as shown in Fig. 3(b), that can be assumed to correspond to single beads. The spread then corresponds to aggregates of beads, either in a linear fashion (larger area, stronger intensity) or rather far apart (same intensity, larger area). Extracting the beads contained in the 3σ elliptical region shows that 74.1% of the measured beads correspond to single ones. To further ensure that no weaker features were missed, measurements were also performed at twice the usual power, which showed no other emitting particles.
3.3 Confocal microscopy
We first demonstrate the use of our proposed method with measurements performed on mouse embryonic fibroblast (MEF) cells stained with the ER-tracker Green and Hoescht 33342 according to the manufacturer’s instructions, which bonds respectively to the endoplasmic reticulum and the nucleus. Cells are then fixed with a 4% paraformaldehyde solution, and observed with a commercial confocal microscope (Olympus FluoView FV10i-DOC). Briefly, the dyes are sequentially excited with 405 and 473 nm lasers, respectively. The fluorescence is collected with a 60x microscope objective (1.35 NA, oil immersion), with an optimal optical zoom of 1.2x, yielding a pixel size of 172.6 nm, focused onto a pinhole (2 Airy units), and spectrally filtered from excitation and other emission bands, before being detected with a photomultiplier tube, at respectively 455 and 511 nm. Cells are observed at different sampling rates, and the reconstructed images are then compared to a measurement performed at full resolution.
We first assess the imaging properties of the microscope by measuring the emission of single fluorescent beads (see Section 3.2). The PSF appears elongated vertically [Fig. 4(a1)], which can be attributed to the linear polarization of the excitation laser, as well as aberrations in the system which induce small distortions. It is possible to see that at the pixel size employed for measurements (172.6 nm), which is slightly smaller than the theoretical resolution at this wavelength (256 nm), the experimental PSF spreads across approximately 4-5 pixels [Fig. 4(a2)].
Fig. 4 Comparison of confocal fluorescence measurements of MEF cells with under-sampling and reconstruction. (a1) Point spread function (square-root dynamics) of the confocal microscope, measured with 100 nm fluorescent beads, which is (a2) re-sampled at the resolution employed for cell samples, showing a spread across approximately 5 pixels. (b) Cells stained in the nucleus (cyan) and ER (green), measured at full resolution, with an inset of ER staining (in gray) showing fine structures. Sub-sampled measurement at 25% and 6.25%, then reconstructed to the full resolution with (c) a straightforward spatial approach, and (d) a method which includes the imaging process modeled by the experimental PSF shown in (a2). (e) Inset obtained through nearest-neighbor interpolation.
The MEF cells are stained both with a nuclear dye (Hoescht, cyan) and a dye bonded to the endoplasmic reticulum (ER, green). A typical cell at full resolution (1024 × 1024) is shown in Fig. 4(b), where several nuclei with multiple nucleoli can be identified, and the ER being distributed throughout the cytoplasm, with some fine structures further away from the nuclei. An inset from another cell with only ER staining is also shown in gray, where finer structures can be identified.
Two additional measurements are taken with the same field of view but at lower resolutions of 512 × 512 and 256 × 256, respectively corresponding to sub-sampling rates of 25% and 6.25%, with all other acquisition parameters kept identical. The images are then reconstructed at the original resolution through the proposed method, either through a simple spatial approach [Fig. 4(c)], or by including the PSF in the reconstruction process [Fig. 4(d)]. As expected from the ill-conditioned matrix used in the simple spatial approach, it provides results which are less sharp than the original image, in particular in the case of a lower sub-sampling, with the images appearing strongly blurred.
On the other hand, the approach including the PSF provides good visual results at both sub-sampling rates, with sharpness comparable to the original image. While there is no visible difference between the full resolution and 25% sub-sampled reconstructed images, it is possible to identify a slight loss of quality in the 6.25% one, especially at the edges of some fine structures, although all features are still resolved compared to the full resolution image. The inset is also shown at the same resolution, obtained through nearest-neighbor interpolation [Fig. 4(e)], hence corresponding to the originally measured points, where the loss of spatial resolution can be identified because of the missing samples, especially at 6.25%. Interestingly, it is also possible to see that the reconstructions in Fig. 4(d) seem to contain less noise than the fully resolved measurement, as it can be seen especially in the constant regions within the nucleus. This can also be attributed to the TV regularization.
To demonstrate another practical benefit of this approach, we also consider the reduction of photo-bleaching by sub-sampling. The bleaching is estimated in situ by continuously measuring the same field of view at the maximum frame rate allowed by a given sampling rate, and extracting the average intensity of the brightest spots in the ER image over time [Fig. 5(a)]. The bleaching rate is then estimated by fitting an inverse exponential function to the intensity curves [Fig. 5(b)]. Logically, sub-sampled images can be acquired at a much higher frame rate [Fig. 5(c)], with similar bleaching rates.
Fig. 5 (a) Fluorescence image of ER staining, where the intensity of the brightest spots (see red squares) are extracted. (b) Average fluorescence intensity showing the exponential decrease (lines, fitted) due to bleaching at full resolution. (c) Average intensity (n = 10) at different sub-sampling rates for a constant amount of images, showing the reduced bleaching when sub-sampling.
To further assess the quality and resolution performances of the reconstructed images in practical conditions, we also consider MEF cells tagged with the MitoTracker Deep Red dye (excitation 635 nm, emission 665 nm), which stains mitochondria, of a size close to the typical resolution limit at this wavelength. An optical zoom of 1x is employed in this case, yielding a pixel size of 207.2 nm. In addition to the sub-sampling rates of 1 in 2 pixels and 1 in 4 pixels, we also adjust the optical zoom of the instruments (1.4x) to reach an effective rate of 1 out of 3, or 11.1%, by changing the field of view while keeping the image size constant (256 px). Reconstruction is only performed with inclusion of the PSF, as pure spatial approaches were shown to yield poor results [Fig. 4(c)]. The resulting images are shown in Figs. 6(c) and 6(d), which can be compared with the full resolution image in Fig. 6(a), with performances comparable with previous measurements, where 25% and 11% results provide similar imaging quality, while several features are not fully recovered in the 6% reconstruction. It is interesting to note that in this case, the reconstruction quality of the blue channel significantly degrades at lower rates. This is due to the larger scanning step, here optimized for red wavelengths, which implies that at 6%, the sampling period becomes greater than the PSF width in the blue. As also shown in the simulations above (see Section 3.1), this essentially implies that the sampling period can be extended up to but not exceeding the approximate width of the PSF before the image shows features that are missing when compared to a fully sampled measurement.
Fig. 6 (a–d) MEF cells measured in fluorescence with nuclear and mitochondrial stains and acquired respectively at 100%, 25%, 11.1% and 6.25% sub-sampling rates and reconstructed with the method including the PSF. A small portion of the image (gray insets, see dashed rectangle) is shown to highlight the resolution performance of each reconstruction. Profiles are extracted on small features (red and green dashed lines), and normal curves are fitted to quantitatively assess the reconstruction resolution. (e) Average width (n = 20) of small features at different locations in the image, showing the gradual loss of resolution with decreasing sub-sampling rates, which can be compared with the (f) averaged images and the expected PSF based on beads measurements.
Insets in gray show a region in the image with a low density of emitters [dashed square in Fig. 6(a)], where the reconstruction performance can be further assessed. It is possible to see that there is a gradual loss of contrast with a decreasing sub-sampling rate, making it more difficult to identify weak features. While the general structure in the inset is maintained up to 11%, several artifacts start to appear at 6%. We can assume that at this rate, the spacing of the measured samples (1 out of 4) becomes too large compared to the extent of the PSF to enable an accurate reconstruction.
To estimate the potential resolution degradation of the reconstruction, we extract profiles from small emitters [red and green dashed lines in Fig. 6], and evaluate the FWHM of Gaussian curves fitted on the measured points, with a fit shown for each image. The procedure is repeated for several emitters in the image (n = 20), to assess repeatability [Fig. 6(e)]. The selected emitters are also averaged together to provide a visual estimation of the resolution [Fig. 6(f)], where the gradual spread can also be identified, and compared with the expected PSF size at this wavelength. It is possible to see that the PSF width at 100% is slightly larger and more circular than the expected PSF shape, which can be attributed to the 3D distribution of emitters, and to the larger size of mitochondria compared to beads.
The measured resolution of the full resolution image is around 430 nm, consistent with the PSF measured from beads. The resolution of the reconstructed profiles is around 490 nm for both 25% and 11%, implying a slight increase compared to the original image. Interestingly, decreasing the amount of samples does not significantly increase the width of the reconstructed emission profiles, although additional loss of contrast is evident. The degradation at 6% can be seen by the larger FWHM values with a wide distribution, along with artifacts in the image. This degradation is comparable with the results obtained in simulations [Fig. 2(d2)].
The results in Figs. 3 and 4 show that fine structures require a higher sampling to be accurately reconstructed, as can be seen for example by comparing the 25% sub-sampling in Figs. 3(d1) and 4(b). This is also confirmed in simulations, where the error is higher at the edges of the pattern. This implies that while the homogenous sub-sampling employed here is easier to implement, particularly in the case of commercial instruments, our approach could potentially be made more efficient by employing a selective sampling focusing on regions with a higher frequency content, as was proposed for instance for Raman confocal spectroscopy [35].
3.4 Raman microspectroscopy
The proposed CS method can be applied to any imaging method based on a laser-scanning acquisition, and we demonstrate its application here to Raman spectroscopic imaging, which could benefit dramatically from an increase in acquisition frame rate. We employ here a laboratory-developed system based on a slit-scanning configuration, described in previous publications [36, 37]. Cells are excited by a continuous wave laser at 532 nm, focused with a 60x microscope objective (NA 1.27, water immersion). The Raman emission is separated from the excitation with a long-pass dichroic mirror. The Raman emission line is then focused back on the slit of a Czerny-Turner spectrometer (500 mm focal length), in which the light is dispersed by a grating (300 lp/mm). The spectrum of the whole line is then detected in parallel by a 2D water-cooled low-noise scientific CMOS detector (Δx = 13 μm, binning 2 × 2), with an exposure time in the range of seconds. A background spectrum, measured without a sample and at the same sampling conditions as the measurements, has been subtracted from all shown images.
There are several differences between confocal microscopy and Raman imaging, so that while the reconstruction strategy is similar, it requires some adjustments to optimize the reconstruction. First, since Raman spectroscopy measures the vibrational spectrum at a given location, the recorded data is a hyperspectral 3D stack instead of a 2D image. This requires the generalization of the reconstruction algorithm to take into account the additional dimension, so that regularization is also applied in the spectral dimension.
Furthermore, as the system employed to measure the Raman signals is based on a slit-scanning configuration, it implies that the spectra of a whole line of the object are measured in parallel on a 2D detector. While this means that the measurements can be performed faster, it also means that there is only one scanning dimension on which sub-sampling can be performed, significantly reducing the potential data reduction that can be provided by CS. To compensate for this shortcoming, we also take advantage of the noise reduction capability of TV regularization, by decreasing the exposure time for the data sets employed for CS reconstruction, so that fully sampled measurements are acquired with an exposure of 6 seconds per line, while sub-sampled ones are performed at half the exposure, at 3 seconds.
The PSF of the Raman microscope is estimated by using the procedure described before, where the signal of the fluorescent beads is resolved by the spectrometer. The PSF has a vertical elongation which can again be attributed to the linear polarization of the excitation laser [Fig. 7(a1)], and the PSF re-sampled at the scanning step extends across 4–7 pixels. It is also possible to identify a large horizontal spread, which is in line with the direction of the spectrometer slit, so that out-of-focus light is not rejected in this axis.
Fig. 7 Raman measurements on MEF cells for fully sampled data compared to sub-sampled data with reconstruction. (a1) Point spread function of the Raman microscope, measured with 100 nm fluorescent beads (square-root dynamics), which is (a2) re-sampled at the resolution employed for measurements, showing a spread across approximately 4–7 pixels. Cells are represented by a C-H stretching band (2925 cm−1, green) and lipids (2853 cm−1, red), (b) measured at full resolution. Sub-sampled measurement, then reconstructed with (c1,2) the method including the imaging process modeled by the experimental PSF shown in (a2).
The sub-sampled measurements are performed by changing the steps on the descanned mirror, in order to measure 1 line out of 3 or 5 samples, leading to 33% or 20% data acquisition, where fixed MEF cells are observed, as in the previous experiments. The PSF model is sequentially applied to all the spectral 2D images. After reconstruction, a C-H stretching band (2925 cm−1) and a lipids band (2853 cm−1) are extracted by averaging 3 spectral images around the bands, and are represented in false colors [Fig. 7(b)]. The cell structure can be identified in green (C-H stretch) and in particular the nucleus, while lipids are highlighted in red, showing structures in yellow.
As previously, the approach including the PSF [Fig. 7(c)] provides a resolution comparable to the original fully sampled image, although very fine structures appear slightly blurred. At lower sub-sampling rates (20%), it is possible to identify some errors in the reconstruction, with some details appearing as horizontally elongated. The fact that the image reconstruction performance appears less good than in the case of the fluorescence imaging can be attributed first to the sub-sampling strategy, which is limited to only one dimension and thus less favorable for efficient reconstruction at low sampling rates, but also to the signal itself, as the Raman images are significantly noisier than the fluorescence ones, due to the low Raman cross-section of biological samples, so that the reconstruction quality is affected. However, this comes at the benefit of acquisition speed, which, combined with the 2 times shorter exposure time and the under-sampling rate of 1 line in 5, corresponds to an increase in the Raman imaging frame rate of 10 times, yielding an effective exposure of 0.6 second per line, or 1.5 ms per pixel.
Nevertheless, the reconstructed images are less noisy than the fully sampled image, as can be estimated by computing the image signal-to-noise ratio (SNR) in different regions, defined as the ratio of the signal mean to its standard deviation. The SNR is estimated in 10 × 10 regions, selected as a nucleolus and in the cytoplasm, on different Raman images, at 2925 cm−1 (lipid bands) and 1330 cm−1, in which the nucleus and in particular the nucleoli are highlighted, which can mostly be attributed to adenine [38]. As shown in Table 1, the reconstructed images have an SNR between 1.5 and 2.5 times higher than the fully sampled measurement, despite the fact that they were measured with only half of the exposure. Interestingly, the lower sampling of 20% provides a higher SNR in some cases, but as shown in Fig. 7, this is at the cost of a loss in spatial resolution.
Table 1. SNR of the sub-sampled reconstructed images [Fig. 7] in different regions, compared to the fully sampled measurement.
The noise characteristics of the Raman images also imply that regularization parameters can be optimized to improve the reconstruction quality. One should also note that the three-dimensional reconstruction significantly contributes to the image quality, due to the TV regularization in the spectral dimension. 3D reconstructions typically provide much better image quality thanks to the inclusion of the spectral dimension than pure 2D reconstructions where the continuity of spectra is not considered. As the regularization is also applied in the spectral domain, the noise rejection of the CS process also significantly improves the SNR of the Raman spectra, as shown in Fig. 8(a). The spectra are also modified by the reconstruction process, where the noise reduction can potentially enable the detection of smaller Raman peaks. In order to assess the possibility of employing these spectra quantitatively, we apply k-means clustering on the hyperspectra, which has been shown to enable the separation of various spatial regions corresponding to cellular structures based on the distances between spectra [39]. K-means clustering and error estimation are performed in the statistics program R, with the packages cclust and fpc. To ensure reproducible results, clustering is performed 10 times per data set, and the result maximizing consistency, estimated by the average silhouette width [40] is selected.
Fig. 8 (a) Normalized Raman spectra of one location in the MEF cell cytoplasm, for the original image at full sampling (6 s exposure) and the reconstruction (3 s exposure) for different sub-sampling rates. Despite the smaller signal at half the exposure, the reconstructed spectra are less noisy than the fully sampled one. (b–d) K-means clustering images computed from the Raman spectra shown in Fig. 7, for (b) the fully sampled measurement, and (c–d) the reconstructed images, at respectively 33% and 20% sub-sampling, showing that the clusters are not significantly changed by the reconstruction process.
The fact of employing 10 separate classes, which is a high number, may force several morphological features of the cell to be divided into several classes. However, this helpfully amplifies small differences between the different image reconstructions. The resulting images are shown in Fig. 8, for the full measurement and reconstructed signals.
The clustering of both the fully sampled measurement [Fig. 7(b)] and of the sub-sampled reconstruction [Figs. 7(c) and 7(d)] present very consistent results, with the same structures highlighted by the different regions. This shows that the regularization which reduces the noise in the signals does not affect the spectral content, and that the spectra reconstructed from the CS procedure can be employed for quantitative analysis.
4. Conclusion
The approach proposed in this article enables compressed sensing for laser scanning microscopy. It makes possible the retrieval of the signal from a significantly reduced amount of samples, although it does not correspond to a standard CS scheme since the measurements are not generated by a specific pseudo-random mask like most CS applications in imaging. The measurement matrix generated here is based purely on the imaging properties of the system, so that it is not possible to choose specific illumination patterns to optimize the recording scheme. On the other hand, it can be used on existing laser scanning microscopes, as it does not require any hardware changes to physically generate the measurement patterns, and is only based on the ability of controlling the steps of the scanning stage. Since the reconstruction relies on modeling the imaging properties of the system by including the PSF in the measurement model to condition the inversion procedure, there exists a trade-off between the PSF size and the possible sub-sampling rate that can be achieved. The largest sampling period which can be allowed should be on the order of the PSF FWHM width, so that unmeasured locations can still be assumed to contribute to the measured data via the physical smoothing induced by the PSF. This however does not reduce the applicability of the method in practice, as microscope systems typically over-sample the PSF to a certain degree to prevent any loss of spatial resolution.
Because the reconstruction method only relies on an imaging model for a scanning-based measurement scheme, it can be applied to any type of microscopy, as shown for instance here for both confocal fluorescence and Raman microspectroscopy. As this approach does not require any encoding of the light pattern and still relies on a tightly focused light excitation, it implies that there is no loss of light power neither at excitation or detection, so that it can be applied both in case of low light applications, such as Raman imaging, or in applications requiring high peak power such as nonlinear microscopy, including two-photon fluorescence, second harmonic generation, or coherent anti-stokes spectroscopy, for example.
The increased acquisition rate provided by our approach would be especially beneficial for live sample imaging, where dynamics measurements are required. While fixed cells were employed here to enable the validation of the method with non-moving samples, the approach can be directly used with live cells or tissues.
It is also possible to employ this method in case of multi-dimensional microscopy, as shown in the case of Raman imaging, which involves both spatial and spectral dimensions. The approach can also be generalized to three-dimensional microscopy in a straightforward way, where the inclusion of a 3D PSF can extend the sub-sampling capabilities, and reach more efficient data reduction rates than demonstrated in this article. For applications such as Raman imaging, the additional dimensions such as the spectral information can contribute to the quality of the reconstruction through regularization.
The sub-sampling strategy enables a faster acquisition of the data, where we showed reductions up to 10–15 times in confocal fluorescence measurements and 5 times reduction in the sampling for Raman spectroscopic imaging with only minor alteration of the image quality. Furthermore, due to the regularization of the data provided by the reconstruction based on total variation minimization, the overall noise can be minimized. In this fashion, we showed a measurement speed up of 10x in Raman measurements, with a 20% data reduction and half the exposure time, which could still outperform standard sampling Raman imaging in terms of image SNR.
Funding
Japan Society for the Promotion of Science (JSPS) (Program for World-Leading Innovative R&D on Science and Technology (FIRST), World Premier International Research Center Initiative Funding Program (WPI), KAKENHI JP15K17473); Japan Science and Technology (PRESTO).
Acknowledgments
The authors would like to thank Dr. T. Saitoh (Tokushima University) and Dr. S. Akira (Osaka University) who donated the spontaneously immortalized wild-type mouse embryonic fibroblast cells, and E. Bostan (Ecole Polytechnique Fédérale de Lausanne) for fruitful discussion about imaging applications of total variation.
References and links
1. D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]
2. E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]
3. E. Candes and T. Tao, “Decoding by linear programming,” IEEE Trans. Inf. Theory 51(12), 4203–4215 (2005). [CrossRef]
4. E. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]
5. M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magn. Reson. Med. 58(6), 1182–1195 (2007). [CrossRef] [PubMed]
6. K. Choi, J. Wang, L. Zhu, T.-S. Suh, S. Boyd, and L. Xing, “Compressed sensing based cone-beam computed tomography reconstruction with a first-order method,” Med. Phys. 37(9), 5113–5125 (2010). [CrossRef] [PubMed]
7. J. Romberg, “Imaging via Compressive Sampling,” IEEE Signal Process. Mag. 25(2), 14–20 (2008). [CrossRef]
8. R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Opt. Eng. 50(7), 072601 (2011). [CrossRef]
9. M. A. Neifeld and J. Ke, “Optical architectures for compressive imaging,” Appl. Opt. 46(22), 5293–5303 (2007). [CrossRef] [PubMed]
10. M. Duarte, M. Davenport, D. Takhar, J. Laska, T. Sun, K. Kelly, and R. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag. 25(2), 83–91 (2008). [CrossRef]
11. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proc. Natl. Acad. Sci. U.S.A. 109(26), E1679–E1687 (2012). [CrossRef] [PubMed]
12. M. E. Gehm, R. John, D. J. Brady, R. M. Willett, and T. J. Schulz, “Single-shot compressive spectral imaging with a dual-disperser architecture,” Opt. Express 15(21), 14013–14027 (2007). [CrossRef] [PubMed]
13. B. M. Davis, A. J. Hemphill, D. C. Maltaş, M. A. Zipper, P. Wang, and D. Ben-Amotz, “Multivariate hyperspectral Raman imaging using compressive detection,” Anal. Chem. 83(13), 5086–5092 (2011). [CrossRef] [PubMed]
14. D. J. Brady, K. Choi, D. L. Marks, R. Horisaki, and S. Lim, “Compressive holography,” Opt. Express 17(15), 13040–13049 (2009). [CrossRef] [PubMed]
15. Y. Rivenson, A. Stern, and B. Javidi, “Compressive Fresnel holography,” J. Disp. Technol. 6(10), 506–509 (2010). [CrossRef]
16. A. Bourquard, N. Pavillon, E. Bostan, C. Depeursinge, and M. Unser, “A practical inverse-problem approach to digital holographic reconstruction,” Opt. Express 21(3), 3417–3433 (2013). [CrossRef] [PubMed]
17. L. Zhu, W. Zhang, D. Elnatan, and B. Huang, “Faster STORM using compressed sensing,” Nat. Methods 9(7), 721–723 (2012). [CrossRef] [PubMed]
18. D. S. Wilcox, G. T. Buzzard, B. J. Lucier, O. G. Rehrauer, P. Wang, and D. Ben-Amotz, “Digital compressive chemical quantitation and hyperspectral imaging,” Analyst (Lond.) 138(17), 4982–4990 (2013). [CrossRef] [PubMed]
19. Y. Wu, P. Ye, I. O. Mirza, G. R. Arce, and D. W. Prather, “Experimental demonstration of an optical-sectioning compressive sensing microscope (CSM),” Opt. Express 18(24), 24565–24578 (2010). [CrossRef] [PubMed]
20. A. Stevens, H. Yang, L. Carin, I. Arslan, and N. D. Browning, “The potential for Bayesian compressive sensing to significantly reduce electron dose in high-resolution STEM images,” Microscopy (Oxf.) 63(1), 41–51 (2014). [CrossRef] [PubMed]
21. Y. Tsaig and D. L. Donoho, “Extensions of compressed sensing,” Signal Process. 86(3), 549–571 (2006). [CrossRef]
22. A. Chambolle and P.-L. Lions, “Image recovery via total variation minimization and related problems,” Numer. Math. 76(2), 167–188 (1997). [CrossRef]
23. F. Guichard and F. Malgouyres, “Total variation based interpolation,” in Proc. Eusipco’98 (1998), pp. 1741–1744.
24. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vis. 20(1/2), 89–97 (2004). [CrossRef]
25. F. Malgouyres and F. Guichard, “Edge direction preserving image zooming: a mathematical and numerical analysis,” SIAM J. Numer. Anal. 39(1), 1–37 (2001). [CrossRef]
26. J. G. McNally, T. Karpova, J. Cooper, and J. A. Conchello, “Three-dimensional imaging by deconvolution microscopy,” Methods 19(3), 373–385 (1999). [CrossRef] [PubMed]
27. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Process. Mag. 23(3), 32–45 (2006). [CrossRef]
28. E. J. Candès and C. Fernandez-Granda, “Towards a mathematical theory of super-resolution,” Commun. Pure Appl. Math. 67(6), 906–956 (2014). [CrossRef]
29. S. Gazit, A. Szameit, Y. C. Eldar, and M. Segev, “Super-resolution and reconstruction of sparse sub-wavelength images,” Opt. Express 17(26), 23920–23946 (2009). [CrossRef] [PubMed]
30. N. Dey, L. Blanc-Feraud, C. Zimmer, P. Roux, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, “Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution,” Microsc. Res. Tech. 69(4), 260–266 (2006). [CrossRef] [PubMed]
31. S. D. Babacan, Z. Wang, M. Do, and G. Popescu, “Cell imaging beyond the diffraction limit using sparse deconvolution spatial light interference microscopy,” Biomed. Opt. Express 2(7), 1815–1827 (2011). [CrossRef] [PubMed]
32. C. Li, W. Yin, H. Jiang, and Y. Zhang, “An efficient augmented Lagrangian method with applications to total variation minimization,” Comput. Optim. Appl. 56(3), 507–530 (2013). [CrossRef]
33. “Tval3 code,” http://www.caam.rice.edu/∼optimization/L1/TVAL3/. Accessed: 2016/12/02.
34. P. Thévenaz, U. E. Ruttimann, and M. Unser, “A pyramid approach to subpixel registration based on intensity,” IEEE Trans. Image Process. 7(1), 27–41 (1998). [CrossRef] [PubMed]
35. C. J. Rowlands, S. Varma, W. Perkins, I. Leach, H. Williams, and I. Notingher, “Rapid acquisition of Raman spectral maps through minimal sampling: applications in tissue imaging,” J. Biophotonics 5(3), 220–229 (2012). [CrossRef] [PubMed]
36. N. Pavillon, A. J. Hobro, and N. I. Smith, “Cell optical density and molecular composition revealed by simultaneous multimodal label-free imaging,” Biophys. J. 105(5), 1123–1132 (2013). [CrossRef] [PubMed]
37. N. Pavillon and N. I. Smith, “Implementation of simultaneous quantitative phase with Raman imaging,” EPJ Tech. Instrum. 2(1), 1–11 (2015). [CrossRef] [PubMed]
38. J. De Gelder, K. De Gussem, P. Vandenabeele, and L. Moens, “Reference database of Raman spectra of biological molecules,” J. Raman Spectrosc. 38(9), 1133–1147 (2007). [CrossRef]
39. M. Hedegaard, C. Matthäus, S. Hassing, C. Krafft, M. Diem, and J. Popp, “Spectral unmixing and clustering algorithms for assessment of single cells by Raman microscopic imaging,” Theor. Chem. Acc. 130(4-6), 1249–1260 (2011). [CrossRef]
40. P. J. Rousseeuw, “Silhouettes: a graphical aid to the interpretation and validation of cluster analysis,” J. Comput. Appl. Math. 20, 53–65 (1987). [CrossRef]
