2.1

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

In this work, the estimated image is forward projected using Siddon’s algorithm,⁵ which is a ray-based technique that can be efficiently implemented on the GPU.³ 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.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.⁷

Figure 1.

Flowchart depicting how raw data is processed to perform tomographic reconstructions with OSC-TV.

2.3

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

Figure 2.

Virtual phantom.

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

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

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.

2.4

Scanning protocols

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.

Figure 3.

Samples: (a) water phantom, (b) granite, (c) chalcopyrite, (d) 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.

2.5

Conversion calibration

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

2.6

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

3.1

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.

Figure 4.

Convergence of virtual phantom. Percent error maps.

3.2

Samples

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

Figure 5.

Water phantom: (left) 100 kVp, (right) 140 kVp. Window: [0.18:0.22]

3.3

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.

Figure 6.

BEC A196-6 100 kVp, reconstructed with OSC-TV applying different techniques: (a) backprojection by bilinear interpolation, (b) ray-driven backprojection and region-of-interest strategy, (c) backprojection by bilinear interpolation and region-of-interest strategy. Window [0.5:0.6] cm^‒1.

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.

Figure 7.

Granite 100 kVp: (a) Siemens algorithm (B30s), (b) OSC-TV and the Siddon backprojection, (c) OSC-TV and backprojection by bilinear interpolation, (d) OSC-TV, backprojection by bilinear interpolation and region-of-interest strategy, (e) plot of the line profiles. Window [1.0:1.8]

Figure 8.

Chalcopyrite 140 kVp: (a) Siemens algorithm (B30s), (b) OSC-TV and the Siddon backprojection, (c) OSC-TV and backprojection by bilinear interpolation, (d) OSC-TV, backprojection by bilinear interpolation and region-of-interest strategy, (e) plot of the line profiles. Window [1.0:1.6]