Autofocusing Algorithm for Pixel-Super-Resolved Lensfree On-Chip Microscopy

In recent years, lensfree on-chip microscopy has developed into a promising and powerful computational optical microscopy technique that allows for wide-field, high-throughput microscopic imaging without using any lenses. However, due to the limited pixel size of the state-of-the-art image sensors, lens-free on-chip microscopy generally suffers from low imaging resolution, which is far from enough to meet the current demand for high-resolution microscopy. Many pixel super-resolution techniques have been developed to solve or at least partially solve this problem by acquiring a series of low-resolution holograms with multiple lateral sub-pixel shifting or axial distances. However, the prerequisite of these pixel super-resolution techniques is that the propagation distance of each low-resolution hologram can be obtained precisely, which faces two major challenges. On the one hand, the captured hologram is inherent pixelated and of low resolution, making it difficult to determine the focal plane by evaluating the image sharpness accurately. On the other hand, the twin-image is superimposed on the backpropagated raw hologram, further exacerbating the difficulties in accurate focal plane determination. In this study, we proposed a high-precision autofocusing algorithm for multi-height pixel-super-resolved lensfree on-chip microscopy. Our approach consists of two major steps: individual preliminary estimation and global precise estimation. First, an improved critical function that combines differential critical function and frequency domain critical function is proposed to obtain the preliminary focus distances of different holograms. Then, the precise focus distances can be determined by further evaluating the global offset of the averaged, low-noise reconstruction from all backpropagated holograms with preliminary focus distances. Simulations and experimental results verified the validity and effectiveness of the proposed algorithm.

In recent years, lensfree on-chip microscopy has developed into a promising and powerful computational optical microscopy technique that allows for wide-field, high-throughput microscopic imaging without using any lenses. However, due to the limited pixel size of the state-of-the-art image sensors, lens-free on-chip microscopy generally suffers from low imaging resolution, which is far from enough to meet the current demand for high-resolution microscopy. Many pixel super-resolution techniques have been developed to solve or at least partially solve this problem by acquiring a series of low-resolution holograms with multiple lateral sub-pixel shifting or axial distances. However, the prerequisite of these pixel super-resolution techniques is that the propagation distance of each low-resolution hologram can be obtained precisely, which faces two major challenges. On the one hand, the captured hologram is inherent pixelated and of low resolution, making it difficult to determine the focal plane by evaluating the image sharpness accurately. On the other hand, the twin-image is superimposed on the backpropagated raw hologram, further exacerbating the difficulties in accurate focal plane determination. In this study, we proposed a high-precision autofocusing algorithm for multi-height pixel-super-resolved lensfree on-chip microscopy. Our approach consists of two major steps: individual preliminary estimation and global precise estimation. First, an improved critical function that combines differential critical function and frequency domain critical function is proposed to obtain the preliminary focus distances of different holograms. Then, the precise focus distances can be determined by further evaluating the global offset of the averaged, low-noise reconstruction from all backpropagated holograms with preliminary focus distances. Simulations and experimental results verified the validity and effectiveness of the proposed algorithm.

