Advanced Denoising for X-ray Ptychography

The success of ptychographic imaging experiments strongly depends on achieving high signal-to-noise ratio. This is particularly important in nanoscale imaging experiments when diffraction signals are very weak and the experiments are accompanied by significant parasitic scattering (background), outliers or correlated noise sources. It is also critical when rare events such as cosmic rays, or bad frames caused by electronic glitches or shutter timing malfunction take place. In this paper, we propose a novel iterative algorithm with rigorous analysis that exploits the direct forward model for parasitic noise and sample smoothness to achieve a thorough characterization and removal of structured and random noise. We present a formal description of the proposed algorithm and prove its convergence under mild conditions. Numerical experiments from simulations and real data (both soft and hard X-ray beamlines) demonstrate that the proposed algorithms produce better results when compared to state-of-the-art methods.


Introduction
Ptychography [1][2][3][4] has become an increasingly popular imaging technique, and it is used nowadays in scientific fields as diverse as condensed matter physics [5], cell biology [6], materials science [7,8], and electronics [4,9], among other areas. Compared to standard lens-based microscopy, the resolution in ptychography is not limited by the size of the illumination probe, but by the wavelength and number of photons used in an experiment. Unfortunately, multiple experimental challenges have to be tackled to achieve a high-quality reconstruction from a ptychographic experiment in practice. Complex experimental systems require nanometer stability while collecting large amounts of diffraction data, which can require several hours depending on the experimental setup. Data can be then contaminated by structured parasitic scattering [10][11][12], detector read-out noise, Poisson noise derived from the photon counting, and different types of outliers. Additionally, the characterization of the illumination source is also commonly considered as a joint problem to the object reconstruction, known as Blind Ptychography (BP) [3,[13][14][15].
During the last decades, researchers have developed several schemes to solve the BP problem. Arguably, the most popular ones are extended Ptychographic Iterative Engine (ePIE) [3], Difference Map [13,16], Maximum Likelihood (ML) method [17], Proximal Splitting algorithm [14], Relaxed Averaged Alternating Reflections (RAAR [18]) based algorithms [19], and generalized Alternating Direction Method of Multipliers (ADMM) [20][21][22]) based BP [15]. Although some implementations of those methods present ad hoc solutions for structural noise removal, there is no formal analysis or in-depth characterization of such experimental problems, even though they can be critical to achieve robust high-resolution images, specially from weakly scattering or low contrast specimens.
Experimental noise overview Direct reconstructions based on raw experimental data often contain visible artifacts. They commonly derive from outliers and structured and randomly distributed uncorrelated noise sources. Below there is a characterization of the main causes of such anomalies.
Structured noise (1) Parasitic scattering: Derives from the scattering from any element along the beam path other than the sample and the optical elements desired harmonic order e.g.: slits, pinholes, lens imperfections, harmonic contamination, air, gas, etc. We refer to this as the background.
(2) Saturation: Occurs when the flux exceeds the detector capacity, the range of the Analog to Digital Converter (ADC) or the maximum photon counting rate.
(3) Dark noise: Produced by the detector dark current, light from position encoders and interferometers, or light from the environment. It differs from the parasitic scattering as it occurs also when the x-ray beam is off and can be measured in advance by turning off the beam and (commonly referred to as dark frame). The constant component can be subtracted from the measured signal before applying a reconstruction solver or alternatively it can be incorporated in the generalized background.

