Alignment methods for nanotomography with deep subpixel accuracy

As the resolution of X-ray tomography improves, the limited long-term stability and accuracy of nanoimaging tools does not allow computing artifact-free three-dimensional (3D) reconstructions without an additional step of numerical alignment of the measured projections. However, the common iterative alignment methods are significantly more computationally demanding than a simple tomographic reconstruction of the acquired volume. Here, we address this issue and present an alignment toolkit, which exploits methods with deep-subpixel accuracy combined with a multi-resolution scheme. This leads to robust and accurate alignment with significantly reduced computational and memory requirements. The performance of the presented methods is demonstrated on simulated and measured datasets for tomography and also laminography acquisition geometries. A GPU accelerated implementation of our alignment framework is publicly available. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Computed tomography (CT) and laminography (CL) methods enable 3D imaging and thus make it possible to acquire statistically relevant information about a sample's internal structure without its destruction. This allows for unique imaging applications broadly applied in material research, medicine, and industry, especially for in situ and ex situ studies. In particular, in combination with high-resolution X-ray microscopy, 3D nano-imaging is important for many applications in biology and materials science [1]. Two-dimensional sub-50 nm spatial resolution can be commonly reached in X-ray microscopy [2][3][4][5][6][7] for sample thickness of up to tens of microns. However, this high 2D resolution can be translated to comparable 3D resolution only if a sufficient number of measured angular 2D projections of the sample provide mutually consistent information. If the measured 2D projections are misaligned, the reconstruction quality deteriorates and the achievable spatial resolution is limited [8,9]. Since reaching the level of accuracy and mechanical stability required for sub-30 nm 3D resolution is currently not achievable and would drastically increase the complexity of the measurement device, alignment of the measured 2D projections during post-processing is an important step.
In the past decades, various methods for alignment of tomography projections were proposed, such as marker-based methods [10][11][12][13], cross-correlation [14][15][16][17], bootstrap methods [18][19][20][21][22][23][24][25][26], common-line [9,[27][28][29], and feature matching [30,31]. However, each of the methods has some drawbacks that hinder their application for high-accuracy alignment of large X-ray datasets. The marker-based methods, often used in electron tomography, require significant effort during the sample preparation and the subsequent marker detection. Furthermore, while for thin samples prepared for electron microscopy a marker is clearly visible in the projection, this is often not the case for majority of X-ray tomography samples. The cross-correlation-based methods suffer from accumulation of subpixel errors resulting in poor sensitivity to misalignment between angularly distanced projections. The common-line methods such as the center-of-mass method [28,29] for estimation of the horizontal shift and the vertical-mass-fluctuation (VMF) method for estimation [28,29] of the vertical shift are limited to the parallel-beam tomography with conventional geometry and cannot be used for more general 3D imaging geometries such as cone-beam tomography, laminography or interior tomography [20,32,33]. Finally, the projection matching alignment (PMA) methods [18][19][20][21][22][23][24][25][26] provide a very general way of projection alignment. However, PMA applications are often limited by their high computational demands and their convergence may depend on quality of the initial guess. Here, we describe a framework for alignment of 3D imaging datasets commonly produced by nanotomography, nanolaminography, and interior nanotomography. Our approach is also well suited for other high-resolution 3D imaging methods such as vector and tensor tomography [34][35][36][37] and non-rigid tomography [38]. The goal of our toolkit is to provide means for fast, robust and fully automated alignment of large datasets, which are collected in modern synchrotron 3D nano-imaging experiments.

Projection alignment strategy
With respect to the broad range of challenges met for high-resolution 3D imaging, the ideal alignment method needs to exhibit high robustness, low computational cost and high alignment accuracy. In order to meet these requirements, we have based our alignment strategy on algorithms with deep subpixel accuracy. The goal is to develop an approach that is able to align projections downsampled to a significantly reduced spatial resolution and simultaneously reach sufficient subpixel accuracy that allows for direct alignment of the full resolution dataset. This approach differs from the conventional pyramidal method [18,19], where the low-resolution alignment is used only as an initial guess in order to improve robustness, but the full resolution alignment step needs to be still performed.
Solving the alignment problem in a lower resolution leads to a significant reduction of the computational time. For example, the algorithmic complexity of alignment of M projections by the cross-correlation-based methods is proportional to O(MN 2 log N) and O(MN 3 ) for the projection matching methods, where N denotes the linear size of the reconstructed volume. Additionally, solving the alignment task in a lower resolution reduces memory requirements, improves angular sampling of the tomographic dataset according to the Crowther criterion [39], and generally increases the average SNR in the downsampled projections. Therefore, as we will demonstrate, our approach provides excellent results even for undersampled datasets with poor signal-to-noise ratio (SNR).
In order to provide maximal alignment robustness under various experimental conditions, our alignment strategy is based on the combined strength of three alignment methods as shown in Fig. 1. After loading the projections, an initial alignment is performed using the cross-correlation alignment (XCA) method. The next alignment step is provided by the vertical-mass-fluctuation (VMF) common-line method. If assumptions of this method are satisfied, it provides a very fast and accurate vertical alignment independent of any horizontal misalignment. If the assumptions are not satisfied, e.g. in laminography, or interior tomography configurations, this prealignment step is skipped. The final alignment is then provided by our implementation of the multi-resolution projection-matching alignment (MR-PMA) method.

Common contrast mechanisms in nanotomography
The combination of good penetration and high achievable resolution makes the hard X-ray wavelength range a common choice for nanotomography imaging. X-ray nano-imaging commonly relies on either of two contrast mechanisms: absorption or phase-shift. Due to broad variation of 3D tomography tasks, our alignment methods were designed to be used with either of these imaging modalities.
A consequence of the good penetration length of the hard X-ray photons is weak absorption contrast. The complex-valued transmission T(x, y) can be approximated for thin samples as where λ is the imaging wavelength, n(x, y, z) denotes the local refractive index n = 1 − δ + i β, and L Z is the sample thickness. Since the material constant ratio δ/β is much larger than one in the hard X-ray regime, phase imaging usually provides a significantly stronger signal, which makes phase usually the better contrast modality for imaging and projection alignment. However, phase-tomography has several major differences compared to the conventional absorption-based approach. Firstly, whereas the attenuation P A (x, y) can be easily calculated from the measured complex-valued projection as the phase is not uniquely determined by the measured complex-valued transmission T(x, y). If the sample is sufficiently weakly phase-shifting, so that the maximal relative phase-shift φ(x, y) is smaller than 2π, φ(x, y) can be recovered directly from the complex-valued projection as However, if the phase-shift exceeds 2π radians, phase unwrapping has to be performed [40] and a unique solution exists only for well-behaved projections, i.e. when the signal to noise ratio is high and the phase evolution is smooth without phase residua or vortices [41]. Additionally, each phase projection φ(x, y) can contain an arbitrary phase offset φ 0 and if the projections are reconstructed by ptychography [42], they also include an arbitrary linear phase ramp [8,9]. As we discuss in Section 8, these degrees of freedom become constrained once the aligned projections are adequately assembled in 3D space and thus they do not affect the final imaging quality. However, this limited reliability motivated us to develop alignment methods that are insensitive to these degrees of freedom and other potential systematic errors in the lowest spatial frequencies.

Initialization by cross-correlation alignment (XCA) method
Alignment of the adjacent projections by the cross-correlation alignment (XCA) method [14][15][16] is commonly used in cryo-electron tomography and a similar idea was tested for alignment of X-ray datasets as well [17,24]. Assuming that the average translation between two ideally aligned projections acquired at adjacent angles should be zero, an estimation of the relative shift between the adjacent projections can be used for alignment.
The XCA method can be used with projections that are real-or complex-valued, for generality we denote here the projections that we seek to align T Θ (x, y). Instead of simply estimating the mutual displacement between the raw projections, our implementation estimates the shift between the local variance fields V Θ (x, y) of the projections defined as Since XCA serves only as an initial guess, we use the computationally cheaper central differences for calculation of the gradients of the measured field T Θ instead of an FFT-based approach.
Since T Θ denotes a generally complex-valued field, our definition of V Θ does not require any phase unwrapping. This definition has also other advantages, for example, it makes the XCA insensitive to phase-residues, additional linear phase ramp, and potential errors in the lowest spatial frequencies of T Θ . Although the XCA method can significantly reduce the positioning error between projections acquired at adjacent angles, it is insensitive to a constant offset in the center of rotation. Assuming that the rotation axis is perpendicular to the illumination beam propagation direction and that the full angular range of 180 degrees is acquired, XCA between projections acquired at 0 and 180 degrees can be used to estimate and correct for such offset of the center of rotation [24].
In order to reduce the computational time and increase robustness against noise, the variance fields V Θ are downsampled to a 4× lower resolution grid using the Gaussian pyramid approach [43]. The total displacement for each projection can be calculated as a cumulative sum of the registered displacements between projections measured at adjacent angles. A drawback of the XCA method is accumulation of registration errors, in particular for datasets with large number of projections.

Common-line method -vertical mass fluctuation (VMF)
Because the XCA prealignment does not use any knowledge about the imaging geometry, the relative shift between angularly distanced projections can be still far from optimal. Therefore, robustness of the entire alignment process can be further improved if an additional prealignment method is used. Assuming the sample rotation around the vertical axis, parallel beam geometry, and the field of view capturing the entire sample along the horizontal axis in each projection, the vertical alignment can be estimated using a method called vertical mass fluctuation (VMF) [9]. VMF is a powerful and computationally cheap approach based on alignment of an integrated conserved quantity that we denote here m Θ (y) The vertical mass distribution m Θ (y) is independent of the rotation angle. In conventional transmission X-ray microscopy, the quantity p Θ (x, y) can be either the attenuation P A (x, y) or the unwrapped phase φ(x, y). A major advantage of the vertical-mass-distribution-based alignment is its independence of the horizontal misalignment. Due to the linearity of Eq. (5), a global offset or a linear slope in the preserved quantity p Θ(x, y) will result in a linear offset in the vertical mass distribution m Θ(y) . The same argument can be also applied for propagation of potential errors in the lowest spatial frequencies of p Θ(x, y) . In order to make our analysis insensitive to these errors, the high-pass filtered vertical mass distribution m Θ(y) by a linear operator is used for alignment, we denote the high-pass filter linear operator F hp . The vertical alignment is then performed as follows: 1. p Θ(x, y) is calculated for each projection, either from the unwrapped phase φ(x, y) or attenuation P A(x, y) .
2. p Θ(x, y) is integrated along the horizontal axis and high-pass filtered to calculate the vertical mass distribution. Note that if the phase is used and the selected high-pass filter is a gradient filter, there is no need for phase unwrapping since the vertical phase-gradient ∂ y φ can be calculated analytically from the complex-valued projections T Θ using where Im{} denotes the imaginary part and T Θ denotes the complex conjugate, u y is the vertical reciprocal-space coordinate and F {} is the Fourier transform. However, this strong high-pass filtering increases sensitivity of the VMF method to noise in the data.

A relative displacement between the one-dimensional VMF vectors is estimated
In order to guarantee a robust search in step 3) and avoid stagnation in local minima, the alignment is performed by subpixel accuracy cross-correlation [44] of the filtered vertical mass fluctuation.
Since an optimal choice of the cross-correlation reference is important, we have selected the most representative F hp {m Θ } vector as the centroid position provided by K-means clustering [45].
If the assumptions of the VMF method are fulfilled, the vertical alignment provided by this step can be considered to be final. Otherwise, the vertical alignment has to be refined in the final alignment step.

