Measuring Identification and Quantification Errors in Spectral CT Material Decomposition

Abstract

Material decomposition methods are used to identify and quantify multiple tissue components in spectral CT but there is no published method to quantify the misidentification of materials. This paper describes a new method for assessing misidentification and mis-quantification in spectral CT. We scanned a phantom containing gadolinium (1, 2, 4, 8 mg/mL), hydroxyapatite (54.3, 211.7, 808.5 mg/mL), water and vegetable oil using a MARS spectral scanner equipped with a poly-energetic X-ray source operated at 118 kVp and a CdTe Medipix3RX camera. Two imaging protocols were used; both with and without 0.375 mm external brass filter. A proprietary material decomposition method identified voxels as gadolinium, hydroxyapatite, lipid or water. Sensitivity and specificity information was used to evaluate material misidentification. Biological samples were also scanned. There were marked differences in identification and quantification between the two protocols even though spectral and linear correlation of gadolinium and hydroxyapatite in the reconstructed images was high and no qualitative segmentation differences in the material decomposed images were observed. At 8 mg/mL, gadolinium was correctly identified for both protocols, but concentration was underestimated by over half for the unfiltered protocol. At 1 mg/mL, gadolinium was misidentified in 38% of voxels for the filtered protocol and 58% of voxels for the unfiltered protocol. Hydroxyapatite was correctly identified at the two higher concentrations for both protocols, but mis-quantified for the unfiltered protocol. Gadolinium concentration as measured in the biological specimen showed a two-fold difference between protocols. In future, this methodology could be used to compare and optimize scanning protocols, image reconstruction methods, and methods for material differentiation in spectral CT.

1. Introduction

Spectral (multi-energy) computed tomography is an extension of dual-energy computed tomography (DECT) or single energy CT. In spectral CT, a single broad X-ray spectrum is compartmentalized into separate energy bins by photon counting detectors to extract energy information to improve contrast discrimination and material decomposition [1]. In recent years, spectral CT with single photon energy sensitive detection systems has offered the potential for a range of CT applications, such as spectral imaging of vulnerable atherosclerotic plaque [2,3], discrimination between iodine and barium [4,5], soft tissue quantification [6,7], beam hardening effects and metal artefact reduction in orthopedic applications [8], quantitative imaging of excised osteoarthritic cartilage [9], and K-edge imaging of high-Z (atomic number) biomedical nanoparticles [10,11,12,13,14,15,16,17,18,19].

The potential to use spectral information for decomposing the linear attenuation coefficient into two or more basis for material discrimination (also called material decomposition MD) was proposed by Alvarez and Macovski in 1976 [20]. They analyzed the material decomposition in the projection image domain (pre-reconstruction). Later in 1977, Brooks [21] successfully demonstrated the dual energy material decomposition method in the reconstruction image domain (post-reconstruction). However, both of these early methods performed decomposition only into photoelectric and Compton components and assumed an absence of K-edges within the energy bins which is not always the case. This method is only valid for low atomic number materials (Z < 20) whose K-edge is not detectable in the human diagnostic energy range.

Elements with more than 40 protons in the nucleus can be referred to as heavy atoms or high Z materials. Atoms heavier than this feature a sharp discontinuity in their X-ray attenuation profile, caused by the binding energy of the innermost (K) shell electrons. Only photons with energies higher than the binding energy have the possibility to knock out the K-shell electrons. Accordingly, photons with energies just above the K-edge are attenuated more than just below. Spectral CT imaging using K-edge energies was first investigated by Riederer and Mistretta in 1977 [22]. Subsequently, the decomposition of the high-Z materials has been evaluated in projection space [15,23,24] and in image space [5,6,25]. The objective of both projection or image based decomposition is to produce basis images with an optimum signal to noise ratio [26]. The main advantage of spectral decomposition is that it allows identification and quantification of multiple high-Z elements in the same scan with high sensitivity and specificity. However, the amount of spectral information can be affected by various factors such as the limited number of independent energy thresholds, incorrect arrangement of the energy bins and their separation size [26], charge sharing [27], charge trapping [28,29,30], K-fluorescence X-rays from high-Z detectors [31], scatter X-rays, semiconductor detector inhomogenities, and pulse pileup [32,33] in the photon counting detectors. These effects are intrinsic to detector properties and some of them can be compensated by flat field correction. However, it is the energy-dependent effects that cannot be corrected by the flat-fielding process, thereby contributing to errors in identification of materials.

