Semi-implicit relaxed Douglas-Rachford algorithm (sDR) for ptychography

Minh Pham; Arjun Rana; Jianwei Miao; Stanley Osher

doi:10.1364/OE.27.031246

1. Introduction

Since the first experimental demonstration in 1999 [1], coherent diffraction imaging (CDI) has been a rapidly growing field due to its broad potential applications in the physical and biological sciences [2–5]. A fundamental problem of CDI is the phase problem, that is, the measured diffraction pattern only contains the Fourier magnitudes without phase information. In the original demonstration of CDI, phase retrieval was performed by measuring the diffraction pattern from a finite object. If the diffraction intensity is sufficiently oversampled [6], the phase information can be directly retrieved by using iterative algorithms [7–11]. Ptychography, a powerful scanning CDI method, relieves the finite object requirement by performing 2D scanning of an extended object with an illumination probe. The object is scanned such that adjacent scan positions are partially overlapping, and a diffraction pattern is collected at each scan position [12,13]. The overlap of the probes not only reduces the oversampling requirement, but also improves the convergence speed of the iterative process. By taking advantage of ever-improving computational power and advanced detectors, ptychography has been applied to study a wide range of samples using both coherent x-rays and electrons [2,5,14–23]. More recently, a time-domain ptychography method was developed by introducing a time-invariant overlapping region as a constraint, allowing the reconstruction of a time series of complex exit wave of dynamic processes with robust and fast convergence [24].

Algorithms for ptychography have been studied exhaustively in theory and practice. The majority follow non-convex optimization approaches [25–27], while a few follow convex relaxation [28]. In recent years, powerful ptychographic algorithms have been developed to handle partial coherence [29], solve for the probe [13,30,31], correct positioning errors [32–34], reduce noise [35,36], and deal with multiple scattering [37,38].

These iterative algorithms can be generally divided into three classes: i) conjugate gradient (CG) [30,34], ii) extended ptychographic iterative engine (ePIE) [31], and iii) difference map (DM) [13], whereas the last two have a close relationship. ePIE is an alternating projection algorithm, while DM is built on both projection and reflection which is believed to provide a momentum to speed up the convergence. The relaxed average alternating reflection (RAAR) method [8] is a relaxation of DM and has been shown to be effective in phase retrieval [39]. All algorithms except ePIE take a global approach, i.e. using the entire collection of diffraction patterns to perform an update of the probe and object in each iteration. In contrast, ePIE goes through the measured data sequentially to refine the probe and object. The advantage of global updating methods is the ability of parallelization which significantly reduces the computation time and allows efficient implementation on graphical processing units (GPUs) and clusters. However, the larger memory requirement is a major drawback of these methods. On the other hand, sequential method ePIE is still being used effectively due to its memory efficiency and simplicity. One drawback of sequential ePIE is that it requires a small update step size for stable convergence. Limiting the convergence speed, the step size restriction may cause divergence if violated. To fix this issue, rPIE was proposed, in which regularization is used for stability [40]. The significant results also show that rPIE obtains a larger field of view (FOV) than ePIE. An attempt to implement DM sequentially is reported [41] but the algorithm is slightly modified by enforcing the Fourier magnitude constraint update twice.

An exceptional ability of the stochastic sequential update method is the power to escape local minima due to stochasticity, which plays a crucial role in non-convex optimization. Especially when measurements are sparse (less overlap), the method has a higher chance to pass local minima than CG and DM. In this light, we motivate the utilization of non-convex optimization tools to improve the robustness and the convergence of sequential updating methods. Our proposed algorithm semi-implicit relaxed Douglas-Rachford (sDR) incorporates two techniques. The first one modifies the update of the probe and object as the algorithm iterates in a semi-implicit fashion known as Proximal Mapping. The second technique is the implementation of relaxed Douglas-Rachford, a generalized version of DM and RAAR, on the local scale. In addition to our formulation, we show that DM and RAAR can also be implemented locally, similarly to ePIE.

2. Proposed algorithm

Given N measured diffraction patterns at N positions, the ptychographic algorithm aims to find a 2D object O and a 2D probe P that satisfy the overlap constraint and the Fourier magnitude constraint

(1)$$|\mathcal{F} (P O_n) | = \sqrt{I_n} \quad \mbox{for} \quad n=1,..,N.$$

Where $O_n$ is the object at the $n^{th}$ scan position. For simplicity, we omit spatial variables and use notations $P, \, O_n$ and $I_n$ for both continuous and discrete cases. The absolute value, multiplication, division, conjugate, and square root operators are applied element-wise on $P, \, O_n$ and $I_n$ which represent 2D complex matrices of the same size in the discrete case. Alternatively, $O_n$ represents the object $O$ restricted to a sub-domain $\Omega _n$. The overlap constraint can be mathematically interpreted as

(2)$$O(x+r_n) = O_n(x) \quad \mbox{if} \quad r_n+x \in \Omega_n \quad \mbox{for} \quad n=1,\ldots,N$$

where $\{r_n\}_{n=1}^N$ are displacement vectors. In short notation, we write $O_n = O|_{\Omega _n}$. We find a better representation of the problem by introducing the exit wave variable $\Psi = PO$. By denoting the Fourier measurement constraint set $\mathcal {T}$ and the overlap object constraint set $\mathcal {S}$, we have

(3)$$\begin{aligned} \mathcal{T} & := {\big\{} \Psi = \{\Psi_n\}_{n=1}^N: | \mathcal{F} \Psi_n| = \sqrt{I_n} \mbox{ for } n=1,\ldots,N {\big\}} \\ \mathcal{S} & : = {\big\{} \Psi = \{\Psi_n\}_{n=1}^N: \exists P, O \mbox{ s.t. } \Psi_n = PO_n \mbox{ for } n=1,\ldots,N \}{\big\}}. \end{aligned}$$

