Full-field Fourier ptychography (FFP): spatially varying pupil modeling and its application for rapid field-dependent aberration metrology

Digital aberration measurement and removal play a prominent role in computational imaging platforms aimed at achieving simple and compact optical arrangements. A recent important class of such platforms is Fourier ptychography, which is geared towards efficiently creating gigapixel images with high resolution and large field of view (FOV). In current FP implementations, pupil aberration is often recovered at each small segment of the entire FOV. This reconstruction strategy fails to consider the field-dependent nature of the optical pupil. Given the power series expansion of the wavefront aberration, the spatially varying pupil can be fully characterized by tens of coefficients over the entire FOV. With this observation, we report a Full-field Fourier Ptychography (FFP) scheme for rapid and robust aberration metrology. The meaning of 'full-field' in FFP is referred to the recovering of the 'full-field' coefficients that govern the field-dependent pupil over the entire FOV. The optimization degrees of freedom are at least two orders of magnitude lower than the previous implementations. We show that the image acquisition process of FFP can be completed in ~1s and the spatially varying aberration of the entire FOV can be recovered in ~35s using a CPU. The reported approach may facilitate the further development of Fourier ptychography. Since no moving part or calibration target is needed in this approach, it may find important applications in aberration metrology. The derivation of the full-field coefficients and its extension for Zernike modes also provide a general tool for analyzing spatially varying aberrations in computational imaging systems.


Introduction
Characterization of spatially varying wavefront aberration is of critical importance in ophthalmology, microscopy, astronomy, photography, and lithography. The knowledge of system aberration allows aberration correction either actively through adaptive optics or passively with post-acquisition aberration removal. Different aberration characterization methods have been reported in past years, including interferometry 1, 2 , wavefront sensing 3-5 , phase retrieval 6,7 , optimization-based spatially-varying pupils recovery 8 , point spread function measurement 9 , calibration masks [10][11][12] , statistical modeling of weakly scattering object 13 , differential phase contrast imaging 14 , and etc. In ptychography, spatially varying probe beam can be recovered using a orthogonal relaxation approach 15 . Among these different implementations, digital aberration measurement and removal play an especially prominent role in computational imaging platforms aimed at achieving simple and compact optical arrangements. A recent important class of such platforms is Fourier ptychography [16][17][18][19][20][21][22][23][24][25][26][27] , which is geared towards efficiently creating gigapixel images with high resolution over a large field-of-view (FOV). In a standard FP system, we illuminate the object with angle-varied plane waves and collect the corresponding images through a low numerical aperture (NA) objective lens. As such, each captured image represents the information of a circular disk in the Fourier domain and its offset from the origin is determined by the illumination angle. We can then apply the phase retrieval process to synthesize the aperture information in the Fourier domain. Unlike the conventional synthetic aperture imaging approach, FP requires no phase information to stitch the apertures in the Fourier domain. Instead, the phase information is recovered thanks to the partial aperture overlap between adjacent acquisitions 28,29 . At the end of the FP reconstruction process, the synthesized information in the Fourier domain generates a higher resolution complex object image that retains the original large FOV set by the low-NA objective. Fig. 1 We use ( , ) to denote the pupil coordinate at the Fourier plane and ( , ) to denote the field coordinate at the image plane. The pupil aberration ( , ) is varying at different field coordinate ( , ) . Here we use off-axis coma aberration as an example with ( , )= · , where 'a' represents the mount of coma aberration. As one expects, the amount of coma is 0 at the center and become larger at the edge of the FOV. The field-dependent coefficient ( , ) follows the form of and we aim to recover the 'full-field' constant 'b' via Fourier ptychography in this work. Such full-field constants provide a higher-level description of the spatially varying aberration across the entire FOV.
In current FP implementations, we need to divide the entire FOV into smaller segments to handle the spatially varying pupil aberrations. At each small segment, the pupil aberration is treated as spatially invariant and can be locally recovered in an optimization process 14,19,30,31 . Such an aberration recovery strategy, however, fails to consider the field-dependent nature of the optical pupil. Given the power series expansion of the wavefront aberration, the spatially varying pupil can be fully characterized using tens of power series coefficients over the entire FOV. As an example, we show the field-dependent coma aberration in Fig. 1, where we use ( , ) to represent the field location of the image plane and ( , ) to represent the wavevector location of the pupil (Fourier) plane. For the coma aberration shown in Fig. 1, the wavefront follows the form of in the pupil plane and the coefficient ( , ) represents the amount of coma aberration at the field location ( , ). To the best of our knowledge, most existing aberration recovery schemes aim to recover the pupil or the coefficient ( , ) for a given field location ( , ). As such, all different segments of the FOV are treated independently and the recovered coefficient is only related to the local segment. As we will discuss in the next section, the coefficient ( , ) in Fig. 1 follows the form of in the image plane, where b is a power series constant across the entire FOV. Therefore, it is better to recover such a 'full-field' constant 'b' that provides a higher-level description of the spatially varying aberration across the entire FOV.
In this work, we report a Full-field Fourier Ptychography (FFP) scheme for rapid and robust aberration metrology. The meaning of 'full-field' in FFP is referred to the recovering of the 'full-field' power series constants that govern the field-dependent pupil over the entire FOV. We show that the image acquisition process of FFP can be completed in ~1s and the spatially varying aberration of the entire FOV can be recovered in ~35s using a regular CPU. We also demonstrate the reported approach for generating a gigapixel image of a biological sample with a ~8-mm FOV and a 0.5 synthetic NA.
Our approach has several advantages over the existing techniques for aberration metrology. First, we recover the 'full-field' constants instead of the pupil wavefront at each small segment. As such, the optimization degrees of freedom are orders of magnitude lower than the previous implementations. The field-dependent power series expansion can also be viewed as a strong constraint for the pupil recovery process, improving the robustness and accelerating the convergence speed. Second, the full-field constants are jointly recovered using segments over the entire FOV. As such, our reconstruction can be viewed as an averaging result across the entire FOV. It is, in general, more accurate than the localized recovery in a regular implementation. Third, we require no mechanical scanning in our approach. The only hardware modification is to place a low-cost LED array under the object. Such a simple and low-cost configuration enables aberration metrology in systems where traditional techniques are difficult to apply. Fourth, no pre-known calibration mask is needed in our approach. In the reconstruction process, we jointly recover the pupil aberration and the unknown object simultaneously. We use a low-cost blood smear slide in our demonstrations and any thin section would work in our approach. This paper is structured as follows. In Section 2, we will discuss the spatially varying pupil modeling and derive the field-dependent expression for both power series polynomials and Zernike modes. In Section 3, we will discuss the forward imaging model and the reconstruction process of the FFP approach. In Section 4, we will present our experimental measurements using the FFP approach and quantify the measurement errors. Finally, we will summarize the results in Section 5. The reported approach may greatly facilitate the development of Fourier ptychographic imaging platforms. Since no moving part or calibration target is needed in this approach, it may also find important applications in aberration metrology. The derivation of the full-field power series coefficients and its extension for Zernike modes also provide a general tool for analyzing spatially varying aberrations in computational imaging systems.

