Microsc. Microanal. 20, 1090–1099, 2014 doi:10.1017/S1431927614000932

© MICROSCOPY SOCIETY OF AMERICA 2014

Fast Deterministic Ptychographic Imaging Using X-Rays Ada W. C. Yan, Adrian J. D’Alfonso, Andrew J. Morgan, Corey T. Putkunz, and Leslie J. Allen* School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia

Abstract: We present a deterministic approach to the ptychographic retrieval of the wave at the exit surface of a specimen of condensed matter illuminated by X-rays. The method is based on the solution of an overdetermined set of linear equations, and is robust to measurement noise. The set of linear equations is efficiently solved using the conjugate gradient least-squares method implemented using fast Fourier transforms. The method is demonstrated using a data set obtained from a gold–chromium nanostructured test object. It is shown that the transmission function retrieved by this linear method is quantitatively comparable with established methods of ptychography, with a large decrease in computational time, and is thus a good candidate for real-time reconstruction. Key words: coherent diffractive imaging, exit-wave reconstruction, ptychography, lensless imaging, phase retrieval, X-ray imaging

INTRODUCTION Determination of the location and type of atoms in a specimen of condensed matter is the basis for understanding the structural and electronic properties of materials, and the mechanisms by which biological structures function. The leading methods of structure determination, X-ray crystallography, and electron microscopy using lenses, have their limits; some biological specimens cannot be crystallized (Ostermeier & Michel, 1997), and lens aberrations limit the resolution of electron microscopy (Scherzer, 1936; Hawkes, 2009). Since 1997, aberration-corrected electron microscopes have been available, increasing the ease of obtaining atomic-resolution images. However, the best resolution achieved so far is about 0.05 nm (Hawkes, 2009), still about 25 times poorer than the theoretical Rayleigh limit, and further correction leads to diminishing returns. This is because higher order aberrations and temporal incoherence of the electron beam begin to dominate. It has recently been shown that magnetic field noise from thermally driven currents in the conductive parts of electron microscopes causes decoherence in the beam, so that having a chromatic aberration corrector to reduce temporal incoherence can actually introduce other aberrations (Uhlemann et al., 2013). This may fundamentally limit the resolution obtained, which motivates the development of lensless electron imaging. Similarly, lensless X-ray imaging is motivated by the difficulty of manufacturing focusing optics for X-rays, such as Fresnel zone plates, to the precision necessary for the specimen to be reconstructed at the Rayleigh limit.

Received December 1, 2013; accepted April 21, 2014 *Corresponding author. [email protected]

Coherent diffractive imaging (CDI) was first demonstrated experimentally in 1999 by Miao et al. (1999), and has the potential to overcome the shortcomings of X-ray crystallography and conventional transmission electron microscopy. It does not require a crystalline specimen, and is lensless, in that lenses are only used to focus the incident beam, and not to image the specimen; the sample structure is retrieved from its diffraction pattern, and so the reconstruction is diffraction limited rather than aberration limited. However, specimen damage is still of concern, especially for biological samples, and it is thus desired that the dosage of illumination be kept as low as possible (Henderson, 1995). Nonlinear iterative retrieval algorithms are commonly used to solve the phase problem. However, they are not deterministic: different starting points may lead to different final reconstructions. A common approach is, therefore, to perform many reconstructions using different starting points to check if the solution reached is unique (Chapman et al., 2006), or to average said reconstructions to eliminate high spatial frequency artifacts (Shapiro et al., 2005). This nonuniqueness is because the equations being solved in the reconstruction process are nonlinear with respect to the specimen transmission function, which describes the interaction of the specimen with the illumination. Hence, when these algorithms search for the global minimum of an error metric in solution space, the error metric is nonconvex, and the reconstruction may converge to a local minimum of the error metric in solution space. It is also not possible to determine beforehand the number of iterations required for convergence, although empirical evidence in imaging and crystallographic applications shows that convergence usually occurs in a reasonable number of iterations (Elser, 2003). Furthermore, a large proportion of nonlinear iterative retrieval algorithms have adjustable parameters related to the

Fast Deterministic Ptychographic Imaging Using X-Rays

convergence rate, which are manually altered for different data sets. There exist noniterative, deterministic phase retrieval methods such as the direct numerical solution of the transport of intensity equation (Gureyev et al., 1995); however, the method is inapplicable when there are phase vortices in the specimen transmission function, and images are taken from the near-field rather than the far-field, which is usually accomplished with lenses. A fast, deterministic algorithm that is robust to noise is desirable, such that a unique solution is guaranteed, and real-time analysis can be conducted in an automated manner. The concept of ptychography was developed about 40 years ago (Hoppe,1969), but it is only recently that practical algorithms have been developed to solve the inverse problem (Rodenburg & Faulkner, 2004; Thibault et al., 2008). In ptychography, multiple diffraction patterns are used to reconstruct a single transmission function. These diffraction patterns may be in the form of a through focal series, or longitudinal phase diversity, where the defocus of the illuminating probe is changed; or transverse phase diversity, where the probe is shifted across the specimen, and diffraction patterns are taken at overlapping probe positions. Ptychography overdetermines the problem to be solved: in other words, there is more information than necessary to reconstruct the transmission function of the specimen, which describes the interaction of the illumination with the specimen. It imposes the constraint that the transmission function must be consistent across all of the diffraction patterns; this makes the reconstruction more robust with respect to measurement noise. It has been shown experimentally that a retrieval using several probe positions with transverse and longitudinal phase diversity is more accurate than that using a single probe position with the equivalent dosage (Putkunz et al., 2011), and that the signal-to-noise ratio and convergence rate increases with the number of diffraction patterns used (Zhang et al., 2007; Bao et al., 2008). Furthermore, since the illumination covers different areas when translated, ptychography can extend the field of view over which the specimen transmission function is reconstructed. It is for these reasons we have extended our fast, deterministic method for reconstruction of the specimen transmission function to ptychography (D’Alfonso et al., 2013).

