X-ray Computed Tomography (CT) is nowadays ubiquitous in medicine for diagnosis, treatment planning and treatment responses purposes. This technology is also increasingly used for non-medical purposes in many fields, providing several advantages such as: (i) non-destructive testing and (ii) high spatial and density resolution.1
In clinical CT scanners, raw acquisition data are stored as sinograms and are typically processed by proprietary methods, notably to reduce beam hardening artifacts in reconstructed images. Raw sinograms are relatively large from a storage perspective and are typically not kept on the long term, as opposed to reconstructed images which are normally sent to a PACS. Even though it is possible to archive sinograms, their proprietary format typically prevents users from reading them. The sinograms could also be preprocessed (e.g. calibration, beam-hardening correction), and in this respect might not truly represent raw CT data, i.e. measures of the attenuation along a ray. In some cases, including for research purposes, it might be desirable to access this raw attenuation data.
Previous works have shown the importance of the forward projection model on image quality.2 The aim of this work was to explore how the reconstruction of small samples with small voxel size (e.g. 512 pixels × 512 pixels, and a pixel size in x and y of 0.00977 cm), are affected by the backprojection technique implemented on a iterative reconstruction algorithm, OSC-TV (Ordered subset convex algorithm with total variation minimization).3
Besides, a technique used to reconstruct regions-of-interest with this class of algorithms was also analyzed,4 as they allow reducing truncation artifacts caused by objects outside the reconstruction matrix.
MATERIALS AND METHODS
We have developed a framework where an in-house iterative algorithm can be used to reconstruct images based on genuinely raw attenuation data. First, the sinogram data in proprietary format are converted through the use of binaries provided by the manufacturer. The conversion generates usable sinogram data, and also provides associated geometry data. These data provided by the vendor is then used to perform the reconstruction with the iterative algorithm OSC-TV. Two backprojection techniques were evaluated in this work with regards to their impact on image features: a Siddon-based (ray-driven)5 and a bilinear interpolation (voxel-driven) approach.
Projection and backprojection techniques
In the OSC-TV algorithm, the estimated image is forward-projected and backprojected several times, depending on the number of iterations and subsets.3
In this work, the estimated image is forward projected using Siddon’s algorithm,5 which is a ray-based technique that can be efficiently implemented on the GPU.3 Backprojection was performed in this work with two techniques: (i) ray-based or (ii) voxel-driven by bilinear interpolation.
In the voxel-driven backprojection by bilinear interpolation (BLI), the center of each voxel is projected on the detector and the reading corresponds to the weighted sum of the four neighboring pixels (the position is calculated relative to the detector pixel centers).
On the other hand, for a ray-driven backprojection, the detector readings are smeared back across the image. The voxels that are incremented in this process depend on the ray path, from the detector reading to the source. Typically, a finite number of rays are defined (e.g. one ray per detector). In this approach, depending on the angle between rays and the voxel size, some voxels might not be traversed by a given ray, or they are under-utilized, for the intersection length is negligible. This problem is not unique to the backprojection, but it is also present in the forward-projection in different techniques (e.g. Siddon, Joseph’s method and bilinear interpolation).2
Working with proprietary format
For this work, a Siemens SOMATOM Definition AS+ 128 scanner was used. This device is installed at the Institut national de la recherche scientifique, in Québec City, Canada. This platform is used for several non-medical applications, including material characterization and custom beam hardening corrections with dual-energy techniques.6, 7
The process of using raw data from this medical CT scanner is illustrated in Fig. 1. In summary, a vendor-provided calibration table is used for acquisitions to cancel any vendor-specific beam hardening corrections (detector calibration is still applied). The raw data is stored in the host system and copied to a different machine for archiving purposes, along with calibration data. These raw data, free of vendor-specific beam hardening corrections, can be read with vendor-provided binaries. These genuinely raw data (except for detector calibration), can thereafter be used in custom reconstruction algorithms designed to handle corrections from first principles, notably through dual-energy approaches.7
Convergence of OSC-TV in numerical simulations
In order to verify the convergence of the reconstructed image, a virtual phantom was defined (see Fig. 2). It is composed of a water cylinder with 25.2 cm of diameter and 9 cylinder rods of distinct materials, each one with a diameter of 2.4 cm. The X-ray absorption properties of the phantom, defined by the linear attenuation coefficient, was retrieved from the NIST XCOM database.8
A noise-free monoenergetic step-and-shoot acquisition of the virtual phantom at 83 keV was simulated by using the geometry of the medical CT scanner Siemens SOMATOM Definition AS+ 128.
The convergence of the reconstruction to an optimal result is tuned by the reconstruction parameters: number of iterations, number of subsets, final number of subsets, and initial image; regularization parameters also play an important role in decreasing the overall noise and controlling the spatial resolution: gradient steps and strength of regularization (rms).3
A high number of subsets, equivalent to half the number of projections and a final number of subsets representing 1/10 of that value, can be used to achieve optimal convergence, as already suggested with few-view acquisitions.9
Different combinations of reconstruction parameters are used to reconstruct the virtual phantom in order to assess the convergence: 5, 9 and 12 iterations; 84, 576 (1/4 projections) and 1152 (1/2 projections) subsets. A total of 9 reconstructions are performed (3×3), with the regularization constant and gradient steps fixed at 0.02 and 20, respectively.
Four samples were imaged in the experimental protocol (see Fig. 3): (a) water phantom, (b,c) small mineral sample (approximately 5 cm) of granite and chalcopyrite, respectively, and (d) a sandstone, with a diameter of 10.0 cm, henceforth called BEC A196-6.
Different scan protocols were used for the samples, and for some cases proprietary beam-hardening correction (BHC) was neutralized (neutral): water phantom and chalcopyrite (neutral) at 140 kVp, BEC A196-6 (neutral) and granite (neutral) at 100 kVp. Tomographic reconstructions performed with the Siemens algorithm, filtered backprojection, used the B30s (smooth) kernel.
Tomographic reconstruction performed by the OSC-TV algorithm are inherently in units of linear attenuation (cm‒1), while the ones obtained with the Siemens algorithm are in Hounsfield units (HU). The relation between these quantities is given by: . Reconstruction of the water phantom, at 100 kVp and 140 kVp, are performed with OSC-TV, and the mean value in a ROI of each case is calculated. The average of the results provides μwater. This value is independent of the tube voltage, for the projections are always normalized for water, so HU is always close to 0 (proprietary preprocessing of raw data).
Reconstruction of raw data using OSC-TV
Following the workflow depicted in Fig. 1, and the numerically-determined optimal reconstruction parameters defined in Section 2.5, images were reconstructed with well-defined parameters of the OSC-TV code.
For the cases where small voxels are used with a regular grid (512 pixels × 512 pixels), a modified iterative reconstruction of a region-of-interest (IR ROI) technique based on the work of Ziegler et al. 20084 was applied. The CT table lies outside the image reconstruction matrix for some cases. This technique allows us to perform reconstruction with a normal grid and small pixel size (e.g. 0.00977 cm), avoiding truncation artifacts caused by objects outside the reconstruction matrix.
Contrariwise, high-resolution reconstructions (e.g. 2048 pixels × 2048 pixels) would have to be made to cover the entire field-of-view (FOV) (e.g. 50 cm), for the table presents an important attenuation, even though a small object (few centimeters) is being reconstructed.
Variations of the OSC-TV algorithm are identified by acronyms, where BLI stands for backprojection by bilinear interpolation (voxel-driven), Siddon for ray-driven backprojection, and ROI for the IR ROI technique
RESULTS AND DISCUSSION
Convergence of OSC-TV in numerical simulations
The visual convergence of the OSC-TV algorithm in terms of number of iterations and subsets is depicted in Fig. 4. As one can notice, 84 subsets is insufficient even when 12 iterations are performed, giving rise to beam-hardening-like artifacts. As those images were reconstructed from monochromatic projections (83 keV), such artifacts were not expected, and so are due to a non-optimal convergence. By increasing the number of subsets (576, or 1/4 of the number of projections), and the number of iterations, such artifacts are decreased. Finally, it is only when a high number of subsets (1152, which is equivalent to half the number of projections) is used that those artifacts are removed, as a result of high rate convergence.
From the measurements in the central slice of the water phantom, at 100 kVp and 140 kVp, it was obtained an average value μwater = 0.1918 cm‒1, which allows conversion from HU to linear attenuation coefficient to be performed (see Fig. 5).
Reconstruction of ROI and small samples: OSC-TV ROI
The importance of the backprojection technique for artifact mitigation and the IR ROI technique for removing truncation artifacts is shown in Fig 6. In (a), the CT table was not included in the reconstruction matrix. It is worth noting that the setup uses a custom table, with much more important absorption properties than a medical one. Secondly, in (b), the image is reconstructed using Siddon backprojection and IR ROI technique, so the table is removed from the projection data. As not all voxels are incremented by the correspondent detector read during backprojection, more artifacts arise, producing a noisy image. Finally, in (c), the IR ROI and the voxel-driven backprojection by bilinear interpolation are combined. These two techniques are capable of mitigating both artifacts: truncation and lack of data in voxel increment during backprojection.
In Figs. 7 and 8, the reconstruction of a small sample of granite and chalcopyrite (approximately 5 cm) is shown with different techniques: (a) Siemens with the B30s kernel (smooth); (b) OSC-TV with Siddon’s backprojection; (c) OSC-TV with backprojection by bilinear interpolation; (d) OSC-TV with IR ROI technique and backprojection by bilinear interpolation. As a smooth kernel was selected for the Siemens reconstruction, its natural its OSC-TV counterpart (d) is less blurry. Image (b) suffers from the inherent problem of the Siddon’s backprojection, where some voxels are not properly incremented by the backprojection and the artifacts are scattered from the center of the image, and the presence of noise. When (c) and (d) are compared, it is clear that the first is smoother. Contrariwise, the edges are smoother in general with the OSC-TV technique. The iterative technique also provides images with less beam-hardening artifacts: at the left edge of the granite in Fig 7, pixels values are steady, while a smooth decrease is observed in the Siemens reconstruction.
In this work, a framework for working with proprietary data from a commercial CT scanner was presented. Binaries provided by the manufacturer allows research to be performed in both the projection and the image space, where beam-hardening preprocessing can also be neutralized, so custom BHC can also be applied. Ray-driven backprojection, with a limited number of rays, is insufficient for small geometries and samples. With a backprojection by bilinear interpolation (voxel-driven), results with less noise and free of reconstruction artifacts were obtained for studied cases. Strategies to reconstruct regions-of-interest were also applied, where truncation artifacts were reduced, and less blurring was observed.
This research was supported by the Sentinel North program of Université Laval, made possible by the Canada First Research Excellence Fund. This work was also supported by a grant from Fonds de recherche du Québec - Nature et technologies [2018-PR-206076]. Karl Stierstorfer (software support), Philippe Letellier and Mathieu des Roches (support in scanning the samples), Stéphanie Larmagnat (providing samples), Jérôme Landry and Pascal Bourgault (algorithm development).
Di Schiavi Trotta, L., Matenine, D., Martini, M., Stierstorfer, K., Lemaréchal, Y., Francus, P., and Després, P., “Beam-hardening corrections through a polychromatic projection model integrated to an iterative reconstruction algorithm,” NDT & E International, 102594 (2021). Google Scholar
Suplee, C., “XCOM: Photon Cross Sections Database,” (2009). Google Scholar