Despite considerable development in material decomposition algorithms made over the last few years [6,24,33,34], almost no work has been done on the quantitative assessment of MD images. Quality of material images is typically assessed through residual images from simulated or phantom data [35,36,37] and their SNR [36]. These techniques give a view of the average quantitative and qualitative image quality, but not an accurate classification on the voxel scale. Thus it is important to know the degree of confidence which can be placed on the identity and quantity of any material provided by spectral CT. The advantage of spectral imaging over current imaging techniques is to differentiate and quantify multiple biomarkers and drug penetration into the target tissue at the same scan. At present, spectral CT imaging requires phantoms containing known quantities of known reference materials in order to calibrate the biological data. However, spectral CT decomposed images that show good correlation over the full range of measured concentrations of the materials may contain unrecognized misidentification errors between the different materials. The ground truth for the identity and concentration of all materials in a complex biological specimen is frequently unknown. Therefore, a need has arisen to both measure and monitor changes in material detectability in the post-MD phantom images before they can be used for calibration purposes in specimens whose composition is unknown.

The aim of this study is to generate a quantitative metric using MD images of a calibration phantom which can be used to quantify correct material identification at various concentrations. We compare two different imaging protocols using one sparse solution material differentiation method for post material decomposition assessment by measuring the sensitivity (true positive rate) and specificity (true negative rate) of individual high-Z material at each concentration in the material image domain.

2. Materials and Methods

2.1. MARS Spectral Scanner

MARS spectral CT uses the Medipix3RX photon-counting X-ray detector technology [14,38,39,40]. The scanner is suitable for preclinical spectral imaging of small animal models and human samples of diseased tissue. It consists of a rotating gantry, detector, X-ray source and sample holder shaft, computer hardware and software as shown in the scanner schematic in Figure 1a. It is the photon processing X-ray detector, Medipix3RX [41], that allows for spectral computed tomography imaging by measuring the number and energy of each X-ray photon, to provide specific identification and quantification (g/cm³) of the component material. Medipix3RX allows simultaneous acquisition of up to eight energy bins at a single X-ray exposure. Energy ranges are user-defined, therefore offering flexibility and customization for a wide range of medical and industrial applications. To reduce the effect of charge sharing [42], the charge summing feature is included in the Medipix3RX ASIC which allows the total charge for a single interaction to be summed and allocated to the closest pixel [43]. To do this, each pixel communicates with its four neighbours to locate the pixel with the highest charge for a co-incident event, the total charge is then allocated to this pixel [42].

The threshold settings for the selection of energy bins have a major impact on spectral image quality in terms of image contrast and noise level. The wider the energy bin width is, the lower the noise level is, but the poorer the reconstructed spectral image contrast is. This is mainly because Medipix3RX detector works as a single photon counting detector using user-specified energy thresholds to gather data for a certain energy range. Figure 1b shows the mass attenuation coefficients of calcium, gadolinium (K-edge energy = 50.2 keV), tissue and adipose and simulated X-ray spectrum with an external 0.375 mm brass filter.

The MARS scanner incorporated a microfocus poly-energetic X-ray source and a Medipix3RX detector within a continuous rotating gantry. We used a 2 mm thick CdTe sensor (14 × 14 mm²), bump-bonded at 110 μm to a Medipix3RX readout chip.

2.2. Experimental SETUP

