3D x-ray imaging of continuous objects beyond the depth of focus limit

1. INTRODUCTION

Over the entire span of the electromagnetic spectrum, x-rays offer a unique combination of nanometer wavelength to enable high-spatial-resolution imaging, large penetration in millimeter-scale specimens, and low values of plural and inelastic scattering to enable straightforward image interpretation and quantization [1,2]. A variety of lens-based and lensless imaging methods have demonstrated 2D spatial resolution better than 20 nm [3,4]. Among these, ptychography [5,6] offers a unique combination of having a spatial resolution determined not by optics but by the maximum detected scattering angle, while displaying robustness in phase retrieval for non-isolated objects [7]. In ptychography, a finite-sized x-ray probe (coherent beam spot) is used to illuminate the specimen at multiple overlapping probe positions, while recording the far-field diffraction intensity corresponding to each probe position. The resulting data redundancy allows one to recover both the object and the probe. One can even reconstruct multiple probe function modes to account for x-ray beam partial coherence [8], sample vibration [9], and continuously moving illumination [10–13]. Ptychography has been used to image thin circuit layers through 300 μm of silicon at 12 nm resolution [14]; sub-10 nm resolution has been achieved with thinner specimens [15–17], and sub-wavelength resolution has been obtained using extended ultraviolet (EUV) light [18].

As the transverse resolution ${\delta }_r$ of x-ray imaging continues to be improved, a challenge lies ahead: the invalidity of the pure projection approximation. When using a circular lens, the depth of focus (DOF) of an image is given by [19]

Normally ptychography reconstructs a single plane, which is the exit wave leaving a 2D object or the pure projection of complex optical modulation by a 3D object. With pure projections of 3D objects, single-slice ptychographic tomography (SSPT) allows one to recover the 3D object by first applying phase unwrapping [21] to the projection images and then using these projections in standard tomographic reconstruction algorithms [22,23]. However, one can deal with objects located at multiple planes along the x-ray beam direction by applying the object’s optical modulation at each plane and propagating the resulting wavefield to the plane of the next object slice before propagating the final exit wave to the detector plane [24]. In the reconstruction process, the probe and object functions are updated at each plane, since the ptychographic update steps work to maximize the separability of object and probe. Using this multislice approach, images have been obtained of a few discrete, separated planes by using a single viewing direction [20,25,26–28] or a limited range of tilt angles [29].

One approach to 3D imaging that builds on the above work was recently demonstrated [30]. This approach, termed multislice ptychographic tomography (MSPT), involved the use of multislice ptychography to reconstruct five depth planes of a 3D object at each viewing direction (with a span over all planes sufficient to encompass the 3D object and a spacing of planes fine enough so that one can ignore to some extent Fresnel diffraction within each plane). Since the object showed primarily phase contrast, the phase projections of the sample at each angle were calculated by pixelwise addition of the phases of the five slices. This again required a phase unwrapping process to yield a pure projection image for standard tomographic reconstruction.

We demonstrate here a different approach to the reconstruction of 3D images of extended, complex objects. Rather than seeking the solution of several object planes from each viewing direction separately before combining calculated pure projections in a standard tomographic approach, we consider the totality of the 3D object in each update step. To do so, we simulate the propagation of probe function illumination waves at various probe positions and incident angles through a present guess of the 3D object so as to produce a set of assumed intensity patterns to be recorded on a far-field detector; we then adjust the guess of the object so as to minimize the difference between the actual detector-plane Fourier magnitudes against our present guess of the same. Such an approach has been used with learning algorithms to guide the object updates [31,32] as well as with the imposition of a sparsity regularizer [33]. In our case, we use a newly developed proximal alternating linearized minimization algorithm for finding the object, as will be discussed below. Our algorithmic approach exploits a parallelized implementation that can address the additional nonlinearities and computational complexity introduced with the fully propagated model and the capabilities of high-performance computing. Our approach requires no phase unwrapping because the phase shift per 3D voxel is always small and because it is based on a forward model rather than backwards propagation of the wavefield through the object. We term our approach multislice optimized object recovery, or MOOR.

The use of multislice propagation to carry out the forward model calculation is well established. First introduced to interpret high-resolution electron microscopy data [34,35], the multislice method has been shown [36] to recreate a wide range of x-ray optical phenomena relevant to nanoscale imaging. These include grazing incidence reflectivity and wave propagation in arbitrarily thick transmission gratings that were previously understood only by using coupled-wave theories for simple, mathematically definable structures [37–40]. Objects that have refractive index boundaries within gridded voxels can be treated by filling the voxel with a weighted sum of materials [36]. Thus, multislice propagation releases one from the limit of considering only those objects for which the Born approximation applies, or objects that are constrained to be within the effective DOF of the imaging scheme employed.