Random noise
(4) Photon-counting noise: Caused by the quantization of the wavefront at the detector.
(5) Read-out and detector noise: Derives from electronic interference when the ADC converts the charge distribution to a digital signal and the random component of the dark current. The interference and parallel ADC read-out electronics can produce ringing and correlation in the digital signal.
Outliers (6) Bad frames: Glitches on the functioning of mechanical components inside an instrument such as shutter timing can cause some measured frames to be corrupted.
(7) Cosmic rays: High energy cosmic particles that penetrate the atmosphere, the building and the instrument enclosure and hit the detector. Although they are rare, those types of signals normally present very high charge.
(8) Bad pixels: Result from fabrication errors may cause dead pixels on a detector, these can be typically measured in advanced and masked out during the reconstruction process.
In this paper we focus on a solution to address (1), (6) and (7), regarding structured noise and outliers, and also propose a technique to deal with a combination of different random noise distributions coming from (4) and (5). The noise from (3) and (8) can be measured in advance, and it is commonly subtracted. Saturation noise (2) can be equally identified and addressed by masking out the saturated signal areas and it is outside the focus of this paper. Random noise sources (4) and (5) have been considered in the literature [15,17,23,24], but there is no solution that can address the combination of two or more random noise distributions. In this paper, we propose a mixed-noise term for ptychography to effectively deal with any combination of Poisson and Gaussian noise sources.
Some of the previous pathologies can be identified inspecting the measured data. An illustrative example of parasitic noise (1) is reported in Fig. 1. A single measured diffraction pattern ( Fig.  1(a)) presents a distinct reverse "L-Shape" artifact, which consistently appears on all other diffraction measurements from the same dataset. If parasitic scattering is not considered, a standard ptychography reconstruction produces very noisy and low contrast images (Figs. 1(d) and 1(e)). The phase part is lost in some areas, concealed by the biased background of the recovered phase, and the resulted intensities ( Fig. 1(b)) contain a slight alteration of the parasitic noise from the measured data. (A more detailed comparison between different methods, with and without structured noise correction, is provided in section 3) Reconstruction results using RAAR (implemented in SHARP [19]) and the proposed algorithm. The dataset corresponds to one of the 3D ptychography experiments from the results presented in [8]. (a) Measured intensities from a single diffraction pattern. Recovered intensities (b), object amplitude (amp.) (d) and phase (e), respectively, using RAAR without structured noise correction. Recovered intensities (c), object amplitude (f) and phase (g), respectively, using the proposed algorithm (alg.). The figures in the first row are shown in scale (·) 0.05 following the colorbar shown on the top right corner; the figures in the second row are shown in linear scale following the colorbars on the left corner (amplitude) and right corner (phase).
Related work A beam stop scheme was employed in [11] to reduce the parasitic noise using intensities from two raster scans (i.e. scanning on a Cartesian 2D grid) with and/or without beamstop. However, the artifacts caused by the structured noise are still noticeable in Fig.  2 of [11], and some fine features still appear undefined. A multiple-mode approach in [25] was proposed to reduce the background. In [12], the data was pre-processed by subtracting an adjustable scalar value multiplied by the dark frame. A method to remove the parasitic background automatically employed a gradient descent method with adaptive stepsizes was proposed in [26] and implemented using a distributed multi-GPU implementation [19] in a real-time streaming production code at a user facility [27]. Several researchers have also studied techniques to remove the random noise from photoncounting or read-out noises by employing a sparse regularization technique to further improve the image quality for conventional ptychography, .e.g. Tikhonov regularization with nonlinear conjugate gradient method [17], total variation regularization with ADMM [28] and dictionary learning method with proximal algorithm [29]. Specially, for the more practical noises, like a mix of Poisson and Gaussian noises, several variational methods have been proposed for linear inverse problems, e.g. generalized Anscombe transform [30], total variation based Shift-Poisson method [31] and weighted least squares method [32]. However, such methods have not been applied to the conventional ptychography yet, which is essentially a nonlinear inverse problem. Beyond algorithmic methods, there are also several contributions from beamline scientists attenuating the background problem via hardware solutions [10,33].