Spatially varying pupil aberration modeling and full-field coefficients
It is well-known that the pupil aberration depends on the field location 32 . In this section, we will derive the field-dependent expression via power series expansions and then extend it for Zernike modes. As shown in Fig. 2, we consider the light wave propagates from the pupil plane with coordinates ( , ) to the image plane with coordinates ( , ). If we fix the field location to be ( , ), the wavefront aberration will be a function of pupil coordinates only, i.e., , , and it represents the regular treatment of spatially-invariant aberration. To obtain the complete information about the system aberrations, we need to consider wavefront aberration at different field locations and it becomes a function of 4 variables ( , , , ). These 4 variables are needed for optical systems without rotational symmetry. With rotational symmetry assumption, we can simplify it to ( + , + , + ), where three variables are enough to describe the system aberrations 32 . Fig. 2 Spatially varying pupil aberration modeling. The wavefront aberration is expanded to a power series, with field-dependence on each term.
In Fig. 2, we expand this aberration expression as a power series in three variables. Based on the power series expansion, we summarize the aberration expressions and its field-dependent function in Fig. 3. The first term, + , represents the defocus aberration and its field dependence is a constant plus a quadratic and a biquadratic term. The second term in Fig. 3, , represents the astigmatism aberration and its field dependent function has a quadratic and a biquadratic term. Similarly, represents coma aberration and its field dependence is a cubic function. We note that, the field-dependent function in Fig. 3 also provide some physical insights on the aberrations. For example, coma aberration is often known as an off-axis aberration and it will vanish at the on-axis position. This point can be verified by the field-dependent function which equals to 0 at the on-axis position. In our FFP recovery scheme, we aim to recover the power series coefficients, , , , etc. We term these coefficients as 'full-field coefficients', which govern the spatially varying aberrations over the entire FOV. Compare to the recovery of the aberrations of a local segment, the full-field coefficients represent a higher-level description of system aberrations. Tens of such coefficients can fully characterize the aberration over the entire FOV. As such, the optimization degrees of freedom are at least two orders of magnitude lower than the previous approaches, enabling a rapid and robust aberration metrology scheme. Fig. 3 Pupil-plane wavefront aberrations and their field-dependence. We aim to recover the power series coefficients ( , , , …) in our FFP scheme. These coefficients are termed 'full-field coefficients', which govern the spatially varying aberrations over the entire FOV. Zernike aberrations are combinations of different Taylor polynomials terms shown in Fig. 3. We also provide the field-dependent functions of Zernike aberrations in Fig. 4, where we ignore the biquadratic terms in the field-dependent function. We also note that, the coefficients in Fig. 4, , , …, are not independent with each other. For example, = and = for a rotational symmetry system. The derivation of their relationship is beyond the scope of this work.