The 720 circular projections over 360° were acquired using a Source-Ray SB-120-350 X-ray tube (Source-Ray Inc., Ronkonkoma, NY, USA) with a tungsten anode having 1.8 mm of aluminium (equivalent) intrinsic filtration; the focal spot size of the X-rays was ~50 μm. Two imaging protocols were used for post-MD assessment: protocol-1 employed a 0.375 mm external brass filter to mimic a clinical CT X-ray energy range and protocol-2 featured no external filtration as often used in pre-clinical X-ray CT imaging (shown in

Table 1. Imaging parameters of two protocols.

Imaging Protocols	Tube Voltage (kVp)	Tube Current (µA)	Exposure Time (ms)	External Filtration (mm)	Intrinsic Filtration (mm)	SDD (mm)	Energy Bins (keV)
Protocol-1	118	24	220/frame	0.375 Brass	1.8 (Al equivalent)	250	30–45, 45–60, 60–78 and 78–118
Protocol-2	118	30	130/frame	-	1.8 (Al equivalent)	250	27–33, 33–49, 49–60 and 60–118

). In this study, four charge summing counters were used for each of the protocols to set low energy thresholds. Vertical dotted lines in Figure 1b is an illustrative example of protocol-1. The source to object (SOD) and source to detector distances (SDD) were 200 and 250 cm, respectively, providing a magnification factor of ~1.25. Five vertical positions of the CdTe camera were used to create a virtual detector of the greater area to cover the 31 mm field of view (FOV). The bias voltage applied to the sensor was −600 V. The MARS camera readout system using Medipix3 has been described previously [44]. Threshold equalization with respect to the noise edge and energy calibration of the detector was performed prior to the measurements. Threshold equalization involves per-pixel calibration to reduce the intrinsic variations in the threshold levels associated with each counter and provides a consistent individual pixel response to a uniform X-ray flux [45]. The energy response of the detector was calibrated against the reference energies such as X-ray fluorescence (XRF) emitted from metallic foils (Mo and Pb) and γ-rays emitted from an Am-241 radioisotope [31].

2.3. Imaging Samples

Two phantoms were used in this study. The first phantom (as shown in Figure 2) was a calibration phantom (Gd/HA phantom) to demonstrate the application of our proposed technique on two different image acquisition protocols. The second phantom with a biological specimen was used to show the use of the technique to exploit the advantages of monitoring changes in detectability of multiple high-Z materials for their correct identification. The calibration phantom and the two biological specimens were scanned with two imaging protocols as described in Table 1. Phantom preparation details are as follows

• Gd/HA phantom: A 31 mm-diameter polymethyl methacrylate (PMMA) phantom was custom-built, containing 6 mm diameter vials for a range of concentrations of gadolinium (1, 2, 4 and 8 mg/mL) and 6 mm diameter solid rods (20 mm long) of calcium hydroxyapatite (HA) (54.3, 211.7 and 808.5 mg/mL) (Quality Assurance in Radiology and Medicine (QRM) GmbH, Moehrendorf, Germany) along with water and a fat surrogate (vegetable oil) as shown in Figure 2a. Hydroxyapatite rods and water serve as bone-like and soft tissue-like material respectively. The gadolinium solutions (1, 2, 4, 8 mg/mL) were prepared by serially diluting Multihance (Gadobenate dimeglumine, Bracco Diagnostics Inc., Princeton, NJ, USA) in water.
• Demonstration with biological specimens: Biological specimens were scanned to demonstrate the effectiveness of our proposed technique for correct material identification and quantification. Two cylindrical cartilage-bone plugs harvested from the bovine stifle joint, each with a diameter of 8 mm and length of 7 mm, were used to quantify gadolinium uptake in articular cartilage. The sample preparation setup was similar to that used in a previous study using iodinated contrast [9]. Samples were incubated for 24 h in Gd contrast agent solution (100% MultiHance) at 37 °C and then rinsed in PBS for 30 s. A non-Gd incubated sample was used as a control. Both incubated and control samples were placed side by side in a 15-mL Falcon tube, with PBS infused cotton wool at both ends, sealed to maintain a humid environment during scanning.

2.4. Image Processing and Reconstruction