Projection-matching alignment (PMA) method
At this point, our dataset is already prealigned by the cross-correlation method and, if the measured dataset satisfies assumptions of the VMF method, the vertical alignment can be considered to be final. However, in order to keep our alignment strategy general, we will assume that both the vertical and the horizontal alignment have to be refined. Our final alignment step is performed by a multi-resolution projection-matching alignment (MR-PMA) method.
The PMA methods [18][19][20][21][22][23][24][25][26], also called "iterative reprojection" or "bootstrap methods", are based on the idea that even misaligned projections can be used for a 3D reconstruction, although the quality can be suboptimal. Reprojections of this initial reconstruction will be closer to the optimally aligned projections since only the consistent features are preserved in these synthetic reprojections. Therefore, shifting of the measured projections to match these reprojections will lead to improved reconstruction quality in the next iteration. Iterative application of this principle leads progressive improvement of the reconstruction quality until convergence is reached.
An application of the multi-resolution approach for projection alignment has been already demonstrated in [18,19]. Here, based on the same idea of multi-resolution processing, we demonstrate that with an optimized implementation, this method can reach alignment with so high accuracy in the coarse resolution levels that alignment at the finest resolution level does not further improve the reconstruction quality.
Solving the PMA task in lower resolution not only reduces the computational cost, but also the multi-resolution approach improves the robustness of the optimization. Essentially, assuming that the projections contain some features at each spatial frequency, the coarser resolution makes the optimization less sensitive to local properties of the measured projections and PMA becomes more robust against stagnation. The ability of MR-PMA to reach deep subpixel accuracy is based on the redundant information between the lowest spatial frequencies of the projections in the 3D Fourier space as indicated in Fig. 2. Since mutual displacement of the measured projections results in a mismatch in these spatial frequencies, this discrepancy can be sensed and used for correction by the PMA method. PMA methods generally minimize the difference between the linearized measurements p Θ and their reprojectionsp Θ defined asp where the variable p denotes a stack of the measurements The linear operator A Θ denotes a forward tomographic projection operator at angle Θ and the inversion operator A −1 denotes generally any common tomography reconstruction method such as filtered backprojection (FBP) or simultaneous algebraic reconstruction technique (SART). In our implementation, we use FBP mainly because of its low computational cost and versatility in non-standard imaging geometries [46]. Since the PMA computation is dominated by the reprojection step, use of Fourier space reconstruction methods such as Gridrec [47] could be beneficial only if the forward tomographic projection A Θ could be implemented efficiently in Fourier space as well.
The displacement between the measurements p Θ and the reprojectionsp Θ is commonly estimated by cross-correlation methods, which can be implemented with deep subpixel accuracy [44,48,49]. However, we have observed that the cross-correlation-based methods may become unstable if the data contains systematic errors. Additionally, the computational complexity of these methods grows significantly with the aimed level of subpixel accuracy [49]. Therefore, we have adopted an optimization approach inspired by the optical flow methods [50], which we found to be more robust and can iteratively reach deep subpixel accuracy. The generalized PMA optimization task can be expressed by the following cost function where W Θ denotes reliability weights, which are equal to zero in missing or unreliable regions. Given the same arguments as in the previous sections, both the measured data p Θ(x, y) and the reprojectionsp Θ(x, y) are high-pass filtered by a linear operator F hp in order to reduce the influence of potential errors in the lowest spatial frequencies. These errors can originate either from imperfections in the projection p Θ or inconsistencies in the reprojected model in Eq. (7), e.g. due to the missing cone in laminography geometry. By applying enough downsampling, we can assume without loss of generality that at the coarsest resolution level the displacement errors are much smaller than 1 pixel. The optimal displacements ∆x (i) Θ , ∆y (i) Θ , which minimize the cost function Eq. (9) in the i-th iteration, can be then estimated analytically as where ∂ xp Θ denote horizontal and vertical gradient of the reprojectionp (i) Θ . The gradients are calculated in Fourier space as where u x and u y denotes the reciprocal-space coordinates. The estimated shifts ∆x (i) Θ , ∆y (i) Θ are used to correct the displacements errors directly in the measured projections The cost function, Eq. (9), needs to be optimized iteratively, because the corrected displacement in p (i+1) Θ will lead to improved quality of the tomogram and therefore of the reprojectionsp (i+1) Θ . The optimization is finished when the update step drops below a certain threshold γ/D where D is the downsampling level and γ is a small constant selected to be 10 −2 for our application. The gradients ∂ xpΘ , ∂ ypΘ are commonly approximated by their finite difference, and displacement of the projections is typically approximated using by bilinear interpolation, which is applied by a change of the projection tomographic geometry [18,19,21,24,25]. However, such strategy is not able to provide the deep subpixel accuracy required to avoid iterative alignment without downsampling. Additionally, it creates numerical noise that hinders convergence [24]. In contrast to other works, we use Fourier transforms to calculate the gradients of the reprojection and to apply the displacements to the measured projections in Eqs. (11) and (12). The discrete Fourier transform is a formally invertible numerical transformation with high numerical accuracy and stability and this is of great importance to achieve deep subpixel alignment accuracy.

