1. Introduction

Phase retrieval has been fundamental to several disciplines, ranging from imaging, microscopy, crystallography and optics to astronomy [1–4]. It aims to recover an object only from its Fourier magnitudes. Without the Fourier phases, the recovery can be achieved via iterative algorithms when the Fourier magnitudes are sampled at a frequency sufficiently finer than the Nyquist interval [5]. In 1972, Gerchberg and Saxton developed an iterative algorithm for phase retrieval, utilizing the magnitude of an image and the Fourier magnitudes as constraints [6]. In 1982, Fienup generalized the Gerchberg-Saxton algorithm by developing two iterative algorithms: error reduction (ER) and hybrid input-output (HIO), which use a support and non-negativity as constraints in real space and the Fourier magnitudes in reciprocal space [7]. In 1998, Miao, Sayre and Chapman proposed, when the number of independently measured Fourier magnitudes is larger than the number of unknown variables associated with a sample, the phases are in principle encoded in the Fourier magnitudes and can be retrieved by iterative algorithms [5]. These developments finally led to the first experimental demonstration of coherent diffractive imaging (CDI) by Miao and collaborators in 1999 [8], which has stimulated wide spread research activities in phase retrieval, CDI, and their applications in the physical and biological sciences ever since [2, 9, 10].

For a finite object, when its Fourier transform is sampled at a frequency finer than the Nyquist interval (i.e. oversampled), mathematically it is equivalent to padding zeros to the object in real space. In another words, when the magnitude of the Fourier transform is oversampled, the correct phases correspond to the zero-density region surrounding the object, which is known as the oversampling theorem [5, 11]. The phase retrieval algorithms iterate between real and reciprocal space using zero-density region and the Fourier magnitudes as dual-space constraints. A support is typically defined to separate the zero-density region from the object. The support constraint, which is also known as the non-negativity constraint, is applied to the density inside the support. In the ER algorithm, the no-density region outside the support and the negative density inside the support are set to zero in each iteration [7]. The HIO algorithm relaxes the ER in the sense that it gradually reduces the densities that violate the support constraint instead of directly forcing them to zero [7]. This relaxation often leads to good reconstructions from noise-free patterns. However, in real experiments, the diffraction intensities, which are proportional to the square of Fourier magnitudes, are corrupted by a combination of Gaussian and Poisson noise and missing data. In the presence of experimental noise and missing data, phase retrieval becomes much more challenging, and the ER and HIO algorithms may only converge to sub-optimal solutions. Simply combining ER and HIO still suffers from stagnation and the iterations can get trapped at local minima [12]. To alleviate these problems, more advanced phase retrieval algorithms have been developed such as the shrink-wrap algorithm and guided HIO (gHIO) [4, 13]. In 2010, a smoothness constraint in real space was first introduced to improve the phase retrieval of noisy data [14]. Later, a noise robust framework was implemented for enhancing the performance of existing algorithms [15]. Recently, Rodriguez et al. proposed to impose the smoothness constraint on the no-density region outside the support by applying Gaussian filters [16]. The resulting oversampling smoothness (OSS) algorithm successfully reduces oscillations in the reconstructed image, and is more robust to noisy data than the existing algorithms.

Since phase retrieval can be cast as a non-convex minimization problem, many efforts have been made to study phase retrieval algorithms from the viewpoint of optimization. For example, Bauschke et al. [17] related HIO to a particular relaxation of the Douglas-Rachford algorithm [18] and introduced the hybrid projection reflection algorithm [17, 19]. Using similar ideas, researchers further proposed several projection algorithms such as iterated difference map [20] and relaxed averaged alternation reflection [21]. In [22], Chen and Fannjiang analyzed a Fourier-domain Douglas-Rachford algorithm for phase retrieval. By taking noise into account, the Wirtinger Flow [23] replaces the amplitude constraint by the intensity constraint which is then relaxed into a least square fidelity term that measures the misfit of measured Fourier intensity data, to which gradient descent is applied. Other methods in this line include alternating direction methods [24–26] that have been widely used in image processing, as well as lifting approaches [27] such as PhaseLift [28, 29] by Candès et al. and its variants [30, 31].

