Three-Dimensional Single-Shot Ptychography

Here we introduce three-dimensional single-shot ptychography (3DSSP). 3DSSP leverages an additional constraint unique to the single-shot geometry to deconvolve multiple 2D planes of a 3D object. Numeric simulations and analytic calculations demonstrate that 3DSSP reconstructs multiple planes in an extended 3D object with a minimum separation consistent with the depth of field for a conventional microscope. We experimentally demonstrate 3DSSP by reconstructing orthogonal hair strands axially separated by 5 mm. Three-dimensional single-shot ptychography provides a pathway towards volumetric imaging of dynamically evolving systems on ultrafast timescales.


Introduction
Coherent Diffractive Imaging (CDI) is a computational microscopy technique that uses a sophisticated algorithm to processes far-field intensity measurements simultaneously returning both the phase and amplitude of the probing illumination and the specimen [1,2,3]. The technique is particularly exciting because it can provide simultaneous amplitude and phase contrast imaging across the EM spectrum in both reflection and transmission modalities [4,5,6]. In ptychography, an advanced implementation of CDI, a sample is scanned transverse to a probe illumination producing diffraction patterns from overlapping regions of illumination. Transverse scanning of a probe illumination limits the minimum collection time but increases the technique's robustness and removes the need for a priori knowledge of the specimen [7,8]. Recently developed single-shot ptychography (SSP) eliminates the need for scanning by introducing a diffractive optical element (DOE), often an array of pinholes, and a 4f imaging system [9,10]. A ray tracing schematic of a conventional SSP imaging system is shown in the bottom panel of Fig. 1. The DOE breaks up incoming illumination into several beamlets. The 4f system collimates the beamlets and then crosses them. The object is slightly offset from the crossover point (which is also the Fourier plane of the 4f imaging system) and the detector is segmented such that the diffraction pattern produced by each beamlet is essentially recorded independently. The system thus provides the same set of diffraction patterns as scanning ptychography but in a single shot. SSP dramatically reduces the acquisition time of a ptychographic dataset and heralds the possibility of time-resolved, simultaneous phase and amplitude contrast imaging, as recently demonstrated experimentally [11,12].
Generally, ptychography depends on the projection approximation, which models the interaction of an illuminating probe with a specimen's complex transmission function as a product. This approximation is accurate in the limit of an optically thin object [13], but it restricts the technique to imaging effectively two-dimensional (2D) samples. Scanning ptychography has been adapted to image in three-dimensions (3D) using a multi-slice approach (3PIE) and ptychographic tomography [14,15,16]. These techniques were applied to image thick samples like biological tissues with relatively simple tabletop experimental systems. While these methods are powerful and robust, their required arXiv:2002.06106v1 [eess.IV] 14 Feb 2020 acquisition times preclude them from measuring dynamic objects. As we probe more complex phenomena at ever shorter time scales the need for metrologies capable of imaging dynamically evolving objects becomes critical. Here we introduce Three-dimensional Single-Shot Ptychography (3DSSP), which uses the same SSP experimental configuration but implements a novel algorithm that leverages an additional object domain constraint to reconstruct multiple 2D planes of a 3D object. Numerical simulations show that 3DSSP can reconstruct multiple planes in an extended 3D object with minimum axial separation consistent with the depth of field for a conventional microscope. We experimentally demonstrate 3DSSP by reconstructing orthogonal hair strands axially separated by 5 mm. This novel imaging technique provides a pathway toward volumetric imaging of dynamically evolving systems on ultrafast timescales, such as the nonlinear response of materials [17,18,19].

Methods
Current single-shot ptychographic experiments depend on reconstruction algorithms developed for conventional ptychography like ePIE [8,10]. Conventional reconstruction algorithms provide sufficient results when applied to 2D single-shot experiments as long as the object is optically thin [13], but thick objects require more sophisticated multi-slice reconstruction algorithms. However, the geometry of the single-shot system prohibits direct application of conventional multi-slice algorithms like 3PIE [16]. In particular, the beamlets are collimated and propagate at an angle through the object. Our 3DSSP technique reconstructs 3D objects based on recorded diffraction patterns from a SSP experiment by taking into account the SSP geometry and applying a multi-slice approach. 3PIE was inspired by electron microscopy techniques [16], and it has been successfully applied to reconstruct phase and amplitude contrast images of continuous 3D objects [15]. Like 3PIE, 3DSSP computes exit waves for each 2D slice (or z-slice) of a 3D object, but where 3PIE propagates the exit wave from one slice to form the illumination function at the next slice, 3DSSP propagates and shifts each exit wave to account for the geometry of a SSP setup.

