Super-resolution microscopy for biological specimens: lensless phase retrieval in noisy conditions

The paper is devoted to a computational super-resolution microscopy. A complexvalued wavefront of a transparent biological cellular specimen is restored from multiple intensity diffraction patterns registered with noise. For this problem, the recently developed lensless super-resolution phase retrieval algorithm [Optica, 4(7), 786 (2017)] is modified and tuned. This algorithm is based on a random phase coding of the wavefront and on a sparse complex-domain approximation of the specimen. It is demonstrated in experiments, that the reliable phase and amplitude imaging of the specimen is achieved for the low signal-to-noise ratio provided a low dynamic range of observations. The filterings in the observation domain and specimen variables are specific features of the applied algorithm. If these filterings are omitted the algorithm becomes a super-resolution version of the standard iterative phase retrieval algorithms. In comparison with this simplified algorithm with no filterings, our algorithm shows a valuable improvement in imaging with much smaller number of observations and shorter exposure time. In this way, presented algorithm demonstrates ability to work in a low radiation photon-limited mode. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Digital holography (DH) is a universal technique applied to a wide range of tasks for objects investigation and control. To name some of them: shape measurement [1,2], vibration control [3], label-free analysis [4], quantitative biological samples phase imaging (QPI) [5], etc. DH techniques have a serious advantage for biological samples because many of them are transparent and, for visualization in conventional microscopy, coloring [6] or fluorescent agents [7] are often used which can be harmful especially for live objects [8]. DH allows a harmless introvision by reconstruction 3D complex-valued wavefront of transparent objects without additional preparations of samples [9].
One of the DH's advantages is an ability to reconstruct the object phase, which brings information about a thickness and refractive index of the object [10]. Estimation of these parameters opens great possibility for science in investigation of live creatures. For example, knowing the sample refractive index and basing on the phase information it is possible to calculate dry mass of the samples [11]. Additionally, based on the phase information, DH allows to plot 3D movement maps of the objects in an investigated media that is a great benefit for ecological aims [12] and behavioral descriptions of live objects [13].
Phase retrieval is another technique to reconstruct a complex-valued wavefront of objects based only on intensity measurements without using reference beams. The iterative Gerchberg-Saxton (GS) algorithms based on alternating projections of wavefronts between the object and sensor planes are prominent instruments for phase retrieval [14,15]. For higher accuracy and fast convergence of iterations the phase retrieval techniques need multiple measurements and their significant diversity. The sufficient diversity can be achieved in many different ways, in particular: by varying a propagation distance between object and registration plane [16]; by using different wavelengths [17]; by defocussing [18]; by introduction of wavefront modulation masks (amplitude [19] or phase [20]); etc. Phase retrieval as compared to DH has an advantage in the areas where DH cannot be applied due to complexity of imaging systems and, therefore, compelled absence of the reference wavefront ( x-rays [21] or low-energy electron microscopy [22]). Another great advantage of phase retrieval is that it could be realized in a lensless configuration contrary to state-of-the-art DH and QPI methods [5,11].
In work with live tissues, it is always need to keep in mind that used radiation could harm an object or change its behavior. Especially it is crucial for X-rays with high-power radiation. Therefore, to minimize radiation influence on the object the lowest possible radiation power should be used for measurements. In turn, a low radiation power always drops down an intensity of the registered holograms, which results in a low signal-to-noise ratio (SNR) and erroneous and noisy reconstructions of objects. In these cases, an intelligent filtering should be applied to eliminate noises without losses in a true signal. Nowadays, the advanced noise suppression techniques are successfully embedded in phase reconstruction algorithms both in DH [23] and phase retrieval [24][25][26].
One of the state-of-the-art noise filtering algorithms is Block Matching 3 Dimensional (BM3D) filtering [27]. It is based on sparse modeling of objects and successfully applied in various optical setups: off-axis DH [28], phase shifting DH [29] and phase retrieval [30]. The recent Super-Resolution Sparse Phase Amplitude Retrieval (SR-SPAR) algorithm is developed for lensless setup [31]. BM3D filtering of phase and amplitude are important components of this algorithm promising for microscopic imaging and multiwavelength absolute phase retrieval [32,33].
The contribution of this paper concerns application and demonstration of high-efficiency of the SR-SPAR algorithm for super-resolution imaging of biological cellular specimens in a noisy conditions. The algorithm is especially modified and tuned to deal with low radiation of laser beams and, respectively, extremely noisy observations. The paper is organized as follows. The phase retrieval problem in noisy conditions is posed in section 2. The SR-SPAR algorithm is presented in section 3. This section is focused on the structure and the details of the algorithm. The experimental setup is a topic of section 4. Section 5 demonstrates application of SR-SPAR for imaging of epithelial cells and erythrocytes. These results are given for different exposure time, i.e. for different levels of noise in observed diffractive patterns. The conclusions are given in section 6.