2. BEYOND THE PURE PROJECTION APPROXIMATION

The Rayleigh resolution criterion [41] uses the position of the first minimum of the Airy intensity distribution as the measure of the transverse resolution ${\delta }_r$, giving us

Ptychography currently allows for the solution of the phase problem from a thin object or from discrete planes. However, the nature of diffraction (that it is jointly contributed to by materials throughout the entire depth of the sample) implies that the exiting wave contains 3D information about the whole object, which inspires an iterative optimization-based solution to the object unknowns. This strategy requires a forward model to propagate the probe wave through the object in each iteration, which can be implemented by using a multislice approach. Multislice propagation decomposes the object into a number of thin layers with thickness $\mathrm{\Delta }z$ along the beam axis. For the $j$th layer, the incident wavefront ${\psi }_j$ is modified by the well-tabulated [43] refractive index $n=1-\delta -i\beta$ of the material in that layer as

3. OPTIMIZATION METHODOLOGY

We discretize the object domain with a 3D regular grid of dimensions $n\times n\times n$ and discretize the probe domain with a 2D regular grid of dimensions $l\times l$; the real-valued measurements are assumed to lie on a 2D regular grid of size $l\times l$. For ease of representation, we employ vector notation and denote the object variable by $\mathbf{y}\in {\mathbb{C}}^{n^3}$, the probe variable by $\mathbf{x}\in {\mathbb{C}}^{l^2}$, and each of $m$ scan data by ${\mathbf{d}}_j\in {\mathbb{R}}_{+}^{l^2}$. The variable $\mathbf{y}$ is related to $\beta$ and $\delta$ by

We let $S_{q,j}\in {\mathbb{R}}^{n^3\times l^2}$ denote the linear operator that samples the $q$th propagation plane in the $j$th scan of the object. In this notation, the forward-model-generating scan $j$ is

We perform the inversion by solving the constrained minimization problem

The first term of the objective function in Eq. (11) of

Because of the nonlinearity inherent to the forward model in Eq. (10), updating the step size parameters in each iteration (e.g., by a line search procedure like in our implementation) can significantly benefit overall performance. The convergence of Algorithm 1 can be further accelerated using an inertial version of the proximal alternating linearized minimization (iPALM) [47], which our implementation exploits. The evaluation of the analytical partial gradients is the main computational bottleneck of Algorithm 1; however, significant computational savings are possible by careful ordering of the terms. Furthermore, the independence of scans, as well as the spatial structure of probe and object functions, can be exploited in a parallel implementation of Algorithm 1 to enable a large-scale solution such as that considered in our demonstration.

Algorithm 1. A proximal alternating linearized minimization algorithm for solving Eq. (11)

View Table

4. DEMONSTRATION

Figure 2 illustrates the setup used for our simulated experiment. One can increase the throughput of x-ray ptychography by using small focused beams and small-pixel-count detector arrays with a high frame rate [48], so this is the configuration we selected for our simulation. We simulated a Gaussian focus spot of 5 keV or $\lambda =0.25\mathrm{nm}$ X rays with a 14 nm full width at half-maxiumim (FWHM) probe size. This was represented with a $l\times l=72\times 72$ pixel array with 1 nm pixel size. We considered an object described below (200 nm in extent, sampled on a 3D grid with $n\times n\times n=256\times 256\times 256$) that is large enough to go beyond the pure projection approximation at the chosen spatial resolution (1 nm voxel size, since sub-wavelength resolution imaging has been demonstrated with ptychography [18]), forcing us to account for beam propagation effects within the specimen. In other words, the “depth of focus” of diffraction effects within the object is found from Eq. (1) to be $5.4{(1\mathrm{nm})}^2/(0.25\mathrm{nm})=22\mathrm{nm}$, while the probe has a DOF of 4,300 nm. With 1 nm voxel size, we were able to use nearest-neighbor sampling rather than interpolation for tomographic object rotation. Assuming Nyquist sampling of discrete Fourier transforms, a $l\times l=72\times 72$ pixel detector array was assumed to collect a scattered signal going out to a maximum spatial frequency of $1/(2·1\mathrm{nm})=500{\mathrm{\mu m}}^{-1}$.

 