In this paper, we propose an optimization-based phase retrieval method, termed generalized proximal smoothing (GPS), which effectively addresses the noise in both real and Fourier spaces. Motivated by the success of OSS [16], GPS incorporates the idea of Moreau-Yosida [32, 33] regularization with heat kernel smoothing, to relax the support constraint into a continuous penalty term. We further relax the magnitude constraint into a least squares fidelity term, for de-noising in Fourier space. To minimize the resulting primal-dual formulation, GPS iterates back and forth between efficient proximal mappings of the two relaxed functions, respectively. Our experimental results using noisy experimental data of biological and inorganic specimens demonstrate that GPS consistently outperforms the state-of-the-art algorithms HIO and OSS in terms of both speed and robustness. We also refer readers to the recent paper [34] about training quantized neural networks, which shows another success of using Moreau-Yosida regularization to relax the hard constraint.

Notations

Let us fix some notations. For any complex-valued vectors u, v ∈ ℂⁿ, $\overline{u}$ is the complex conjugate of u, whereas $u^{*}:={\overline{u}}^{\top }$ is the Hermitian transpose. Re(u) and Im(u) are the real and imaginary parts of u, respectively.

$〈u,v〉:=u^{*}v={\displaystyle {\sum }_{i=1}^n\overline{u_i}{v}_i}$ is the Hermitian inner product of u and v. u ⊙ v is the element-wise (Hadamard) product of u and v given by (u ⊙ v)_i = u_i v_i. $\Vert u\Vert :=\sqrt{〈u,u〉}$ denotes the ℓ2 norm of u. Given any Hermitian positive semi-definite matrix K ∈ ℂ^n×n, we define ${\Vert u\Vert }_K:=\sqrt{〈u,Ku〉}$.

$\arg (u):={\tan }^{-1}\left(\frac{\mathrm{Im}(u)}{\mathrm{Re}(u)}\right)$ denotes the argument (or phase) of a non complex-valued vector u = Re(u) + i · Im(u). For convenience we set arg(0) = 1.

${ℐ}_{\mathcal{X}}$ is the indicator function of a closed set $\mathcal{X}\subset {ℂ}^n$ given by ${ℐ}_{\mathcal{X}}(x)=0$ if $x\in \mathcal{X}$ and ${ℐ}_{\mathcal{X}}(x)=\infty$ otherwise

${\text{proj}}_{\mathcal{X}}(u):=\arg {\min }_{v\in \mathcal{X}}\Vert v-u\Vert$ is the projection of u onto $\mathcal{X}$,

$\text{pro}{x}_f(u):=\arg {\min }_v\left\{f(v)+\frac{1}{2}{\Vert v-u\Vert }^2\right\}$ is the proximal mapping of the function f.

$h^{*}(y):=\underset{u\in \mathcal{C}}{\sup }\left\{\mathrm{Re}〈y,u〉-h(u)\right\}$ is the extended Legendre transformation to the complex plane of a real-valued function h.

2. The proposed model

We first fix some settings of the phase retrieval problem. Consider the reconstruction of a 2D image u defined on a discrete lattice Ω := {(i, j) :1 ≤ i ≤ n₁, 1 ≤ j ≤ n₂}.

For simplicity, we represent u in terms of a vector in ℝⁿ by the lexicographical order with n = n₁ × n₂. Then u_i represents the density of image at the i-th pixel. Due to oversampling, the object densities reside in a sub-domain S ⊂ Ω known as the support, and u is supposed to be zero outside S. Throughout the paper, we assume that the support S is centered around the domain Ω. The support constraint is $\mathcal{S}:=\{u\in {ℝ}^n:u_i\ge 0\text{if}i\in S,u_i=0\text{otherwise}\}$.

The Fourier magnitude data is obtained as b = |ℱu|, where $ℱ:{ℝ}^{n_1\times n_2}\to {ℂ}^{n_1\times n_2}$ is the discrete Fourier transform (DFT). We denote the magnitude constraint by $\mathcal{T}:=\{u\in {ℝ}^n:|ℱu|=b\}$.

In the absence of noise, phase retrieval (PR) problem is simply to

This amounts to the following composite minimization problem

2.1. A new noise-removal model

