X-ray ptychography with extended depth of field

Esther H. R. Tsai; Ivan Usov; Ana Diaz; Andreas Menzel; Manuel Guizar-Sicairos

doi:10.1364/OE.24.029089

1. Introduction

The depth of field (DOF) refers to the distance by which a sample can be translated in the direction of beam propagation while remaining in the focus of an imaging system. It limits the regions within which high-resolution diffraction-limited imaging can be achieved. To extend the DOF, focus stacking has been used in photography and optical microscopy by merging images taken at different focus through image processing techniques [1]. Application of this method to scanning transmission electron microscopy has also been demonstrated [2]. Light field cameras [3] utilize microlenses placed in front of the image sensor to obtain directional light ray information, allowing the image to be refocused numerically from a single light field image. Together with deconvolution, the method is also applied to optical microscopy to form tomograms [4]. Similarly, for coherent imaging reconstructed wavefields, for example measured by holography, can be propagated to different planes to form a composite image with large DOF using focus stacking. A different strategy that involves modification of the optical system is wave-front coding [5]. The method extends the DOF of an incoherent optical system by using a phase mask to create a modulation transfer function that is insensitive to defocus and subsequently deconvolving the measurement to recover a sharp image of the sample.

Complementary to extending the DOF, some methods remove from the image light arising from out-of-focus planes. These methods allow for sectioning, i.e. separating a specific axial slice, and therefore enable 3D imaging. Confocal microscopy [6,7] achieves optical sectioning by scanning a coherent point illumination at an axial slice of the sample that is in focus while an aperture is placed at the confocal plane before the detector to allow only in-focus signals to be captured. With conventional microscopy, an illumination with structured light can also be used to give a sectioning effectiveness that depends on the frequency spectrum of the imaged structure [8].

The depth-of-field limit also exists in ptychography, a method that is on a fast route to being an established X-ray microscopy technique due to its ability to overcome limitations on lens fabrication by computational image reconstruction. In ptychography [9, 10], the sample is illuminated with a spatially confined coherent beam at overlapping regions while for each illumination a diffraction pattern is recorded. These patterns are used in an iterative phase retrieval algorithm to form a reconstruction with image resolution that can be orders of magnitude better than the scanning step size and the diameter of the illumination beam. Its extension to 3D imaging, ptychographic X-ray computed tomography (PXCT) [11], offers quantitative contrast [12] and is now routinely used for characterizing a wide variety of materials. There has also been great progress in developing instrumentation that gives an isotropic 3D image resolution better than 20 nm [13]. However, the limited DOF, which is ubiquitous to any imaging technique due to diffraction and multiple scattering effects, leads to a compromise between the sample thickness and the achievable resolution. This compromise, for example, hinders using PXCT for the 3D characterization of hierarchical structures in life and materials science where a representative volume needs to be imaged with high resolution.

For ptychography, the DOF limitation can be lifted by incorporating the simultaneous reconstruction of multiple axial slices of the object along the propagation direction. In this manner, a ptychographic reconstruction algorithm addresses propagation effects including multiple scattering and diffraction within the sample. This multi-slice solution enables ptychographic imaging beyond the DOF limit using a conventional 2D ptychography data set. The ptychographical iterative engine (PIE) [14] was generalized for this problem and called 3PIE [15]. The principle has been demonstrated in the optical regime for 3D sectioning without sample rotation [16] and in the X-ray wavelengths with proof-of-principle examples [17,18]. Due to the large numerical aperture feasible in the optical wavelengths regime, it is possible to achieve an almost isotropic resolution, i.e. an axial section sampling that is comparable to the transverse resolution. However, for X-ray wavelengths with photon energies exceeding a few keV, the weak wave-matter interaction results in small diffraction angles at which the data can be practically collected. This leads to a transverse resolution that can be two orders of magnitude worse than the wavelength and, correspondingly, an axial resolution that is also much worse than the transverse resolution. The low axial resolution is insufficient to directly provide 3D X-ray imaging. However, the larger depth of field compared to that for the optical wavelengths also means that, for X-ray wavelengths, only a few object slices in the multi-slice method are required to achieve high-resolution projections, which is the main interest of this work.

Here, we study in detail the dependence of the reconstruction quality versus object thickness for the case of conventional single-slice ptychography, and discuss the popularly recognized theoretical DOF limit. We introduce multi-slice reconstruction algorithms, which are based on iterated projections and maximum likelihood, and provide numerical and experimental demonstrations. We further quantify the benefits gained by multi-slice reconstructions, giving insight into when and how the method should be applied.

2. Depth of field in ptychography

The DOF limit arises in ptychography due to the assumption that the exit X-ray wavefield can be expressed as the product of a 2D object and the incident illumination, given by

ψ_{j, r} = 𝒫_{r - r_{j}} O_{r},

where ψ_j,r is the exit wavefield for the j-th scanning position, r is the 2D set of transverse coordinates (x, y), P_{r−r_j} is the complex-valued incident wavefield centered at position r_j, and O_r is the object transmissivity. In this assumption, the object is presumed to be well represented by a 2D function, usually the integral of the object transmissivity along the propagation direction. This simplification breaks down when the illuminating beam changes significantly during propagation through the object [19] or when diffraction or multiple scattering effects are non-negligible. The limit was considered in the supplementary material of [19], where it was found that

T ≪ a (δ r) / λ,

where T is the sample thickness in the beam propagation direction, a is the focus size of the illuminating beam, δr is the image resolution, and λ is the illumination wavelength. It is shown that Eq. (2) is the sufficient condition under which the object and the incident illumination are independent of each other while satisfying the inhomogeneous wave equation [19]. The sample thickness limit was also analyzed in terms of Ewald sphere constructions for ptychography [20], arriving at the expression that the sample thickness should satisfy

T \leq 2 {(δ r)}^{2} / λ

if the Ewald sphere departure is less than the speckle width 1/T. A similar expression was obtained in the context of coherent diffraction imaging [21], using the stricter condition of an Ewald sphere departure of 1/(4T), giving then T ≤ (δr)²/(2λ). As the effect is mainly driven by diffraction, it is not surprising that this expression and Eq. (3) are similar to that for the DOF of conventional imaging or for the depth of focus of a lens, which is given by 4(δr)²/λ. We have tested Eq. (3) through a numerical study, described in the supplementary information of [13], and suggested a minor refinement, given by

T \leq 5.2 {(δ r)}^{2} / λ .

The tomogram in [13] has a 16 nm resolution on a sample of 6 μm thickness measured at 6.2 keV photon energy, which comes very close to this limit, T ≤ 6.7 μm. Equations (2)–(4) are all inversely proportional to the wavelength, showing that the thickness limit is even stricter in the soft X-ray regime. At 710 eV, a nanoplate with 5 nm structures has been measured [22] and, according to Eq. (4), the limit for 5 nm resolution is T ≤ 74 nm. If either Eq. (2) or Eq. (4) is not satisfied, the resolution of the reconstruction could be affected by propagation or multiple scattering effects.

The idea behind multi-slice propagation is to represent the object by a finite number of axial slices that individually satisfy the multiplicative approximation in Eq. (1) [15]. In the beam propagation method [23, 24], the incident beam is multiplied by the integrated transmissivity of the first slice and then propagated to the next object slice assuming free-space propagation. Figure 1 shows the object slices, $O_{r}^{(1)}$ to $O_{r}^{(N)}$ , in the beam propagation direction, i.e. the positive z direction, with N the total number of slices and Δz_n,n+1 the separation between slices n and n + 1. As a result, the exit wavefield is given by

ψ_{j, r}^{(N)} = 𝒫_{Δ z_{N - 1, N}} {𝒫_{Δ z_{N - 2, N - 1}} {\dots 𝒫_{Δ z_{2, 3}} {𝒫_{Δ z_{1, 2}} {P_{r - r_{j}} O_{r}^{(1)}} O_{r}^{(2)}} O_{r}^{(3)} \dots} O_{r}^{(N - 1)}} O_{r}^{(N)},

where 𝒫_Δz {·} is the free-space propagation over distance Δz in the z direction. In this manuscript, we refer to the virtual discretization of objects in Fig. 1 as object ’slices’ and the separation between the slices ’slice separation’, while ’thickness’ refers to the entire extent of the sample in the beam propagation direction comprising all the object slices.

Fig. 1 Illustration of multi-slice ptychography, which reconstructs multiple axial slices of the object along the propagation direction to account for propagation effects including multiple scattering, diffraction, and propagation within the sample.