INTRODUCTION
In recent years, high-resolution wide-field optical imaging has become a valuable tool in various biomedical applications such as cell counting [1], cell morphology measurement [2], and optical microscopy techniques that require low-cost and compact imaging systems [3][4][5]. Nevertheless, it is difficult for traditional microscopic imaging technology to satisfy the large field-of-view and high-resolution at the same time. The recently developed lensfree on-chip microscopy provides a new opportunity to process high-resolution and large field-of-view images, which has the advantages of small size, low cost, and excellent performance [6][7][8][9][10][11][12][13][14]. Lensfree on-chip microscopy reconstructed the image of a specimen based on the interference patterns created by the diffracted object. To achieve wide-field imaging, the samples are placed as close as possible to the imaging sensor [13,14]. In this case, the resolution of the reconstructed image is limited by the pixel size of the sensor instead of the system numerical aperture, which results in the system resolution does not meet the needs of current applications. In order to solve this problem, superresolution algorithms that process a series of low-resolution holograms have been applied to lensfree imaging [15][16][17][18]. Many related studies have been proposed and can be classified as: twodimensional movement of light sources by mechanical structures [19], tightly distributed fiber arrays [14], multi-height superresolution algorithm [11], different illumination wavelengths [18], different illumination angles [20], and so on.
Furthermore, to improve the data processing efficiency of super-resolution and twin-image removal, related research has been proposed. Luo et al. [12] proposed a propagation phasor approach, which combines pixel super-resolution and phase retrieval techniques. Zhang et al. [8] proposed an adaptive pixelsuper-resolved lensfree imaging (APLI) method based on the Zaxis scanning holograms. Both of these methods require a series of holograms of different heights. In the process, the sampleto-sensor distances must be accurately determined, otherwise, it will directly affect the final super-resolution image resolution. For example, an error of 1 µm in the reconstruction distance directly results in the absence of information corresponding to group 9 of 1951 USAF resolution test target for typical lensfree onchip systems [8]. Therefore, accurate autofocusing is an essential step in the actual experimental measurement of lensfree superresolution microscopy.
In the past decades, various autofocusing methods have been proposed for digital holography [21][22][23][24][25][26][27][28][29][30][31]. Automatic evaluation function is the primary consideration of image definition, which is mainly determined from the spatial domain and the frequency domain. Compared with the defocused blurred image in the spatial domain, the focused image has obvious sharpened edges. Based on the concept, the Sobel operator [26], Laplace operator [27], and the total sum of gradient [23] have been used as evaluation function. In the frequency domain, the focused image has more high frequency components than the blurred image, so the high frequency component value can be used as the focus evaluation standard, such as bandpass-filtered power spectrum [27]. In addition, evaluation functions based on other models have appeared, such as differential critical function (DIF) [23], self-entropy [28], spectral L1 norms [29] etc. However, the above algorithms are difficult to directly apply to the autofocusing of multi-height pixel-super-resolved lensfree on-chip imaging. The main problems are as follows: (1) Due to the use of lensfree microscopy, the sharp edge of the focused image is not obvious, so it is difficult to judge the focus directly. (2) The lensfree on-chip microscopy is inevitably interfered by twin-image, which affects the accuracy of evaluation function. (3) In superresolution imaging, each hologram of different heights needs to be accurately reconstructed, which requires high accuracy of the focusing distance of each plane. Currently, there is no specific and effective technology to achieve high-precision autofocusing under the influence of twin-image for lensfree on-chip imaging.
In this paper, a high-precision autofocusing method is proposed based on the multi-height pixel-super-resolved lensfree on-chip microscopy. The proposed method can be divided into two steps: individual preliminary estimation and global precise estimation. First, an improved differential critical function (SDIF) has been proposed to obtain the approximate focus distances of different height planes, which combines differential critical function and frequency domain autofocus. In order to solve the problem that the DIF value will remain flat near the peak or valley, we introduced frequency domain autofocus to improve the sensitivity of DIF. Then, the intensity or phase image of different heights is synthesized into a high-resolution reference without twin-image. we regard a series of holograms of different heights as a whole, and a high-resolution reference image can be obtained by consistently varying the reconstructed distances of the holograms at different heights. we transform the problem of focusing on different planes into a problem of solving the optimal solution of the overall system. Last, by searching for the extremum of the synthesized image, we can determine the overall offset and get the final reconstruction distance sequence. Through the proposed method, we can more accurately determine the focus distances of different plane heights and improve the accuracy of super-resolution based on the multi-height pixel-super-resolved lensfree on-chip imaging. The structure of the paper is as follows: Basic introduction of lensfree on-chip holography imaging and theoretical knowledge of proposed methods are carried out in section 2. Simulations are carried out in section 3. The section 4 shows the experimental setup and experimental results for amplitude-contrast or phasecontrast objects. Conclusions are given in the in section 5.

