Automatic marker-free estimation methods for the axis of rotation in sub-micron X-ray computed tomography
Abstract
Misalignment of the rotation axis causes severe artifacts in X-ray computed tomography. Calibration of this parameter is often insufficient for sub-micron resolution measurements and needs to be corrected during the post-processing. This correction can be accelerated by various automatic methods. These vary in mechanisms and performance, making them suitable for different use-cases. This work summarizes existing automatic methods for estimating the rotation axis in X-ray computed tomography, with a focus on sub-micron applications. Some of the methods are implemented and compared in the context of a laboratory sub-micron scanner to demonstrate practical considerations of this task.
Full text
Tomography of Materials and Structures 1 (2023) 100002 Contents lists available at ScienceDirect Tomography of Materials and Structures journal homepage: www.journals.elsevier.com/tomography-of-materials-and-structures Automatic marker-free estimation methods for the axis of rotation in sub-micron X-ray computed tomography Marek Zemek a,1 , Jakub Šalplachta a , Tomáš Zikmund a,⁎ , Kazuhiko Omote b , Yoshihiro Takeda b , Peter Oberta c,d , Jozef Kaiser a a Laboratory of X-ray micro and nano computed tomography, Central European Institute of Technology, Brno University of Technology, Purkyňova 656/123, Brno 612 00, Czechia b Rigaku Corporation, 3-9-12, Matsubara-cho, Akishima-shi, Tokyo 196-8666, Japan c Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, Prague 182 21, Czechia d Rigaku Innovative Technologies Europe s.r.o., Za Radnicí 868, Dolní Brežany 252 41, Czechia ARTICLE INFO Keywords: Computed tomography Rotation Axis Tuning-fork artifact Automatic ABSTRACT Misalignment of the rotation axis causes severe artifacts in X-ray computed tomography. Calibration of this parameter is often insufficient for sub-micron resolution measurements and needs to be corrected during the post-processing. This correction can be accelerated by various automatic methods. These vary in mechanisms and performance, making them suitable for different use-cases. This work summarizes existing automatic methods for estimating the rotation axis in X-ray computed tomography, with a focus on sub-micron applications. Some of the methods are implemented and compared in the context of a laboratory sub-micron scanner to demonstrate practical considerations of this task. 1. Introduction X-ray computed tomography (CT) is a tool for the non-destructive imaging of internal structures of samples [1–4,5]. A CT measurement consists of scanning a series of projections over a range of angles and processing the acquired data by tomographic reconstruction to yield cross-sectional images (tomograms) of scanned objects [6]. The raw data produced by a modern laboratory CT scan can be viewed either as individual, usually two-dimensional projections or as sinograms, which display a single row of projection data at all acquired angles (Fig. 1) [6]. Contemporary scanners are capable of reaching sub-micrometer (sub-micron CT) [7] or even higher (nanoCT) [8] resolutions. The geometric alignment of scanner components (Fig. 1) is a significant factor affecting the quality and uncertainty of CT measurements [9]. The most critical alignment parameter is the position of the axis of rotation (AoR) [10,11] and its projection onto the detector. Tomograms are severely degraded by characteristic tuning fork streak artifacts in half-scan data [6,10] or double edges in full-scan data [6] when the physical position of the AoR does not coincide with the AoR assumed during the reconstruction. These artifacts decrease in magnitude as the assumed AoR position approaches its true location (Fig. 2) but even relatively small errors can cause a sufficient disturbance to render the resulting tomograms unusable for any further analysis. Such artifacts should thus be reduced as much as possible to minimize their negative impact on the diagnostic potential of images. The AoR can generally be calibrated using specialized phantoms and procedures before performing a scan [11]. The phantoms are usually made of metal wires [11] or other high-contrast objects. The AoR is localized by evaluating the symmetry of the phantom as it rotates on the sample stage [11]. Correction of the AoR is then done by hardware or software tools like adjusting the rotational stage or shifting the acquired data, respectively. However, this calibration is often insufficient at the high resolutions of sub-micron CT because of the limited precision of the sample stage [14,15] and a reintroduction of the AoR misalignment while mounting the sample or due to thermal expansion [16,17]. The position of the AoR must therefore be localized and corrected after scanning by analyzing and processing the acquired projection data. The AoR can also be estimated post-scan in a manual or automatic fashion. The estimation process can be enhanced using fiducial markers placed on the sample [18]. Typical fiducial markers are small particles (comparable in size to a single pixel [18]) with a high contrast [19,20], but other objects like thin wires [18] may also be used. Features of the sample stage assembly within the field of view (FoV) can also serve as references for the AoR position in some cases [21]. The disadvantage of https://doi.org/10.1016/j.tmater.2022.100002 Received 6 August 2022; Received in revised form 9 December 2022; Accepted 18 December 2022 Available online 26 December 2022 2949-673X/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ]]]] ]]]]]] ⁎ Corresponding author. E-mail addresses: [email protected] (M. Zemek), [email protected] (T. Zikmund). 1 0000-0002-3236-4111
fiducial markers is that they can obscure sample structures and complicate both the sample preparation and the data processing [18,22,14,20] This means that it is preferable to estimate the AoR using unaltered projection data whenever possible. A manual AoR estimation consists of an operator interactively changing the assumed position of the AoR in an effort to minimize artifacts in the reconstructed tomograms (Fig. 2) [16,17]. This is a laborious process that is prone to errors, but a human operator can also leverage their experience and intuition to deal with complex cases [17]. The automatic estimation mimics this process by scoring quantifiable effects of AoR shifts in images using objective metrics. This can increase the data throughput, lower the operator’s workload, and provide objective and repeatable results. An automatic AoR estimation is often implemented as a stand-alone processing step or as a part of a more general scheme for correcting geometric misalignments. An estimate of the AoR position can also be obtained as a by-product of some algorithms for a per-projection motion correction [18,22]. This is because the AoR misalignment is in essence a special case of the general motion problem in which the sample and the stage are misaligned by a constant horizontal shift in all projections. This article focuses on methods aimed specifically or primarily at the automatic AoR estimation by processing the scanned datasets without the aid of dedicated phantoms or fiducial markers. The circular scan trajectory (the source and the detector midpoint both rotate in a single plane) is the most widely used trajectory in most laboratory CT applications, and it will be referred to in this work unless stated otherwise. Similar methods are also documented for related modalities such as laminography [23] or optical projection tomography (OPT) [24], but these are beyond the scope of this review. AoR estimation methods generally fall into four distinct categories as summarized in Table 1. The first are the methods based on the Fig. 1. Illustration of a typical CT scan geometry, showing key scanner components (bold), basic geometric parameters (italicized), and different views of raw and processed data, which are relevant for AoR estimation. Fig. 2. Simulated image reconstructed from projections from a fullscan (360 ∘ ) and halfscan (180 ∘ ) range, showing the effects of double edges and tuning fork artifacts, respectively. The 257 × 257-pixel test image was created in Matlab by generating, smoothing, and thresholding uniformly distributed random values. Both projection data (scanned in 0. 5 ∘ increments) and tomograms were created using the ASTRA toolbox [12,13]. Table 1 Overview of automatic AoR estimation methods grouped by category. Center of mass Opposite projection Sinogram symmetry Tomogram evaluation registration evaluation Azevedo et al.[16] (1990) Olander[25] (1994) Liu, Malcolm[26] (2006) Brunetti, De Carlo[27] (2004) Hogan et al.[18] (1993) Pan et al.[17] (2012) Patel et al.[28] (2008) Walls et al.[24] (2005) Jun, Yoon[29] (2017) Yang et al.[30] (2015) Liu[31] (2009) Donath et al.[32,33] (2006) Li et al.[34] (2010) Dong et al.[35,36] (2013) Yang et al.[37] (2012) Yang et al.[38] (2017) Yang et al.[23] (2013) Cheng et al.[39] (2018) Vo et al.[40] (2014) Zhou et al.[21] (2021) Meng, Wu[41] (2017) Lin et al.[42] (2019) Ma et al.[43] (2020) Vo et al.[44] (2021) Vacek, Jacobsen[45] (2022) M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 2
center-of-mass or related concepts in the projections, where a welldefined point is identified in the projection data, tracked throughout the scan, and aligned with the AoR position. The second cateegory consists of opposite projection registration approaches in which two projections taken at opposite views are aligned and their displacement is used to estimate the AoR. Sinogram symmetry evaluation methods also use redundancies in the acquired data but include additional processing steps in order to estimate the AoR in a wider range of geometries. The last category includes tomogram evaluation methods which evaluate reconstructed CT slices with various assumed AoR positions. Each of these categories is discussed in section 2 and category-specific limits and conditions are also compiled and outlined. A practical example of implementing automatic AoR estimation methods is also included in section 3. This example is presented in the context of a laboratory sub-micron CT scanner with a quasi-parallel geometry. 2. Theory 2.1. Center-of-mass methods Center-of-mass (CoM) methods include some of the oldest AoR estimation procedures in the literature [16], with several authors citing their relative robustness to noise [16,29] and fast computation times [16] as their main strong points. Azevedo et al. [16] are among some of the first authors to describe a marker-free AoR alignment method. The center of mass in their approach is calculated as a weighted mean of pixel positions of a projection, with the weights being the absorption values recorded by the corresponding pixels [16]. It is important to perform this calculation after applying a logarithmic transform [18,29] to linearize the raw X-ray intensity data acquired by the CT scanner according to the Lambert-Beer law. The resulting point is then treated as a fiducial marker which can be tracked throughout the scan. A sine curve is fitted onto the detected points and its offset can be used to estimate the AoR [16]. Hogan et al. [18] present a similar approach for correcting the per-projection movement, which is then adapted and adjusted by Olander [25] to specifically estimate the AoR misalignment. Jun and Yoon [29] published a more recent method where they used the center of attenuation (CoA) to detect a fixed point in parallel-beam projection data and align this point from each projection to a common line [29]. The CoA is a related concept to the center of mass, but it is derived specifically for objects expressed as functions of the mass attenuation coefficient rather than relying solely on the measured X-ray absorption values [29]. CoM methods are only suitable for parallel-beam (or quasi-parallel) scans with a circular trajectory [31]. They cannot be used in the presence of lateral data truncation where the sample is not entirely contained within the scan FoV and extends past the left or right edge of the projection data [16,29,30], which is a common occurrence in applications such as region-of-interest tomography [6]. A heterogeneous detector response [25,27,32], X-ray beam instability [27], and nonlinear effects such as beam hardening [16,40] will also negatively impact the results of methods in this group. The performance of CoM is also impaired if diffraction and refraction are not negligible [21,40]. Systematic single-pixel errors such as defective pixels can also skew the results [16]. Most of these methods are additionally hindered by low contrast [16,27] and some sources [25,40,27,32] have called the robustness to noise of this category into question. 2.2. Opposite projection registration Opposite projection registration (OPR) approaches date almost as far back [25] as the CoM category in the literature. These methods assume that two projections taken 180 ∘ apart are mirror images which can be aligned when one of them is flipped along the direction perpendicular to the AoR. The speed of these methods is among their biggest strengths [25]. The OPR methods are also relatively robust to noise according to some sources [25,30]. Olander [25] compares two parallel-beam projections using the mean square difference (MSD) metric, which can be implemented in both the spatial domain and the phase component of the Fourier spectrum. An optimization scheme is applied to find the minimum MSD which presumably corresponds to the optimal AoR position [25]. Pan et al. [17] use cross-correlation (in the spatial or Fourier domain) to estimate the shift between two opposite projections. They then compensate the dataset by half the found shift to adjust the AoR position [17]. Yang et al. [30] suggest several other AoR estimation methods in the OPR category, including the phase correlation or registration using image keypoint detection and description algorithms [30]. The principles of OPR can also be applied to other than parallelbeam geometries with some additional conditions and processing. Fanbeam projection datasets can be rebinned (rearranged) to parallel-beam projections [25], but this is only possible in the central slice of the conebeam data scanned along a circular trajectory [46]. More advanced concepts can also be applied to accommodate the cone-beam data and specialized trajectories such as helical CT [47]. These concepts include differences along PI-lines [47] (parametric interval lines [48]) or epipolar consistency conditions [49] in the projection data. All of these concepts generally require more than a single pair of projections to find matching opposite datapoints due to the cone-beam geometry. This means that although they can be used to create a synthetic pair of projections for use in an OPR alignment, they are more often applied in the SSE approaches discussed in section 2.3 to further leverage information in the data at all scan angles. A large amount of methods related to these concepts has been used to solve the more general problem of motion correction [47,49]. The works of Patel et al. [28] and Liu [31] apply principles very close to the OPR in the fan-beam and cone-beam geometries. Both are based on a direct comparison of individual opposite datapoints [31] or entire projections [28]. However, these methods operate primarily on sinograms and use data at all available angles instead of only one or several projection pairs. The overall characteristics of these approaches correspond the best with the SSE category and both methods are thus further described in section 2.3. Most OPR methods are limited to parallel-beam CT geometries [34] and geometries where the sample size is negligible compared to its distance from the radiation source [25]. Typical cone-beam geometries may cause a noticeable difference in the geometric magnification between sample features closer to the source and features further away from it. This can cause two opposed projections to be different enough so that their registration fails which leads to an incorrectly estimated AoR. Using only a subset of the projection data (usually only two projections) makes the OPR methods vulnerable to sample movement [17,30,40] and to fixed-pattern optical defects such as unresponsive pixels [40]. Insufficient contrast [40] or a lack of distinct features [32] can also negatively influence the results of the OPR. Some sources also dispute the robustness of OPR to noise [38]. The accuracy of the results also suffers if the two registered images are not exactly 180 ∘ apart, which is a real possibility depending on the scan setup [11]. Nevertheless, most AoR estimates tend to be located at most within several pixels of the true AoR position unless scan conditions are extremely adverse [17,30]. This makes the OPR methods suitable as a quick and rough first estimate before applying a more sophisticated algorithm [25,17,30]. 2.3. Sinogram symmetry evaluation Sinogram symmetry evaluation (SSE) methods exploit the periodic nature of the projection data through comparisons between symmetrical/redundant data points. This category expands on the concepts introduced by the OPR methods by employing additional processing in the sinogram domain rather than operating on individual projections. This makes the SSE much more suitable for fan-beam and cone-beam geometries than the OPR. M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 3
This is the most diverse of the four groups as it features methods which operate in both the spatial and frequency domains. The spatialdomain SSE can further be split based on whether an algorithm employs any identification of points or areas of interest in the sinogram. This leads to a total of three sub-categories: segmentation-based SSE, comparison-based SSE, and frequency domain SSE. Methods in the SSE category are generally fast [31,34,40] and many of them are inherently well-equipped to deal with noise and other fluctuations [31,37] by using data from all available views during the AoR estimation. Segmentation-based SSE methods simplify the AoR estimation problem by transforming the input sinogram into a reduced binary representation, which leads to simple and easy-to-implement AoR estimation algorithms. The performance of these methods suffers in lowcontrast data [26,31] and in the presence of substantial unsharpness [26] or lateral data truncation [31]. These methods also require the input data to be full-scan and they will generally fail when applied to shorter scan protocols [26,34]. Liu and Malcolm [26] locate the leftmost and rightmost edges of a sinogram using edge detection and estimate the AoR location as the midpoint between these two extrema. Li et al. [34] take an approach in which the sinogram is binarized using an algorithm such as Otsu’s method [50]. The mean of ray positions above the threshold is then taken as the AoR estimate [34]. Ma et al. [43] combine the segmentation-based and the comparison-based SSE into a hybrid method. A rough AoR estimation is accomplished by segmenting the input sinogram using the mean sinogram value as the threshold and localizing the center of the segmented data corresponding to the scanned sample. The input sinogram is then split into two 180 ∘ sections, which are correlated to find a refined estimate of the AoR [43]. Approaches that can be categorized as the comparison-based SSE are summarized below. These methods process the input sinogram data itself rather than converting them to a simplified representation. They generally evaluate differences of symmetric datapoints for a given AoR using appropriate similarity measures. The correct AoR position is estimated by optimizing these measures. This subcategory also assumes full-scan input data [37,23,40,42] and it can be sensitive to misalignment of the components of the scanner [23]. It tends to offer better performance in low-contrast [31] or noisy [37] data. Lateral truncation is also not an issue here. In fact, comparison-based SSE methods can operate only on a narrow strip of the sinogram as long as the AoR position is contained within this strip [31,42]. However, these methods can still fail when there are not enough prominent features in the sinogram data [37]. The SSE methods are also sensitive to systematic errors such as ring artifacts in practice. This is because erroneous differences between corresponding data points reduce the similarity of these points, which comparison-based SSE relies on. Liu [31] notes the redundancy of measurements in fan-beam full scans and derives a method related to the PI-line principle [47,48] mentioned in section 2.2. An X-ray takes the same path through the sample twice in a fan-beam full scan, only in opposite directions. The positions of these corresponding measurements in the projections can be calculated for a certain AoR position. A sum of squared differences of all these redundant measurements is then computed to evaluate the suitability of each AoR estimate. Minimizing this score with respect to the AoR position leads to an estimate of the true AoR [31]. Patel et al. [28] compare opposite projections in a cone-beam geometry. They note that the opposite projection pairs will show differences due to the cone-beam magnification, but posit that this should not affect the estimation [28]. The method is based on selecting a single row (sinogram) of a cone-beam dataset and comparing all opposite pairs of projection vectors within it. The comparison is performed by flipping one of the vectors, shifting it by different amounts, and finding the shift that minimizes the root-mean-square difference between the pair of vectors. An average of shifts calculated for all opposite pairs is assumed as the AoR estimate [28]. This process can be repeated for a range of rows in the data, which the authors use to additionally estimate the AoR tilt [28]. Yang et al. [37] opt to rebin the full-scan fan-beam sinogram data to a parallel-beam geometry. They then split the rebinned sinogram into two 180 ∘ halves, mirror the second half-sinogram to reverse the projections within, and apply the cross-correlation to find shifts between corresponding projections of the two halves. The mean value of these shifts is then divided by two, leading to an estimate of the distance between the current midpoint of the dataset and the estimated AoR [37]. Yang et al. [23] further streamline this approach by omitting the rebinning step and summing all the projection data in a fan-beam sinogram to produce a single vector of values. This vector is ideally symmetric around the AoR position. The autocorrelation of this vector reveals the amount of shift needed to align the dataset’s midpoint with the AoR [23]. Meng and Wu [41] investigate the fan-beam geometry and come to a similar conclusion as Liu [31] in terms of the redundancy of the fullscan data. A pair of redundant rays at a certain distance has a welldefined angular interval between them. This interval is exactly 180 ∘ only for the rays intersecting the AoR. The sinogram is thus split into two 180 ∘ halves and the cross-correlation coefficient of each corresponding projection location from both halves is calculated [41]. Applying the cross-correlation on columns of the sinogram instead of its rows allows the authors to avoid the rebinning step necessary in the method of Yang et al. [37]. Maximizing the value of the correlation coefficient yields the estimated AoR position [41]. Comparison-based SSE methods are also suitable for specialized scanning protocols such as offset scans [42], where the detector or the AoR is deliberately displaced to one side of the projection data and additional processing is applied to the dataset to extend the FoV [42]. Lin et al. [42] present a method for such offset-scan data with a displaced AoR. The central ray (the ray intersecting the AoR) is no longer perpendicular to the detector in such scans. This causes most of the existing AoR estimation methods to fail [42]. The method bears resemblance to the authors’ earlier work [23], but it adds a transform of the data onto a virtual detector that is perpendicular to the current AoR estimate. The transformed data on both sides of the AoR are summed to produce two vectors of values that are expected to be mirrored around the AoR. The variance of the difference between these vectors serves as a cost function which is minimized to estimate the correct AoR position. Vo et al. [44] employ an approach that is slightly different to the previous ones, as it is tailored to the parallel-beam offset-scan data. The full-scan sinogram is divided into 180 ∘ halves and the second half is flipped. The overlap between both halves is then determined by optimizing an image correlation metric of narrow regions in both half-sinograms. The half-sinograms are then stitched into one wide half-scan sinogram with the AoR positioned in the middle of the resulting dataset [44]. The frequency domain SSE is a small subcategory which nevertheless offers some unique properties compared to the previous categories. These properties include possibility of use with half-scan data [40] and an effective performance in noisy or otherwise adverse conditions [40,45]. The method of Vo et al. [40] is designed specifically for the parallel-beam half-scan data. A half-scan sinogram is converted to an artificial full-scan one by copying and flipping the sinogram and concatenating the original and the flipped copy. The AoR misalignments will cause discontinuities in the stacked sinogram and corresponding frequency components in its two-dimensional frequency spectrum. These components can be easily distinguished from the frequency coefficients belonging to the genuine structures of the sinogram based on their location within the spectrum. Minimizing the total magnitude of frequency components in these locations by shifting the two sinogram halves with respect to each other will lead to an estimate of the AoR position [40]. Vacek and Jacobsen [45] sum together opposite projections in a full-scan sinogram to produce even vectors symmetric around the AoR. The misalignment of the AoR forms a ramp in the phase component, which can be converted back to a shift value of M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 4
the AoR in the spatial domain. The influence of noise on the algorithm is minimized by examining only the low-frequency components, but the method can fail in the presence of low-frequency fluctuations ir if the projections are truncated [45]. 2.4. Tomogram evaluation Tomogram evaluation (TE) methods model the behavior of a human operator during manual AoR estimation. They score reconstructed cross-sectional images using metrics which evaluate changes in the data quality caused by the AoR shifts. Optimization of such metrics ideally yields the true AoR position. Operating on the CT slices makes the TE methods independent on the used acquisition geometry [32]. The TE methods also tend to be robust to fluctuations in the projection intensity [32]. Brunetti and De Carlo [27] observe that the AoR misalignment smears intensities of datapoints over neighboring pixels. They use the number of non-zero (non-background) points in the tomogram as a metric to minimize in order to estimate the AoR. Walls et al. [24] and Dong et al. [35,36] view the AoR misalignment as a cause of blur which decreases the variance of image values. Their method is thus based on maximizing the variance of values in the reconstructed tomogram as a function of the AoR position [24]. This metric is sensitive to truncation artifacts on the edges of the FoV [24], streaks caused by limited views [21], and other image artifacts which artificially increase the image variance. Zhou et al. [21] expand upon the use of this metric by first segmenting the image into background and foreground and only computing the metric in the latter. The authors report that this improves the performance of the metric in measurements with a poor signal-to-noise ratio and a limited number of projections [21]. Donath et al. [32,33] suggest three distinct metrics which all yield an estimate of the true AoR position when minimized. The discrete versions these metrics are: the sum of absolute values of pixels in the tomogram, the sum of negative values, and the entropy of the grayscale histogram of the tomogram. The first two metrics are supported by mathematical proofs in the original work [32] and their use is limited to cases where only positive values are expected in the tomogram data. The third entropy-based metric is a heuristic, but it can also be applied to tomograms with negative values [32]. The robustness of these metrics to noise can be increased by averaging several neighboring sinograms before reconstruction, but the method may still exhibit lower precision if the signal-to-noise ratio in tomograms is too low [33]. Yang et al. [38] consider the AoR estimation as an image classification problem. The task is to label tomograms as aligned or misaligned based on whether they were reconstructed with the true AoR position or not. The authors train a convolutional neural network on a database of patches of manually classified tomograms [38]. Tomograms reconstructed at a range of AoR positions are then presented to the trained network, which predicts the images with features most indicative of an accurate alignment [38]. One limit of this method is that it may have suboptimal results for datasets which have very different characteristics and features from the training data [38]. Cheng et al. [39] evaluate the quality of tomograms using their total variation (TV). Artifacts caused by a misaligned AoR increase the TV, which in turn means that minimization of this value yields an estimate of the correct AoR. The influence of noise on the TV metric is reduced by applying a smoothing filter to the tomogram [39]. The authors also note that the image should contain enough information (quantified using entropy) for this metric to provide reliable results [39]. TE methods can be enhanced by any artifact reduction or denoising algorithms or the use of advanced reconstruction algorithms. However, this further increases the complexity and computation time relative to other AoR estimation categories. The high time-cost of these methods can be especially concerning when dealing with large datasets [34,30,40,42,45], but researchers have used various techniques to reduce the total calculation time. These techniques include the use of optimization algorithms [39] and coarse-to-fine approaches [35,39]. A prealignment using a coarse but fast method can also be used to acquire an initial AoR estimate and a more precise TE method can then be applied only in a narrow range around this estimate [35,43,21]. 2.5. Methods with AoR estimation as a secondary trait Some of the works referenced in sections 2.1 to 2.4 combine the AoR estimation with the estimation and/or correction of other parameters such as the AoR tilt [28,30], angular errors [39], motion of the sample [18,28], or other geometrical parameters [32]. Some authors use properties of the acquired datasets for even more complex and holistic alignment approaches. These algorithms often fall into one of the categories outlined in sections 2.1 to 2.4, but they perform a multivariate optimization of the scan geometry instead of localizing only the AoR position. It is common for such methods to estimate three and more parameters such as the source-object distance, slant of the detector, and vertical misalignments. Examples of such algorithms include works by Viskoe [51] (most related to the CoM and TE approaches), Kyriakou et al. [52] (TE), Panetta et al. [53] (SSE), and Kingston et al. [54] (TE). Reprojection-alignment algorithms were first used in electron tomography [55] and later applied to CT [56]. This group of methods is closely related to the AoR estimation, but distinct in several aspects. These approaches iteratively correct the per-projection drift and the sample movement [14] or misalignments in the scan geometry [57]. There is a large variety of different approaches in this group [22], and their primary application is tomography with resolutions in the order of tens of nanometers and less [15,22], where the sample can exhibit nonnegligible motion throughout a scan. The latest of these approaches jointly reconstruct, reproject, and compare data to the original projections to estimate and correct the drift of samples [22,58] and in some cases even their deformation [8] in a single calculation framework. These reprojection methods can also correct systematic shifts of the AoR, essentially a special case of sample drift, as a by-product of the iterative process [22]. Another notable alignment method presented by [59] uses comparison of image features in opposite projections for estimating the misalignment of the AoR pitch, roll, and position. The scale-invariant feature transform [60] is used in this work, but other feature descriptors and extractors [61] can also be used [59]. Features extracted from pairs of the opposite projections are matched using a brute force comparison. The AoR position and other misalignments can then be deduced from spatial relationships of matching features [59]. This method differs from the post-scan methods in section 2.2 because it is used prior to a scan to iteratively adjust the sample stage and ensure the alignment of the hardware itself. 2.6. Methods aimed primarily at other modalities Yang et al. [62] applied an AoR estimation method in computed laminography using an approach similar to the authors’ previous work aimed at CT [23]. This approach is not strictly limited to computed laminography, but the oblique scan geometry requires a slightly different theoretical approach when compared to CT [62]. All projection data are summed to create a single image whose rows form even vectors symmetrical around the AoR. The autocorrelation of the rows yields an estimate of the AoR position [62]. The AoR estimation is also essential in transmission optical projection tomography (OPT), an imaging modality closely related to CT which uses visible light instead of X-rays [24]. Key differences between CT and OPT include the much more significant effect of refraction and scattering and a shallower depth of field in the latter [24]. Several automatic AoR estimation methods aimed primarily at OPT have been published [24,63,36,35,43]. Many of these are also applicable in CT, and so they are included in the overview in sections 2.3 and 2.4. M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 5
3. Implementation of axis-of-rotation estimation Different automatic AoR estimation methods vary in robustness, processing time, and requirements placed on the input data. The choice of an appropriate method to be implemented in a CT data processing workflow must be considered carefully and in context. This example shows how AoR estimation methods may be implemented in the context of a laboratory sub-micron CT scanner with a quasi-parallel geometry. This intermediate geometry combines traits of both parallel-beam and cone-beam CT. It is a relatively under-represented geometry in the existing literature on the automatic AoR estimation, so this example provides a useful framework for providing additional and complementary insight to previous publications. The Rigaku nano3DX scanner [64] used in this example utilizes a quasi-parallel geometry with a powerful source at a large distance from the sample and a relatively short sample-detector distance. The resulting X-ray beam has a very shallow cone angle and a low geometric magnification of around 1.02 × for most regular measurements. X-rays are detected by the XSight Micron LC X-ray CCD camera. Rays are converted to visible light by a scintillator and then magnified onto the detector array using microscope optics. The scanner is equipped with a MicroMax-007 HF X-ray source with interchangeable copper (Cu) and molybdenum (Mo) targets operating at 40 kV/30 mA and 50 kV/24 mA, respectively. Measurements with the Mo target use a 0.1 mm Al filter and no filtering is applied in Cu measurements. 3.1. Model dataset selection Five model cases (Table 2) were chosen to cover a range of typical applications of the used scanner. The collection features datasets with low contrast and high noise (case 1), homogeneous objects (case 2), beam hardening (case 3), significantly truncated projections (case 4), as well as a dataset without any significant complicating characteristics (case 5) (Fig. 3). There is also the possibility of sample movement that is noticeable over the duration of a scan but negligible in-between projections. All test datasets contain 800 projections with 1648 × 1250 pixels each, which are scanned over a half-scan range (180 ∘ ). 3.2. Positional accuracy of the AoR estimate AoR estimation methods need to be accurate, which means that the estimate of an effective method must be close to the true value. However, there is no universally accepted level of the AoR accuracy. Some sources regard the accuracy of 0.5 pixels as sufficient [17], while others cite a shift of as little as 0.4 pixels unacceptable [23,37] and aim for a maximum deviation of several hundredths of a pixel [32]. This might be because the magnitude of tuning fork artifacts is fundamentally influenced by aspects such as noise, contrast, sample shape and movement, reconstruction algorithm settings, and more. An appropriate accuracy level must therefore be determined on a case-by-case basis. For this example – a sub-micron CT with a quasi-parallel geometry and relatively low-contrast data - a maximum error of 0.5 pixels was deemed appropriate based on practical tests (Fig. 4) as well as some of the literature [17,27,24]. 3.3. Tested estimation methods The selection of suitable AoR estimation methods was first narrowed down based on the characteristics of the used CT scanner. The high probability of truncation due to the small FoV of sub-micron CT made the CoM category impractical. Most SSE methods were also excluded due to their incompatibility with half-scan data. Four appropriate methods were finally selected for further evaluation: . •M1, an OPR method which uses the gradient cross-correlation [40,65]. •M2, the SSE method of Vo et al. [40] in which a sinogram is duplicated, flipped, and stacked onto the original with a varying AoR. Spectral analysis of this stack reveals discontinuities between the two sinograms which are ideally minimal at the true AoR. Table 2 Overview of model cases for testing selected automatic AoR estimation methods in the case of the nano3DX CT scanner. The used X-ray source has fixed settings for each target material, which are mentioned in section 3. Case Sample Target material Exposure Voxel size FoV width Cone angle 1 Fiber-reinforced polymer Mo 8 s 0.53 μm 873 μm 0.09 ∘ 2 Ruby ball Mo 10 s 1.03 μm 1697 μm 0.18 ∘ 3 Li-ion battery cathode Mo 40 s 0.53 μm 873 μm 0.09 ∘ 4 Toothpick Cu 23 s 0.53 μm 873 μm 0.09 ∘ 5 Foraminifera micro-fossil Mo 20 s 0.53 μm 873 μm 0.09 ∘ Fig. 3. Sample CT slices of model datasets shown with enhanced contrast. Case 1 and its detail is shown in (A) and (F), case 2 in (B) and (G), case 3 in (C) and (H), case 4 in (D) and (I), and case 5 in (E) and (J). M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 6
•M3, a TE method based on the sum of absolute values [32], which iteratively reconstructs CT slices in order to minimize a cost function. The sum of absolute value was chosen as the cost function over others in [32] based on the authors’ recommendation. •M4, a TE method based on the total variation [39] which operates in a similar manner to M3. Other TE metrics were omitted as they are based on similar concepts as M3 or M4. Methods M1 and M2 are suited for the parallel-beam geometry, but their performance in a quasi-parallel geometry is unclear. The methods are still expected to output estimates close to true AoR positions, but their precision may be impacted. The results of M2 can also be influenced by the amount of truncation in case 4 in particular. Methods M3 and M4 are not expected to be impacted by the geometry, but they may be affected by some of the artifacts present in the test cases due in part to the lack of any additional processing applied to the data. These methods are also expected to take significantly longer to compute than M2 and especially M1. Method M1 requires only a single-pass calculation, while methods M2, M3, and M4 optimize their result iteratively. A simple grid search in 0.5-pixel intervals in a sufficiently wide AoR range was used for the iterative methods. Data were flat-field corrected in the nano3DX acquisition software and all further processing was performed using Matlab R2020a. All estimations were run on a Windows 10 machine with a six-core AMD Ryzen 5 2600X CPU, an nVidia GeForce GTX 1050 Ti GPU, and 64 GB of RAM. Tomograms were reconstructed via the Feldkamp-Davis-Kress algorithm [66] implemented in the ASTRA toolbox 1.9 [12,13] with a cosine filter cut off at 0.85-times the maximum frequency. No other processing (such as artifact reduction) was applied on the tomograms. Minor adjustments were made in the implementation of the selected methods compared to their original publications in order to make the testing seamless. These and other details are described in appendix A. Reference AoR values for each case were obtained as median values of a manual AoR estimation performed by ten experienced operators who were checking with the 0.5-pixel accuracy level determined above. Estimates of the individual operators are listed in appendix B. Results of the tested automatic methods were then compared to the reference values while considering the target accuracy level. The manual estimation and all automatic methods except M1 were performed on the sinogram closest to the center row of the detector array [37]. This sinogram contained adequate information in all test cases, so a more sophisticated choice of slices such as the one described in [36] was not necessary. 3.4. Estimation results and discussion Experimental results (Fig. 5) show that none of the tested methods were robust enough to handle all five test cases with a sufficient accuracy. M1 and M4 exceeded the 0.5-pixel threshold in case 5 and M2 and M3 did not reach sufficient accuracy in case 4. M2 also deviated from the reference AoR in case 1. Figure 6shows results of the tested methods for case 4. Results of M1 and M4 are within the ± 0.5-pixel accuracy window and show no obvious differences from the reference. The outcome of M2 is close to the reference and contains no noticeable tuning fork artifacts, but changes to some of the sample structures can be observed. The result of M3 is severely distorted. The results reveal some interesting and not immediately obvious properties of the tested methods. The truncation in case 4 caused M2 to deviate from the reference. This was expected as the performance of the method is tied to the sample size [40]. However, M2 also reached the second largest overall error in case 1, presumably due to the low contrast and high noise. The performance of M2 and M1 in the rest of the cases shows that the methods are indeed compatible with quasi-parallel geometries. Tomographic artifacts were a general concern for M3 and M4. The truncation in case 4 had the most impact on M3 in particular. This was surprising at first, but it does in fact seem to be consistent with the original publication describing M3 [32]. The authors mention that the Fig. 4. Influence of a misaligned AoR in case 2. Effects of a positive and negative AoR misalignment are highlighted by red arrows. Only very faint streaks with minimal impact on subjective image quality are present with an AoR shift of 0.5 pixels, but at a shift of 1 pixel streaks are already noticeable. Fig. 5. Errors of each tested method in all test cases. A bar ending in an arrow indicates the value is outside the displayed range. Green lines indicate the window of values that can be assumed as accurate enough ( ± 0.5 pixels). M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 7
used cost function is suited only for measurements with strictly positive attenuation coefficients [32]. This renders it ineffective for measurements of samples in a surrounding medium like water. The results shown here and in appendix Csuggest that truncation in the projection data may have a similar effect. The differences in the computation time between the four methods were as expected. The two TE methods were computed simultaneously, taking 26.8 s to process when averaged over five runs. Most of this time was taken up by the reconstruction. The more lightweight M2 ran for an average of 13.38 s. Computation of the single-pass alignment used in M1 was the fastest at an average of 0.64 s. While these run-times are indicative of the relative time demands of each method, they can still be significantly decreased through mathematical and code optimization. 4. Outlook The results suggest that a general and fully autonomous application of any tested method is still limited, as the performance of these methods depends strongly on the character of the input data. Different methods vary in their strengths and weaknesses. A compromise between the automation, processing speed, and reliability can thus be reached by a hybrid AoR estimation approach using multiple methods. This has also been suggested in the prior literature [17,30]. Another possible approach may be to base the final AoR estimate on a consensus between the results of all applied methods. Depending on the required level of automation and processing throughput, a manual check can be added to leave the final decision up to a human operator. For a laboratory setting with a lower volume of processed samples, this is an acceptable compromise between the convenience and objectivity of automatic methods and the control of a manual estimation. The number of AoR estimation methods in the recent literature skews heavily towards the SSE and TE categories. Advances in the computer hardware mean that even demanding operations like tomographic reconstruction can be performed relatively quickly. This makes practical use of the demanding TE methods feasible. At the same time, the SSE and TE categories tend to be the most versatile in terms of scan geometry and the character of the scanned data, so it is likely that they will see more use in the future. The TE methods in particular can employ advanced image quality metrics to better replicate the evaluation of CT data by a human observer. The use of artificial intelligence may also become more popular for AoR estimation, as its is a major trend in most other areas of image processing. An example of this is the convolutional neural network used by Yang et al. [38] (mentioned in section 2.4). The increased capabilities of modern computers also favor methods for motion correction and algorithms for multi-parameter geometric alignment 2.5. These methods tend to be more demanding than the specialized AoR estimation, but they also offer much more flexibility. This is especially true in non-standard scan geometries and high-resolution applications where the sample motion is a major issue. 5. Conclusion The alignment of the axis of rotation in a CT scan is essential for a high-quality tomographic reconstruction. The estimation of this position directly from the scanned data is often the most expedient approach to achieve this alignment. The laborious task of estimating the AoR position can be automated using a range of published methods which fall into four distinct categories with specific strengths and limits. These categories include center-of-mass methods, opposite projection registration, and sinogram symmetry evaluation, which all operate on the projection data or sinograms. The final category is tomogram evaluation. The choice of an appropriate AoR estimation method for a particular application should be carefully considered based on the characteristics of potential methods. Most published approaches assume a specific scan geometry or angular range. The tomogram evaluation category is unique in this regard as it avoids any such assumptions by operating on the reconstructed data. However, it is also the most computationally demanding category. The registration of opposite projections appears to be a very quick and robust approach for parallel-beam and quasi-parallel geometries. The sinogram symmetry evaluation and center-of-mass methods are more specialized, making them more susceptible to fail when assumptions about the input data are not fulfilled. A truly universal and robust automatic AoR estimation method does not seem to exist at the moment, so the combined estimation of multiple different methods can be leveraged for an increased robustness. A basic AoR estimation may also be inadequate in non-standard CT applications. For instance, high-resolution nanoCT scans show significant amounts of movement between consecutive projections. A more comprehensive correction of any motion of the sample and geometric misalignments, such as one of the reprojection-alignment methods, might be more appropriate in such cases, despite the added computational complexity. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 6. Reconstruction of the dataset in case 4. The same cut-out area as in Fig. 3 (H) is shown here in (A). Details in (B) to (E) of the area marked by a red square show different AoR positions estimated manually (B) and by the tested methods. Arrows point to areas where differences between C-E are the most prominent. M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 8
Acknowledgments Funding: We acknowledge CzechNanoLab Research Infrastructure supported by MEYS CR [grant number LM2018110]; and the Brno University of Technology [grant numbers FSI-S-20-6353, CEITEC VUTK-21-6956]. We thank members of the Laboratory of X-ray micro and nano computed tomography for providing the test datasets for the purposes of this work, and for lending their time to manually estimate the AoR positions of the test datasets. We also thank Martina Pořízková and Eva Zikmundová for proofreading various stages of the manuscript. Appendix A. Implementation Details for AoR Estimation Methods This section summarizes details of the implementation of automatic AoR estimation methods used in this study. The basic structure of the implementations is summarized below in the form of pseudocode. There is a number of changes from the original descriptions of these methods in order to ensure proper functionality in the context of the subsequent tests. These are highlighted in bold and briefly described. Processing and registration are flexible in method M1 (algorithm A.1). In this implementation, median filtration is applied to increase robustness to noise, and a 2D Tukey window is applied to make the subsequent cross-correlation more reliable. Gradient cross-correlation [40,65] was chosen as a fast and reliable registration method. Algorithm A.1. for method M1. For method M2 (algorithm A.2), steps were taken to increase the robustness of the algorithm for nanoCT data. A different σparameter was chosen for the Gaussian filter than in the original publication. A row normalization step was also added; each row of the sinogram is normalized to reduce the effect of a changing projection intensity throughout the scan, which can produce false discontinuities in the stacked sinogram. Algorithm A.2. for method M2. A simple grid search optimization was used for M3 and M4 in algorithm A.3 to simplify the code compared to the original publications. Grid search with a step of 0.5 pixels was sufficiently precise for the purposes of testing, not overly time-demanding, and potentially more robust to local optima than more sophisticated optimization methods. Compared to the original description in [32], the metric used for the method M3 was simplified in a way that does not affect the results. A simple FBP-type algorithm with a Ram-Lak filter was used for the reconstruction. M. Zemek, J. Šalplachta, T. Zikmund et al. Tomography of Materials and Structures 1 (2023) 100002 9