Forward imaging model and the optimization problem
In the proposed FFP approach, the forward imaging model is similar to the original FP paper 16 . We use J different angle-varied plane waves to illuminate complex object ( , ) and acquire the corresponding intensity images ( , ) (j = 1,2,3…J) via a low-NA system. In our implementation, we divide the captured image ( , ) into M different small segments and obtain ( , ), where m = 1,2,3…M. For each small segment, the wavefront aberration is treated as a spatially invariant pupil and we can obtain the following forward imaging model for the m th image segment: where '*' denotes convolution, ( , ) denotes the m th small segment of the object, and ( , ) denotes the spatially invariant point spread function (PSF) at the m th small segment.
The convolution operation in Eq. (1) is typically calculated in the Fourier domain via fast Fourier transform (FFT). We can define the following complex exit wave , at the pupil plane using the pupil function , as follows: where , is the Fourier transform of ( , ) and , is the Fourier transform of ( , ). Based on the definition in Eq. (2), the forward imaging model of Eq. (1) can be rewritten as where ℱ denotes inverse Fourier transform. Assuming the central position of the m th image segment is ( , ), the pupil function for the m th image segment in Eq. (2) can be expressed as follows where ' ' denotes a circular mask with a radius of ⋅ and is the wavelength of the light wave.
In Eq. The goal of the FFP approach is to recover the 13 full-field coefficients { , , … , } in 17 fielddependent functions in Fig. 3 (several full-field functions share the same full-field coefficients). We define the following cost function for the m th image segment and the j th incident angle: The full-field coefficients can be recovered by solving the following optimization problem: Compared with other pupil recovery schemes, the key innovation here is to impose the constraints of the full-field functions. As such, we only need to recover 13 parameters over the entire FOV. The optimization degrees of freedom are even lower than the spatially invariant blind deconvolution of a single image segment. In our implementation, we divide the entire FOV into ~320 smaller segments and each segment has 256 by 256 pixels.

Recovering full-field coefficients via gradient descent
We provide two solutions for solving the optimization problem in Eq. (6). The first solution is an extension of the EPRY (or ePIE) approach 19,33 . In the recovery process, we update the object ( , ) via the EPRY approach and then update the full-field coefficients { , , … , } via gradient descent. Different from previous single-segment implementations, the gradient update for the full-field coefficients is based on the average from M segments over the entire FOV. Fig. 5 Recovering full-field coefficients via gradient descent. Figure 5 summarizes the recovery process, where the subscript 'm' denotes different image segments over the entire FOV, the subscript 'j' denotes the j th incident angle, and the subscript 't' denotes the t th aberration modes in Eq. (4). In this process, we first initialize the object spectrum , and the 13 full-field coefficients { , , … , }. For each iteration n, the images ( , ) are addressed according to different segments and different incident angles. We perform the Fourier magnitude projection as follows: The updated complex object is then used to generate the updated spectrum of the exit wave at the pupil plane: , = ℱ ( , ) . We then update the object spectrum using the following equation: Instead of updating the pupil function as in EPRY, we directly update the full-field coefficients { , , … , }. The gradient of the loss function with respect to , are given by: where ( , ) and ( , ) represent the gradient of ( , ) with respective to and , ' ' represents complex conjugate, and ' ' means taking the imaginary part. The gradient of the loss function with respect to other coefficients can be obtained in a similar manner. For each full-field coefficient, we calculate the update (Δ , Δ , …) by taking the gradient average of all M small segments: The gradient descent method for updating the full-field coefficients is, thus, given by: where α is the step-size for the updating process. This process is repeated for n iterations until the solution converges. In a typical implementation, we use 10-50 loops. With the full-field coefficients, the pupil function at any location of the entire FOV can be obtained via Eq. (4).