M ATERIALS AND M ETHODS Algorithms In this section we will summarize the linear approach to ptychography used in this paper, global ptychographic iterative linear retrieval using Fourier transforms (GPILRUFT), as first presented by D’Alfonso et al. (2014) and based on the original linear phase reconstruction method by Martin & Allen (2008), and compare this to established nonlinear methods, scanning X-ray diffraction microscopy (SXDM) (Thibault et al., 2008), and the ptychographic iterative engine (PIE) (Rodenburg & Faulkner, 2004).

1091

In ptychographic CDI it is assumed that the wave at the exit surface of the specimen for each position of the illuminating probe can be written in the form (1) ψ exit;j ðrÞ ¼ ψ ill;j ðrÞTðrÞ; where r is the real space coordinates in the object plane, ψill (r) the illumination in real space, j the probe position number, and the specimen transmission function T(r) may conveniently be expressed as TðrÞ ¼ 1 + OðrÞ:

(2)

The validity of the wave factorization assumption in equation (1) is expected to be good for X-ray scattering and is discussed further in Thibault et al. (2008) and Rodenburg & Bates (1992). The transmission function T(r) describes the interaction of the illumination with the specimen. For overlapping probe positions, the reconstructed transmission function must be consistent with all of the reconstructed exit surface waves in the region of overlap, and this consistency is enforced as a constraint, overdetermining the problem. GPILRUFT The linear retrieval approach of Martin & Allen (2008) uses generalized holography to reconstruct the wave at the exit surface of the specimen from its autocorrelation. In conventional holography, the unknown exit surface wave is interfered with a separate reference wave: a pinhole for Fourier holography (Stroke et al., 1965), or a plane wave for off-axis holography (Leith & Upatnieks, 1962). In the generalization by Martin and Allen, the unscattered part of the illumination is used as an effective reference wave. The autocorrelation of a given exit surface wave, which by the Wiener–Khinchin theorem is the inverse Fourier transform of the intensity of the diffraction pattern (Wiener, 1930), can be written as Z -1 F ½I ðqÞ¼ ½ψ ill ðr+r0 ÞOðr+r0 Þψ *ill ðr0 ÞO* ðr0 Þ +ψ ill ðr+r0 ÞOðr+r0 Þψ *ill ðr0 Þ+ψ ill ðr+r0 Þψ *ill ðr0 ÞO* ðr0 Þ +ψ ill ðr+r0 Þψ *ill ðr0 Þdr0  ðψ ill OÞ ðψ ill OÞ +ðψ ill OÞψ ill +ψ ill  ðψ ill OÞ+ψ ill ψ ill ;

ð3Þ

where q is the reciprocal space coordinates in the detector plane, I(q) the intensity of the far-field diffraction pattern, and ⊗ a cross-correlation: Z AB ½A* ðr0 ÞBðr+r0 Þdr0 : (4) Assuming the illumination is well characterized, we can subtract the autocorrelation of the illumination [the last term in equation (3)] to obtain ψ ill  ðψ ill OÞ + ðψ ill OÞ  ψ ill + ðψ ill OÞ  ðψ ill OÞ ¼ F - 1 ½I  - ψ ill  ψ ill :

(5)

1092

Ada W. C. Yan et al.

This is a set of equations for O; for each pixel in the autocorrelation of the exit surface wave, we can construct a complex equation for the complex pixels of O. Finding the least-squares solution to a set of simultaneous equations is equivalent to finding the global minimum of some error metric. For nonlinear equations, the error metric to be minimized has at least one global minimum in solution space, but, in general, also many local minima, such that search algorithms can stagnate in local minima instead of reaching the true solution. For linear equations of the form Ax = b where x is our trial solution, a metric that can be minimized is the residual: ϵr ¼k Ax - b k2 ;

(6)

where the L2-norm is used, and it turns out that minimizing the residual is equivalent to solving the normal equation (Hestenes & Stiefel, 1952): AT Ax ¼ AT b:

(7)

As such, the metric is a quadratic function in solution space; as there is only one minimum in solution space, stagnation in a local minimum does not occur. This means that solving a linear equation is deterministic; there is no need to seed the algorithm with different starting points to check the uniqueness of the solution, or to average retrieved transmission functions to obtain the final solution. This motivates solving the phase problem using linear equations. To make the set of equations linear, we wish to exclude the third term on the left-hand side of equation (5), which has a second-order contribution in O. One approach is to only use equations from the area in the autocorrelation of the exit surface wave that contains only terms linear in O (Martin & Allen, 2008). However, in practice, much of the intensity of the autocorrelation of the exit surface wave is concentrated in the area that contains these nonlinear terms. As measurement error can be modeled as Poissonian, the autocorrelation of the exit surface wave will have approximately uniform noise across the autocorrelation plane, so the area of the autocorrelation of the exit surface wave that contains nonlinear terms is also the area least affected by measurement noise. Thus, we would like to include this area when constructing the linear equations. If the interaction between the illumination and the specimen is sufficiently weak, we can simply ignore the nonlinear term ðψ ill OÞ  ðψ ill OÞ in equation (5) and obtain an equation for an approximate solution (Morgan et al., 2013): ψ ill  ðψ ill OÞ + ðψ ill OÞ  ψ ill  F - 1 ½I  - ψ ill  ψ ill : (8) In discrete notation, equation (8) is equivalent to X

