Extending the dynamic range of biomedical micro-computed tomography for application to geomaterials

2.1 Theoretical basis of the correction method

The Beer-Lambert law, which relates the absorbance of photons and a material’s attenuation properties, states: I = I 0 e - μ x , where I is the measured intensity at the detector, I _o is the incident intensity, μ is the energy-dependent linear attenuation coefficient of the material, and x is the path length traveled by the x rays though the material. Based on this relationship, the absorbance A can be easily calculated as:A = ln ( I 0 / I ) = μ x .(1)

From Equation 1, the relationship between absorbance and material thickness is linear for mono-energetic x-ray beams, in the absence of scatter. Deviations from this relationship arise with poly-energetic x-ray beams (as are typically used in micro-CT scanners). When poly-energetic beams are used, the low-energy photons are absorbed preferentially, resulting in a “harder beam” (i.e. the mean energy of the transmitted x-ray spectrum is higher than that of the incident x-ray spectrum). As a result, relatively more photons will pass through the sample than if a mono-energetic x-ray beam (of an effective mean energy and effective flux) was used. In addition, scattering of x-rays within the material also results in an increase of detected photons. By “mimicking” decreased absorbance, both beam hardening and scatter result in an underestimation of the material’s attenuation coefficient, particularly for longer path lengths through the material. The artifacts following the CT reconstruction process, appear as “cupping” in the CT image (i.e. decrease of measured attenuation coefficient in the center of the image, compared to the edges of the object). Figure 1 is a representative plot of the measured absorbance as a function of material thickness, demonstrating the expected non-linear relationship.

The artifacts resulting from beam hardening and scatter can be reduced by implementing a correction that is based on a curve fit of empirically acquired absorbance data of objects with known materials and thicknesses. Previous authors have used point-to-point linear [10] and polynomial [8] fits to the empirically derived absorbance data; Herman’s original beam hardening correction showed negligible improvement in artifact reduction when increasing the polynomial order above two (i.e. A = ax ² + bx, where a and b are constants). For this study, and based on the non-linear shape observed in Fig. 1, we also considered a functional form that increases monotonically to a plateau value, i.e. inverse exponential decay as described below:

OPEN IN VIEWER

Fig.1

Idealized plot of absorbance vs. material thickness to demonstrate the effect of beam hardening and scatter. The straight line represents the response expected from Equation 1. This non-linear curve (dashed line) represents the measured absorbance.

A = P ( 1 - e ( Kx ) ) .(2) where in Equation 2, P is the absorbance at infinite thickness (the plateau) and K is a rate constant.

The fitted functions can then be used to linearize the data by rearranging the formulae to solve for the thickness from the recorded absorbance. This effectively remaps the observed absorbance data to linearized (i.e. corrected) values on a pixel-by-pixel basis, for each projection view. The linearization of a second-order polynomial results in the equation:x = - b + b 2 + 4 aA 2 a(3) while the linearization of the inverse exponential decay results in the equation:x = ln ( P - A ) P / - K .(4)

For the empirical determination of absorbance vs. material thickness curves, two calibration phantoms (i.e. calibration objects) were designed and built. Each phantom was fabricated of a single material and covered a range of thicknesses appropriate for the range of anticipated applications. The phantoms comprised eight coaxial cylindrical plates with known diameters ranging from 4.5 mm to 60 mm arranged in in an alternating pattern to mimic the general shape of geomaterial specimens, as shown in Fig. 2. One calibration phantom (Fig. 2a) was manufactured using poly[methylmethacrylate] (PMMA, Lucite), which mimics soft, water-like materials typically used in biomedical applications [10]; the PMMA plate thickness was 4.8 mm. The second calibration phantom (Fig. 2b) was manufactured using aluminum (6061-T6) to mimic dense objects like bone, silicate minerals, and other geomaterials with similar electron densities; the Al plate thickness was 6.4 mm. To ensure larger areas of uniform thickness on each plate, two opposing sides were flattened prior to imaging. Using these calibration phantoms, a single x-ray projection image yields the eight absorbance measurements (Fig. 2c, d) required to determine characterize the relationship between absorbance and material thickness.

OPEN IN VIEWER

Fig.2