Flat field measurements (1440 for each camera position and energy bin) were taken after the sample scan to correct for fixed-pattern variations in individual pixel responses. The raw data in DICOM format was transferred from the scanner’s server to the scanner’s inbuilt patient archive and communication system (PACS), where an automated image processing was performed.

The image processing has three stages: pre-reconstruction processing, reconstruction into attenuation volumes, and decomposition into material volumes (see Section 2.5). In this case, the image processing produced slices with 0.1 mm² voxels on a 480 × 480 matrix.

The pre-reconstruction processing consists of three steps; the first step is pixel masking to remove any bad pixels; the second step is flat field correction [32,46]; the final step is ring filtration. The ring filter is based on the work by Jan Sijbers and Andrei Postnov [47], but applied to the projection data and adapted into 3D to exploit more information provided by the detector (x, y, and theta).

The reconstruction uses a low-resolution version of the polychromatic form of the Beer–Lambert Law. What this means is that the overlapping, low-threshold counters that measure the data are processed simultaneously to produce attenuation volumes representing non-overlapping energy bins across the measured spectrum. For example, four counters representing the energy range 30–120 keV, 45–120 keV, 60–120 keV and 78–120 keV would simultaneously produce four attenuation volumes representing the energy ranges 30–45 keV, 45–60 keV, 60–78 keV and 78–120 keV.

Processing the counters simultaneously exploits all photons available in the scan. Consider the attenuation volume for 78–120 keV. All four counters contribute to this volume. The counter for 60–120 keV provides specific information about the energy range but has fewer photons (higher statistical noise). The counter for 30–120 keV is less specific but contains all the photons in the scan (lower statistical noise). Given that the counter’s measurements are acquired simultaneously, we know each counter maps the same structures. This reconstruction process uses all photons available in the wider counters to reconstruct clean boundaries and structures, while using the narrower counters to guarantee that the correct energy response is obtained.

The reconstruction algorithm itself is a statistical iterative technique. It uses multiplicative correction terms similar to OSEM [48], or MART [49]. The statistical element is introduced by weighting the corrections to the volume by a normalized Poisson distribution to slow down convergence when approaching the solution, thereby reducing noise. The reconstruction algorithm also adopts a multi-stage approach where it initially reconstructs voxels that are eight times larger than requested. Later on, this is repeatedly subdivided for a total of four stages until the requested voxel size is reached [50]. This approach allows for the reconstruction to proceed quickly. It is also a weak form of a sparsity constraint as a large voxel is the same as a set of small voxels with the same value. Lastly, the larger the voxel, the more pixels from the projection images will contribute to it. This reduces the effect of dead regions. This is particularly useful during the initial reconstruction stages where the effects of dead regions are the most significant.

To evaluate the spectroscopic response, CT slices corresponding to a thickness of ~0.5 mm were averaged and then analyzed. The CT numbers were quantified in spectral Hounsfield Units [51] by normalizing CT attenuation (commonly known as linear attenuation coefficient µ) to water and air to account for the differing attenuation at different energies as defined by Equation (1)

$$H{U}_{Energybin\left(E\right)}=1000\frac{{\mathsf{\mu }}_{mat\left(E\right)}-{\mathsf{\mu }}_{water\left(E\right)}}{{\mathsf{\mu }}_{water\left(E\right)}-{\mathsf{\mu }}_{air\left(E\right)}}$$

(1)

where ${\mathsf{\mu }}_{mat}$ is the attenuation measurement of the given material in the given energy bin, i.e., $E=1,2,3\&4$.

The relationship between CT number and energy bins was analyzed graphically by manually selecting ROIs comprising ~1409 voxels (N) over each solution insert. The standard error was calculated by the normalized standard deviation $(stddev/\sqrt{N})$. The relationship between signal intensity and material concentration was also determined for all four energy bins by linear regression. This was achieved by plotting CT number as a function of concentration.

2.5. Material Decomposition

