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 three-dimensional (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–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–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 $\mu$ 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 $\mu (x,y,z)$ with the vectors ${\widehat{\mathbf{e}}}_{\mathbf{x}}$ and ${\widehat{\mathbf{e}}}_{\mathbf{y}}$ spanning the plane normal to the axis of rotation and ${\widehat{\mathbf{e}}}_{\mathbf{z}}$ pointing parallel to the axis of rotation. We start with the radon transform of the local mass attenuation coefficient $\mu (x,y,z)$:

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\times 500\times 500$ voxels. A projection of the sphere phantom is given in Fig. 1(a), and an exemplary slice is given in Fig. 1(b).

 

Fig. 1. (a) Projection of the sphere phantom used for testing the registration. (b) shows an exemplary slice of the phantom, taken from the center of the phantom (slice no. 250). (c) shows a graph of the vertical displacement of each projection.

Download Full Size | PDF

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_{\theta }(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_{\theta }(z)$. Using the profile $M_0$ of the first projection as a reference, all other projection profiles $M_{\theta ,i}$ have been shifted to minimize the root mean square (RMS) deviation: $\mathrm{RMS}(M_0,M_{\theta ,i})$ has been calculated for pixel-wise shifts of each projection’s profile $M_{\theta ,i}$. Close to the minimum, the $\mathrm{RMS}(M_0,M_{\theta ,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 ${\sigma }_{\mathrm{RMS}}=0.0076\mathrm{pixels}$ with the largest deviation ${\mathrm{\Delta }}_{\max }=0.015\mathrm{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.

 

Fig. 2. $M_{\theta }$ profiles. (a) shows the profiles after applying the random movement and drift. (b) shows the same profiles after applying the correction.

Download Full Size | PDF

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 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.

 

Fig. 3. Visualization of the reconstructed backshifted data: (a) 3D ortho-slice, (b) xy-plane, (c) xz-plane, and (d) yz-plane. (e)–(h) show the difference between the original and the backshifted reconstruction: (e) 3D ortho-slice, (f) xy-plane, (g) xz-plane, and (h) yz-plane. The double bilinear interpolation creates rings in the xy-plane of the backshifted data.

Download Full Size | PDF

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:

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 ${\sigma }_{\mathrm{RMS}}(256\mathrm{ADU})=0.154$ is well below one pixel, and the maximum discrepancy between fit and simulated drift is ${\mathrm{\Delta }}_{\max }(256\mathrm{ADU})=0.469$. The difference decreases to ${\sigma }_{\mathrm{RMS}}(65536\mathrm{ADU})=0.005$ and ${\mathrm{\Delta }}_{\max }(65536\mathrm{ADU})=0.027$. Figure 4 shows the resulting differences in the fits from the simulated drift.

 

Fig. 4. Difference between fitted and simulated drift values for different ADU levels: both the RMS variance ${\sigma }_{\mathrm{RMS}}$ and the maximum deviation ${\mathrm{\Delta }}_{\max }$ asymptotically approach a value $ϵ$ close to zero.

Download Full Size | PDF

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 ${\sigma }_{\mathrm{RMS}}(256\mathrm{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 $\sigma =2\mathrm{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.

 

Fig. 5. Metrics for comparing the reconstruction from the original and backshifted data. The structural similarity index (MS-SSIM) of the original and backshifted data score very similar values at all noise levels. The integral of absolute value $Q_{\mathrm{IA}}$, integral of negativity $Q_{\mathrm{IN}}$, and histogram entropy $Q_H$ are measures for the quality of reconstruction and show an actual improvement of the reconstruction for the backshifted data which can be attributed to a statistical averaging of noise between adjacent slices.

Download Full Size | PDF

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_{\mathrm{IA}}$, integral of negativity $Q_{\mathrm{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.

 

Fig. 6. (a) Projection of the test sample. (b) Deviations between induced and calculated shifts. Possible movements of the beam in the flat-field and optics imperfections lead to a slightly larger deviation than for perfect simulated data.

Download Full Size | PDF

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.