Calibration phantoms for beam hardening correction in CT. Phantoms were fabricated from a) PMMA and b) aluminum. Note the flattened areas on each plate, machined to provide a larger measurement area. The scale cube (1 cm) provides a size reference. Representative absorbance images of the PMMA (c) and aluminum (d) phantoms are also shown.

Two biomedical laboratory micro-CT scanners were used in this study to validate the versatility of the methods. Both scanners utilize cone-beam geometry, where the x-ray source and detector rotate around the specimen.

The GE eXplore Ultra (GE Healthcare, London Canada) is a large-bore (20 cm) slip-ring gantry scanner designed for rapid and fast throughput imaging. The operation and characterization of the scanner have been described in detail previously [11]. Briefly, 1000 projection views, of matrix size 1024×680, are obtained over a single 16-second scan covering 360°. Corrections for bad detector elements, pixel-to-pixel non-uniformity, and dark current were performed prior to cone-beam CT reconstruction. The resulting three-dimensional (3D) images have a transaxial field of view of 14 cm and axial extent of 10.4 cm and an isotropic voxel spacing of 154μm [11]. For all experiments in this study the x-ray beam parameters used on the Ultra scanner were 120 kVp and 20 mA.

The second scanner was a GE eXplore speCZT (GE Healthcare, London Canada) scanner, which is a step-and-shoot device with an 89 mm bore size. For each scan, 900 views were obtained in 0.4° increments, with an exposure interval of 16 ms per frame; the total scan time was 5 minutes. The reconstructed CT images have a transaxial field of view of 70 mm and cover an axial extent of 5.4 cm; the voxel spacing on this scanner, with this acquisition protocol, is 50μm, isotropic. For all experiments in this study the x-ray beam parameters used on the speCZT scanner were 90 kVp and current of 40 mA. All data was reconstructed into signed 16-bit data; this extends the dynamic range of the CT data in Hounsfield units to a maximum of 32,767 HU.

2.4 Acquisition and analysis of calibration data

For each scanner, calibration was performed using both the aluminum and PMMA phantoms. Each phantom was placed on the scanner bed and the source-detector system was oriented perpendicular to the flat edges of the calibration phantom. A single x-ray image of each phantom was acquired and was sufficient to generate measurements of absorbance for all thicknesses. The images were corrected for detector-element non-uniformity, then converted to absorbance images (Fig. 2c, d) prior to analysis. ImageJ software (National Institutes of Health, Bethesda, MD) was utilized for analyzing the absorbance images. Regions of interest (ROI) – 4×4 pixels – were placed on each of the flat regions of each plate (representing a different thickness) of each phantom and the mean and standard deviation of the absorbance were recorded.

The absorbance-thickness relationships were fitted using non-linear Levenberg-Marquardt fitting in Prism (GraphPad Software, Inc., La Jolla, CA). As described above, both a second-order polynomial [8, 12] and the one-phase exponential decay function (Equation 2) were fitted. Calibration fits were obtained using both fitting functions for each scanner and both the aluminum and PMMA phantoms.

2.5 Implementation of beam hardening correction

To incorporate the beam hardening (and scatter) correction prior to reconstruction, custom software was written to linearize each absorbance value on a pixel-by-pixel basis for each projection image acquired. Linearization was performed using Equations 3 and 4 for the polynomial and inverse exponential decay fitted functions. It is important to note that the linearization process results in values of absorbance that are floating point numbers ranging up to approximately 5, for geomaterials. While some reconstruction software can accept the floating-point absorbances (e.g., the reconstruction software used with the Ultra scanner data) others (e.g. the speCZT scanner) requires absorbance data be saved as 16-bit integer values. In the latter, the linearized absorbances must be rescaled prior to reconstruction, while taking care that the scaled values do not exceed the maximum 16-bit integer range.

2.6 Evaluation of beam hardening correction using single-material test objects

2.6.2 Analysis of reconstructed images (cupping)

Reconstructed images were analyzed using ImageJ. For each test-rod image, nine 4×4 voxel ROIs were drawn – one at the center, four at the edges of the rod, and four in the air surrounding the rods, as shown in Fig. 3. From these measurements, the extent of cupping was calculated as follows:cupping = 1 - ( C - B ) E - B(5)

OPEN IN VIEWER

Fig.3

Locations of the nine regions of interest used to calculate percent cupping for the test rods. The image shown in this image is from the 63.5 mm PMMA cylinder scanned on the speCZT scanner.