In real experiments, the Fourier data are contaminated by experimental noise. Moreover, the densities outside the support are not exactly equal to zero. In the noisy case, the image to be reconstructed no longer fulfills either the Fourier magnitudes or the support constraint. The ER algorithm (a.k.a. alternating projection), which performs alternatively projects onto these two constraints, apparently does not take care of the noise. The HIO “relaxes” the support constraint on densities wherever it is violated. This relaxation only helps in the noiseless case. In the presence of noise, the feasible set $\mathcal{S}\cap \mathcal{T}$ can be empty, and alternating projection method and its relaxation, like ER and HIO, may fail to converge and keep oscillating. The OSS [16] improves the HIO by applying extra Gaussian filters to smooth the densities outside the support at different stages of the iterative process. None of them, however, seems to properly address the corruption of the Fourier magnitudes.

Introducing the splitting variable z = ℱu ∈ ℂⁿ, we reformulate (2) as

In the presence of noise, we seek to relax the indicator functions f and g that enforce hard constraints into soft constraints. To this end, we first relax g into a least square fidelity, which is the sum of squared errors as follows

This fidelity term has been considered in the literature by assuming the measurements being corrupted by i.i.d. Gaussian noise; see [24] for example. In practice, we observe that it works well even with a combination of Gaussian and Poisson noise. This relaxation is related to the Moreau-Yoshida regularization [32, 33], which will be explained in detail in the appendix.

Following this line, we further relax $f={ℐ}_{\mathcal{S}}$ into

This constraint optimization can be solved by a splitting method, such as alternating direction method of multipliers (ADMM) [35–37] or primal dual hybrid gradient (PDHG) algorithm [38–41].

2.2. A primal-dual formulation

We solve the constraint optimization problem (6) by introducing the Lagrangian

For notational convenience, we denote ${\mathcal{S}}^{*}:=\{y\in {\mathcal{S}}^n:\mathrm{Re}(y)\le 0\text{on}S\}$ the dual constraint set of S. Our finalized min-max problem is thus given by

3. Generalized proximal smoothing (GPS) algorithm

We carry out the minimization of the saddle point problem (9) or (12) by a generalized primal dual hybrid gradient (PDHG) algorithm [38–41], which iterates

Moreover, we need to find the proximal mapping of $s{f}_G^{*}$ for updating y^k+1, which reduces to

Eq. (15) has closed-form solution only when G is a diagonal matrix. In the other case, we provide an approximation.

We devise two versions of GPS algorithm based on different choices of G.

3.1. Real space smoothing

One choice of G is G = γ D^TD. Here D is the discrete gradient operator, and then D^TD is the negative of discrete Laplacian. In this case,

The linear inversion is in fact the Laplacian smoothing that follows the projection to ensure the smoothness of the reconstructed image after each iteration. The real space smoothing is related to the diffusion process and can be approximated by a low-pass filter. Recall ${\mathcal{G}}_t(x):=\frac{1}{{(4\pi t)}^{n/2}}\exp \left(-\frac{{\Vert x\Vert }^2}{4t}\right)$ is a (heat) Gaussian kernel and its Fourier transform W_t (ξ) = exp(−4π²t|ξ|²) is a normalized Gaussian function. Then replacing the linear inversion by the Gaussian kernel and convolution give a fast approximated implementation of Eq. (17) for small γ

The projection on ${\mathcal{S}}^{*}$ can be computed as follows

Here x⁻ := min(0, x) and x⁺ := max(0, x) for x ∈ ℝ. The remaining convolution can be done via the efficient DFT by computing

Using Eq. (14), Eq. (19), and Eq. (20), we arrive at our first algorithm

Algorithm 1: GPS-R features low-pass filters for smoothing. Here we abuse notation γ to imply the filter. Inspired by OSS [16], we choose a sequence of increasing spatial frequency filters {γ_l} (a sequence of finer filters). In our experiments, we do 1000 iterations with a total of 10 stages. Each stage contains 100 iterations, in which the filter frequency is held constant. We monitor the R-factor (relative error) during the iterative process, which is defined as

The reconstruction with minimal R_F at each stage is fed into the next stage. By applying the smoothing on the entire domain, GPS-R can remove noise in real space and obtain the spatial resolution with fine features.

Algorithm 1. GPS-R: GPS with smoothing in real space.

View Table | View all tables in this article

3.2. Fourier space smoothing