Problem formulation
In the computational lensless phase microscopy, it is commonly assumed that the whole information about object scattered coherent wave passed through the object is needed to reconstruct the object image. It implies the reconstruction not only amplitude but the phase as well. Contrary, the phase retrieval is a problem of complex-valued object's wavefront (u = b · exp( jϕ), where b is an amplitude, and ϕ is a phase) reconstruction provided that only intensities of light radiation are registered. In our wavefront manipulation modeling, we assume that the object is flat and the planes of the object and sensor are parallel, the laser beam is perpendicular to the object plane.
Following to [31], we formalize the measured intensities in the form: where: u(x, y, 0)∈ C N ×N is an N × N complex-valued object's wavefront; x, y are spatial coordinates, 0 means that the free propagation distance d equals to zero, that corresponds to the object plane; P s : C N ×N −→ C M×M is a complex-valued forward propagation operator from object to sensor planes, y s (x, y, d)∈ R M×M It is useful to note, that the object's wavefront u(x, y, 0) in this notation is the wavefront just after the object plane. It coincides with the object transfer function provided that the input light beam is coherent monochromatic parallel to the optical axis of the system and with the magnitude equal to 1. In our setup as in [31], the object wavefront is modulated by phase masks, such that u(x, y, 0) is changed for M s (x, y) · u(x, y, 0). Then the forward propagation from the object to the sensor plane can be represented as Here u s (x, y, d) is a complex-valued wavefront at the sensor plane at distance d from the object plane, M s ∈ C N ×N are phase-only modulation masks, M s (x, y) = exp( jφ s (x, y)). Then, Eq. (1) takes the form: The phases φ s (x, y) in M s (x, y) are generated as random. It results in the observations known as "coded diffraction patterns" (e.g. [34]). The phase modulation is able to change vastly the diffraction pattern of P{u(x, y, 0)} enabling redistribution of observed intensities from low to higher frequencies.
For noisy observations Eq. (3) is changed for where G stands for a generator of random variables (observations). The wavefront intensity registration process in optics amounts to count the photons hitting the sensor's elements and is well modeled by independent Poisson random variables (e.g. [35]). The mean and the variance of Poisson random variable z s are equal and given by y s , i.e., E {z s } = var{z s } = y s , here E {} is a mathematical expectation, var{} is a variance. Defining the observation signal-to-noise ratio (SNR) as the ratio between the square of the mean and the variance of z s , we have SN R = E 2 {z s }/var{z s } = y s . In real experiments the intensities y s are scaled by the exposure time t ex p , which essentially controls noisiness of observations. Smaller values of t ex p lead to smaller SN R, i.e. to noisier observations z s (x).
For the wavefront propagation P from the object to the sensor plane, we use the Rayleigh-Sommerfeld model with the transfer function defined in the Fourier domain through the angular spectrum (AS) [31].

