Fourier ptychographic reconstruction using Wirtinger flow optimization

Recently Fourier Ptychography (FP) has attracted great attention, due to its marked effectiveness in leveraging snapshot numbers for spatial resolution in large field-of-view imaging. To acquire high signal-to-noise-ratio (SNR) images under angularly varying illuminations for subsequent reconstruction, FP requires long exposure time, which largely limits its practical applications. In this paper, based on the recently reported Wirtinger flow algorithm, we propose an iterative optimization framework incorporating phase retrieval and noise relaxation together, to realize FP reconstruction using low SNR images captured under short exposure time. Experiments on both synthetic and real captured data validate the effectiveness of the proposed reconstruction method. Specifically, the proposed technique could save around 80% exposure time to achieve similar retrieval accuracy compared to the conventional FP. Besides, we have released our source code for non-commercial use.


Introduction
Fourier Ptychography (FP) is a newly reported technique for large field-of-view (FOV) and high-resolution (HR) imaging [1][2][3]. This technique sequentially captures a set of lowresolution (LR) images describing different spatial spectrum bands of the sample, and then stitches these spectrum bands together in the Fourier domain to reconstruct the entire HR spatial spectrum, including both amplitudes and phases. This HR spectrum can be transformed to the spatial domain to recover the HR image of the sample. Mathematically, FP reconstruction could be treated as a typical phase retrieval problem, which needs to recover a plural function given the magnitude measurements of its Fourier transform. Specifically, we only obtain the magnitudes of images corresponding to the sub-bands of the HR spatial spectrum, and intend to retrieve the plural HR spectrum. So far, all the FP applications [1,[4][5][6] utilize the alternating projection (AP) algorithm [7,8], a widely used method for phase retrieval, to implement the reconstruction process.
Recently, FP technique has been successfully applied to microscopic imaging and brings Fourier ptychography microscopy (FPM) [1]. FPM assumes plane-wave illuminations from the LED array. Therefore, by sequentially lightening LEDs at different positions in the illumination plane, according to [9], we can obtain different shifted versions of the sample's spatial spectrum. Since the whole light field is filtered by the microscope's objective lens, the changing of incident angles results in a set of LR images carrying different sub-bands of the entire spectrum. Finally, by utilizing the FP reconstruction technique (the AP algorithm), FPM can achieve large FOV and HR microscopic imaging. As stated in [1], the synthetic NA of the FPM setup is ∼0.5, and the FOV could reach ∼120 mm 2 , which greatly improves the throughput performance of existing microscopes. However, the AP algorithm utilized in the reconstruction process is sensitive to the input noise [10], and thus long exposure time is required for capturing high signal-to-noise-ratio (SNR) inputs. As reported in [1], the FPM setup needs around 3 minutes to acquire high SNR images under 137 angularly varying illuminations, for subsequent reconstruction of a gigapixel grayscale image. This largely limits the practical applications of FPM and other FP applications. To shorten the acquisition time, this paper proposes a new phase retrieval method for FP reconstruction, which is able to deal with low SNR input images.
As for existing phase retrieval methods, typically they can be classified into two categories, namely alternating projection algorithms and semi-definite programming (SDP) based algorithms [10]. The former kind alternately operates in spatial space and Fourier space, imposing corresponding spatial-plane and Fourier-plane constraints to the retrieved plural function. These constraints include consistency with measured magnitudes [11], magnitudes' nonnegativity [12], signal support constraint [7], and so on. This kind of methods are efficient but at the risk of non-convergence and reaching to a local optimum [10]. The latter kind of approaches rely on the observation that the quadratic equations in the phase retrieval problem can be rewritten as linear equations in a higher dimensional space [13]. Typical SDP algorithms include PhaseLift [14,15] and PhaseCut [16], and PhaseLift has been successfully applied to the common phase retrieval task from captured coded diffraction patterns [17]. Such kind of methods could converge to a global optimum by a series of convex relaxations. However, they require matrix lifting to work in a higher space, and thus causes high computation cost, which makes them less competitive.
Recently, based on Wirtinger derivatives [18,19], Candes et al. [20] develop a non-convex formulation of the phase retrieval problem, and utilize the gradient descent scheme to derive a computation-cost-saving solution, termed Wirtinger flow (WF) algorithm. As stated in [20], although the quadratic model is non-convex, Wirtinger flow algorithm can rigorously retrieve exact phase information from a nearly minimal number of random measurements, by starting with a relatively accurate initialization. Besides, the algorithm converges at a geometric rate in the case of random Gaussian sampling mode, and at a linear rate in the case of coded diffraction pattern mode. In this paper, we apply the Wirtinger flow scheme to FP, and further introduce a noise relaxation constraint for a new FP reconstruction framework. The proposed framework is termed WFP (Wirtinger flow optimization for Fourier Ptychographic). The advantages of the proposed WFP are threefold: • Compared to the existing AP algorithm for FP, WFP can better handle the detector noise, and thus largely reduce the requisite long exposure time; • Compared to the SDP based algorithms such as PhaseLift and PhaseCut, WFP doesn't need matrix lifting and largely decreases the computation cost; • WFP is a general optimization framework being able to incorporate various priors and constraints, and hence can be extended to save costs and increase the retrieval accuracy further.
The remainder of this paper is organized as follows: modeling and derivation of the optimization algorithm are explained in Sec. 2. Then, we conduct a series of experiments on both synthetic and real captured data to validate the proposed approach in Sec. 3. Finally, we conclude this paper with some summaries and discussions in Sec. 4.

Optimization framework
In this section, we first review the Wirtinger flow algorithm [20], and then introduce the Wirtinger flow formulation into the FP reconstruction process, and incorporate the noise constraint. Finally, we derive the WFP reconstruction algorithm robust to the capturing noise.

Review of the Wirtinger flow algorithm
Wirtinger flow algorithm [20] is a recently reported technique to solve the standard phase retrieval problem. Specifically, it retrieves a plural signal x ∈ C n from a series of its real sampling measurements b ∈ R m , with the measurement formation defined as b = |Ax| 2 = (Ax) * Ax.
Here A ∈ C m×n is a linear sampling matrix, and stands for the dot product.
Based on the quadratic loss function, the Wirtinger flow algorithm transforms the phase retrieval task into a minimization problem as Here || · || F is the Frobenius norm, and is calculated as ||X|| F = ∑ i, j X 2 i j . Such an optimization model can be solved in an iterative manner, utilizing the gradient descent scheme. According to [18,19], the derivative of the complex quadratic cost function with respect to x * , i.e. ∂ f ∂ x * , is necessary for updating x in each iteration. In implementation, x is updated in a gradient descending manner as [20] x Here ∆ is the gradient descent step size set by users, and ∂ f ∂ x * can be easily calculated according to the Wirtinger derivatives as With the above derivations, the Wirtinger flow algorithm is summarized as Alg. 1.

Output
: Retrieved plural signal x.

Wirtinger flow optimization for Fourier Ptychographic-WFP
In terms of the FP reconstruction, the target is to recover the HR spatial spectrum from a series of LR images captured in spatial space. The relation between the HR reconstruction and the LR observations corresponds to two sequential linear operations: (i) down-sampling caused by the object aperture, and (ii) inverse Fourier transform to the LR spectrum bands caused by the microscope imaging system. We treat these two operations as a whole and use A ∈ C m×n to represent this combinational sampling process. Specifically, denoting the HR spatial spectrum as a plural vector x ∈ C n , the corresponding sampling matrix A = FS is composed of two components: inverse Fourier transform F and down-sampling S.
In addition, we also consider the capturing noise explicitly, and the measurement formation model becomes where n ∈ R m denotes the capturing noise, we assume which to be Gaussian. We use σ ∈ R to represent the standard deviation of the noise. The three sigma rule tells that, nearly all (99.73%) the samples of a random variable lie within 3 times standard deviation from its mean. Therefore, we can approximatively formulate our noise constraint as |n| 3σ .
Introducing a relaxation vector ε ε ε ∈ R m , we can transform the above inequality into an equality Combining the above noise constraint (Eq. 6) with the measurement formation (Eq. 4), we can get the optimization model for FP reconstruction as In the following, according to the Wirtinger flow algorithm, we derive the WFP optimization algorithm to solve the above model. First, we introduce a weighting parameter µ to incorporate the noise constraint into the objective function, and the model becomes This is similar to the augmented Lagrangian function [21], which can be solved by sequentially updating each variable [22], while keeping the other variables constant. The updating can be conducted either by assigning zero to the function's partial derivative with respect to the updating variable, or by the gradient descent technique. Here we utilize a similar scheme, and sequentially update the optimization variables, i.e., x, n and ε ε ε, in Eq. 8. For x, by calculating the partial derivative of f to x * , we can get its updating rule using the gradient descent technique as with ∆ x being the gradient descent step size of x. Similarly, we can set the step size of n as ∆ n , and update n following = n (k) − ∆ n (|Az|. 2 + n − b) + µ(n n − 9σ 2 + ε ε ε ε ε ε) 2n | n=n (k) .
Based on the above derivations, the proposed WFP algorithm is summarized in Alg. 2. As for the initialization, we set x (0) as the spatial spectrum of the up-sampled version of the LR image, which is captured under the normal incident light.

end
As for the parameters settings of ∆ x and ∆ n , similar to WF, we assign ∆ (k) ||σ σ σ || 2 , where θ (k) = min 1 − e −k/k 0 , θ max . As stated in [20], k 0 = 330 and θ max = 0.4 work well, so we also use these settings in our algorithm. We have released our source code for non-commercial use, which can be downloaded here.

Experiments
In this section, we conduct a series of experiments on both synthetic and real captured data to validate the proposed WFP algorithm.

Algorithms for comparison:
To demonstrate the performance and advantages of the proposed WFP algorithm, here we run the WFP as well as the AP algorithm on simulated FP data. Besides, to investigate WFP's ability in addressing acquisition noise, we further compare its results with those produced by applying denoising before or after AP reconstruction. Here we use BM3D [23] for denoising considering its promising results and high efficiency, as stated in [24]. For simplicity, we respectively refer to these two methods including the denoising operation as "BM3D+AP" and "AP+BM3D". Criterion: Besides the visual results, we also utilize two quantitative criteria to assess the recovery performance of the above methods. The first one is the peak signal to noise ratio (PSNR), which has traditionally been widely used to assess the quality of processed images compared to benchmark. PSNR intuitively describes the intensity difference between two images, and would be smaller for lower quality recovery. Another adopted criterion is the structure similarity (SSIM) [25] that measures the spatial structural closeness between two images, and thus consists with human perception better than PSNR. The SSIM score ranges from 0 to 1, and is higher when two images are of more similar structural information. Note that here both PSNR and SSIM are calculated on the intensity images.

Experiment parameter settings:
The convergence experiment in [20] shows that the Wirtinger flow algorithm works successfully when measurements are more than 6 times of the signal entries to be recovered. In terms of the FP problem, assuming that the overlap ratio is ξ , the LR images are of m × m pixels and we capture k × k LR images, the sampling ratio between measurements and signal entries can be calculated as Similar to [1], by setting ξ = 65%, we can get the sampling ratio as around 8. The ratio is higher than the minimum convergence requisition (∼6) in [20]. Therefore, we adopt the above experiment settings in our simulation experiment, namely ξ = 65%, k = 15, m = 100.

Results:
Based on the above specifications, the captured image volume in our simulation experiment is synthesized by following three steps: 1) we apply FFT to the original HR image, and select subregions corresponding to different incident angles, by multiplying the HR spectrum with an ideal pupil function (all ones in the pupil circle and zeros outside). 2) We shift these sub spatial spectra to the origin location, and do inverse Fourier transform to recover the LR plural images in the spatial domain. 3) We retain only the intensity of these LR plural images, and add Gaussian white noise to obtain the simulated captured noisy images. In our implementation, we use the 'Lena' and the 'Map' image (512 × 512 pixels) from the USC-SIPI image database [26] as the HR intensity and phase image, respectively. The LR images' pixel numbers are set to be one tenth of the HR image along both directions.
First, we apply the proposed WFP on the simulated data with varying noise levels to study the algorithm's performance. Specifically, the standard derivation σ of the additive noise ranges from 0.002 to 0.008 with a 0.002 interval. By algorithm testing, 500 iterations are enough for WFP to converge, and hence we set the iteration number of WFP as 500. The visual and quantitative results are respectively shown in Fig. 1 and Tab. 1. From the results we can see that WPF works well to reconstruct both intensity and phase information. Besides, as the noise level increases, the reconstruction quality does not degenerate much. This illustrates the robustness of WFP to different noise levels, and thus a wider applicability.
Then, we compare WFP with the above mentioned three other methods, i.e., AP, BM3D+AP, and AP+BM3D, to show their pros and cons. Here the noise level is fixed at σ = 0.004. The iteration number of conventional AP is set to 50 to ensure convergence. The simulated acquisition image under the normal illumination is shown in Fig. 2(a). The quantitative and visual reconstruction results are respectively shown in Tab. 2 and Fig. 2(b)-(d).
Due to the noise corruption, the SNR of captured images is very low, especially for the images corresponding to high spatial frequencies. As a result, the reconstruction intensity and phase image of AP in Fig. 2(b) are very noisy. When BM3D is applied before the AP reconstruction, a lot of high frequency information are filtered out, thus there exist serious artifacts in the final recovery, as shown in Fig. 2(c). If we apply BM3D after the AP reconstruction, though most of the noise is removed, many crucial image details are filtered out as well (see the areas of hat tassels in Fig .2(d)). Differently, the proposed WFP incorporates noise suppression into the reconstruction framework, and conducts these two operations jointly. This largely avoids error accumulation in successive processing and achieves higher performance. Consequently, WFP could obtain satisfying reconstruction results with less noise and more details. In a nutshell, WFP largely outperforms the other three methods, on both visual and quantitative metrics.  The advantageous performance is at the expense of high computational cost. We implement all the four different methods using Matlab on an Intel i7 3.6 GHz CPU computer, with 16G RAM and 64 bit Windows 7 system. The running time of different methods are listed in the bottom row of Tab. 2. Obviously, WFP takes longer time than the other methods.

