Adaptive pixel-super-resolved lensfree in-line digital holography for wide-field on-chip microscopy

High-resolution wide field-of-view (FOV) microscopic imaging plays an essential role in various fields of biomedicine, engineering, and physical sciences. As an alternative to conventional lens-based scanning techniques, lensfree holography provides a new way to effectively bypass the intrinsical trade-off between the spatial resolution and FOV of conventional microscopes. Unfortunately, due to the limited sensor pixel-size, unpredictable disturbance during image acquisition, and sub-optimum solution to the phase retrieval problem, typical lensfree microscopes only produce compromised imaging quality in terms of lateral resolution and signal-to-noise ratio (SNR). Here, we propose an adaptive pixel-super-resolved lensfree imaging (APLI) method which can solve, or at least partially alleviate these limitations. Our approach addresses the pixel aliasing problem by Z-scanning only, without resorting to subpixel shifting or beam-angle manipulation. Automatic positional error correction algorithm and adaptive relaxation strategy are introduced to enhance the robustness and SNR of reconstruction significantly. Based on APLI, we perform full-FOV reconstruction of a USAF resolution target (~29.85 mm2) and achieve half-pitch lateral resolution of 770 nm, surpassing 2.17 times of the theoretical Nyquist–Shannon sampling resolution limit imposed by the sensor pixel-size (1.67µm). Full-FOV imaging result of a typical dicot root is also provided to demonstrate its promising potential applications in biologic imaging.


Materials and Methods
Experimental Setup. Figure 1(a) depicts the configuration of the lensfree imaging setup. The coherent or partially coherent light irradiates the specimen, and then the scattered light and the transmitted light co-propagate in the same direction, finally forming interference fringes on the imaging device. In order to make the emitted light impinging on the object plane to be considered as a plane wave, the distance between the source and samples should be by far larger than that between the sensor and samples. Furthermore, the length scale of the distance between samples and the imaging device typically is on the order of submillimeter 34 , and finally this structure can achieve the whole active area of the imaging sensor as the FOV and reduce the unacceptable resolution loss of the reconstructed images while the magnification (F) approaches unit [see Fig. 1(a), Z 2 >> Z 1 and F = (Z 1 + Z 2 )/Z 2 ≈ 1]. Based on this, the scale of the FOV will be directly restricted by the number of pixels and the pixel-size of the imaging device, and at the meantime the latter will be the main limiting factor of improvement of the spatial resolution. Unfortunately, during the process of Z-scanning which is implemented to achieve super-resolution, whether electric or manual adjustment of the sample-to-sensor distance will result in the tiny lateral positional error inevitably and the co-propagated light beam carrying the error will be sampled by the imaging device. As shown in Fig. 1(b), the tiny error will result in the subpixel shift among the intensity images on the different planes.
As depicted in Fig. 1(a), our lensfree in-line digital holographic imaging system mainly contains three parts: a single mode fiber-coupled light source (LP660-SF20, Thorlabs, the United States), a monochrome imaging device (DMM 27UJ003-ML, the imaging source, Germany), and thin specimen placed above the imaging device. In our experimental system, the optical fiber is put at ~20 cm over the samples. At the meantime, the imaging sensor is placed ~400-900 μm away from the sample which is attached to a piezo-driven positioning stage (MAX301, Thorlabs, the United States). The stage holds the specimen with the self-designed 3D-printed support and can move vertically to change the distance between the sample and the image sensor. The camera in our lensfree in-line digital holographic imaging system has 1.67 μm pixel-size and 10.7 megapixels, and it is used for the acquisition of ten holograms with different sample-to-sensor distances Z 1 . Theoretically, the FOV can reach ~29.85 mm 2 and simultaneously the resolution is up to 2.5 fold of camera resolution. Finally, the spatial resolution of the experimental results is enhanced to the 2.17 times of the theoretical Nyquist-Shannon sampling resolution and the space-bandwidth product (Megapixel) is increased from 10.7 to 50.34.

