1. Introduction

Ptychography [1,2] has grown into a mature technique for x-ray, extreme ultraviolet (EUV), and electron microscopy. It has revolutionized synchrotron-based x-ray microscopy, where it improves upon previously existing scanning transmission x-ray microscopy (STXM) data analysis techniques [3–6]. Three major benefits of ptychography over STXM are: (1) decoupling of the illumination spot size from the achievable lateral resolution, (2) quantitative amplitude and phase contrast, and (3) access to wavefront diagnostics [7–11]. Similar benefits have subsequently been demonstrated for scanning transmission electron microscopes (STEMs) [12–14], where it recently produced micrographs at record-breaking resolution [15,16]. A parallel line of development is EUV laboratory-scale microscopy, where ptychography is a promising candidate for actinic inline metrology for lithography applications [17–19] and a tool for chemically-resolved microscopy [20–22]. In ptychography, a specimen is laterally scanned through a localized illumination, referred to as probe. A detector downstream of the specimen records a sequence of diffraction patterns. These observations lack phase information preventing direct inversion. Ptychography solves this problem by recording data from laterally overlapping specimen regions of interest. This acquisition scheme opens up the possibility for phase retrieval and simultaneous deconvolution of illumination and specimen information. Beyond operation with x-ray and electron radiation, ptychography has been implemented with extreme ultraviolet, visible, near-infrared, and terahertz radiation [19,23–26].

Fourier ptychography [27] follows a similar operational principle as (conventional) ptychography, denoted FP and CP, respectively, throughout this paper. In FP, a specimen is illuminated from different directions, typically steered by means of an LED array, which serves as a controllable condenser. A sequence of low-resolution bright and dark field images is recorded in a lens-based microscope. Changing the illumination direction amounts to shifting the object spectrum with respect to the pupil of the optical system. If the illumination direction is changed in such a way that two recorded images share information in the Fourier domain, phase retrieval techniques can be applied to separately reconstruct the object spectrum and the pupil of the optical system. Thus FP has three attractive features: (1) The low-resolution data can be stitched together to a large synthetic numerical aperture (NA), resulting in both a large field of view and high resolution. In contrast to most wide-field systems, FP thus does not trade-off resolution and field of view. (2) after conversion to real-space, the recovered object spectrum gives quantitative amplitude and phase maps of the sample; (3) the reconstructed pupil function enables aberration diagnostics of the optical system at hand [28–30]. While FP has mostly found applications in the visible domain, recent implementations using infrared radiation and x-rays have been reported [31,32].

1.1 Contribution

In both CP and FP, the recorded data jointly sample real and reciprocal space [33–37]. In CP, the probe serves as a real space window that selects local spatial frequency content. In FP, the pupil selects a low-resolution real space image from a localized Fourier space bandpass. In fact, the forward models of CP and FP are mathematically equivalent and the measured data cubes may be regarded as rotations of one another in phase space [36]. Although this equivalence is well-known, CP and FP have evolved into two separate communities with different algorithmic approaches and self-calibration techniques. Here, we report on a numerical data analysis toolbox, named PtyLab (code available online [38]), which places the equivalence of CP and FP at the center of its logical structure, resulting in three main contributions of this work:

In addition, PtyLab features a variety of algorithmic extensions as compared to currently available ptychography software packages. Some of these were previously reported by us and are now provided open source. This includes axial (zPIE) [39] as well as angle (aPIE) [40] correction engines, code to analyze ptychographic optical coherence tomography (POCT) data, efficient wave propagation algorithms both valid for high NA and polychromatic radiation, and detector subsampling (sPIE) [25,39–42]. Other novelties are reported here for the first time such as external reference wave ptychography. In addition, previously reported algorithms developed by other groups are included such as the extended ptychographic iterative engine (ePIE) [43], multislice (e3PIE) [44], mixed states [45], information multiplexing (PIM) [46], momentum acceleration (mPIE) [47], Tikhonov and total variation regularization [20,48,49], correlation-based lateral position correction [50], and orthogonal probe relaxation (OPR) [51]. In writing this manuscript, we pursued the goal of providing a concise overview of the various engines available to date in ptychography.

1.2 Related work

Most CP packages reported to date have focused on high performance computing, which is key for day-to-day user operation at large scale facilities, where beamtime is scarce and experimental feedback is needed quickly [52–59]. Another line of research has investigated the capabilities opened up by modern automatic differentiation (AD) and machine learning (ML) toolboxes [60,61]. AD frameworks offer flexibility as they simply require the user to specify a desired forward model, typically specified in the form of a suitable loss function and possible regularizers. This approach is convenient for the user as it dispenses with the need to derive challenging gradient expressions, the latter of which oftentimes involve complex-valued (Wirtinger) derivatives [62,63]. It is thus for instance straightforward to switch from one regularizer to another without analytically deriving and programming the underlying gradient expressions into the underlying software. ML approaches, in particular those based on neural networks, have been used to significantly speed up the reconstruction process, lower the sampling requirements in the raw data, and to embed denoising priors [64,65]. However, neural network approaches need to be trained based on data sets that have already been solved for the corresponding real space images. In Ref. [64] the neural network was trained based on the solution of an iterative solver. Moreover, training a neural network capable of solving ptychography problems is a memory-consumptive, large-scale computational challenge, that cannot be performed on small hardware architecture. We thus believe the need for memory-efficient but possibly slower iterative algorithms remains, despite the exciting possibilities opened up by neural networks [66]. From the above referenced work, we shortly describe some of the features of two prominent code projects, namely PtychoShelves [57] and PtyPy [54], to illustrate some of the different design choices made in PtyLab.

PtychoShelves [57] is a Matlab-based software package for ptychography, designed with large-scale synchrotron facilities in mind. Shelves refer to the modular coding framework representing bookshelves, from which desired books (e.g., detector module, reconstruction engine etc.) can be taken out and inserted into the processing pipeline. To provide data handling across synchrotron facilities worldwide, PtychoShelves supports commonly used X-ray detectors as well as instrument control software. Reconstructions are highly-optimized for Matlab-based GPU acceleration as well as CPU processing through reconstruction engines written in binary C++ code and parallelized through the OpenMP multiprocessing interface. The C++ code supports Difference Map [8] and Maximum Likelihood [48] engines, together with other features such as mixed state [45] and multi-slice ptychography [44]. A wider range of reconstruction features are available through Matlab-based GPU engines, including an iterative least-squares solver for maximum-likelihood ptychography [55], orthogonal-probe-relaxation [51], near-field ptychography [67], and position correction [68].