Interpolation by offset discrete Fourier transform (DFT)
The multi-resolution optimization approach requires a way to upsample and downsample the projections to different interpolation grids. The selection of the optimal interpolation methods is critical to preserve the alignment quality between various interpolation grids and to preserve the deep subpixel accuracy acquired in the coarse resolution levels. In our implementation, we resample of the original signal x(n) of length N to a new signalx(n) of length M using the discrete Fourier interpolation with an offset of half of the sampling step wherex(n) denotes the resampled signal of length M, and the discrete inverse Fourier transform of x(n) is computed also with an offset, namely, X(m) is given by The main motivation for this offset is to preserve the pixels centers after the interpolation, i.e. there is no shift of the bulk properties such as center of mass of the projections. The computational cost of the offset Fourier interpolation is comparable to the conventional discrete Fourier interpolation, because we implement this offset interpolation using two fast Fourier transformations with a linear phase ramp multiplication in the Fourier domain.

Ambiguities of 3D phase-imaging using ptychography
The global phase offset and linear phase ramp are inherent degrees of freedom of ptychography [9]. Although the discussed alignment methods exhibit high robustness against various lowspatial-frequency errors, the quantitative contrast of the final reconstruction can be guaranteed only if these ambiguities are adequately constrained. The conventional approach is to select regions in each frame that do not contain any sample contribution [9,24] and subtract the phase ramp in each projection separately. However, this approach may result in inaccuracies if the sample-free regions are not well detected, e.g. when contrast of the sample is low, or if the sample projections do not contain enough empty regions around the sample. Our approach here is based on the fact that not all possible phase ramps are consistent with the 3D reconstruction, in other words tomographic consistency can be used to deal with these ambiguities. Using the projection matching idea in Eq. (7), the optimal phase ramp φ Θ (x, y) = β Θ, x x + β Θ, y y + β Θ,0 at angle Θ can be found by minimization of the following least-squares optimization task min β Θ, x ,β Θ, y ,β Θ,0 where p Θ denotes the unwrapped phase and W Θ are reliability weights. Without additional boundary conditions, the conventional parallel beam tomography allows for two degrees of freedom per tomogram: a shared phase ramp along the rotation axis direction and a shared phase offset. The laminography geometry reduces the number of the degrees of freedom to a single one: the shared phase offset. Although these maximally two scalar variables per tomogram still need to be calibrated by some prior knowledge, robustness of the phase ramp removal process is significantly improved because this knowledge does not need to be available for each projection. An identical approach can be used for calibration of the absorption tomogram, for which ptychography has an ambiguity in the overall amplitude of the reconstructed complex-valued transmission.