The Basic Principle of Lensfree On-Chip Holography
First, we introduce the basic recording and numerical reconstruction principles of lensfree on-chip holography, which is essentially in-line digital holographic imaging. When the light source with sufficient spatial and temporal coherence illuminates the sample, the hologram H(x, y) is formed by the interference between the scattered light O(x, y) on the sample and the unscattered background light R(x, y), which can be expressed as: Where, r(x, y) and o(x, y) are the amplitude information of the reference light and the object light, respectively, and φ(x, y) is the phase information of the object light wave. Through the angular spectrum algorithm, the light intensity and phase information of the object light wave can be obtained, which can be expressed as: Where F and F −1 denote the Fourier transform and the inverse Fourier transform, respectively. G(f x , f y , z) is the optical transfer function in the frequency domain. λ is the wavelength of the light source, z 1 is the distance from the sample to the sensor plane. For lensfree on-chip microscopy, in order to expand the field of view and improve the resolution, the distance between the sample and the sensor is very small (z 1 ≪ z 2 ) as shown in Figure 1A. The typical distance of sample-to-sensor is 400-1,000 µm. Since the system does not have a microscope objective, it is necessary to propagate the hologram back to the focal plane. In actual operation, the distance between the sample and the sensor is usually obtained by the autofocusing algorithm. Based on that, the maximum field of view is similar to the working area of the image sensor, and the maximum image resolution is approximately the resolution of the image sensor. Due to the limitation of the pixel size of the sensor, the resolution of lensfree on-chip holography is insufficient. Reconstruction based on multiple heights is an effective super-resolution reconstruction method, which can also remove twin-image at the same time. The triangle, ellipse, and pentagram in Figure 1B represent grid points of captured holograms on different planes, respectively, and each of the low-resolution images contains similar and complementary information. The super-resolution reconstruction combines all these new information into one image to obtain a high-resolution image. The most common matrix notations used to formulate the general super-resolution model in the pixel domain can be defined as [32]: Where, the matrix F k is the geometric motion operator between the high-resolution (HR) frame X and the low-resolution (LR) frame Y k . The camera's point spread function (PSF) is modeled by the blur matrix H k , and matrix D k represents the decimation operator. V k is the system noise and N is the number of available LR frames. Based on this, the problem of obtaining a highresolution image from a set of low-resolution images can be attributed to the solution of the following equation: In the process of super-resolution, we need to get a series of low-resolution holograms with different distances of sample-tosensor. However, the reconstruction distance of each hologram is difficult to determine accurately in actual experiments, which directly affects the resolution of the final super-resolved image. Therefore, accurate autofocusing is an essential step in the lensfree on-chip super-resolution imaging.

The Proposed Critical Function
First, we initially determine the reconstruction distances for different height holograms. In the autofocusing algorithm, the key issue is to find a high-precision critical function, which directly affects the subsequent resolution accuracy problem. Considering the influence of twin-image, the differential critical function could be more suitable as the critical function. By making a difference between adjacent reconstructed images, the influence of noise such as twinimage can be reduced to a certain extent. The plane in focus can be determined by searching for the minimum value of differential critical function for amplitude-contrast objects or the maximum value for amplitude-contrast objects. Therefore, the reconstruction distance can be preliminarily determined by searching the extreme value in the initial interval [z 1 , z 2 ] with the step size , as shown in the following formula: where l = 0, 1, ..., (z 2 − z 1 )/ − 1. Ū is the light intensity or phase depending on the type of object. It should be emphasized that for traditional the differential critical function, usually only rely on the intensity. The initial reconstruction distance can be determined by finding the minimum value for the amplitude-contrast object or the maximum value for the phasecontrast object. When the need for focusing accuracy is not very high, the position in focus can be approximated by the above method. However, when the need for focusing accuracy is high, it is more accurate to use intensity as the evaluation data for amplitude-contrast objects, and phase as evaluation data for phase-contrast objects. Although noise interference can be reduced to a certain extent through the critical function of differential, this discriminant factor is insufficient for highprecision focusing sensitivity. Especially when the step size is very small, it is easy to have multiple extreme values. The corresponding simulations and experiments are in the following sections.
In order to solve this problem, the frequency domain critical factor is introduced. The main reasons can be attributed to the following aspects: The critical factor in the frequency domain is more sensitive to high-frequency information. However, if the hologram is directly reconstructed, the high-frequency information generated by the noise will make the frequency domain critical factor unstable or make errors. When the interference of twin-image is weakened, the frequency domain focusing has higher accuracy and stability. Actual simulations and experiments below verify the reliability and stability of the proposed critical factor. Therefore, we can define the criteria function SDIF as shown in the following formula: Where DCT is the Fast Cosine Fourier Transform. M,N are the size pixels of the image. By calculating the value of SDIF at different position intervals, we can preliminarily determine the focus distances of different holograms based on the extreme value of SDIF.