Material decomposition is the process of converting spectral attenuation (energy information) into information about the constituent materials contributing to that attenuation. This typically involves obtaining material properties inverting the mass attenuation equation for one of its variants (i.e., volume fraction [52], Compton-photoelectric [20], or $\rho Z$ [53,54]). The weighted average X-ray attenuation of a compound or mixture of the constituent materials is given by

$${\mathsf{\mu }}_E={\displaystyle \sum }_{i=1}^m{\rho }_m{\left(\frac{\mathsf{\mu }}{\rho }\right)}_{mE}$$

(2)

In this equation, ${\mathsf{\mu }}_E$ is the linear attenuation of some composite material for a given energy range $E$. The material properties desired are the densities of the constituent materials (indexed by $m$). These are proportionally connected to the linear attenuation through the mass attenuation of the respective materials ${\left(\mathsf{\mu }/\rho \right)}_{mE}$ for the given energy ranges. Since attenuation varies differently over energy for different materials, such material properties can be deduced when using multiple energies.

The MARS small animal spectral CT scanner incorporates a heuristic image space material decomposition algorithm (MARS-MD) for converting reconstructed energy bins into sparse material images [55,56,57]. This algorithm follows a three-step process for assigning materials to each voxel: (1) each voxel classified as either high, low, or insignificantly attenuating; (2) a calculation of uniquely constrained least squares decomposition for each feasible pair of materials falling into the respective attenuation category; and (3) a classification of the result which has the smallest regression error. The output of this process is that every voxel is assigned the 1–2 materials with the most significant contribution to the voxel’s attenuation. The algorithm was calibrated using mass attenuation coefficients estimated from the reconstructed phantom data. For more details on this algorithm we direct the reader to our draft paper [57]. Using this algorithm, the spectral data analyzed in this study was decomposed at a voxel level into sparse combinations of gadolinium, hydroxyapatite, fat, and water. We are exploiting the single high-Z property of the MD algorithm to generate a quantitative metric in decomposed material images which can be used to both measure and monitor changes in detectability and correct identification of multiple high-Z materials over a wide range of concentrations.

To measure the amount of misidentification between different materials, a metric was generated for quantitative evaluation of post-MD volumes of the calibration phantom which provides sensitivity (true positive rate) and specificity (true negative rate) for correct material identification at various material concentrations. As the composite material contains high-Z material in solution form (i.e., in water), all voxels of high-Z material density images should also be identified as water on the water density image. Therefore, water density image was taken as a reference for our metrics. In addition, voxels in all material images with value >0 were masked for further processing.

In our metric, sensitivity refers to the test’s ability to quantify successful voxel identification in the target material’s region of interest (ROI). Manual selection of a ROI comprising 1409 voxels (N) over each of the given solution inserts of the target material image was compared with the corresponding voxels in the water density image, as shown in Figure 3. Corresponding voxels with values equal to 1 in both images (target material and water density) were identified as successful voxels (voxels that are identified as the correct material).

The percentage of successful voxels identification in the target material’s ROI provides the true positive voxels (TP) such that

$$TP=\left\{v:vinR,TM\_v=W\_v\right\}$$

(3)

Here $v$ is a voxel and R shows the ROI, TM_v is the target material image voxel, and W_v is the water material image voxel. Similarly, the percentage of unsuccessful voxels provide the false negative voxels (FN) in the target material’s ROI such that

$$FN=\left\{v:vinR,TM\_v\ne W\_v\right\}$$

(4)

The sensitivity (identification as a percentage) of every voxel of each material at each of the given concentration was calculated by Equation (5).

$$Sensitivity=\frac{TP}{TP+FN}\times 100$$

(5)

The specificity refers to the test’s ability to quantify voxels which are correctly identified as not the target material in the control ROIs. The specificity describes that the control material will be identified in the target material image when misidentification occurs.

Manual selection of ROIs in the target material image comprising 1409 voxels (N) over each of the given solution inserts of the control material was compared with the corresponding voxels in the water density image as demonstrated in Figure 4. Corresponding voxels in control solutions in the target material image with values not equal to 1 were identified as successful voxels (correctly identified as not containing the target material).