Then we write the ptychography problem in a minimization fashion

(4)$$\min_{ \Psi} \quad i_{\mathcal{S}}(\Psi) + i_{\mathcal{T}}(\Psi)$$

where $i_{\mathcal {S}}(\Psi )$ and $i_{\mathcal {T}}(\Psi )$ are the indicator functions of sets $\mathcal {S}$ and $\mathcal {T}$ respectively, defined as

(5)$$i_{\mathcal{S}} (\Psi) = \begin{cases} 0 & \Psi \in \mathcal{S} \\ \infty & \textrm{otherwise} \end{cases} $$

To solve Eq. (4), an alternating projection method is proposed. At each iteration, we select a random position $n$ and update $\Psi _n$

(6)$$\Psi_n' = \Pi_{\mathcal{T}} ( \Psi_n^k ) = \mathcal{F}^{{-}1} {\big(} \sqrt{I_n} \arg (\mathcal{F} \Psi_n^k) {\big)} $$

(7)$$ \{ P^{k+1},O_n^{k+1} \} = \mathop{\mathrm{argmin}}_{P,O_n} \frac{1}{2} \| PO_n - \Psi_n' \|^2 $$

(8)$$\Psi_n^{k+1} = P^{k+1} O_n^{k+1} $$

The Frobenius norm is used in this minimization problem and entire paper unless a different norm is specified. The minimization of Eq. (7) is difficult due to instability. One way to solve this non-convex problem is to minimize each variable while fixing the other ones.

(9)$$\begin{aligned} O_n^{k+1} & = \mathop{\mathrm{argmin}}_{O_n} \frac{1}{2} \| P^k O_n - \Psi_n' \|^2 = \frac{\Psi_n'}{P^k} \\ P^{k+1} & = \mathop{\mathrm{argmin}}_{P} \frac{1}{2} \| P O_n^{k+1} - \Psi_n' \|^2 = \frac{\Psi_n'}{O_n^{k+1}} \end{aligned}$$

This approach is unstable because of the division. A cut-off method is used to avoid the divergence and zero-division. A modification is recommended by adding a penalizing least square error term (i.e. regularizer)

(10)$$\{ P^{k+1},O_n^{k+1} \} = \mathop{\mathrm{argmin}}_{P,O_n} \frac{1}{2} \| P O_n - \Psi_n' \|^2 + \frac{1}{2 s} \| P-P^k \|^2 +\frac{1}{2t} \|O_n-O_n^k\|^2$$

The idea of regularization appears throughout the literature such as proximal algorithms [42,43]. Equation (10) is more reliable to solve than Eq. (7) but is still very expensive since the variables are coupled. $P^{k+1}$ and $O_n^{k+1}$ can be solved via a Backward-Euler system derived from Eq. (10).

(11)$$\begin{aligned} O_n^{k+1} & = O_n^k -t \overline{P^{k+1}} {\big(} P^{k+1} O_n^{k+1} - \Psi_n' {\big)} \\ P^{k+1} & = P^k - s \overline{O_n^{k+1}} {\big(} P^{k+1} O_n^{k+1} - \Psi_n' {\big)} \end{aligned}$$

ePIE proposes a simple approximation by linearizing the system so that it can be solved sequentially.

(12)$$\begin{aligned} \quad O_n^{k+1} & = O_n^k -t \overline{P^k} {\big(} P^k O_n^k - \Psi_n' {\big)} \\ \quad P^{k+1} & = P^k - s \overline{O_n^{k+1}} {\big(} P^k O_n^{k+1} - \Psi_n' {\big)} \end{aligned}$$

Indeed, ePIE is gradient descent method applied on the minimization of Eq. (7). The system is solved by alternating direction methods (ADM) [44]. The remaining part is to choose appropriate step sizes $t$ and $s$ to ensure stability. ePIE suggests $t = \beta _O/ \|P^k\|_{\max }^2$ and $s = \beta _P/ \|O^{k+1}\|_{\max }^2$ where $\beta _O, \; \beta _P \in (0,1)$ are normalized step sizes. The max matrix norm is the element-wise norm, taking the maximum in absolute values of all elements in the matrix. The final version of ePIE is given by

(13)$$\begin{aligned} O_n^{k+1} & = O_n^k - \left. \beta_O { \overline{P^k} {\big(} P^k O_n^k - \Psi_n' {\big)} } \middle/ {\|P^k\|_{\max}^2} \right. \\ P^{k+1} & = P^k - \left. \beta_P { \overline{O_n^{k+1}} {\big(} P^k O_n^{k+1} - \Psi_n' {\big)} } \middle/ {\|O_n^{k+1}\|_{\max}^2} \right. \end{aligned}$$

We will exploit the structure of Eq. (11) to give a better approximation.

2.1 Semi-implicit algorithm

We replace the minimization of Eq. (10) by two steps

(14)$$\begin{aligned} \mbox{Step 1:} \quad O_n^{k+1} & = \mathop{\mathrm{argmin}}_{O_n} \frac{1}{2} \| P^k O_n - \Psi_n' \|^2 + \frac{1}{2t} \| O_n - O_n^k \|^2 \\ \mbox{Step 2:} \quad P^{k+1} & = \mathop{\mathrm{argmin}}_{P} \frac{1}{2} \| P O_n^{k+1} - \Psi_n' \|^2 + \frac{1}{2s} \| P - P^k \|^2 \end{aligned}$$

This results in a better approximation to the linearized system of Eq. (12) and is simpler than the Backward-Euler Eq. (11)

(15)$$\begin{aligned} O_n^{k+1} & = O_n^k -t \overline{P^k} {\big(} P^k O_n^{k+1} - \Psi_n' {\big)} \\ P^{k+1} & = P^k - s \overline{O_n^{k+1}} {\big(} P^{k+1} O_n^{k+1} - \Psi_n' {\big)} \end{aligned}$$

