Correction of axial position uncertainty and systematic detector errors in ptychographic diffraction imaging

Abstract. Ptychography is a diffraction imaging method that allows one to solve inverse problems in microscopy with the ability to retrieve information about and correct for systematic errors. Here, we propose techniques to correct for axial position uncertainty, detector point spread, and inhomogeneous detector response using ptychography’s inherent self-calibration capabilities. The proposed methods are tested with visible light and x-ray experimental data. We believe that the results are important for precise calibration of ptychographic experimental setups and rigorous quantification of partially coherent beams by means of ptychography.


Introduction
Ptychographic coherent diffraction imaging (PCDI) is a scanning microscopy technique that has gained wide popularity during the last decade. [1][2][3] Its ability to recover illumination and object information with phase contrast sensitivity without the need for sophisticated specimen preparation makes it suitable for both beam diagnostics 4,5 and quantitative imaging [6][7][8][9] in a variety of spectral ranges and experimental geometries. [10][11][12][13] However, PCDI relies on accurate experimental calibration, and incorrect modeling of the data can lead to reconstruction artifacts. Here, we report on correction schemes to mitigate detector-sided systematic errors, namely axial position uncertainty, point spread, and inhomegeneous response. We derive theoretical models for the aforementioned systematic errors and present correction strategies to improve reconstruction quality. One result is that the effect of axial position uncertainty is to scale the estimated coordinates of each scan position. An algorithm is presented to detect and correct for axial position errors based on lateral position correction. A second result is that both point spread and inhomegeneous response of the detector can be mitigated by mixed state ptychography. 14 The paper is organized as follows: Sec. 2 reviews details regarding PCDI and gives a theoretical description of the aforementioned systematic errors. Section 3 demonstrates experimental results on the classification and correction of axial position uncertainty and detector point spread. Concluding remarks are given in Sec. 4.

Ptychography
Under the thin object approximation, 2 the exit wave ψ downstream an object is modeled as the product of illumination P and object transmission O: where r ∈ R 2 denotes object coordinates, t j ∈ R 2 is the translation vector of the j'th scan position, and J is the total number of scan positions. The far-field diffraction intensities I are given as the modulus squared of the Fourier transform of the exit wave: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 4 1 1 where q ∈ R 2 denotes detector coordinates. In the case of partially coherent illumination, the model is generalized to incorporate coherence state mixtures, 14 involving a set of coherent modes: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 3 3 4 ψ m;j ðrÞ ¼ P m ðrÞOðr − t j Þ; m∈ f1; _; Mg; which are propagated independently through the optical system and added incoherently upon detection, i.e., E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 2 8 0 where the summation index corresponds to the m'th mode in the orthogonal mode expansion of the illumination. 15 Standard ptychography algorithms iterate between two constraints in reciprocal and real space. The reciprocal space constraint forces the iterate to comply with the measured intensities for each scan position, i.e., where n denotes iteration, ϵ > 0 is a small constant preventing division by zero, andψ n m;j ðqÞ is the Fourier transform of the m'th exit surface mode at the j'th object position. Scanning the object in overlapping regions provides ptychography with phase contrast sensitivity. 2 The overlap in scan positions together with the factorization assumption in Eq. (1) gives the real space constraint, which is typically imposed iteratively via nonlinear optimization algorithms. 2,[16][17][18] In this report, we use the real space updates: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 7 0 8 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 6 1 7 where β ¼ 0.5 is a feedback parameter, P nþ1 m is the updated m'th probe mode estimate, and O nþ1 is the updated object estimate. The latter updates are adapted versions of the extended ptychographic iterative engine (ePIE) for the case of partially coherent illumination. 14,19 Assuming quasimonochromatic radiation, the reconstructed coherent modes of the probe relate to the mutual intensity, Jðr 1 ; r 2 Þ, via 14 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 4 7 5 Jðr 1 ; The mutual intensity can be spectrally decomposed into 15 which is used here as a measure for the degree of spatial coherence of the beam. The reconstructed probe modes are constrained to be proportional to the orthonormal modes in the spectral decomposition of the mutual intensity. Combining and discretizing Eqs. (8) and (9) lead to the matrix equation: 21 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 1 9 1 PP Ã ¼ ΦΛΦ Ã ; where P ¼ ½P 1 ; _; P M ∈ C N×M contains the probe modes along its columns, Φ ¼ ½ϕ 1 ; _; ϕ M ∈ C N×M is orthonormal, Λ ∈ R M×M is diagonal containing the eigenvalues λ 1 ; _; λ M ≥ 0, and N is the total number of samples in the discretized beam. In practice, calculation of PP Ã ∈ C N×N is computationally expensive and should be avoided. This can be achieved using the truncated singular value decomposition for the orthogonalization of the probe: 22 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 7 5 2 where U ∈ C N×M , S ∈ R M×M , and V ∈ C M×M . Comparing E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 7 1 9 and Eq. (12) leads to the identities Φ ¼ U and Λ ¼ S 2 . A set of orthogonalized probe modes P ⊥ is derived from the set of nonorthogonalized probe modes P using the transformation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 6 5 5 P → P ⊥ ¼ US: The latter step allows one to identify the orthogonalized set of probe modes with the orthogonal modes and their corresponding weights in the spectral decomposition of the mutual intensity, namely P ⊥ ¼ ΦΛ 1∕2 . The orthogonalization of the probe modes [Eq. (15)] does not have to be calculated at every iteration of the ptychographic algorithm. In the software implementation used here, the orthogonalization step is carried out every 10th iteration to save computational resources. The next sections discuss departures from the ideal experimental model discussed above and correction strategies.