*ill * *ill ill ψ ill n + m ψ n On + ψ nm ψ n On

n

¼ F - 1 ½I m -

X n

ill ψ *ill nm ψ n ;

(9)

where the subscripts in equation (8) have been raised to superscripts to facilitate the discrete notation. Equation (9) can be written as a matrix multiplication Ax ¼ b;

(10)

where A is constructed using the known illumination, x is the vectorial form of the unknown object wave function O and b is constructed using the measured diffraction pattern and known illumination. The operation of A on x is the sum of the two cross-correlations on the left-hand side of equation (8). The matrix equation can be solved efficiently using the conjugate gradient least-squares method (CGLS) (Hestenes & Stiefel, 1952) implemented using fast Fourier transforms (D’Alfonso et al., 2012). The advantage of the CGLS algorithm is that A need not be constructed explicitly; we only need to know the results of multiplying A and AT by an arbitrary vector. Explicitly, the results of multiplying A and AT by an arbitrary vector can be calculated using (D’Alfonso et al., 2012) n o Ap ¼ F - 1 ψ^ ill ½F ðψ ill ϕÞ* + ψ^ *ill ½F ðψ ill ϕÞ   AT d ¼ ψ *ill F - 1 ψ^ ill^γ * + ψ^ ill^γ ;

ð11Þ

where ψ^ is the Fourier transform of ψ, ϕ the function representation of p and γ the functional representation of d. Furthermore, a CGLS reconstruction exhibits semiconvergence; in other words, the components of the solution fitted first are those least affected by errors on b, such as those caused by the second-order contribution. The components of the solution affected by the error on b appear at later iterations. Hence, we can truncate the CGLS reconstruction to obtain an optimal solution (Hansen, 2010). In general, the solution of equation (8) will be different from the true object O that is the solution of equation (5), but it will be a good approximation if the second-order contribution is small and we truncate the CGLS reconstruction such that the later iterations, which are affected by the second-order contribution, are not included (Morgan et al., 2013). We can correct for the second-order contribution, as discussed in the Regularization in GPILRUFT section. In GPILRUFT, the A-matrices associated with all probe positions and the b-vectors formed from each diffraction pattern are concatenated (D’Alfonso et al., 2013). For example, if the object is entirely covered by every illumination, equation (10) becomes 2 3 2 3 ½A1  ½b1  6 ½A2  7 6 ½b2  7 6 7 6 7 (12) 6 .. 7½x ¼ 6 .. 7; 4 . 5 4 . 5 ½AJ  ½bJ  where J is the number of probe positions. The generalisation of Eq. 12 to extended objects is trivial. The algorithm can be extended to update the illumination wave functions and the retrieved transmission function simultaneously (D’Alfonso et al., 2013). However, this requires that the illumination (up to a transverse and/or longitudinal

Fast Deterministic Ptychographic Imaging Using X-Rays

translation) be the same for each probe position. In the work that follows, this will not be the case. Regularization in GPILRUFT In contrast to previous work by D’Alfonso et al. (2013) and Morgan et al. (2013),  we consider  the case where the nonlinear term ψ ill;j O  ψ ill;j O is significant compared with the linear terms, so a correction scheme must be used to take into account its contribution to the equations. If we have ~ for the unknown object, we can construct the an estimate, O, ~ can be obtained by first nonlinear term with this estimate. O solving (9)  equation   without correction, i.e., assuming ψ ill;j O  ψ ill;j O ¼ 0, then truncating the reconstruction after a certain number of CGLS iterations; the optimum cutoff number of iterations will be discussed later on in this subsection. ~ ðrÞ + OD ðrÞ, where OD ðrÞ is the difWe can let OðrÞ ¼ O ference between the estimated and true object wave functions, and expand the nonlinear term in equation (5) as follows         ~  ψ ill;j O ~ ψ ill;j O  ψ ill;j O ¼ ψ ill;j O         ~  ψ ill;j OD + ψ ill;j OD  ψ ill;j O ~ + ψ ill;j O     + ψ ill;j OD  ψ ill;j OD : ð13Þ Equation (5) can then be written as       ψ ill;j  ψ ill;j O + ψ ill;j O  ψ ill;j ¼ F - 1 Ij     ~  ψ ill;j O ~ - ψ ill;j  ψ ill;j - ψ ill;j O         ~  ψ ill;j OD - ψ ill;j OD  ψ ill;j O ~ - ψ ill;j O    - ψ ill;j OD ψ ill;j OD :

ð14Þ

If OD is small, the last three terms in equation (14) can be ignored. Following the approach of equation (8), we can ignore these terms and solve the equation       ψ ill;j  ψ ill;j O + ψ ill;j O  ψ ill;j ¼ F - 1 Ij     ~  ψ ill;j O ~ : - ψ ill;j  ψ ill;j - ψ ill;j O ð15Þ Compared with equation (8), the degree of approximation required to equate the two sides is reduced. This is the approach used by Morgan et al. (2013), where it was found that the nonlinear term was small for their particular data set, and so the correction produced little change in the retrieved transmission function. It is possible to improve on the approach in equation (15). It has been shown that a linear retrieval can be applied to the ptychographic problem in a serial way, using the part of the exit surface wave retrieved by one probe position as a reference wave to retrieve the unknown part of the exit surface wave at a different probe position (D’Alfonso et al., 2012).

1093