PtyPy [54] is an open-source ptychography framework written in Python. It follows the Python coding style and is therefore modularized and object-oriented. The physical models are abstracted, which results in readable and concise code interfaces at a user level. The key element is the so-called POD class (probe, object, diffraction), which holds the access rule for a single position, object mode, and probe mode. For parallelization of the reconstructions, PtyPy uses a message passing interface (MPI), which allows for flexible usage on standard office PCs, workstations, and clusters. MPI favors non-sequential reconstruction algorithms that can be parallelized (e.g. Difference Map [8] and Maximum Likelihood [48]). So far, a broad range of forward models are implemented (e.g. mixed state ptychography [45], near-field ptychography [67], and lateral position correction [68]). The PtyPy framework is actively developed and novel features (e.g. GPU acceleration) are constantly added.

Both PtychoShelves and PtyPy are powerful ptychography data analysis platforms, but their design for high-performance computing poses an entry barrier for simple, one-off reconstructions in an academic-lab setting. In such cases, rapid code prototyping and ease-of-use can be more desirable than highly-optimized data handling and reconstructions.

Unlike CP, FP has not seen the same wide-spread use within research institutions that resulted in well-developed and maintained coding platforms. The existing open-source FP codes in [69–74] exist mainly to supplement publications by providing only minimal working examples with limited functionality intended. Recently an attempt has been made to provide a Matlab-based FP reconstruction platform [75], which among other features provides raw data denoising, GPU processing, and LED misalignment correction.

Our goal here is to bridge the gap between CP and FP, thereby allowing for a cross-pollination of the two domains and a unified algorithmic framework. As compared to the above highlighted software packages PtyLab is less focused on high performance and distributed computing, but puts emphasis on providing an ptychography ecosystem for researchers interested in rapid prototyping and exchanging algorithms - across modalities and programming languages.

1.3 Outline

In Section 2 we revisit the idea of reciprocity, which formalizes the equivalence between CP and FP - the central idea for the unified software design in PtyLab. Section 3 details language-specific implementation details in Matlab, Python, and Julia. Section 4 serves as a comprehensive overview of the available forward models and the corresponding optimization algorithms. Practical features for scan grid optimization are described in section 6. PtyLab is released with various data sets and hands-on tutorials, which are described in section 7.

2. Implications of reciprocity for ptychography

One may think of the data sets recorded in ptychography in analogy to a musical score, where frequency information is prescribed at particular signatures in time. Once such a time-frequency, or phase-space, representation is given in form of a musical score, we can convert this information into either the time or the frequency domain. For example, we can digitally record a concert and Fourier transform the resulting signal. These processing steps would involve the temporal waveform and its frequency spectrum, respectively. Likewise, ptychography jointly samples real and reciprocal space representations of a signal, where for simplicity we ignore the additional complication of phase retrieval for the moment. The goal of ptychography is to convert partial phase-space information of a signal into a pure space or a pure spatial frequency representation. Physically, the phase-space description of ptychography [33] is intimately connected to the principle of reciprocity [77], which states that by interchanging the illumination and detection direction in an optical system identical data sets can be observed.

We would like to distinguish two types of reciprocity in ptychography. Type-I reciprocity refers to the ability to design separate CP and FP optical systems, both of which produce 4D data cubes which are essentially related by a phase space rotation [36]. In addition, we define type-II reciprocity, which refers to the ability to algorithmically convert a 4D data cube from one domain to the other. Thus type-I reciprocity is essentially a statement about the ability to design different CP and FP hardware embodiments producing the same 4D experimental data cube. Type-II reciprocity is a matter of data processing: once a 4D data cube is measured in either a CP or a FP representation, it can be converted into the other respective domain and subsequently reconstructed.

Figure 1 illustrates this idea of reciprocity in the context of ptychography. In CP (Fig. 1(a)), an object is laterally translated against a typically focused beam. This probe localizes the origin of the measured signal in real-space (scan coordinates $\boldsymbol {s}$). A sequence of diffraction patterns is captured on a pixelated detector, which is assumed here to be located in the far field (spatial frequency $\boldsymbol {u}$). Hence the data cube in CP consists of a sequence of angular scattering maps, each corresponding to a particular real-space specimen region localized to the spatial extent of the incident probe (Fig. 1(b)). In FP, an object is held at a fixed location under variation of the illumination angle (Fig. 1(d)). Thus the angular spectrum emanating from the sample is shifted over the finite-sized lens aperture. The pupil serves to bandpass filter the object’s spatial frequency spectrum, resulting in a data cube consisting of dark and bright field images (Fig. 1(c)). Thus the data cube in FP consists of a sequence of real-space images, each corresponding to a particular reciprocal space portion of the object spectrum localized to the passband admitted by the pupil of the imaging system. In summary, both flavors of ptychography sample real and reciprocal space.

 

Fig. 1. Illustration of the operation principle and equivalence of conventional and Fourier ptychography. (a) An object is laterally translated against a localized illumination profile. (b) In CP, the recorded data cube is a sequence of diffraction patterns, providing spatial frequency ($\boldsymbol {u}$) information for each scan position ($\boldsymbol {s}$). Each detector pixel alone contains a sequence real-space information that may be reshaped into a low-resolution real-space image of the sample. (c) In FP, the recorded data cube is a sequence of low-resolution bright and dark field image plane ($\boldsymbol {s}$) data corresponding to bandpass-filtered versions of the object spectrum. The shifts with respect to the pupil are controlled by shifting the illumination direction ($\boldsymbol {u}$). (d) Single-lens FP experimental configuration. Data in panel (b) and (c) from [76]

Download Full Size | PDF

We illustrate the aforementioned types of reciprocity in two ways: First, consider replacing the detector in the FP setup (Fig. 1(d)) with a pixelated light source. Turning on a single point source at a time and placing the detector in the far field of the sample, we can record a CP data set by scanning the point source location, as pointed out in a recent review article [2]. Of course, this hardware conversion ability faces practical limits imposed by the lens, which may cause a space-variant probe, but this complication is ignored here. Thus via hardware modification we can convert a CP experimental setup into and a FP system, which is a type-I reciprocity. Type-II reciprocity concerns the captured data itself and does not require any hardware modifications. It is possible to convert the measured data cube from one modality to the other. Suppose in Fig. 1(a) we scan the sample (for conceptual simplicity) on a regular raster grid and record a sequence of diffraction patterns. Then each pixel of the detector can be regarded as a traditional (single-pixel) scanning microscope data set. The data on each individual pixel may directly be reshaped into a two-dimensional real-space image, simply by permuting its dimension in correspondence with the scan trajectory. Practically, aperiodic translation trajectories and high NA effects require interpolation techniques. However, at low NA and using raster scan grids the data reshaping operation can be implemented with a single line of code (namely, a permute operation), converting for instance a CP data set into a sequence of low-resolution bright and dark field images, the latter of which constitutes the raw data for FP. While we described type-II reciprocity phenomenologically in this section, a mathematical proof of this conversion ability is provided in the appendix. The mathematical details also elucidate the correspondence between reconstructed quantities in CP and FP. We provide online tutorials that illustrate the conversion between CP and FP [38].

