Virtual single-pixel imaging-based deconvolution method for spatial resolution improvement in wide-field fluorescence microscopy

: Deconvolution technique has been widely used in ﬂuorescence microscopy to restore ﬁne structures of biological samples. However, conventional deconvolution methods usually achieve little contrast enhancement in dense structures that have the intervals close to the Rayleigh criterion. Herein, we developed a novel deconvolution method, termed virtual single-pixel imaging ( v -SPI). Diﬀering from existing deconvolution methods, v -SPI aims to retrieve the less blurred image directly, not the sample distribution which cannot be actually obtained. And the result can be retrieved simply by solving a linear matrix in spatial domain. In addition, the proposed method has no requirement of calibrating parameters of microscope system. Simulation and experimental results demonstrated that the proposed v -SPI method can enhance the contrast of dense structures signiﬁcantly and acquire a 24% increase in resolution.


Introduction
The light diffraction limits the spatial resolution of optical microscopy and thus renders distinguishing subcellular structures at length-scales smaller than the Rayleigh criterion [1][2][3][4]. Various super-resolution (SR) microscopies have been developed to improve spatial resolution. These SR methods narrow the shape of the point spread function (PSF) mainly through changing the microscope design. Although these SR approaches, such as confocal [5], stimulated emission depletion (STED) [6] and structure illumination microscopy (SIM) [7], can effectively break the diffraction limitation, they however should build technical microscope setups, which require careful calibration and stable maintenance. Otherwise, artifacts will be generated in the reconstructed images.
Theoretically, the blurred fluorescence image can be taken as the convolution result of samples and microscope PSF [8][9][10]. It thus is possible to retrieve images numerically with fine details and enhanced contrast. Comparatively, numeral processing method (so-called deconvolution method) is more flexible, and can be applied to the simplest microscope system or taken as the post-processing tool for other SR methods [11][12][13]. Because it is useful in many contexts, various deconvolution algorithms have been proposed [14][15][16][17][18][19], such as the Wiener filtering, the Gold algorithm, the Tihkonov regularization method and the Richardson-Lucy (RL) method, etc. Most of the existing deconvolution methods estimate the sample distribution over multiple iterations of creating a blurred image with system PSF and comparing it to the original. It not only requires precise calibration of system PSF, but also has the difficulty in deciding when to end the iteration procedure as the 'real' sample distribution actually cannot be obtained. Moreover, existing deconvolution methods mainly increase the image contrast for sparse structures (the minimum distance between two objects is higher than the Rayleigh criterion) but shows little contrast enhancement to the area of denser structures [20]. With more rigorous resolution definition (relating to the minimum distance that two objects can be distinguished from one another), this structure-dependence performance of existing deconvolution methods just makes the structures shrink but with little spatial resolution improvement.
To this end, a novel deconvolution method termed virtual single-pixel imaging (v-SPI) is presented in this work. After combining an introduced equivalence relationship of integral value and SPI algorithm, it retrieves results by simply solving a linear matrix in spatial domain. Besides, we also altered our estimation objective from the sample distribution to the less blurred sample distribution (blurred with a smaller PSF) which is possible to be obtained in practice. It is verified by the simulation and experimental results that, with all these efforts, the proposed v-SPI method can improve the spatial resolution effectively by enhancing the contrast of dense structures and without calibrating system parameters.

Principle
In v-SPI method, the blurred fluorescence image is post-modulated with digital patterns. We first deduced an equivalence relationship, i.e., the integral value of the post-modulated image is equal to that of the dot production result of sample distribution and convoluted pattern (with system PSF). The reverse problem of deconvolution was then converted into solving a linear matrix with SPI algorithm. In the linear matrix, convoluted patterns and the corresponding integral values are known, the sample distribution is desired.

Equivalence relationship of the integral value
The schematic diagram of equivalence relationship of the integral values is shown in Fig. 1. From the diffraction theory, the blurred fluorescence image I can be defined as the convolution result of sample density S and microscope PSF h, as follows: where, x and y indicate the coordinates in both directions on the image plane, and x' and y' the coordinates in both directions on the specimen plane. For simplicity, the theoretical model for imaging formation is described here in two-dimensional space. We then post-modulate it with digital patterns P i : where M i is the modulated image with respect to the i th digital pattern P i , and m denotes the number of digital patterns. The integral value of the modulated image M i may be calculated as: By exchanging the order of the two integrals [21][22][23] in Eq. (3), we rewrite this expression as: We may now define the convoluted patterns P i as By substituting Eq. (5) into Eq. (4), we acquire the following expression From Eqs. (3)-(6), we can draw the equivalence relationship of integral value, i.e., the integral value of modulated image M i is equal to that of the dot production result of sample distribution S and convoluted pattern P i , which can be simply expressed as where * represents convolution operation and • denotes dot product operation.