A variation of this is to use our current estimate of the exit surface wave:   ~ ðrÞ ; (16) ψ ref ;j ðrÞ  ψ ill;j ðrÞ 1 + O as the reference wave to recast the set of linear equations. If we cross-correlate ψ ill;j ðrÞOðrÞwith ψ ref ;j ðrÞ, we obtain       ψ ref ;j  ψ ill;j O + ψ ill;j O  ψ ref ;j ¼ ψ ill;j  ψ ill;j O       ~  ψ ill;j O + ψ ill;j O  ψ ill;j + ψ ill;j O     ~ ð17aÞ + ψ ill;j O  ψ ill;j O         ¼ ψ ill;j  ψ ill;j O + ψ ill;j O  ψ ill;j + ψ ill;j O  ψ ill;j O         ~  ψ ill;j O ~ - ψ ill;j OD  ψ ill;j OD + ψ ill;j O ð17bÞ       ~  ψ ill;j O ~ ¼ F - 1 Ij - ψ ill;j  ψ ill;j + ψ ill;j O     - ψ ill;j OD  ψ ill;j OD :

ð17cÞ

Equation (17b) is obtained from equation (17a) by substitution of equation (13); subsequent substitution of equation (5) gives equation (17c). Thus, we have constructed another set of (almost) linear equations for O using the known data. The only nonlinear term is the last term of equation (17c), and if the object estimate is close to the true object, this nonlinear term will be small. Following the approach of equation (8), we can ignore this term and solve the equation     ψ ref;j  ψ ill;j O + ψ ill;j O  ψ ref ;j       ~  ψ ill;j O ~ : ð18Þ ¼ F - 1 Ij - ψ ill;j  ψ ill;j + ψ ill;j O   ~  The approximation  (15) is that    in equation  ψ ill;jO ~ + ψ ill;j OD  ψ ill;j OD ¼ 0, ψ ill;j OD + ψ ill;j OD  ψ ill;j O in equation (18) is that  while the approximation  ψ ill;j OD  ψ ill;j OD ¼ 0. Since the terms discarded in the approximation in equation (15) are larger than those in equation (18), equation (18) is more exact than equation (14) given the same diffraction data. In either case, a correction to b is applied by either adding  or subtracting   the concatenated vec~  ψ ill;j O ~ . In addition, when tors that represent ψ ill;j O constructing the matrix equation Ax = b from equation (18), the matrix products become Ap ¼ t  t1kktJ n    * *   o ¼F -1 ψ^ ref ;1 F ψ ill;1 ϕ +^ ψ ref ;1 F ψ ill;1 ϕ   * *   o ψ ref;J F ψ ill;J ϕ k¼kF -1 fψ^ ref;1 F ψ ill;J ϕ +^   X ψ ref ;j^γ : AT d¼ ψ *ill;j F -1 ψ^ ref ;j^γ *+^

ð19Þ

j

The question remains as to how to truncate the uncor~ The point rected CGLS solution of equation (9) to obtain O.

1094

Ada W. C. Yan et al.

albeit slowly, whereas methods such as the difference map algorithm have a faster convergence rate, but diverge when the diffraction patterns are affected by measurement noise. Hence, the modulus projection is chosen for this data set. The progress of the reconstruction may be monitored using the modulus error metric ϵm ¼

2 X  Pm ψ exit;j - ψ exit;j  :

(22)

j

Figure 1. Flowchart for correction of the nonlinear term in GPILRUFT. GPILRUFT, global ptychographic iterative linear retrieval using Fourier transforms.

at which to truncate and correct should be chosen to maxi~ mize the convergence rate—if truncation occurs too soon, O will not have sufficient detail to be a useful estimate, but if it ~ will contain artifacts due to noise and occurs too late, O ~ will not be small, nonlinear terms. In both cases, OD ¼ O - O and the correction will be ineffective. D’Alfonso et al. (2013) used the L-curve (Hansen, 1992) to truncate the uncorrected CGLS solution, but this can be ineffective if the nonlinear term is large. In this case, an iterative correction scheme where n correction runs are performed, with the ith correction run truncating at i CGLS iterations then correcting, is found to be effective. The retrieved object from each correction run, xi, is retained as the starting guess x0 for CGLS in the next correction run. A flowchart for this procedure is given in Figure 1. SXDM SXDM (Thibault et al., 2008) is a nonlinear ptychographic iterative algorithm, where the transmission function is updated using the exit surface waves at all probe positions simultaneously. At each iteration, the transmission function is updated according to P j ψ ill;j ðrÞψ exit;i;j ðrÞ Ti + 1 ðrÞ ¼ ; (20) 2 P   ψ ð r Þ   j ill;j where j denotes the probe position and i the iteration number. New exit surface waves are constructed by multiplying the new transmission function by the illumination at each probe position. Then each exit surface wave is updated via the modulus projection: " # pffiffi ψ_ exit -1 Pm ψ exit ¼ F I _  ; (21) ψ exit  or using alternative approaches, such as the hybrid input– output (Fienup, 1978) or difference map (Elser, 2003) algorithms. This constitutes one complete iteration. For our data set, simulation shows that a simple modulus projection as per equation (21) enables the reconstruction to converge,

Like GPILRUFT, SXDM can also be implemented such that the illumination is updated along with the transmission function. Once again, this requires that the illumination be the same for each probe position, aside from translation and defocus, which is not the case for the data presented in this paper (for details, see Thibault et al., 2008). PIE Unlike SXDM, the PIE updates the transmission function by addressing each probe position serially (Rodenburg & Faulkner, 2004). Each iteration of PIE consists of J updates, where J is the number of probe positions. For each update, a probe position is selected at random and the following update made: ψ *ill;j ðrÞ Ti + 1 ðrÞ ¼ Ti ðrÞ + α  2   ψ ill;j ðrÞ

h

Pm ψ exit;i;j ðrÞ - ψ exit;i ðrÞ

i

max

where ψ exit;i ðrÞ ¼ ψ ill;j ðrÞTi ðrÞ:

ð23Þ