The ability to convert CP and FP data has a bearing on the computational complexity of inversion algorithms underlying ptychography. Suppose we are given a CP data cube consisting of diffraction patterns with $U^2$ pixels at $S^2$ scan points. A single iteration of a parallelized ptychography solver (for example difference map [8]) requires us to numerically propagate exit waves from all scan positions to the detector plane and back. The Fourier transform operations involved will have a computational complexity of $\mathcal {O}\left [S^{2}\cdot U^{2}\cdot \log \left (U\right )\right ]$ if we work in the CP domain, while it scales with $\mathcal {O}\left [U^{2}\cdot S^{2}\cdot \log \left (S\right )\right ]$ if we convert the data into the FP domain. The difference in the log-terms can result in a practical speed up, provided that the number of detector pixels per dimension $U$ and the number of scan positions per dimension $S$ is largely different.

In summary, utilizing type-II reciprocity is the central motivation for the design of PtyLab: CP and FP data can be converted into each other. A unified data analysis framework thus allows to migrate between the two modalities. A benefit of this data conversion ability is the applicability of diverse inversion algorithms and self-calibration routines in the domain in which they are most conveniently applied. Another benefit of reciprocity is the trade-off in computational complexity.

3. Code structure

In this section we describe structural workflow in PtyLab. Our overall goal is to provide a code that enables flow of algorithmic ideas and rapid prototyping beyond the boundaries of modality (CP/FP) and programming language (Matlab/Python/Julia). Thus collaborators with different programming language preferences or from different communities (e.g. synchrotron-based CP versus visible light FP) can easily exchange code, without being perfectly literate in the other programming language. This approach comes at the benefit of a unified structure and naming convention but at times at the expense of certain language-specific programming conventions. The following subsection describes the common structure independent of programming language. Subsequently, we address differences in the Matlab, Python, and Julia implementations.

3.1 Platform- and modality-independent workflow

A high-level overview of the PtyLab’s workflow is illustrated in Fig. 2. Assuming CP or FP data is stored in a local folder on the user’s computer, the first step is preprocessing the data (see Fig. 2). Preprocessing converts the raw data into a PtyLab class object. The user specifies physical input (illumination wavelength), hardware properties (detector pixel pitch, binning), and geometric parameters (sample-detector distance [CP], lens magnification [FP], sample scan trajectory [CP], illumination angles [FP]). In addition, the user specifies a forward model that describes the propagation from end to end (CP: probe to detector, FP: pupil to detector). The preprocessing pipeline then writes a single or multiple PtyLab class objects into a hdf5 file (see Fig. 3) [78].

 

Fig. 2. Workflow in PtyLab. Experimental data is converted into a preprocessed hdf5 data set. The remaining parameters controlling algorithmic and monitoring behavior and required for reconstruction are set by an initialization routine. Various reconstruction engines and forward models can be chosen to analyze the preprocessed data. After the reconstruction is finished the reconstructed data is written into a final hdf5 file.

Download Full Size | PDF

 

Fig. 3. PtyLab HDF file structure of preprocessed (a) and reconstructed (b) data. The blue boxes (top) show the mandatory fields, with specific differences between CP (yellow) and FP (green) data. The grey box in panel (a) shows optional fields, that are nevertheless recommended. Both CP and FP reconstructions struggle to converge when background is not appropriately subtracted or accounted for in the forward model, subtraction being the easier route. Initial estimates for the probe diameter in CP and the pupil diameter in FP are recommended to be specified (both referred to as entrancePupilDiamater). This can aid initial convergence in CP. Moreover, the circle fitting routine for position calibration in FP [79], which is used in PtyLab, requires an estimate of the pupil diameter. Arrays indicated with (*) are specified in SI units.

Download Full Size | PDF

Second, the reconstruction script loads the preprocessed data. An initialization function generates uniform or radomized starting estimates for the probe (CP) or pupil (FP) and object (CP) or object spectrum (FP). In each the probe/pupil, object/object spectrum, and detector planes, meshgrids are calculated. These meshgrids depend on the specified physical, hardware, and geometrical parameters, as well as a forward model (propagator) that describes the mapping between probe (CP) or pupil (FP) and detector planes. A variety of propagation models can be specified, including angular spectrum (AS), scaled angular spectrum (SAS), Fresnel (Fresnel), Fraunhofer (Fraunhofer), and tilted Fresnel diffraction. The latter is relevant for non-coplanar reflection geometries and is typically performed only once on the raw diffraction data stack, provided no angle correction is performed (see subsection 5.7 for further details). Fraunhofer and Fresnel diffraction distinguish each other by an incorporation of a quadratic phase inside the propagation model. This quadratic phase may be absorbed into the probe/pupil function and can be compensated for post-reconstruction when a quantitative analysis of the reconstructed wavefront (CP) or pupil (FP) is of interest to the user.

3.2 Matlab structure

The Matlab code structure is shown in Fig. 4(a). Here an object of class PtyLab is generated. Its first-order properties (obj."firstOrder") contain the physics as well as the geometry of the experiment. In addition, there are second-order properties (obj."firstOrder"obj.params."secondOrder"), which are algorithmic parameters (obj."firstOrder"obj.params), monitoring control (obj.monitor), propagators as part of the forward model (obj.propagator), and file export parameters (obj.export). Certain notational conventions are noteworthy: the diffraction or image data is contained in obj.ptychogram, a term borrowed from time-frequency analysis where the raw data is often referred to as spectrogram [80]. The dimensions of (obj.ptychogram) are $({\tt y, x, numFrames})$, which is different from the Python convention (see subsection 3.3). The order along the first two dimensions stems from Matlab’s most convenient use when adhering to row-column convention. Similarly, obj.positions follows row-column convention.

 

Fig. 4. The Matlab code structure (left) comprises a single class, which contains all field relevant for ptychographic data analysis. The Matlab class is organized into first- and second-order properties. First-order properties contain physical information (e.g. wavelength) and geometrical parameters of the experiment. Second-order properties are mainly found in params, which contains algorithmic properties (step sizes, number of iterations, etc.) that are optimized during data analysis. Other second-order properties comprise monitoring behaviour (monitor), specification of the wave propagation model (propagator), and input-output control (export). The Python code (right) consists of five separate classes: ExperimentalData, Reconstruction, Params, Monitor, and Engines. The Julia implementation consists of 4 main abstract types called ExperimentalData, Reconstruction, Params, and Engines.