Download Full Size | PDF

The illumination incident on the first object slice, P_{r−r_j}, is common for all positions. The illumination incident on the n-th slice is in general affected by the previous slice transmissivity and is denoted by $P_{r}^{(n)}$ . The wavefield after interaction with the n-th object slice is given by $ψ_{j, r}^{(n)}$ , such that $ψ_{j, r}^{(n)} = P_{j, r}^{(n)} O_{r}^{(n)}$ . The propagation between slices results in the relation $P_{j, r}^{(n + 1)} = 𝒫_{Δ z_{n, n + 1}} {ψ_{j, r}^{(n)}}$ . After propagation through the object, the far-field diffraction patterns are given by

I_{j, q} = {| {\tilde{ψ}}_{j, q}^{(N)} |}^{2} = {| ℱ {ψ_{j, r}^{(N)}} |}^{2},

where ℱ {·} is a 2D Fourier transform, q = (q_x, q_y) is the reciprocal space Cartesian coordinate vector, and

{\tilde{ψ}}_{j, q}^{(N)}

is the Fourier transform of

ψ_{j, r}^{(N)}

.

3. Reconstruction algorithms

The general interaction model between the object and illumination in Eq. (5) requires the reconstruction of all the object slices, $O_{r}^{(n)}$ , and the incident illumination, P_r. This reconstruction was first incorporated into ptychography in [15] with the introduction of the 3PIE algorithm.

For conventional ptychographic reconstructions with single-slice, projection-based algorithms are used for their robustness and fast convergence [19], followed by a maximum-likelihood refinement that is used to provide the final solution accounting for the noise statistics [25,26]. Motivated by this, we introduce here a projection-based ptychography algorithm as well as an optimization-based algorithm that incorporates the multi-slice method. These algorithms use all measurements simultaneously for updating the object slices and the illumination, in contrast to 3PIE where each diffraction pattern is used sequentially for updating the reconstruction.

3.1. Projection-based method

Algorithms designed to find the intersection of two projection sets and robust to stagnation have been developed for phase retrieval and ptychography [27–29]. Here the sets of constraints and their associated projection operators are generalized to incorporate the multi-slice method. As customary for phase retrieval, the Fourier projection is obtained by replacing the magnitudes of the exit view Fourier transform by $\sqrt{I_{j, q}}$ while keeping their phase unchanged, namely

Π_{F} {ψ_{j, r}^{(N)}} : ψ_{j, r}^{(N)} \to ψ_{j, r}^{(N) F} = ℱ^{- 1} (\frac{{\tilde{ψ}}_{j, q}^{(N)} \sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) .

For conventional ptychography, the object space overlap constraint is satisfied when the smallest modification is made to the object and the illumination to ensure that the exit waves can be factored as Eq. (1) by minimizing an error function [29]. In the case of multi-slice ptychography, the overlap projection must be extended so that this factorization between object and illumination is possible at all slices of the object, resulting in the form

Π_{O} {ψ_{j, r}^{(n)}} : ψ_{j, r}^{(n)} \to ψ_{j, r}^{(n) O} = {\hat{P}}_{j, r}^{(n)} {\hat{O}}_{r}^{(n)}, \forall n,

where

P_{j, r}^{(1)} = P_{r - r_{j}}

. Extending the derivation in [29] to the multi-slice case, the projection can be calculated by minimizing the error function

ℰ = \sum_{n} \sum_{j} \sum_{m_{r}} {| ψ_{j, r}^{n} - {\hat{P}}_{j, r}^{(n)} {\hat{O}}_{r}^{(n)} |}^{2},

where indices m_r = (m_x, m_y) = (x/Δ_x, y/Δ_y) and Δ_x and Δ_y are the discretization sizes along the x and y directions, respectively. Calculating the partial derivatives of Eq. (9) with respect to

O_{r}^{(n)}

and

P_{j, r}^{(n)}

and setting them equal to zero give the projection in the form of a system of equations,

{\hat{P}}_{r} = \frac{\sum_{j} O_{r + r_{j}}^{(1) *} ψ_{j, r}^{(1)}}{\sum_{j} {| O_{r + r_{j}}^{(1) *} |}^{2}}, for n = 1,

{\hat{P}}_{j, r}^{(n)} = \frac{O_{r}^{(n) *} ψ_{j, r}^{(n)}}{{| O_{r}^{(n)} |}^{2}}, for 1 \leq n \leq N,

{\hat{O}}_{r}^{(n)} = \frac{\sum_{j} P_{j, r}^{(n) *} ψ_{j, r}^{(n)}}{\sum_{j} {| P_{j, r}^{(n)} |}^{2}}, \forall n,

where (*) denotes the complex conjugate. For free-space propagation between slices, the angular spectrum method is used, namely

𝒫_{Δ_{z}} {ψ_{r}} = ℱ^{- 1} {ℱ {ψ_{r}} exp (i k Δ z \sqrt{1 - {| λ q |}^{2}})},

where ψ_r is an arbitrary wavefield and k = 2π/λ is the wave number. This system of equations given in Eq. (10) can be solved iteratively in a manner analogous to [29]. In this work, one inner iteration within Eq. (10) was used for computation simplicity as all exit waves need to be propagated through all object slices for each additional iteration, which adds significantly to the computation time. In other words, instead of a simultaneous update, the object slices were updated sequentially, for which Eqs. (10b) and (10c) were applied as an update step starting from the exit wave at n = N and ending with n = 1 at which point the incident illumination was updated with Eq. (10a). As the reconstruction progresses and better estimates of the object and illumination become available, this simplification approaches the solution to Eq. (10). The update of multiple object slices and the Fourier transforms in the propagation leads to an increase in the computation time that scales approximately linearly with the number of slices.

Based on the conventional algorithm, difference map (DM) [29], we introduce the multi-slice ptychography algorithm described above as 3D difference map (3DM). The views for the k-th iteration are therefore updated as

ψ_{j, r}^{(N) [k + 1]} = ψ_{j, r}^{(N) [k]} + Π_{F} {2 Π_{O} {ψ_{j, r}^{(N) [k]}} - ψ_{j, r}^{(N) [k]}} - Π_{O} {ψ_{j, r}^{(N) [k]}} .

3.2. Maximum Likelihood

We also introduce a multi-slice gradient-based optimization method, which we call 3D maximum-likelihood (3ML). Generalizing the method in [25,26] for multi-slice ptychography, we define an error metric that approximates to the first order the negative log-likelihood for Poisson distributed noise [26]. The error metric is given by

ℒ = \sum_{j} \sum_{m_{q}} {(| ℱ {P_{j, r}^{(N)} O_{r}^{(N)}} | - \sqrt{I_{j, q}})}^{2},

where indices m_q = (m_{q_x}, m_{q_y}) = (q_x/Δ_{q_x}, q_y/Δ_{q_y}) and Δ_{q_x} and Δ_{q_y} are the discretization sizes in the reciprocal space. The goal is to iteratively minimize Eq. (13) by providing improved estimates of the object slices and the incident illumination. Due to the large number of parameters to optimize, an analytical expression for the gradient is necessary for the optimization routine. Using Wirtinger derivatives [30,31], the gradient of the error metric with respect to the parameters to be reconstructed is given by

\frac{\partial ℒ}{\partial O_{r}^{(n)}} = {\begin{array}{l} \sum_{j} P_{j, r}^{(n)} {(𝒫_{- Δ z_{n, n + 1}} {\frac{\partial ℒ}{\partial P_{j, r}^{(n + 1) *}}})}^{*} & ; 1 \leq n \leq N \\ \sum_{j} P_{j, r}^{(N)} χ_{j, r}^{(N) *} & ; n = N \end{array}

\frac{\partial ℒ}{\partial P_{j, r}^{(n)}} = {\begin{array}{l} O_{r}^{(n)} {(𝒫_{- Δ z_{n, n + 1}} {\frac{\partial ℒ}{\partial P_{j, r}^{(n + 1) *}}})}^{*} & ; 1 \leq n \leq N \\ O_{r}^{(N)} χ_{j, r}^{(N) *} & ; n = N \end{array}

χ_{j, r}^{(N)} = ℱ^{- 1} {(1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) {\tilde{P}}_{j, q}^{(N)}} .

To compute the gradients, first the gradients at slice N are calculated and subsequently the gradients for slice n to n − 1 are calculated recursively. Finally for the incident illumination, P_r, the gradient is expressed by

\frac{\partial ℒ}{\partial P_{r}} = \sum_{j} \frac{\partial ℒ}{\partial P_{j, r + r_{j}}^{(1)}} .

Optimization algorithms are well suited to refine other experimental uncertainties. Additional to the reconstruction of the object and illumination, the slice separation can be retrieved. The gradient with respect to the slice separation is given by

\begin{array}{l} \frac{\partial ℒ}{\partial Δ z_{n, n + 1}} = & 2 \sum_{j} \sum_{m_{q}} (1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \\ \times ℱ {𝒫_{Δ z_{N - 1, N}} {\dots 𝒫_{Δ z_{n + 1, n + 2}} {ℱ^{- 1} {ℱ {P_{j, r}^{(n)} O_{r}^{(n)}} \\ \times i k \sqrt{1 - {| λ q |}^{2}} exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}})} \\ \times O_{r}^{(n + 1)}} O_{r}^{(n + 2)} \dots} O_{r}^{(N)}}}, \end{array}

