Analytical registration of vertical image drifts in parallel beam tomographic data

Reconstructing tomographic images of high resolution, as in x-ray microscopy or transmission electron microscopy, is often limited by the stability of the stages or sample drifts, which requires an image alignment prior to reconstruction. Feature-based image registration is routinely used to align images, but this technique relies on strong features in the sample or the application of gold tracer particles, for example. In this Letter, we present an analytic approach for achieving the vertical registration based on the inherent properties of the data acquired for tomographic reconstruction. It is computationally cheap to implement and can be easily integrated into existing reconstruction pipelines.

Reconstructing tomographic images of high resolution, as in x-ray microscopy or transmission electron microscopy, is often limited by the stability of the stages or sample drifts, which requires an image alignment prior to reconstruction. Feature-based image registration is routinely used to align images, but this technique relies on strong features in the sample or the application of gold tracer particles, for example. In this Letter, we present an analytic approach for achieving the vertical registration based on the inherent properties of the data acquired for tomographic reconstruction. It is computationally cheap to implement and can be easily integrated into existing reconstruction pipelines. Image registration is often required for high-resolution tomographic imaging, as the resolution approaches or surpasses the mechanical or thermal stability of the setup, for example, in x-ray ptychography, full-field transmission microscopy (TXM), and transmission electron microscopy (TEM).
While image registration algorithms are readily available with sub-pixel precision [1,2], in the case of tomographic data, the image changes for every rotation step. Since the sample is rotated between any two images, a direct image registration is not possible. Each tomographic projection can be registered to the preceding and/or following projection, but this method can only compensate for small-scale random movements. Drifts cannot be easily corrected in this manner, and an image drift can accumulate over a stack of several hundreds or thousands of images, even if each individual registration has sub-pixel accuracy. A better route is using information from the threedimensional (3D) structure of the sample itself.
The routine technique of image registration in TEM and soft x-ray microscopy employs known markers (fiducials), most often spherical gold nanoparticles, which are suspended in a solution added to the sample prior to cryo-fixation. The fiducials then disperse in the sample which gives a relative uniform distribution of these particles. These markers give a strong signal and are easy to track in the images to register the latter [3][4][5].
If the use of markers is not possible or desired, the image data are analyzed directly, either by using image corners (image regions with a large local contrast gradient) and tracking their positions [6][7][8] or by using the image cross-correlation [9,10], of either the whole image or of small subsets.
While these current techniques work well, they only use very distinct features of the projections. The information about the sample is contained in the full images, even if the signal strength (e.g., the local absorption) is rather weak compared to markers. However, as witnessed by the tomographic reconstruction, it carries all the information required. In this Letter, we present an analytical approach to the vertical image registration in parallel beam tomography based on the tomography geometry itself. The approach is tested against different levels of noise, and the quality of the reconstruction is compared.
Many full-field imaging techniques such as TXM can make use of this novel algorithm to correct for vertical movements. However, any correction of the horizontal jitter would need to be performed by other means. Furthermore, it could also be used as a computationally cheap check for confirming sample stability in applications such as microtomography if inexplicable reconstruction artifacts appear and sample movement is suspected.
Our approach requires that the sample be completely included in the image and does not move out of the field-of-view. This limitation prohibits the applicability to techniques such as TEM, where an extended thin sample is measured and the whole sample volume is not captured in each image. In addition, the application of our approach requires that images are flat-field corrected and that the flat-field correction does not induce any systematic artifacts in the image, for example, from beam drifts.
We follow here the nomenclature of μ for x-ray absorption data, but any other normalized signal such as the phase shift can be used too. The sample is described by its mass attenuation coefficient μx; y; z with the vectorsê x andê y spanning the plane normal to the axis of rotation andê z pointing parallel to the axis of rotation. We start with the radon transform of the local mass attenuation coefficient μx; y; z: The coordinate transforms t x cos θ y sin θ and u − sin θ y cos θ have been used in Eq. (1). Expressed in variables x; y; z, Eq. (1) is We know that μx; y; z ≠ 0 only for x ∈ x 0 ; x 1 , y ∈ y 0 ; y 1 with the limits x 0 ; x 1 ; y 0 ; y 1 defined by the sample outer dimensions. Therefore, the integration boundaries of x and y in Eq. (2) can be limited to the intervals given above. Furthermore, consider an integration of p θ t; z over t: Substituting for p θ t; z with Eq. (2), we get Here, the integration order has been changed, and the integral over the Dirac delta function has been executed. The resulting integral of μx; y; z over x and y is a constant for any given z. The value M θ; z M θ z is independent of θ, i.e., a constant for every angular projection. Therefore, the function M θ z can be computed for every angular projection and has a distinct profile with the height. If the detector image shifts over the course of a measurement, the profile M θ z will shift, too. For discrete detectors used in tomography setups, the integration simply changes to a sum over each detector row. Using discrete data does not impose any limitations, as the measured intensity in each discrete pixel is the integrated intensity over the pixel area. The vertical data are sampled only at discrete positions, corresponding to the pixels.
For registration of an image sequence, the profile can be tracked and registered to a reference, for example, the first image of a series.
A simulation has been performed with a phantom of spheres generated in XDesign [11]. The sample has been generated from 2500 spheres with radii of 2-30 pixels in a cylinder in a volume of 500 × 500 × 500 voxels. A projection of the sphere phantom is given in Fig. 1(a), and an exemplary slice is given in Fig. 1(b).
For testing, each image has been shifted vertically by the combined value of a random movement and a global drift of the sample position with a maximum displacement of approximately 25 pixels. As the test structures are spheres of up to 30 pixel radius, this operation will render the data set unusable. Figure 1(c) shows the vertical displacement for each image. Sub-pixel displacement has been performed by bilinear interpolation of the original image. Reconstructions of the phantom were performed using GridRec [12] and the TomoPy package [13].
The profile M θ z has been calculated for the shifted data set. Small deviations in the local mass density for each slice lead to a variation of M θ z. Using the profile M 0 of the first projection as a reference, all other projection profiles M θ;i have been shifted to minimize the root mean square (RMS) deviation: RMSM 0 ; M θ;i has been calculated for pixel-wise shifts of each projection's profile M θ;i . Close to the minimum, the RMSM 0 ; M θ;i curve can be well approximated by two linear functions. The intersection of these two curves gives the sub-pixel value for the vertical shift. This simple method is quite robust and yields sufficiently good results. However, more sophisticated methods to calculate the difference between two one-dimensional curves could be used, if desired. The mean RMS deviation between the fit and the simulated drift was σ RMS 0.0076 pixels with the largest deviation Δ max 0.015 pixels. The profiles of the randomly shifted and the corrected data are shown in Fig. 2. Again, the data have been shifted to sub-pixel accuracy with bilinear interpolation.
The reconstruction of the backshifted data given in Fig. 3(a) shows that the drift can be corrected and an artifact-free reconstruction is possible. Note that the double bilinear interpolation of the data does introduce rings around some features, as shown in Figs. 3(e) and 3(f ): these occur when the  Letter diameter of the features in the xy-plane changes strongly for changes in z, i.e., at the north and south caps of the spheres.