Super-resolution algorithm for phase retrieval
The SR-SPAR algorithm derived from the variational formulation of the phase retrieval problem for optimal reconstruction of the object from noisy data is presented in [31]. In this section, we go through the structure of this algorithm in descriptive manner following the successive operations step-by-step. For details of the algorithm, we refer to [31].
The computational wavefront reconstruction is going from the continuous domain wavefront propagation to the corresponding discrete models based on pixelation of the object's wavefront distribution u(x, y, 0), thus, the discrete modeling of the system at hand links the discrete values of sensor output (observations) with the discrete values of the object, where the integral Fourier Transform (FT) is replaced by Fast Fourier Transform (FFT). This discretization is necessary both for digital data processing as well as for modeling of wavefront propagation and image formation. Contrary to the sensor pixels ∆ S size defined by the corresponding optical-electronic devices, the object pixels ∆ o are computational ∆ o = ∆ c , which maybe be taken arbitrary small.
In the pixel-wise phase retrieval ∆ c = ∆ S , therefore the resolution of the object u(x, y, 0) reconstruction from the observations z s (x, y, d) is dictated by the pixel size of the sensor ∆ S . In Input: {z s }, {M s }, s = 1, ..., L Initialization: super-resolution ∆ c < ∆ S . It is conventional to apply an integer up-sampling factor r S ≥ 1 as ∆ S = r S · ∆ c . Thus, a larger r S means a smaller computational pixel ∆ c and a larger order of the super-resolution. All wavefronts in our algorithm are calculated with the discretization given by ∆ c .

Object wavefront update:
To avoid aliasing effects in the AS propagation modeling, the zero-padding of the wavefronts is used of the size defined from the inequality [36]: which binds the number of pixels in one dimension of the zero-padded object N and propagation distance d, for given values of the computation pixel size ∆ c and the wavelength λ. Figure 1 shows the block-scheme of the SR-SPAR algorithm [31]. The inputsz s in this algorithm are the observations z s upsampled by factor r s , and M s are known phase masks introduced into the setup. We use the zero-order upsampling givingz s with piece-wise invariant values for computational pixels corresponding to each of the larger size pixels of the sensor. At the initialization step the initial estimate of the object wavefront x 0 is created. It is computed as the backward propagation of the wavefronts which amplitudes equal to the square root of the upsampled observationsz s and phases of these wavefronts are zeros. Then, the phases of the masks M s are subtracted from the propagated wavefronts and x 0 is obtained as the average of the estimates obtained for each of the L observations, the '•' stands for the entry-wise product of two matrices. Here M * s means a complex conjugate M s . At Step 1 (forward propagation), the object wavefront estimate x t is multiplied by the phase mask M s and propagates by the operator P to the sensor plane, the result is denoted as v t s . These wavefronts are calculated for the computational diffraction area N × N denoted as S D . At Step 2 (Poissonian noise suppression), the wavefronts v t s are updated to the variables u t s by filtering the amplitudes of v t s . The filtered amplitudes b t s are obtained according to the algorithm optimal for Poissonian observationsz s [30]: The filtering at Step 4 (sparse phase and amplitude filtering) is one of the key elements of the SR-SPAR algorithm. It follows from the sparse representations of the object phase and amplitude. In this modification of the algorithm the sparsification of the object is implemented for the wrapped phase. As a result the BM3D phase filtering is produced for the wrapped phase, what makes a difference with the SR-SPAR algorithm in [31], where the sparsification and BM3D filtering are produced for the absolute phase. The step of phase unwrapping is excluded in presented version of the SR-SPAR algorithm [31]. At Step 5 (object wavefront update) the iteration of estimate is updated. SR-SPAR is the iterative algorithm, that needed a number of iterations to produce a reliable result. The number of iterations T depends on the object structure and on the level of the noise in the observations, but usually, T = 30 is enough. After T iterations at the output of the algorithm, final reconstructions for phase and amplitude are obtained.
The success of the sparsity modeling and the corresponding filtering depend on how rich and redundant are dictionaries (sets of basis functions) used for complex-valued image analysis and synthesis. In the SR-SPAR algorithm, the BM3D sparsity is achieved due to a quite artificial procedure, where similar blocks of the images are grouped together and processed jointly [27]. The nonlocal block-wise sparsity imposed on phase and amplitude results in the separate BM3D filtering of phase and amplitude as it is shown in Step 4: Hereφ andb are sparse approximations of ϕ and b; phase and ampl as indices of BM3D are used in order to emphasize that the parameters of BM3D can be different for the phase and amplitude; T h ϕ and T h b are thresholding parameters of the algorithms [27]. Let the filtering at Steps 2 and 4 are dropped. Then, the SR-SPAR algorithm can be treated as a super-resolution version of the GS algorithms using the conventional alternating projections of wavefronts between the object and sensor planes and our interpolation at the sensor plane for the pixels in the area S D \S S . In order to demonstrate the importance of this filtering, in what follows, we will also present results obtained with no filtering. To make a difference between the algorithms with and without filtering we name the latter one as the Super-Resolution Gerchberg-Saxton (SR-GS) algorithm.