Download Full Size | PDF

3.3 Python structure

The Python structure is similar in idea to the Matlab structure, but is designed with a stronger emphasis on modularity. As shown in Fig. 4(b), the Python implementation contains five classes: ExperimentalData, Reconstruction, Monitor, Params, and Engines. Most but not all of these classes reflect second-order properties in the Matlab structure. The ExperimentalData class imports data from a preprocessed .hdf5 file, checks if all required parameters for a ptychographic reconstruction are included, and saves them into an instance that is immutable.

The Reconstruction class takes the ExperimentalData instance as the input, and creates a mutable instance containing attributes that are optimized during a reconstruction process, e.g. the probe/pupil, and the object, as well as attributes that are related to the optimizable parameters, e.g. the error, the coordinates and meshgrids. Note that in the Python implementation, the probe/pupil and the object are set as 6D arrays with the fixed axes [nlambda, nosm, npsm, nslice, row, col], which are the number of wavelength, object state mixtures, probe state mixtures, slices (for multislice ptychography), and rows as well as columns.

The Monitor class is used to display a reconstruction process. Equivalent to its Matlab counterpart, one or two figures are created depending on the verbosity level set by users. A default figure shows the updated object, probe, and reconstruction error. An optional figure shows the comparison of the measured and estimated diffraction patterns. The update frequency of the plots can also be controlled by users (Monitor.figureUpdateFrequency).

The Params class holds parameters that determine how a reconstruction is performed, for instance whether a reconstruction is carried out on a CPU or a GPU, the propagator type such as Fraunhofer, Fresnel, (scaled) angular spectrum (ASP, scaledASP), etc., and whether the order of position iterations is random or sequential. Switches and parameters of various regularization types are also included in the Params instance, for example controlling how frequenctly orthogonalization is applied in the context of a mixed-states reconstruction.

The Engine class consists of a BaseEngine as a parent class, and other child engine classes, for instance ePIE [43], mPIE [47], zPIE [39], aPIE [40], and qNewton [81]. All four instances of ExperimentalData, Reconstruction, Params, and Monitor are taken as inputs for a chosen engine, then get modified/updated by the engine, and can be passed to a different engine easily. Each engine stores its own attributes such as the number of iterations (numIteration), and the update step sizes for the probe/pupil (betaProbe) and object (betaObject).

3.4 Julia structure

PtyLab.jl is the most recent translation of PtyLab to Julia. Due to the differences in Julia to Matlab and Python consequently small differences exist in the implementation but most of the common principles still hold. The amount of features is less than in the other two packages since its focus was on performance first. The basis are four main types ExperimentalData, Reconstruction, Params, and Engines. Engines is an abstract type which is subtyped (indicated by <:) by specific algorithms (such as ePIE <: Engines) as composite types. Via this mechanism, generic functions can be commonly used by all Engines solvers. However, Julia’s multiple dispatch allows that different functionalities can be specified if they belong to a ePIE or a zPIE solver. Params is a composite type storing second-order-properties. Further, ExperimentalDataCPM <: ExperimentalData exists and similarly ReconstructionCPM <: Reconstruction to store the experimental data and the reconstruction data. During the iteration of the algorithms, the fields of ReconstructionCPM are allowed to change. Julia’s language features allow for a functional style of programming implying that memory buffers are not explicitly exposed to the user but instead are implicitly stored via closures.

4. Inverse modeling

The inverse modeling workflow in PtyLab consists of several modular processing steps, which are shown in Fig. 5. All optimization algorithms in PtyLab iterate between the object (CP) / object spectrum (FP) plane (orange) and the detector plane (green), where the two planes are linked via a suitable propagation model (yellow). We subdivide this section into several parts, describing the individual modules that the user can stack up to build customized data analysis pipelines.

 

Fig. 5. Optimization workflow. The schematic illustrates the building blocks of a user-defined reconstruction routine in PtyLab. In the object plane, the forward exit wave model as well as the inverse model for the object and probe gradients, subject to several regularization options, are specified. In the detector plane, the underlying noise model and various regularization options lead to an optimal update for the estimated detector wave. After initialization of the object and probe, the reconstruction engine iterates between the object and detector plane until the reconstruction error is low or other stopping criteria are satisfied.

Download Full Size | PDF

4.1 Forward model

In CP and FP the goal is to retrieve a wide-field high-resolution reconstruction of a sample of interest. In addition, the probe (CP) and pupil (FP) of the imaging system are recovered. A forward model links the estimated detector intensity $I$ to the object $O$ and probe $P$ (CP) or object spectrum $\tilde {O}$ and pupil $\tilde {P}$ (FP),

Here $\mathcal {D}$ describes wave propagation between the sample (CP)/pupil (FP) and the detector plane, $\boldsymbol {r}$ refers to spatial coordinates, and $\boldsymbol {q}$ refers to reciprocal space coordinates. The index $j=1,\ldots,J$ denotes scan position. For simplicity, we drop the coordinate dependence and use the CP notation throughout; the conversion to FP of all results discussed below is straightforward. The symbol $O$ refers to a particular object box of equal size as the probe FOV (compare red region in Fig. 2). The entire object field of view is denoted by $O_{\textrm {FOV}}$ (compare blue region in Fig. 2).

In the presence of noise in the observed signal, for instance caused by photoelectric conversion and read-out, we cannot expect to find a combination of sample and prope/pupil that exactly matches the recorded data. In what follows we therefore assume the measured data $m$ to arise as a probabilistic response to the true intensity $I$ incident on the detector and discuss several maximum likelihood estimation (MLE) models [48,55,81–84]. These MLE models aim at estimating the most likely combination of object and probe, a viewpoint extended by the addition of maximum a posteriori (MAP) estimation, which originates from a Bayesian viewpoint and enables flexible embedding of regularization into the forward model [48,84,85]. Before going into the details of inverse models, we briefly review the continuous and discrete viewpoints on optimization, which are both encountered in the ptychography optimization literature [45,47,48,55,81–84,86–88].

4.2 Inverse modeling

In this section we review various forward models implemented in PtyLab. A summary of these forward models is given in Fig. 6. We first describe general techniques to tackle the inverse problem underlying ptychography. Subsequently, we detail the individual solvers that allow the user to build and invert modular forward models.

 