Motivations and contributions
In a ptychography experimental setup it is impossible to isolate parasitic noise from the measured data. Furthermore, the data is often contaminated by outliers and an additional mix of random noise types. In this context, an automatic denoising algorithm with convergence guarantee to jointly remove both structured and random noise can be very valuable. In this paper, we propose ADMM Denoising for Phase retrieval (ADP), a new algorithm that addresses all the main sources of experimental noise in a blind X-ray ptychography experimental setting. The proposed algorithm achieves similar convergence speed as state-of-the-art blind ptychography methods and it is constructed on the framework presented in [15]. Even if additional hardware or pre-processing techniques are in place, ADP can be used to further enhance the reconstruction quality of the retrieved object. The main contributions of this work are listed below: i. A new algorithm for automatic structured and random mixed noise removal with outliers correction. The algorithm is based on the forward physical model of the parasitic scattering noise (additive frame-invariant background) and on the maximum a posteriori (MAP) estimation combined with the shift-Poisson method [31] for mixed types of noise. By assuming the piecewise smoothness of sample patches, We propose an additional framewise regularization term to further enhance the denoising process and improve the image quality, instead of using a simple sparsity term of the sample [28]. To the best of our knowledge, it is the first time an iterative algorithm with rigorous analysis for background removal in ptychographic imaging is proposed.
ii. By formulating the parasitic noise variable with positive constraint as a quadratic term, fast ADMM algorithms are designed for the proposed model, with and without regularization. The convergence guarantee of the method is proved under very mild conditions. Computationally, each subproblem can be solved very efficiently, using pointwise operations or simple linear algebra solvers, requiring few arithmetic operations per iteration and exposing high parallelism.
iii. Multiple numerical experiments on both simulated and experimental data from different light sources demonstrate the enhanced quality of the proposed algorithm. Reconstruction results using the proposed algorithm can be found in Fig. 1 and in the experimental results section.

Proposed iterative algorithm
Mathematical formula for standard ptychography In a standard ptychography experiment, a localized coherent X-ray probe (or illumination) ω scans through an image (or sample) u, while the detector collects a sequence of phaseless intensities in the far field. Throughout this paper, we consider the following discrete setting: The variable u ∈ C n (C n is the complex Euclidean space) corresponds to a 2D sample with √ n × √ n pixels, and ω ∈ Cm is a localized 2D probe with √m × √m pixels, both u and ω written as a vector by a lexicographical order. A stack of phaseless measurementsĨ j ∈ Rm + ∀0 ≤ j ≤ J − 1 (Rm + is the Euclidean space with non-negative elements) is collected withĨ j = |F (ω • S j u)| 2 , where | · | represents the element-wise absolute value of a vector, • denotes the element-wise multiplication, and F denotes the normalized discrete Fourier transformation. S j ∈ Rm ×n is a binary matrix that specifies a small window with the index j and sizem over the entire sample u. The BP problem can then be expressed as follows: To find ω ∈ Cm and u ∈ C n , such that |A(ω, u) where bilinear operators A : Cm × C n → C m and A j :