Inter-slice Propagation
In SSP an object is placed a distance δ away from the crossover point. Beamlets created by a DOE are collimated, crossed, and then separate geometrically, as shown in Fig. 1. A critical input to any ptychographic reconstruction algorithm is the position of the probes on the object. In the single-shot geometry, these positions change as a function of δ. To account for transverse shifts between z-slices, we define the following Inter-Slice Propagator (ISP) as: (1) Here F{·} is the two-dimensional discrete Fourier transform, ψ es is the exit surface wave at slice s = {1, ..., N s }, N s is the total number of slices, ψ is+1 is the incident surface wave at slice s + 1, f ≡ (f x , f y ) are spatial frequency grids for the object domain grids r ≡ (x, y), k 0 is the vacuum angular wave number, and the transfer function H is given by: The transfer function in Eqn. 2 is simply a product of the free-space transfer function [20] with a linear phase that accounts for the shift between planes including sub-pixel contributions. The magnitudes of the shifts are give by, ∆x = ∆zX/f and ∆y = ∆zY /f . Where X and Y are the position of a beamlet on the detector, ∆z is the propagation distance between z-slices for a specific beamlet, which in the small angle approximation is the same for all beamlets, and f is the focal length of the second lens. The positions of beamlets at the first object slice are determined from an image of the DOE, δ, and the focal length of the second lens [10]. The inverse inter-slice propagator is then defined as: where * is the complex conjugate. Accounting for the position change of the beamlets at each object slice enables 3D reconstructions via a multi-slice approach for the single-shot geometry. In comparison to 3PIE, 3DSSP places an additional constraint on the problem because the positions at each slice shift by a known amount. However, because the beamlets separate from each other in the SSP geometry, the total 3D volume that can be imaged by 3DSSP is limited.

The 3DSSP Algorithm
In this section, we describe one iteration of the algorithm for a single beamlet and highlight deviations from other imaging methods that make three-dimensional imaging possible in the single-shot geometry. A flow diagram of the 3DSSP algorithm is shown Fig. 2 for reference. For a full reconstruction, the steps presented below are repeated for each beamlet and then multiple overall iterations until some prescribed number of total iterations are completed or until the error reduces below some threshold [21,22]. The 3DSSP algorithm forms the exit wave at the first z-slice as the product between the incident probe and first slice of the object as ψ e1 = p o 1 . We initialize the algorithm with a guess for the probe by taking the average of the inverse Fourier transforms of all of the segmented DOE beamlet images without the object present. The object is initialized as free space at each slice. The exit waves at subsequent z-slices are then calculated using: with ψ i1 ≡ p is the input probe illumination. This exit wave is then propagated and shifted to form the incident wave for the next slice as: Eqns. 4 and 5 are applied repetitively until the exit wave for the last slice is formed i.e., s = N s . The exit wave at the last slice is updated with the measured data according to: where ψ e Ns is the updated exit wave at slice N s ,ψ e Ns = F{ψ e Ns } and I is the measured diffraction pattern intensity from the section of the segmented detector that corresponds to the current beamlet. The modulus constraint in Eqn. 6 is a common feature of all CDI algorithms however, in 3DSSP the measured intensity for a given position is a segment of the full detector. The exit waves are back propagated through the object using the inverse ISP in Eqn. 3 and the object and incident waves are updated using the corresponding exit surface wave according to Eqn. 3 from [16].

