Registration-based multi-orientation tomography

Wepropose a combination of an experimental approach and a reconstruction technique that leads to reduction of artefacts in X-ray computer tomography of strongly attenuating objects. Through fully automatic data alignment, data generated in multiple experiments with varying object orientations are combined. Simulations and experiments show that the solutions computed using algebraic methods based on multiple acquisitions can achieve a dramatic improvement in the reconstruction quality, even when each acquisition generates a reduced number of projections. The approach does not require any advanced setup components making it ideal for laboratorybased X-ray tomography. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
X-ray computed tomography (CT) is employed in a wide range of fields including clinical imaging, biomedical research, materials sciences, industrial non-destructive testing and cultural heritage [1][2][3][4][5].In some cases, specialized and highly optimized CT systems can deliver accurate and reproducible three-dimensional reconstructions.Nonetheless, the available CT scanners often cannot be optimized for scanning particular types of specimens, producing sub-optimal reconstructed data.Typical issues in laboratory-based X-ray CT are, for instance, beam hardening [6], scattering [7] and photon starvation artefacts [8] observed in the presence of high density regions inside the specimen.
To improve the quality of CT reconstruction based on sub-optimal data, an ever increasing number of advanced computational approaches are being developed.Such methods include: model based and empirical beam hardening corrections [9][10][11][12], statistical noise models [13], regularized reconstructions [14], machine learning-based approaches [15,16] and approaches based on hardware calibration [17].
In addition to data processing methods, a number of approaches exploit a principle of multi-axis acquisition to improve the quality of reconstruction.For instance, multi-axis approaches for electron microscopy [18,19], X-ray CT [20,21], medical interventional X-ray CT scanners [22,23] accompanied by a number of reconstruction techniques designed for scanning orbits of an arbitrary shape [24][25][26][27].
Similarly to multi-axis approaches, in the current paper we address the issue of accurate CT reconstruction avoiding the complexities and ambiguities present in the processing-based approaches mentioned above.Namely, avoiding the need for algorithm parametrization, prior knowledge of the specimen's properties and correct statistical models of the detection process, spectrum calibration procedures and advanced hardware.We also aim at avoiding the necessity in complex hardware such as gimbal-type sample holders or other high-precision multi-axis rotation devices used in the other multi-axis-based approaches.Alternatively, we use a volume registration approach that allows to acquire multi-axis data through manual orientation of the object inside the scanner.Instead of explicitly addressing the nonlinear nature of image formation in laboratory X-ray CT systems, we propose to explore an experimental approach based on the acquisition of additional tomographic datasets, where each dataset corresponds to a different orientation of the object (i.e.rotating around a different axis).Each additional dataset can have a reduced number of projections or shorter exposure times to make sure that the total dose stays the same as in the case of a single acquisition.This approach is coupled with a robust volume registration method and a parameter-free algebraic reconstruction approach.The latter is based on penalizing data with higher degree of non-linearity or statistical error.

Tomographic acquisition
In most applications, X-ray CT reconstruction algorithms are based on variants of the well known Radon transform.While that model is accurate in case of monochromatic illumination, applying it to the experimental data generated by a typical polychromatic X-ray CT scanner is associated with both systematic and random errors.Instead of relying on the known pre-and post-reconstruction error-correction methods, in the current work, we will explore a more data-driven, linear approach.We will, however, assume that a significant fraction of measurements conform well with the linear model, while a minor part of measurements are associated with larger systematic and random errors and should be penalized (see Fig. 2(a)).Several CT scans with different orientations of the rotation axis should be performed to increase the number of available measurements that conform with the linear model (see Fig. 2(b)).We will assume that the object is rigid (it is not deforming during the mount) and it is possible to change the orientation of the rotation axis by simply changing the objects orientation inside the scanner.
To avoid a need for complicated gimbal-type sample holders or multi-axis rotation stages, in each scan, the object will be reoriented manually.The exact coordinate transformation that corresponds to different orientations of the object can be found by a volume registration algorithm described in the Section 2.3.After the volume registration step, a reconstruction based on all available data can be computed using either an explicit approach based on Filtered-Back Projection (FBP), Feldkamp, Davis, and Kress (FDK) or algebraic approach based Least Squares minimization.A flow chart of this approach is shown in Fig. 1.