The percentage of successful voxels which are correctly identified as not the target material in the control ROI provides the true negative voxels (TN) such that

$$TN=\left\{v:vinR,CM\_v=W\_v\right\}$$

(6)

where $CM\_v$ shows control material image voxel. Similarly, the percentage of unsuccessful voxels (incorrectly identified as containing the target material) provide the true false positive voxels (FP) such that

$$FP=\left\{v:vinR,CM\_v\ne W\_v\right\}$$

(7)

Specificity of every voxel of each material at each of the given concentration is calculated by Equation (8)

$$Specificity=\left(1-\frac{TN}{TN+FP}\right)\times 100$$

(8)

It is important to mention that $TP+FN$ (Equation (5)) and $TN+FP$ (Equation (8)) should sum to the total number of the selected voxels. However, the specificity of target material to water cannot be determined by Equation (8) as the water basis is used as a reference material. Therefore, percentage of voxels in the composite materials which are only classified as water are calculated as

$$SpecificityofTMtowater=\frac{FN-FP}{Totalnoofvoxels}$$

(9)

Finally, the relationship between the known and measured concentrations of gadolinium and hydroxyapatite was analyzed in the material domain images by manually selecting ROIs comprising ~1409 voxels (N) over each solution insert.

3. Results

To calibrate biological data against the phantom containing known concentrations of known materials, calibration phantom images were evaluated in the energy domain (spectroscopic images) as well as the material domain (material segmentation and quantification) before they could be used for calibration purposes.

The overall image quality produced by the MARS scanner was adequate for the analysis task, i.e., no major artefacts (beam hardening, partial volume effect, photon starvation or undersampling) were observed in the images. A few ring artefacts were observed towards the center of the reconstructed images. Figure 5 shows the spectral response of the transverse slices of the Gd/HA phantom with 4 energy bins for two imaging protocols.

Figure 6 shows the relationship between CT number and energy bins. The CT number of the gadolinium and hydroxyapatite decreases with energy due to the decreasing influence of the photoelectric effect and increasing Compton scattering, however, for both protocols, a higher value for CT numbers is noticeable due to the influence of k-edge at gadolinium’s k-edge energy range. Moreover, the CT number of the lipid (effective Z ≈ 6.5) increases uniformly with energy. The spectral separation between different materials is the key information required for the segmentation task of the material decomposition.

The linear fitted lines in Figure 7 indicate the system’s linearity response over a range of concentration and energy bins ($R^2\ge 0.99)$. Linearity determines the ability of a system to detect the presence or absence of any material. This information is important for material quantification and is therefore used to correlate between the known (contrast material) and unknowns (suspected regions). It can be affected by random noise, beam hardening, and photon starvation [58]. The non-zero CT number at zero concentration can be related to possible imperfections of the photon counting detectors or contributions of random (quantum) noise and systematic effects of the system.

Figure 8 and Figure 9 show identification of hydroxyapatite, gadolinium, lipid and water in Gd/HA phantom with protocol-1 and protocol-2 respectively. Some misidentification in material voxels (crosstalk) is observed between different materials, especially at the lower concentrations of gadolinium and hydroxyapatite decomposed images. It is important to highlight that despite clear material separation in attenuation profiles of gadolinium and hydroxyapatite in energy images (see Figure 6 and Figure 7), the misidentification of material type and its concentration is difficult to assess visually in the images in Figure 8 and Figure 9. It is also observed that the phantom body (PMMA) was significantly misidentified as gadolinium with protocol-2 (see Figure 9c), hence the need for a voxel by voxel assessment method for quantifying the degree of misidentification.

For post material decomposition evaluation, images corresponding to a thickness of ~2 mm were analyzed slice by slice by using Equation (5) for sensitivity, Equation (8) for the specificity of the target material with high-Z materials and Equation (9) for the specificity of the target material with water. Results were then averaged and plotted in the MD misidentification charts as shown in Figure 10 and Figure 11.