Results for simulated tomography datasets
Our alignment framework is implemented in Matlab using the ASTRA framework [21,51] with additional modifications [52]. The reconstructions were performed with CPU Intel Xeon E5-2690 and GPU Nvidia P100 and 512 GB RAM. Performance of the proposed alignment methods was first verified using a synthetic dataset. In order evaluate our methods for a common class of samples in nanotomography, we have created a phantom of 800×800×800 pixels of a porous structure with two material phases, shown in Fig. 3(a), the structure was inspired in the nanostructure of catalyst materials. This phantom was used to produce 1260 projections, which provide full angular sampling according to the Crowther criterion [39]. The projections were displaced by a combination of a normally distributed random offset with amplitude of 20 pixels and an angularly dependent sinusoidal offset with period of 90 degrees and amplitude of 20 pixels, as shown in Fig. 3(b). Three test datasets were created using these projections. Dataset 1 was fully angularly sampled with noiseless projections, dataset 2 was 8× angularly undersampled, i.e. it contained only 161 noiseless projections with equal angular sampling, and finally dataset 3 adds 10% of Gaussian noise to the dataset 2.
The XCA method is able to significantly reduce the initial positioning errors. As shown in Fig. 3(b), the position variance between adjacent projections is significantly reduced. However, the overall alignment of the projections, in particular between those with large angular distance, is suboptimal, as expected. The remaining positioning errors for the three methods are detailed in Table 1. The main advantage of the XCA method as a prealignment step is its low computational cost. In this case, the variation calculation, shown in Eq. (4), and downsampling operation required 7 s and the cross-correlation between all 1260 projections required additional 2 s. The alignment quality, summarized in Table 1, is comparable for the first two tested datasets, but the root mean square (RMS) error was increased for Dataset 3 due to noise-sensitivity of our variance projection preprocessing in Eq. (4).  , the result at previous downsampling level is provided as the initial guess for the next higher resolution level. At downsampling of 4×, 2× and 1× only one iteration is needed due to the deep subpixel accuracy at higher downsampling levels.
The following optional prealignment method is VMF. In order to provide a more fair comparison, we have not used the cross-correlation alignment from the previous step as the initial guess, and we have applied VMF method directly on the initially displaced projections, as it is used in [9,22,53]. The alignment quality in Table 1 shows that our implementation of the VMF method can reach deep subpixel accuracy, resulting in a displacement accuracy better than 0.2 pixel RMS for all the test datasets. The computational time is comparable with the XCA method, the 2D FFT-based unwrapping [40] of all 1260 projections required 20 s and the cross-correlation alignment of the vertical mass distributions needed only 0.7 s.
The last tested method was our implementation of the MR-PMA method. Similarly, to the previous test, we have not used the advantage of an initial guess provided neither by the XCA nor by the VMF methods. Instead, we perform the reconstruction on multiple resolution levels. Once the MR-PMA method converges according to the criterion in Eq. (13) on the coarser resolution level, the result serves as the initial guess for the progressively finer resolution levels. Our MR-PMA method reached sub-0.2 pixel RMS accuracy for all tested configurations in the final resolution level as shown in Table 1. However, the remarkable feature of the deep-subpixel alignment is the option to finish the alignment already in one of the coarser resolution levels without need to iteratively compute the tomogram and projections at full resolution. This is demonstrated in Fig. 4(a), where we show the evolution of the RMS displacement error in different undersampling levels. The results are very consistent between the tested datasets and only a weak dependence on the angular sampling and the noise in the data can be observed. Note that the errors in Fig. 4(a) refer to the corresponding residual errors at full resolution. Alignment at the 32× downsampled resolution level was already sufficient to reach the RMS projection displacement error below 0.2 of full resolution pixels, i.e. 0.00625 downsampled pixels, for all the tested synthetic datasets as shown in Fig. 4(a). The level of detail in the sample obtained by the FBP method for different downsampling ratios can be seen in Fig. 5.
The per-iteration computational time of our PMA implementation is summarized in Fig. 4(b). In the case of highly downsampled resolution, the computational cost is limited by a constant overhead. For downsampling factor below 8, the per-iteration computational cost grows roughly 8× with each resolution level. On the other hand, the accurate initial guess from previous steps allows reaching the convergence threshold provided by Eq. (13) in a single iteration as show in Fig. 4(b).

