Suppression of Cone-Beam Artefacts with Direct Iterative Reconstruction Computed Tomography Trajectories (DIRECTT)

The reconstruction of cone-beam computed tomography data using filtered back-projection algorithms unavoidably results in severe artefacts. We describe how the Direct Iterative Reconstruction of Computed Tomography Trajectories (DIRECTT) algorithm can be combined with a model of the artefacts for the reconstruction of such data. The implementation of DIRECTT results in reconstructed volumes of superior quality compared to the conventional algorithms.


Introduction
Computed tomography (CT) systems for non-destructive testing and material analysis generally use a cone beam on a sample that rotates in a circular orbit [1], with cylindrical samples being among the most common [1][2][3][4]. The exact reconstruction of data acquired during such a measurement is not possible because the geometry does not satisfy Tuy's sufficiency condition [5]. This is demonstrated in Figure 1 with the reconstruction of a concrete rod by the commonly used algorithm developed by Feldkamp, Davis and Kress (FDK) [6]. For higher cone angles, there is a decrease of the grey values, which represent the attenuation coefficient µ, and of the image quality in the direction of the rotation axis (z-axis in Figure 1).
Several algorithms have been proposed to reduce such artefacts. Hsieh proposed a two-pass algorithm that estimates the cone-beam artefacts from the segmented highdensity material and then subtracts them from the FDK reconstruction [7]. Han and Baek went further by devising a multi-pass approach that they tested for larger cone angles and different material densities [8]. Maaß et al. proposed an iterative algorithm that also subtracts the estimated artefacts from the FDK reconstruction without requiring segmentation [9].
Here, we will describe how we have adjusted the Direct Iterative Reconstruction of Computed Tomography Trajectories (DIRECTT) algorithm [10][11][12] to estimate such artefacts and compensate for them.

Sample Images
A set of 3000 cone-beam projections of the concrete rod of Figure 1 was acquired over 360° on an in-house GE v|tome|x L 300 scanner. A 2024 × 2024 PerkinElmer detector with a pixel size of 0.2 mm was used. Source-object and source-detector distances of 81 mm and 1018 mm, respectively, resulted in a magnification of 12.5 for a voxel size of 0.016 mm. The voltage and current settings of the source were set to 140 kV and 80 μA, respectively. A 0.5 mm Cu prefilter was used. The acquisition time per projection was 6 s.
The geometry of the CT scan of the cylindrical sample is represented, not to scale, by Figure 2. The orange cone represents the field of view (FoV), while the blue dashed lines represent rays that traverse the front and rear edges of the sample. Near the lower edge of the FoV, an inverse conical area of the sample is defined by the solid and dashed orange lines. This is the part of the sample that lies within the FoV during only some of the projections and, therefore, is not fully reconstructible by FDK [13]. We consider the case that dimensions of the sample are not known precisely. It is possible to determine them from the projections. The total height h of the sample within the FoV is the sum of its parts h1 and h2 that extend respectively above and below the plane

Sample Images
A set of 3000 cone-beam projections of the concrete rod of Figure 1 was acquired over 360 • on an in-house GE v|tome|x L 300 scanner. A 2024 × 2024 PerkinElmer detector with a pixel size of 0.2 mm was used. Source-object and source-detector distances of 81 mm and 1018 mm, respectively, resulted in a magnification of 12.5 for a voxel size of 0.016 mm. The voltage and current settings of the source were set to 140 kV and 80 µA, respectively. A 0.5 mm Cu prefilter was used. The acquisition time per projection was 6 s.
The geometry of the CT scan of the cylindrical sample is represented, not to scale, by Figure 2. The orange cone represents the field of view (FoV), while the blue dashed lines represent rays that traverse the front and rear edges of the sample. Near the lower edge of the FoV, an inverse conical area of the sample is defined by the solid and dashed orange lines. This is the part of the sample that lies within the FoV during only some of the projections and, therefore, is not fully reconstructible by FDK [13].