Another simple choice of G is the diagonal matrix G = γ Diag(r ⊙ r) where r ∈ ℝⁿ and r_i is the distance in the original 2D lattice between the i-th pixel and the center of image. Note that G is not invertible since r_i = 0 for the pixel at the center. By Eq. (41),

So f_G(u) is a weighted sum of squares penalty on u. The weight is inversely proportional to the squared radius, which is infinity for density in the center. The further the density off the center, the smaller the penalty for the support constraint being violated.

By Parseval’s identity, for a square-integrable function u, we have

Therefore, by Eq. (10),

This means that we are smoothing f^*(y) by regularizing with the ℓ₂ gradient of Fourier coefficients of y. We thus refer to it as Fourier space smoothing.

Since G is diagonal, Eq. (15) has the closed-form solution

Define $W_{\frac{s\gamma }{4{\pi }^2}}:=\exp (-s\gamma r^2)$. The above quantity can be approximated by a direct multiplication for small γ. We update y^k+1 as

GPS with smoothing in Fourier space (GPS-F) is summarized in Algorithm 2.

Algorithm 2. GPS-F: GPS with smoothing in Fourier space.

View Table | View all tables in this article

We illustrate the core of GPS algorithm as a flowchart in Fig. 1

 

Fig. 1: Flowchart of GPS algorithm for stage 1. For simplicity, we denote $ϵ=\frac{\sigma }{\sigma +1}$ in Eq. (14). The input z⁰ is initialized with random phases and magnitudes are equal the measurements, while y⁰ is a zero vector. The ${\text{proj}}_{{\mathcal{S}}^{*}}$ and the convolution ∗ follows Eq. (19) and Eq. (20) respectively. z^k and y^k with the smallest R_F will be fed as an input into the next stage where a Gaussian kernel ${\mathcal{G}}_{\gamma }$ with larger γ is applied.

Download Full Size | PDF

3.3. Incomplete measurements

In practice, not all diffraction intensities can be experimentally measured. For example, to prevent a detector from being damaged by an intense X-ray beam, either a beamstop has to be placed in front of the detector to block the direct beam or a hole has to be made at the center of the detector, resulting in missing data at the center [42]. Furthermore, missing data may also be present in the form of gaps between detector panels. For incomplete data, the alternating projection algorithms skip the projection onto the magnitude constraint in this region. Similarly, we only apply the relaxation $g_{\sigma }=\frac{1}{2\sigma }{{\displaystyle {\sum }_i\left|\left|z_i\right|-b_i\right|}}^2$ on the known data for GPS. A simple exercise shows that

4. Experimental results

4.1. Reconstruction from simulated data

We expect GPS to be a reliable PR algorithm in the reconstruction of the images of weakly scattering objects, in particular biological specimens, which have become more popular [43]. Since OSS has been shown to consistently outperform ER, HIO, ER-HIO, NR-HIO [16], we perform both quantitative and qualitative comparisons between GPS and OSS.

Due to the discrete nature of photon counting, experimentally measured data inherently contain Poisson noise that is directly related to the incident photon flux on the object. In addition to Poisson noise, the data is also corrupted with zero-mean Gaussian noise to simulate readout from a CCD. Any resulting negative values are set to zero. Therefore, an accurate simulation of b_i can be calculated as

For simulation, the Fourier magnitudes of a vesicle model are corrupted with Poisson and Gaussian noises by Eq. (29). The corresponding relative error is approximated R_noise ≈ 6%.

In some cases, the reconstructed image yields a small relative error R_F but has low quality. This is the issue of over-fitting, an example of which can be seen in certain reconstructions using ER-HIO [16]. Smoothing is a technique to avoid data over-fitting. To validate results and show that our algorithm does not develop over-fitting, we measure the difference between the reconstructed image and the model by

In addition, we also look at the residual $\text{Res}=\left|\right|ℱu|-b||$ which measures the difference between the Fourier magnitudes of the reconstructed image and the experimental measurements. The residual can validate the noise-removal model, telling how much noise is removed. Figure 2 shows the reconstruction of vesicle model from simulated noisy data using HIO, OSS, GPS-R, and GPS-F. GPS-F and GPS-R obtain lower R_F and R_real than OSS. Moreover, GPS-F can get very close to the ground truth with R_F =5.90% and R_real = 0.7%. In addition to lower R values, GPS-R and GPS-F converge to zero outside the support. They both obtain lower residuals than OSS, specifically GPS-F produces the least residual. If we choose larger parameter for the ℓ₂ gradient regularizer in Fourier space, we will get a smoother residual. Overall for realistic Gaussian and Poisson noise in Fourier space, GPS-F is a suitable noise-removal model.

 