Results for simulated laminography datasets
Our implementation of the MR-PMA method can be applied for different imaging geometries, e.g. laminography. Laminography is a generalization of conventional tomography with the only difference that the sample rotation axis is not perpendicular to the illumination beam direction. This geometry can be used to image planar 3D objects that would be difficult or impossible to image by conventional tomography in a nondestructive way [54]. However, the laminography geometry brings additional challenges to the projection alignment. Most importantly, the vertical and horizontal alignment cannot be solved separately. Additionally, the tilted axis of rotation implies that a meaningful reconstruction can only be obtained by calculating the full volume or a significant part of it, which leads to a significant increase in computational demands. Furthermore, there are no conserved quantities, such as VMF or center of mass that can be leveraged for alignment.
We have used the phantom from the previous section to generate 1260 projections in parallel beam laminography geometry with a 60 degrees angle between the rotation axis and beam propagation direction. These projections were used to create three datasets in the same way as the tomography datasets. The displacement error after alignment by our PMA method for different downsampling levels is shown in Fig. 6(a). Laminography alignment is also a more complex task due to the missing cone information and typically unlimited support of the imaged sample, which may result into inconsistencies between the measured projections and their reprojections. Despite these limitations, the alignment threshold of 0.2 pixel RMS was reached already for 8× downsampled resolution for all three test datasets.