Sample preparation. A standard 2′′ × 2′′ positive 1951 USAF resolution test target (Edmund Scientific
Corporation, Barrington, New Jersey, USA) is used to quantitatively demonstrate resolution improvement. Besides the resolution target, the typical dicot root (Carolina Biological Supply Company, Burlington, North Carolina, USA) stained with fast green and the counterstain safranin is a representative sample for the study of the internal structure of plants, and applied to demonstrate universality of the proposed method.
Adaptive pixel-super-resolved lensfree imaging (APLI). In order to recover a super-resolution image of the complex object field based on a series of the out-of-focus low-resolution holograms in the spatial domain, the overview flowchart of our method is shown in Fig. 2 which is mainly composed of the following three stages. Stage 1: Generation of an initial guess. A stack of the holograms (e.g, the pixel dimension of the hologram is m × n) is captured on different sample-to-sensor planes and the first plane should get as close as possible to the sensor. After capturing the raw images, up-sampling will be carried to all holograms with the nearest neighbor interpolation which coincides with the imaging theory of cameras, and all the up-sampling images [e.g, the pixel dimension of the hologram is M × N with the interpolation weight k (M × N = km × kn)] will back-propagate to the object plane with the auto-focusing algorithm 35,36 . All the up-sampling intensity images are superimposed together to acquire a good initial guess which will be used as the input of Stage 2.
Although the single back-propagated up-sampling hologram can be regarded as the initial guess, simply summing up all back-propagated up-sampling holograms can significantly suppress the twin image noise, aliasing signals and up-sampling related artifacts 31,32 . Furthermore, with the same set of raw data, this initialization method can have resolution improvement of the initial guess compared to the previous initialization method 14 and then have faster convergent rate in the Stage 2.
Stage 2: Iterative multi-height images reconstruction. The whole iteration process is a procedure for phase retrieval, and it is essentially a process of solving the inverse problem which is very common in the computational imaging. To solve this problem, we need to build a precise forward model and reconstruct the super-resolution intensity images and the phase map from the captured discretized intensity images. The above-mentioned initialization is input into the precise forward model to obtain the estimated captured images. If the estimated captured images cannot match the corresponding captured images, the iterative loop in the model will continue. The loop will terminate until the estimated captured images can accord with the actual captured images, and the current assumptions will be regarded as the actual super-resolution intensity and phase images.
Thus, the following two key elements must be taken seriously to make the model to be consistent with actual physical process. Firstly, in order to obtain the precise imaging model, the limited sensor pixel-size and the unpredictable disturbance during image acquisition should be taken into account. The pixel binning as the process of recording the image is a down-sampling procedure which can be regarded as the spatial averaging. On the other hand, to acquire the diffraction patterns on the distinct planes, the longitudinal shift of the samples is inevitable which will lead to the accidental lateral displacement. Hence, the lateral positional error must be absorbed in the model. (More details are given in the Section Automatic lateral positional error correction). Secondly, the stability and the robustness of the solution to the phase retrieval problem must be improved, and the algorithm should be able to converge to a desired optimal solution that can be considered as the optimization of phase recovery based on multi-height measurements. Solving the optimization problem is deemed to make the current estimate close fit the input actual captured images as a whole, and the quantification is given by the real-space error as described in the following equation where . is the Euclidean norm, I i is the amplitude of i th captured image, | |g i is the current down-sampling estimated amplitude corresponding to the i th measurement, which is the output after consideration of system uncertainties. The solution process is an incremental gradient optimization process, which unfortunately will provide a relatively correct solution in early iterations, but then overshoot. This problem is often attributed to the non-convex nature of phase retrieval, but we find the reason for this is more closely related to the choice of the relaxation factor based on the analysis of the similar problem in pre-work 37 , and the relaxation factor needs to be gradually diminishing for convergence even in the convex case. Thus, the adaptive relaxation factor should be introduced into the model to achieve the improvement in the stability and robustness of the reconstruction towards noise. (Details can be referred to the Section Adaptive relaxation factor).
The step of updating amplitude is crucial and depends on a correction coefficient matrix. This matrix is determined by the product of the adaptive relaxation factor α and the proportional relation matrix between the up-sampling captured images and the prior estimated intensity images. The specific process can be seen in Fig. 3 and described as the following three steps.
Step 1, the (i-1) th estimated complex amplitude − O i j 1 is forth-propagated to the next height (O i j ) and then the i th captured images are up-sampled (j represents the current index of the iteration cycle and i represents the index of the out-of-focus plane.). Next, the estimated intensity image | | O i j 2 is registered with the up-sampling captured images I upsample , and the positional error is represented by (x shift , y shift ). Then we shift | | O i j 2 in place by the amount of (−x shift , −y shift ) and the original estimated intensity image is substituted by the refined intensity termed (the amplitude can be denoted as A).
Step 2, we implement down-sampling to the refined estimated super-resolution intensity images | | O _