Once the material separation was quantified, the relationship between the known and measured concentrations of gadolinium and hydroxyapatite was determined from material images. Results indicate that despite good spectral and linear correlation of gadolinium in energy images with protocol-2 (Figure 6c and Figure 7c), the measured concentration of gadolinium is underestimated for protocol-2 (Figure 12a) whereas, protocol-1 indicates a good correlation between the known and measured concentrations over the full range of gadolinium $(R^2\ge 0.99).$ Hydroxyapatite concentrations were observed to be linear with both protocols. Material misidentification and mis-quantification results from Figure 10, Figure 11 and Figure 12 are summarized in

Table 2. Qualitative and quantitative results summary of the Gd/HA calibration phantoms with two imaging protocols.

Known Material (mg mL−1)	Protocol-1 Measured Material (mg mL−1)	Protocol-2 Measured Material (mg mL−1)	Protocol-1 Measured Misidentification (%)	Protocol-2 Measured Misidentification (%)
Gd08	8.19	3.15	0	7.22
Gd04	4.35	1.29	1.52	31.68
Gd02	2.31	0.59	8.93	41.45
Gd01	1.18	0	37.57	57.66
HA808.5	818.9	852.9	0	0
HA211.7	215.4	246.4	0	0
HA54.3	63.78	68.09	4.66	0

. These results demonstrate that the methodology can be used to compare how different scanning protocols using one material differentiation method affect the measurement of the concentration of materials.

Spectral CT and material component images of ex-vivo bovine cartilage-bone plugs are shown in Figure 13. The bone in the bovine plugs is represented in the HA component images (see Figure 13b,f). The layer of articular cartilage of the incubated knee specimen is represented by the gadolinium density images (see Figure 13c,g) whereas no gadolinium signal from control specimen was detected with either protocol. However, misidentification of gadolinium as water is observed with protocol-2.

Figure 14a,b shows fused images of gadolinium and hydroxyapatite densities with protocol-1 and protocol-2 respectively. Quantitative variations in the gadolinium concentration are plotted (Figure 14c) along the depth of the cartilage. Contrast quantification is crucial in a biological specimen especially when radiologists determine disease burden and contrast delivery by using high-Z contrast materials. Figure 14c shows that despite good qualitative gadolinium decomposition, measured gadolinium concentration with protocol-2 is over 50% lower than with protocol-1. Which measurement of these two concentrations is more reliable? Protocol-2 underestimates gadolinium concentration in comparison to Protocol-1 in the phantom (Figure 12a), by a similar amount.

4. Discussion

Material misidentification affects the overall sensitivity and specificity of the measurements which is the key requirement of specific identification and quantification of multiple tissue components for spectral molecular imaging. The MARS MD algorithm assumes that each voxel consists of primarily one high-Z material. Therefore, in the segmentation task, once a material is identified in any given voxel, the value for this voxel in other basis materials is assumed to be zero. Therefore, this particular analysis method was designed using a sparse solution MD as a reference, only taking into account whether the solution is zero/non-zero in each material and therefore exploits this “binary property” and quantifies correct material identification at various material concentrations. However, for MD algorithms that have noise in all material channels (i.e., every voxel contains a non-zero value for all material channels), modification like low-value thresholding or smoothing using techniques such as compressed sensing can be applied. Figure 10 and Figure 11 show that sensitivity of any given material in the post-MD images is density dependent whereas specificity is not, and therefore specificity ideally should be 100% for every voxel of each material at each of the given concentrations.