Fig. 6. Selection of forward models implemented in PtyLab: (a) The basic coherent diffraction model assumes the thin element approximation (TEA), where the exit wave $\psi _j$ is modeled as a product of probe $P$ and object box $O_j$ a scan position $j$. The exit wave is propagated into the detector plane via a suitable propagator $D$. (b) In mixed state ptychography the object interacts with $k$ mutually incoherent probe modes, giving rise to independent exit waves $\psi _{j,k}$. These exit waves are propagated into the detector plane and incoherently added to form the intensity forward model. (c) Multispectral ptychography. Here a polychromatic probe interacts with a dispersive object, both of which are functions of wavelength $\Lambda$. (d) In multislice ptychography the exit wave is modeled using the beam propagation method (BPM), which models a three-dimensional object to consist of several two-dimensional slices (index $s$). Inside each slice the TEA is used, while the propagation between slices is carried out via angular spectrum propagation $A$. (e) Orthogonal probe relaxation can model scan position dependent probes $P_j$ as a linear combination of mutually coherent orthogonal basis modes $U_k$. (f) A coherent external reference wave can be added to the forward model.

Download Full Size | PDF

4.2.1 Continuous viewpoint

In the continuous viewpoint, we aim to minimize a cost functional

If Eq. (5) is not amenable to a direct solution, we can iteratively solve it by seeking the steady state solution of the diffusion equation [91,92]

4.2.2 Discrete viewpoint

The discrete viewpoint is used when considering inverse problems over sampled functions. In this case, we oftentimes wish to minimize the sum of squares cost function

A non-iterative solution is formally obtained by setting the gradient in Eq. (9) to zero and solving for $f,$

4.3 Maximum likelihood (MLE) estimation

We now discuss models for the detector noise commonly used in ptychography. Two particularly prominent models that have been addressed [48,55,81–83] are the Poisson likelihood

Computing the Wirtinger derivatives of the negative log-likelihood $\mathcal {L}=-\sum _{\boldsymbol {q}}\log \left (p\left [m\left |I\right.\right ]\right )$ of Eqs. (12) and (13) results in [48,55,81,83]

It appears to us that the Anscombe forward model is used by the vast majority in the ptychography literature. Although the Poisson distribution works in practice, we have observed that the Anscombe model is more robust in practical data analysis. One has to keep in mind that the Poisson model assumes that the photoelectric counting distribution’s mean equals its variance and is only valid for shot-noise limited data - a somewhat restrictive assumption considering the manifold fluctuations that are present in typical experiments, including partial spatial and temporal coherence effects as well as detector read out. We note that other models for the statistics of photoelectic counting distributions have been proposed in the literature, albeit to our knowledge they not yet have been used for ptychography. Noteworthy is the negative binomial distribution

4.4 Maximum a posteriori (MAP) estimation

The MLE approach in the previous subsection can be extended by MAP estimation, which introduces prior knowledge into the reconstruction process. In MAP the detector intensity is regarded as a random variable with underlying probability density $p\left [I\right ]$. Since the detector intensity $I=\tilde {\psi }^{*}\tilde {\psi }$ is a function of the real space object and probe, namely $\psi \left (\boldsymbol {r}\right )=P\left (\boldsymbol {r}\right )\cdot O\left (\boldsymbol {r}-\boldsymbol {r_{j}}\right )$ and $\tilde {\psi }\left (\boldsymbol {q}\right )=\mathcal {F}_{\boldsymbol {r}\rightarrow \boldsymbol {q}}\left [\psi \left (\boldsymbol {r}\right )\right ]$, MAP opens up a convenient way to formulate and impose constraints in the inverse problem underlying ptychography. The optimization problem is then

4.4.1 Proximal detector updates

The detector update can be restricted to small step sizes by introducing the proximal prior

These updates are referred to as proximal. A large value of the tuning parameter $\lambda$ forces the updated wave $\tilde {\psi }_{n+1}$ to remain in the proximity of the previous estimate $\tilde {\psi }_{n}$. Intuitively, the updates along the gradient directions in Eqs. (14) and (15) enforce the magnitude of the updated wave to be equal to the measured data, either in intensity or modulus for the Poisson and Anscombe gradients, respectively. However, due to noise, sequential update schemes can only be as certain as the noise in a single diffraction pattern permits. Proximal gradients incorporate the memory of previous updates and do not naively accept the update suggested by the data. The gradient direction suggested by the data is followed in case that the deviation from the current estimate is small. It is conjectured here that this incorporates dose fractionation effects into ptychography, resulting in superior signal-to-noise in the reconstruction. This is supported by previous reports that observed improved reconstruction quality using proximal gradients for Gerchberg-Saxton type phase retrieval [104] and ptychography [84,87]. The cost functions above (Eqs. (22) and (23)) result in the update steps

We note that the proximal Poisson update is different from the result in [104], since here the prior is Gaussian in modulus, while in the related work the prior is Gaussian in intensity. The approach reported here avoids the need to solve a cubic polynomial to compute the proximal Poisson update.

4.4.2 Proximal probe and object updates via (e)PIE

While in the previous section a proximal update step has been discussed in the detector update step, we may impose a similar type of regularization for the probe and object update, which has been shown to result in the ptychographic iterative engine (ePIE) [45]. This derivation is reviewed here from the discrete viewpoint outlined above. Considering the cost function

The particular choice

Later the extended ptychographic iterative engine (ePIE) was proposed [43], which uses the regularization matrix [45]

An intuition of the difference between PIE and ePIE can be obtained in the limit of $\epsilon \rightarrow 0$ and $\beta =1$, for which we get

$\Gamma _{\textrm {PIE}}$ is small, when $\left |P\right |$ approaches $\max \left (\left |P\right |\right )$ or $0$. Thus PIE allows object updates for locations where the probe amplitude is small or large. In contrast, $\Gamma _{\textrm {ePIE}}$ is small, only when $\left |P\right |$ approaches $\max \left (\left |P\right |\right )$. Thus ePIE allows object updates only for locations where the probe amplitude is large. The derivation of the probe update is similar, resulting in a joint optimization of $P$ and $O$. The robustness of ePIE is attributed to the penalized step size at low probe intensities. However, in FP the large dynamic range of the object spectrum can cause problems in conjunction with ePIE. The pupil would be updated only in the center of k-space, where the object spectrum exhibits values close to its maximum amplitude. High illumination angles produce dark field images, which have a reduced signal-to-noise ratio as compared to bright field images from lower illumination angles. In CP the illuminating beam is always aligned with the detector resulting in images with similar signal-to-noise ratio as compared between scan positions. For this reason, we use the PIE-type update by default in PtyLab for FP and the ePIE-type update for CP data analysis. Other choices of regularization can be embedded into the reconstruction routine by changing the weight matrix $\Gamma$. We note in passing that the PIE-type update rule has been identified as a quasi-Newton algorithm [81].