Experimental setup
For the biological samples investigation, we modified the experimental setup that was presented in [31]. The main reason of this modification is a shortening propagation distance d between Fig. 2. Experimental setup. The laser is a light source with wavelength λ = 532 nm, L 1 , L 2 are beam expanding lenses, BS is a beamsplitter, SLM stands for Spatial Light Modulator, L 3 , L 4 are lenses in a 4f-telescopic system, Object is an investigated specimen, CMOS is a registration camera. object and registration plane to fulfill the condition of Eq. (5). Figure 2 demonstrates the modified setup: laser LASOS GLK 3250 TS with the wavelength λ = 532 nm, digital 8-bit CMOS camera EVS 3264 × 2448 with pixel-pitch ∆ S = 1.4 µm, and phase only reflective SLM Holoeye LETO 1920 × 1080 with pixel-pitch ∆ SL M = 6.4 µm and fill factor 93%. The light beam emitted by the laser source expands in the two lenses L 1 and L 2 collimator and goes through the beamsplitter (BS) to SLM, where the light wavefront obtains the desired phase retardation corresponding to the phase masks M s . Next, the wavefront is transfered to the object plane by 4f system (lenses L 3 and L 4 ), where it interacts with an object and after the wavefront freely propagates to CMOS camera for registration. In comparison with previous configuration [31], we added 4f system (lenses L 3 and L 4 ) for transferring masks to the object plane, thereby the optical distance between object plane and CMOS is shortened to 1.44 mm, instead of 21.55 mm. According to [31], we implement the phase modulation masks M s on SLM with the large super-pixels of the size of 8 × 8 of the SLM's pixels. The super-resolution up-sampling factor r s is equal to 4, therefore, computational pixel size ∆ c = 0.35 µm. According to Eq. (5) for the given values of ∆ c , d, and λ the zero-padded object wavefront should have as a support 6200 × 6200 pixels.

Samples preparation
Epithelial cells used in the experiments are taken from intact oral cavity mucosa of a 33 year old volunteer having no signs of inflammation. Clean cotton swab was used to scrap the intact oral cavity mucosa of the volunteer mouth, then the swab was placed on a coverslip, where after its removal, oral cavity environment containing epithelial cells left. After it was dried in a room conditions for 30 minutes, it was placed under investigation. No additional preparation was performed on the epithelial cells.
The blood sample with erythrocytes was taken from the peripheral vein of the same volunteer by dropper and placed on the coverslip. It was dried in a room conditions for 30 minutes and placed under investigation. No additional preparation was performed on the blood sample.