where ℛ{·} gives the real part of a complex-valued array. Derivations for Eqs. (14)–(16) and Eq. (18) are given in the Appendix. Compared to 3PIE and 3DM, the gradient-based optimization has the advantage of using all data for updating all slices and all variables simultaneously and is therefore not affected by the choices on the order of slice update. Moreover, the method allows the incorporation of an appropriate noise model, provides a natural stopping criterion as the algorithm converges, and offers the unique capability to refine the axial position of the slice. The slice position refinement could be helpful when the default definition of equally spaced slices is not optimal for the representation of the object, e.g. when there are strongly scattering features sparsely distributed with unknown positions inside a sample.

4. Experiment

To demonstrate and validate these algorithms, we performed experiments at the cSAXS beamline (X12SA), Swiss Light Source, Paul Scherrer Institut, Switzerland. For the X-ray illumination, monochromatic radiation of 6.2 keV photon energy, i.e. wavelength λ ≈ 0.2 nm, was defined with a double crystal Si(111) monochromator. A Fresnel zone plate (FZP) made of Au [32] with an outer-most zone width of d = 60 nm and a diameter of D = 250 μm was used, giving a focal length of f = Dd/λ = 75 mm and a depth of focus of 4d²/λ = 72 μm. A central beam stop was placed upstream the FZP and an order sorting aperture with 20 μm diameter was mounted downstream near the focus. Both the FZP and the central stop were fabricated in the Laboratory for Micro and Nanotechnology, Paul Scherrer Institut, Switzerland. The sample consisted of 2 μm-diameter SiO₂ microspheres dispersed on two layers of Si₃N₄ membranes, which were separated along the direction of beam propagation by a nominal layer separation of 200 μm. The sample was placed 1.5 mm downstream of the focus, where the illumination had about 4 μm diameter, and scanned over a 20×20 μm² field of view (FOV) following a Fermat spiral pattern [33] with a 1 μm average step size, giving 400 scanning points. At each scanning point, a diffraction pattern was recorded 7.2 m downstream of the object using a Pilatus 2M detector [34] with 172×172 μm² pixel size and 0.15 s exposure time, corresponding to 5.1×10⁷ photons per pattern. With a diffraction pattern size of 192×192 pixels, the pixel size of the reconstructed object was 43.9 nm. According to Eq. (4), for a thickness of 200 μm, the resolution would be limited to 88 nm.

For comparison with conventional methods, the projected object transmissivity is obtained from the reconstructed slices by a product, namely

O_{r} = \prod_{n = 1}^{N} O_{r}^{(n)} .

The reconstruction of the microspheres using the conventional single-slice algorithm, 3DM, and 3ML, are shown in Figs. 2(a)–2(c), respectively. The multi-slice images in Figs. 2(b) and 2(c) are the phase of projected object, i.e. O_r in Eq. (19). For the multi-slice algorithms, reconstructions were performed with N = 2 and with a fixed Δz_1,2 of 200 μm. Figures 2(d)–(g) show the reconstruction of individual slices,

O_{r}^{(n)}

. Figures 2(h) and 2(i) show the comparison of conventional ptychography and 3ML. The microspheres on membrane layer 2 are blurry in Fig. 2(h) but are clearly resolved in 2(i). Moreover, the microspheres on layer 1 have dark spot artifacts in the center, indicated by the pink arrowheads in Fig. 2(h), but not in 2(i). Besides the microspheres, impurities deposited during the sample preparation are clearly reconstructed on slice 2 in Figs. 2(d) and 2(e), indicated by the yellow arrowheads.

Fig. 2 Experimental reconstructions of microspheres on membranes at 6.2 keV photon energy. Comparison of (a) conventional difference map algorithm, (b) 3DM, and (c) 3ML. In (b) and (c), individual slices are combined according to Eq. (19). Reconstructions of individual slices for 3DM are shown in (d) and (f), and for 3ML in (e) and (g). Yellow arrowheads point to impurities in the sample, clearly reconstructed on slice 2 only. Insets of DM and 3ML are shown in (h) and (i), respectively. Pink arrowheads point to particular artifacts of the central area of particles in the DM reconstruction.

Download Full Size | PDF

To quantify the reconstruction quality, we use Fourier ring correlation (FRC) [35,36]. FRC is a commonly used technique for determining image resolution by computing the normalized cross-correlation coefficient between two independently acquired images at a given thin shell in reciprocal space, giving a measure of the consistency of the image signal as a function of spatial frequency that monotonically decreases with decreasing signal-to-noise ratio. We determine the resolution by the intersection of the FRC curve and the threshold. However, from two independent measurements, the FRC computed for images reconstructed with the conventional algorithm provided high resolution estimations of tens of nm that were not consistent with the observed image quality, shown by the low-resolution reconstruction of the 2 μm microspheres in Figs. 2(a) and 2(h). We attribute this to the fact that the artifacts introduced by exceeding the ptychography DOF limit are systematic, producing highly reproducible features that boost the FRC curve in a manner that is not representative of the image sharpness, rendering the FRC resolution estimation unsuitable in this case. To circumvent this issue, we used the fact that the sample is known and modelled portions of the image in Figs. 2(a) and 2(c) with the transmissivity of a perfect sphere convolved with a Gaussian point spread function. A similar approach was used in [37] where the model was based on a series of perfect spheres whose sizes and positions were deduced from a higher resolution scanning electron microscope image. This model was convolved with a Gaussian kernel to quantitatively match the reconstruction. The only fitting parameter in [37] was the Gaussian full-width-half-maximum (FWHM), which was also their resolution estimate.

In this work, for a microsphere in the reconstructed projection, a non-linear optimization routine was used to minimize an error metric that compares a fitted model with the reconstructed microsphere. The model parameters that were optimized were the radius and refractive index of the microsphere, its transverse position, the global phase offset and the linear phase terms, which are inherent degrees of freedom for ptychography [38], and the standard deviation of the Gaussian kernel. The FRC between the reconstructed image and the fitted model image were calculated after subpixel image registration [39] and tapering of the image edges. The half-period resolution was then estimated with the 1-bit FRC threshold curve [35], which corresponded well to the Gaussian FWHM, giving confidence in the resolution estimates. This technique allows us to fairly assess and compare the resolution for both the conventional method and the multi-slice methods. The conventional method in Fig. 2(a) gives a resolution estimate of 90 nm for the spheres in focus and 290 nm for the others. With 3ML, Fig. 2(c) gives a resolution estimate of 75 nm for spheres on both layers, showing a significant improvement in resolution compared to the conventional method as well as consistency in the resolution at different axial slices.

Inspection of the slice reconstructions, e.g. slice 2 in Figs. 2(f) and 2(g), reveals faint dark shadows of the microspheres on slice 1 seeping through and appearing on the reconstruction of slice 2, and vice versa. These shadows, together with the bright halos that can be seen around individual microspheres in Figs. 2(d)–(g), can be explained by the dependence of multi-slice effectiveness on spatial frequencies, which will be discussed in Section 5.1. As shown in Figs. 2(b) and 2(c), the artifacts were reduced when the slices were combined for the purpose of generating high-resolution projections.

A unique feature of the gradient-based optimization is the ab initio determination of the slice separation, Δz, during the reconstruction. Figures 3(a) and 3(c) show the slices reconstructed using 3ML with a fixed wrong separation of 100 μm. The wrong slice separation prohibits the algorithm to separate the slices without artifacts on both slices.