38
, where a h is the gray value of the super-resolution intensity images, k is a down-sampling factor.
Step 3, after down-sampling, the estimated low-resolution intensity image has the same dimension with the original captured image on the corresponding i th sensor-to-sample plane. Then we up-sample the estimated low-resolution intensity image and the corresponding captured image (the amplitude can be referred as the matrix B and C respectively.) with the nearest neighbor interpolation. To acquire the correction coefficient matrix, we multiply the adaptive relaxation factor α with proportional value between the matrix C and B, and the expression will be regarded as the updated i th estimated amplitude. The relaxation factor α in above expression is a diminishing value differing from the traditional fixed value ~0.5, and the guided filter is taken into account to further eliminate the influence of noise as well. At last, the complex amplitude containing the new updated estimated amplitude and the previous unchanged phase is forth-propagated to the next height using the Angular Spectrum Method 39 .
The process of the Step 1-Step 3 is repeated until all the sample-to-sensor distances are gone through. That is to say, all the raw measurements are used for once, and it will be considered as one iteration cycle.
Stage 3: Reconstruction on the object plane. After some iterations, we will achieve the complex amplitude on the plane closest to the imaging device and then back-propagate the complex amplitude to the object plane as shown in Fig. 2.
Physical modeling of the pixel binning. Blurring may be caused by an optical system (inherent noise of the camera, diffraction limit, etc.), and the PSF of the imaging device. The former can be modeled as linear space invariant while the the latter is considered as linear space variant 38 . It is difficult to obtain the exact information about the linear space invariant, so it is usually compensated by the specific algorithms or avoided as much as possible. Besides the linear space invariant, in the process of image reconstruction, the PSF of the imaging device (which can also be regarded as the finiteness of the physical pixel-size) is an important factor for blur, which should be incorporated into the reconstruction procedure. As a complementary interpretation, there is a natural loss of spatial resolution caused by the insufficient sensor density and noise that occurs within the sensor or during transmission. As shown in Fig. 4(b), the spectrum loss will be more serious while the decimation factor increases. Figure 4(b) indicates again that the pixel-size is the main limiting factor of the systems which will determine whether it can directly record the high frequency fringes corresponding to the super-resolution of the samples.
In traditional multi-height reconstruction method 23,33 , many efforts are made to implement subpixel shift to achieve the super-resolution, but the pixel binning is not taken into account, which does not accord with actual physical process. Thus, involving the process of recording digital images in the reconstruction procedure has drawn attention, and the enhancement in resolution has been validated 31 . The process of recording digital images is a down-sampling process which is usually modeled as a spatial averaging operator where a h is the gray value of the super-resolution intensity images, k is a decimation factor] as shown in Fig. 4(a). In the iterative process of the reconstruction, we convolve the estimated intensity of the field | | O i j 2 with the PSF of the image sensor 38 , and then it has the same dimension as the raw measurement. Automatic lateral positional error correction. In many lensfree systems, the different sample-to-sensor distances are provided by the mechanical movement, so the mechanical lateral error will be generated as shown in Fig. 1(b). In order to eliminate the unavoidable error at subpixel-scale, the traditional method is registering images before reconstruction 23,24 called beforehand lateral positional error correction (BLPEC). However, in the actual imaging process, the light illuminates the specimen which has tiny lateral movement because the positioning stage will lead into the lateral mechanical error while it moves longitudinally, and then forth-propagates to the imaging plane carrying the information of the object. Thus, the lateral positional error appears before camera sampling and the captured images carry the error signals. Based on this, the BLPEC can only correct the lateral positional error cursorily, and the accuracy of the correction will decrease due to existence of the artifacts and aliasing. Additionally, this method has incapacity to rectify registration error in the later process which will affect the final quality of reconstruction.
In order to solve the problems existing in traditional BLPEC, we introduce automatic lateral positional error correction (ALPEC) into our method. The crux of solving the general problem of subpixel image registration is computing the cross correlation between the image to register and a reference image by means of a fast Fourier transform (FFT), and locating its peak 40 . The cross correlation of the captured image f(x, y) and its corresponding estimated image g(x, y) is defined by: where M and N are the image dimensions, (*) represents complex conjugation, F and G denote the discrete Fourier transform of the f and g respectively. The expression of and there is a similar expression for G (u, v). It is important to determine accurately the peak of the cross-correlation function r fg (x shift , y shift ), and relax the limitation on computational speed and memory caused by the FFT. Thus, the refined initial estimate method 40,41 is used, which uses aid by the existence of analytic expressions for the derivatives of r fg (x shift , y shift ) with respect to x shift and y shift , and the algorithm iteratively searches for the image displacement (x shift , y shift ) that maximizes r fg (x shift , y shift ). At last, it can achieve registration precision to within an arbitrary fraction of a pixel at a fast rate.
Adaptive relaxation factor. In most cases, the propagation phasor approach 31 is an effective solution to the pixellation problem and it gives a unified mathematical framework combining phase retrieval and pixel super-resolution. Nevertheless, in practical operation, the stability and reconstruction quality of the method may be significantly degraded due to the existence of nonnegligible noise during the sampling process. This problem is often attributed to the non-convex nature of phase retrieval and the ill-condition process of the super-resolution reconstruction. Although numerous super-resolution algorithms have been proposed in the literature 11,23,25,[42][43][44] , the super-resolution image reconstruction remains extremely ill-posed 43,44 . Moreover, the noise effect will accumulate as the iterations augment.
The choice of the relaxation factor can suppress noise to a certain extent, typically the relaxation factor α = 0.5 in traditional phase retrieval methods. The relaxation factor will be utilized to update amplitude, and the corrected super-resolution intensity images substitutes for the earlier estimated super-resolution intensity images incompletely. In other words, the corrected super-resolution intensity images will occupy a part in the new updated estimated intensity (as shown in Fig. 3), which have a close relationship with the captured intensity images containing noise. So the new updated estimated intensity images suffer from the noise because the captured intensity images impose ill-condition restrictions with the fixed relaxation factor.
The stability and reconstruction quality may be significantly degraded when non-negligible noise is present in the captured images, and the same problem is encountered in FPM. We find that the reason for the phenomenon in this field is the non-convex nature of phase retrieval and more closely related to the choice of the step-size, so the adaptive step-size strategy is introduced to successfully solve this problem 37 . Considering that the problems of the iterative method in lensfree imaging are also attributed to the non-convex nature of phase retrieval, instead of the traditional fixed relaxation factor, the adaptive relaxation factor which diminishes to an infinitely small value will be used to improve the performance of the incremental solutions. So the critical issue in practical application will be how to determine a suitable relaxation factor sequence α iter to get close to a solution within fewer iterations. The α iter must satisfy the two conditions that are shrinking the relaxation factor to zero and making the diminishing speed not be too fast. Because if the relaxation factor shrinks too fast, the estimated object field may converge to a point that is not a minimum especially when the initial point is sufficiently far from the optimum. So the relaxation factor should not be reduce too fast and then the algorithm can travel infinitely far. Thus, in this paper, we give an alteration of the relaxation factor when the global error ε(O ite−1 ) and ε(O iter ) obtained in consecutive cycles satisfies the following criterion:  The super-resolution image needed to recover in theory, (b) the low-resolution image captured by the camera on the object plane theoretically, (c,d) the images captured by the camera on the planes of distinct sample-to-sensor distances with noise and not respectively, (e,f) the reconstructed super-resolution images using a fixed relaxation factor (α = 0.5) under the noise-free and noisy circumstances separately, (g,h) the reconstructed super-resolution images using the adaptive relaxation factor under the noise-free and noisy circumstances separately. system uncertainties such as the lateral positional error. Finally, the algorithm will converge to the stationary point when the relaxation factor reach a pre-specified minimum.