Principle of v-SPI
Based on the equivalence relationship, the integral value of the modulated image M i is calculated, it can be taken as that of the dot product of sample distribution S and convoluted pattern P i . Since the convoluted pattern P i is determined with digital pattern P i and the PSF h. We thus can retrieve the desired sample density S by solving the following equation: where A ∈ R m×n indicates the modulation matrix (m convoluted patterns P i and each pattern consists of n pixels), X ∈ R n×1 denotes the reshaped result of the desired sample density S (aligned as a vector), and b ∈ R m×1 is a vector that consists of m integral values. The schematic diagram of v-SPI method is shown in Fig. 2. By referring to the SPI technique [24,25], we find that the Eq. (8) is similar to the measurement formation in SPI technique. From Eqs. (3) and (5), we can determine the integral values and convoluted patterns with the blurred fluorescence image I, digital patterns P i and PSF h. Thus, the desired sample distribution S can be reconstructed using the SPI algorithm. In addition, a set of specific digital patterns is applied here, in which each pixel will be sequentially set to a value of 1, and the rest of the pixels are set to zero. That means the number of digital patterns is equal to the pixel number of the blurred fluorescence image. With the specific digital patterns, this reconstruction involves solving a huge matrix equation (16384*16384 matrix for reconstructing 128*128 image with 16384 patterns). By referring to the comparison of different SPI algorithms in [24] and taking the computer memory limit into account, a linear iterative algorithm (Alternating Projection, AP) is selected in this paper to retrieve the desired sample distribution S.

Practical usage of v-SPI
Although we provided an ideal numerical retrieving strategy with Eq. (8), however the 'real' sample distribution S cannot be actually obtained, like in existing deconvolution methods. Usually, retrieved results are just less blurred images than the originals. Based on this fact, we changed the estimating goal of v-SPI from retrieving the sample distribution to the less blurred distribution. As the 2D system PSF can be accurately described with Gaussian model [26], we then can regard the system PSF h (2D Gaussian approximations) as the convolution result of two smaller Gaussian PSFs (h 1 and h 2 ), which have smaller half full-width at half-maximum (FWHM) than that of system PSF h. By substituting it to Eq. (7), we can acquire: where S = S * h 1 denotes the less blurred distribution. It is demonstrated from Eq. (9) that the less blurred deconvolution result S' can be retrieved by v-SPI with the smaller PSF h 2 . Compared with the unreachable objective (retrieving sample distribution) of existing deconvolution methods, our changed estimation objective in v-SPI is possible to be achieved and thus can yield more reliable results. As any smaller PSF h 2 can be applied in v-SPI to retrieve the less blurred distribution S', the proposed usage of v-SPI actually eliminates the requirement of accurate PSF calibration, which is necessary for existing deconvolution methods.