Sample Images
A set of 3000 cone-beam projections of the concrete rod of Figure 1 was acquired over 360° on an in-house GE v|tome|x L 300 scanner. A 2024 × 2024 PerkinElmer detector with a pixel size of 0.2 mm was used. Source-object and source-detector distances of 81 mm and 1018 mm, respectively, resulted in a magnification of 12.5 for a voxel size of 0.016 mm. The voltage and current settings of the source were set to 140 kV and 80 μA, respectively. A 0.5 mm Cu prefilter was used. The acquisition time per projection was 6 s.
The geometry of the CT scan of the cylindrical sample is represented, not to scale, by Figure 2. The orange cone represents the field of view (FoV), while the blue dashed lines represent rays that traverse the front and rear edges of the sample. Near the lower edge of the FoV, an inverse conical area of the sample is defined by the solid and dashed orange lines. This is the part of the sample that lies within the FoV during only some of the projections and, therefore, is not fully reconstructible by FDK [13]. We consider the case that dimensions of the sample are not known precisely. It is possible to determine them from the projections. The total height h of the sample within the FoV is the sum of its parts h1 and h2 that extend respectively above and below the plane We consider the case that dimensions of the sample are not known precisely. It is possible to determine them from the projections. The total height h of the sample within the FoV is the sum of its parts h 1 and h 2 that extend respectively above and below the plane SOD. The plane is defined by the source (S), the centre of rotation (O) and the centre of the detector (D). The plane SOD is assumed to be perpendicular to the detector.
We can see from Figure 2 that, and, where r is the radius of the sample, and DF and DR the distance between the central detector row and the detector rows where the front (F) and rear (R) edge of the concrete rod are respectively projected. Dividing Equation (1) by Equation (2) and rearranging, we obtain: The length h 1 can be calculated now from either Equation (1) or Equation (2), while the length h 2 is: where H is the height of the detector. The part of the sample that, as mentioned above, does not always lie within the FoV has been accounted for through the inclusion of the radius r in Equation (4).

The Direct Iterative Reconstruction of Computed Tomography Trajectories (DIRECTT) Algorithm
The DIRECTT algorithm was first proposed for the reconstruction of two-dimensional (2D) images by Lange et al. and, in a previous article [12], we introduced a new, more efficient, and fully 3D version. The algorithm operates on finding the best solution possible by mimicking the actual physical projection process, instead of directly solving the inverse problem. It only reconstructs certain voxels during each iteration, simulating the projection of the partial reconstruction, and repeating the workflow for the difference between measured and simulated projections until this difference is sufficiently close to zero [10][11][12]. Although the concrete rod of Figure 1 was one of the two datasets that were used to showcase the performance of DIRECTT, the algorithm was implemented only on the slice that corresponds to the cross section of the sample with the plane SOD [12]. That slice will be hereafter referred to as the central slice. Attempting to implement DIRECTT on the whole dataset does not lead to an improvement over the FDK reconstruction of Figure 1. On the contrary, it results in severe artefacts and missing data because the algorithm fails to predict the decreasing grey values along the z-axis. However, these artefacts can be reduced by modelling them based on the shape of the sample.

Software
For this work, the forward-and back-projection operations involved in DIRECTT were performed using the Python programming language and the open-source ASTRA (All Scale Tomographic Reconstruction Antwerp) toolbox [14]. Via ASTRA, computationally demanding operations are offloaded to a graphics processing unit using the CUDA (Compute Unified Device Architecture) language. The toolbox also includes several reconstruction algorithms, such as the FDK, the simultaneous iterative reconstruction technique (SIRT) [15], and a conjugate gradient (CG) method based on the Krylov subspace [16], that can run with little input from the user [14].