Here α is a scale factor related to the step size in a steepest descent method, set to unity in this work; i denotes the update number. The probe position for the next update is then chosen randomly from the subset of probe positions that have not already been considered. Once all probe positions have been addressed we consider this a complete iteration. The progress of the reconstruction can be monitored using the modulus error metric defined in equation (22). The extended PIE (Maiden & Rodenburg, 2009) extends the PIE algorithm to correct the illumination; once again, this requires that the illumination be the same for each probe position, which is not the case for the data presented in this paper. For further details, see Maiden & Rodenburg (2009).

Experimental Set-Up Diffraction data for a gold–chromium (Au–Cr) nanostructured test object was taken using a dedicated Fresnel imaging end station situated at beam line 2-ID-B of the Advanced Photon Source, Argonne National Laboratory. This data has previously been presented in Putkunz et al. (2011), where the test object’s transmission function was reconstructed using a variant of SXDM. The test object was mounted in vacuo on a silicon nitride (Si3N4) window. A 160 µm diameter Fresnel zone plate was used to longitudinally shift the focus of the 2.5 keV X-ray illumination incident on the specimen. Longitudinal

Fast Deterministic Ptychographic Imaging Using X-Rays a

c

b

a

b

c

d

e

f

1095

d

Figure 2. a: The intensity of the illumination at the object plane at the first probe position. The object area is outlined in red. b: The phase of the illumination at the detector plane in the absence of a specimen, where the exponential phase factor due to defocus has been removed for clarity. c: An SEM image of the nanostructured Au–Cr object. d: The diffraction pattern for the first probe position and associated defocus. SEM, scanning electron microscopy; Au–Cr, gold–chromium.

and transverse phase-diverse data were obtained for seven probe positions. At each probe position, a diffraction pattern is taken at a single defocus; there are five distinct defoci across the seven probe positions, ranging from 1.8 overfocus to 2.2 mm overfocus in steps of 100 µm. The diffraction patterns are recorded using a 2,048 × 2,048 direct detection cooled CCD with pixels 13.5 µm wide. Measurements of the unperturbed illumination were taken for each position and reconstructed using the method of Quiney et al. (2006). The real space illumination intensity at the first probe position is shown in Figure 2a, and the phase of the illumination at the detector in the absence of a specimen is shown in Figure 2b, where the exponential factor owing to defocus has been removed for clarity. The assumed object area is outlined in red. Figure 2c shows a scanning electron microscopy image of the nanostructured Au–Cr object, and Figure 2d the diffraction pattern for the first probe position and associated defocus. For more details on the experimental set-up, refer to Putkunz et al. (2011).

RESULTS AND D ISCUSSION Figure 3a shows the best (minimum standard deviation) reconstructed phase of the transmission function using GPILRUFT without correction for the nonlinear term, at 50 CGLS iterations. The retrieval is similar to that presented in Putkunz et al. (2011), showing the feasibility of GPILRUFT as an alternative reconstruction method. However, it is known

Figure 3. a: GPILRUFT reconstruction of the transmission function phase with no correction (50 iterations). b: The three regions in the specimen transmission function that are known to be of uniform phase. The line scans in Figures 5, 8, and 9 are taken along the red arrow. The GPILRUFT reconstruction in (a) is included to show the relative position between the specimen and the three regions. c, d: GPILRUFT reconstruction of the transmission function phase with (c) correction by subtracting the estimated nonlinear term [equation (15)]; (d) correction using the reference wave approach [equation (18)] at 20 correction runs (210 iterations) each. e: SXDM reconstruction (4,000 iterations). f: PIE reconstruction (2,000 iterations). For the implementations in this paper, each iteration of CGLS, SXDM, and PIE takes ∼1 s, while either subtracting the nonlinear term or constructing the reference wave takes ∼4 s for each iteration where correction takes place. Because SXDM and PIE retrieve the phase up to an arbitrary constant, the phase in (e) and (f) is adjusted by a constant such that the mean phase within the object area is equal to that in (d). GPILRUFT, global ptychographic iterative linear retrieval using Fourier transforms; SXDM, scanning X-ray diffraction microscopy; PIE, ptychographic iterative engine; CGLS, conjugate gradient least-squares method.

that the three regions of the sample shown in Figure 3b are of roughly uniform phase, and this uniformity has not been reproduced—we can see “rings” that are artifacts of the reconstruction, as indicated by the arrow. Figure 3c shows the reconstructed phase with correction by subtracting the estimated nonlinear term [equation (15)]; (d) uses the reference wave correction in equation (18). Reconstructions (c) and (d) are truncated after 20 correction runs (a total of 210 CGLS iterations). For all of the reconstruction methods, the number of iterations at which to truncate is based on the knowledge that the three regions outlined in Figure 3b should be uniform in phase. When both the mean and standard deviation of the phase of the retrieved transmission function in each of the three regions

1096

Ada W. C. Yan et al.

Figure 4. Line scans of the retrieved phase of the specimen transmission function taken along the red arrow in Figure 3b, for GPILRUFT without correction (Fig. 3a), with correction by subtracting the nonlinear term (Fig. 3c), and using the reference wave (Fig. 3d).