There could be various factors that likely to have an effect on sensitivity, specificity and material mis-quantification. The statistical noise of the system could be one of the factors affecting the sensitivity and specificity of the given material. Noise not only affects the ability to detect low-Z materials but also affects the differentiation of low concentrations of high-Z materials [59]. Another noise related factor affecting material decomposition is pulse pileup [33] which contributes negatively to the energy calibration. It has its origin in photon statistics and also depends on the read-out speed of the detector electronics [60,61]. Pileup occurs when two or more photons arrive within the same detector region and within the same integration time of the readout electronics. These events appear to be only one thus resulting in an information loss. It also leads to events in low energy ranges being detected which exceeds the maximum possible photon energy which also complicates the process of energy calibration and results in a loss of energy resolution [31]. Low X-ray flux was used with the brass filter protocol (protocol-1) to reduce the probability of the pileup effect whilst providing adequate statistics in a reasonable exposure time as shown in our previous studies [32,62]. This might be the reason for less material misidentification in protocol-1 than protocol-2. The underestimation of the gadolinium quantity with protocol-2 can be associated with higher gadolinium misidentification as hydroxyapatite. The threshold energy settings for the selection of energy bins could also have a major impact on spectral image quality in terms of image contrast and noise level. Wide energy bins have lower noise but suffer from poor contrast in the reconstructed images. In single photon counting detectors, a photon is only counted if the signal pulse is higher than the user specified energy threshold. However, the signal created by a single X-ray photon could be registered simultaneously by all energy thresholds depending on photon energy. Therefore, the selection of energy ranges needs to be balanced with the requirements for achieving desirable contrast and material’s K-edge. Moreover, CT reconstruction for spectral CT imaging needs to be performed in subtracted energy bins to reconstruct the difference of the photon counts between the two subsequent counters.

It is difficult to pinpoint the exact cause of misidentification but factors such as X-ray scatter, Compton scattering, K-fluorescence X-rays from high-Z semiconductor materials and imperfections of the CdTe detector (pulse pileup, charge trapping, semiconductor detector inhomogenities, and crosstalk within the semiconductor material) are likely to play a part. All these effects are intrinsic to detector properties and not all of them can be compensated by flat field correction, especially energy-dependent effects. They contribute to misidentification of materials in material decomposition. Image related artefacts also have negative effects on the material decomposition.

Other groups [25,63,64] have reported noisy regions with incorrect material identification of high-Z materials but they claimed their ability to detect low concentrations of high-Z materials with high accuracy by measuring the concentrations in decomposed images and comparing it to the true concentration value. Our results show that good correlation over the full range of the known and measured concentrations fails to provide information about the misidentification between different materials and the extent of such misidentification in decomposed material images. Results from the misidentification chart can be used to choose contrast materials best suited to a spectral CT imaging task in relation to the anticipated concentration of those materials. We can use the binary information to quantitatively compute sensitivity and specificity performance of the (data acquisition and data processing) system for target materials of different concentrations and composition and scanned with different imaging protocols.

The materials selected for this study are commonly used for radiology related studies due to their attenuation properties. A shortcoming of the study could be the evaluation of a limited number of high-Z materials and limited image acquisition protocols. Despite these limitations, we showed that the technique can be applied to various concentrations of high-Z materials simultaneously and can be used (1) as a quality assurance tool to set up a threshold for lower detectable limits of any given material in the post-MD domain; (2) for optimization of the spectral CT imaging parameters and acquisition protocol to maximize the sensitivity and specificity of the given materials; and (3) to establish a baseline of misidentification measurement for daily QA of spectral CT scanners.

The ground truth for the identity and concentration of all materials in a complex biological specimen is frequently unknown. There are many potential sources of error in spectral CT analysis of compositionally complex specimens, including the MD method itself. This paper describes a new method that measures the net effect of all these sources of error on the ability to identify and quantify materials at spectral CT. We have successfully tested the method in both a phantom and a biological sample. We have demonstrated its utility for one material differentiation method applied to two scanning protocols. In future, this methodology could be used to compare and optimize scanning protocols, image reconstruction methods, and methods for material differentiation in spectral CT.

5. Conclusions

In summary, we developed and validated a method for quantifying misclassification of the identity and quantity of gadolinium and hydroxyapatite in spectral CT material images. The method can be applied to any chosen material to measure the image quality of the final result. Although developed with multi-energy CT, it is equally applicable to dual-energy CT. These results demonstrate that the methodology can be used to compare how different scanning protocols using one material differentiation method affect material identification and measurement of the concentration of that material for clinically relevant contrast agents and tissues.