High-speed label-free multimode-fiber-based compressive imaging beyond the diffraction limit

: Glass fibers are miniature optical components that serve as ultra-narrow endoscopy probes. Ideally, one would want to perform imaging through a fiber at the highest achievable resolution and speed. State-of-the-art super-resolution techniques have shattered the diffraction limit, but more than twofold improvement requires fluorescent labeling and a long acquisition time. Moreover, it is challenging to implement super-resolution microscopy in a fiber format. Here we present fiber-based label-free video-rate imaging at more than 2-fold higher resolution than the diffraction limit. Our work paves the way to rapid, sub-wavelength endo-microscopy in unlabeled live specimens.


Introduction
Optical microscopy is an essential tool in various fields of science and technology. For a long time, far-field optical imaging was held back by the Abbe limit: the resolution is confined by λ/(2NA), where λ is the wavelength of light and NA is the numerical aperture of the optical system [1]. The diffraction limit has been bypassed by the help of fluorescent molecules [2]. Many biological discoveries would not have been possible without super-resolution techniques, such as stimulated emission-depletion microscopy (STED) [3], photoactivation localization microscopy (PALM) [4] and stochastic optical reconstruction microscopy (STORM) [5].
Nowadays, optical nanoscopy is widely used worldwide. However, to significantly (more than twice) overcome the diffraction limit, a specific fluorescent labeling is required [6]. Moreover, the acquisition speed of super-resolution microscopy is rather slow. For a 40 × 40 µm 2 image, the acquisition time varies from several seconds per frame for STED and confocal microscopies to several minutes per frame for single-molecule localization methods [7,8]. Fluorescent labelling is not always possible, because it can change the sample properties or even have a destructive influence. The important application area where additional fluorescent markers are not desirable includes wafer defect inspection. Microcracks reduce the fracture strength of Si wafers and lead to wafer breakage during manufacturing, which reduces manufacturing yield. The microcracks can be detected by near-infrared-bright-field imaging [9]. However near-infrared cameras are costly and require temperature stabilization [10]. There is still a high demand for a fast super-resolution technique that does not rely on fluorescent labeling.
The endoscopic probe is often put inside the region of interest to overcome the limitations posed by scattering of light in the tissues [11]. To make the imaging process minimally invasive, the probe should be as compact as possible. Multimode fibers (MMFs) -miniature flexible optical waveguides that can carry thousands of modes -have a massive potential in endoscopic imaging [12]. Recently, compressive imaging through a MMF has been introduced [13]. Multiple speckle patterns generated in a MMF are projected on the sample and a bucket detector records the fluorescence intensity. The sparsity constraint allows the sample to be computationally reconstructed from a small set of speckle patterns combined with the fluorescent signal [14]. The robustness of the approach to fiber bending was demonstrated [15]. The method can be merged with acoustic imaging to successfully produce a hybrid technique that permits imaging with two contrast mechanisms of fluorescence and nonradiative absorption [16].
Here we demonstrate the new method of label-free MMF-based compressive imaging (MMFCI). The advantages of MMFCI over ghost imaging and conventional raster-scan microscopy has been demonstrated by imaging with a spatial resolution 2.5 times better than the diffraction limit at a speed of 5 fps (for 400 × 400-pixel frame) without any fluorescent labeling. Imaging speed up to 250 fps can be easily achieved by implementing fast beam scanning system, such as acousto-optic deflectors. The MMFCI approach is different from traditional methods of MMF imaging, where fields on the fiber output facet are pre-optimized to take the form of diffraction-limited foci and the number of single-pixel measurements is equal to the number of pixels in the final image. In contrast, in MMFCI, the computationally reconstructed image of a sample contains much more pixels than the number of single-pixel measurements.
Fiber-format and demonstrated long-term stability of our approach preface rapid, subwavelength endo-microscopy in unlabelled live specimens.