In Equation 5, E, C and B are the average voxel values in the edge, center and background air, respectively. Averages and standard deviations of cupping (percent) were calculated for five slices of each reconstructed image of a PMMA or Al test rod. Additionally, lines were drawn through the center of each image and the image intensities were plotted to enable visual evaluation of the effectiveness of the corrections.

The ultimate application of the proposed semi-empirical correction technique is to improve the uniformity and accuracy of scans of geomaterials, such as meteorites. As a proof-of-principle example, we have applied the derived calibration curves to correct a representative geomaterial sample, in this case a 2.5 cm³ sample of the fresh Bassikounou meteorite (H5 chondrite), which fell on Oct 16, 2006 [13]. The specimen was scanned using the speCZT scanner, with the same acquisition parameters as all other test objects (90 kVp and 40 mA). Absorbance data were corrected using the poly-Al calibration and 3D images were reconstructed using both corrected and un-corrected. Line profiles were drawn through the reconstructed images to demonstrate the reduction in cupping artifact. Furthermore, volumetric analysis was performed to quantify the metallic fraction (by volume) of the meteorite.

3.1 Calibration

Both non-linear fit functions fit the experimental data acquired from both scanners and both calibration phantoms with an R² >0.99; details are provided in Fig. 4 and Table 1. As expected, both scanners yielded similar results, but different fit parameter values. The differences are attributed to the fact that the images were acquired with different x-ray energies and using different detectors. In addition, the scanners store absorbance values in different data formats and scale the data differently.

OPEN IN VIEWER

Fig.4

Plots of absorbance vs. thickness for the aluminum and PMMA calibration phantoms scanned in the Ultra (a) and specCZT (b) scanners. The symbols represent the measured data (error bars are plotted, but in most cases are smaller than the symbol size). The lines represent the polynomial (poly) and inverse exponential (exp) fits as indicated in the legend.

4.1 Calibration phantoms

In this study we have shown that biomedical micro-CT scanners can be used to analyze geomaterial specimens following correction for beam hardening and scatter. The correction step requires multiple thicknesses of a material, with x-ray absorption properties similar to those of the samples to be scanned, to be imaged with a fixed protocol. While an estimate of the change in absorbance with thickness can be obtained by sequentially imaging samples of different thicknesses, the phantoms presented in this study enable calibration to be performed in a single step. Furthermore, the shape of the phantoms – involving plates of varying diameter arranged in an alternating manner on either side of the largest plate – mimics not only the beam hardening effect, but also mimics the effect of scattered photons. It is important to note that the calibrators provide estimates of the two effects and cannot mimic either perfectly for a range of applications.

In this study we used an aluminum calibration phantom to estimate the parameters for the beam-hardening correction methods. We expect that for tailored applications, phantoms can be made of different materials – e.g. a concrete phantom of similar shape to those shown in Fig. 2 could be used to calibrate for scans of concrete cylinders; phantoms could also be molded of clay or machined out of granite for further examination of those materials.

4.2 Correction for cupping artifacts

The results in Fig. 5 through 8 clearly demonstrate that the proposed correction methods effectively reduce the cupping artifacts. Linearization of the projection data prior to reconstruction effectively extends the dynamic range, enabling the reconstruction of denser and thicker objects. In the present study, effective correction was achieved on both scanners for both PMMA and aluminum test rods. Overall, the correction using the aluminum calibration phantom and the second order polynomial fit to the calibration data (poly-Al) yielded the best results overall, with |percent cupping|<5.9% for all PMMA test rods and <4.1% for all aluminum rods in diameter smaller than 38 mm.

For the large diameter (50.8 mm) Al-test rod the Ultra scan, corrected with poly-Al yielded good correction (|percent cupping|=8.5%), however the scan of the same rod on the speCZT scanner could not be corrected effectively (|percent cupping|=45%), which was in fact worse than the uncorrected result. This result is attributed to the fact that a higher-energy spectrum was (120 kVp) was used on the Ultra scanner, compared to the 90 kVp used on the speCZT scanner. At the lower energy, very few x-ray photons penetrate the full thickness of the Al rod, resulting in “under ranging” – the signal (I) at the detector is nearly zero and it is impossible to differentiate between different thicknesses based on I. In this case the correction step only amplifies the noise in the image and results in even worse non-uniformity across the sample. Therefore, it is important to select a scanning protocol with the highest possible x-ray energy spectrum and highest current and to check for the sufficient photon flux though the highest attenuating path prior to applying correction. (Note, the speCZT scanner is capable of generating x-ray spectra with energies as high as 110 kVp. However, for this study we chose to use 90 kVp in order to demonstrate the ability to correct dense materials, including geomaterials, even at this low energy.)