Simulation results
The performance of v-SPI was first evaluated using the numerically simulated image of a thin spoke-like sample placed at the objective focal plane (with numerical aperture NA=1.45) [ Fig. 3(a)]. In the numerical simulation, the sample was uniformly illuminated and the wavelength of fluorescent light is 510nm (λ em =510nm). The image pixel size was set to 30 nm. At first, noise was not introduced in the blurred image to evaluate the performance under the ideal condition. With these simulation settings, a 2D Gaussian PSF model h is utilized to blur the thin spoke-like sample and the blurred result is shown in Fig. 3(b). Although any smaller PSF h 2 can be used in v-SPI method to retrieve results, they however may lead different results of the reconstruction and we need to seek the most suitable PSF in practice. In this paper, the 2D PSF is described with Gaussian model. And the comprehensive standard deviation (CSD,σ = σ/ d,σ is the standard deviation of Gaussian model and d is the pixel size) is defined to determine the size of PSF. In this simulation, the CSD of nominal system PSF is 2.97 (σ = 2.97). The v-SPI results [Figs. 3(c)-3(i)] were reconstructed using Gaussian PSF model h 2 in different sizes (σ 2 =0.9, 1.2. . . , 3.0), respectively. By comparing the v-SPI results [Figs. 3(c)-3(i)], it can be found that the dense strips in the central zone gradually become visible with increased σ 2 . However, the oversized PSFs (σ 2 >1.5) may have little contribution to the resolution improvement and cause artifacts in the reconstructed results. Therefore, the smaller 2D PSF h 2 with CSD of 1.5 (σ 2 = 1.5) was used in v-SPI to retrieve results in this simulation. It should be noted that the suitable CSD of the PSF h 2 (σ 2 ) may vary with respect to the changed CSD of system PSF (σ ). The selected CSD (σ 2 = 1.5) only can be used for the case simulated here.
Then, three different deconvolution methods, including Regularized inverse filtering (RIF) method, Richardson-Lucy (RL) deconvolution method and RL with total-variation regularization (RL-TV) method, are adopted for performance comparison. As seen in Figs. 4(c)-4(f), it exhibits the retrieved images of v-SPI, RIF, RL and RL-TV, respectively. The result of v-SPI was obtained by using MATLAB (R2015b) on a computer equipped with 8 CPU (Intel Xeon E5-1620 v4 @3.5 GHz). Results of other three deconvolution methods are computed with a free software package (DeconvolutionLab2) [8]. Parameters of used deconvolution methods are set as follows: regularization parameter λ RIF =0.01, iterative number N RL =100, regularization parameter λ RL−TV = 1E-12 and iterative number N RL−TV = 20. The simulated 2D system PSF is used by three comparative methods (RIF, RL and RL-TV). To estimate the ability of contrast enhancement, we calculated the contrast C(p) from the blurred image and results of four numerical processing methods (v-SPI, RIF, RL and RL-TV) as a function of the period (in wavelength, λ em =500nm): where f p represents the intensities corresponding to the period p, f max and f min are maximum and minimum intensity of the whole image. The calculated contrast curves of the blurred image and results of four deconvolution methods are shown in Fig. 4(g). The black dot dash line denotes the contrast about 0.1, i.e., the contrast is still 'eye-visible' [20]. This figure illustrates the fact that, there still exist areas in the blurred fluorescence image with small contrast but 'eye-invisible' when the periods are smaller than the Rayleigh criterion (0.61λ em / NA ≈ 0.42λ em ).
As shown in Fig. 4(g), the RIF method and the RL-TV method mainly increase the contrast for sparse structures (periods bigger than the Rayleigh criterion), while the RL method can work for the more dense structures (0.37λ em , slightly smaller than the Rayleigh criterion). In comparison, the uniform and significant contrast enhancement can be seen in the result of v-SPI even when the periods smaller than Rayleigh criterion (till to 0.3λ em ). It means that, by greatly enhancing the contrast of the areas which originally have weak contrast but 'eye-invisible', v-SPI method acquired the resolution improvement of 0.12λ em . Comparing the radius of yellow circles in Figs. 4(b)-4(f), we can intuitively see the resolution improvement of v-SPI than the blurred image and results of other three deconvolution methods.
To assess the noise sensitivity of the proposed method, the simulation fluorescence image [ Fig. 5(a)] was contaminated by Poisson noise and then processed by using v-SPI, RIF, RL and RL-TV respectively, as shown in Figs. 5(b)-5(e). The signal to noise ratio (SNR) was calculated as follows to quantitatively assess the noise sensitivity [27], where,Ī signal andĪ noise are the mean values of signal and noise components, respectively. As shown in Fig. 5, tiny Poisson noise in the blurred fluorescence image nearly have no influence on the results of RIF, RL and RL-TV, but cause obvious noise in the v-SPI result (SNR was reduced from 36 dB to 19 dB). It is verified that, the v-SPI method is more sensitive to noise than other three deconvolution methods.