Discussion and Results
The comparison between the adaptive and fixed factor. Figure 5 shows the influence of the noise on the system and emphasizes the important role played by the relaxation factor in iterative process. A theoretical super-resolution image needed to reconstruct is shown in Fig. 5(a), and Fig. 5(b) shows the low-resolution image captured by the camera on object plane in theory. Figures 5(d) and (c) describe a set of emulational images captured by the camera on the planes of distinct sample-to-sensor distances with Gaussian noise and not respectively. Figures 5(e) and (f) depict the reconstructed super-resolution images using the fixed relaxation factor ~0.5 under the noise-free and noisy circumstances separately. The two yellow curves in the sub-graphs of red region convey the information that if we update the earlier estimated super-resolution intensity images using a fixed relaxation factor (α = 0.5) as the new estimated amplitude in the above-mentioned Section APLI, the reconstructed result under the noisy condition [ Fig. 5(f)] is much worse in respect of the resolution and background after the same iterations compared to the result without noise [ Fig. 5(e)]. To demonstrate that an adaptive relaxation factor can effectively solve the problem having a close relationship with the over-amplification noise, we test our method under the noise-free and noisy circumstances separately. The results can tell apart the densest line and give relatively clean background under the two different conditions as shown in Figs. 5(g) and (h) respectively. Furthermore, the quantitative comparison of reconstruction accuracy versus intensity noise among using the adaptive and fixed relaxation factors (α = 0.5,0.01), as well as the rate of convergence is shown in Fig. 6. Figure 6(c) depicts the reconstructed results corresponding to iterations labelled in Fig. 6(b) under the condition of the adaptive and fixed relaxation factors respectively. Figure 6(a) shows the curves of the intensity error following the iterations increasing with different relaxation factors and Fig. 6(b) shows the local enlarged drawing of Fig. 6(a) (shaded region). Among the curves, the purple bight represents that using the perfect initial guess and a very small fixed relaxation factor ~0.01, the intensity error still accumulates as the iterations increase. This offers an explanation for the following phenomenon that with the fixed relaxation factors, the reconstruction error will have convergent tendency at the outset and then get worse after reaching their respective minima which can be obviously seen in the green curve of Fig. 6(a). The same goes for the small relaxation factor corresponding to the red curve, but it is not obvious to observe the the turning point of the curve (the red curve reaches the minimum in the about 700 iterations and then overshoots), because the speed of convergence is extremely slow and the curve rises at a glacial pace after reaching the minimum. The iteration should be suspended when the curves reach their respective minima due to the overshooting of the curve in the later period. The cause of the overshooting is that even if the reconstruction converges to a true value, the captured images still provide the ill intensity constraints as before. Non-convergence is a disadvantage for the iterative methods, and suspending the iteration when reaching the minimum will result in loss of image details or taking a long time. Comparing the orange curve with the green one or the red one, we can find that using adaptive relaxation factor can obtain the converged reconstruction and effectively prevent the overshooting. Meanwhile the introduction of the adaptive relaxation factor into our method can retain the relatively fast initial convergence speed and it is seen that this method decreases more rapidly than fixed relaxation factor methods (α = 0.01) and converges in the early 20 iterations.
Although the adaptive relaxation factor has the anti-noise capability to a certain extent, the noise will cause the estimated intensity image to have no tendency to be consistent one with another on the next plane and the reconstructed image to deviate from the theoretically calculated values. To avoid the over-amplification of noise we take into account nonlinear denoising algorithm termed guided filter 45 . It is essentially equivalent to add the relevant transcendental knowledge that objects are piecewise smooth. The introduction of the guided filter will further restrain the noise and the smooth regions of the reconstruction results will tend to the ideal value. However, the reconstruction results are slightly flawed at edge, because the guide filter cannot distinguish whether it is noise or jump edges of the object and preserves the edges during the reconstruction process. The combination of the adaptive relaxation factor and the guided filter can effectively suppress the noise and achieve better reconstruction results corresponding to the blue curve. From these results, we can safely conclude that the adaptive relaxation factor method outperforms the fixed relaxation factor methods, with both faster convergence rate and lower mis-adjustment error simultaneously achieved. Figure 7 shows that the mechanical lateral positional error has great effect on the reconstruction results and lateral positional error correction significantly improves reconstruction performance. The theoretical super-resolution image needed to reconstruct is shown in Fig. 5(a).  (1) The down-sampling by a factor of four is implemented which can also be regarded as the spatial averaging operating weight of the sensor (the average value of sixteen pixels in the super-resolution image is equivalent to the value of the corresponding pixel in the low-resolution image). The actual physical phenomena and process is that the super-resolution image propagates in the free space, and then the two-dimensional continuous intensity distribution of the hologram (diffraction patterns) is discretized into a matrix through the two-dimensional convolution of the hologram and a pixel unit in an imaging array, which results in the low-resolution image. (2) In each group of simulations, variable-controlling approach (for instance, the number of iteration remains unchanged) is used to make the conclusion more convincing. (3) In every simulation, at least 32 raw low-resolution images (the number of pixels is m × n) are utilized to reconstruct super-resolution images [the number of pixels is M × N (M = k × m, N = k × n, k = 4)]. These captured diffraction patterns are enforced as object constraints, gradually converging to the missing two-dimensional phase information 46 and the corresponding super-resolution amplitude. For a complex intensity object function, to obtain the super-resolution intensity, the recovery problem becomes undetermined by a factor of 2 since there are 2 × M × N pixels defining the object function (M × N pixels for the real part and M × N pixels for the imaginary part), whereas there are only m × n pixels in the measurement matrix 47,48 . In order to solve this underdetermined problem, more information about the object function needs to be acquired and incorporated as a constraint on the solution space, so at least 2 × k × k raw low-resolution images are needed in theory 47 .