Fig. 1. Depth of focus as a function of x-ray photon energy for a variety of transverse resolution ${\delta }_r$ values. Also shown is the $\exp [-1]$ penetration depth ${\mu }^{-1}$ for amorphous ice (for frozen hydrated biological specimens) and silicon (for microelectronics specimens) as proxies for the thickness range of x-ray imaging as a function of photon energy. As the transverse resolution ${\delta }_r$ in x-ray microscopy is improved to finer values, the DOF decreases with the square of the resolution improvement [Eq. (1)], leading to a decrease in the size of a specimen that can be imaged within the projection approximation required by standard tomography.

Download Full Size | PPT Slide | PDF

 

Fig. 2. Experimental geometry used for our simulated experiment. A lens was assumed to produce a Gaussian coherent illumination probe of size 14 nm FWHM through which the object is scanned at each object rotation angle. A far-field diffraction pattern is then captured in each scan.

Download Full Size | PPT Slide | PDF

Figure 3 shows the simulated object shape and the incident illumination wavefront. The focus spot is considered to be at the center of the object, so that the 2D illumination function shown in Fig. 3 is the exit wave of the probe backpropagated by half of the 3D object size, although this effect is small given the relatively large DOF of the probe compared to the object size.

 

Fig. 3. Left: an isosurface rendering of the true object with 200 nm length along the cone’s axis. The simulated object was a conical glass capillary tube with embedded Ti nanospheres. Right: the as-designed Gaussian probe function with 14 nm FWHM size, as represented at the midpoint of the object region; brightness indicates the amplitude of the wave, and hue indicates the phase (see color wheel inset).

Download Full Size | PPT Slide | PDF

The object was designed to resemble the thin capillary tubes that are widely used in laboratories as a carrier for cryo-microscopy of cells in liquid suspension [49] and for applications such as electrophysiological probing or photocatalysis [50]. Composite structures like these are interesting subjects for x-ray 3D imaging because of the presence of features ranging from nanospheres with diameters of tens of nanometers to quantum dots several nanometers in size, close to the wavelength of a soft x-ray probe beam. Further challenges in imaging are contributed by the capillary tube substrate, whose diameter can be orders of magnitude higher than that of the finer features, rendering the projection approximation no longer valid. The object was designed by using the open-source package XDesign [51], using x-ray refractive index data retrieved from a database within Xraylib [52]. It is composed of the tip region of a capillary tube made of Si (with $\delta =1.98\times {10}^{-5}$ and $\beta =1.13\times {10}^{-6}$ at 5 keV in the complex refractive index $n=1-\delta -i\beta$). Nanospheres of ${\mathrm{TiO}}_2$ ($\delta =3.01\times {10}^{-5}$, $\beta =3.58\times {10}^{-6}$ at 5 keV), with diameter ranging from 8 to 20 nm, were computationally “deposited” on the outer surface. The cone-shaped capillary tip has a top diameter of 20 nm, a bottom diameter of 160 nm, and a length of 200 nm. The tube wall thickness varies linearly from 5 nm at the top to 15 nm at the bottom.

This demonstration involved a complex object array with $n^3={256}^3$ voxels, which implies that we are inverting for approximately 17 million complex variables. We simulated the data acquisition scans in the following way. First, data from ${23}^2$ probe positions from a single large plane were generated. The center of each probe position on the plane was at $[12j,12k]$ for $j,k\in 0,\dots ,22$. We tested two rotation sampling schemes: one with 90 specimen rotations at 4° increments about a single axis over an 360° range, and another with 180 rotations at 2° increments over a 180° range. The total number of far-field diffraction patterns recorded for both sampling schemes was thus $m=90·23·23=47610$. We note that the 2° rotation increment is not fine enough to meet the Crowther criterion for fine angular spacing in standard tomography [53]. Finally, while our simulation did not include the effects of noise due to finite photon statistics, we note the following:

Therefore, we expect that our MOOR approach should work within the limits of the resolution that is achievable given a particular feature contrast and photon fluence.