Individual Preliminary Estimation and Global Precise Estimation
According to sections 2.1, 2.2, the reconstruction distances for different height holograms can be initially obtained. However, the resulting reconstructed image still has the effect of the twin-image and it is difficult to further determine the focal distance by a single evaluation factor. Therefore, to ensure the overall resolution of the reconstructed images, we synthesize the reconstructed images of different heights into a single reconstructed image. Based on this, the problem of focusing on different planes is transformed into the problem of solving the optimal solution of the whole system. There are two main issues that need to be addressed. Firstly, when combining reconstructed images of different heights into a whole, consideration must be given to eliminating the effects of twin-image and increasing the resolution of the reconstructed image. Secondly, the current process is not the same as the super-resolution algorithm. It is simply a pre-processing process for super-resolution to obtain the exact reconstructed distance. Therefore, this process should not take too much time. In this paper, we use the weighted average method, which can be expressed as: Where K is the number of hologram sequences with different heights and d k is the corresponding initial reconstruction distance based on sections 2.1, 2.2. In the actual operation of multi-height lensfree super-resolution, the defocus range of the reconstruction distance is very small, and the highfrequency detail information of the reconstructed image will change significantly. On the whole, by summing and averaging the reconstructed images of different heights, the noise such as twin-image can be quickly reduced and the high-frequency information of the reconstructed image can also be improved. In actual operation, there are two problems that need special attention. On the one hand, it should be emphasized that in the actual operation of lensfree on-chip holography, the smaller the distance between the sample and the sensor, the more highfrequency information of the object. Therefore, in the summation process, α is introduced as a weighting factor. In this article, α is set to 1/k for simplicity, and theoretical relationships are more complex. On the other hand, when the sample or sensor is displaced in the axial direction, lateral displacement will inevitably occur in the X-Y plane, which will cause serious spatial domain aliasing. In this article, the cross-correlation sub-pixel shift method [33] is applied to correct lateral position errors. The reconstructed image closest to the sensor is used as the initial position reference. The reconstructed image after correction U cor (x shift , y shift ) can be expressed as: (12) Where, (x shift , y shift ) is the sub-pixel shift calculated by the crosscorrelation method.
(f x , f y ) is the frequency domain coordinates.
Therefore, Equation (11) can be rewritten as: Us 1 =Ū(x, y, d 1 ), a = 1/k k = 1, 2, 3, ......, K − 1 So far, we have obtained a preliminary set of reconstructed distance sequences of holograms and synthetic reconstructed images with less noise and more detailed information. We will then evaluate the synthetic reconstructed images to further determine more accurate reconstructed distances. The reconstructed distance sequence [d 1 , d 2 , ......, d k ] will also be processed as a whole. Now we can further narrow the interval to [z 3 , z 4 ](z 1 < z 3 < z 4 < z 2 )with a smaller step size 1

SIMULATIONS
In this section, we will perform simulation experiments to verify the principles of the proposed algorithm. According to the proposed algorithm as shown in Figure 2, the whole simulation is divided into two parts: the individual preliminary estimation of the reconstructed distance sequence and the global precise estimation. The simulation procedure is described as follows. First, we simulate an object with amplitude contrast which has a size of 520×520 pixels. The object is propagated to 10 planes by an angular spectrum algorithm. The closest distance between the hologram and the sensor is 500 µm, and the distance interval of each hologram is set to 20 µm. In the actual lensfree onchip holography, the resolution is limited by the pixel size, so in the simulation experiment, the object light wave is also down-sampled [34] with a down-sampling factor of 4 as shown in Figure 3.