Tomographic reconstruction
Let us start by describing the CT reconstruction problem for multiple datasets that correspond to different orientations of the object.We describe a projected attenuation p as a function of the object linear attenuation µ in case of a cone-beam geometry through the X-ray transform: here the projected attenuation p is defined for each source location x ∈ Γ, where Γ is the trajectory of the source and θ is a unit vector pointing towards a detector pixel (see Fig. 2).According to the Beer-Lambert law, projected linear attenuation p( x, θ) is computed from the transmitted and incoming intensity (I out and I in respectively): Once the coordinate transformation describing the change in the object orientation inside the scanner is known, a new equation for relating the measured data to the object attenuation can be defined as: here T is a rotation transform and s is a translation vector.In a discrete form, the linear transformation defined by Eq. (3) can be described by a matrix or evaluated on the fly.Assuming that all measurements are described by a discrete vector p and all points of the reconstructed image correspond to a vector µ, Eq. (3) can be written in a short matrix notation p = Rµ.This problem is commonly solved using a L2 minimization: To make sure that the measurements with high error magnitude are penalized, we will consider two minimization approaches: the so called Penalized Weighted Least Squares (PWLS) [28] and an approach based on Student's-t penalty (SP) [29].The PWLS approach can be expressed as a minimization problem of the following form: where W is a diagonal weight matrix with elements proportional to the error variance.We will assume that both the magnitude of the systematic error due to non-linearity (i.e.beam hardening, detector response) and the variance of the random error (i.e.shot noise) are highest in the measurements with high attenuation [21].So, as a first approximation, we will assume that the elements of the weight matrix W are inversely proportional to intensity, and In the approach based on Student's-t penalty, the minimization problem expressed by Eq. ( 4) is replaced by the following minimization problem: where summation is carried out for all elements of the residual x; σ is an error variance parameter which is estimated in [29] using the following linear search: where M is the number of elements in x.
For the sake of comparison we included a third reconstruction approach: an explicit inversion approach first presented in [21].It does not rely on solving a minimization problem but instead fuses reconstructed volumes computed by the FDK algorithm using weighted summation.The weights are computed by back-projecting the estimated error variance that is inversely proportional to the measured intensity (similarly to our work).

Volume registration
In a CT system where the orientation of the rotation axis can be controlled, an optimal combination of axes orientations can be pre-computed using optimization methods similar to the one described in [30].In the current work we address the case of a standard laboratory setup, where the rotation axis cannot be adjusted and the user has to simply reposition the specimen on its mount by hand.Typically, two or three scans can be carried out with object orientations close to mutually orthogonal.If it is possible to assess the directions of dominant attenuation (based on previous scans or based on the shape of the object), it is sensible to align the rotation axes parallel with such directions.This would lower the apparent attenuation and thus reduce the amount of measurements associated with higher non-linearity.
Assuming that the object does not change between consecutive scans, a rigid transformation can be computed for each volume reconstruction, allowing to register all of the available data in a common coordinate system.While marker-based registration methods (such as [31]) may be considered more robust than intensity-based ones [32], marker-based registration complicates the data acquisition step.To avoid that, we developed an intensity-based volume registration method that, in our experience, delivers reliable results when applied to X-ray CT data.
To make sure that only the relevant object features influence the registration process, some preprocessing can be applied to the initial reconstructions.To avoid influence of low intensity noise with a large spatial footprint, a simple binary threshold can be applied (e.g.Otsu threshold [33]).If the mount is clearly visible in the reconstructed dataset, it may need to be truncated.
The registration step can be subdivided in two parts: initial alignment using image moments [34] and fine alignment using gradient descent-type optimization [35].The first step of the registration is needed to make sure that the gradient-descent optimizer starts near a global minimum of the objective function.Major and minor axes of the spatial distribution of the object attenuation (proportionally, the mass) can be computed using eigenvectors of the covariance matrix constructed from the second order central moments of the reconstructed volume.Central moments are defined as follows: here x 0 , y 0 and z 0 are the center of mass (attenuation) of the reconstructed volume.Subsequently, the normalized moments are: mijk = m i jk /m 000 .The covariance matrix is constructed from the second order normalized moments: The rotation matrix that corresponds to the orientation of the object M q , can be obtained by stacking of the eigenvectors of the matrix C. Subsequently, the relative rotation between two datasets is defined as M 0 M T q , where M 0 is the rotation matrix of the "master" volume.Image moments provide information about the shape of the mass distribution but not the direction (they can be flipped) and can be especially ambiguous in symmetrical objects.To avoid ambiguity, we generate sets of four rotations, each time over an angle of π /2 around one of the three moment vectors, and select one that corresponds to the minimum of a simple objective function: where µ 0 is the "master" volume to which all subsequent volumes µ q are registered; T k is the coordinate transformation that corresponds to the initial moments alignment, plus one of the π /2 rotations.Further registration is done using one of the gradient descent algorithms based on an objective function similar to Eq. ( 12).The volume registration method described above may certainly be replaced with other techniques (see [36] for instance) that are sufficiently stable when applied to CT reconstruction volumes oriented nearly orthogonally or having an arbitrary orientation.