Fig. 2: First row: vesicle model with log-scaled diffraction pattern (left), zoom-in image (center) and residual (right). Second row: HIO: R_F = 12.87%, R_real = 21.14%. Third row: OSS: R_F = 6.08%, R_real = 3.59%. Fourth row: GPS-R R_F = 5.90%, R_real = 2.85%. Fifth row: GPS-F R_F = 5.89% and R_real = 0.7%. Third column: residual.

Download Full Size | PDF

Figure 3 shows the histogram and the convergence of R_F and R_real on 100 independent, randomly seeded runs using HIO, OSS and GPS on the simulated vesicle data. The histogram shows that GPS is more consistent and robust than OSS. It has a higher chance to converge to good minima with lower R_F and R_real than OSS. Furthermore, R_F and R_real of OSS scatter widely due to the initial value dependency. In contrast, GPS is more selective and less dependent on initial values. R_F and R_real of GPS are seen at a lower mean minimum with less variance.

 

Fig. 3: Histogram (first row), and convergence (second row) of R_F and R_real on Vesicle data using HIO, OSS, GPS. GPS consistently produces smaller R_F and R_real than HIO or OSS. Moreover, GPS converges fastest with the fewest oscillations.

Download Full Size | PDF

Similar to HIO and ER, OSS keeps oscillating until a finer low-pass filter is applied. In contrast, GPS converges faster and is less oscillatory than OSS. In the presence of noise, alternating projection methods (ER, HIO, OSS) keep oscillating but do not converge. Applying smoothing and replacing the measurement constraint by the least squares fidelity term $g_{\sigma }(z)=\frac{1}{2\sigma }{\Vert \left|z\right|-b\Vert }^2$ helps to reduce the oscillations; hence, the method can converge to a stationary point. Note that larger σ reduces more oscillations, but also decreases the chance to escape from local minima. Alternating projection methods have σ = 0 since they impose measurement constraints. GPS obtains both smaller R_F, R_real, and lower variance. Even though R_F are close to each other, R_real of GPS is much smaller than OSS. This means GPS recovers the vesicle cell with higher quality than OSS.

Next, we test the consistency of GPS algorithm for vesicle model. Since OSS has been shown to outperform HIO, ER-HIO and HR-HIO [16], it suffices to carry out the comparison among HIO, OSS and GPS. We look for the mean and variance of the best five of 500 independent reconstructions by HIO, OSS, GPS-F and GPS-R (lowest R_F). Figure 4 shows that GPS has much smaller variance than OSS and HIO. These tests shows that GPS is more reliable and consistent than OSS and HIO.

 

Fig. 4: The (scaled) means and variances of the best five of 500 independent reconstructions by HIO, OSS, GPS-F and GPS-R

Download Full Size | PDF

4.2. Reconstructions from experimental data

4.2.1. S. pombe yeast spore

To demonstrate the applicability of GPS to biological samples, we do phase retrieval on the diffraction pattern in Fig. 5 taken of a S. Pombe yeast spore from an experiment done using beamline BL29 at SPring-8 [45]. We do 500 independent, randomly seeded reconstructions with each algorithm and record R_F, excluding the missing center. We choose default parameters for these experiments: t = 1, s = 0.9. The sequence of low-pass filters are chosen to be the same as in OSS [16]. For the first 400 iterations, σ = 0.01, then is increased to σ = 0.1 for the remaining 600 iterations. The left column of Fig. 5 is the mean of the best 5 reconstructions obtained by the respective algorithm. The right column shows the variance of the same 5 reconstructions. It is evident from the variance that GPS achieves more consistent reconstructions. Figure 6 shows the histogram and convergence of R_F. We can conclude that not only are GPS-R results the most consistent, but also the most faithful to the experimental data.

 