Fast implementation via alternating projection
The second solution for solving the optimization problem in Eq. (6) is based on alternating projection. In this case, we first rewrite the pupil function as where represents the aberration coefficient of the t th aberration mode on the m th image segment. In this scheme, we update the coefficients and then fit these coefficients to the full-field functions ( , ).
Based on the gradient in Eq. (14), we then update the aberration coefficient via With the updated , we fit it to the full-field functions ( , ) to update the full-field coefficients. For example, given = 1, ( = 1,2, … ) are fitted to ( , ) = + ( + ) + ( + ) . The fitting process can be considered as the following optimization problem: Based on Eq. (16), we can get the full-field functions with the updated full-field coefficients , , . We then update the aberration coefficient again by projecting it to the full-field functions, i.e., = ( , ). This process is repeated for n iterations until the solution converges. The full-field functions can be viewed as constraints for the optimization problem. The updating process of can be viewed as a projection operation. This scheme is ideal for parallel processing since different image segments can be processed at the same time. In a typical implementation, we use 10-100 loops for reaching a converged solution.

Field sampling and other considerations
There are several considerations for the reported FFP approach. First, we do not need to use all image segments in our implementation. Where to choose these image segments is an important consideration. In Fig. 7, we discuss different sampling strategy 34 for choosing the image segments: 1) uniform sampling over the entire FOV, 2) a higher sampling density at the center, and 3) a higher sampling density at the edge of the FOV. The ground truth aberration is generated using 320 image segments with 256 by 256 pixels each and with 25 incident angles. To test the three cases listed above, we use 81 segments with 64 by 64 pixels each and with 15 incident angles. We quantify the results via root mean squared error (RMSE). We can see that a higher sampling density at the center give a higher convergence speed. Other considerations for the reported approach include the number of image segments, the number of pixels for each segment, and the number of incident angle for each segment. In our implementation, we use 81 image segments, 64 by 64 pixels per segment, and 15 incident angles. This choice is a good compromise between the convergence speed and processing time per iteration. In particular, the acquisition time of the FFP approach is less than 1s (we acquire ~20 brightfield images in total). The spatially varying aberration of the entire FOV can be recovered in ~35s using the alternating projection scheme with parallel processing of an Intel i7 CPU.