The experimental results of the USAF resolution test target. A standard 1951 USAF resolution test
target as the experimental samples is utilized to prove that our method has the universality and stability during the actual measurements. In order to test our method, we acquire 10 raw holograms at different sample-to-sensor distances (~547-577 μm) with the standard 1951 USAF resolution test target and each raw hologram is digitized by the imaging device with 1.67 μm pixel-size. Figure 8(a) shows a full FOV (~29.85 mm 2 ) low-resolution hologram which is captured by the camera directly. The inset shows local enlarged drawing of the dashed rectangular area in Fig. 8(a) which corresponds to the full FOV of 10X objective lens. Due to the relatively large pixel-size resulting in down-sampling, our method is applied to diminish the effective pixel-size namely improving the resolution. During the process of reconstruction (the image-processing steps are depicted in the Section APLI), the raw holograms are used as the intensity constraints and the recovered super-resolution intensity image is shown in Figs. 8(b) and (c). In Fig. 8(c), we can deduce that the smallest resolved half-pitch can reach 0.77 μm, which Although there is a non-negligible phenomenon that the imaging results with a 10X objective and Kohler illumination (condenser aperture wide open) outperforms the proposed method, the space-bandwidth product of the reported method has been increased by nearly 100 fold. Furthermore, the whole system requires no lenses, which provides the possibility of miniaturization and low cost. In Supplementary Video 1, we show a zooming video of the full-FOV reconstructed images of a USAF target with our method and the traditional method 23,33 (with BLPEC) respectively. To process the data in parallel, the large format raw image (3872 × 2764 raw pixels) is divided into 35 portions (700 × 700 raw pixels) for computation. Here, the blocks at the end of each row or each column do not have the same pixels as others and the adjacent portions have 200 pixels overlap with each other. For the reconstructed image (2800 × 2800 pixels), we cut away 200 pixels at the edge and the adjacent portions introduce a certain degree of redundancy (400 pixels) into our stitching. Thus, no observable boundary is present in the stitched region and the blending comes at a small computational cost of redundantly processing the overlapping regions twice. Figure 9 describes additional experimental work to address the significance of ALPEC, and intuitively shows the comparison between reconstructed results based on the adaptive and the fixed relaxation factor. Figure 9(a) presents the FOV of the USAF resolution target recorded by the camera directly. Figures 9(a1) and (a2) show the local enlarge enlargements of the rectangular areas in Fig. 9(a) and the respective reconstruction results are shown in Figs. 9(b1)-(c6). To illustrate great effects of ALPEC on experimental results, we have carried out two groups of experiments with the adaptive and the fixed relaxation factor respectively. Figures 9(b1) and (c1) show the reconstructed results of enlargements of two different small segments in Fig. 9(a) based on the single raw image, and as shown in them, the reconstructed results are blurry. Figures 9(b2) and (c2) show the reconstructed results without positional error correction. From Fig. 9(b2) we can deduce that the tiny lateral positional error has little effect on low frequency of the reconstructed results, but exerts a tremendous influence on high frequency which corresponds to the super-resolution [ Fig. 9(c2)]. The same data set is employed to recover the super-resolution image for each segment using our adaptive method with BLPEC and ALPEC respectively. In addition, the same iterations are conducted in two reconstruction methods, the only difference between the methods utilized in this paper is that the latter involves ALPEC while the former puts into effect positional error correction in advance. Figures 9(b3) and (c3) present the recovered super-resolution intensity images with BLPEC corresponding to the same segment of Figs. 9(b1) and (c1) respectively. It can be seen that the silhouette of lines in Fig. 9(b3) is unsharp but recognizable because the corrected images using BLPEC have removed the relatively large lateral positional error, but still remain tiny the positional misalignment. Even worse, due to the remnant lateral positional error, the higher frequency of the object is still unable to recover as shown in Fig. 9(c3), although the resolution of the reconstructed results has been improved compared to Fig. 9(c2). With the help of ALPEC, high-quality recovered intensity distributions are obtained, as shown in Figs. 9(b4) and (c4). The blur in Figs. 9(b3) and (c3) is eliminated completely and the clarity of image is increased, meanwhile better outlines are given in Figs. 9(b4) and (c4). Similarly, with the fixed relaxation factor (α = 0.5), the reconstructed results are shown in Figs. 9(b5) and (c5) as well as Figs. 9(b6) and (c6) using the BLPEC and ALPEC separately. We can also come to the same conclusion that the ALPEC can bring benefits to improve the resolution with the adaptive relaxation or fixed factor.
In order to experimentally illustrate that the adaptive relaxation factor can improve the stability and robustness of the reconstruction towards noise, the comparison between the adaptive relaxation and the fixed factor based reconstructed results is also shown in Fig. 9. Figures 9(b6) and (c6) show the reconstructed results with fixed relaxation factor (α = 0.5) and ALPEC, while the reconstructed results using our adaptive relaxation-factor method and ALPEC are shown in Figs. 9(b4) and (c4). It is obvious that using the adaptive relaxation-factor method can suppress the over-amplification noise without extra auxiliary information.
Imaging of the typical dicot root. Another experiment was demonstrated that our method can also be used for the dense sample such as plant slice, which can be seen in Fig. 10. Figure 10(a) shows the full FOV of the typical dicot root (~29.85 mm 2 ), and the whole sample can be captured which is challenging for the traditional high-magnification lens microscope. From upper left enlarged region of the orange dashed box in raw full FOV low-resolution hologram [see the inset of Fig. 10(a)], we can find that the details in the typical dicot root are hard to be observed because they are submerged in the diffraction fringes. Figures 10(b) and (c) show the reconstructed intensity images based on the traditional method 23,33 (α = 0.5, BLPEC) and our method respectively. The selected area [ Fig. 10(c)] occupy only 1% of the FOV of Fig. 10(a), which corresponds to the whole FOV with the 10X objective as shown in Fig. 10(d). From Fig. 10(c), it is easy to distinguish endodermis, pericycle, primary phloem, primary xylem, and parenchyma cell, which are extremely important for botanical studies. To observe the details inside the amyloplasts, we further select a small area in Fig. 10(c) which is the full FOV with 60X objective [Fig. 10(d1)] and we can find that the unit magnification lensfree systems greatly expand the FOV. Figures 10(b1) and (c1) are the local enlarge enlargements of the rectangular areas in Figs. 10(b) and (c) separately. As shown in the enlargements [ Fig. 10(c1)], the grains in amyloplasts are distinguishable and the sharp improvements are noticed in the image contrast compared to the results shown in Fig. 10(b1). Figures 10(b2)-(d2) show the line profiles along the respective arrow, and two particles in the middle cortex are easily to distinguish which is impossible in the traditional method. However, there are many horizontal and perpendicular lines in the enlargements [ Fig. 10(c1)], and blur is obvious in Fig. 10(c2). The reasons are mainly that guide filter is sensitive to the smooth background, but this experimental sample is not piecewise smooth and the guided filter brings the aberration to the reconstructed results. Furthermore, the test object has a certain thickness and the diffraction patterns of the non-target objects in the vertical direction of objects on focal plane will influence imaging reconstructed results. In Supplementary Video 2, we show a zooming video of the full-FOV reconstruction result of a typical dicot root with our method and the traditional method 23,33 (α = 0.5, BLPEC) respectively.