Experiment on real captured data
To further validate the effectiveness of WFP, we build a FPM setup to capture LR images as inputs for FP reconstruction. Similar to [27], an upright microscope with a 2× (NA = 0.1) objective is used in the platform. The LED array is placed around 8cm under the specimen, and the lateral distance between two adjacent LEDs is 4 mm. The central wavelength of incident light is 632nm. The pixel size of the captured raw images is ∼1.85µm.
First, we capture the LR images corresponding to the 15×15 LED positions, with the exposure time for each LED as 1ms, and apply AP as well as WFP to the image set for performance comprison. Fig. 3(a) shows the LR images of the USAF chart and the blood smear sample captured under the normal illumination. The HR reconstruction results of AP and WFP are respectively presented in Fig. 3(b) and Fig. 3(c), where the left columns show recovered amplitudes, and the right columns present recovered phases. From the results we can conclude three advantages of WFP over conventional AP. First, the reconstruction results of WFP own higher resolution than those of AP, and thus contain more details (see the close-ups for clearer comparison). Second, WFP could effectively suppress noise (see the smooth regions in the USAF chart). Third, WFP could reconstruct much more accurate phase than AP. Note that there exists some phase jumps recovered by WFP in the feature areas of the USAF chart, where the magnitudes are close to zero. This is due to that in these areas, the phase can be any value and their assignments would not affect successful magnitude recovery.
Then, to quantitatively evaluate the advantage of WFP over conventional AP in the view of exposure time, we increase the exposure time, and apply AP to corresponding acquisition data under longer exposure time. See the amplitude reconstruction results in Fig. 4. The proposed WFP could resolve a comparable resolution under 1ms exposure time, to that of the conventional AP under 5ms exposure time. This indicates that WFP could save round 80% exposure time than AP to achieve the same reconstruction accuracy.
In all, WFP offers a feasible way for the FP technique to reconstruct highly accurate results when there exists non-ignorable capturing noise, such as the case of low exposure time, or when the hardware is not so precise. This will do lot of help in real applications.