Algorithms and implementation
To implement volume registration, we used the Python package numpy to compute image moments and the Powell optimization with an objective function based on mutual information (see [37]) implemented in the Insight Segmentation and Registration Toolkit (ITK) [38].

Algorithm 1 Volume registration
1: for Each dataset p q do 2: Compute a preliminary reconstruction µ q (FBP / FDK) 3: Apply a binary threshold (Otsu) 4: Compute second order central moments mijk 5: Construct the covariance matrix C and compute its eigenvalues and eigenvectors 6: Construct the rotation matrix from the eigenvectors 7: Find the minimum µ 0 − T k µ q 2 2 for a number of π /2 rotations (initial + 3 per axis) 8: Run Gradient-Descent optimizer for fine registration A description of the regularized version of the PWLS algorithm is given in [28].We implemented an unregularized version of PWLS with the inverse of the estimated variance proportional to the intensity exp(−p).

Algorithm 2 Unregularized PWLS
Pre-compute the weights: In each iteration of the multi-orientation algorithm, an update of the reconstructed attenuation µ is computed for each available dataset p q and the total update is calculated as an average.X-ray transforms R q in shifted and rotated coordinates (estimated for each dataset p q ) were implemented using the ASTRA tomography toolbox [39] which allows to apply an arbitrary transformation to coordinates of any element of the tomographic setup.

Algorithm 3 Multi-orientation PWLS
Pre-compute forward-projection and back-projection weights: for Each dataset p q do ω q = exp(−p q ) s q = R T q ω q R q Iterate: while Termination condition not reached do for Each dataset p q do Update solution: µ k+1 = µ k − 1 s q N R T q ω q (R q µ − p q ) Apply non-negativity constraint: µ k+1 = max(0, µ k+1 ) The multi-orientation Student's-t approach was implemented by replacing the L2 misfit function of the Least Squares by the Student's-t misfit function.The error variance parameter was estimated for each iteration using a downhill simplex algorithm implemented in scipy Python package.

Algorithm 4 Multi-orientation Student's-t
Pre-compute projection weights: for Each dataset p q do s q = R T q R q Iterate: while Termination condition not reached do for Each dataset p q do Compute residual: x q = R q µ − p q Estimate variance: σq = arg min σ M log(πσ) Apply non-negativity constraint: µ k+1 = max(0, µ k+1 ) The source code for the multi-orientation tomography can be found on GitHub at [40].

Data
Data demonstrated in this publication were both modeled and acquired experimentally.Modeled data were produced using the Monte Carlo simulation platform GEANT 4 [41] (see Fig. 3).The developed model uses components of the Geant 4 toolkit to define a geometry of the experimental setup, simulate passage of photons through the matter of the object and simulate the detector response.Physical models in this simulation take into account the photoelectric effect and Compton scattering to calculate direct interaction between X-rays and object matter.Thus, beam hardening, scattering, and photon starvation artefacts can be reproduced.Photons diverge as a cone from a point source.Their energy is distributed according to one of the measured spectra (shown in Fig. 3(b)).In the simulated projections, pixel values are equal to the number of registered photons and do not depend on the photon energy or momentum.Thus, the statistics of the simulated data is only to some extent representative of the statistics of a conventional integrating detector.The following geometry was considered: the source-to-object distance and the object-todetector distance 100 mm; detector pixel size 0.2 mm; 3e8 photons generated per projection.
Flat-field images are generated with 1.5e9 photons per image (five times the photon count of a single projection).The phantom is defined as a Nylon cylinder of 9 mm radius and 12 mm height with an aluminium block 12 × 4 × 12 mm and two iron rods or 0.6 mm radius and 12 mm height (see Fig. 3(a)).
The reconstruction results for 4 simulated datasets are shown in Figs.Monte Carlo simulation allows to take into account all relevant physical effects, for example, a change of the photon trajectory due to the Compton scattering.Thus, pixel values can not be computed using linear integrals along straight lines (i.e.X-ray transform).Monte Carlo methods usually require significant time for the calculations, limiting the statistical accuracy of the data.Due to limited photon count, that falls from around 2000 photons per pixel outside of the object to around 200 photons behind metal absorbers, an error, according to the Poisson statistics, can be roughly estimated as 3 -7 % for a separate projection (depending on the area of projection; photon count decreases behind iron rods in comparison with nylon) and 1.5 % for the flat-field image.
Experimental data were acquired at the FlexRay lab [42], CWI, Amsterdam for three objects: a printed circuit board (PCB), an incandescent light bulb, and a nine-volt alkaline battery (see Fig. 5).All three objects were placed on a polyester sponge (ensuring the mount transparency) in two orientations close to orthogonal.The scan settings can be found in Table 1.As we are dealing with a circular orbit cone-beam data, the FDK algorithm [43] was employed for preliminary reconstructions which were used in the volume registration step.The registration step was followed by multi-axis PWLS with a fixed number of iterations.In order to accelerate the convergence of PWLS, an ordered subset version of the algorithm was used [44].For all three reconstructions 50 iterations with 25 subsets were carried out.The total number of projections used for the reconstruction was kept constant.The displayed intensity in Figs.3-5 represents an apparent linear attenuation coefficient in mm −1 (e.g. as reconstructed by FDK).