Conclusion
In this work, a method termed as APLI is proposed to mitigate the artifacts and simultaneously obtain super-resolution images only with Z-scanning. According to the method more than double pixel resolution of camera is successfully achieved, and there is no extra embedding medium between the object and sensor, like the refractive index matching oil. Here we emphasize that this super-resolution technique does not require lateral displacements, wavelength changing, and illumination angles scanning. Throughout the experiment, an imaging sensor with pixel-count of 10.7 million and pixel-size of 1.67 μm provides a large FOV (~29.85 mm 2 ), and the samples are moved vertically to generate ten out-of-focus undersampling intensity images with artifacts. Instead of the traditional fixed-step, an adaptive relaxation factor strategy has been firstly introduced into our method to suppress the over-amplification noise and retain the convergence speed under noisy conditions. Furthermore, we introduce an ALPEC method into our method which tallies with actual physical process, and it can avoid the misalignment effectively and improve the reconstruction stability. APLI offers a way to exploit the resolution potential of lensfree microscopy and achieves smallest resolved half-pitch of 770 nm, surpassing 2.17 times of the theoretical Nyquist-Shannon sampling resolution limit.
We believe that our method will broadly benefit the lensfree imaging microscope and acquire higher resolution with the same amount of data comparing to the traditional reconstruction methods. In addition, our method can vastly not only remove the adverse impact of alteration in multiple systematic parameters on the reconstructed results, but also reduce the complexity of the actual operation. The results of the resolution target and botanical samples demonstrate that the proposed reconstructed method can offer a new way to make the lensfree microscope to be a competitive and promising tool for the medical care in remote areas in future.
However, some issues still deserve further consideration. Although the number of the captured images may influence the reconstructed results, the artifacts cannot be removed completely due to the trade-off between the resolution and the artifacts. Specifically, if we need to weaken the influence of the artifacts, the sample-to-sensor distance must be increased, but at the meantime the long-distance will result in failure of acquisition of the high-frequency patterns because many patterns become denser suffering from more severe pixel aliasing or exceed the sensor area limits. On the other hand, in the actual situation, the resolution cannot be further improved while the number of raw images increases. We consider the reason for this phenomenon is that during the whole process, the longitudinal error is not taken seriously. In future work, we will make effort to correct the tiny longitudinal error automatically.