stagnate, we consider the reconstruction converged and truncate the reconstruction. Hence, the reconstructions in Figures 3c–3f are truncated at 210, 210, 4,000, and 2,000 iterations, respectively. Further iterations did not change the reconstructions appreciably. Comparing Figure 3c with Figure 3d, using the reference wave for correction removes the artifacts in (a) more effectively than subtracting the nonlinear terms. The uniformity of the phase of the reconstructed transmission function can be quantified using line scans along the red arrow in Figure 3b. Figure 4 shows line scans for GPILRUFT without correction (Fig. 3a), with correction by subtracting the nonlinear term (Fig. 3c), and with correction using the reference wave (Fig. 3d). We see that in the three regions where the phase is known to be uniform, correction using the reference wave gives the most uniform phase. We can compare the reconstruction using GPILRUFT to reconstructions using nonlinear algorithms. Figure 3e shows the phase reconstructed using SXDM at 4,000 iterations. This is qualitatively comparable to Figure 3d; however, the number of iterations required is an order of magnitude larger. Figure 3f shows the reconstructed phase of the transmission function at 2,000 iterations of PIE. Once again, the number of iterations required is an order of magnitude larger. Furthermore, artifacts such as truncation of the object on the left-hand side are visible (as indicated by the arrow), and the background is not uniform as it should be, given this area corresponds to vacuum. Figure 5 shows line scans of the retrieved phase of the specimen transmission function taken along the direction that is indicated by the red line in Figure 3b, for GPILRUFT using reference wave correction (Fig. 3d), SXDM (Fig. 3e), and PIE (Fig. 3f). We see that within the regions of uniform phase, the phase of the transmission function reconstructed using GPILRUFT varies much more rapidly than that reconstructed using SXDM and PIE. However, the

Figure 5. Line scans of the retrieved phase of the specimen transmission function taken along the red arrow in Figure 3b, for GPILRUFT with reference wave correction (Fig. 3d), SXDM (Fig. 3e), and PIE (Fig. 3f). SXDM, scanning X-ray diffraction microscopy; PIE, ptychographic iterative engine.

magnitude of the fluctuations is small, and comparable with the magnitude of fluctuations for SXDM and PIE within the uniform regions. Also, we can see that while the phase profiles are generally similar across the different methods, for PIE the magnitude of the phase increases much more slowly on the left-hand edge of the specimen, as indicated by the green arrow. This corresponds to the observed truncation of the left edge of the specimen indicated by the arrow in Figure 3f. Furthermore, the phase difference between the three regions reconstructed using PIE is different from that reconstructed using GPILRUFT and SXDM. We contend that the high-frequency artifacts in the GPILRUFT reconstruction in Figure 3d are mainly due to discrepancies between the reconstructed illumination and the actual illumination. In the absence of a specimen, the reciprocal space illumination for all probe positions in the data set should be the same except for a phase factor due to defocus and/or translation, but it actually changes slightly, due to a windowed air gap ∼5 m upstream of the Fresnel zone plate that produces Fresnel diffraction fringes at the zone plate. The motion of this window during data collection causes inconsistency in the illumination between different probe positions. Hence, we cannot update the illumination using any of the three algorithms. For the same reason, the assumed illumination, which is constructed from white-field images of the illumination where the specimen was absent using the method of Quiney et al. (2006), does not accurately reflect the incident illumination when the specimen was present. Measurement noise on the images from which the illumination is reconstructed also leads to difficulty in characterizing the illumination. We use a simulated specimen transmission function to test this hypothesis. The simulated transmission intensity and phase are constructed to mimic the experimentally reconstructed specimen transmission function, and are shown in Figures 6a and 6b, respectively. Diffraction patterns are

Fast Deterministic Ptychographic Imaging Using X-Rays a

1097

b

Figure 6. (a) Intensity and (b) phase of simulated transmission function, where the phase in the three uniform regions is indicated. a

b

c

d

e

f

Figure 7. Phase of transmission function reconstructed from simulated data: (a) GPILRUFT (with reference wave correction) (b) SXDM, (c) PIE, each at 2,000 iterations. Phase of transmission function reconstructed using assumed illumination that differs from the illumination used to form the diffraction patterns: (d) GPILRUFT, (e) SXDM, (f) PIE at 400, 4,000, and 2,000 iterations, respectively. The phase for SXDM and PIE is adjusted by a constant such that the mean phase within the object area is equal to that of the simulated transmission function. GPILRUFT, global ptychographic iterative linear retrieval using Fourier transforms; SXDM, scanning X-ray diffraction microscopy; PIE, ptychographic iterative engine.

constructed using the recorded experimental illumination and simulated transmission function. GPILRUFT, SXDM, and PIE are then used to reconstruct the specimen transmission function from the simulated diffraction patterns. Figures 7a–7c show the transmission function phase reconstructed using GPILRUFT (with reference wave correction for the nonlinear terms), SXDM, and PIE, respectively, with the simulated transmission function as input, at 2,000 iterations. It is found that although the mean phase of the transmission function stagnates as the number of iterations increases, its standard deviation continues to decrease

Figure 8. Line scans of the retrieved phase of the specimen transmission function taken in the direction that is indicated by the red line in Figure 3b, for GPILRUFT with reference wave correction (Fig. 7a), SXDM (Fig. 7b), and PIE (Fig. 7c). A line scan of the input transmission function is shown for reference. GPILRUFT, global ptychographic iterative linear retrieval using Fourier transforms; SXDM, scanning X-ray diffraction microscopy; PIE, ptychographic iterative engine.