Axial Position Correction
The first imperfection considered is axial misalignment of the object-detector distance z, which has been reported to cause nonuniqueness and artifacts in the reconstruction. 23,24 The effect of axial displacement of the detector is to scale the real space coordinates t j attributed to each object frame Oðr − t j Þ. Assuming far-field diffraction, the real space pixel size Δx and field of view L are given by Δx ¼ λz∕D and L ¼ λz∕Δq, respectively, where D is the detector size and Δq is the detector pixel size. On an integer grid, a point at coordinate x is converted into dimensionless pixel coordinates N ¼ x∕Δx, where rounding is neglected. Then, the difference between the true pixel location, N t , and the measured pixel location, N m , is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 3 2 6 ; 3 3 3 From this, it follows that the coordinate error is proportional to the distance from the center of the coordinate system x and to the ratio Δz∕ðz þ ΔzÞ. The latter quantity can be approximated by Δz∕z assuming jΔzj ≪ z. Equation (16) is illustrated in Fig. 1. The sign of Δz in Eq. (16) leads to contraction or inflation of the scan grid: If the measured object-detector distance z m is larger than the ground truth object-detector distance z t , i.e., Δz ¼ z m − z t > 0, then ΔN < 0 and the estimated lateral positions are contracted with respect to the actual scan grid. In the other case, when z m < z t , i.e., Δz ¼ z m − z t < 0, then ΔN > 0 and the estimated scan grid is inflated. The former case is shown on the left and the latter case is shown on the right for a concentric 25 and a Fermat scan grid, 26 respectively. The proportionality of ΔN and x in Eq. (16) causes variable lateral position error that is most pronounced at the extremal points of the scan grid. Therefore, detector position uncertainty may be detected by the use of lateral position correction algorithms. [27][28][29] In this report, we use a random walk lateral position correction scheme similar to algorithms reported previously. 29,30 In our scheme, no underlying annealing schedule or drift model is assumed. At every iteration, the random walk position correction performs the following steps: (1) For the j'th position obtain unshifted and shifted exit wave modes, where the shift of the latter is given by AE1 pixel in both lateral directions. (2) Update exit wave modes according to Eq. (5) and calculate the error metric: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; 3 9 3 e j ¼ for the j'th position. Equation (3) shows that if e j;s < ce j;u , where e j;s and e j;u are the errors of the shifted and unshifted exit wave modes, respectively, then the j'th position is replaced by the shifted position. The constant c < 1 prevents premature position updates. In this report, we use c ¼ 0.95.

Partial Spatial Coherence as a Convolution Operation
Assuming partial spatial coherence (PSC), the far-field diffraction intensity is given by where Jðr 1 ; r 2 ; tÞ is the mutual intensity in the object plane. 31 For a Schell model field, the mutual intensity can be written as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 6 3 ; 1 4 1 Jðr 1 ; r 2 ; tÞ ¼ ψðr 1 ; tÞψ Ã ðr 2 ; tÞμðr 2 − r 1 Þ; where μ is the complex coherence factor. 32 Changing to centered coordinates r ¼ ðr 1 þ r 2 Þ∕2 and Δr ¼ r 2 − r 1 , the partially coherent diffraction intensity is given by [33][34][35] E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 3 2 6 ; 4 7 0 Using the autocorrelation and convolution theorems for Fourier transforms, this is equivalent to where ⊗ denotes convolution and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 3 2 6 ; 3 4 6 is the the coherent diffraction intensity. Equation (21) is Schell's theorem, which states that, given a complex coherence factor that is only a function of distance, the partially coherent diffraction intensity is equal to the coherent diffraction intensity convolved with the Fourier transform of the complex coherence factor. 35

