Phase retrieval for object and probe using a series of defocus near-ﬁeld images

: Full ﬁeld x-ray propagation imaging can be severely de-teriorated by wave front aberrations. Here we present an extension of ptychographic phase retrieval with simultaneous probe and object reconstruction suitable for the near-ﬁeld diffractive imaging setting. Update equations used to iteratively solve the phase problem from a set of near-ﬁeld images in view of reconstruction both object and probe are derived. The algorithm is tested based on numerical simulations including photon shot noise. The results indicate that the approach provides an efﬁcient way to overcome restrictive idealizations of the illumination wave in the near-ﬁeld (propagation) imaging.


Introduction
Lens-less x-ray coherent diffractive imaging (CDI) has provided a manifold of possible approaches to solving the classical phase problem in x-ray diffraction, with the goal of nanoscale imaging beyond the limits imposed by lens fabrication [1,2]. In its first demonstrations [3], based on the Error Reduction or Hybrid-Input-Output algorithms [4], CDI was limited to samples of compact and known support, a requirement needed in order to constrain the underdetermined data set of a far-field intensity pattern [5,6]. This restriction was lifted by so-called ptychographic CDI, replacing the support constraint with the overlap constraint [7,8]. Ptychography uses the redundancy of data collected by scanning a compact beam over an extended object with partial overlap between successive scan points. Proper sampling (or oversampling) conditions are achieved not by a compact object o, but a compact and stationary illumination wave field denoted as the probe p [9]. The technique is now applied in nanoscale imaging for example of biological specimens [10][11][12][13][14][15]. More important than its compatibility with extended samples is a second intriguing aspect of ptychography, namely the capability to reconstruct not only an unknown object o but also the complex-valued field of an unknown probe p [16], based on the Difference Map (DM), ePIE algorithm [17], or conjugate gradient implementations [18,19]. Ptychography has thus overcome the ubiquitous approximation of plane wave illumination in conventional x-ray diffraction. Coherent diffractive imaging is thus tractable also using distorted beams, i.e. wave fields with strong aberrations and phase front errors. Simultaneous characterization of beams and optical systems has thus also become possible [20][21][22]. It is important to discuss the restrictions applicable to the illumination wave front. For example, in most applications great care is undertaken to fulfill oversampling conditions limiting the allowable size of the illumination probe on the object by the Nyquist frequency of detector sampling. Ptychography is to date therefore primarily realized based on compact beams, as defined by pinholes or focusing optics. This condition of a 'compact probe' has recently been shown to be overly restrictive by Edo et al. [23], according to them the size of the illumination does not matter, provided that overlap and lateral diversity of the probe are sufficient. In practice, this means that an extended nearly plane wave illumination which exhibits little lateral diversity can often not be reconstructed unless the wave front is randomized as shown in [24]. To this end, Stockmar et al. have used a diffuser in x-ray near-field imaging [25]. By, lateral scanning of the object in an extended but structured probe they recover the object and the extended near-field wave front. They find that the incoming illumination must differ strongly from a uniform wave front, a condition which can in some cases be already satisfied by wave front distortions due to imperfect optics. It is directly understandable that efforts to simultaneously reconstruct o and p in the case of an extended illumination wave front by lateral scanning must fail if the wave fronts are close to plane waves, since the lack of phase diversity implies that all diffraction patterns are essentially equal. The goal of the present work is therefore to extend the seminal emancipation of phase reconstruction from the illuminating beam to the case of extended and arbitrary wave fronts, including both the case of nearly perfect (weakly phase shifting) wave fronts as well as strongly distorted illuminations, based on a full and simultaneous quantification of o and p. As shown below, phase diversity is created by translation of the object along the optical axis, i.e. in different defocus positions [26,27]. This data is then combined with two empty beam recordings, i.e. intensity recordings of the illuminating wave field without specimen, at different detector positions. From this set of data, a simultaneous reconstruction of object and probe can be achieved even for quasi unlimited sample and probe, as demonstrated by numerical simulations, using generalized ptychographic update equations for o and p. In contrast to the approach by Stockmar and coworkers [25], we thus use longitudinal rather than lateral scanning, rendering any need of a diffuser or wave front randomizier unnecessary.
Before we describe the reconstruction algorithm in the next section, followed by the presentation of the simulation results, let us briefly comment on the relevance of near-field x-ray imaging in extended beams. Near-field x-ray imaging using extended wave fronts in the Fresnel regime is a full field imaging technique similar to classical radiography but requiring partial coherence and treatment in correct wave optical terms [1,2,28]. Phase contrast emerges based on coherent free space propagation, leading to in-line holograms recorded at one or several defocus distances z j behind the object. Near-field imaging experiments can be carried out in quasi-parallel beams without magnification and hence macroscopic resolution, or using divergent (curved) illumination wave fronts [29-31], resulting in geometric magnification and microscopic resolution down to scales below 50nm. While higher resolution remains an asset of far-field CDI, a larger field of view (FOV) combined with still fairly high resolution provides an advantage of the near-field setting for many applications. In particular, nanoscale tomography of objects with diameters on the range of 100µm or larger [32, 33] is easily possible by adapting resolution and FOV in a divergent beam (defocus position). For small objects such as single biological cells, which are also amenable to high resolution CDI [34-36], near-field tomography is attractive because of its dose efficiency [37]. With respect to far-field imaging, near-field (propagation) imaging can thus combine high FOV and microscopic resolution down to below 50nm [10,37]. It also tolerates relaxed coherence and monochromaticity conditions, is dose efficient and leads to a more evenly distributed signal on the detector.
Significant efforts have been devoted to solve the phase problem in the near-field setting of propagation imaging using either the framework of the transport-of-intensity-equation [2,[38][39][40][41][42][43], an inversion scheme based on the oscillatory contrast transfer functions for phase and amplitude [28], or iterative algorithms such as the seminal Gerchberg and Saxton algorithm [44] and more elaborate iterative algorithms involving either a single [30,31,45] or several defocus distances [46]. Further work includes gradient based optimization tools again using a single [47,48] or several [27] near-field diffraction patterns. While previous near-field reconstruction techniques have included probe retrieval [49,50], the probe retrieval always required knowledge about the pupil function and had to be carried out prior to the imaging experiments.
In this work we present a ptychographic algorithm for simultaneous (intertwined) phase retrieval for o and p in the near-field setting of extended beams. Lateral scanning is replaced by longitudinal scanning along the optical axis. No additional information on the optical system (such as support in the pupil plane) is required. In section 2 we present the reconstruction algorithm, with some derivations placed in the appendix, followed by the simulation results in section 3, and a brief conclusion.

Reconstruction algorithm
The near-field dataset is obtained by illuminating the object o in N different defocus planes {1, ..., j, ..., N} along the optical axis. We assume a disturbed parallel beam illumination as the probe p which impedes an aberration free imaging. In each plane the complex valued field of the exit wave behind the object ψ j = o · p j is calculated by a point wise multiplication of object o and probe p j in plane j, in other words we assume the projection approximation to hold [31].
At each position along the optical axis, the object o is the same, whereas the probe p j is a propagated version of the probe from the preceding or following position. The optical field in plane j thus can be written as D ∆ j−1, j denotes the Fresnel near-field propagator written as (see [1]) where ∆ a,b is the distance between plane a and b, k = 2π λ is the wave number with wavelength λ , (x, y) are real space coordinates while (u, v) are reciprocal space coordinates. F is the Fourier transform and F −1 marks the inverse Fourier transform. Without loss of generality constant phase factors were set to one. A detector is placed at a distance z (defined from the first plane j = 1) to record the near-field intensity patterns The distance ∆ j,det is the distance from plane j to the detector. In addition two intensity patterns M p a , M p b of the empty beam are recorded in two different planes by translating the detector. The task is thus to reconstruct both o and p using the entire defocus series and the two empty beam recordings. For this purpose, two constraint sets and corresponding projectors are defined. The magnitude constraint P M projects the reconstructed near-field intensities in the detector plane onto the measured ones: with and with ε ≪ 1. For numerical stability this possibility to formulate the magnitude constraint was used. It can be found in [51]. The second constraint set should separate object and probe. (1) intensities (2) original and reconstructed object (3) original and reconstructed probe Gerchberg-Saxton on probe initial guess for object initial guess for probe p o * * j Fig. 1. Schematic of the phase retrieval algorithm. In the outer loop, a reference plane for the probe update is chosen. In the shown simulations we define the first plane to be this reference plane. In general any plane can be chosen. Three different error metrics are calculated (see Eqs. (11), (12) and (13)). The update plane j for the object is selected randomly. A Gerchberg-Saxton step using two empty beam measurements at two detector positions a and b is needed to predefine the probe. The complex valued field of the exit wave behind the object is calculated as the product of the object and the current, propagated probe. Next comes a modulus constraint according to the respective image recorded with the object in plane j, i.e. propagated intensities are set to the measured data, followed by back propagation. Object o and probe p are separated using to Eqs. (7) and (8). Following [16], the difference between the guessed ψ j and the separableψ j =ô ·p j needs to be minimized, i.e. the function Ideally, all S j should be zero and ψ j would be identical to the productψ j =ô · D ∆ 1, j [p 1 ]. To achieve this, we set the partial derivatives of all S j with respect toô andp j to zero. The analytical derivation (shown in the appendix) results in two symmetric and rather intuitive expressions forô andp 1 :ô as well asp Here * denotes the complex conjugated field. Of course, any of the defocus planes can be chosen as reference plane where the probe is updated, if the propagation operators are adapted accordingly. As the equations are coupled, they have to be applied sequentially. They compare well with the update equations given in [16] for separating object and probe in far-field ptychographic reconstruction. However, the two additional empty beam intensity recordings M p a , M p b used for a Gerchberg-Saxton step [44] nested into the update scheme are found to be indispensable to properly constrain the probe.
Putting all parts together, the algorithm for simultaneous reconstruction of object and probe as depicted in Fig. 1 can shortly be described as follows: An initial guess for object and probe (in plane one) is created. Usually homogeneous amplitude distributions for object and probe with zero phases are chosen as initial guess. A reference plane for the probe update is selected. For our simulation this reference plane is plane one and it is kept constant during the whole performance of the algorithm. A Gerchberg-Saxton step serves to predefine the probe. Next, a plane j of the object update is selected randomly. Care is taken that each plane can be chosen only once during one evaluation of the inner loop in Fig. 1. The exit wave ψ j in the update plane of the object is calculated by multiplication of o and p j . The magnitude constraint is applied to ψ j resulting in ψ ′ j = P M [ψ j ]. Now, the probe in plane one is updated followed by an update of the object in plane j These updates are similar to those used in [17]. The parameters α and β allow for an adjustable weighted sum of the current iterate with respect to the previous iterate. Setting α or β to zero is equal to not changing the initial guess, setting α or β to one does not take the solution of the iteration before into account. This entire scheme is iterated for several times. Three different kinds of errors are evaluated to monitor convergence. First the measured holograms are compared to the reconstructed ones: N is the number of holograms; K, L are the numbers of pixels in horizontal and vertical direction. The pixel size of squared pixels in one dimension is ∆d. Second the original object is compared to the reconstructed object: Finally the original probe is compared to the reconstructed probe: |p rec (k ·∆d, l · ∆d)| − |p orig (k · ∆d, l · ∆d)| |p orig (k ·∆d, l · ∆d)| 2 . (13) The terms 'original probe' and 'original object' refer to the original images used for object and probe as shown in Figs. 2(a) and 2(b). Of course, such error metrics only make sense when evaluating simulations. For real experiments the original object and probe are unknown and the error metrics should be calculated considering the values of the previous iterate for object and probe. To validate the algorithm and to check its performance by numerical simulation, object and probe were modeled using four photographs representing amplitude and phase of object ( Fig.  2(a)) and probe in plane one (Fig. 2(b)), respectively. To generate a data set, six near-field intensity patterns were simulated, three of them shown in Fig. 2(c) as well as two empty beam intensity recordings for two different detector positions (see Fig. 2(d)). For the first empty beam image, the detector was assumed to be at a distance z as for all defocus measurements, for the second empty image at position z/2; both values are with respect to the first plane. The Fresnel numbers of the holograms (see Table 1) are given for ten pixels each:

Simulation results
∆d is the pixel size in one dimension. As squared pixels are chosen for the simulations, only the Fresnel number in one direction is needed. The choice of six near-field intensity patterns was made for the purpose to insure convergence. We did observe that the algorithm can converge using a minimum of three images, provided sufficiently differing Fresnel numbers. All simulation parameters can be found in Table 2. At first we verified that the algorithm is capable to reconstruct the object having full information about the probe. For this purpose the probe was set to the original probe as shown in Fig. 2(b), and the reconstruction was run without any updates of p, i.e. neither the update using the separability constraint nor the Gerchberg-Saxton step. As initial guess for the object's amplitude and phase an array of ones (amplitude) and zeros (phase) was chosen. As demonstrated in Fig. 3(a), the known probe together with the intensity measurements of the defocus series is sufficient to reconstruct both amplitude and phase of the object. Contrarily, it was found that for the given data, the combination of known object (no update on o) and unknown probe is insufficient to reconstruct the probe, if no additional Gerchberg-Saxton cycle is used. In fact this observation prompted us to extend the scheme to incorporate two empty beam measurements. The failure of probe retrieval from the separability constraint is illustrated in Fig. 3(b). While the phase of the probe becomes visible, the amplitude seems to be disturbed by the phase information such that the original amplitudes appear only partially. This result indicates that the problem is under-determined concerning the probe, and that probe and object do not enter symmetrically in the problem as one may mistakenly assume. As expected from the above, simultaneous retrieval of an unknown object and probe is of course also impossible without the additional constraint of two empty images (Gerchberg-Saxton step). In fact, Fig. 3(c) shows the result of a reconstruction after 50 iterations (without the Gerchberg-Saxton step), initializing both o and p with ones (zeros for phase) and setting α = 0.3, β = 0.2 in Eqs. (9) and (10). While the phases of o and p are in general reconstructed in a semi-quantitative manner, amplitude reconstruction is strongly corrupted, as concluded from many similar reconstruction results which are compared to the respective input fields in Figs. 2(a) and 2(b). To compensate for the loss of phase information and to overcome under-specified data, the data recording scheme must therefore be extended to incorporate two independent empty beam intensity measurements at two different detector positions. Phasing of the probe in the reference plane is then significantly enhanced by a corresponding Gerchberg-Saxton step which acts in addition to the separability constraint. From many simulations, we find that this is fully sufficient to reconstruct both unknown object and probe. This finding is illustrated in Fig. 4(a): Object and probe can be reconstructed in amplitude and phase using a defocus series and two empty beam recordings. Fig.  4(b) shows the corresponding reconstruction of o and p for noisy intensity recordings. For this purpose, Poisson noise was simulated for 3.8 · 10 4 incident photons per pixel. Object and probe are very well reconstructed despite the high noise level. A somewhat grainy structure can be visually identified as a noise induced artifact. However, o and p can clearly be separated from each other. Convergence is illustrated in Fig. 4(c) for the ideal data set and in Fig. 4(d) for the noisy one. In both cases the algorithm converges with respect to the calculated error metrics.

Conclusion
In this article a novel concept of 'longitudinal ptychography' (i.e. exploiting diversity along the optical axis) with simultaneous probe and object retrieval in near-field imaging was presented. A procedure similar to the ePIE of [17] was used to reconstruct both object and probe in amplitude and phase from a set of near-field images and two empty beam recordings with different detector positions. Via minimization of the distance between the guessed exit wave and a complex valued function, which can be written as a product of an object transmission function and a propagated version of the probe, it is possible to derive Eq. (7) and Eq. (8). This then leads to an efficient routine of simultaneous reconstruction of object and probe. As the system is found to be under-determined else, two empty beam recordings treated by an additional Setting these derivatives to zero results in which is equation (7) andp Shiftingp j with the near-field propagator to plane one we arrive at Eq. (8).