Figures 3(g) and 3(h) show renderings in the xz∕yz-direction
which show that the rings become more pronounced at the caps and disappear at the equator of the spheres.
Noisy data have been simulated by applying a Poisson noise to otherwise ideal data. In accordance with a typical tomography measurement scheme, the projection has been calculated from sample transmission and flat field with Poisson noise: Here, Pλ denotes a random Poisson sample of the mean λ.
The theoretical sample transmission is given by T, and the statistics is determined by the number of analogue to digital units (ADUs), i.e., the detector counts. The sample transmission T is scaled to a minimum T min exp−2 ≈ 0.135 to use the optimum dynamic range of the detector [14]. ADU values from 2 8 256 up to 2 16 65536 have been tested, as these correspond to typical numbers of detector counts. The quality of the fit improves with reduced noise levels, but even at very low count rates, a reasonable result is achieved. At 256 ADUs, the variance σ RMS 256 ADU 0.154 is well below one pixel, and the maximum discrepancy between fit and simulated drift is Δ max 256 ADU 0.469. The difference decreases to σ RMS 65536 ADU 0.005 and Δ max 65536 ADU 0.027. Figure 4 shows the resulting differences in the fits from the simulated drift.
Reconstructions have been performed for data with different noise levels to test the robustness of the vertical registration. Even at 256 ADUs with an average discrepancy of σ RMS 256 ADU 0.154, no artifacts stemming from misalignment of the data are visible. However, the reconstruction is dominated by noise. The quality of the reconstructions has been assessed by applying various metrics. The similarity of the reconstructed images and the original phantom has been compared using the multi-scale structural similarity index (MS-SSIM) [15] at different scales (2,4,8,16, and 32 pixels). The similarity of the reconstructions and the original phantom is scale-dependent, but both the original reconstructions and those from the backshifted data are virtually identical; the largest discrepancy from the noise-less data and the data sampled with 256 ADU is 0.31%. Figure 5 shows the data for the MS-SSIM metric at σ 2 pixels. Slight variations occur only at low ADU levels, i.e., with very noisy data. Interestingly, the backshifted data achieve slightly higher results than the original data. This can be explained by the noise levels: the backshifted data were twice interpolated through the bilinear image shifts,   leading to an effective averaging over neighboring slices. Therefore, the metrics of the backshifted data correspond very well with the ones from the original data with an ADU level increased by a factor of two. However, these metrics only assess the in-plane quality of the reconstruction. The obvious problem of the interpolation is a reduced fidelity in the z-direction, as shown in Fig. 3(b), and it is not advantageous to perform unnecessary interpolation steps to achieve nominally enhanced metrics.
The quality of the reconstructions was estimated using metrics described by Donath et al. [16]. While these metrics were originally developed to test the position of the center of reconstruction, they are generally valid for comparing the quality of different reconstructions. The integral of absolute value Q IA , integral of negativity Q IN , and histogram entropy Q H were calculated for all reconstructions. Figure 5 shows the resulting metric values. Again, the backshifted data achieve slightly better results which are comparable with the original data of twice as many ADUs, but at the cost of reduced quality in the z-direction.
A further test of this scheme was performed by applying random shifts to microtomography data acquired at the Diamond Manchester Imaging Branchline I13-2. A series of five tomographies of a limestone sample with different statistics was acquired. A sample projection is given in Fig. 6(a). A slow drift and a random jitter have been applied to these microtomography projections, and the shift has been calculated using the scheme presented in this Letter. The results are shown in Fig. 6(b). For high statistics, the achieved fit quality is very similar to the simulated data, but the deviation is larger for lower count rates. In this case, lower count rates have been achieved by shorter exposure times, which increase the effect of fast beam movements for flat-field corrections. However, even for the lowest statistics, a fit accuracy of well below 1 pixel has been achieved.
An analytic scheme for tomographic data registration in the vertical direction was presented. While this scheme cannot be used in all geometries, it is valid for full-field parallel beam tomography as used in synchrotron sources. Implementation is computationally cheap and can be easily integrated into different reconstruction packages as it is a standalone step. For example, prealigning an image stack in height using the scheme presented in this Letter would give a better starting point for iteratively aligning projections by cycling between reconstructions and projecting the reconstructed data.
The required interpolation of image data poses a constraint on the visibility of the smallest features, but it has been shown that the general effect of finite statistics in the tomographic dataset is limiting the quality more than the image interpolation. In addition, interpolation of the data is intrinsic in correcting for any drifts or movements and not specific to this method. In any application, the image would only be shifted once and not twice, as was required in the case of simulated noise, reducing the effect of image interpolation even further.