Although we have found the aluminum calibrated polynomial fitted beam-hardening correction to be the most effective method in the present study, it may be important to evaluate different fitting functions for different scanners and protocols. In addition, even for the scanners used in this study the calibration step must be repeated for each protocol (various kVp) that is to be used.

Furthermore, it is important to note that the presented correction method can be applied retrospectively to existing data – provided the projection images are available – offering the ability to re-evaluate previously scanned specimens.

A goal of this study was to determine scanning and beam hardening correction protocols that will allow micro-CT scanners designed for biological applications (e.g. imaging of experimental rodents) to be applied to geomaterials for non-destructive testing. In general, biomedical micro-CT scanners are less expensive, compared to specialized industrial CT scanners, and are ubiquitous at academic institutions. By demonstrating that the same beam-hardening correction method (poly-Al) can be used to correct both dense objects (aluminum rods) and objects mimicking soft biological tissues (PMMA rods) the scanners can implement a single correction for both biomedical and geological applications.

In addition, biomedical scanners also offer the advantage that a majority employ rotating source-detector geometry, allowing the sample to remain stationary during the scanning process. By minimizing the disruption of sensitive samples, the use of biomedical scanners for geophysical samples may extend the use of CT to an increasing number of applications, such as pore and flow analysis of soils.

4.4 The presented correction in the context of prior beam-hardening correction algorithms

Semi-empirical corrections for beam hardening have been described and adopted in the past, including the approach originally described in 1976 [14] and later expanded by Herman [8], who demonstrated excellent correction using a quadratic polynomial fit. Despite the introduction of correction algorithms almost as soon as CT was invented, work in this area continues to this day – even in the medical field. Beam hardening can be reduced at the cost of photon flux by using the filtration and has been shown to reduce cupping in water phantoms [15]. A post-reconstruction method using re-projections of CT objects and corrections using curves based on aluminum has been shown to reduce beam hardening in aluminum up to 14% [10]. While these methods can also be applied to geomaterials, they employ intermediate stage corrections performed on the original or re-projected CT data, and often require knowledge of the CT reconstruction process to apply the correction. The simple correction methods presented here can – in general – be applied to the projection data off line, then returned to the scanner for CT reconstruction.

4.5 Application to a meteorite sample

Corrected images of the scanned Bassikounou (H5 chondrite) meteorite sample (Fig. 9) clearly demonstrate that the beam-hardening artifact (cupping) was minimized following correction. If a single threshold value is selected as a cut-off for differentiating metallic inclusions from background material the majority of the inclusions in the centre of the sample are not counted (see profiles in Fig. 9c). In fact, even the brightest (densest) inclusion seen in the middle of the selected profile is barely included in the cumulative volume count in the uncorrected images; however, it is clearly above the threshold when the image is corrected. Earlier destructive studies, using thin sections, studying the composition of another sample from the same meteorite have reported metallic content at 8% by volume, with an additional 6.6% by volume of polycrystalline troilite [13] – both of which would appear as radiodense inclusions in the CT scans. The agreement between this prior report (14.6% by volume as the sum of metal plus troilite) and the results from the corrected micro-CT image (Fig. 9b, Table 2) of 13.7% by volume is excellent and suggests that beam-hardening corrected micro-CT using a biomedical scanner can be used for non-destructive analysis of meteorite samples. These encouraging results also suggest that additional micro-CT analysis – e.g. using dual energy approaches [16, 17] – could be used for even more detailed composition analysis. The methodology presented here has been used to correct reconstructions of a range of materials, from spinal tissue [18] to meteoritic samples [19, 20]. The restoration of grayscale values more representative of the material’s given x-ray properties enables effective quantitative analysis that is normally not feasible in CT data of radiodense materials scanned with low-energy polychromatic x-ray sources.