beyond machine precision, such that for all three reconstruction methods, the quality of retrieval continues to increase with the number of iterations. Since there is no point of stagnation, it is natural to compare reconstructions for each of the methods at the same number of iterations. We see that when the illumination is well characterized, all three reconstructions converge to the simulated transmission function, and are free of artifacts. The line scan in Figure 8 shows that the phase is uniform in all three regions in Figure 3b for GPILRUFT and PIE, but much less uniform for SXDM. Figures 7d–7f show the phase of the transmission function reconstructed using GPILRUFT, SXDM, and PIE, respectively, where the assumed illumination differs from the illumination with which the diffraction patterns are formed. The assumed illumination is simulated by randomly swapping the reciprocal space illumination (without defocus or translation) between probe positions, then normalizing and applying the correct defocus and translations to construct a new illumination wave function at each probe position that has the same position and summed intensity as the original illumination. We note that the discrepancy between assumed and actual illumination would be much smaller in experiment. The same diffraction patterns for Figures 7a–7c are used, but the assumed illumination is used to solve for the transmission function from the diffraction patterns. In this case, the GPILRUFT and SXDM reconstructions converge at 400 and 4,000 iterations, respectively, and Figures 7d and 7e show the phase of the transmission function at these iterations. For the PIE reconstruction, the mean phase of the specimen transmission function in region 2 diverges as the number of iterations increases, so we truncate the reconstruction when the mean and standard deviation of

1098

Ada W. C. Yan et al.

Figure 9. Line scans of the retrieved phase of the specimen transmission function taken in the direction that is indicated by the red line in Figure 3b, for GPILRUFT with reference wave correction (Fig. 7d), SXDM (Fig. 7e), and PIE (Fig. 7f), reconstructed using assumed illumination that differs from the illumination used to form the diffraction patterns. GPILRUFT, global ptychographic iterative linear retrieval using Fourier transforms; SXDM, scanning X-ray diffraction microscopy; PIE, ptychographic iterative engine.

the phase of the specimen transmission function in the other two regions stagnates at 2,000 iterations. We see that highfrequency “ring” artifacts are present in the GPILRUFT reconstruction, but not in the SXDM and PIE reconstructions, which have more slowly varying artifacts. In addition, the SXDM and PIE backgrounds become nonuniform. This matches what we see in experiment, and is an artifact because we know that the transmission function phase should be zero outside the object area. Figure 9 shows line scans of the retrieved phase of the specimen transmission function taken in the direction that is indicated by the red arrow in Figure 3b, for GPILRUFT using reference wave correction (Fig. 7d), SXDM (Fig. 7e), and PIE (Fig. 7f), where the assumed illumination differs from the illumination with which the diffraction patterns are formed. We see that for SXDM and PIE, the fluctuations in phase within the three uniform regions are of lower frequency than for GPILRUFT. However, for PIE the phase differences between the three regions is not reproduced. Also, the phase within each region is not uniform, especially in region 3. We conclude that the high-frequency artifacts in the GPILRUFT reconstruction are due to discrepancies between the assumed and actual illumination, which causes the reconstruction to stagnate prematurely, such that the reconstructed phase is less uniform in the three regions. The inconsistency between assumed and actual illumination affects the SXDM and PIE algorithms less because the essential step, modulus projection using equation (21), involves only the exit surface waves. The illumination is only used as a weighting factor to extract the transmission function from the exit surface waves using equations (20) and (23), and to construct trial exit surface waves for the modulus projection. On the other hand, knowledge of the illumination is required

for both matrix multiplication steps of the CGLS algorithm, i.e., equation (11) or, if using reference wave correction, equation (19). The correction step, either forming the secondorder nonlinear term to be subtracted or forming the reference wave, also requires knowledge of the illumination. Therefore, if the illumination is not known to have high accuracy, the nonlinear terms cannot be taken into account accurately. If the illumination is the same across all probe positions, one can update the illumination to ameliorate the effect of inaccurately determined illumination; alternatively, one can use more sets of phase-diverse data and further overdetermine the problem. Another approach is to use GPILRUFT to quickly obtain an estimate for the specimen transmission function, then use this as a starting guess for SXDM or PIE. This drastically decreases computing time while allowing for further refinement of the solution.

CONCLUSIONS We have applied a fast, deterministic ptychographic reconstruction algorithm, dubbed GPILRUFT, to X-ray diffraction data. It can be applied to the same data as nonlinear ptychographic methods; they all share the advantage that the longitudinal and transverse phase diversity of the probe positions allows for a better reconstruction for the same dose of illumination, in comparison with single-shot retrievals. However, unlike nonlinear algorithms, there is no need to seed GPILRUFT with different starting points to check the uniqueness of the solution, or to average retrieved transmission functions to obtain the final solution. Because each run solves a set of linear equations, which is equivalent to minimizing a parabolic function, there is a single minimum and stagnation in a local minimum does not arise. We have demonstrated GPILRUFT on experimental X-ray data, by reconstructing the transmission function of an Au–Cr nanostructured object. We corrected for the nonlinear contribution to the autocorrelation of the exit surface waves to obtain a result quantitatively comparable with results from nonlinear ptychographic reconstruction algorithms, using an order of magnitude less computing time. We demonstrated that reconstruction artifacts were due to discrepancies between the assumed and actual illumination. Overall, GPILRUFT is shown to be a fast, deterministic ptychographic phase retrieval algorithm that works on experimental X-ray data.

ACKNOWLEDGMENTS This research was supported under the Australian Research Councils Discovery Projects funding scheme (Project No. DP1096025) and its DECRA funding scheme (Project No. DE130100739).

REFERENCES BAO, P., ZHANG, F., PEDRINI, G. & OSTEN, W. (2008). Phase retrieval using multiple illumination wavelengths. Opt Lett 33, 309–311.