Fig. 3 Reconstruction using 3ML with fixed wrong slice separation, shown by (a), (c), and (e), and with separation determination, shown by (b), (d), and (f). (a)–(b) Reconstructed phase for slice 1. (c)–(d) Reconstructed phase for slice 2. (e)–(f) The slice separation, Δz, over the iterations. The separation was fixed at 100 μm in (e) while it converged to 220 μm in (f).

Download Full Size | PDF

To optimize the separation, Δz was fixed at an initial guess of 100 μm in the first 30 iterations while only the illumination and the object slices were updated. This number of iterations was chosen to allow the algorithm to arrive at a rough reconstruction before starting the optimization of Δz. The separation Δz was then updated together with the illumination and the object slices for the following 30 iterations, after which the distance was fixed for another 30 iterations to allow the algorithm to focus on improving the estimate of object and illumination based on the latest estimate of the separation. Subsequently this sequence of alternating between 30 iterations of updating all variables and 30 iterations of no separation update was repeated. We do not expect the overall reconstruction to be sensitive to a particular choice of the number of iterations in this cycle, although in extreme cases it may affect the overall convergence speed. The slice separation as a function of iteration is shown in Fig. 3(f). After 3000 iterations, the separation converged to 220 μm and the reconstructed slices, shown in Figs. 3(b) and 3(d), exhibited improved image quality with significantly reduced artifacts compared to Figs. 3(a) and 3(c), respectively. For a wide range of initial guesses for the slice separation, e.g. 1 μm and 300 μm, the algorithm converged to separations between 219 μm and 222 μm, in around 9000 iterations and 3000 iterations, respectively. Measurement of the membrane layer separation with an optical microscope later verified a separation of 220±3 μm, which closely matches the Δz obtained by our algorithm.

5. Numerical studies

We have validated experimentally in Section 4 both the projection-based method, 3DM, and the optimization-based method, 3ML. In this section, we use numerical simulations to quantify the benefits, constraints, and boundary conditions for the practical application of multi-slice ptychography. Additionally, we study the ability to achieve sectioning and the effectiveness of our methods in terms of improving transverse resolution by extending the DOF.

For incident illumination, we simulated a spatially and temporally coherent beam focused by a FZP with 170 μm diameter and 50 mm focal length, giving an outer zone width of 59 nm at λ = 0.2 nm and a depth of focus of 70 μm. The beam was propagated numerically 1 mm downstream from the focus, giving a diameter of 3.4 μm incident on the object. The object consisted of three layers, each taken from a different mouse histology image [40,41], as shown by the first column of Fig. 4. To simulate a realistic X-ray scattering, the greyscale images were normalized and converted to a height map with thickness ranging from zero to 1 μm, and assigned a refractive index (1 − δ + iβ) with δ = 1.19 × 10⁻⁵ and β = 3.36 × 10⁻⁸. These values correspond to the index of refraction of carbon at 6.2 keV photon energy and provide each slice transmissivity with a maximum phase shift of −0.12π rad.

Fig. 4 Ptychography simulation at 6.2 keV, imaging objects with three layers and thicknesses of 40 μm, 200 μm, and 2 mm. With increasing layer separation, the single-slice method gives worse reconstructions, shown by the first row, while the multi-slice method gives visibly a constant reconstruction quality, shown by the projections in the second row. The light blue lines indicate the FRC curves for the case of 40 μm-thickness, the blue for 200 μm, dark blue for 2 mm, and the red lines give the 1-bit threshold. Rows three to five show the individual slices. The sectioning effectiveness improves with increasing separation and, with small separation, low spatial frequency features cannot be effectively sectioned, as shown by the reconstructed images and pointed by the yellow arrowheads at the FRC curves.

Download Full Size | PDF

The object was scanned over a FOV of 20×20 μm² on concentric rings with 5n scanning positions on the nth ring [42] and a scanning step s = 1.5 μm. The detector was 7.2 m downstream of the object with pixel size 172×172 μm². A diffraction pattern of size 512×512 was generated at each scanning point, with added Poisson noise corresponding to an average of 10⁸ photons per diffraction pattern, realistic in synchrotron experiments. We have considered these high signal-to-noise conditions in order to focus on the influence of object thickness on the image quality. The pixel size of all reconstructed objects in the numerical studies was 16.4 nm, giving a thickness limit of 7 μm according to Eq. (4) if the resolution equals the pixel size. In this section, the FRC with the 1-bit criterion [35] was used for resolution estimation based on the true image of the simulation that is available and the reconstructed image.

5.1. Image sectioning

For tomography applications, we are primarily interested in using multi-slice ptychography to obtain high-resolution projections. However, in some applications, e.g. when multi-slice methods are used instead of tomography to achieve 3D imaging in the optical regime, the axial information is of interest, hence we also study the performance of axial sectioning. In Fig. 4, we show simulations with three equally-spaced object layers with a total object thickness of 40 μm, 200 μm, and 2 mm. To probe exclusively the effect of sample thickness and disentangle it from a possible enhanced image quality due to larger integrated phase shift, in these simulations we kept the total phase shift introduced by the sample constant by using the same object layers and only varied the separation to simulate different sample thicknesses. Reconstructions were performed using 3ML with N = 3 and with the true separations as prior knowledge, i.e. Δz_1,2 = Δz_2,3 = 20 μm in the simulation for a thickness of 40 μm, and similarly for thicknesses of 200 μm and 2 mm. The FRC curve for the case of 40 μm-thickness is shown in Fig. 4 in light blue, the 200 μm case in blue, the 2 mm case in dark blue, and the FRC 1-bit threshold in red. Using maximum likelihood with a single slice, the reconstruction quality degrades with increasing thickness, which can be observed from the reconstructed images on the first row in Fig. 4 and also quantitatively by the FRC curves. On the second row of Fig. 4, the high correlation of the FRC curves in the low frequencies indicates that these frequency components were reconstructed rather accurately, therefore there was almost no visible difference in the combined reconstructed slices, giving confidence in the faithful reconstruction of projections for all thicknesses. Moreover, the resolution determined by FRC between the true image and the reconstructed image showed that better resolution was obtained with larger separation.

Reconstructed images on rows three to five in Fig. 4 clearly show that, with smaller separations, the quality of the individual slices degrades. The reconstructed object slices are affected by features from neighboring slices seeping through, i.e. the sectioning ability of the algorithm is reduced. While the exit wavefield, $ψ_{r}^{(N)}$ , is expressed as