MMF-based compressive imaging
In MMFCI, the process of projecting different speckle patterns and measuring the transmitted intensity can be formulated as the following linear problem where A is the measurement matrix (M × N 2 ), x is the sample vector (N 2 × 1) and y is the measurement result (M × 1), N 2 is the total number of pixels, M is the number of measurements, n is noise. Fig. 1(a) shows the experimental setup and Fig. 1(b) represents the simplified procedure of MMFCI. Speckle patterns from a MMF are sequentially projected on a sample while a single-pixel detector records total transmitted intensities. One measurement event includes the projection of a single speckle pattern and recording the total signal from a sample. To construct measurement matrix A, each N × N illumination pattern from the MMF, I i (x, y), (x, y are the transverse coordinates) recorded by the camera was transformed into a 1 × N 2 vector. The vectors were stacked one by one to form measurement matrix A. The signal from the sample recorded with a single-pixel detector forms measurement vector y. The sample image is reconstructed by solving a linear problem Eq. (1). Original sample x and noise n are unknown. The problem is ill-posed: in the case of M ≪ N 2 , there are multiple x that satisfy Eq. (1). However a prior information about x can be used to solve Eq. (1) exactly [17]. Compressive sensing uses a sparsity constraint -the knowledge that in some domain x contains a lot of zeros. Sparsity helps to solve the problem by: where ||.|| 1 is the l 1 norm defined as ||x|| 1 = ∑︁ i |x i |. Several well-developed algorithms provide robust l 1 -norm minimisation [18]. Taking into account noise in the measurement system the optimisation problem can be formulated as: where ||.|| 2 denotes the Euclidian norm and τ ≥ 0 is the nonnegative parameter.  [18]. The advantage is that the l 1 norm is replaced by a linear function (functional), making the objective function smooth. Therefore the gradient method of minimization can be applied. For multimode-fiber-based imaging, we used the GPSR algorithm MATLAB code developed for compressed sensing and other inverse problems in signal processing and statistics [19]. Reconstructed vector x (N 2 × 1) was transformed into N × N sample image O(x, y).

Ghost imaging
Ghost imaging uses the same measurement procedure as described by Eq. (1): different speckle patterns are projected on the sample, and the transmitted intensity is recorded. Sample image O(x, y) is obtained by correlating the transmitted intensity with the corresponding illumination pattern intensity distribution [20]: Signal-to-noise ratio (SNR) of ghost imaging scales as the square root of the number of realizations, and in practise M ≫ N 2 patterns are required to achieve SNR ≫ 1 [20]. The spatial resolution of the system is bound by the diffraction limit [21]. For reconstruction by ghost imaging, we used the same experimental data as for MMFCI.

Diffraction-limited raster-scan imaging
Conventional raster-scan imaging was simulated by taking into account the theoretical diffraction limit of the MMF and the known geometrical dimensions of the samples. To obtain O(x, y), ground truth image S(x, y) was filtered in Fourier domain with the 2D incoherent optical transfer function [1]: where S(x, y) is the ground-truth image of the sample, * is the convolution, PSF(x, y) is the incoherent point spread function of the optical system, ν x , ν y are the spatial frequencies, OTF(ν x , ν y ) is the incoherent optical transfer function, which coincides with the theoretical power spectrum of speckles defined by Eq. (7), S F (ν x , ν y ) is the Fourier transform of the ground truth image of S(x, y), F −1 denotes the inverse Fourier transform, |.| denotes the absolute value.

Experimental setup
The experimental setup is presented in Fig. 1 . The cut-off frequency of the magnified speckle pattern was calculated as ν 0 = 2NA/(λ · MG). The total intensity transmitted through the sample was measured by the avalanche photodiode (Thorlabs APD410A2). The incident laser power was adjusted by the half-waveplate and the polarizing beamsplitter cube (CCM1-PBS251/M, Thorlabs) to be 2 -200 µW. The photodiode and the camera were synchronized. In general terms, the resolution can be defined as the minimum separation distance at which two objects can be sufficiently distinguished. As a sample we used round defects with a known size at a known distance from each other (as shown in Fig. 1(c)). The samples were specifically designed for these experiments and made in the cleanroom environment using CMOS-compatible technologies and standard photolithography procedures without any fluorescent labelling. Maskless ultraviolet photolithography, followed by lift-off of sputtered aluminum patterns has been used to fabricate a sample that represents a simple model of the defects on a wafer surface or a microcrack in a silicon solar cell. Samples were accurately characterized by a high-NA optical microscope. The results are summarized in the table in Fig. 1(c). The geometrical parameters of the samples are normalized by the 54× magnified diffraction limit of the MMF. Each sample was placed between the beam splitter and avalanche photodiode in the image plane of the tube lens. Samples were illuminated with the magnified image of the output facet of the MMF. The same image was recorded by the camera.