Detector Point Spread as a Convolution Operation
If the detector is subject to point spread, the contrast in the diffraction pattern is reduced. In the ideal case, the detector discretely samples the diffraction intensity. In practice, the measured data on each pixel is integrated over a finite area. It is now shown that the integration over a finite pixel area can be modeled as a convolution operation, where only the one-dimensional case is derived. The extension to two dimensions is straightforward. The integrated intensity I int is given by where Δq is the pixel size of the detector. Fourier transformation of the latter expression with respect to q and changing orders of integration leads to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 6 3 ; 6 8 1Ĩ where f q is the Fourier conjugate variable to the detector coordinate q and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 6 3 ; 5 9 8PSF is the transfer function of the detector. Inverse Fourier transform Eq. (24) to get E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 6 3 ; 4 9 8 I int ðqÞ ¼ IðqÞ ⊗ PSFðqÞ; where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 6 3 ; 4 5 5 PSFðqÞ ¼ rect is the point spread function of the detector. The latter model for the point spread includes only the effect of area integration of the intensity incident of the detector. Other effects may additionally broaden the point spread function but depend on the electronic architecture of the detector at hand. Charge-coupled devices, for instance, may exhibit a broadened point spread due to tunneling between neighboring potential wells and charge diffusion during readout. 36,37 The main point of this section is that the detector point spread may be modeled as a convolution operation on a diffraction intensity, as described by Eq. 26.

Partial Spatial Coherence versus Detector Point Spread Ambiguity
If both partial coherence and detector point spread are nonnegligible, the observed diffraction intensity is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 6 3 ; 2 3 4 Iðq; tÞ ¼ I c ðq; tÞ ⊗μðqÞ ⊗ PSFðqÞ: (28) Under the above discussed approximations (Schell model beam and space-invariant detector point spread), ptychography's capability to reconstruct individual terms in the orthogonal mode decomposition has principally no means to distinguish the effects of partial coherence and detector point spread. However, it can be tested whether the terms in the orthogonal mode decomposition are due to partial coherence or detector point spread by changing the coherence defining aperture of the optical system. This is demonstrated in Sec. 3.2.

Mixed-State Ptychography in the Presence of Inhomogeneous Detector Response
In Ss. 2.2 and 2.4, it was discussed that the effect of PSC (for the case of a Schell model fields) and a space invariant detector point spread are mathematically modeled through a convolution operation with an a priori unknown kernel. Instead of estimating the kernel, the mixed state formalism allows one to incorporate convolution effects such as partial coherence, 14,38 sample vibrations, 39,40 stage movement during exposure 41,42 (fly scan effects), and point spread of the detector 14,40 into a mixed state probe. Mathematically, this is justified by Mercer's theorem, which states that a non-negative definite, hermitian, and square integrable kernel may be decomposed into a series of orthogonal modes. 43 A general cross-spectral density satisfies the latter conditions, 15 including the special case of Schell model fields. In addition, it was observed that the mixed state algorithm can mitigate static detector imperfections such as imhomogeneous response. 38 It is easily seen that an inhomogeneous detector response obeys the conditions for Mercer's theorem to apply. Let a detector responce be given by a real-valued function: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 3 2 6 ; 5 1 7 The hermiticity Tðq; q 0 Þ ¼ T Ã ðq 0 ; qÞ follows from the fact that T is real valued and diagonal. Further T is square integrable since it is bounded by 1 and nonzero only over a finite domain. Finally, it is non-negative definite since E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 3 2 6 ; 4 3 0 for arbitrary functions fðqÞ. This suggests that Mercer's expansion may be used to correct for static detector imperfections.