4.4.3 Tikhonov regularization

A popular idea in solving inverse problems is Tikhonov regularization. The general idea is to add an additional term to the cost function, which penalizes variations in the object,

We emphasize that the regularization term in Eq. (26) in the previous subsection penalizes fluctuations at for a fixed object pixel between successive iterations. In this subsection the regularization term in Eq. (34) penalizes fluctuations between neighboring object pixels. Applying functional gradient descent, as described by Eq. (7), to the cost in Eq. (34) gives

4.4.4 Total variation regularization

Total variation (TV) regularization penalizes changes in the image while to a certain degree preserving edge features [91]. The corresponding cost function can be approximated by

Applying functional gradient descent (Eq. (7)) gives

4.5 Momentum acceleration (mPIE and mqNewton)

The momentum-accelerated ptychographic iterative engine (mPIE) [47] is the standard solver used for CP in PtyLab. In mPIE a predefined number of ePIE iterations $T$ is carried out, after which the search direction is complemented by a momentum term $\nu$ updating the entire object field of view $O_{\textrm {oFOV}}$,

Here $\eta$ is a damping term that is set to 0.7 by default [47]. Similar to conjugate gradient solvers [48,55,106], the momentum term accelerates the search direction and prevents zigzag motion towards the optimum. We emphasize that $O_{n+1,\textrm {oFOV}}$ in this subsection denotes the entire probe field of view, while in other subsection $O$ is an object box of the same size as the probe window.

While addition of momentum is typically done for another regularized version of the PIE-type family of algorithms (rPIE, see [47]), it can complement any of the existing reconstruction engines including PIE. To avoid naming ambiguities, addition of momentum to PIE will be referred to as momentum accelerated quasi-Newton (mqNewton), which we often use as an FP solver.

5. Robust inverse models

In the foregoing section, we have reviewed the basic inverse models underlying ptychography. However, oftentimes a variety of systematic errors are present in the experimental data that requires more robust inverse models. This is the case when the data is corrupted by for example partial spatial as well as partial temporal coherence, when the illumination wavefront profile is unstable throughout the scan, and when the scan positions are imprecisely known. In what follows, we discuss robust forward models that account for and mitigate the aforementioned sources of error.

5.1 Mixed States

In mixed state ptychography [45] the intensity in the detector plane is modeled as

The real space cost function for mixed state ptychography is

5.2 Multispectral ptychography

In multispectral ptychography [41,46,110–115] a light source of multiple individual spectral lines or a continuous spectrum is used. Because different colors are mutually incoherent, the detector update is identical to mixed state ptychography (compare Eq. (42)), but the index $k$ now denotes wavelength instead of spatial mode. The differences between mixed state and multispectral ptychography lie in the real space updates and in the propagator between the sample and the detector plane. In PtyLab we minimize the following cost function

These cost functions result in the updates

The latter update also results for the special case $\mu =0$. In this case the spectral channels are only coupled by the incoherent model for the detector intensity (Eq. (41)). In the presence of spectral regularization ($\mu \neq 0$) we do not need a priori knowledge about the spectral weights of the incident beam [41]. In the original work proposing multispectral ptychography no spectral regularization of adjacent channels was used. Instead the spectrum of the incident polychromatic probe was known a priori and used as a constraint in the optimization routine [46].

The second difference between mixed state and multispectral ptychography is the propagation model. Other code projects, for instance PtyPy [54], use zero padding to model the wavelength dependence of far-field wave propagation. In PtyLab the pixel size of a monochromatic wave can be scaled by using two-step propagators (scaledASP), which omits the need for spectrally dependent zero padding of the exit wave. For details the reader is referred to the supplementary information of [41].

5.3 Multislice ptychography (e3PIE)

In multislice CP [44] the specimen is modeled by a stack of 2D slices. The beam propagation method (BPM) [116] is used as a forward model for the exit wave. In each of the slices the thin element approximation is assumed to be valid. The cascade of multiplication with each slice and subsequent propagation enables the BPM to model multiple forward scattering effects. A basic version of multislice CP (termed e3PIE) is implemented in PtyLab. For details, the reader is referred to the original work by Maiden et al. [44] and subsequent work [117,118]. We have not yet implemented multislice FP [119] in the current version of PtyLab, although such an enigine may come in future releases.

5.4 Orthogonal probe relaxation (OPR)

A basic version of orthogonal probe relaxation (OPR) [51] is implemented in PtyLab. Instead of sharing the same probe across all scan positions in CP, as done for example in simple engines such as ePIE, OPR relaxes the requirement for a stable probe. The exit waves from different scan positions are used to estimate a low-rank basis, which seeks to model probe variations that occurred during a full scan, for example caused by pointing instability of the source. For details, the reader is referred to the original work by Odstrcil et al. [51]. We note that OPR can be combined with mixed states, as recently described by Eschen et al. [21].

With regard to FP, to our knowledge OPR has not been applied, although this may be an interesting approach to effectively model space-variant pupil functions. The latter is typically achieved by partitioning a larger field of view into a set of smaller sub-regions, each of which may be subject to different pupil aberrations. This approach, known as embedded pupil recovery (EPRY) [28], has to date essentially remained the only model for space-variant pupil aberrations in FP. However, because EPRY requires the reconstruction of a separate pupil for each sub-region, the model requires many degrees of freedom and ignores that adjacent sub-regions are unlikely to have strongly differing pupil aberrations. OPR could be a promising candidate to robustify EPRY in future FP applications for spatially varying aberration reconstruction. An example of the use of OPR in FP is shown in Fig. 7. The image FOV was split into small segments, each with its own unique pupil function and OPR was used to impose a low-rank consistency constraint on all the pupil functions. In applying OPR to FP adjacent pupils share information and poor convergence of some isolated field segments can be avoided (compare highlighted red boxes in Fig. 7). For further implementation details the reader is referred to [120].

 

Fig. 7. (a) Orthogonal probe relaxation scheme in FP. In step 1, the reconstructed pupils at a given iteration are factorized using SVD to produce an orthogonal full-field aberration basis. In step 2, the low-rank representation is obtained by eliminating modes with low-contributions to actual aberrations. This way, noise and errors are eliminated. In step 3, the low-rank full-field aberration reconstruction is performed from the low-rank basis. In this example, it can be seen that noisy pupils (top right corner) were replaced with a better pupil estimate. The whole process ensures that pupil aberrations remain well conditioned. (b) Experimental validation of pupil initialisation [120], where the use of OPR resulted in stable pupil reconstruction at the corners (indicated by red boxes).