Speckle-based imaging system
The quality of the reconstruction in compressive sensing largely depends on the properties of matrix A: it is desirable to have a linearly independent set of speckle patterns. Otherwise, some measurements do not add any new information. The absence of correlation between vectors guarantees their linear independence [22]. Here we analyse how the correlation between speckle patterns and their linear independence can be controlled by the scanning strategy (number of points and the distance between points) for a square scanning grid.
Optimal compressive imaging in a diffraction limited system requires the number of measurements to be less or equal to the number of modes, M ≤ M modes [23]. We estimated the number of modes by the formula for an ideal cylindrical optical waveguide [24]: where d is the diameter of the fiber core. Our fiber supports about 1050 modes for a single polarization. The output speckle pattern is the result of interference of propagated MMF modes. An uncorrelated set of speckle patterns requires an optimal scanning grid. We consider a coupling coefficient between the i-th linearly polarized input field E in i (x, y) and the j-th mode Since the NA of the coupling objective is greater than the NA of the fiber, we can approximate the input field as a delta function: , where x 0i and y 0i are the center coordinates of the i-th input beam. If the two input beams are closer than the diffraction limit of the fiber, their coupling coefficients are approximately the same. As a result, they excite two similar sets of fiber modes and the produced output speckle patterns are linearly dependent and correlated. The distance between input scanning points should be not smaller than the diffraction limit of the MMF. However, we can also excite another set of linearly independent speckle patterns by using orthogonal polarization on the input [25].
In our experiments, the distance between adjacent points was set to 1.14 µm that is equal to the diffraction limit of the MMF. The illumination patterns were formed by scanning the focused beam over the input facet on a 31 × 31 points equidistant square grid with a size of 35 × 35 µm ( Fig. 2(a)). The total number of incident speckle patterns was M = 961. The typical speckle pattern produced by the MMF is presented in Fig. 2(c).
In the first set of experiments, we characterized the measured speckle patterns, I i (x, y). We calculated the mean power spectrum density (PS), as PS exp (ν ′ x , ν ′ y ) = <|F (I i (x, y)| 2 >, where ν ′ x = ν x /ν 0 , ν ′ y = ν y /ν 0 are the normalized spatial frequencies. The theoretical PS is [26,27]: The cross sections along ν ′ x for the experimentally measured and theoretical power spectra are presented in Fig. 2(b) by the orange open circles and the solid green line, respectively. We see a good match, and as expected, the experimentally measured PS is limited by a cut-off frequency ν 0 . However, PS exp has a peak at low frequencies and some high frequency components at |ν ′ |>1. To explain the difference, we calculated the PS of the average background as PS bnd The results, presented in Fig. 2(b) by the solid blue line, show that the low-frequency peak in PS exp corresponds to the slowly varying background. High-frequency components can be explained by noise in the background and are several orders of magnitudes lower than the amplitude of the low-frequency part. Complicated effects of light propagation in MMF lead to the existence of a correlation between incident and transmitted fields or optical memory effects [28]. For example, in MMF angular memory effect, the rotation of the input field around the fiber axis results in the rotation of the output field by the same angle. These effects were not taken into account in our simple analysis of scanning strategy optimization and can lead to unexpected correlations of speckle patterns. To ensure that experimentally measured patterns are both linear independent and uncorrelated we calculated the probability distribution of correlation coefficients and singular values. To get rid of high-frequency noise in measured speckle patterns, we filtered the recorded images in the Fourier We characterized the correlation between different speckle filtered patterns by the Pearsson correlation coefficient: where I ′ i (x, y), I ′ j (x, y) are the intensity distributions of two filtered speckle patterns and I ′ i , I ′ j are their mean values. The probability density function of correlation values is presented in Fig. 2(d) by blue dots. The mean value of correlation µ = 0.02 was calculated by fitting experimental The fit is presented in Fig. 2(d) by orange line. Such a low µ means that, in average, speckle patterns are uncorrelated.
We studied the linear independence of the measured speckle patterns. We reshaped each To calculate the rank of matrix A, we performed singular value decomposition. All the singular values are far from zero ranging from 1.5 × 10 3 to 2.5 × 10 5 with a well-pronounced peak around 3000 (Fig. 2(e)). It indicates that there are no values that can be neglected. The rank of matrix A coincides with the number of measured speckle patterns meaning that they are linearly independent. Therefore, we ensured that the diffraction-limited speckle patterns are optimal for MMFCI.

Fiber-based label-free super-resolution imaging
In the second set of measurements, we experimentally investigated the resolution limits for fiber-based label-free MMFCI. The samples were reconstructed using the GPSR algorithm with τ opt (see Appendix). Figure 3(b) represents the experimental results. Columns first till last corresponds to samples 1-4. The high-resolution bright-field images of the samples are presented in Fig. 3(a). We defined the resolution improvement factor as the theoretical diffraction limit divided by the minimum feature size of the sample resolved by MMFCI. The smallest features of samples 1-4 were 1.4, 1.8, 2.5, and 4 times smaller than the diffraction limit, respectively. Our experiments demonstrated that samples 1-3 are nicely resolved by the proposed approach, meaning that the resolution of MMFCI is at least 2.5 times better than the diffraction limit. However, we were not able to experimentally resolve features that are four times smaller than the diffraction limit. Whereas classical compressive sensing theory does not put any restrictions on the resolution, it has been shown that the successful reconstruction cannot always be achieved for a diffraction limited system [23]. The theoretical limit of the spatial resolution improvement factor calculated for our samples was more than an order of magnitude larger than the experimentally achieved values. Unsuccessful reconstruction of sample 4 is not a result of the fundamental limit of MMFCI and can be explained by the presence of noise in matrix A and measurement vector y [23].
We compared the performance of the proposed MMFCI with conventional approaches. Figures 3(c) and 3(d) present the results for raster-scan microscopy and ghost imaging, respectively. Columns first to last correspond to samples 1-4. These methods can resolve only sample 1 because its feature size is comparable to the diffraction limit.
To analyze imaging resolution in more detail, we plotted the cross-sections along the x axis passed through the centers of the dots for the MMFCI (orange line), conventional microscopy (blue line), ghost imaging (white line), and ground truth (green line) in Fig. 3. Ghost imaging and raster-scan approaches can resolve the sample only if the feature size is close to the diffraction limit. Two distinctive peaks in the image of Sample 1 can be easily identified by all the methods (Fig. 3, first row). However, the contrast (the difference between the peak and the gap intensities) is close to one for MMFCI and only about 0.5 for conventional microscopy and ghost imaging techniques. MMFCI also provides an almost perfect contrast for sample 2, whereas it is barely resolvable by conventional imaging and not resolvable by ghost imaging (Fig. 3, second row). While two peaks in sample 3 are not distinguishable neither by ghost imaging nor by conventional microscopy, MMFCI clearly resolves them (Fig. 3, third row). Our results demonstrate the supremacy of MMFCI over classical microscopy techniques in terms of spatial resolution and imaging contrast.
In the next set of measurements, we studied the stability of the MMFCI. We repeated the measurements of sample 3 using the same signal, y, but different matrices A constructed from speckle patterns measured in 2 hours and 2 days. The results are presented in Fig. 4. Even though the resolution slightly degrades in comparison with the measurements performed synchronously, the two dots are well resolved even if 2 days between the actual sample measurements and the speckle pattern pre-calibration have passed. The experimental setup and the proposed approach of MMFCI are stable enough to perform long-term super-resolution measurements with a single pre-calibration episode.

Image compression
We are working in a highly compressive regime by reconstructing a 400 × 400-pixel image from only 961 single-point measurements. Moreover, our computational approach allows for forming final images of almost any size. To characterize the real compression ratio, we take into account the amount of spatial information that we reconstruct. We imaged 1.9 × 1.9 mm sample resolving two objects placed at a distance of 26 µm from each other. To be able to image this sample with a conventional approach of similar quality, we would need to measure at least 150 × 150 pixel. This results in a real compression ratio of 150 2 /961 = 23.
In the final set of experiments, we investigated the minimal number of measurements required to resolve the sub-diffraction features with the proposed approach. The correlation between the reconstruction results and the ground truth as a function of M for samples 1-3 is presented in Fig. 5(e). For M ≥ 48, correlation R O,S >0.3, meaning that the reconstructed image and the ground truth are correlated even for such a very low number of measurements. The correlation increases with the number of measurements followed by the saturation, which is marked by a dotted blue line in the Fig. 5(e) and corresponds to M sat = 192. At this moment, two dots can be resolved, as shown in Fig. 5(b). To ensure the good reconstruction quality, we chose a twice higher number of measurements despite the fact that the correlation does not show rapid growth with the number of measurements for M ≥ M sat . We show that M = 432 provides to R O,S >0.8 for samples 1-3. Thus, about 432 measurements are enough to get a high-quality image with MMFCI approach. As a result, the compression ratio can be increased up to 150 2 /432 ≈ 52.
High compression allows to improve the imaging speed significantly. Because pre-calibration can be done several hours in advance (as we showed in the previous section), the measurement speed is only limited by a beam scanner. Imaging speed of 200 ms per frame (5 fps) is available with the standard galvo mirror system.

Conclusion
Here we report on the successful reconstruction of sub-diffraction features (super-resolution) by fast fiber-based label-free compressive imaging. The prior knowledge that the sample is sparse allows achieving superior image quality and speed compared to the conventional techniques. We believe that sparsity-based imaging can be applied to any sample because all natural objects have a sparse representation in an appropriate basis, which is supported by many popular image compression techniques, such as JPEG.
We have experimentally demonstrated the imaging speed of 5 fps and the spatial resolution 2.5 times better than the diffraction limit without using any fluorescent labeling. In our experiments, the imaging speed was limited by the time that galvo mirrors switch between adjacent points of the scanning grid. Acousto-or electro-optical deflectors, which do not contain any mechanical moving parts, can improve scanning speed to approximately 10 µs per point [29]. As a result, high-speed super-resolution imaging with 4 ms per frame or 250 fps for 1.6-megapixel images can be achieved by using non-mechanical beam scanning system. It is three orders of magnitude faster than state-of-the-art imaging through the same multimode fiber with a diffraction limited resolution [8].
Currently, the results are largely influenced by choice of the reconstruction algorithm. As we show, the best performance can be achieved with the proper selection of regularisation parameter. The future direction would be the development of algorithm optimization methods that do not rely on the prior information about the sample.
To summarize, the proposed MMFCI approach merges high-speed super-resolution optical label-free imaging with a potentially compact optical probe. We have shown that MMFCI is stable enough to be used for long-term (more than a day-long) ultra-high-speed super-resolution imaging. We experimentally demonstrated the resolution improvement better than could be achieved by state-of-the-art label-free techniques, such as structured illumination microscopy. Super-resolution fluorescence methods such as STED, PALM and STORM can theoretically achieve infinitely small resolution. However this comes at a cost of long acquisition time of about several minutes per frame [7]. Therefore, the demonstrated imaging speed of MMFCI is significantly higher than the imaging speed of fluorescence super-resolution microscopy. The fiber nature of the proposed approach paves the way toward high-speed and high-quality endoscopy.

Appendix: Algorithm optimization for a MMF-based super-resolution imaging
We experimentally investigated the optimal parameters of the GPSR reconstruction algorithm for MMFCI. We illuminated the sample with M = 961 speckle patterns produced by the MMF and recorded the total transmitted intensity for each illumination pattern. The intensities formed measurement vector y and the speckle patterns filtered in the Fourier domain formed measurement matrix A.
In Eq. (3), regularisation parameter τ defines the balance between the l 1 norm and noise. The reconstruction results are greatly influenced by τ: too large values of τ leads to overestimation of sample sparsity. For example, for τ>τ max = |A T y| ∞ , x = 0 and the reconstruction fails [30].
To study the influence of the regularisation parameter on the image quality, we performed the reconstruction with different values of τ. The reconstructed images for Sample 1 in comparison with ground truth are presented in Fig. 6(a). For τ<τ max , the algorithm always reproduced the number and the positions of dots. To quantify the image quality, we calculated the Pearson correlation coefficient R O,S between the reconstructed image and the ground truth. The correlation coefficient as a function of normalized regularization parameter τ/τ max for all the samples is presented in Fig. 6(c). We defined τ opt as the regularization parameter for which R O,S is maximal. We see that for τ opt <τ<τ max , the feature size is smaller and for τ<τ opt is bigger than the ground truth ( Fig. 6(a)).
To explain this behavior, we performed numerical experiments with the measured set of speckle patterns and corresponding A. We simulated the photodiode signal as y sim = Ax sim +n, where x sim is the 1D vector (N 2 × 1) reshaped from the ground truth image and n = max(Ax sim ) · r · ϕ/100%, where max(Ax sim ) is the maximum value of noise-free photodiode signal, ϕ is a noise ratio in percent, and r is a vector (M × 1), which entities are randomly and uniformly distributed numbers in the range (0,1).
We simulated the imaging procedure of samples 1-4 with different noise ratio ϕ and different values of τ. The examples of Sample 1 reconstructions for 5 % noise level are presented in Fig. 6(b). The correlation coefficients as functions of τ/τ max are presented in Fig. 6(e). Different columns correspond to different samples, and violet, yellow, red and blue colors correspond to ϕ = 50%, 10%, 5% and 0%, respectively. In the noise-free case, the correlation function for all the samples has a plateau -a range of τ values where the correlation is equally high (Fig. 6(e), blue areas). However in the presence of noise, the R O,S curves have peaks that define τ opt . For τ<τ opt , the algorithm overestimates the feature size, for τ opt <τ<τ max , the size is underestimated (Fig. 6(b)) similar to the experiment (Fig. 6(a)). Higher noise level corresponds to the larger value of τ opt . It is in agreement with Eq. (3), where the larger value of τ makes the sparsity term more significant compared to the difference between measurement result y and Ax. As (c) Correlation between the measured image and the ground truth as a function of τ/τ max for sample 1 (violet), sample 2 (yellow), sample 3 (orange) and sample 4 (blue). (d) Optimal regularisation parameter τ opt /τ max as a function of noise ratio in simulations for sample 1 (solid violet line), sample 2 (dashed yellow line), sample 3 (solid orange line), sample 4 (solid blue line). (e) Correlation between reconstruction result and ground truth as a function of regularization parameter τ/τ max for the noise-free case (blue), 5 % noise ratio (red), 10 % noise ratio (yellow), 50 % noise ratio (violet). From top to bottom: sample 1, sample 2, sample 3, sample 4 the noise in the measurement increases, the value of the term ||y − Ax|| 2 in problem Eq. (3) increases. As a result, the sparsity term ||x|| 1 is overestimated and the size of the sample features is underestimated. To reconstruct sample vector x correctly the significance of the sparsity term ||x|| 1 should be increased by increasing τ. The higher noise level corresponds to the smaller value of R O,S or the degradation of reconstruction (Fig. 6(e)).
The numerical result qualitatively reproduces experimental behavior. To achieve the highest quality of the MMFCI image, GPSR algorithm requires the optimized value of regularisation parameter τ opt . In case of τ ≠ τ opt , the feature size is reconstructed incorrectly. For τ<τ opt the resolution of the image degrades, for τ opt <τ<τ max the feature size is underestimated and for τ ≥ τ max the reconstruction fails. The value of optimal regularisation parameter τ opt is specific to both the sample and the noise level. To overcome this problem either an empiric adaptive choice or an automated selection of regularisation parameter can be used [31,32]. In further experiments, we use the experimentally determined optimal values of τ.