Axial Resolution and Reconstruction Volume
An important parameter in any 3D imaging system is the axial resolution. The 3DSSP algorithm separates planes in the object because the position of the beamlets change at each slice. Thus, the minimum resolvable separation of planes should be achieved when the position of the outer most beamlet shifts by one transverse resolution unit between two planes separated by the axial resolution, δz. Based on this premise we calculate the expected axial resolution as: where λ 0 is the vacuum wavelength, f is the focal length of the second lens in the 4f imaging system, and X det , X dif are the size of the full detector and the size of one recorded diffraction pattern, respectively. Recall that these sizes are not the same in the single-shot geometry because the detector is segmented to record diffraction patterns independently. Using the parameters from our experimental single-shot system: f = 5 cm, λ 0 = 532 nm, X det = 10.85 mm, X dif = 1.59 mm, we calculate an axial resolution of δz ∼ 150 µm. It is worth noting that if X det = X dif Eqn. 7 gives the paraxial depth of field for a microscope with magnification much greater than one.
We now calculate a value for the imaging volume in 3DSSP by first noting that the overlap fraction between beamlets, η ≡ A/πr 2 , changes as a function of the distance from the crossover point, z. Here, A is the overlap area shared between two beamlets and r is the radius of each beamlet. It is common in ptychography to use overlap fractions between η = 0.9 and η = 0.6, which we use to set the minimum and maximum z values, respectively. Our task is to write an expression, z(η), that provides bounds on the axial extent of the imaging volume.
The specific relationship for z(η) depends on the shape of DOE elements and their relative arrangement. Here we derive the relation for a DOE similar to our experiment which is an array of circular pinholes of radius, r p ∼ 32µm, arranged in a Fermat spiral pattern [23]. Beamlets passing between the two lenses in the single-shot configuration have a cross-section given by the 2D Fourier transform of a DOE element. This means that the shape of each beamlet when passing between the 4f lenses is the Airy pattern [24], which has a central Airy disk that we take to be the illumination function. The Airy disk is surrounded by a circular boundary where the field amplitude drops to zero, which defines the beam radius. The area of overlap between two Airy disks may be approximated for small center-to-center distances, d as: A ≈ πr 2 − 2rd. In this approximation, the value for the area becomes more accurate with increasing overlap and reaches 10% error when η ∼ 0.35. The Vogel model for the Fermat spiral gives the radial distance to each pinhole as X det /2 × n/N p [25], where n ∈ {0, ..., N p − 1} is the index number of each pinhole and N p is the total number of pinholes. The angular position of each pinhole is given by n θ g , where 2π/θ g = ϕ, is the golden ratio.
We approximate the distance between nearest neighbor pinholes as the distance between the n = 0 and n = 1 pinhole. This approximation is only valid near the center of the Fermat pattern but this is the region where the reconstruction will have the highest fidelity in any case. The center-to-center distance between Airy disks as they propagate from the crossover point to the second lens is then given by Writing the approximate overlap area as A = ηπr 2 , substituting d, and solving for z gives: where we have used r = 0.61λ 0 f /r p to relate the radius of the Airy disk to the pinhole radius. Using Eqn. 8, we calculate an axial extent, z(η = 0.6) − z(η = 0.9) ∼ 1.4 cm. This extent sets the height of a truncated cone that fills a volume of ∼ 75 mm 3 .