The Individually Preliminary Estimation
In this section, we need to initially calculate the reconstruction distances for holograms of different heights. To verify the accuracy of the proposed critical functions, we used the closest hologram to the sensor to calculate the reconstruction distance (d 1 = 500 µm). We calculated the critical functions for DIF and SDIF separately. The initial reconstruction distance range is [485,515] µm with the interval of 2 = 0.5 µm as shown in Figure 4A. By calculating the minimum value of critical functions, we can see that the reconstruction distance is about 500 µm for the critical functions of DIF and the SDIF. In order to determine the reconstruction distance more accurately, we shorten the interval and step size with the reconstruction distance range of [498, 503] µm, the interval of 2 = 0.5 µm. As can be seen from Figure 4B, critical functions of DIF changes smoothly in the range of [499,501] µm, and has no obvious minimum value. However, the critical functions of SDIF can also clearly determine the minimum value of 500.5 µm in the range of [500, 501] µm.

The Global Precise Estimation
For lensfree on-chip holography, the reconstruction distance accuracy directly affects the resolution of high-frequency information. Figure 4C shows the reconstructed intensity at different reconstruction distances. The system has no lens, and the reconstructed image is affected by the twin-image. At this point, it is difficult to directly determine the focus distance of the high frequency information based on a single low-resolution hologram. Based on section 3.1, we can get the reconstruction distance sequence [500.5, 520.5, 540.5, 560.5, 580.5, 600.5, 620.5, 640.5, 660.5, 680.5] µm corresponding to different height holograms. Then, the synthesized intensity can be obtained by Equations (12) and (13) as shown in Figure 5A. It can be clearly seen that the synthesized intensity has clearer detailed information compared with the low-resolution light intensity. Every time the entire distance sequence moves d, a new synthesized intensity image can be obtained. In this case, the proposed method transforms the problem of focusing on different planes into a problem of solving the optimal solution of the overall system. By calculating the minimum value of SDIF for synthesized intensity images, the overall displacement distance d can be obtained as shown in Figure 5B. In the global precise estimation, the range of d is

Experimental Setup
To identify the robustness and adaptability of the proposed method, we respectively conducted experiments on an amplitude-contrast object and a phase-contrast object, respectively, based on the lensfree on-chip holography system as shown in Figure 6. The system mainly contains four parts: a single mode fiber-coupled light source with the wavelength of 661.5 nm (LP660-SF20, Thorlabs, the United States), a monochrome imaging device with 3,872×2,764, and pixel size of 1.67 µm (DMM 27UJ003-ML, the imaging source, Germany), displacement platform for sample and the motor drive hardware. In the system, the distance between the light source and the sample is about 20 cm. The distance between the sensor and the sample is between 500 and 1,000 µm, and the distance can be adjusted by the motor drive hardware. M1, M2, M3 are fixed reflecting prisms and M4 is a reflecting prism that can be moved around the x and y axes. The sensor is fixed to a displacement table that can be moved precisely vertically by means of a motor drive. However, in actual experiments, when the z-axis distance is adjusted by the motor drive hardware, there is usually a deviation between the preset distance and the actual distance. In addition, during the movement of the z-axis, there will usually be horizontal displacement due to mechanical vibration. It should be emphasized that only the z-axis height needs to be adjusted in our experiments, and there is no need to change the x-axis and y-axis distances.In the following experiment, we will set the hardware program to raise the sample stage by 20 µm in turn and collect 10 low-resolution holograms as experimental samples.