Proposed model
Experimentally, not all of the parasitic scattering (background) is the same on each frame (here, "frame" means one image from the CCD detector). However, the frame-invariant component is one of the key sources of the artifacts in standard reconstruction results due to the redundancy of ptychography [11,12]. Essentially, it is an additive components of signals, apart from those originating from the sample and the beam hitting it. Hence, we start by assuming that the background of each frame is approximately unchanged [26], i.e., it is frame-invariant. It is important to note that parasitic noise can also be interpreted as an independent structured noise source per frame. The frame-invariant approximation is critical for proper modeling and the problem is actually not well defined if this approximation is not considered. This is due to the fact that the independent noise per frame alternative has a trivial solution where the background is equal to the residual between the measured and reconstructed intensities from the iterative solution of the sample and probe. The frame-invariant approximation is, besides simple, very efficient at removing structural noise to produce high-contrast recovery images. Formally, the true intensity dataĨ is contaminated by non-negative parasitic structured noisê φ ∈ Rm + , such that the measured dataÎ = (Î T 0 ,Î T 1 , · · · ,Î T J−1 ) T ∈ R m + withÎ j ∈ Rm + , ∀0 ≤ j ≤ J − 1 is recorded for each frame as follows: In practice, the data is contaminated by noise aŝ where Noi denotes the degrading operator caused by different noise sources. Without loss of generality, we consider the Poisson and Gaussian mixed noise model with parasitic noise as: where i.i.d. is short for "Independent and Identically Distributed", n j denotes the white Gaussian noise with zero mean and variance σ > 0. It is not difficult to verify that var(Î j (t) + σ 2 ) = var(Î j (t)) = |A j (ω, u)(t)| 2 +φ + σ 2 , and mean(Î j (t) + σ 2 ) = |A j (ω, u)(t)| 2 +φ + σ 2 . Hence, I j (t) :=Î j (t) + σ 2 has the same variance and mean. By further estimating I j by a Poisson distribution with the shift-Poisson technique [31] and using KL divergence derived by maximum likelihood estimation of Poisson noise, the following optimization problem is derived: where 1m is a vector of all ones, φ :=φ + σ 2 , ·, · denotes the inner product as and positivity constraint set is defined as E := φ ∈ Rm + : φ ≥ 0 . Here, the indicator function I E is defined as: However, when this metric is used to measure the distance between the collected and recovered intensities, the corresponding first-order iterative algorithm is significantly slow [15]. This is because the residual slowly decreases as the iterations go, and the subproblem with respect to the auxiliary variable is related with quartic equations, which are solved with high computational cost [34].
In order to deal with the misfit caused by experimental outliers, a weight vector C := To further handle the noise or other artifacts, we consider the total variation regularization on the split patches of the sample, instead of the regularization of the entire sample in [28]. As a result, we propose the following Kullback-Leibler (KL) divergence-based problem with framewise sparsity and parasitic noise retrieval: with a penalization parameter ε > 0, where TV denotes the total variation and λ is a positive parameter to balance the regularization and data fitting terms.

ADMM denoising for phase retrieval (ADP)
Generally speaking, the ADMM is a variant of the augmented Lagrangian method [20,21] that uses partial updates for the dual variable, and each subproblem can be easily solved compared with the original optimization problem. Due to its scalability and flexibility, it has been a popular algorithm for large-scale optimization problems arising in computer vision, statistics, machine learning, and other related areas, and see more details in the review paper [22] and references therein. Since the proposed model is non-convex and non-smooth, ADMM will be adopted to solve it, which allows for bigger stepsizes by avoiding directly calculating the gradient of the objective functional and therefore has fast convergence with low computation cost per iteration.
We formulate below the proposed ADP algorithm to solve Eq. (7). First, we consider the problem without regularization by setting λ = 0. If directly following [35,36], the subproblem with respect to the variable φ does not have a closed form solution: A gradient descent or proximal type algorithm reported in [15] to directly solve this problem presents very slow convergence. In order to develop a more efficient algorithm for each subproblem, we first rewrite Eq. (7) with λ = 0 as min ω,u,φ whereφ := √ φ. Note that the constraint for φ is automatically removed in Eq. (9). Also, the variables A j (ω, u) andφ play now the same role in the objective function and adding auxiliary variables can derive a more efficient algorithm than directly solving the previous problem.
Introducing the auxiliary variables z j = A j (ω, u) and µ j =φ ∀0 ≤ j ≤ J − 1, produces the following equivalent constraint optimization problem: with The corresponding augmented Lagrangian for Eq. (10) reads: where denotes the real part of the complex-valued number. The proposed ADP algorithm is formulated as follows: with The solution to each subproblems above are reported in Algorithm 1, and further details can be found in the Appendix B. When setting µ 0 j =φ 0 = 0, and Λ 0 2 = 0, following the above algorithm we have that µ k j = 0 ∀k = 1, 2, · · · , since µ = 0 is a stationary point of the proposed model Eq. (7). Therefore, we use the following warm-start scheme to initialize µ: first solve a problem without structured noise correction by setting µ 0 j = 0, and after J 0 iterations reset µ k j using µ k+1 j := max{0, 1 J j (I j − |A(ω k+1 , u k+1 )| 2 )}. The weight function {C j } can be simply fixed to be one. It can also be dynamically changed if the measured data contains outliers, such that the ADP algorithm is slightly modified with additional update of the weight function between Step 3 and Step 4, and please see more details for update scheme of this function in Remark 2.2.
Step 5. If k = J 0 : reinitialize the average of µ k+1 by Else if (satisfying some stopping condition): output u k+1 as the final reconstructed result Else: set k := k + 1, and goto Step 1