Simulated Experiments
In the simulations presented here we calculate ideal diffraction patterns using parameters similar to our experimental system and the forward model in Eqns. 4 and 5. The simulated 4f system has unit magnification with f = 5 cm focal lengths for each lens. The detector has 2048×2048 square pixels with a dX = 5.3 µm pitch and the input wavelength was λ 0 = 532 nm. This produces a collimated beam diameter between the lenses given by D ∼ 0.5 mm. With these parameters the oversampling is calculated as σ = f λ 0 /DdX, which corresponds to a value of 5. In addition to satisfying oversampling, 3DSSP requires the distance from the crossover point to each slice in the reconstruction. In the simulated and experimental results presented here, we specify the distance to each object slice as the distance from the crossover point to the first slice, δ, plus the distance between slices, ∆z.
The object was scaled to fit within the theoretical field of view at each plane. Then we simulated the input illumination using parameters of our experimental DOE and 4f imaging system. For the simulations presented here, the DOE is a Fermat spiral of N p = 40 identical pinholes. A Fermat spiral DOE lends itself nicely to SSP experiments because it gives the most uniformly overlapped probes on the object [23]. We Fourier transform the DOE to calculate the total illumination at the crossover point and then use the free space transfer function to propagate the illumination a distance δ to the first object slice.
We model the interaction of this illumination with the first 2D slice of the object by multiplying it by the complex transfer function of that slice, and define the result as the first exit wave. Then we propagate the first exit wave to the second slice and multiply it by the complex transfer function of that slice. We repeat this process N s times to calculate the final exit wave. A two-dimensional discrete Fourier transform of the final exit wave gives a simulated diffraction pattern. Finally, the detector is segmented into N p sections by applying a centroidal Voronoi tessellation, which gives the maximum area for each beamlet's diffraction pattern. Together the N p simulated diffraction patterns make a ptychographic data set, which we feed into the 3DSSP reconstruction algorithm. Fig. 1 presents the results from our first simulated experiment where the 3D object is a continuous broken loop structure. We simulated and reconstructed three slices taken from the continuous 3D object each separated by 3 mm. The 3DSSP algorithm successfully reproduced the slices of the object. This object is particularly interesting because it is inherently anisotropic and standard 2D projection-based imaging techniques like Schlieren imaging or interferometric imaging would have incorrectly produced an image of a ring. Attempts to use Abel inversion to produce a 3D image from 2D projections would also fail.
The first simulation shows the power of 3DSSP as the algorithm successfully reconstructs the object where other techniques fail. However, the object in that simulation is relatively simple as it does not contain fine details and its size does not change significantly. To investigate the quality of reconstructions from the 3DSSP algorithm, we performed a second simulated experiment in which the object was two slices from a MATLAB example MRI scan separated by 1.5 mm. We chose the MRI scan slices as an object because they contain small features and the size of the images continuously vary. The two slices that we used to simulate diffraction patterns are shown in a) and b) of Fig. 3 above, while the reconstructed slices are shown in c) and d). The left panel of Fig. 3 shows the reconstructed slices overlaid on the head. While the scale of a human head is much larger than the field of view of our current imaging system, this simulation shows that the 3DSSP algorithm is capable of capturing fine details from continuously varying 3D samples.
To test the axial resolution in Eqn. 7, we performed a series of simulated experiments. We simulated two slices of another example MRI data set from MATLAB and varied the separation between the slices. To quantify the accuracy of the reconstructions we developed two error metrics. The object error complement (OEC) is calculated as one minus the RMS difference between the absolute value of each pixel in each slice of the simulated and reconstructed object. Similarly, the diffraction error complement (DEC) is calculated as one minus the RMS difference between the absolute value of the intensity pattern produced by the simulated and reconstructed slices. The results of these simulations are shown in Fig. 4. They show that our calculated axial resolution is accurate because both error metrics exhibit dramatically different behavior above and below 150 µm.
The OEC is effectively constant until the separation between the slices reduces below the axial resolution, where it decreases monotonically. The diffraction error compliment data presented in the top right panel of Fig. 4 is normalized such that 1 is the maximum DEC and 0 is the minimum. Like the OEC, the DEC decreases dramatically as the slice separation decreases below then axial resolution limit. However, as that separation continues to decrease the diffraction error compliment increases again. This unexpected result can be explained by considering that once the separation between the slices is so small that the propagation of light between them becomes negligible, the diffraction pattern from two separated slices is effectively the same as the diffraction pattern from one slice comprised of the product of the two individual slices. The reconstruction of two slices with no separation (∆z = 0) is shown in Fig. 4. There are no discernible details of either slice in the reconstructions, but the product of the reconstructed slices matches the product of the ground truth slices. This corroborates the above explanation of the behavior of the DEC below the axial resolution limit, and justifies the claim that the microscope is diffraction limited in the axial direction. Moreover, the simulated axial resolution matches the theoretical value, which gives confidence in Eqn. 7. There are a few parameters in Eqn. 7 that can be optimized for better resolution in the single-shot system. For instance, we could theoretically push the resolution of the system to less than 10 µm by using lenses in the 4f imaging system with half the focal length and doubling the size of the detector, which would also double the size of the diffraction patterns. Making these changes would give axial resolution 16 times smaller than the current system achieves. Axial resolution on the order of 10µm is ideal for many applications, including plasma imaging.

Experimental Results
We experimentally verified 3DSSP using our existing SSP system. The DOE is a Fermat spiral of 40 circular pinholes of 55 µm diameter which is illuminated by a large isotropic 532 nm collimated beam. The 4f imaging system has unit magnification and f = 5 cm lenses, and the detector is a Thorlabs 8051M-USB 8 megapixel monochrome CCD camera with the face-plate and wedge window removed. The detector has 3296x2472 square pixels with a dX = 5.5 µm pitch. The DOE slightly under-filled the detector such that it only has a square effective area of about 2000 pixels. The experimental setup is shown on the left of Fig. 5.
To create an object for testing 3DSSP, we used two orthogonally oriented 50 µm thick strands of hair axially separated by 5 mm in a cross-like pattern. The discrete two plane nature of the object was chosen because it makes identification of a successful reconstruction trivial. We determined the crossover point and mounted the object on a micrometer such that the first hair was at the crossover plane. Then we used the micrometer to move the object such that the first hair was 5 cm downstream of the crossover point, δ = 5 cm, where the overlap fraction calculated by Eqn. 8 is η = 0.66.
Diffraction patterns were collected using multiple exposure times that were stitched together using a high dynamic range (HDR) algorithm, which improved the fidelity of the reconstructions for this proof of principle. After collecting the diffraction pattern data, we applied the centroidal Voronoi tessellation to break up the detector into 40 segments of 220x220 pixels. We calculated the initial positions of the beamlets on the first object slice. These were determined by removing the object, imaging the DOE directly to determine beamlet positions on the detector, and then calculating the transverse shifts for each beamlet as described above.
The diffraction patterns and calculated positions on the first slice were fed into the 3DSSP algorithm and the reconstructions are shown on the right of Fig. 5. As shown in the figure, the algorithm successfully deconvolves the first slice which shows the vertical bar in the figure form the second slice which shows a horizontal bar. The modulations in the reconstruction are a result of imperfect knowledge of the distance between strands, δ, and a limited number of iterations. The reconstruction of this relatively simple object demonstrates that experimental 3DSSP can successfully deconvolve multiple 2D slices of a 3D object. Future work will expand on this experimental proof of concept by increasing the number of reconstructed slices, determining the resolution of the system experimentally and establishing the minimum transmission and/or maximum extent of an object that can be imaged.