Amplitude-Contrast Objects
In the first experiment, a standard 2 ′′ × 2 ′′ positive 1951 USAF resolution test target is regarded as amplitude-contrast objects to verify our proposed method. Through the hardware drive setting 10 times in the Z-axis direction, 10 raw holograms at different sample-to-sensor distances can be obtained. The first hologram captured by the camera as shown in Figure 7A. In order to verify that the algorithm automatically focuses on highfrequency information, we mainly analyze the information of the rectangular frame as shown in Figures 7B,C. We will compare the proposed SDIF with three typically existing critical functions, namely DIF, Sobel, SPEC as shown in Figure 7D. It is clearly seen that four methods almost locate the focal plane range of [585, 600] µm. In the focal plane, the critical functions of Sobel and SPEC search for the maximum value of the critical function, while functions DIF and SDIF search for the minimum value of the critical function. However, the critical functions of Sobel, SPEC, and DIF have multiple extremes near the focal plane, which has not a good sensitivity. In high frequency information, reconstruction distance errors of a few microns will directly affect the focusing accuracy. As can be seen in Figures 7B,C, when the reconstruction distance is in the range of [580,610] µm, the resolution is essentially the same for group 7, However, the reconstruction distance has a large impact on the resolution of groups 8 and 9. Thus, when the resolution is not high, the method described above allows a better focus plane to be determined; but when the resolution is high, especially when pixel super-resolution is performed, the focus distance must be determined accurately. We can see that the SDIF has one of the sharpest valleys, indicating a better sensitivity than the traditional method. Based on the critical function of SDIF, the sequence of   reconstructed distances for the different planes can be uniquely determined, as shown in Figure 7E.
To determine the reconstructed distances in each plane more accurately, the 10 reconstructed intensity images were transformed into one synthetic intensity image according to Equations (12) and (13), as shown in Figure 7F. It can be clearly seen that the synthesized intensity image has fewer twin-image and clearer detail information. Then, by finding the minimum value of the SDIF critical function, the optimal distance for the overall shift of the reconstructed distance sequence can be determined, as shown in Figure 7F. Therefore Figure 7G shows the high-resolution intensity image obtained by the APLI algorithm, which proves the effectiveness of the proposed algorithm.

Phase-Contrast Objects
In order to evaluate the capability of the proposed method to detect the pure phase object as shown in Figure 8A, we have taken the second experiment. The sample is a benchmark quantitative phase microscopy target, which is completely transparent and can be considered as the phase-contrast sample. The first hologram captured by the camera as shown in Figure 7A. In order to verify that the algorithm automatically focuses on high-frequency information, we mainly analyze the information of the rectangular frame as shown in Figures 8B,C.  Figure 8D shows the criterion curve of SDIF, DIF, Sobel, and SPEC. It is clearly seen that the critical function of DIF is the form of approximate cosine oscillation, which is difficult to determine the focus position by searching extreme values. The critical functions of Sobel, SPEC, and SDIF could locate the focal plane range of [425,435] µm by searching the minimum value of the critical functions. For phase-contrast objects, the focus position corresponds to the maximum point of the SDIF critical function. Based on the critical function of SDIF, the initial reconstruction distance sequences of different planes can be uniquely determined [427,448,469,486,506,528,550,568,586,606] µm as shown in Figure 8E. By calculating the maximum value of SDIF for synthesized phase images with a series of overall displacement distance as shown in Figure 8F, Figure 8G shows the high-resolution phase image obtained by the APLI algorithm. It can be clearly seen that the resolution of the phase image has been significantly improved, which proves the effectiveness of the proposed algorithm.

CONCLUSIONS
In this paper, we present an autofocusing method for multiheight pixel super-resolution lensless on-chip holography. In order to eliminate the influence of twin-image and improve the sensitivity of the critical function, the proposed method is divided into two steps: individual preliminary estimation and global accurate estimation. First, an improved critical function, named SDIF, is proposed to obtain approximate focal distances in different height planes by combining the differential critical function and frequency domain autofocus. To further obtain more accurate reconstructed distance sequences, the intensity or phase images of different heights are synthesized into a high-resolution reference image without twin-image. We transform the focusing problem on different planes into the problem of solving the optimal solution of the whole system, where we consider a series of holograms of different heights as a whole. By shifting the whole of the distance sequence, we can obtain a high-resolution reference image. With the proposed method, we can determine the focal distance at different plane heights more accurately and improve the supersolution accuracy of the imaging. Compared with conventional methods, the proposed method is more sensitive to the accurate reconstruction of distances for high frequency information and is less disturbed by twin-image. Simulations and practical experiments demonstrate the correctness and stability of the proposed algorithm for amplitude-contrast and phase-contrast objects. Considering that the proposed algorithm performs the estimation of individual reconstruction distances and the estimation of overall reconstruction distances, the time spent by the algorithm will increase. In the experiments, the time consumption for 10 images was around 50 s, and how to further reduce computational complexity is the next step of our research.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
YW conceived the idea. LL and ZL did the experiments. JZ provided Super-resolution algorithm support. CZ supervised the project and provided Technical Support. All the authors contributed to discussion on the results for this manuscript.