In this uncoupled system, we can derive a closed form solution for each sub-problem.

(16)$$\begin{aligned} O_n^{k+1} & = \left. {\big(} O_n^k + t \overline{P^k} \Psi_n' {\big)} \middle/ {\big(}1 + t |P^k|^2 {\big)} \right.\\ P^{k+1} & = \left. {\big(} P^k + s \overline{O_n^{k+1}} \Psi_n' {\big)} \middle/ {\big(} 1 + s|O_n^{k+1}|^2 {\big)} \right. \end{aligned}$$

By choosing the step sizes $s$ and $t$ as in the ePIE algorithm and normalizing the parameters $\beta _O$ and $\beta _P$, we obtain

(17)$$\begin{aligned} O_n^{k+1} & = \frac{(1-\beta_O)\|P^k\|_{\max}^2 O_n^k + \beta_O \overline{P^k} \Psi_n' }{(1-\beta_O)\|P^k\|_{\max}^2 + \beta_O |P^k|^2} \\ P^{k+1} & = \frac{(1-\beta_P)\|O_n^{k+1}\|_{\max}^2 P^k + \beta_P \overline{O_n^{k+1}} \Psi_n' }{(1-\beta_P)\|O_n^{k+1}\|_{\max}^2 + \beta_P |O_n^{k+1}|^2} \end{aligned}$$

$O_n^{k+1}$ is updated as a weighted average between the previous value $O_n^k$ and $\displaystyle \frac {\Psi _n^k}{P^k}$. The object update can be rewritten as