Results and discussion
Figures 4 and 5 illustrate the results of application of different reconstruction algorithms to the simulated and real experimental data.In Fig. 4, one can see a comparison of the standard FDK approach [43], and the FDK-based fusion approach that was recently described in [21] and three algebraic algorithms: Simultaneous Iterative Reconstruction Technique (SIRT), Penalized Weighed Least Squares (PWLS) and an algorithm based on Student's-T distribution (see [29]).Axis of rotation of the additional data used in each reconstruction is indicated by a white arrow (axis of rotation) and a yellow circle (orbit of rotation).Since the simulated data are polychromatic, instead of computing an intensity-based error metric (e.g.mean squared error), the error is computed as a fraction of misclassified voxels after application of segmentation based on three-level Otsu threshold [33].
Our goal is to compare a relatively fast approach of intensity weighted fusion of several FDK reconstructions with slower iterative Least Squares -type methods.We expect the well known SIRT implementation to lead to a solution where errors from all scans are averaged.The PWLS approach should lead to a better solution, giving a lower weight to the measurements with low intensity and, supposedly, high error.The Student's-T approach [29] should lead to a comparable result as it makes outliers less significant.
The first column of the Fig. 4 shows how the error in segmentation decreases from 16% to 5% when more projections are used in the FDK fusion approach.In this method, a separate reconstruction is computed for each scan, after which, a weighted sum of all reconstructions is calculated using weights back-projected into the reconstruction volume.The second to fourth columns show the results of the algebraic methods.In an attempt to provide comparable data, we fixed the number of iterations to 100 for SIRT and the PWLS method, while the Student's-T algorithm was run for 200 iterations since we observed a slower convergence rate for that method.In both methods the values of the initial guess were set to zero.
Additionally, to keep the amount of input data constant, we decimated the number of projections per dataset by a factor of two and three respectively in cases when two and three axes were used.This was not done for the fusion approach, as the reduction of the projection number seriously deteriorates individual FDK reconstructions, making this method inapplicable to undersampled data.It can be seen that PWLS and Student's-T methods provide the best results as the number of datasets used in the reconstruction is increased, even when the total number of projections used in the reconstruction stays constant.Furthermore, if the input datasets are not decimated, due to a lower statistical error, these approaches lead to a further decrease in the segmentation error (1% error for PWLS with compete three-axes data).In the Student's-T approach, the error variance parameter is estimated globally using a misfit function of the projections, as a result, high density regions (metal rods) seem to be treated as outliers during the start of the algorithm and their density is underestimated.This effect is counteracted by using a larger number of iterations.
Reconstruction based on the PWLS approach was tested on several physical objects.Fig. 5 shows slices and projections of a reconstructed volume computed from two scans of a printed circuit board (PCB), an incandescent light bulb, and a nine-volt alkaline battery.The PCB provides a great example of an object where a standard FDK reconstruction fails due to strong attenuation in the plane of the board and around parts that contain metal.The multi-axis PWLS method allows a significant increase in the quality of the reconstruction by combining data from two scans with orthogonal rotation axes.
Reconstruction artefacts are less apparent in the light bulb data, however, in the first scan, the projections are slightly truncated, so part of the object is not reconstructed.In this case, the volume registration algorithm finds a correct orientation of the object even for truncated data.It can be seen that the overall intensity scale of the PWLS-based light bulb reconstruction differs from that of the FDK-based reconstruction, however, that should likely be attributed to the normalization that is applied to the reconstructions.
Scans performed for two different orientations of the nine-volt battery produce the most different results with respect to the FDK approach.For instance, it can be seen that in the second dataset, metal electrodes inside the anode (nail shaped) are not visible due to strong absorption, as well as the space between the anodes and cathodes and the overall intensity of the anodes is heterogeneous (likely due to beam hardening).The multi-axis PWLS reconstruction provides a clearer image of the layer between the anodes and cathodes, the porosity of anodes and reduces the "halo" effect produced by beam hardening, however, the intensity of the inner parts of the metal electrodes seems to be underestimated.This effect may be reduced by further optimization of the weights used in PWLS scheme.However, such optimization lies outside the scope of the While there is no ground-truth data available for these objects, and so the reconstruction error cannot be estimated, the improved quality of reconstructions based on two axes of rotation is apparent.Namely regions of the object that are over-attenuated by metallic parts are more homogeneous, there are less artefacts visible outside of the object (including parts of the mount that are averaged out).These effects are achieved through breaking the symmetry of a single-axis scan and not by introducing any type of explicit beam-hardening correction.In the situation where the bulk of a strongly attenuating object causes a significant beam-hardening effect, the multi-axis approach described in this publication can be combined with one of the beam-hardening approaches [9][10][11][12].
Considering the accuracy of the reconstructed data, FDK fusion may provide a suboptimal result comparing to the PWLS approach (respectively 9% of misclassified pixels comparing to 2% for simulated dual-axis data).However, FDK-based methods should be considered in applications where long computation times represent a critical issue.Reconstruction of the experimental data presented in this manuscript (1450 projections, 972 x 768 pixels each) based on the FDK algorithm takes roughly one minute per dataset on our workstation (Intel Xeon CPU E5-1620 v3, GeForce GTX 970), while computing a PWLS-based reconstruction takes on the order of one to two hours (depending on the parametrization of the PWLS solver).