\begin{array}{l} ψ_{r}^{(N)} & = P_{r}^{(N)} O_{r}^{(N)} \\ = 𝒫_{Δ z_{N - 1, N}} {ψ_{r^{'}}^{N - 1}} O_{r}^{(N)} \\ = 𝒫_{Δ z_{N - 1, N}} {P_{r^{'}}^{(N - 1)} O_{r^{'}}^{(N - 1)}} O_{r}^{(N)}, \end{array}

in the limit of a vanishingly small propagation distance it is reduced to

ψ_{r}^{(N)} \approx P_{r}^{(N - 1)} O_{r}^{(N - 1)} O_{r}^{(N)} .

The latter expression shows that features in

O_{r}^{(N - 1)}

and

O_{r}^{(N)}

cannot be separated in this case, resulting in the object slice ambiguity in which features from one object slice can seamlessly be attributed to a different slice during reconstruction.

The dependence of slicing effectiveness versus spatial frequency can be examined in the FRC curves from rows three to five in Fig. 4, which give the sectioning fidelity versus spatial frequency for each separation and each slice. The different true image slices, shown in the first column of Fig. 4, have similar content in the very low frequencies, hence spatial frequencies below 10% of the Nyquist frequency have always high correlation that are not attributed to a correct sectioning. For the case of 40 μm-thickness, the correlation of each slice takes a pronounced dip in the low frequencies, indicated by the yellow arrowheads, while the corresponding FRC of the projection (on the second row) shows a high correlation in these frequencies. This is consistent with the observation that low-spatial frequency features are accurately reconstructed for the projection but these features seep through neighboring slices, i.e. that low spatial frequencies are more difficult to be axially sectioned. This effect is akin to that which is present in optical microscopes and is described by the 3D optical transfer function [43], in which low spatial frequencies require larger slice separations to be effectively sectioned. For multi-slice ptychography in the paraxial, or small-angle, approximation we can consider the coherent paraxial propagation transfer function [44],

H_{q} = exp (i k Δ z) exp (- i π λ Δ z {| q |}^{2}),

which explicitly shows that the effect of free-space propagation is caused by a phase shift that depends on the spatial frequency. If a significant propagation effect, and thus the ability for sectioning, occurs at a phase difference of π with respect to the zero frequency component, it gives

Δ z = \frac{1}{λ {| q |}^{2}},

which suggests an inverse relationship between the minimal sectioning distance and the square of the spatial frequency. Hence, for a given slice separation, sectioning has decreased performance for lower spatial frequencies, which correspond to broad features and to the shadows that seep through to adjacent slices, as also observed in [16]. It should be emphasized that, for tomography, the sample rotation provides a significant redundancy of the low and middle spatial frequencies, and thus effective sectioning of those frequencies is not required. Therefore, we expect this issue to pose no significant hurdle for the method’s application to tomographic synthesis.

5.2. Extended depth of field

For application of multi-slice ptychography to nanotomography, a high resolution projection is sought with the aim to overcome traditional limitations on resolution due to the sample thickness. Therefore, for each multi-slice reconstruction, the reconstructed slices are combined to give a numerical object projection, given by Eq. (19). For practical applications, the conditions under which a multi-slice reconstruction is beneficial are of particular interest. To elucidate this matter, conventional ptychography and multi-slice reconstruction are compared at different sample thicknesses and different scanning steps.

We performed numerical simulations with different scanning steps, ranging from 1.5 μm to 2.0 μm. This gave, e.g. 141 scanning points for a scanning step s = 1.5 μm and 83 points for s = 2.0 μm. Twelve different sample thicknesses were simulated using three equally-spaced object layers, with a thickness that ranged from 5 μm to 200 μm. The projected objects, obtained through Eq. (19), were compared between different algorithms: the conventional method and multi-slice with two and three slices. Figure 5 shows the FRC resolution of the phase of the projected object versus object thickness. The blue curve a₁ in Fig. 5 gives the result for conventional difference map at s = 1.5 μm, a₂ for 3DM reconstructed with N = 2 slices, and a₃ for 3DM reconstructed with N = 3 slices. Similarly for s = 2.0 μm, the results are shown by the red curves, b₁ for single-slice, b₂ for two, and b₃ for three. Reconstructions were performed using the correct slice separations. Similar curves were obtained using 3ML. It should be noted that while the true object has three layers, reconstructions with two slices were also successful, determined by the comparison of the resolutions of projected objects. This shows the robustness of our algorithms to the mismatch between the object layers and the reconstruction slices, which is important in practice when imaging 3D axially contiguous objects.

Fig. 5 Numerical simulation for resolution versus object thickness. With s = 1.5 μm, the blue curve a₁ gives the result for conventional method, a₂ for 3DM reconstructed with two slices, and a₃ for 3DM reconstructed with three slices. Similarly for b₁, b₂, and b₃, with s = 2.0 μm. The intersections of DM and 3DM with two slices indicate the point where multi-slice ptychography becomes beneficial and are indicated by circles: blue circles for s = 1.5 μm, red ones for s = 2.0 μm, and green for intermediate scanning steps. Gray circles are the intersecting points between DM and 3DM with three slices. For simplicity, full curves are not shown for these intermediate scanning step data sets. Equations (3) and (4) are shown as gray and black lines, respectively.

Download Full Size | PDF

Figure 5 shows that by using multi-slice ptychography, the resolution can be improved significantly if the resolution was limited by the object thickness, effectively removing the thickness limit. For example, for the case with a step size s = 1.5 μm, shown in Fig. 5, the resolution is improved from 47 nm for a₁ to 20 nm for a₂ and a₃ on the same data set, i.e. 200 μm-thickness, by simply using a multi-slice reconstruction.

Comparing the conventional reconstructions with different scanning steps, i.e. a₁ and b₁, it is notable that the relation between the resolution and the object thickness is dependent on the ptychography scanning step. While Eq. (4) predicts a resolution limit of 85 nm for Δz = 200 μm, the resolution obtained for s = 2.0 μm is 70 nm and is further improved to 47 nm by merely decreasing the step size to s = 1.5 μm. Projection resolution and depth of field were improved from b₂ to a₂ and also from b₃ to a₃, as shown in Fig. 5. The sectioning was also more effective for the smaller step size. Ptychography is considered to be a well-constrained problem from which even additional multiplexed functions can be extracted, e.g. mutually incoherent modes [45,46]. Although the scanning step is the central sampling criterion for ptychographic data sets [47,48], the concept that this could also be used to alleviate the thickness limited resolution is to our knowledge a novel finding.

The thickness limit, given by Eq. (3) or (4), is widely regarded as a required trade-off between resolution and object thickness that spans across imaging techniques. These predictions are shown respectively as gray and black solid lines in Fig. 5. However, all the curves in Fig. 5, either conventional for a₁ and b₁ or multi-slice for the rest, show a larger depth of field than what is theoretically predicted by Eq. (4). The intersection points between conventional and multi-slice reconstructions, shown as circles in Fig. 5, indicate the scenario in which it becomes beneficial to apply multi-slice ptychography. We notice that the trend of these intersections points follows Eq. (4). For example, the blue circles are around thickness 10 μm for s = 1.5 μm, red ones are around thickness 35 μm for s = 2.0 μm. In other words, although Eq. (4) and its variations are often taken to be the absolute compromise between transverse resolution and sample thickness, numerical simulations show that this limit can be overcome even by a conventional ptychography reconstruction. Moreover, here we show that Eq. (4) instead provides a guideline for the axial sampling criterion for multi-slice ptychography.

These simulations were performed with a fixed average number of photons per diffraction pattern. To rule out the possibility of having the effect of improved photon statistics with the smaller scanning steps, we also performed simulations where the dose on the sample was fixed by having a total of 10¹⁰ photons for the whole scan and obtained similar results to Fig. 5. This is expected as we are performing these simulations in a high SNR regime where the photon statistics are not expected to be the limiting factor for resolution.

One criterion for multi-slice is that the data contains enough information for the algorithm to reconstruct multiple slices. The curve for reconstructions using three slices with s = 1.5 μm, a₃, follows closely the two-slice reconstructions a₂ for the entire thickness range. However, the same curves with s = 2.0 μm, b₂ and b₃, follow each other closely but, for large thicknesses, the reconstruction of the projection worsens for three-slice reconstruction. This result is counter-intuitive as we expect that, as the thickness of the object increases, we should be able to better represent the object with more slices. We attribute this result to the effect of the scanning step on the ptychography data sets. A large step size can result in a data set that fails to sufficiently constraint the reconstruction as the number of parameters to reconstruct increases. In the small object thickness regime, i.e. if the resolution is not limited by the thickness and the object is well-represented by a 2D function, using multi-slice degrades the resolution compared to single-slice methods because the number of degrees of freedom increases linearly with the number of slices and reconstruction becomes poorly constrained due to the object slice ambiguity discussed in Section 5.1. This degraded performance is even more significant with s = 2.0 μm than with s = 1.5 μm, observed at a thickness of 5 μm in Fig. 5.

Using a single virtual slice as an approximation for nominally thick objects works well for ptychography under a compromise between sample thickness and image resolution. With multi-slice, a large object is virtually divided into several of these slices, as shown in Fig. 1. The multi-slice method has already shown to work well in thick continuous samples at optical wavelengths [16], and it should likewise work for X-rays. Yet until now there has been little study on what is the optimal axial sampling criterion for multi-slice ptychography. We show in Fig. 5 that Eq. (4) offers a guideline for such sampling, shedding a light into the subject. Understanding the relation between slice separation, resolution, and the step size for multi-slice reconstructions facilitates the widespread application of the method in the future.

6. Conclusions

We have introduced two novel multi-slice ptychography algorithms, a projection-based method and an optimization method, with the latter providing the unique capability to determine the object thickness by optimizing the slice separation. We have experimentally validated both methods and show the ab initio determination of the thickness. Numerical studies shown in this work reveal that the theoretical DOF prediction in Eq. (4) or its variations do not serve as a resolution or thickness limit even for conventional ptychography. Instead, Eq. (4) gives the condition under which multi-slice becomes beneficial compared to conventional algorithms for the purpose of obtaining an object projection. For a given scanning step size and incident flux, a multi-slice reconstruction pushes the resolution to the limit by removing the thickness limitation. We have discussed the axial sectioning in Section 5.1 and the axial sampling criterion in Section 5.2. We also show that the DOF for ptychography has a dependence on the scanning step. This combined knowledge will be instrumental for the routine implementation of multi-slice reconstructions in X-rays, optical regime, or electrons [49], especially in the latter two, where multiple scattering effects are particularly important. With suitable adjustment of the propagator operators, multi-slice methods are also applicable to the growing field of Fourier ptychography [50, 51]. This work further serves as the first step to applying multi-slice to tomography to allow for high-resolution measurements on volumes larger than what would be possible with the DOF limit. Additionally, the axial information obtained with multi-slice via sectioning can potentially be built into the tomography reconstruction to allow for coarser tomography angular sampling and reduce significantly the measurement time.

Appendix - Derivation for multi-slice optimization-based ptychography

Here we show the derivation of analytical expressions of the gradients of object slices, the incident illumination, and the slice separation as needed for the 3ML algorithm.

For refinement of the object slice, we calculate the gradient of the error metric by taking the partial derivative with respect to the real and imaginary parts of the complex-valued object slices, $O_{r}^{(n)} = O_{r}^{(n) ℛ} + i O_{r}^{(n) ℐ} = ℛ {O_{r}^{(n)}} + i ℐ {O_{r}^{(n)}}$ ,

\frac{\partial ℒ}{\partial O_{r}^{(n) ℛ}} = 2 \sum_{j} \sum_{m_{q}} (| {\tilde{ψ}}_{j, q}^{(N)} | - \sqrt{I_{j, q}}) \frac{\partial | {\tilde{ψ}}_{j, q}^{(N)} |}{\partial O_{r}^{(n) ℛ}},

where

\begin{array}{l} \frac{\partial | {\tilde{ψ}}_{j, q}^{(N)} |}{\partial O_{r}^{(n) ℛ}} & = \frac{\partial \sqrt{{\tilde{ψ}}_{j, q}^{(N)} {\tilde{ψ}}_{j, q}^{(N) *}}}{\partial O_{r}^{(n) ℛ}} \\ = \frac{1}{2 | {\tilde{ψ}}_{j, q}^{(N)} |} (\frac{\partial {\tilde{ψ}}_{j, q}^{(N)}}{\partial O_{r}^{(n) ℛ}} {\tilde{ψ}}_{j, q}^{(N) *} + \frac{\partial {\tilde{ψ}}_{j, q}^{(N) *}}{\partial O_{r}^{(n) ℛ}} {\tilde{ψ}}_{j, q}^{(N)}) \\ = \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {\frac{\partial {\tilde{ψ}}_{j, q}^{(N)}}{\partial O_{r}^{(n) ℛ}} {\tilde{ψ}}_{j, q}^{(N) *}} \\ = \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {\frac{\partial [ℱ {P_{j, r}^{(N)} O_{r}^{(N)}}]}{\partial O_{r}^{(n) ℛ}} {\tilde{ψ}}_{j, q}^{(N) *}} . \end{array}

Using the discrete Fourier transform, Eq. (25) becomes

\frac{\partial | {\tilde{ψ}}_{j, q}^{(N)} |}{\partial O_{r}^{(n) ℛ}} = \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \frac{\partial}{\partial O_{r}^{(n) ℛ}} [\sum_{m_{r^{'}}} \frac{1}{\sqrt{M}} 𝒫_{j, r^{'}}^{(N)} O_{r^{'}}^{(N)} exp (- i 2 π q \cdot r^{'})]},

where q · r = (q_x, q_y) · (x, y) = q_x x + q_y y = (m_{q_x} Δ_{q_x} m_xΔ_x) + (m_{q_y} Δ_{q_y} m_yΔ_y) and M = M_x M_y = [1/(Δ_{q_x} Δ_x)][1/(Δ_{q_y} Δ_y)]. For object slice n = N, Eq. (26) becomes

\begin{array}{l} \frac{\partial | {\tilde{ψ}}_{j, q}^{(N)} |}{\partial O_{r}^{(N) ℛ}} & = \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {{\tilde{ψ}}_{j, q}^{(N) *} [\sum_{m_{r^{'}}} \frac{1}{\sqrt{M}} 𝒫_{j, r^{'}}^{(N)} δ (r - r^{'}) exp (- i 2 π q \cdot r^{'})]}, \\ = \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {\frac{1}{\sqrt{M}} {\tilde{ψ}}_{j, q}^{(N) *} P_{j, r}^{(N)} exp (- i 2 π q \cdot r)} \end{array}

Insert Eq. (27) to Eq. (24) gives

\begin{array}{l} \frac{\partial ℒ}{\partial O_{r}^{(N) ℛ}} & = 2 \sum_{j} \sum_{m_{q}} (| {\tilde{ψ}}_{j, q}^{(N)} | - \sqrt{I_{j, q}}) \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {\frac{1}{\sqrt{M}} {\tilde{ψ}}_{j, q}^{(N) *} P_{j, r}^{(N)} exp (- i 2 π q \cdot r)} \\ = 2 \sum_{j} ℛ {χ_{j, r}^{(N) *} P_{j, r}^{(N)}} \end{array}

where

χ_{j, r}^{(N)}

is defined in Eq. (16). Analogous to Eqs. (27) and (28), the gradient for the imaginary part of the object slice is derived as

\frac{\partial | {\tilde{ψ}}_{j, q}^{(N)} |}{\partial O_{r}^{(N) ℐ}} = \frac{1}{| {\tilde{ψ}}_{j, q}^{(N)} |} ℛ {i \frac{1}{\sqrt{M}} {\tilde{ψ}}_{j, q}^{(N) *} P_{j, r}^{(N)} exp (- i 2 π q \cdot r)}

\frac{\partial ℒ}{\partial O_{r}^{(N) ℐ}} = - 2 \sum_{j} ℐ {χ_{j, r}^{(N) *} P_{j, r}^{(N)}},

thus giving the Wirtinger derivative for object slice N

\frac{\partial ℒ}{\partial O_{r}^{(N)}} = \sum_{j} χ_{j, r}^{(N) *} P_{j, r}^{(N)} .

The gradient of

P_{j, r}^{(N)}

can be derived similarly and expressed as

\frac{\partial ℒ}{\partial P_{j, r}^{(N)}} = χ_{j, r}^{(N) *} O_{j, r}^{(N)} .

The explicit dependence of object slice transmissivity

O_{r}^{(n)}

is given by the equation

P_{r}^{(n + 1)} = 𝒫_{Δ z_{n, n + 1}} {P_{r}^{(n)} O_{r}^{(n)}}

. Using the chain rule and assuming that different object slices are independent from each other and that the slice transmissivity is independent from the illumination incident on said slice, i.e.

O_{r}^{(n)}

is independent from

P_{r}^{(n)}

, we form the gradients of

O_{r}^{(n)}

in terms of the gradient of

P_{r}^{(n + 1)}

, namely

\begin{array}{l} \frac{\partial ℒ}{\partial O_{r}^{(n) ℛ}} & = \sum_{j} \sum_{m_{r^{'}}} (\frac{\partial ℒ}{\partial P_{j, r^{'}}^{(n + 1) ℛ}} \frac{\partial P_{j, r^{'}}^{(n + 1) ℛ}}{\partial O_{r}^{(n) ℛ}} + \frac{\partial ℒ}{\partial 𝒫_{j, r^{'}}^{(n + 1) ℐ}} \frac{\partial P_{j, r^{'}}^{(n + 1) ℐ}}{\partial O_{r}^{(n) ℛ}}) \\ = 2 \sum_{j} \sum_{m_{r^{'}}} (ℛ {\frac{\partial ℒ}{\partial P_{r^{'}}^{(n + 1)}}} ℐ {\frac{\partial P_{j, r^{'}}^{(n + 1)}}{\partial O_{r}^{(n) ℛ}}} - ℐ {\frac{\partial ℒ}{\partial P_{j, r^{'}}^{(n + 1)}}} ℐ {\frac{\partial P_{j, r^{'}}^{(n + 1)}}{\partial O_{r}^{(n) ℛ}}}) \\ = 2 \sum_{j} \sum_{m_{r^{'}}} ℛ {\frac{\partial ℒ}{\partial P_{r^{'}}^{(n + 1)}} \frac{\partial P_{j, r^{'}}^{(n + 1)}}{\partial O_{r}^{(n) ℛ}}} \\ = 2 \sum_{j} \sum_{m_{r^{'}}} ℛ {\frac{\partial ℒ}{\partial P_{r^{'}}^{(n + 1)}} \frac{\partial}{\partial O_{r}^{(n) ℛ}} 𝒫_{Δ z_{n, n + 1}} {P_{j, r^{'}}^{(n)} O_{r^{'}}^{(n)}}} . \end{array}

By approximating the angular spectrum calculation using discrete Fourier transforms, Eq. (33) becomes

\begin{array}{l} \frac{\partial ℒ}{\partial O_{r}^{(n) ℛ}} = & \frac{2}{M} \sum_{j} \sum_{m_{r^{'}}} ℛ {\frac{\partial ℒ}{\partial P_{r^{'}}^{(n + 1)}} \frac{\partial}{\partial O_{r}^{(n) ℛ}} \sum_{m_{q}} \sum_{m_{s}} P_{j, s}^{(n)} O_{s}^{(n)} exp (- i 2 π q \cdot s) \\ \times exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}}) exp (i 2 π q \cdot r^{'})} \\ \frac{2}{M} \sum_{j} \sum_{m_{r^{'}}} ℛ {\frac{\partial ℒ}{\partial P_{r^{'}}^{(n + 1)}} \sum_{m_{q}} \sum_{m_{s}} P_{j, s}^{(n)} δ (r - s) exp (- i 2 π q \cdot s) \\ \times exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}}) exp (i 2 π q \cdot r^{'})} \\ = \frac{2}{M} \sum_{j} \sum_{m_{r^{'}}} ℛ {\frac{\partial ℒ}{\partial P_{r^{'}}^{(n + 1)}} \sum_{m_{q}} P_{j, r}^{(n)} exp (- i 2 π q \cdot r) \\ \times exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}}) exp (i 2 π q \cdot r^{'})} \\ = \frac{2}{M} \sum_{j} ℛ {P_{j, r}^{(n)} \sum_{m_{q}} [\sum_{m_{r^{'}}} \frac{\partial ℒ}{\partial P_{j, r^{'}}^{(n + 1)}} exp (i 2 π q \cdot r^{'})] \\ \times exp (- i 2 π q \cdot r) exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}})} \\ = 2 \sum_{j} ℛ {P_{j, r}^{(n)} {(𝒫_{- Δ z_{n, n + 1}} {\frac{\partial ℒ}{\partial P_{j, r}^{(n + 1) *}}})}^{*}} \end{array}