Alignment of experimental datasets
In order to demonstrate the full potential of the developed methods, we present here the alignment of a ptychographic X-ray laminography (PyXL) dataset. Details of the measurement, acquired at the cSAXS beamline, Paul Scherrer Institut (PSI), Switzerland, can be found in [55]. The measurements were acquired in the laminography nano-imaging instrument (LamNI) [56] equipped with optimized illumination optics [57]. The sample was an integrated circuit (IC) printed using 16 nm node technology. The substrate was mechanically polished to a thickness of about 20 µm, however all the scattering features are confined within a 4 µm thickness. We have collected 2704 ptychographic projections with a circular field of view of 40 µm diameter in the plane of the sample. The dimensions of the reconstructed ptychographic projections were 2549×4155 pixels with a pixel size of 13.0 nm. The corresponding reconstruction volume contains 3800×3800×600 elements.
This dataset presents several challenges, in particular the size of the reconstructed volume. As discussed beforehand, laminography alignment does not allow the vertical and horizontal alignment to be solved separately and, depending on the sample thickness, the entire imaged volume or at least a significant part of it needs to be calculated in each iteration. Since the data volume of the reconstructed ptychographic complex-valued projections is 230 GB in singleprecision (32-bit), the alignment task becomes extremely demanding, both in terms of memory and computation. Another complication for this IC sample is that it contains quasi-periodic features, which makes alignment methods more prone to convergence into local minima. The combination of these challenges makes this dataset ideal for demonstration of the unique capabilities of our alignment framework.
As in the previous sections, the alignment was solved in several resolution levels. The coarsest resolution level was 64× downsampled and the finest was 2× downsampled. Due to memory limitations, we were not able to perform the alignment step at the finest resolution, therefore the ability to reach alignment with subpixel accuracy was critical for this dataset. Since the ground truth alignment is unknown, in Fig. 7(a) we report the alignment progress in each resolution level using the displacement RMS value.
The displacements per iteration reported in Fig. 7(a) indicate that the 4× downsampled resolution level already provides sufficient alignment precision, because the update RMS magnitude in the following step is below 0.1 pixel. The required computational time in the 2× downsampled dataset was 450 s per iteration as shown in Fig. 7(b). The total displacement for each projection angle is shown in Fig. 7(c). The projection displacements are caused among others by tracking the selected field of view in our sample by coarse motors upon rotation, imperfections of the rotation stage, and long-term thermal drifts. The maximal displacement error was up to 4 µm, i.e. 10% of the targeted field of view.
The laminography reconstruction after alignment is shown in Fig. 8. The reconstruction was performed by a filtered back-projection method modified for laminography geometry [58]. The corresponding spatial resolution of 18.5 nm was estimated from the intersection of the Fourier shell correlation (FSC) with the 1/2-bit threshold [59] shown Fig. 8(c). The FSC calculation excluded the Fourier-domain missing cone. This spatial resolution is sufficient to recognize the