Results
A single measured projection of the concrete and the corresponding simulated projection of a virtual homogeneous cylinder of height h and radius r are shown in Figure 3a,b, respectively. The two horizontal lines in the former indicate the detector rows that correspond to points F and R of Figure 2. The two rows are identified automatically from the projections. Specifically, F is the lowermost detector row containing exclusively values that correspond to the background. Similarly, R is the lowermost row containing exclusively values lower than the mode of all absorption projections, which roughly corresponds to the absorption of rays that penetrate the sample perpendicularly.
correspond to points F and R of Figure 2. The two rows are identified automatically fr the projections. Specifically, F is the lowermost detector row containing exclusively valu that correspond to the background. Similarly, R is the lowermost row containing exc sively values lower than the mode of all absorption projections, which roughly cor sponds to the absorption of rays that penetrate the sample perpendicularly.
By simulating the geometry of the scan, projecting and back-projecting the virt cylinder, and normalizing for the central slice, a 3D model M of the artefacts that ar during the reconstruction is computed. Normalizing for the central slice ensures t hardly any artefacts are considered for the parts of the volume that always lie within FoV and can be accurately reconstructed by FDK. Note that the model M needs to be co puted just once. The slices of M that correspond to those in Figure 1 are shown in the t row of Figure 4. The volume resulting from the back-projection of the unfiltered proj tions of the concrete rod, the first step during each iteration of DIRECTT [12], is shown the bottom row of Figure 4. Both volumes of Figure 4 comprise more slices along the axis than that of Figure 1. The extra slices have been included because it is essential the implementation of any iterative reconstruction algorithm, such as DIRECTT, that parts of the sample that do not lie within the FoV for every projection are reconstruc too.  By simulating the geometry of the scan, projecting and back-projecting the virtual cylinder, and normalizing for the central slice, a 3D model M of the artefacts that arise during the reconstruction is computed. Normalizing for the central slice ensures that hardly any artefacts are considered for the parts of the volume that always lie within the FoV and can be accurately reconstructed by FDK. Note that the model M needs to be computed just once. The slices of M that correspond to those in Figure 1 are shown in the top row of Figure 4. The volume resulting from the back-projection of the unfiltered projections of the concrete rod, the first step during each iteration of DIRECTT [12], is shown in the bottom row of Figure 4. Both volumes of Figure 4 comprise more slices along the z-axis than that of Figure 1. The extra slices have been included because it is essential for the implementation of any iterative reconstruction algorithm, such as DIRECTT, that the parts of the sample that do not lie within the FoV for every projection are reconstructed too.
There is an obvious qualitative relation between the two rows of Figure 4. The pixelby-pixel division of the reconstructed volume by the model M results in a "corrected" volume, on which DIRECTT can be successfully implemented. The threshold values, based on which the voxels to be reconstructed during each iteration are selected, are calculated from the central slice as described in [12] but are applied simultaneously on the whole volume. The algorithm terminates when the projections of the reconstruction array match the measured ones. The final reconstructed volume is shown in Figure 5. There is an obvious qualitative relation between the two rows of Figure 4. The pixelby-pixel division of the reconstructed volume by the model M results in a "corrected" volume, on which DIRECTT can be successfully implemented. The threshold values, based on which the voxels to be reconstructed during each iteration are selected, are calculated from the central slice as described in [12] but are applied simultaneously on the whole volume. The algorithm terminates when the projections of the reconstruction array match the measured ones. The final reconstructed volume is shown in Figure 5.

Discussion and Conclusions
The reconstruction by DIRECTT is an improvement over that by FDK in Figure 1 as there are no artefacts linked to the increasing cone angle. Lacking the ground truth for the reconstructed volume, the evaluation of the two algorithms using a full-reference metric is meaningless. Nevertheless, a quantitative evaluation of the algorithms can be  There is an obvious qualitative relation between the two rows of Figure 4. The pixelby-pixel division of the reconstructed volume by the model M results in a "corrected" volume, on which DIRECTT can be successfully implemented. The threshold values, based on which the voxels to be reconstructed during each iteration are selected, are calculated from the central slice as described in [12] but are applied simultaneously on the whole volume. The algorithm terminates when the projections of the reconstruction array match the measured ones. The final reconstructed volume is shown in Figure 5.

Discussion and Conclusions
The reconstruction by DIRECTT is an improvement over that by FDK in Figure 1 as there are no artefacts linked to the increasing cone angle. Lacking the ground truth for the reconstructed volume, the evaluation of the two algorithms using a full-reference metric is meaningless. Nevertheless, a quantitative evaluation of the algorithms can be