Results and discussion
Two biological objects were investigated in the paper: erythrocyte and cell of human oral cavity epithelium, under different exposure times (t ex p ) of the camera. In the Poissonian noise distribution, the noise variance is proportional to the exposure time, therefore, variation of the  Figure 3 shows fragments of the registered intensity distributions with the corresponding histograms of different exposure times for the epithelial cell sample. Top row -observed raw intensity distributions obtained from a single experiment (s = 1), from left to rightt ex p =1 ms, 5 ms and 25 ms (pay your attention to colorbar scales, which demonstrate intensities of different levels); bottom row -the corresponding histograms. Case of t ex p = 25 ms corresponds to the highest SN R level in the observations, SN R 25ms = 16.8 dB, and the lowest SN R is almost two times smaller with SN R 1ms = 8.8 dB for t ex p = 1 ms case. Histograms show used dynamic range of the camera, and as it can be seen from Fig. 3, there are full range 8 bit, less then 7-bit, and less then 6-bit for 25 ms, 5 ms, and 1 ms exposure cases, respectively. Figure 4 illustrates fragmented amplitudes (top row) and phases (bottom row) of the human cheek epithelial cell reconstructed by SR-SPAR from experiments with different exposure times t ex p = 1 ms, t ex p = 5 ms, and t ex p = 25 ms, from left to right, respectively. For all reconstructions the borders of the cell and its inner structure with nucleus are clearly seen both in amplitudes and phases. Phase distribution maps show that phase delay in a region of cell nucleus is significantly higher than in other parts of the cell, that corresponds to higher refractive index of the nucleus, knowing that thickness of the cell changes slightly [10]. The reconstruction results for amplitudes for all exposure times are similar, but for phases there are differences between t ex p = 1 ms and others cases. This difference is cosed by the low value of SN R 1ms = 8.8 dB, but even in such noisy case SR-SPAR produce reliable result that can be seen from middle cross-sections plotted in Fig. 5. These cross-sections illustrate phase distribution in reconstructed  phases for experiments with different exposures: t ex p = 1 ms -blue circles line, t ex p = 5 msred diamonds line, and t ex p = 25 ms -green squares line. For all cases cross-sections behavior is almost the same, and the peak values for nucleus are nearly identical.

Epithelial cell
A serious difference in the width of the peaks concerns only the case t ex p = 1 ms, where the peak width is essentially smaller. This difference can be explained by the data quantization in intensity registration by the sensor. As it can be concluded from the histogram images in Fig. 3 the dynamic range of the sensor is about 6 bit while it is much larger for the larger exposure time. The smaller dynamic range results in a narrower boundary imaging. The ability of SR-SPAR to enable high quality imaging for the low-bit sensors is an interesting feature of this algorithm. For instance, such small errors in phase retrieval with lens optical setup in paper [26] is obtained  only with the camera dynamic range equal or more than 12 bits.
Nevertheless, since the peak value for noisy case of t ex p = 1 ms coincides with values for less noisy cases, it is possible to use the SR-SPAR algorithm even in such noisy conditions for cells visualization and refractive index and/or thickness of the cell estimation.
In order to illustrate image improvement due to filtering in the SR-SPAR algorithm, we present the results obtained by the SR-GS algorithm where as it already discussed above the filtering is omitted. Comparison of the phase reconstructions is produced for the nosiest case with t ex p = 1 ms. Figure 6 shows SR-SPAR and SR-GS phase reconstructions obtained for the number of experiments L = 6 and L = 18 and plot of middle cross-sections of these reconstructions. Two images for L = 6 (first and third) are obtained by SR-GS and SR-SPAR algorithms, respectively, where the parameters of the filters for SP-SPAR are optimized for the best results. Comparing these two images we may realize that indeed the filtering results in a dramatical improvement of imaging providing clear information on the general shape and its small variations. The first image in Fig. 6 is given by SR-GS for L = 6 without filtering. The quality of this image is quite poor: the inner structure of the cell and the phase variations are corrupted to the range higher than 4 radians, while this range in the second image is lower. The second image is obtained by SR-GS but for much larger number of experiments L = 18. Comparing this image with the third one, that was made by SR-SPAR algorithm with L = 6, we may note that they are nearly identical with clear boundaries and the same range for the phase. Corresponding phase cross-sections in Fig. 6 demonstrate erroneous reconstruction of the SR-GS method with number of experiments L = 6 (red crosses line) in comparison with other two reconstructions by SR-GS L = 18 (green diamonds line) and SR-SPAR L = 6 (blue triangles line). Therefore, the filtering in the SR-SPAR algorithm allows to decrease number of experiments to L = 6 without losses in quality of reconstructions. Furthermore, less number of experiments decreases the illumination