(18)$$O_n^{k+1} = O^k + \beta_O \frac{ \overline{P^k} {\big(} \Psi_n' - \Psi_n^k {\big)} } {(1-\beta_O)\|P^k\|_{\max}^2 + \beta_O |P^k|^2}$$

Recall that the object update by rPIE has a similar form [40]

(19)$$O_n^{k+1} = O^k + \frac{ \overline{P^k} {\big(} \Psi_n' - \Psi_n^k {\big)} } {\alpha \|P^k\|_{\max}^2 + (1-\alpha) |P^k|^2}$$

It is clearly to see that $\alpha = 1-\beta _O \in (0,1)$ and rPIE does not have the parameter $\beta _O$ in front of the fraction. i.e. rPIE has a larger step size for the object update than our proposed algorithm sDR which will be discussed later. Consequently rPIE converges faster. However, it is still unable to escape local minima when being trapped. This “weighted average” regularization is more mathematically correct and enhances the algorithm’s stability.

In the next section, we apply the Douglas-Rachford (DR) algorithm to solve for the exit wave $\Psi$ and improve the ability to escape local minima.

2.2 Relaxed Douglas-Rachford algorithm

The Douglas-Rachford algorithm was originally proposed to solve the heat conduction problem [45], which represents a composite minimization problem

(20)$$\min_{\Psi} \quad f(\Psi) + g(\Psi)$$

The iteration consists of

(21)$$\Psi^{k+1} = \Psi^k + \mathrm{prox}_{tf} {\big(} 2\; \mathrm{prox}_{tg} (\Psi^k) - \Psi^k {\big)} - \mathrm{prox}_{tg} (\Psi^k)$$

Over the past decades, this accelerated convex optimization algorithm has been exhaustively studied in both theory and practice with many applications [46–51]. Here we apply the algorithm to ptychographic phase retrieval. Note that the Douglas-Rachford algorithm reduces to Difference Map (DM) when $f=i_{\mathcal {T}}$ and $g=i_{\mathcal {S}}$ are characteristic functions of constraint sets $\mathcal {T}$ and $\mathcal {S}$, respectively

(22)$$\Psi^{k+1} = \Psi^k + \Pi_{\mathcal{T}} {\big(} 2\; \Pi_{\mathcal{S}} (\Psi^k) - \Psi^k {\big)} - \Pi_{\mathcal{S}} (\Psi^k)$$

The reflection operator $2\Pi _{\mathcal {S}}-I$ accelerates the convergence speed in the convex case and helps to escape local minima in the non-convex case. Momentum also compensates for the step size of the object update in Eq. (18) and speeds up convergence. However this momentum, caused by reflection, might be too large and causes divergence or oscillation. Especially in the case of high sparsity (less overlap), momentum can also cause over-fitting since the ptychography problem is no longer over-constrained. Therefore, we relax the reflection by introducing a relaxed reflection parameter $\sigma \in [0,1]$

(23)$$\Psi^{k+1} = \mathrm{prox}_{tf} \Big( (1+\sigma)\; \mathrm{prox}_{tg} (\Psi^k) - \sigma\Psi^k \Big) + \sigma \Big(\Psi^k - \mathrm{prox}_{tg} (\Psi^k) \Big)$$

Since the experimental measurements are contaminated by noise, a direct projection of the measurement constraint is not an appropriate approach. We thus relax the Fourier magnitude constraint by a least square penalty

(24)$$\min_{ \Psi_n } \quad \sum_{n=1}^N \| |\mathcal{F} \Psi_n | - \sqrt{I_n} \|^2 +i_{\mathcal{S}}(\Psi_n)$$

Recall that $\mathrm {prox}_{tf}(\Psi )$ has a closed form solution

(25)$$\begin{aligned} \mathrm{prox}_{tf}(\Psi^k) & = \mathop{\mathrm{argmin}}_{\Psi} \frac{1}{2} \| | \mathcal{F} \Psi| - \sqrt{I_n} \|^2 + \frac{1}{2t} \| \Psi - \Psi^k \|^2 \\ & = \left. {\Psi^k + t \mathcal{F}^{{-}1} \Big[\sqrt{I_n} \arg(\mathcal{F}\Psi) \Big] } \middle/ {(1 + t)} \right. = (1-\tau) \Psi^k + \tau \mathcal{F}^{{-}1} \Big[\sqrt{I_n} \arg(\mathcal{F}\Psi) \Big] \end{aligned}$$

where $\tau = t/(1+t) \in (0,1]$ is the normalized step size. We call $\tau$ the relaxed modulus constraint. Combining this result with DM, we obtain

(26)$$\begin{aligned}\Psi^{k+1} &= (1-\tau) \Big( (1+\sigma)\Pi_{\mathcal{S}} (\Psi^k) - \sigma \Psi^k \Big) +\tau \Pi_{\mathcal{T}} \Big( (1+\sigma)\Pi_{\mathcal{S}} (\Psi^k) - \sigma\Psi^k \Big) + \sigma \Big( \Psi^k - \Pi_{\mathcal{S}}(\Psi^k) \Big)\\ &= \tau \Big( \sigma \Psi^k + \Pi_{\mathcal{T}} {\big(} (1+\sigma)\Pi_{\mathcal{S}} (\Psi^k) - \sigma \Psi^k {\big)} \Big) + \Big(1-\tau(1+\sigma) \Big) \Pi_{\mathcal{S}}(\Psi^k) \end{aligned}$$

When we let $\beta = 1-\tau$ and $\sigma =1$, the update reduces to RAAR. To end this section, we conclude that relaxed Douglas-Rachford is a generalized version of RAAR. We now move to our main algorithm.

2.3 Semi-implicit relaxed Douglas-Rachford algorithm (sDR)

In a combination of the semi-implicit algorithm and relaxed Douglas-Rachford algorithm, we propose the sDR algorithm, shown in Algorithm 1 and presented as a flow chart in Fig. 1.

Fig. 1. Flow chart of the sDR algorithm.

Download Full Size | PDF

In most ptychography experiments, prior information about the probe is given and can be used as an initial guess. This knowledge is very helpful in non-convex optimization. In light of stochastic gradient descent methods, small step size $\beta _P$ performs efficiently and better since the probe starts from a good initial guess and gets updated more often than the object. The performance of forward Euler and semi-backward Euler methods is not very different because of this small step size. To keep simplicity and efficiency, we only apply the semi-implicit method on $O_n^k$ while $P^k$ can be integrated with the regular gradient descent method. At each iteration, sDR processes all diffraction patterns in a random sequence and updates the probe and object.

Comparing memory usage with ePIE, sDR stores an extra variable $\{Z_n\}_{n=1}^N$ which has the same size as the diffraction patterns. sDR uses more computations than ePIE due to the introduction of momentum. These extra computations are simple element-wise additions and divisions. sDR approximately spends 20% more computation time than ePIE.

For parameter selections, we choose a small relaxed modulus constraint $\tau$. Specifically $\tau = 0.1$ in all experiments. The values of the relaxed reflection parameter $\sigma$ depends on the specific problem. In many cases, $\sigma =1$ works very well (full reflection). However in noisy or low overlap data, large $\sigma$ might cause divergence. A smaller $\sigma$ is then preferred in these cases such as $\sigma =0.5$. Similar to the choice of $\sigma$, we select large $\beta _O$ in many cases, for example $\beta _O=0.9$, and smaller $\beta _O$ in noisy data. On the other hand, $\beta _P$ is selected differently. It needs not to be large since the probe gets updated more often than the object. A value of $\beta _P=0.1$ can start the reconstruction and decrease as a function of iteration. In our case, we choose a square root decreasing function

(27)$$\beta_O^k = \beta_O^0 \sqrt{ \left. k \middle/ (\textrm{K}-k) \right. }$$

where $\beta _O^0$ is the initial $\beta _O$. This adaptive step-size has been introduced as a strategy for noise-robust Fourier ptychography [52].

3. Experimental results

3.1 Reconstruction from simulated data

To examine the sDR algorithm, we simulate a complex object of $128 \times 128$ pixels with a cameraman and a pepper images representing the amplitude and the phase, respectively (Fig. 2). The circular aperture is chosen as probe with a radius of 50 pixels. We raster scan the aperture over the object with a step size of 35 pixels, resulting in 4x4 scan positions. The overlap is therefore 56.4% which is approximately the lower limit for ePIE to work in this simulated data test. Poisson noise is added to the diffraction patterns with a flux of $1\times 10^8$ photons per scan position. We use $R_{noise}$ to quantify the relative error with respect to the noise-free diffraction patterns

(28)$$R_{noise} = \frac{1}{N} \sum_{n=1}^N \left. \| |\mathcal{F} ( P^0 O^0_n ) | - \sqrt{I_n} \|_{1,1} \middle/ \| \sqrt{I_n} \|_{1,1} \right.$$

where $P^0$ and $O^0$ is the noise-free model and the $L_{1,1}$ matrix norm represents the sum of all elements in absolute value of the matrix. The above flux results in $R_{noise} = 3.73\%$. Figure 3 shows that three algorithms (ePIE, rPIE and sDR) all successfully reconstruct the object in the case where the overlap between adjacent positions is high and the noise level is low.

Fig. 2. A simulated complex object with the amplitude being a camera man image shown in (a) and the phase being a pepper image shown in (b).

Download Full Size | PDF

Fig. 3. The reconstructions of ePIE, rPIE and sDR of a complex object consisting of $128\times 128$ pixels, a scan step size of 35 pixels and $4\times 4$ diffraction patterns. Poisson noise was added to the diffraction patterns with $R_{noise} = 3.73\%$. (a-c) The amplitude and (d-f) the phase of ePIE, rPIE and sDR reconstructions, respectively.

Download Full Size | PDF

Next, we apply the three algorithms to the reconstruction of sparse data, which is centrally important to reducing computation time, data storage requirements and data acquisition time. We increase the scan step size to 50 pixels while keeping the same field of view, which reduces the number of diffraction patterns to $3\times 3$. Consequently, the overlap is reduced to 39.1%. Not only is the overlap between adjacent positions low, but the total number of measurements is also small, creating a challenging data set for conventional ptychographic algorithms. $\beta _O = 0.7, \; \beta _P^0 = 0.01$ and $\sigma = 0.9$ are used for sDR in this test. Figure 4 show that sDR can work well with sparse data, while ePIE and rPIE fail to reconstruct the object faithfully.

Fig. 4. Ptychographic reconstructions of sparse data by ePIE, rPIE and sDR. The data consist of $3\times 3$ diffraction patterns with a scan step of 50 pixels. Poisson noise was added to the diffraction patterns with $R_{noise} = 3.73\%$. (a-c) The amplitude and (d-f) the phase of ePIE, rPIE and sDR reconstructions, respectively. For this sparse data, ePIE and rPIE fail to converge in the reconstructions no matter how many iterations are used, but sDR converges to a high quality image. Furthermore, amplitudes and phases in the ePIE and rPIE reconstructions are mixed and could not be separated.

Download Full Size | PDF

3.2 Reconstruction from experimental data

3.2.1 Optical laser data

As an initial test of sDR with experimental data, we collect diffraction patterns from an USAF resolution pattern using a green laser with a wavelength of 543 nm. The incident illumination is created by a $150\, \mu m$ diameter pinhole. The pinhole is placed approximately 6 mm in front of the sample, creating a illumination wavefront on the sample plane that can be approximated by Fresnel propagation. The detector is positioned 26 cm downstream of the sample to collect far-field diffraction patterns. We raster scan across the sample with a step size of $50\, \mu m$ and acquire 169 diffraction patterns. We perform a sparsity test by randomly choosing 85 diffraction patterns (50% density) and run ePIE, rPIE, and sDR on this subset with 300 iterations. If we assume the probe diameter is to where the intensity falls to 10% of the maximum, then the overlaps are 73% and 46.4% for the full and sparsity sets respectively. The parameters used for the sDR test are $\beta _O = 0.9, \; \beta _P^0 = 0.5$ and $\sigma = 0.3$. Figure 5 shows that rPIE and sDR obtain a larger FOV than ePIE as both use regularization. Furthermore, sDR removes noise more effectively and obtains a flatter background than ePIE and rPIE. We monitor the R-factor (relative error) to quantify the reconstruction, defined as

(29)$$R_F = \frac{1}{N} \sum_{n=1}^N \| |\mathcal{F} ( P O_n ) | - \sqrt{I_n} \|_{1,1} \Big/ \| \sqrt{I_n} \|_{1,1}$$

$R_F$ is 16.94%, 13.95% and 13.28% for the ePIE, rPIE and sDR reconstructions, respectively.

Fig. 5. The ePIE (a), rPIE (b) and sDR (c) reconstructions respectively of a sparse data with 300 iterations, where sDR obtains a better quality reconstruction than ePIE and rPIE. Both sDR and rPIE produce a larger FOV than ePIE. Scale bar $200 \mu m$.

Download Full Size | PDF

3.2.2 Synchrotron radiation data

To demonstrate the applicability of sDR to more realistic experimental data, we reconstruct a ptychographic data set collected from the Advanced Light Source [16]. In this experiment, 710 eV soft x-rays are focused onto a sample using a zone plate and the far-field diffraction patterns are collected by a detector. A 2D scan, which imparts an estimated total dose of $1.17 \times 10^3$ Gy on 7,500 scan positions, spans approximately $10\times 4\, \mu m$. The sample is a portion of a HeLa cell labeled with nanoparticles, which is supported on a graphene-oxide layer. To compare the three algorithms, we choose a subdomain of a $3.70\times 3.70\, \mu m$ region, consisting of 2,450 diffraction patterns. Under the same assumption, the overlap is computed to be 79.5%.

First, we run the three algorithms with many iterations to assure the convergence of reconstructions (results do not change after such many iterations). Figure 6 displays the ePIE, rPIE, and sDR reconstructions after 10000 iterations, respectively. All three algorithms converge to images with high quality. When reducing the number of iterations to 100, we observe that sDR converges faster and reconstruct a larger FOV than both ePIE and rPIE. The individual nanoparticles, which serve as a resolution benchmark, are better resolved in the sDR reconstruction than reconstructions obtained by ePIE and rPIE. Furthermore, the reconstruction by ePIE at 100 iterations still contains artifacts such as faint square grids and halo-artifacts, which have been removed by rPIE and sDR.

Fig. 6. The ePIE (a), rPIE (b) and sDR (c) reconstructions of a $3.70 \times 3.70 \mu m$ region of the HeLa cell after 10000 iterations with $R_F = 15.80\%, \, 15.84\%$ and $14.40\%$, respectively. (d-f) Magnified regions ($1.66 \times 1.66 \mu m$) in (a-c), respectively. (g-l) The ePIE, rPIE, and sDR reconstructions and their magnified regions after 100 iterations with $R_F = 18.69\%, \, 16.60\%$ and $14.92\%$, respectively. sDR converges fastest among the three algorithms. Scale bar $500 nm$ and $200nm$ respectively.

Download Full Size | PDF

We next perform a sparsity test by randomly picking 980 out of 2,450 diffraction patterns, i.e. a reduction of data by 60%. The corresponding overlap of the sparsity set is 50.8%. Since fewer data is selected, the problem becomes less over-constrained. As a result, not only the quality but also the R-factor decreases. Figure 7 shows the reconstructions by ePIE, rPIE, and sDR with 10000 and 300 iterations. sDR is shown to perform better than ePIE and rPIE in the sense of quality and convergence rate. At 10000 iterations, the features are resolved better by sDR where more distinguishable individual nanoparticles are obtained. Especially, sDR converges faster than ePIE and rPIE which encounter slow convergence. At 300 iterations, sDR reconstructs a much better qualitative image with a larger FOV while ePIE and sDR cannot resolve the nanoparticles and still contain artifacts. The sDR result at 300 iterations looks very similar to the one at 10000 iterations. Additionally, the R-factor of sDR at 300 iterations is very close to the one at 10000 iterations while the R-factors of ePIE and rPIE are still very far from their final values. In both full and sparsity cases, $\beta _O=0.6, \; \beta _P^0=0.02$ and $\sigma = 0.7$ are used for the sDR tests.

Fig. 7. (a-c) The ePIE, rPIE, and sDR reconstructions of a sparse dataset and (d-f) their magnified regions after 10000 iterations with $R_F = 15.62\%, \,15.72\%$ and $13.89\%$, respectively. (g-l) The ePIE, rPIE and sDR reconstructions and their magnified regions after 300 iterations with $R_F = 19.52\%, \,16.58\%$ and $14.26\%$, respectively. To create the sparse data, we randomly pick 980 out of 2,450 diffraction patterns from the HeLa cell dataset. sDR reproduces a more distinguishable features than ePIE and rPIE. In particular, ePIE and rPIE encounter slow convergence while sDR converges faster with fewer iterations.

Download Full Size | PDF

To validate the quality of the sparsity results, we use Fourier Ring Correlation (FRC) [53] which measures the normalized cross-correlation coefficient between two 2-dimensional (or 3-dimensional) images over corresponding rings (or shells) in Fourier space (as a function of spatial frequency correlation). We assume that the ePIE reconstruction from the full dataset with 10000 iterations serves as the ground truth while all ePIE, rPIE and sDR reconstructions from 40% data set with 10000 iterations are the subjects. Figure 8 shows the FRC curves of ePIE, rPIE and sDR on a normalized frequency domain (% of Nyquist). The high correlation to high spatial frequency of all three algorithms confirms that their reconstructions are comparable to the full dataset result. FRC indicates that the reconstruction obtained by sDR is more highly correlated to the ground truth as compared to the other reconstructions. We now move to discussion section and future potential applications.

Fig. 8. Fourier ring correlation (FRC) of ePIE, rPIE, sDR reconstructions of the 40% dataset. FRC shows that the sDR reconstruction is more highly correlated to the ground truth(reconstruction using full dataset) than the ePIE and rPIE reconstructions.

Download Full Size | PDF

4. Discussion

The sDR test with these simulated and experimental data has been shown to perform well with fewer scan positions (low overlap ratio) covering the same FOV while regular methods such as ePIE and rPIE fail to converge or experience slow convergence. As the overlap is reduced, the reconstruction not only requires more iterations for convergence but the quality also decreases. This result is consistent with the “dose fractionation theorem” of Hegerl and Hoppe about the dose needed for a given resolution [54,55]. The high FRC results shown in Fig. 8 demonstrate that all algorithms can produce sustainable images with low overlap given enough time and iterations. However, the slow convergence makes ePIE and rPIE reconstructions infeasible in terms of computation time. Thanks to sDR, the reconstruction is accelerated and can be obtained with better quality in practical application.

Depending data and noise, the lower limit of overlap ratio with which sDR works can vary, but is roughly 40% in our experiments. For example in the synchrotron data, 60% of data reduction corresponding to 50.8% overlap, can still produce an acceptable reconstruction by sDR. This implies that sDR may have useful application where beam time is scarce and sample positioning is a significant fraction of the total acquisition time. Additionally, the effect of systematic errors such as instability in the illumination and sample drift, is reduced. sDR hence allows users to effectively collect more data per given beam time, which is important in this scared beam time case.

5. Conclusion

In this work, we have developed a fast and robust ptychographic algorithm, termed sDR, which incorporates two techniques relaxed Douglas-Rachford and semi-implicit scheme (semi-Backward Euler). Using both simulated and experimental data, we have demonstrated that sDR outperforms ePIE and rPIE in both quality and convergence speed. Especially in sparse data, sDR has been evidenced to have capability to defeat local minima, significantly reduce computation time, and obtain ptychographic reconstructions without losing much quality. We believe that sDR provides a powerful ptychography algorithm to the X-ray community. To promote the use of sDR algorithm, we publish our Matlab source code of sDR on our website www.physics.ucla.edu/research/imaging/sDR/index.html while the experimental data used in this paper can be found at zenodo website https://zenodo.org/record/3408463#.XX2A0ChKhdg.

Funding

National Science Foundation (DMR-1548924); U.S. Department of Energy (DE-SC0010378).

References

1. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400(6742), 342–344 (1999). [CrossRef]

2. J. Miao, T. Ishikawa, I. K. Robinson, and M. M. Murnane, “Beyond crystallography: Diffractive imaging using coherent x-ray light sources,” Science 348(6234), 530–535 (2015). [CrossRef]

3. K. J. Gaffney and H. N. Chapman, “Imaging atomic structure and dynamics with ultrafast x-ray scattering,” Science 316(5830), 1444–1448 (2007). [CrossRef]

4. I. Robinson and R. Harder, “Coherent x-ray diffraction imaging of strain at the nanoscale,” Nat. Mater. 8(4), 291–298 (2009). [CrossRef]

5. F. Pfeiffer, “X-ray ptychography,” Nat. Photonics 12(1), 9–17 (2018). [CrossRef]

6. J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15(6), 1662–1669 (1998). [CrossRef]

7. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef]

8. D. R. Luke, “Relaxed averaged alternating reflections for diffraction imaging,” Inverse Probl. 21(1), 37–50 (2005). [CrossRef]

9. J. A. Rodriguez, R. Xu, C.-C. Chen, Y. Zou, and J. Miao, “Oversampling smoothness: an effective algorithm for phase retrieval of noisy diffraction intensities,” J. Appl. Crystallogr. 46(2), 312–318 (2013). [CrossRef]

10. S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78(1), 011301 (2007). [CrossRef]

11. Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging: A contemporary overview,” IEEE Signal Process. Mag. 32(3), 87–109 (2015). [CrossRef]

12. J. M. Rodenburg, A. C. Hurst, A. G. Cullis, B. R. Dobson, F. Pfeiffer, O. Bunk, C. David, K. Jefimovs, and I. Johnson, “Hard-x-ray lensless imaging of extended objects,” Phys. Rev. Lett. 98(3), 034801 (2007). [CrossRef]

13. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-resolution scanning x-ray diffraction microscopy,” Science 321(5887), 379–382 (2008). [CrossRef]

14. M. Dierolf, A. Menzel, P. Thibault, P. Schneider, C. M. Kewish, R. Wepf, O. Bunk, and F. Pfeiffer, “Ptychographic x-ray computed tomography at the nanoscale,” Nature 467(7314), 436–439 (2010). [CrossRef]

15. M. Holler, M. Guizar-Sicairos, E. H. R. Tsai, R. Dinapoli, E. Müller, O. Bunk, J. Raabe, and G. Aeppli, “High-resolution non-destructive three-dimensional imaging of integrated circuits,” Nature 543(7645), 402–406 (2017). [CrossRef]

16. M. Gallagher-Jones, C. S. B. Dias, A. Pryor, K. Bouchmella, L. Zhao, Y. H. Lo, M. B. Cardoso, D. Shapiro, J. Rodriguez, and J. Miao, “Correlative cellular ptychography with functionalized nanoparticles at the fe l-edge,” Sci. Rep. 7(1), 4757 (2017). [CrossRef]

17. A. M. Maiden, M. C. Sarahan, M. D. Stagg, S. M. Schramm, and M. J. Humphry, “Quantitative electron phase imaging with high sensitivity and an unlimited field of view,” Sci. Rep. 5(1), 14690 (2015). [CrossRef]

18. D. A. Shapiro, Y.-S. Yu, T. Tyliszczak, J. Cabana, R. Celestre, W. Chao, K. Kaznatcheev, A. L. D. Kilcoyne, F. Maia, S. Marchesini, Y. S. Meng, T. Warwick, L. L. Yang, and H. A. Padmore, “Chemical composition mapping with nanometre resolution by soft x-ray microscopy,” Nat. Photonics 8(10), 765–769 (2014). [CrossRef]

19. D. F. Gardner, M. Tanksalvala, E. R. Shanblatt, X. Zhang, B. R. Galloway, C. L. Porter, R. Karl Jr, C. Bevis, D. E. Adams, H. C. Kapteyn, M. M. Murnane, and G. F. Mancini, “Subwavelength coherent imaging of periodic samples using a 13.5 nm tabletop high-harmonic light source,” Nat. Photonics 11(4), 259–263 (2017). [CrossRef]

20. J. Marrison, L. Räty, P. Marriott, and P. O’Toole, “Ptychography – a label free, high-contrast imaging technique for live cells using quantitative phase information,” Sci. Rep. 3(1), 2369 (2013). [CrossRef]

21. J. Deng, Y. H. Lo, M. Gallagher-Jones, S. Chen, A. Pryor, Q. Jin, Y. P. Hong, Y. S. G. Nashed, S. Vogt, J. Miao, and C. Jacobsen, “Correlative 3d x-ray fluorescence and ptychographic tomography of frozen-hydrated green algae,” Sci. Adv. 4(11), eaau4548 (2018). [CrossRef]

22. S. Gao, P. Wang, F. Zhang, G. T. Martinez, P. D. Nellist, X. Pan, and A. I. Kirkland, “Electron ptychographic microscopy for three-dimensional imaging,” Nat. Commun. 8(1), 163 (2017). [CrossRef]

23. Y. Jiang, Z. C. Chen, Y. Han, P. Deb, H. Gao, S. Xie, P. Purohit, M. W. Tate, J. Park, S. M. Gruner, V. Elser, and D. A. Muller, “Electron ptychography of 2d materials to deep sub-ångström resolution,” Nature 559(7714), 343–349 (2018). [CrossRef]

24. Y. H. Lo, L. Zhao, M. Gallagher-Jones, A. Rana, J. J. Lodico, W. Xiao, B. C. Regan, and J. Miao, “In situ coherent diffractive imaging,” Nat. Commun. 9(1), 1826 (2018). [CrossRef]

25. R. Hesse, D. R. Luke, S. Sabach, and M. K. Tam, “Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging,” SIAM J. Imaging Sci. 8(1), 426–457 (2015). [CrossRef]

26. A. J. D’Alfonso, A. J. Morgan, A. W. C. Yan, P. Wang, H. Sawada, A. I. Kirkland, and L. J. Allen, “Deterministic electron ptychography at atomic resolution,” Phys. Rev. B 89(6), 064101 (2014). [CrossRef]

27. L. Bian, J. Suo, G. Zheng, K. Guo, F. Chen, and Q. Dai, “Fourier ptychographic reconstruction using wirtinger flow optimization,” Opt. Express 23(4), 4856–4866 (2015). [CrossRef]

28. R. Horstmeyer, R. Y. Chen, X. Ou, B. Ames, J. A. Tropp, and C. Yang, “Solving ptychography with a convex relaxation,” New J. Phys. 17(5), 053044 (2015). [CrossRef]

29. P. Thibault and A. Menzel, “Reconstructing state mixtures from diffraction measurements,” Nature 494(7435), 68–71 (2013). [CrossRef]

30. M. Guizar-Sicairos and J. R. Fienup, “Phase retrieval with transverse translation diversity: a nonlinear optimization approach,” Opt. Express 16(10), 7264–7278 (2008). [CrossRef]

31. A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109(10), 1256–1262 (2009). [CrossRef]

32. A. M. Maiden, M. J. Humphry, M. C. Sarahan, B. Kraus, and J. M. Rodenburg, “An annealing algorithm to correct positioning errors in ptychography,” Ultramicroscopy 120, 64–72 (2012). [CrossRef]

33. F. Zhang, I. Peterson, J. Vila-Comamala, A. Diaz, F. Berenguer, R. Bean, B. Chen, A. Menzel, I. K. Robinson, and J. M. Rodenburg, “Translation position determination in ptychographic coherent diffraction imaging,” Opt. Express 21(11), 13592–13606 (2013). [CrossRef]

34. A. Tripathi, I. McNulty, and O. G. Shpyrko, “Ptychographic overlap constraint errors and the limits of their numerical recovery using conjugate gradient descent methods,” Opt. Express 22(2), 1452–1466 (2014). [CrossRef]

35. P. Thibault and M. Guizar-Sicairos, “Maximum-likelihood refinement for coherent diffractive imaging,” New J. Phys. 14(6), 063004 (2012). [CrossRef]

36. P. Godard, M. Allain, V. Chamard, and J. Rodenburg, “Noise models for low counting rate coherent diffraction imaging,” Opt. Express 20(23), 25914–25934 (2012). [CrossRef]

37. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29(8), 1606–1614 (2012). [CrossRef]

38. E. H. R. Tsai, I. Usov, A. Diaz, A. Menzel, and M. Guizar-Sicairos, “X-ray ptychography with extended depth of field,” Opt. Express 24(25), 29089–29108 (2016). [CrossRef]

39. S. Marchesini, H. Krishnan, B. J. Daurer, D. A. Shapiro, T. Perciano, J. A. Sethian, and F. R. N. C. Maia, “SHARP: a distributed GPU-based ptychographic solver,” J. Appl. Crystallogr. 49(4), 1245–1252 (2016). [CrossRef]

40. A. Maiden, D. Johnson, and P. Li, “Further improvements to the ptychographical iterative engine,” Optica 4(7), 736–745 (2017). [CrossRef]

41. M. Rose, T. Senkbeil, A. R. von Gundlach, S. Stuhr, C. Rumancev, D. Dzhigaev, I. Besedin, P. Skopintsev, L. Loetgering, J. Viefhaus, A. Rosenhahn, and I. A. Vartanyants, “Quantitative ptychographic bio-imaging in the water window,” Opt. Express 26(2), 1237–1254 (2018). [CrossRef]

42. D. P. Bertsekas, “Incremental gradient, subgradient, and proximal methods for convex optimization: A survey,” in Optimization for Machine Learning, vol. 3S. Sra, S. Nowozin, and S. J. Wright, eds. (MIT press, 2012).

43. N. Parikh and S. Boyd, “Proximal algorithms,” FNT in Optimization 1(3), 127–239 (2014). [CrossRef]

44. Z. Wen, C. Yang, X. Liu, and S. Marchesini, “Alternating direction methods for classical and ptychographic phase retrieval,” Inverse Probl. 28(11), 115010 (2012). [CrossRef]

45. J. Douglas and H. H. Rachford, “On the numerical solution of heat conduction problems in two and three space variables,” Trans. Amer. Math. Soc. 82(2), 421 (1956). [CrossRef]

46. P. Tseng, “Applications of a splitting algorithm to decomposition in convex programming and variational inequalities,” SIAM J. Control Optim. 29(1), 119–138 (1991). [CrossRef]

47. P.-L. Lions and B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM J. Numer. Anal. 16(6), 964–979 (1979). [CrossRef]

48. J. Eckstein and D. P. Bertsekas, “On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators,” Math. Program. 55(1-3), 293–318 (1992). [CrossRef]

49. Y. Nesterov, Introductory lectures on convex optimization: A basic course, vol. 137 (Springer, 2018).

50. P. Tseng, “On accelerated proximal gradient methods for convex-concave optimization,” submitted to SIAM J. on Optim.2, 3 (2008).

51. S. Bubeck, “Convex optimization: Algorithms and complexity,” FNT in Machine Learning 8(3-4), 231–357 (2015). [CrossRef]

52. C. Zuo, J. Sun, and Q. Chen, “Adaptive step-size strategy for noise-robust fourier ptychographic microscopy,” Opt. Express 24(18), 20724–20744 (2016). [CrossRef]

53. M. Van Heel and M. Schatz, “Fourier shell correlation threshold criteria,” J. Struct. Biol. 151(3), 250–262 (2005). [CrossRef]

54. R. Hegerl and W. Hoppe, “Influence of electron noise on three-dimensional image reconstruction,” Zeitschrift für Naturforschung A 31(12), 1717–1721 (1976). [CrossRef]

55. M. Howells, T. Beetz, H. Chapman, C. Cui, J. Holton, C. Jacobsen, J. Kirz, E. Lima, S. Marchesini, H. Miao, D. Sayre, D. Shapiro, J. Spence, and D. Starodub, “An assessment of the resolution limitation due to radiation-damage in x-ray diffraction microscopy,” J. Electron Spectrosc. Relat. Phenom. 170(1-3), 4–12 (2009). [CrossRef]

Semi-implicit relaxed Douglas-Rachford algorithm (sDR) for ptychography

Abstract

1. Introduction

2. Proposed algorithm

2.1 Semi-implicit algorithm

2.2 Relaxed Douglas-Rachford algorithm

2.3 Semi-implicit relaxed Douglas-Rachford algorithm (sDR)

3. Experimental results

3.1 Reconstruction from simulated data

3.2 Reconstruction from experimental data

3.2.1 Optical laser data

3.2.2 Synchrotron radiation data

4. Discussion

5. Conclusion

Funding

References

Cited By

Figures (8)

Equations (29)

Optics Express