Download Full Size | PDF

5.5 Subsampling (sPIE)

In some applications it is challenging to sufficiently sample the captured detector signal. For example, in EUV ptychography generating a highly focused probe is oftentimes restricted by the available hardware [19]. In other applications, such as near-field ptychography, the detector pixel size can be a limiting factor [121]. In both situations one may attempt to solve for the probe and object in CP (or the pupil and object spectrum in FP) from undersampled measurements, being too coarse to oversample the diffraction data. In principle, the detrimental effect of undersampling can be compensated by high overlap. In the context of CP this technique is known as reciprocal space upsampling and was first demonstrated by Batey et al. [122] (where the algorithm was named sPIE). In the latter work, the captured ptychography measurements were deliberately undersampled by means of binning, but the original oversampled data could be recovered thanks to the high scan grid overlap in the captured data. We later generalized this principle to arbitrary sensing matrices that are not necessarily a result of an operation equivalent to binning, but that could result from any sensing architecture that compresses multiple, not necessarily neighboring pixels into a smaller data cube [42]. In such situations one seeks to minimize the cost function

 

Fig. 8. Illustration of sPIE in CP. Assuming a far-field diffraction geometry, the real and reciprocal space sampling conditions are inversely proportional (indicated by the dashed lines): A small probe field of view (pFOV) in CP requires only coarse detector pixels $\Delta \boldsymbol {q}$. Conversely, if the physical probe wavefront extends over a larger region, the observed data is undersampled (top row). In such situations, the detector pixels need to be sub-divided into smaller pixels $\Delta \boldsymbol {q}'$ (bottom row). This allows to extend the numerical pFOV to be larger than the physical probe size. The resulting constraint in the forward model is that the sum over the intensities over a set $\boldsymbol {S}$ of sub-sampled pixels equals the corresponding observation over the same region in the observed data with a coarse sampling grid.

Download Full Size | PDF

5.6 Lateral position correction (pcPIE)

In ptychography the scan positions may not be accurately known, for example in the case of a low-precision scanning stage [68,86,123,124]. The wrong estimation of the scan positions will cause errors and artifacts during the stitching of the object patches into the large object field of view. In PtyLab a momentum-accelerated version of a cross-correlation-based lateral position correction algorithm is used [50]. The rationale of this position correction is based on the observation that, at iteration $n+1$ of the reconstruction procedure, the object patch estimate at scan position $j$ is slightly shifted towards its true position. This shift is detected and used to update the scan grid by maximizing the cross correlation

The updated scan position at iteration $n+1$ is

Default values are $\alpha = 250$, $\beta = 0.9$, and $d_0=0$.

5.7 Reflection ptychography with angle calibration (aPIE)

In reflection ptychography the sample plane and the detector are non-coplanar. Assuming far-field diffraction for simplicity, the captured data is related to the specimen exit wave by a Fourier transformation plus an additional coordinate transformation [17,125]. An inverse coordinate transformation can be applied to the captured raw data to simplify the forward model. However, this operation requires accurate knowledge of the angle between the optical axis and the specimen surface normal. If this angle is not calibrated within a fraction of a degree, the reconstruction quality can suffer notably. We have recently presented an algorithm for angular auto-calibration in reflection ptychography (aPIE), which is part of PtyLab. aPIE uses a heuristic strategy to estimate the unknown angle within an iteratively shrinking search interval. For details the reader is referred to [40].

5.8 Axial position correction (zPIE)

Similarly to pcPIE (position correction) and aPIE (angle correction in reflection ptychography) another self-calibration algorithm provided as part of PtyLab is zPIE, which can be used to estimate the sample-detector distance [39]. The main idea is that when the sample-detector distance is miscalibrated, the reconstructed object oftentimes exhibits characteristics of a slightly defocused inline hologram, including ringing at edges. An autofocus metric based on total variation (TV) is then used to calibrate the correct sample detector distance. We observed that TV-based autofocusing performs best in the near-field on binary specimens, although it can also be used on biological specimens. Other choices of autofocusing metrics can easily be implemented by the user, if the TV-based sharpness metric fails [126].

5.9 Ptychography combined with an external reference beam

We recently reported ptychographic optical coherence tomography, which combines full field frequency-domain OCT with ptychography [25]. In the latter work there was no need for an external reference wave, as common in OCT applications. Instead the reference was provided from a direct surface reflection of the sample itself. Thus the technique can principally be applied to the short-wavelength regime, where providing an external reference comes with extra experimental challenges. However, in the visible and infrared spectral range a reference wave is readily provided and can make POCT more convenient. Providing an external reference wave in ptychography requires adjustments to the forward model. In this case, we seek to minimize the cost function density

We note that the mathematical structure of external reference beam ptychography opens up a trivial ambiguity. Suppose that the triplet of probe $P$, object $O$, and reference $\tilde {\rho }$ yields the observed intenisty $I$, i.e.

Then it immediately follows that the triplet of probe $P$, object $-O$, and reference $\tilde {P}+\tilde {\rho }$ is also a solution, since

6. Scan grid optimization

In CP, a certain amount of consideration is needed to optimize the scan trajectory. To date, the majority of CP setups employ mechanical scanners, although variants exist where the beam is rapidly steered over the sample by means of galvo mirrors [127]. The latter offers advantages in terms of speed and overall cost of the experimental setup, but the isoplanatic illumination patch of such mirror systems is finite and thus limits the field of view over which the probe wavefront can be assumed to be stable, thus compromising one of the very benefits of CP. Hence mechanical scanners are still the preferred option. For such systems it is important to minimize the total scan distance in order to reduce scan time and prevent mechanical wear. In addition, unlike other scanning microscopy systems ptychography requires non-periodic scan grids to avoid ambiguities in the reconstruction [128]. A popular choice are Fermat scan grids [129], as depicted in Fig. 9(a). This type of scan grid is conveniently described in polar coordinates, where its trajectory assumes the form of a spiral. Minimizing to total travel path can be done using a solver for the traveling salesman problem (TSP). In PtyLab, we use a 2opt [130] TSP heuristic solver, which offers a good compromise between optimality and optimization time. Figure 9(b) shows an example of a distance-optimized scan trajectory, where the color scale indicates the start (blue) and end position (red).

 