Field-dependent aberration metrology and full-field FP imaging
In this section, we demonstrate the use of the FFP approach for aberration metrology and full-field FP imaging. In our imaging setup, we illuminate the object using an LED array and acquire images using a Nikon microscope with a 2X, 0.1-NA objective lens. The object is a low-cost blood smear ($6, Carolina, human blood film slide), which has rich features over the entire FOV.
In the first experiment, we use power series (Taylor) polynomials in Fig. 3 as the aberration modes. Figure 8 shows the recovered full-field functions for 7 aberration modes (out of 17 modes in total). Similarly, we can also use Zernike polynomials in Fig. 4 as the aberration modes. The recovered full-field functions for Zernike modes are shown in Fig. 9. In this case, we have not considered the relationship between different full-field coefficients in Fig. 4. Therefore, we recover all 26 parameters listed in Fig. 4 and assume they are independent with each other. Fig. 8 Spatially varying aberrations of the microscopy imaging system with a 2X 0.1 NA objective lens. (a)-(g) The recovered field-dependent functions for 7 Taylor-polynomial aberration modes, where ( , ) represents the field coordinates in millimeter and the z axis represents the amount of the aberration mode. Fig. 9 The recovered field-dependent functions for 7 Zernike aberration modes, where ( , ) represents the field coordinates in millimeter and the z axis represents the amount of the aberration mode.
In a regular FP implementation, one need to recover both the object and the pupil aberration for each segment. Figure 10 shows the comparison between the regular FP implementations and the reported FFP implementations. Figure 10(a) shows the raw image segment of the blood smear. We choose a region close to the edge of the FOV in this demonstration. Figure 10(b1) and 10(b2) show the recovered intensity and phase using the EPRY-FPM approach 19,33 . Figure 10(c) shows the recovered results using the newly reported 'rPIE+Momentum' approach 35 . For both approaches, the recoveries cannot converge to a good solution due to the large pupil aberration presented in the system. To address this problem, a conventional FP implementation needs to recover the segments from the center to the edge sequentially. The recovered pupil from the central segment will be used as the initial guess for the adjacent segments. In the reported FFP, we can first recover the full-field functions as shown in Figs. [8][9]. We can then obtain the pupil aberration for each segment. Finally, we can recover the high-resolution object without updating the pupil function. Different image segments can be processed in parallel since sequential pupil updating is not needed. We also note that, both the gradient-descent and the alternating projection schemes can reach the converged solution in Fig. 10(d) and 10(e). The alternating projection scheme is, in general, preferred because it allows efficient parallel processing of different image segments at the same time. Based on the pupil recovered from the FFP approach, we also generate a full FOV, high-resolution image of the blood smear sample in Fig. 11. In this experiment, the entire area is divided into ~320 segments. We first recover the field-dependent pupil function for each segment via the FFP approach. We then recover the object without updating the pupil function in the reconstruction process. Figure 11(a) shows the recovered full-FOV gigapixel image of the blood smear, with a synthetic NA of 0.5. The insets in Fig. 11(a) show the wavefront aberration of three regions. The aberration at the edge is much more severe compared to that at the center. Figure 11(b)-(d) show the magnified views of recovered intensity, recovered phase, and the raw images in three locations, where we can see significant resolution improvement and clear phase profiles of the blood cells. Finally, we quantify the accuracy of the recovered pupil aberration of the reported FFP approach. In this experiment, we introduce known defocus aberrations by setting the sample to different defocus positions: z = 0 m, 1 m, 2 m, 5 m, 10 m, 15 m, and 20 m. For each position, we acquire an FFP dataset and recover the full-field functions. Figure 12(a) shows two recovered full-field functions for the defocus aberration mode (i.e., the + term in the pupil plane). The bottom and top surfaces in Fig.   12(a) correspond to the cases by placing the object at the 0-m and 20-m z-positions, respectively. As such, we have the recovered pupils at both the 0-m and 20-m positions. To quantify the accuracy of the recovered pupil aberrations, we first generate a standard pupil wavefront with 20-m defocus distance. We then add this generated 20-m defocus pupil to the recovered 0-m pupil. The result is served as the ground-truth for 20-m defocus pupil. The difference between the ground-truth and the recovered 20-m pupil is the wavefront error in the reconstruction process. This wavefront error can be quantified at the different regions of the entire FOV, as shown in Fig. 12(b1)-12(b4). We can see that the wavefront errors are on the nanometer scale and there is no significant difference between different regions. We further quantify the root-mean-square (RMS) wavefront errors with different defocus distances in Fig. 12(c), where the red dots represent the average wavefront errors in all regions and the error bars represent the standard deviation at different regions of the entire FOV. Again, the RMS wavefront errors are on the nanometer scale, validating the accuracy of the recovered pupils.

Summary
In summary, we report a full-field Fourier ptychography (FFP) scheme for rapid and robust aberration metrology. The power series aberration modes and their field-dependent functions are given for modelling the spatially varying aberration in our optimization scheme. Compared to the regular aberration recovery approaches, we recover the 'full-field' coefficients instead of the pupil wavefront at each small segment. As such, the optimization degrees of freedom are orders of magnitude lower than the previous implementations. We show that the image acquisition process of FFP can be completed in ~1s and the spatially varying aberration of the entire FOV can be recovered in ~35s using a regular CPU. We also demonstrate the reported approach for generating a gigapixel image of a biological sample with a ~8-mm FOV and a 0.5 synthetic numerical aperture. The full-field functions employed in our scheme can be viewed as a strong constraint for the pupil recovery process. It can accelerate the optimization process and improve the robustness of the reconstruction. Future directions of the reported approach include the testing of other advanced optimization schemes for reducing the number of acquired images. Our on-going effort is to implement the full-field model for conventional blind deconvolution problems.