Fast Deterministic Ptychographic Imaging Using X-Rays CHAPMAN, H.N., BARTY, A., BOGAN, M.J., BOUTET, S., FRANK, M., HAU-RIEGE, S.P., MARCHESINI, S., WOODS, B.W., BAJT, S., BENNER, W.H., LONDON, R.A., PLÖNJES, E., KUHLMANN, M., TREUSCH, R., DÜSTERER, S., TSCHENTSCHER, T., SCHNEIDER, J.R., SPILLER, E., MÖLLER, T., BOSTEDT, C., HOENER, M., SHAPIRO, D.A., HODGSON, K.O., VAN DER SPOEL, D., BURMEISTER, F., BERGH, M., CALEMAN, C., HULDT, G., SEIBERT, M.M., MAIA, F.R.N.C., LEE, R.W., SZÖKE, A., TIMNEANU, N. & HAJDU, J. (2006). Femtosecond diffractive imaging with a soft-X-ray free-electron laser. Nat Phys 2, 839–843. D’ALFONSO, A.J., MORGAN, A.J., MARTIN, A.V., QUINEY, H.M. & ALLEN, L.J. (2012). Fast deterministic approach to exit-wave reconstruction. Phys Rev A 85, 013816. D’ALFONSO, A.J., MORGAN, A.J., YAN, A.W.C., WANG, P., SAWADA, H., KIRKLAND, A.I. & ALLEN, L.J. (2013). Deterministic electron ptychography at atomic resolution. Phys Rev B 89, 064101. ELSER, V. (2003). Phase retrieval by iterated projections. J Opt Soc Am A 20, 40–55. FIENUP, J.R. (1978). Reconstruction of an object from the modulus of its Fourier transform. Opt Lett 3, 27–29. GUREYEV, T.E., ROBERTS, A. & NUGENT, K.A. (1995). Phase retrieval with the transport-of-intensity equation: Matrix solution with use of Zernike polynomials. J Opt Soc Am A 12, 1932–1942. HANSEN, P.C. (1992). Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 34, 561–580. HANSEN, P.C. (2010). Discrete Inverse Problems: Insight and Algorithms. Philadelphia: SIAM. HAWKES, P.W. (2009). Aberration correction past and present. Phil Trans R Soc A 367, 3637–3664. HENDERSON, R. (1995). The potential and limitations of neutrons, electrons and X-rays for atomic resolution microscopy of unstained biological molecules. Q Rev Biophys 28, 171–193. HESTENES, M.R. & STIEFEL, E. (1952). Methods of conjugate gradients for solving linear systems. J Res Nat Bur Stand 49, 409–436. HOPPE, W. (1969). Diffraction in inhomogeneous primary wave fields.1. Principle of phase determination from electron diffraction interference. Acta Crystallogr A 25, 495–501. LEITH, E. & UPATNIEKS, J. (1962). Reconstructed wavefronts and communication theory. J Opt Soc Am 52, 1123–1128. MAIDEN, A.M. & RODENBURG, J.M. (2009). An improved ptychographical phase retrieval algorithm for diffractive imaging. Ultramicroscopy 109, 1256–1262. MARTIN, A.V. & ALLEN, L.J. (2008). Direct retrieval of a complex wave from its diffraction pattern. Opt Commun 281, 5114–5121.

1099

MIAO, J., CHARALAMBOUS, P., KIRZ, J. & SAYRE, D. (1999). Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens. Nature 400, 342–344. MORGAN, A.J., D’ALFONSO, A.J., WANG, P., SAWADA, H., KIRKLAND, A.I. & ALLEN, L.J. (2013). Fast deterministic single-exposure coherent diffractive imaging at sub-Ångström resolution. Phys Rev B 87, 094115. OSTERMEIER, C. & MICHEL, H. (1997). Crystallization of membrane proteins. Curr Opin Struct Biol 7, 697–701. PUTKUNZ, C.T., CLARK, J.N., VINE, D.J., WILLIAMS, G.J., PFEIFER, M.A., BALAUR, E., MCNULTY, I., NUGENT, K.A. & PEELE, A.G. (2011). Phase-diverse coherent diffractive imaging: High sensitivity with low dose. Phys Rev Lett 106, 013903. QUINEY, H.M., PEELE, A.G., CAI, Z., PATERSON, D. & NUGENT, K.A. (2006). Diffractive imaging of highly focused X-ray fields. Nat Phys 2, 101–104. RODENBURG, J.M. & BATES, R.H.T. (1992). The theory of superresolution electron microscopy via Wigner-distribution deconvolution. Phil Trans R Soc Lond A 339, 521–553. RODENBURG, J.M. & FAULKNER, H.M.L. (2004). A phase retrieval algorithm for shifting illumination. Appl Phys Lett 85, 4795–4797. SCHERZER, O. (1936). Some defects of electron lenses. Zeitschrift für Physik 101, 593–603. SHAPIRO, D., THIBAULT, P., BEETZ, T., ELSER, V., HOWELLS, M., JACOBSEN, C., KIRZ, J., LIMA, E., MIAO, H., NEIMAN, A.M. & SAYRE, D. (2005). Biological imaging by soft X-ray diffraction microscopy. Proc Natl Acad Sci USA 102, 15343–15346. STROKE, G.W., BRUMM, D. & FUNKHOUSER, A. (1965). Threedimensional holography with “lensless” Fourier-transform holograms and coarse P/N polaroid film. J Opt Soc Am 55, 1327–1328. THIBAULT, P., DIEROLF, M., MENZEL, A., BUNK, O., DAVID, C. & PFEIFFER, F. (2008). High-resolution scanning X-ray diffraction microscopy. Science 321, 379–382. UHLEMANN, S., MÜLLER, H., HARTEL, P., ZACH, J. & HAIDER, M. (2013). Thermal magnetic field noise limits resolution in transmission electron microscopy. Phys Rev Lett 111, 046101. WIENER, N. (1930). Generalized harmonic analysis. Acta Math 55, 117–258. ZHANG, F., PEDRINI, G. & OSTEN, W. (2007). Phase retrieval of arbitrary complex-valued fields through aperture-plane modulation. Phys Rev A 75, 043805.