Similarly for the imaginary part,

\frac{\partial ℒ}{\partial O_{r}^{(n) ℐ}} = - 2 \sum_{j} ℐ {P_{j, r}^{(n)} {(𝒫_{- Δ z_{n, n + 1}} {\frac{\partial ℒ}{\partial P_{j, r}^{(n + 1) *}}})}^{*}},

giving the Wirtinger derivative for object slice n < N,

\frac{\partial ℒ}{\partial O_{r}^{(n)}} = \sum_{j} P_{j, r}^{(n)} {(𝒫_{- Δ z_{n, n + 1}} {\frac{\partial ℒ}{\partial P_{j, r}^{(n + 1) *}}})}^{*} .

Similarly for

P_{j, r}^{(n)}

with n < N, the gradient is expressed as

\frac{\partial ℒ}{\partial P_{j, r}^{(n)}} = O_{r}^{(n)} {(𝒫_{- Δ z_{n, n + 1}} {\frac{\partial ℒ}{\partial P_{j, r}^{(n + 1) *}}})}^{*} .

Using Eqs. (36) and (37), we can recursively calculate the gradients of the object slices from slice n = (N − 1) to n = 1. At the first slice, i.e. n = 1, the gradient of the incident illumination is sum of Eq. (37) over j, as shown in Eq. (17).