Convergence analysis
We assess the convergence of ADP for the proposed model with λ = 0 in the following theorem: Theorem 1. By denoting Y k := (ω k , u k , z k , µ k ,φ k , Λ k 1 , Λ k 2 ), any limit point of {Y k } produced by ADP is the stationary point of Eq. (7) with λ = 0 if the stepsize r is sufficiently large and ε > 0.
The proof of the above theorem can be found in Appendix C. The referred proof demonstrates that any limit point of the iterative sequence is the stationary point of the proposed model. In order to prove the convergence of the whole sequence we need more constraints for the object and illumination, and also a careful selection scheme of the proximal parameters α 1 , α 2 [15]. The introduction of ε is also needed for the convergence analysis. Numerically, such penalization does not produce any obvious improvement on either convergence speed or reconstruction quality. Hence, for simplicity, we only show the performance of proposed algorithm with ε = 0.

ADP with regularization (ADPr)
In this subsection we consider the proposed model Eq. (7) with λ > 0. Introducing the auxiliary variables p j = ∇S j u, z j = A j (ω, u) and µ j =φ ∀0 ≤ j ≤ J − 1 produces the following equivalent constraint optimization problem: The corresponding augmented Lagrangian reads We propose the following generalization of ADP to solve the problem above, referred to as ADP with regularization (ADPr): Details of solvers for these subproblems are reported in Appendix B. We only summarize the overall algorithm in Algorithm 2 below.
In a similar manner as with ADP, one can derive the convergence of ADPr. We provide the following theorem omitting the proof details in this case: Eq. (7) if the stepsizes r, β are sufficiently large and ε > 0.

Remark 2.2.
In order to handle outliers, we propose the following adaptive weights {C j }. They are designed to produce bigger values when the residual between the measured data and recovered intensities are smaller. The weight function C j is given based on the inverse of the residual as follows: with γ ≤ 2. The weight function can be updated between Step 3 and Step 4 in ADP algorithm, or Step 4 and Step 5 in ADPr algorithm.

Remark 2.3.
Raster grid scanning can cause visible periodical artifacts, and therefore a constraint of the lens can be enforced as [19]. Furthermore, a detector mask can be used as well.

Experimental results
The experimental results of this section are generated using both simulated and experimental ptychography data. In most experiments, the performance of the proposed algorithm is compared with that of SHARP (using RAAR) [19], a ptychography software solution that implements state-of-the-art blind ptychography algorithms and background retrieval techniques [26] (further details can be found in Appendix A). In some test we will evaluate SHARP without and with background retrieval, the later will be referred to as SHARP-B. All experimental results in this section employ raster scans. Regarding the parameter initialization of the proposed method, we set J 0 = 5 (after J 0 iterations reinitialize the background) for both ADP and ADPr. We use 5 iterations for Step 2 of ADPr. We introduce two different criteria to evaluate performance. For simulation data, the ground truth can be used to compute the signal-to-noise ratio (SNR) of u k as: where u g corresponds to the ground truth image. The error is computed up to the translation T * , and the phase shift and scaling factor ζ * are determined by: For experiment data, the R-factor for ADP and ADPr is used, defined as: The R-factor of SHARP is defined by settingφ = 0, whereas the R-factor of SHARP-B is defined in the same way as with ADP and ADPr.