Detection and Correction of Axial Detector Position
To test the effect of axial detector position uncertainty, a visible light ptychographic scan was acquired. In the experiment, a He-Ne laser beam (λ ¼ 682.8 nm) was focused onto a 12 bit CMOS detector (IDS UI-3370CP-M-GL, 2048 × 2048 pixel, 5.5 μm pixel size) by a lens with a focal length of 100 mm, as depicted in Fig. 2(a). The object was placed a distance of 41 mm downstream the lens and z t ¼ 59 mm upstream the detector. Under the paraxial approximation in this configuration, the detector measures the scaled Fourier transform of the exit wave behind the object emulating far-field diffraction. 44 The probe was approximately Gaussian with a standard deviation of σ ¼ 110 μm as calculated after reconstruction by 45 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 ; 3 2 6 ; 1 2 8 corresponding to a full-width at half-maximum (FWHM) of 259 μm. Reconstructions were carried out using an ePIE algorithm 19 with the random walk position correction scheme described in Sec. 2.2 for distances z m ¼ 50; 59; 62 mm, as shown in Fig. 3. It is seen that for z m < z t (panel a) and for z m > z t (panel c), the initial scan grids are inflated and contracted, respectively, as compared to the corrected scan grids. The scaling of the corrected scan grids with respect to the initial scan grid is in agreement with Eq. (16) and may be used to correct for axial position uncertainty. This step may be automated or carried out manually as done here.

Detector Point Spread and Static Detector Imperfections
To test whether the coherent mode structure of the illuminating beam can be attributed to PSC or detector point spread, experiments were carried out at the MAXYMUS end station at the UE46-PGM2 beam line at the BESSY II synchrotron radiation facility. 46 A kinoform spiral zone plate of 32 μm diameter and 400 nm outer zone width was placed 3 m downstream crossed exit slits to generate a charge one vortex beam with a spot size of 1.5 μm to critically illuminate a binary test target at a photon energy of 800 eV. The experimental setup is depicted in Fig. 2(b). The region of interest on the test target was ∼6 μm in each lateral dimension with a scan step size of 150 nm. The high linear overlap 47 of 90% ensured stable recovery of the higher coherent mode structure (m > 1) of the beam. A total of 1600 diffraction patterns were recorded on a CCD (cropped to 128 × 128 pixel, 48-μm pixel size) placed 15 cm downstream the object resulting in a real space pixel size and field of view per object patch of Δx ¼ 38 nm and L ¼ 4.8 μm, respectively. The ptychographic reconstructions are shown in Fig. 4 for exit slit sizes of 10 μm × 10 μm (top) and 20 μm × 20 μm (bottom). For both exit slit openings, 1000 iterations were carried out with one (a, d) and nine (b, c, e, f) probe modes, four of which are shown (c, f). For the single mode reconstruction, the objects show artifacts at the outer region indicating that the tails of the probe were not reconstructed properly. For the multimode reconstruction, both objects show stronger similarity than for the single-mode reconstruction. Closer inspection of the multimode probe structure reveals that the degree of coherence did not significantly change between the two scans. The beam purities for the 10 μm × 10 μm and 20 μm × 20 μm slits are 63.9% and 60.5%, respectively. If partial coherence had been the cause for the coherent mode structure of the probe, then increasing the slit size by a factor of 2 would have resulted approximately in a twofold decrease in beam purity. From this, we conclude that the probe mode structure can mainly be attributed to detector point spread and static detector imperfections.

Conclusion
We have discussed two common detector-sided errors relevant in ptychographic diffraction imaging. The first error, axial misalignment, causes scaling of the correct lateral scan positions. It was shown that lateral position correction algorithms can be used to detect and correct for axial position uncertainty. The second detector-sided error discussed was detector point spread. We showed that this error can be identified by changing the spatial coherence defining element in the optical system, here the exit slit of a synchrotron beam line. If a change in exit slit size causes no response in the coherent mode structure of the illumination beam, the decreased diffraction pattern contrast can be attributed to detector point spread and static detector imperfections rather than PSC. However, while detector point spread is not attributable to decoherence in the beam, spectral mode decomposition, typically used for the representation of partially coherent beams, can be used to increase ptychographic reconstruction quality in the presence of detector point spread. Not shown in this work are initial tests on deconvolution strategies that we found to be less robust as compared to the mixed state algorithm. We believe that this is due to the dependence of the particular deconvolution ansatz on the specific underlying model. By contrast, the mixed state algorithm is conveniently applied since no a priori model of the detector error is required. Work on deconvolution of ptychographic data with explicit recovery of the detector point spread may be found elsewhere. 48 We believe that the results presented are important for improving reconstructions as well as rigorous quantification of partially coherent beams by means of ptychography.

Disclosures
The authors declare no conflicts of interest.