Conclusion
By exploiting a feature of the SSP setup, that the position of the beamlets change as a function of distance from the crossover point, we developed a novel algorithm for three-dimensional single-shot ptychographic imaging. While 3D ptychography techniques exist, the geometry of single-shot ptychography setups lends itself well to create highly constrained systems for reconstructing objects at multiple z-slices, and allows for time resolved imaging. Through a number of simulations, we showed the power of the 3DSSP algorithm to image 3D objects with high reconstruction quality. The axial resolution of this technique was derived and found to be consistent with the depth of field of a standard microscope. Simulations support the axial resolution of 150 µm and a reconstruction volume of 75 mm 3 given our setup parameters. The axial resolution can be improved down to around 10 µm with reasonable modifications to the setup. The 3DSSP technique was experimentally realized by imaging orthogonally oriented hairs separated by 5 mm. The reconstructions of each hair were deconvolved as predicted by our simulations. Through both simulated and experimental reconstructions, we have shown that 3DSSP is a reliable method for imaging a multitude of 3D objects. This is currently the only method to image transient, non-reproducible objects in three-dimensions in both phase and amplitude and with temporal resolution set by the illumination. Metrologies such as this will be imperative for studying phenomena such as dynamically evolving plasma.

Funding Information
The authors gratefully acknowledge funding from the Air Force through AFOSR FA9550-18-1-0089 and Los Alamos National Laboratory through contract number 501188. consists of incoming plane-wave illumination, a DOE, 4F imaging system and detector. The object is placed a distance δ away from the crossover point, and the inter-plane spacing is ∆z. The detector is segmented to record individual diffraction patterns, which are fed into the 3DSSP algorithm along with the initial positions of the beamlets on the first object slice. The middle panel of the figure shows a simulated example reconstruction of three planes of a continuous 3D object. The reconstructed object planes match up well with the continuous 3D object. Interestingly, the broken loop object that is simulated here is not recoverable by standard 2D projection-based imaging techniques.  show the simulated slices used to produce ideal diffraction patterns while c) and d) are the corresponding reconstructions. Note that the fine details of the simulated slices are well represented by the reconstructed slices. The left panel shows the reconstructed slices superimposed on the head geometry. The reconstructed slices fit as expected showing that the reconstructed slices also accurately account for the changing size of each slices. This promising result suggests that the 3DSSP technique can be applied to continuously varying 3D objects with fine detail and good reconstruction quality. Figure 4: Results of the simulated experiments to determine the axial resolution of the 3DSPP system. The final OEC as a function of the separation between the two slices is shown in the top left graph, while the final normalized DEC is shown in the top right. These error metrics are calculated as described in the text and the dashed red line in each represents the calculated axial resolution, 150 µm. The images below show the slices used to simulate ideal diffraction patterns on the left, which we call the ground truth, and the resultant reconstructed slices at the specified axial separations. These results were selected because they exemplify the general behavior of the reconstructions at 0 separation, just above and well above the axial resolution. The bottom right most panels show the product of the recosntructed slices for ∆z = 0 (above), and the product of the ground truth slices (below). On the left, a diagram shows the experimental setup of the 3DSSP system. The lenses (5cm focal length) make a 4-f imaging system that images the DOE (pinhole pattern) to the detector. The object is comprised of 2 hairs that are orthogonally oriented and separated by 5mm. The first hair was placed 5mm out of the focus of the first lens. Images of the reconstructions of the two slices are shown on the right. The top image shows the amplitude (brightness) and phase (color) of the first z-slice, while the bottom image shows the amplitude and phase of the second z-slice. The algorithm was able to deconvolve the two slices such that the first slice only shows a vertical hair and the second slice only shows the horizontal hair.