Erythrocyte
In microscopic scales, the epithelial cell is quite a big sample with size of about 60 µm and the super-resolution advantage is not obvious. For a more descriptive super-resolution demonstration, we took a 8 times smaller biological object -erythrocyte, which typical size is about 8 µm. During reconstruction strong halo effect arose and for it suppression we used Hilbert transform approach for filtering it out [37]. Figure 7 shows fragmented amplitudes (first column) and phase maps (second column) of the erythrocyte reconstructed by SR-SPAR with the pixel resolution ∆ c = 1.4 µm (top row) and the super-resolution ∆ c = 0.35 µm (bottom row). In the pixel resolution, it is hard to see details of the erythrocyte inner structure, what could result in errors to distinguish the types of erythrocytes (e.g. discocyte, spherocyte or echinocyte) or the pressure conditions (orhypotonic and hypertonic) [38]. Presented type of the erythrocyte is hypertonic discocyte as it has a disk shape with a zone of central smaller thickness and corresponding deformations usually observed for hypertonic cases [38]. The super-resolution allows to achieve a much better imaging.

Parameters of the SR-SPAR algorithm
The performance of the SR-SPAR algorithm essentially depends on tuning of its parameters. The following values are selected in our experiments. The image patches in BM3D are square 8 × 8 pixels. The group size is limited by 39 patches. The step size between the neighboring patches is equal to 3 pixels. The discrete cosine (for patches) and Haar (for the group length) transforms are used for 3D group data processing in BM3D [31]. The thresholding parameters in BM3D are taken very large T h b = 100, T h ϕ = 80, as compared with the standards used in [27]. The parameters defining the iterations of the algorithm are as follows: γ is adoptable and proportional to exposure time, number of iterations T = 30. For reconstruction we use MATLAB R2016b and the computer with the processor Intel® Core™ i7-3770 CPU @ 3.40 GHz, 32 Gb RAM. We characterize the complexity of the algorithm by the time required for processing. For the number of experiments L = 6, the super-resolution image in the computational pixels 4896 × 4896 and this image zero-padded to 6200 × 6200 , this time for a single iteration is equal to 75 sec. Thus computationally, the algorithm is quite demanded manly due to high memory size required for calculations.

Conclusion
Application of the SR-SPAR algorithm as a tool for computational super-resolution microscopy is demonstrated. A complex-valued transparent biological specimen can be restored from multiple intensity diffraction patterns. The algorithm is developed for lensless optical setup with multiple exposure experiments and random phase-modulation of wavefronts at the object plane.
The filterings in the observation domain and in the specimen variables are important features of the algorithm. If these filterings are omitted the algorithm becomes a super-resolution version of the standard iterative phase retrieval algorithms. In comparison with this simplified algorithm with no filtering, our algorithm demonstrates a valuable improvement in imaging with smaller number of observations and shorter exposure time. In particular, it is shown for for epithelial cell specimen, that this filtering allows to decrease a number of experiments from 18 to 6 with the same quality of imaging. It is demonstrated also for epithelial cells that the proposed method provides good reconstruction results even from intensities distribution with high noise level 8.8 dB of SN R and for dynamic range of the camera less than 6 bit.Therefore, in such conditions, illumination radiation could be reduced for minimization of radiation influence on the biological cellular specimen.
The effects of the super-resolution are demonstrated for erythrocyte phase and amplitude imaging. With the super-resolution up-sampling factor 4 the computation resolution is ∆ c = 0.35 µm, what is four time better than the sensor size resolution. In comparison with the wavelength λ = 532 nm, we way note that the achieved resolution better than the wavelength resolution.
SR-SPAR algorithm reconstructed high quality amplitude and phase images for both epithelial cell and erythrocyte. Comparison of super and pixel resolutions is presented. Hence, lensless phase retrieval SR-SPAR algorithm was successfully applied for biological samples investigation and can be considered as an alternative tool to traditional microscopic systems.

Disclosures
The authors declare that there are no conflicts of interest related to this article.