The slice separation can also be optimized. For multiple slices in general, the gradient of Δz_n,n+1 is given by

\begin{array}{l} \frac{\partial ℒ}{\partial Δ z_{n, n + 1}} \\ = 2 \sum_{j} \sum_{m_{q}} (1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \frac{\partial {\tilde{ψ}}_{j, q}^{(N)}}{\partial Δ z_{n, n + 1}}} \\ = 2 \sum_{j} \sum_{m_{q}} (1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \frac{\partial}{\partial Δ z_{n, n + 1}} [\frac{1}{\sqrt{M}} \sum_{m_{r}} P_{j, r}^{(N)} O_{r}^{(N)} exp (- i 2 π q \cdot r)]} \\ = 2 \sum_{j} \sum_{m_{q}} (1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \\ \times \frac{\partial}{\partial Δ z_{n, n + 1}} [\frac{1}{\sqrt{M}} \sum_{m_{r}} 𝒫_{Δ z_{N - 1, N}} {\dots 𝒫_{Δ z_{n + 1, n + 2}} {𝒫_{Δ z_{n, n + 1}} {𝒫_{j, r}^{(n)} O_{r}^{(n)}} \\ \times O_{r}^{(n + 1)}} O_{r}^{(n + 2)} \dots} O_{r}^{(N)} exp (- i 2 π q \cdot r)]} \\ = 2 \sum_{j} \sum_{m_{q}} (1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \\ \times \frac{1}{\sqrt{M}} \sum_{m_{r}} 𝒫_{Δ z_{N - 1, N}} {\dots 𝒫_{Δ z_{n + 1, n + 2}} {\frac{1}{\sqrt{M}} \sum_{m_{q}} ℱ {P_{j, r}^{(n)} O_{r}^{(n)}} \\ \times i k \sqrt{1 - {| λ q |}^{2}} exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}}) exp (i 2 π q \cdot r) O_{r}^{(n + 1)}} \\ \times O_{r}^{(n + 2)} \dots} O_{r}^{(N)} exp (- i 2 π q \cdot r)} \\ = 2 \sum_{j} \sum_{m_{q}} (1 - \frac{\sqrt{I_{j, q}}}{| {\tilde{ψ}}_{j, q}^{(N)} |}) ℛ {{\tilde{ψ}}_{j, q}^{(N) *} \\ \times ℱ {𝒫_{Δ z_{N - 1, N}} {\dots 𝒫_{Δ z_{n + 1, n + 2}} {ℱ^{- 1} {ℱ {P_{j, r}^{(n)} O_{r}^{(n)}} \\ \times i k \sqrt{1 - {| λ q |}^{2}} exp (i k Δ z_{n, n + 1} \sqrt{1 - {| λ q |}^{2}})} O_{r}^{(n + 1)}} O_{r}^{(n + 1)} \dots} O_{r}^{(n + 1)}}} . \end{array}

Funding

Swiss National Science Foundation (SNSF) grant number 200021_152554

Acknowledgments

We acknowledge support from the Data Analysis Service (142-004) project of the swissuniversities SUC P-2 program. We also thank Marianne Liebi and Xavier Donath for their kind help on sample preparation and Oliver Bunk for his valuable input.

References and links

1. K. E. Herkenhoff, S. W. Squyres, R. Anderson, B. A. Archinal, R. E. Arvidson, J. M. Barrett, K. J. Becker, J. F. Bell III, C. Budney, N. A. Cabrol, M. G. Chapman, D. Cook, B. L. Ehlmann, J. Farmer, B. Franklin, L. R. Gaddis, D. M. Galuszka, P. A. Garcia, T. M. Hare, E. Howington-Kraus, J. R. Johnson, S. Johnson, K. Kinch, R. L. Kirk, E. M. Lee, C. Leff, M. Lemmon, M. B. Madsen, J. N. Maki, K. F. Mullins, B. L. Redding, L. Richter, M. R. Rosiek, M. H. Sims, L. A. Soderblom, N. Spanovich, R. Springer, R. M. Sucharski, T. Sucharski, R. Sullivan, J. M. Torson, and A. Yen, “Overview of the microscopic imager investigation during spirit’s first 450 sols in gusev crater,” J. Geophys. Res.-Planet 111, E02S04 (2006). [CrossRef]

2. R. Hovden, H. L. Xin, and D. A. Muller, “Extended depth of field for high-resolution scanning transmission electron microscopy,” Microsc. Microanal. 17, 75–80 (2011). [CrossRef]

3. R. Ng, “Digital light field photography,” Ph.D. thesis, Stanford University (2006).

4. M. Levoy, R. Ng, A. Adams, M. Footer, and M. Horowitz, “Light field microscopy,” ACM Trans. Graph. 25, 924–934 (2006). [CrossRef]