Fig. 9. A variety of scan grids can be generated and optimized in PtyLab. (a) The typical workflow is to generate an aperiodic scan grid in polar coordinates, here a Fermat grid, and subsequently preprocessing steps on it. (b) The total path of the scan trajectory is minimized by solving the traveling salesman problem. In some cases, morphological operations to scan grids are useful, such as non-uniform scaling (c) and rectification (d). (e) Another useful technique is checkpointing, where the same scan point is revisited during a long scan. In panel (e) a Fermat grid with 200 scan points plus 20 checkpoints is shown, which are equally spaced in time. (f) For large scan grids an overlapping k-means (OKM) algorithm can be used to partition the scan grid into overlapping clusters and subsequently process each cluster separately. The overlap between clusters is required to synchronize phase information, which can otherwise differ by a global offset, and for stitching a large-field-of-view image.

Download Full Size | PDF

Moreover, several operations are available to transform scan grids, including non-uniform scaling and rectification as shown Fig. 9(c) and (d), respectively. The former allows for non-uniform spatial sampling, adjusted such that the sampling is higher in regions that are challenging to resolve, while the latter clips the field of view to a rectangular (here square) region. Another practically useful strategy is checkpointing (see Fig. 9(e)), which alters a given scan grid such that it revisits a certain reference point throughout the scan. Deviations in the diffraction data at the checkpoints allow for identifying sources of error in the experimental setup, including position drift, flux instability, and illumination wavefront variations. The checkpoints are equi-spaced in time.

The aforementioned techniques are primarily scan grid preprocessing techniques, meaning that the scan grid is optimized prior to the actual experiment. After the data acquisition, scan grid postprocessing techniques may be required. For example, large scan grids can be partitioned to prevent memory limitations. In this way large data sets may be chopped up into smaller pieces which are then individually processed [131]. It is important that the scan partitions spatially overlap, so that adjacent regions can be phase synchronized and stitched up post-reconstruction. To ensure overlap between the clusters, an overlapping k-means (OKM) algorithm can be used [132,133]. Figure 9(f) shows an example of a scan grid partitioned into four overlapping clusters (filled circles/triangles, unfilled diamonds/squares), each containing 160 scan points. In the middle, the clusters overlap. A second reason for scan grid partitioning can be to define batch gradients which speed up convergence and robustness [55].

As a note on FP scan grids, at first sight it appears surprising that the technique does not exhibit raster scan artefacts although the most commonly employed LED arrays are typically regularly spaced. However, the typically regular LED spatial arrangement still corresponds to a non-periodic spacing in angle, which explains that it is not subject to the aforementioned raster scan artefacts. While most of the aforementioned preprocessing steps are not required for FP due to the absence of mechanical movement, checkpointing and partitioning may still be used for monitoring stability and distributed data analysis, respectively.

7. Open experimental data and tutorials

We publish a variety of CP and FP data sets and tutorials with the aim to introduce users to the functionality of PtyLab. Figure 10 depicts two such data sets. The top row shows a soft x-ray ($\lambda =2.48$nm) data set collected at a synchrotron (experimental details in [76]). The bottom row depicts a visible light ($\lambda =625$nm) FP data set of lung carcinoma. For both data sets we show from left to right a single frame of the raw data, the recovered quantitative phase image (QPI) of the object (resolution test target and lung carcinoma histology slide), and reconstructed probe/pupil for the case of CP (top) and FP (bottom). Hue and brightness depict the phase and amplitude, respectively, of the complex-valued reconstructed quantities. A variety of additional data sets are published alongside PtyLab, which are summarized in Table 1. Each of these data sets comes with an online tutorial explaining suitable data analysis approaches, including self-calibration and regularization.

 

Fig. 10. Examples of CP and FP experimental data analyses using PtyLab. Top row: synchrotron-based soft x-ray CP (a) raw data, reconstructed (b) object QPI and (c) probe wavefront. Bottom row: visible light FP (d) raw data, reconstructed (e) object QPI and (f) pupil. Amplitude and phase are depicted as brightness and hue; experimental details in [76,134]. Figure adapted from [135].

Download Full Size | PDF

Table 1. Overview of open data sets and PtyLab tutorials. $\dagger$: previously unpublished.

View Table

8. Discussion and conclusion

PtyLab is a versatile ptychography software which we hope will aid researchers to explore the capabilities of CP and FP. Nevertheless, despite our excitement about this endeavor, we should mention some of its shortcomings: (1) Researchers with large-scale and high-throughput data analysis tasks (e.g. beamline scientists at synchrotrons) may be better off with one of the currently available high-performance ptychography packages mentioned in the introduction. However, we believe increased performance comes at the cost of flexibility in algorithm prototyping. (2) PtyLab currently does not support tomographic reconstruction with specimen rotation. It is to be noted that external CT toolboxes, such as Astra [138] or Tigre [139], can principally be used for ptychographic computed tomography once a sequence of 2D reconstructions at different angles is available. However, some ptychotomographic software embed specialized regularization techniques within the reconstruction routine [140], which are not available in standard CT packages. In contrast, ptychographic optical coherence tomography (POCT) [25] does not require angle diversity and 3D reconstructions can be obtained simply by performing a Fourier transform along the wavelength dimension in a multispectral reconstruction object stack. The latter is readily performed in PtyLab.

We have built PtyLab in an object-oriented and highly modular fashion. This enables the user to assemble reconstruction algorithms from elementary algorithmic building blocks (e.g. exit wave formation, propagation model, detector intensity and noise and model, partial coherence, experimental setup parameters), which interact seamlessly and concurrently. The parameters of a forward model can be modified during a reconstruction session.

We have designed PtyLab based on the principle of reciprocity. An interesting implication of the conversion between CP and FP is performance, as it provides the freedom to choose whether the computational complexity of the Fourier transform operation, used in typical inversion algorithms, scales logarithmically with the number of pixels in the detector or with the number of scan positions in a given data cube - numbers which can be orders of magnitude apart so that even on a logarithmic scale, practical speed ups can be achieved. However, in order to take full advantage of reciprocity, interpolation techniques to non-equidistantly sampled geometries are required, which will be explored in future work.

In summary, we have presented PtyLab, a cross-platform, open source inverse modeling toolbox for CP and FP. We believe PtyLab’s major strengths lie in (1) the uniform framework for CP and FP enabling cross-pollination between the two domains, (2) the availability in three widely used programming languages (Matlab, Python, and Julia), making it easy for researchers with different programming backgrounds to exchange and benchmark code snippets and data anlyses, and (3) its versatile code architecture suited both for beginners and experts interested in rapid ptychographic algorithm prototyping. In addition, a plethora of self-calibration features (e.g. aPIE, zPIE) and algorithmic novelties (e.g. conversion between CP and FP, POCT, CP with external reference beam, sPIE) are available that to our knowledge have previously not been featured in open access ptychography code. Various functions for scan grid generation help the user to optimize data acquisition and postprocessing. For further information the reader is referred to the GitHub website with its accompanying tutorials as well as the open data provided along with it [38].