Discussion and Conclusions
The reconstruction by DIRECTT is an improvement over that by FDK in Figure 1 as there are no artefacts linked to the increasing cone angle. Lacking the ground truth for the reconstructed volume, the evaluation of the two algorithms using a full-reference metric is meaningless. Nevertheless, a quantitative evaluation of the algorithms can be undertaken by calculating the histogram entropy of the respective volumes. The histogram entropy (HE) is a global metric defined according to the relation: where p(µ) is the distribution function of the grey values. The value of the HE increases if homogeneously distributed noise is present in the image, and decreases if sharp edges are present [17]. Therefore, a low value is an indication of a good balance between noise and blur. The FDK-and DIRECTT-reconstructed volumes (Figures 1 and 5, respectively) have histogram entropies of 2.86 and 1.89, respectively. An alternative way to evaluate the results of the algorithms is to simulate the projection of the reconstructed volumes and calculate the Pearson correlation coefficient [12,18] between these projections and the measured ones. The Pearson correlation coefficient (PCC), the value of which can range between 1 for total linear correlation and −1 for total linear anti-correlation, is calculated according to the relation: where σ M,P is the covariance and σ M and σ P are the standard deviation of the measured and simulated projections, respectively. Such coefficients are calculated for each detector pixel and are plotted in Figure 6. While the PCC values decrease for large cone angles in the case of FDK, they remain near 1 regardless of the angle in the case of DIRECTT.
J. Imaging 2021, 7, x FOR PEER REVIEW 6 o undertaken by calculating the histogram entropy of the respective volumes. The his gram entropy (HE) is a global metric defined according to the relation: where p(μ) is the distribution function of the grey values. The value of the HE increase homogeneously distributed noise is present in the image, and decreases if sharp edges present [17]. Therefore, a low value is an indication of a good balance between noise a blur. The FDK-and DIRECTT-reconstructed volumes (Figures 1 and 5, respectively) ha histogram entropies of 2.86 and 1.89, respectively. An alternative way to evaluate the results of the algorithms is to simulate the proj tion of the reconstructed volumes and calculate the Pearson correlation coefficient [12, between these projections and the measured ones. The Pearson correlation coeffici (PCC), the value of which can range between 1 for total linear correlation and −1 for to linear anti-correlation, is calculated according to the relation: where σM,P is the covariance and σM and σP are the standard deviation of the measu and simulated projections, respectively. Such coefficients are calculated for each detec pixel and are plotted in Figure 6. While the PCC values decrease for large cone angles the case of FDK, they remain near 1 regardless of the angle in the case of DIRECTT. The comparison between the two algorithms can also be done locally on parts of volume that are affected in a greater degree by the cone-beam artefacts. For instance, profiles along the x-axis through the centres of the xy-slices in both Figures 1 and 5 plotted in Figure 7. It is evident that, in the case of FDK, the attenuation values away fr the centre of the sample are underestimated. The comparison between the two algorithms can also be done locally on parts of the volume that are affected in a greater degree by the cone-beam artefacts. For instance, the profiles along the x-axis through the centres of the xy-slices in both Figures 1 and 5 are plotted in Figure 7. It is evident that, in the case of FDK, the attenuation values away from the centre of the sample are underestimated.
Finally, the reconstruction of the volume has also been performed using SIRT and CG, which were programmed to perform a fixed number of 600 iterations. The results are shown in Figure 8. The values of the metrics for each reconstruction algorithm are listed in Table 1. While the simulated projections in the case of both SIRT and CG have a high correlation to the measured ones, the respective histogram entropy values are significantly higher than that of DIRECTT. Moreover, although CG has, for the better part, suppressed the artefacts that are linked to the increasing cone angle, it has, at the same time, resulted in ring-like artefacts exactly where the regions that lie fully within the FoV border the regions that are only partially within it. Finally, the reconstruction of the volume has also been performed using SIRT and CG, which were programmed to perform a fixed number of 600 iterations. The results are shown in Figure 8. The values of the metrics for each reconstruction algorithm are listed in Table 1. While the simulated projections in the case of both SIRT and CG have a high correlation to the measured ones, the respective histogram entropy values are significantly higher than that of DIRECTT. Moreover, although CG has, for the better part, suppressed the artefacts that are linked to the increasing cone angle, it has, at the same time, resulted in ring-like artefacts exactly where the regions that lie fully within the FoV border the regions that are only partially within it.
To sum up, we have described how the DIRECTT algorithm is adjusted for the reconstruction of cone-beam CT data. The artefacts that normally arise during such an operation have been suppressed resulting in a reconstruction that is a clear improvement over that computed by FDK. The performance of DIRECTT has been quantified by both global and local metrics and compared to the performance of other iterative reconstruction algorithms.  Figure 1 through omission of its upper and lower slices. The position of the slices shown corresponds precisely to that of the slices in Figure 1.  Figure 8. Orthogonal slices through the volume of a concrete rod as reconstructed by simultaneous iterative reconstruction technique (SIRT, top row) and conjugate gradient (CG, bottom row). The volume has the same size as the volume shown in Figure 1 through omission of its upper and lower slices. The position of the slices shown corresponds precisely to that of the slices in Figure 1.
To sum up, we have described how the DIRECTT algorithm is adjusted for the reconstruction of cone-beam CT data. The artefacts that normally arise during such an operation have been suppressed resulting in a reconstruction that is a clear improvement J. Imaging 2021, 7, 147 8 of 9 over that computed by FDK. The performance of DIRECTT has been quantified by both global and local metrics and compared to the performance of other iterative reconstruction algorithms.  Informed Consent Statement: Not applicable.

Data Availability Statement:
Restrictions apply to the availability of these data. Data were obtained from BAM and are available from the authors with the permission of BAM.