5. E. R. Dowski and W. T. Cathey, “Extended depth of field through wave-front coding,” Appl. Opt. 34, 1859–1866 (1995). [CrossRef] [PubMed]

6. R. H. Webb, “Confocal optical microscopy,” Rep. Prog. Phys. 59, 427 (1996). [CrossRef]

7. J.-A. Conchello and J. W. Lichtman, “Optical sectioning microscopy,” Nat. Methods 2, 920–931 (2005). [CrossRef] [PubMed]

8. M. Neil, R. Juškaitis, and T. Wilson, “Method of obtaining optical sectioning by using structured light in a conventional microscope,” Opt. Lett. 22, 1905–1907 (1997). [CrossRef]

9. H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: a novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). [CrossRef] [PubMed]

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

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

12. A. Diaz, P. Trtik, M. Guizar-Sicairos, A. Menzel, P. Thibault, and O. Bunk, “Quantitative x-ray phase nanotomography,” Phys. Rev. B 85, 020104 (2012). [CrossRef]

13. M. Holler, A. Diaz, M. Guizar-Sicairos, P. Karvinen, E. Färm, E. Härkönen, M. Ritala, A. Menzel, J. Raabe, and O. Bunk, “X-ray ptychographic computed tomography at 16 nm isotropic 3D resolution,” Sci. Rep. 4, 3857 (2014). [CrossRef] [PubMed]

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

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

16. T. M. Godden, R. Suman, M. J. Humphry, J. M. Rodenburg, and A. M. Maiden, “Ptychographic microscope for three-dimensional imaging,” Opt. Express 22, 12513–12523 (2014). [CrossRef] [PubMed]

17. A. Suzuki, S. Furutaku, K. Shimomura, K. Yamauchi, Y. Kohmura, T. Ishikawa, and Y. Takahashi, “High-resolution multislice x-ray ptychography of extended thick objects,” Phys. Rev. Lett. 112, 053903 (2014). [CrossRef] [PubMed]

18. K. Shimomura, A. Suzuki, M. Hirose, and Y. Takahashi, “Precession x-ray ptychography with multislice approach,” Phys. Rev. B 91, 214114 (2015). [CrossRef]

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

20. J. M. Rodenburg and R. H. T. Bates, “The theory of super-resolution electron microscopy via Wigner-distribution deconvolution,” Philos. Trans. R. Soc. London, Ser. A 339, 521–553 (1992). [CrossRef]

21. H. N. Chapman, A. Barty, S. Marchesini, A. Noy, S. P. Hau-Riege, C. Cui, M. R. Howells, R. Rosen, H. He, J. C. H. Spence, U. Weierstall, T. Beetz, C. Jacobsen, and D. Shapiro, “High-resolution ab initio three-dimensional x-ray diffraction microscopy,” J. Opt. Soc. Am. A 23, 1179–1200 (2006). [CrossRef]

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

23. J. M. Cowley and A. F. Moodie, “The scattering of electrons by atoms and crystals. I. A new theoretical approach,” Acta Crystallogr. 10, 609–619 (1957). [CrossRef]

24. P. Goodman and A. Moodie, “Numerical evaluations of n-beam wave functions in electron scattering by the multi-slice method,” Acta Crystallogr. Sect. A 30, 280–290 (1974). [CrossRef]

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

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

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

28. V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]

29. P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, “Probe retrieval in ptychographic coherent diffractive imaging,” Ultramicroscopy 109, 338–343 (2009). [CrossRef] [PubMed]

30. R. C. Gunning and H. Rossi, Analytic functions of several complex variables, vol. 368 (American Mathematical Soc., 1965).

31. L. Kaup and B. Kaup, Holomorphic functions of several variables: an introduction to the fundamental theory, vol. 3 (Walter de Gruyter, 1983). [CrossRef]

32. S. Gorelick, J. Vila-Comamala, V. A. Guzenko, R. Barrett, M. Salomé, and C. David, “High-efficiency Fresnel zone plates for hard x-rays by 100 kev e-beam lithography and electroplating,” J. Synchrotron Radiat. 18, 442–446 (2011). [CrossRef] [PubMed]

33. X. Huang, H. Yan, R. Harder, Y. Hwu, I. K. Robinson, and Y. S. Chu, “Optimization of overlap uniformness for ptychography,” Opt. Express 22, 12634–12644 (2014). [CrossRef] [PubMed]

34. B. Henrich, A. Bergamaschi, C. Broennimann, R. Dinapoli, E. F. Eikenberry, I. Johnson, M. Kobas, P. Kraft, A. Mozzanica, and B. Schmitt, “Pilatus: A single photon counting pixel detector for x-ray applications,” Nucl. Instrum. Methods Phys. Res. A 607, 247–249 (2009). [CrossRef]

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

36. J. Vila-Comamala, A. Diaz, M. Guizar-Sicairos, A. Mantion, C. M. Kewish, A. Menzel, O. Bunk, and C. David, “Characterization of high-resolution diffractive x-ray optics by ptychographic coherent diffractive imaging,” Opt. Express 19, 21333–21344 (2011). [CrossRef] [PubMed]

37. R. Hoppe, J. Reinhardt, G. Hofmann, J. Patommel, J.-D. Grunwaldt, C. D. Damsgaard, G. Wellenreuther, G. Falkenberg, and C. Schroer, “High-resolution chemical imaging of gold nanoparticles using hard x-ray ptychography,” Appl. Phys. Lett. 102, 203104 (2013). [CrossRef]

38. M. Guizar-Sicairos, A. Diaz, M. Holler, M. S. Lucas, A. Menzel, R. A. Wepf, and O. Bunk, “Phase tomography from x-ray coherent diffractive imaging projections,” Opt. Express 19, 21345–21357 (2011). [CrossRef] [PubMed]

39. M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. 33, 156–158 (2008). [CrossRef] [PubMed]

40. S. Mikula, I. Trotts, J. M. Stone, and E. G. Jones, “Internet-enabled high-resolution brain mapping and virtual microscopy,” Neuroimage 35, 9–15 (2007). [CrossRef] [PubMed]

41. “Brainmaps: An interactive multiresolution brain atlas,” http://brainmaps.org (2015).

42. M. Dierolf, P. Thibault, A. Menzel, C. M. Kewish, K. Jefimovs, I. Schlichting, K. Von Koenig, O. Bunk, and F. Pfeiffer, “Ptychographic coherent diffractive imaging of weakly scattering specimens,” New J. Phys. 12, 035017 (2010). [CrossRef]

43. B. R. Frieden, “Optical transfer of the three-dimensional object,” J. Opt. Soc. Am. 57, 56–66 (1967). [CrossRef]

44. J. W. Goodman, Introduction to Fourier optics (Roberts and Company Publishers, 2005).

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

46. D. J. Batey, D. Claus, and J. M. Rodenburg, “Information multiplexing in ptychography,” Ultramicroscopy 138, 13–21 (2014). [CrossRef] [PubMed]

47. T. B. Edo, D. J. Batey, A. M. Maiden, C. Rau, U. Wagner, Z. D. Pešić, T. A. Waigh, and J. M. Rodenburg, “Sampling in x-ray ptychography,” Phys. Rev. A 87, 053850 (2013). [CrossRef]

48. J. C. da Silva and A. Menzel, “Elementary signals in ptychography,” Opt. Express 23, 33812–33821 (2015). [CrossRef]

49. M. J. Humphry, B. Kraus, A. C. Hurst, A. M. Maiden, and J. M. Rodenburg, “Ptychographic electron microscopy using high-angle dark-field scattering for sub-nanometre resolution imaging,” Nat. Commun. 3, 730 (2012). [CrossRef] [PubMed]

50. L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an led array microscope,” Optica 2, 104–111 (2015). [CrossRef]

51. P. Li, D. J. Batey, T. B. Edo, and J. M. Rodenburg, “Separation of three-dimensional scattering effects in tilt-series fourier ptychography,” Ultramicroscopy 158, 1–7 (2015). [CrossRef] [PubMed]

X-ray ptychography with extended depth of field

Abstract

1. Introduction

2. Depth of field in ptychography

3. Reconstruction algorithms

3.1. Projection-based method

3.2. Maximum Likelihood

4. Experiment

5. Numerical studies

5.1. Image sectioning

5.2. Extended depth of field

6. Conclusions

Appendix - Derivation for multi-slice optimization-based ptychography

Funding

Acknowledgments

References and links

Cited By

Figures (5)

Equations (40)

Optics Express