Fourier ptychographic microscopy (FPM)^[1–3] is a recently developed optical microscopy for wide-field, high-resolution, and quantitative phase imaging. FPM combines the concepts of ptychography^[4–7], synthetic aperture^[8–10], and phase retrieval^[11–14]. In a typical FPM setup, a programmable light emitting diode (LED) array serves as a partially coherent light source to provide angularly varying oblique plane wave illuminations. At each illumination angle, the image sensor records a low-resolution intensity image of the specimen. The captured low-resolution images contain different information in the Fourier domain and can be synthesized to recover a high-resolution image. The resolution of the recovered image is determined by the sum of the objective lens NA and the illumination NA^[15]. Similar to ptychography, FPM recovers the ‘lost’ phase information by applying phase retrieval methods. With no directly measured phase information or mechanical scanning, the construction of FPM is compact and low-cost. Within a few years, FPM has been applied in many fields, such as hematology^[16], pathology^[17,18], and quantitative phase imaging^[2]. These applications show the great potential of FPM in biomedical observation and clinical diagnosis.

FPM improves the space-bandwidth product (SBP) of a microscope system by collecting images with high-frequency information. As the images are orderly captured, the temporal resolution of the system will decrease^[19]. To reduce the acquisition time, some modifications have been proposed, such as high-brightness illuminations and the multiplexed method^[19–21]. As a result, it becomes possible to scan a living tissue with a proper speed. But the reconstruction process is still very slow compared to the capture process and usually consumes dozens of minutes. Low-reconstruction speed brings difficulties to dynamic observation in practical situations.

In this Letter, we propose a modified global reconstruction method for FPM termed gradational FPM (gFPM), which significantly improves the reconstruction speed. Our method applies a gradational strategy to realize iteration from low-frequency to high-frequency and use low-frequency reconstruction as complex initialization. In addition, a graphic processing unit (GPU) accelerated version of gFPM (GPU-gFPM) is implemented. We evaluate the noise performance and speed of the gFPM method, the GPU-gFPM method, the sequential method, and the global method. Both simulations and experiments demonstrate that a faster and more robust reconstruction is achievable utilizing the GPU-gFPM method.

In FPM, the bio-optical property of a thin sample can be represented by its transmission function o(r), where r=(x,y) represents the lateral coordinates at the sample plane. Assuming the specimen is illuminated by an oblique plane wave with wave vector ul(l=1,…,Nimg), the complex field on the sample plane is defined by exp(i2πulr). The exit field passing through the sample is expressed as e(r)=o(r)exp(i2πulr). By applying Fourier transform, the exit field is F{e(r)}=O(u−ul), where u=(fx,fy) represents the 2D frequency coordinates. After being low-pass filtered by the objective pupil function P(u), the estimated complex fields at the image plane can be written as gle(r)=F−1{P(u)O(u−ul)},(1)which is the forward model of FPM. The reconstruction can be performed by iteratively replacing the amplitudes of estimated complex fields with the captured intensities, which can be formulated as Ou(u−ul)=P(u)F{Ilc(r)|gle(r)|gle(r)},(2)where Ou(u−ul) is the updated sub-spectrum. By transforming the spectrum back to the spatial domain, the high-resolution intensity and phase images are extracted. To simplify the expression, the vectorization notations can be adopted. The FPM forward model can be formulated as gle=F−1(POl),(3)and the reconstruction can be expressed as Ou=∑lPF(Ilc|gle|gle).(4)

There are mainly two classes of strategies for FPM reconstruction: sequential and global. Sequential methods are widely used for their flexibility and efficiency, while the global methods perform better reconstruction at the cost of low reconstruction speed. The reconstruction process with the global method is shown in Fig. 1.

In gFPM, we introduce the concept of ‘gradation,’ which corresponds to different regions in the Fourier domain. The low gradation covers the low-frequency region in the Fourier domain, and upper gradations cover the higher frequency region. The captured intensity images are divided into several sets, corresponding to different gradations. Figure 2 shows a flow chart of gFPM, and the detailed process is as follows.