Fig. 5: S. pombe yeast spore log-scaled diffraction pattern, size 500 × 500 (top). Means (first row) and variance (second row) are computed from the best 5 of 500 independent reconstructions. HIO: R_F = 15.697% ± 0.526%, OSS: R_F = 9.775% ± 0.202%, and GPS: R_F = 8.672% ± 0.025%.

Download Full Size | PDF

 

Fig. 6: Top: histogram of R_F in 500 independent runs (top). Bottom: the convergence curves of a singe construction of R_F on S. pombe yeast spore by GPS-RF, OSS, and HIO.

Download Full Size | PDF

4.2.2. Nanorice

To demonstrate the generality of GPS, we also do testing with experimental data of inorganic samples. The diffraction patterns shown in the top row of Fig 7 from ellipsoidal iron oxide nanoparticles (referred to as ‘nanorice1’and ‘nanorice2’) were taken at the AMO instrument at LCLS at an X-ray energy of 1.2keV [46]. This data is freely available online on the CXIDB [47]. We choose default parameters for these experiments: t = 1, s = 0.9. The sequence of low-pass filters are chosen to be the same as in OSS [16]. The fidelity parameter σ is chosen small for the first 800 iterations, specifically σ ∈ [0, 0.01], to produce oscillations which is necessary for the algorithm to skip bad local minima. Once the reconstruction arrives at a good local minimum region, we increase σ to reduce oscillations. This later value of σ depends on noise level and data. We test different values of σ and σ = 1 has been found to be the optimal for both nanorice data. Figure 7 shows OSS (second row) and GPS-F(third row) reconstructions of the two nanorice particles. Figure 8 shows again that GPS obtains more consistent and faithful reconstructions as compared to those obtained by OSS. GPS-F with σ = 1 converges to lower relative error than OSS at all times. OSS cannot get lower relative error because σ = 0 does not work for this case. In general, alternating projection methods do not treat noise correctly. For example, in this case, HIO keeps oscillating but does not converge. Therefore, its results are omitted here. The better approach, OSS model, can reduce oscillations by smoothing but this is not enough. In contrast, the least squares $g_{\sigma }(z)=\frac{1}{2\sigma }{\Vert \left|z\right|-b\Vert }^2$ of GPS works for noise removal since relaxing the constraints allows GPS to reach lower relative error. The values of σ depend on noise level and type. To optimize the convergence of GPS-F on nanorice2, we apply σ = 0.01 for the first 400 iterations, σ = 0.1 for the next 300 iterations, and σ = 1 for the last 300 iterations. This test shows the effect of σ on the convergence. OSS (σ = 0) oscillates the most. GPS with σ = 0.01, 0.1, 1 oscillates less and less. As σ increases, GPS also gets to lower R_F. The algorithm finally reaches a stable minimum as σ goes up to 1. Continuing to increase σ does not help with R_F. Choosing large σ in the beginning may reduce oscillations but also limit the mobility to skip local minima. We recommend start with small σ and then gradually increase it until the iterative process reaches a stable minimum.

 

Fig. 7: Diffraction pattern of nanorice1 and nanorice2 253 × 253 (first column) and reconstructions using OSS: R_F = 18.23%, 16.32% and GPS-F: R_F = 17.40%, 15.83% respectively. GPS obtains less noise on the boundary and lower relative error R_F.

Download Full Size | PDF

 

Fig. 8: Histogram (first row) and convergence (second row) of OSS and GPS-F on nanorice1 (first column) and nanorice2 (second column). The results of HIO are omitted due to lack of convergence.

Download Full Size | PDF

5. Conclusion

In this work, we have developed a fast and robust phase retrieval algorithm GPS for the reconstruction of images from noisy diffraction intensities. Similar to [34], the Moreau-Yosida regularization was used to relax the hard constraints considered in the noiseless model. GPS utilizes a primal-dual algorithm and a noise-removal technique, in which the [2 gradient smoothing is effectively applied on either real or Fourier space of the dual variable. GPS shows more reliable and consistent results than OSS, HIO for the reconstruction of weakly scattered objects such as biological specimens. Looking forward, we aim to explore the role of dual variables in non-convex optimization. Smoothing the dual variable, which is equivalent to smoothing the gradient of convex conjugate, represents a new and effective technique that can in principle be applied to other non-smooth, non-convex problems.