Synthetic data
Given a true illumination and sample ω and u , respectively, and φ as the parasitic noise component, the simulated intensities are generated as follows: where n j denotes the white Gaussian noise (the variance is set to 1.0 × 10 −6 × 1 m t I poi (t), with I poi (t) := Poisson |A j (ω , u )| 2 (t) + φ (t) ). We also consider a further corruption of the intensities by outliers. They are simulated as patches of bad measurements appearing in 10% of the frames, with size ∼ 50% of the frame size. The corrupted intensity values are set to 0. We use raster grid scan with stepsizes of 16 pixels, so an additional constraint of the illumination (a.k.a. support of the lens or illumination Fourier mask) is used to prevent potential periodical artifacts. In the experiments below, the frame size employed is 64 × 64 pixels and all algorithms stop after 1000 iterations. Fig. 2 presents reconstruction results employing two different simulation datasets. The first experiment (second and fourth rows) uses data contaminated only by parasitic noise (φ ). The second experiment (third and fifth rows) reconstructs data with a mix of white Gaussian and Poisson noises as in Eq. (32) and the outliers described above as well. The ground truth images for the sample, illumination and background are depicted in Figs. 2(c), 2(k), Figs. 2(a)-2(b), and Fig. 3(a), respectively. For each of the previous reconstruction experiments, we report the retrieved background in Fig. 3 for SHARP-B, ADP and ADPr. When no mix of noises and outliers are considered, ADP produces much higher quality reconstructions (Figs. 2(f) and 2(n)) than SHARP and SHARP-B (Figs. 2(d)-2(e) and Figs. 2(l)-2(m)). The retrieved background, for both SHARP-B and ADP (Figs. 3(b)-3(c)), is consistent with those results. When additional outliers are introduced, ADP still generates better result with a dramatic reduction in outliers-related artifacts (Figs. 2(i) and 2(q)). In this case, ADPr (Figs. 2(j) and 2(r)) helps further removing the noise from the baseline ADP reconstruction. For the second experiment, the proposed algorithms also retrieve the background (Figs. 3(e)-3(f)) more precisely than SHARP-B (Fig. 3(d)). The SNRs of recovery results are reported: SNR = 13.2, 15.7, 76.7 by SHARP, SHARP-B, and ADP for the case of single parasitic noise; SNR = 9.8, 10.2, 18.5, 27.9 by SHARP, SHARP-B, ADP and ADPr for the case with additional mixture noises and outliers. Again, one can readily see the advantages of proposed ADP and ADPr.