Algorithm 1 was implemented in C and parallelized across $m$ diffraction patterns in Eq. (11) using the Message Passing Interface (MPI). In particular, the reconstructed object defined by 17 million complex variables shown in Figs. 4 and 5 was obtained by running the MOOR algorithm on 2880 cores for a total of 350,000 core hours. We set the initial guess of the variables ${\mathbf{y}}^0=0$ and ${\mathbf{z}}_j^0={\mathcal{F}}^{-1}\{{\mathbf{d}}_j\}$; the magnitude of the probe variable ${\mathbf{x}}^0$ is initialized with $|{\mathbf{x}}^0|={\mathcal{F}}^{-1}\{\frac{1}{m}\sum _{j=1}^m{\mathbf{d}}_j\}$, and its phase is set to a Gaussian function. The step size parameters ${\alpha }_x$ and ${\alpha }_z$ were fixed to $4m$ and 2, respectively, throughout the algorithms. The step size ${\alpha }_y$ for the object variable was initialized to 100 and updated using a line search at each iteration. These parameters were chosen to ensure the convergence of Algorithm 1. The constraint variables ${\nu }_x$ and ${\nu }_y$ were both set to 1, although none of the constraints associated with these variables are active at the reconstructed image. We set the regularization parameter $\kappa =0.01$ and ran Algorithm 1 for a time budget that corresponded to approximately 600 iterations. This algorithm yielded successful reconstructions that agree well with the designed object, as shown in the images of Figs. 4 and 5 and also in the quantitative values of the reconstructed phase $\delta$ as shown in the histograms of Fig. 6.

 

Fig. 4. Comparisons of two cuts (zonal, top row; and meridonal, bottom row) of the reconstructed phase-shifting part $\delta$ of the x-ray refractive index $n=1-\delta -i\beta$ for different methods and data collections. These cuts show the reconstructed value of $\delta$ in the voxels at the selected planes. The true object is shown at left, followed by the MOOR reconstruction using 360° and then 180° rotation axes. Finally, a reconstruction from the standard PPA to ptychographic tomography is shown. As noted in the text, the pure projection approximation does not properly reproduce the object, and it also suffers from regular artifacts due to insufficient probe overlap for the reconstruction method used [59]. For this reason, we also show a column of “PPA+filtering” images with post-processing. These figures show that 360° MOOR gives a reconstructed image that represents the true object with a high degree of fidelity.

Download Full Size | PPT Slide | PDF

 

Fig. 5. Comparisons of 3D isosurface renderings of the true object and MOOR reconstructions using 360° and 180° object rotation. Rotation over a full 360° range gives improved results. The 3D object reconstructed using the PPA is not shown, as its errors (Fig. 4) do not make it possible to obtain a clean isosurface rendering.

Download Full Size | PPT Slide | PDF

 

Fig. 6. Histograms of the distribution of absolute error in the per-voxel values of the phase-shifting part $\delta$ of the x-ray refractive index $n=1-\delta -i\beta$ for both our MOOR reconstruction approach and for the filtered pure projection approximation (PPA, and PPA + filtering) reconstructions. At left are shown the histogram for object-containing voxels, while at right are shown the histograms for object-absent voxels. Also indicated are the values of $\delta$ for Si ($1.98\times {10}^{-5}$ at 5 keV) and ${\mathrm{TiO}}_2$ ($3.01\times {10}^{-5}$ at 5 keV). As can be seen for the object-containing voxels (left), the MOOR reconstruction approach gives results with absolute errors that are small compared to the actual values of $\delta$, while for PPA the absolute errors are almost as large as the expected values of $\delta$. For the object-absent voxels (right), again the errors in the MOOR reconstruction are very small compared to the expected values of $\delta$ in the object-present regions, while the PPA reconstructions have errors in the object-absent voxels that approach the true values expected in the object-present voxels.

Download Full Size | PPT Slide | PDF

For comparison with the (erroneous in this case) assumption that the object can be described by the pure projection approximation (PPA), we also reconstructed the simulated experimental data set using the single-slice ptychographic tomography (SSPT) approach [22,23] discussed in the Introduction. We first reconstructed a set of 2D complex-valued projections using 2000 iterations of a GPU implementation [60] of the extended ptychographic iterative engine (ePIE) algorithm [61]. To align all phase projections and account for linear and constant phase terms that typically plague independent ptychographic reconstructions of the same sample, we subtracted a best-fit plane from each reconstructed phase [23]. The aligned phase and magnitude projections were then tomographically reconstructed by using 200 iterations of the simultaneous iterative reconstruction tomography (SIRT) method as implemented in the TomoPy package [62]. We note that the employed ptychography and tomography algorithms assume a pure projection forward model, so their poor performance shown in Fig. 4 is expected in this case where the pure projection approximation does not apply. For example, propagation fringes at the boundaries of the sample are clearly visible, along with periodic artifacts caused by the relatively low number of projections (which can be suppressed by Fourier filtering [59]) and ePIE’s inability to retrieve the correct illumination function in the circumstances of our demonstration. In order to provide a fair comparison between the PPA and our MOOR approach, we post-processed the result of the PPA reconstruction (we emphasize that no post-processing steps were used for the MOOR algorithm results) to give a secondary result that we have labeled as “PPA + filtering.” First, we used median filtering with a window length of 10 in order to suppress the line artifacts that can be observed in the PPA column of Fig. 4; we then thresholded the values of $\delta$ of the filtered image with a threshold of ${10}^{-5}$.