Conclusion
We have demonstrated here that use of alignment methods for 3D nano-imaging implemented with deep subpixel accuracy can be achieved by a computationally efficient approach for tomography and laminography datasets. The subpixel accuracy of our multi-resolution projection matching method (MR-PMA) makes it possible to avoid sample reconstruction and reprojection at the full spatial resolution during alignment and thus helps to significantly reduce the computational and memory requirements.
Additionally, as we have demonstrated, our alignment methods perform well for noisy and undersampled datasets. The difficulties related to undersampling are alleviated since sufficient downsampling of the projections makes it always possible to satisfy the Crowther angular sampling criterion. Therefore, we have not observed any improvement of the reconstruction performance when iterative reconstruction methods such as SART were used for reconstruction instead of the conventional FBP method.
The combination of all three presented alignment methods allows us to perform the tomography processing in most cases in an unsupervised automatic manner. This helps to avoid the need of user interaction during the alignment process which is a significant bottleneck in nanotomography and nanolaminography imaging.
Our GPU-accelerated MATLAB-based alignment and reconstruction framework is publicly available in 10.5281/zenodo.3360818 [52] including a simple example of an artificial dataset. A high-resolution ptychotomography measurement of nanoporous glass [53] is provided in 10.5281/zenodo.3539513 [60]. The aligned full-resolution complex-valued projections from our experimental laminography dataset and the corresponding reconstruction are provided in 10.5281/zenodo.2657340 [61].

Disclosures
The authors declare no conflicts of interest.