Soft X-ray dataset from the ALS, Lawrence Berkeley Lab
The data used in the following experiment was published in [8]. It corresponds to a 3D ptychography imaging experiment of battery cells at 708 eV, obtained from different projection angles.
The recovered amplitude and phase are reported in Fig. 4 and zoom-in view in Fig. 5 when using SHARP, SHARP-B, ADP and ADPr. The recovered intensities of a single frame and the recovered background for the same experiment are shown in Fig. 6. When no background retrieval is used, the recovery results from SHARP (Figs. 4(a) and 4(e), and zoom-in view in Figs. 5(a) and 5(e) ) are significantly noisy, with some areas presenting vague or inappreciable features. SHARP-B was specially designed for structured noise from ALS soft X-ray sources and the results in Figs The R-factors of the previous experiment are 0.2399, 0.1112, 0.0957, and 0.0969, when using SHARP, SHARP-B, ADP and ADPr, respectively. This demonstrates how the proposed algorithm achieves smaller residuals thus producing higher accuracy results. We also provide cutline values results from the previous experiment in Fig. 8, selecting a line from the center of the reconstructed image. It can be seen from the figure how the proposed algorithm generates a higher contrast reconstruction than SHARP-B, specially in the retrieved phase. Convergence results are presented in Fig. 7, showing the R-factor as the iterations go. Those results illustrate how the R-factors of the proposed algorithm steadily decrease, providing additional evidence on the stability of the method.
Hard X-ray dataset from PETRA III at DESY We also report experimental results using a dataset measured at beamline P06 at PETRA III, DESY [11] at a photon energy of 11919 eV. This dataset consists of 2 ptychographic measurements, one without and one with beamstop. The following results are generated using (1) the data without beamstop (noBS), (2) the data with beamstop (BS), and (3) the merged data from BS and noBS (the low frequencies from noBS plus the high frequency from BS). The datasets are reconstructed using ePIE (extended Ptychographic Iterative Engine [3]), which corresponds to the results presented in [11], and with the proposed ADP algorithm. For a fair comparison, all results are produced without position retrieval. The amount of noise of this dataset is significantly higher than the one reported in the   previous experiments: in this case, the phase contrast produced by the hard X-ray beam is much weaker than in the soft X-ray dataset. Because of this, more iterations are required to achieve reasonable reconstructions. The following results use at most 2000 iterations (we use early-stop for ePIE since it will blow up finally) for both algorithms.
The reconstruction results are shown in Fig. 9. The reconstructed phase parts using ePIE is significantly noisy, with lots of blurred out features ( Fig. 9 (a) and (b)). Visually, ADP generates a much cleaner reconstructed phase parts with much better defined features (especially noticeable in the features around the color circles). ADPr in this case produces similar quality results as ADP with no regularization.
We also report an additional reconstruction result using only the BS dataset with ADP ( Fig. 9  (e)). It is important to note that for this dataset low-frequency information is almost completely lost. We can see how very sharp features are well recovered in the reconstructed phase, while producing a very clean background. This illustrates the robustness of the proposed algorithm even when the measured data is incomplete or contaminated by heavy noise. Some features from the BS experiment with ADP seem to be enhanced, compared with the merged data (specially noticeable again in the color circle areas). We remark that to make our algorithm work well, a good initialization by the illumination from the non-beamstop dataset is adopted. As a future work, we will explore more efficient strategy for the initialization.
Since our proposed algorithms are implemented by python, which is not specially optimized (a) ePIE (noBS) Fig. 9. Hard X-ray experimental results from the data presented in [11]. First row: recovered phase using ePIE, with the noBS dataset (a), and with the merged dataset (b). Second row: recovered phase using ADP, with noBS data (c), merged data ( compared with SHARP (implemented by highly optimized CUDA code). Therefore, we only give the estimate of computational complexity of proposed ADP as about 63m + 2J × FFT per iteration, while about 78m + 2J × FFT for SHARP-B, that demonstrates the proposed ADP needs a bit less computational cost than SHARP-B, where FFT means the computational complexity of 2D fast Fourier Transform for the matrix with size of √m × √m (same size of the illumination).

Conclusion
This paper proposes a novel algorithm for (blind) ptychography reconstruction with integral experimental denoising. All main experimental sources of noise are addressed and characterized in the proposed solution as a mixture of (1) structured parasitic noise, (2) a mix of random noise sources (both Gaussian or Poison), and (3) data outliers. The algorithm is formulated to be robust and efficient, requiring few arithmetic operations and being inherently parallel in most of the steps. The convergence of the method is also proved under mild conditions. Experimental results analyze a variety of datasets with different experimental noise sources and demonstrate how the proposed algorithm achieves superior reconstruction results and denoising than state-of-the-art solutions. As a future work, we will consider partial coherence [37] and position retrieval [26] in the model, and also explore additional sparsity techniques [38], deep-learning [39], and dictionary-learning methods [29] to further enhance the reconstruction results. and lim k→+∞ j S T j ((ω n k ) * • F * (z n k −1 j + Λ n k −1 1, j )) = j S T j ((ω ) * • F * (z j + Λ 1, j )).