3D ptychographic experiments typically only sample 180° around the object. This is a valid procedure under the pure projection approximation; in this setup, two samples taken from symmetrically opposed locations would, in theory, produce the exact complement of the projection image and corresponding diffraction pattern. However, situations beyond the pure projection approximation (thus requiring the multislice approach) break this symmetry, and our numerical results show that in this case a sampling around all 360° is necessary. Figure 4 shows the comparison of four ptychographic reconstructions with the true object: MOOR using 180° sampling, MOOR using 360° sampling, a reconstruction from a pure projection approximation (PPA), and a post-processed PPA reconstruction. In addition, Fig. 5 shows isosurface renderings of the true object and the 360° and 180° MOOR reconstructions. Together, these results show that the “hidden” side of the object is poorly resolved when using MOOR with only 180° rotation sampling.

5. DISCUSSION AND SUMMARY

SSPT is extremely successful [22,23,63,64] at obtaining high-quality 3D x-ray reconstructions of objects to which the pure projection approximation applies. For objects with larger extent, MSPT approach [30] has shown promise for treating the object as being represented by several planes along each viewing direction, which can be reconstructed and added to yield approximations of a pure projection. With MOOR, we take the approach of allowing for multiple slice propagation through thick objects along each viewing direction; an optimization approach is then used to recover a 3D object that is consistent with all diffraction measurements and which requires no phase unwrapping. Although we have not made an explicit comparison with the MSPT approach, the fact that one obtains an improved reconstruction with a 360° rotation relative to a 180° rotation means that a simple summation of separate reconstruction planes may not be sufficient to accurately reconstruct objects with increasing sizes beyond that to which the pure projection approximation applies.

In this first demonstration of the MOOR approach, a proximal alternating linearized minimization algorithm is used to obtain rapid convergence for the case of ptychography (where a small coherent probe is scanned across the entire specimen projection field of view at each rotation angle). This approach is consistent with the expectations of SSPT and MSPT, where the iterative rules alternate between updates of the object and the probe [61,65]. However, our approach employs the capabilities of high-performance computing to carry out a 3D calculation with isotropic voxel size and no need for phase unwrapping. It could also be used for other situations, such as probe overlap during rotation rather than at each viewing angle [66]. Our approach might also be useful in situations where the probe is larger than the object, and hence, the overlapping probe feature of ptychography is absent [67].

Although our approach is successful in recovering a computer-synthesized 3D object to which the pure projection approximation does not apply, there is room for future development. In this first demonstration, data parallelism was used so that each computing node owned and updated not only a subset of the variables ${\mathbf{z}}_j$ (and associated scans ${\mathbf{d}}_j$) but also a copy of the full 3D object, which, in turn, created a synchronization point in the parallel computation for merging those copies after every update step. This allowed us to show that the approach works but in an inefficient manner, where data communication (between all the nodes, with each one calculating a different probe position, and the 3D data) limited computational throughput. This limitation could be overcome by exploiting domain decomposition (i.e., parallelism in the object/reconstruction domain) for partitioning the 3D object based on the scanning pattern and the associated sparse intersections of localized illumination probes with the overall object. One could then have each node hold the probe function, a local subregion of the object, and the measured Fourier magnitudes for rapid computation with periodic synchronization happening only where the subregions overlap (we have used such an approach in standard ptychography reconstructions [60,68]). In addition, a stochastic asynchronous version of the proximal alternating linearized minimization algorithm [69] is likely to reduce the computational cost of the algorithm as well as the synchronization cost. We will also examine the capability of our algorithm with ptychography experiments for nanoscale objects including nanofabricated devices and subcellular structures of eukaryons. Technical challenges such as probe alignment (already addressed in 2D ptychography [70–72]) could be addressed when applying the algorithm to experimental results.