At first, the upsampled center low-resolution intensity and zero phase are used as the initial guess of the object field, labelled as O0. The pupil function is set as a circular-shape, zero-phase, and low-pass filter. The radius of the pass band is 2π·NA/λ, where NA is the numerical aperture of the objective lens, and λ is the wavelength.

Secondly, the ‘gradation-1’ iteration is carried out using O0 as the initialization. We based it on the global method in consideration of its superiority and potential. In gFPM, only 3×3 images in the center of the array are employed in ‘gradation-1’. After K iterations, a low-gradation object spectrum with an acceptable error acts as a good starting point for upper gradation iterations.

Thirdly, the ‘gradation-1’ result is used as the initial of the upper gradation iterations. Another set of images is applied to update the high-frequency region of the spectrum. This iteration process is repeated at each gradation. With a complex initialization, upper gradation iterations of gFPM converge much more quickly than the global method with the intensity-only initialization.

Finally, after the last gradation is carried out, the result will converge to a high space-bandwidth spectrum. The whole reconstruction process is expressed as {gl1=F−1(POl0),O1=W1∑lPF(Il1|gl1|gl1)gl2=F−1(POl1),O2=W2∑lPF(Il2|gl2|gl2)⋮⋮glj=F−1(POl(j−1)),Oj=Wj∑lPF(Ilj|glj|glj),(5)where Oj(1≤j≤J) denotes the recovered object spectrum in ‘gradation-j,’ and Olj denotes the shifted spectrum corresponding to the lth LED. glj is the estimated complex field corresponding to the lth LED in ‘gradation-j,’ and Ilj is the captured intensity image under the same condition. The subscript l has a different range in each gradation because different image sets are applied to each gradation.

The reconstruction results of different gradations of the gFPM method are shown in Fig. 3. The covered region enlarges, and the recovered images become clearer as the gradation grows. The gFPM converges much faster because of the gradational strategy and the calculation reduction of each iteration. As a modified global method, the gFPM method can be further accelerated by using a GPU.

We first validate gFPM with simulations. The parameters are carefully chosen to match a real system, with a 629 nm wavelength, 6.5 μm pixel size, and a 0.09 NA 4× objective. The 9×9 intensity images are collected to reconstruct the complex field. The overlap of the sub-spectrums is set as 0.6 because this value well balances the reconstruction quality and speed. We compare the reconstruction result of gFPM with the conventional sequential method and the global method. The criterion we adopt is the structure similarity (SSIM) for its better consistency with human perception than root-mean-squared error (RMSE), which is commonly used in previous studies. The SSIM is higher when two images are more similar in structure and ranges from zero to one. The computer we use is equipped with an Intel i7-6700K CPU, 16GB DDR4 memory and an NVidia GTX 1080 graphics card. The software we use is Matlab 2016b installed on a Windows 10 operating system.

The FPM reconstruction results are deteriorated by noise, because the noise badly decreases the quality of dark-field images. The tolerance of noise intuitively shows the robustness of the algorithm. As the reconstruction of gFPM is the same with GPU-gFPM, we use the sequential method, global method, and GPU-gFPM method to recover the images under different levels of Gaussian noise and compare the results.

Figure 4 shows the SSIM scores between the reconstruction results and the ground truth with the standard deviation of noise σ varying from 1×10−4 to 1×10−3. The intensity and phase reconstruction of the sequential method, global method, and GPU-gFPM when σ is 7×10−4 are also shown. The simulations show that the results of GPU-gFPM are better than the sequential method and global method in different levels of Gaussian noise, which indicates that the GPU-gFPM method is more robust when non-negligible noise exists.

The reconstruction speed is another performance we care most about. We compare the reconstruction speed of the sequential method, global method, and GPU-gFPM method under the same condition. Two metrics we used are the number of iterations and total runtime. The SSIM score is used as a convergence index to evaluate the reconstruction process. The iteration stops when the differences of three continuous SSIMs are less than 1×10−4. As for gFPM method, the iteration turns to a higher level when the differences between continuous SSIMs are less than 1×10−3 and stops when the differences of three continuous SSIMs are less than 1×10−4 in the last level.

Figure 5 shows the phase SSIM curves with the overlap set as 0.4, 0.5, 0.6, and 0.7 and the reconstruction results of the three methods with 0.6 overlap. It can be reached from Fig. 5 that the gFPM strategy significantly increases the reconstruction speed and reaches better results. Besides, the gFPM strategy can reduce the overlap requirements to a certain extent.

Table 1 shows that the runtime of gFPM is much less than the sequential method and the global method. Besides, the GPU computing contributes several times acceleration to gFPM. With the data size growing larger, the GPU acceleration can be even more efficient. In the reconstruction using 9×9 intensity images of 512×512 size, the GPU-gFPM method can be five times faster than the sequential method and 23 times faster than the global method.

We verify the effectiveness of the GPU-gFPM method on a real FPM platform, as shown in Fig. 6. The microscope is equipped with a four times objective (0.13 NA) and a scientific CMOS camera with 2560×2160pixels (6.5 μm pixel size). A customized 13×13 LED array is carefully calibrated for angularly varying illumination (629 nm central wavelength). The LED array is controlled using an Arduino circuit board. The distance between adjacent LEDs is 8 mm, and the distance between the sample and the LED array is 100 mm. All images are recorded with 100 ms exposure time, and a noise suppression operation is applied to all images.

Figure 7(a) presents the whole field of view (FOV) of a 1951 USAF resolution target, and Fig. 7(b1) shows a magnified region of interest (ROI). Figures 7(b2)–7(b4) are recovered intensity images of the region with the sequential method, global method, and GPU-gFPM method, respectively. Figures 7(c1)–7(c4) show the line traces of the resolution target image corresponding to Figs. 7(b1)–7(b4). All three methods are carried out with enough iterations for the convergence of results. The recovered images have higher contrast than the raw image, which indicates that the recovered image has a higher spatial resolution. Besides, the global method and GPU-gFPM method outperform the sequential method.

We also test GPU-gFPM on an animal testis tissue, and compare the GPU-gFPM with the sequential method and global method. Figure 8(a) presents the FOV of the specimen, and Fig. 8(b) shows a magnified ROI. Figures 8(c1)–8(c3) are recovered intensity images of the ROI with the sequential method, global method, and GPU-gFPM method, respectively. Figures 8(d1)–8(d3) show the recovered phase images of the ROI with the three methods. It can be observed that Fig. 8(c1) contains fewer details and more distortions comparing with Figs. 8(d2) and 8(d3). In summary, the GPU-gFPM method outperforms the sequential method and has similar accuracy with the global method.

To evaluate the reconstruction speed of the four methods under actual conditions, the FOV is divided into 42 sections of 400×400pixels. The runtime of the reconstruction of the four methods are recorded in Table 2. It shows that the GPU-gFPM method accomplished the reconstruction at a much higher speed.

In conclusion, we propose a modified global method for FPM reconstruction termed gFPM, which is efficient on time and noise robust. We also validate the efficiency of gFPM through simulations and experiments. Different from previous methods that utilize the entire image set in each iteration, gFPM divides the image set into several parts that correspond to different gradations of iterations. By applying the image set gradation by gradation, the high-frequency information is gradually introduced into the recovered spectrum. For this reason, we call this strategy the gradational strategy. Simulations show that the gFPM performs better under different noise levels and converges faster compared to conventional sequential methods and global methods. The gFPM method can also reduce the overlap requirements of a system. The high parallelism of gFPM makes it possible for us to further accelerate gFPM with the GPU. The gFPM strategy and GPU acceleration greatly release the potential of global methods, which offers an approach to reach dynamic observation through real-time FPM reconstruction.

Although gFPM performs better than conventional sequential methods and global methods, it cannot currently recover the pupil aberration. The combination of gFPM and the gradient decent algorithm may overcome this weakness. Besides, the CUDA programming rather than Matlab built-in functions may utilize the computing power of the graphics card better. It will be a subject of future work to further enhance the performance and efficiency of gFPM and achieves real-time FPM reconstruction.