Experimental results
After numerical study, the proposed v-SPI method was then tested with experimental images obtained using a wide-field microscope (Nikon Ti-2A) in epi-illumination mode. A commercial sample (Fluocells prepared slide#1, Invitrogen) with stained BPAE cells was used, in which F-actin was stained using green-fluorescent Alexa Fluor 488 phalloidin. The F-actin in the fixed BPAE cells was uniformly illuminated through a high-numerical aperture objective (NA = 1.45, ×100, Plan Apo λ, Nikon) at an excitation wavelength of 488 nm. Emitted fluorescent light was filtered using a band-pass filter (ZET488/640m, Chroma, λ em = 510 nm). The fluorescence light was then imaged by a Scientific CMOS camera (ORCA-Flash 4.0, HAMAMATSU) with pixel size of 65 nm (after magnification). In this experiment, the v-SPI results were reconstructed with the Gaussian PSF models h 2 which has the CSD of 1.1 (σ 2 = 1.1). In comparison, the fine structure of F-actin which is difficult to distinguish in the captured wide-field image [ Fig. 6(a)], was retrieved in the v-SPI image [ Fig. 6(b)]. It implies the obvious resolution enhancement of v-SPI method. As shown in the intensity profiles [ Fig. 6(g)], two close F-actin which have a distance of 167 nm between each other, can be distinguished in the v-SPI image. And the distance is 47 nm shorter than the Rayleigh criterion (0.42 λ em , 215nm). In order to verify the factuality of v-SPI result, the sample was also imaged with an established SR technique [11] (multifocal structured illumination microscopy, MSIM), as shown in Fig. 6(f). From the retrieved images and the intensity profiles, we can see similar structures of F-actin both in v-SPI result and MISM image. By comparing the retrieved images of four deconvolution methods (v-SPI, RIF, RL and RL-TV) and the intensity profiles, we can clearly see that much larger enhancement in image contrast can be obtained by v-SPI method.
To verify the applicability of v-SPI in imaging live cells, we also performed dynamic imaging of mitochondria in living COS-7 cells, as shown in Fig. 7 and Visualization 1. The sample was stained with commercial MitoTracker (MitoTraker Green FM probe, Thermo Fisher Scientific). The mitochondrial dynamics was recorded using a wide-field microscope. In this experiment, the sample was uniformly illuminated through a TIRF objective lens (NA = 1.49, ×60, ApoTIRF, Nikon) at an excitation wavelength of 488 nm. The wide-field images were captured with 1.5-fold interpolation (using the microscope). During the image acquisition process, the exposure time was set to 500 ms, to improve the SNR of the captured images. From the captured wide-field image of mitochondria [Figs. 7(a) and 7(c)], we find the difficulty in observing the clear structure and dynamic movements of mitochondrial cristae in living cells. However, after processed with v-SPI method, the complicated structures of mitochondria in the living cells become visible [see Figs. 6(b)-6(c), Visualization 1] due to contrast enhancement and improved spatial resolution. By using v-SPI method, we could even identify the fusion process of cristae in a single mitochondrion [ Fig. 6(b)].

Discussions
Although the excellent performance of proposed v-SPI method has been demonstrated with simulation and experimental results, many considerations also need to be taken into account in practice. At first, the v-SPI method belongs to a 2D deconvolution method, thin-layer sample imaging is required to avoid the influence of defocus signals. However, we still can image thick sample by combining v-SPI method with optical-sectioning techniques, such as confocal, SIM and SD method [28].
When a suitable PSF h 2 is used in v-SPI method, we can obtain the reconstructed results with better contrast enhancement and less artifacts. Actually, as the wide-field images may be captured with various microscope systems, the suitable CSD of the PSF h 2 (σ 2 ) should be determined according to the CSD of system PSF (σ ). Here, to facilitate the application of v-SPI method in practice, we provide some empirical values (σ 2 ) which can be selected based on the CSD of nominal system PSF (σ ) and shown as following: In view of the Nyquist criterion, the pixel size d should be less than the half FWHM of system PSF (d ≤ FWHM/2). And considering the relationship between FWHM and standard deviation of system PSF (FWHM=2.355σ), small CSDs (σ <0.85) thus do not listed in Eq. (12). For a microscope system with too large CSD (like in the simulation,σ = 2.97), the selected CSD of PSF h 2 should be less than 1.5 (σ 2 ≤ 1.5) to avoid artifacts.
Finally, the boundary effect should be noted in application of the proposed v-SPI method. The fluorophores that are outside the fluorophore density x and close to the borders also contribute to the intensity of captured wide-field image I. However, the influence of these fluorophores is neglected in the reconstruction procedure, which will cause artifacts in the edge areas of the reconstructed image. In order to address the boundary effect, we can simply cut the false edge areas out from the reconstructed image. Usually, the width of removed edge areas can be set as the FWHM of system PSF.

Conclusions
In this paper, a novel deconvolution method is proposed to improve the spatial resolution, which does not require calibration of system parameters. From simulation and experimental results, it is revealed that an increase in resolution of 24% (0.12/0.5) is achieved by v-SPI method. Furthermore, it may be applied as a post-processing approach for other diffraction-limited or super-resolution techniques (such as confocal, two-photon, and STED methods) to further improve spatial resolution.