Conclusion
Our results show that typical tomographic reconstruction artefacts can be dramatically reduced without the use of complex nonlinear models or knowledge about the properties of the object or the scanner.The proposed approach relies on automatic registration of two or more independent tomographic volumes computed from scans done with different orientations of the rotation axis and a subsequent joint reconstruction that uses all available data.The joint reconstruction algorithm can be based on one of the minimization approaches that take into account estimated error variance.In this work, we have explored the performance of the PWLS solver and a solver based on the Student's-T misfit function.An important feature that was observed in our experiments is that the number of projections in each of the scans can be reduced to keep the total radiation dose same as in one full scan, still allowing to reduce reconstruction artefacts.The fact that the proposed approach does not require utilization of advanced rotation stages and gimbal-type sample holders, makes it suitable for a wide variety of applications in laboratorybased X-ray tomography.Further investigation should be pursued in future for optimization of the rotation axis orientation as well as for the reconstruction methods that allow introduction of additional regularization terms or faster convergence.

Fig. 2 .
Fig. 2. Tomographic Reconstruction.(A) Rotation around an object with high density regions: zone adjacent to high-density regions is poorly reconstructed.(B) Addition of 'redundant' data: reconstruction of the over-attenuated zone is now based on more observations

Fig. 4 .
Fig. 4. Comparison of the reconstruction methods: FDK fusion, SIRT, PWLS, Student's-T.Top two rows -single rotation axis, second row -two rotation axes, third row -three rotation axes.The colorbar is given for the reconstructed linear attenuation coefficient in mm −1 .
3(c)-3(f).Corresponding axes of rotation are schematically indicated by a white arrow (axis of rotation) and a yellow circle (orbit of rotation).A reconstruction of the monochromatic data Fig. 3(c) is showing the influence of a limited number of photons (sharp streaks) while the reconstruction of the polychromatic data Figs.3(d)-3(f) is mainly showing the effects of beam hardening.The polychromatic data correspond to different axes of rotation of the object.

Fig. 5 .
Fig. 5. Comparison of reconstruction methods applied to experimental data.Columns (left to right) correspond to FDK reconstructions of data acquired for two rotation axes and a PWLS reconstruction applied to the combined data.The top two rows show slices and maximum projections of a reconstructed PCB.The middle two rows show reconstructions of a light bulb.The bottom two rows show two orthogonal slices of a reconstructed alkaline battery.Colorbars are given for the reconstructed linear attenuation coefficient in mm −1 .