Conclusions and discussions
This paper proposes a reconstruction framework termed Wirtinger flow optimization for Fourier Ptychography (WFP). Based on the recently reported Wirtinger flow algorithm, WFP formulates the FP recovery as a quadratic optimization problem, and presents a solution utilizing the gradient descent scheme. By incorporating priors on the capturing noise, WFP can save around 80% of the exposure time for the present FP technique, while without obvious performance degeneration. Results on both synthetic and real captured data validate the effectiveness of WFP.
One extension of WFP is to handle non-uniform noise. This can be easily realized by treating the standard derivation of noise σ in Eq. (5) as spatial non-uniform, namely by changing σ ∈ R to σ σ σ ∈ R m . Besides, as a flexible optimization framework, WFP can also be easily extended by introducing other priors and constraints. For example, we can incorporate the sparsity of the latent HR image [28] into our framework, which may further reduce the snapshot number and thus the acquisition time. We can also introduce the total variance prior [29,30] into WFP to further suppress noise in the reconstruction results. What's more, in the optimization model in Eqs. (7), the sampling matrix can be composed of any kinds of linear operations (downsampling and inverse Fourier transform in conventional FP). Therefore, WFP is applicable for different variants of conventional FP, such as multiplexed FP [31,32] and extended FP for fluorescence imaging [5].
Although advantages in multiple ways over the conventional FP algorithm, WFP is limited in running efficiency, i.e., WFP needs more running time than AP. Therefore, shortening the running time of WFP is one of our future work. Utilizing accelerated gradient descent methods and introducing parallel computation techniques are two promising speeding up options.