2.1

ME-CBCT setup

The experiment in this work was based on the SmART preclinical radiation platform shown in Fig.1, which included a rotating C-arm gantry and an small animal couch assembled in a self-shielded cabinet. The gantry can perform CBCT scanning using an x-ray tube (Comet iVario 225 kV, JME Ltd, Suffolk, UK) mounted on the gantry and an amorphous silicon flat-panel detector (FPD) (XRD 0820 AN3-ES, Perkin-Elmer, Wiesbaden, Germany). An XC-Thor Cadmium telluride (CdTe) PCD (Direct Conversion AB, Sweden) was mounted on an aluminum frame attached to the gantry. The PCD has an active sensor area of 51.2×100.0 mm² with a 0.1×0.1 mm² pixel size. The source to isocenter distance was 30.5 cm and source to detector distance was 44.2 cm. This PCD was able to output the number of photon counts for each pixel above a user-specified energy threshold T for a time interval Δt.

Figure 1.

Installation of the PCD on the SmART preclinical radiation platform (left) and the front view of of the PCD (right).

We used a 100 kVp x-ray beam for CBCT data acquisition, with a rotation time of 1 min. The data acquisition of the PCD was controlled via the Application Programming Interface provided by the vendor. Specifically, we triggered the PCD to acquire projection data at a frequency of 30 frame/sec and each data frame receive photons for Δt = 20 msec. The remaining ~ 13 msec time interval was allocated for data processing and readout. The energy threshold T was set to sequentially cycle through three energy thresholds T^l, l = 1, …, 3. PCD operated under anti-coincidence high sensitivity mode for all acquisitions. For each projection, the projection angle was read out from the encoder of the gantry motor. With this setting, we acquired 600 projection images for each energy threshold.

The workflow covered by this study is summarized in Fig. 2. We will present details of key steps in subsequent sections.

Figure 2.

Flowchart of PCD-based ME-CBCT imaging procedures.

2.2

Selection of energy thresholds

The lowest energy threshold T₁ was set to be 26 keV to avoid electronic noises.⁹ The other two energy thresholds were determined to balance signal-to-noise ratio (SNR) in log-transformed projection images of all three energy channels, as the SNRs are directly related to those of the subsequently reconstructed images from each channel. As such, we acquired air projection images at 0 degree gantry angle under 3 mA for N_t = 42 thresholds T_i, i = 1, …, N_t ranged from 28 to 110 keV, with increment of 2 keV. Note that the thresholds were nominal input values to the PCD, therefore our search range went beyond the 100 kVp used. At every threshold T_i, we acquired N_f = 200 projections and denote the count image where x is the index of each pixel, k = 1, 2, …, N_f.

For a given threshold combination T^l, l = 1, 2, 3 selected among possible T_i, i = 1, …, N_t values, the count images at the three energy channels E₁ =[T¹, T²] keV, E₂ = (T²,T³] keV and E₃ = (T³,100] keV, were for l = 1, 2 and . Note that we wrote the upper bound 100 keV in the third energy channel because of the 100 kVp x-ray beam, yet in practice this channel represents photons counts with energy above the threshold T³. Following the same procedure, we acquired projection images of a CT calibration phantom, denoted as C_{l, k} (x) for the energy channel l and frame k.

We then computed the log-transformed projection image , where is the averaged counts of the air projection images over all frames. Using the repeatedly acquired 200 projections, we calculated the mean and standard deviation of _gl,k(x), denoted as ḡl(x) and σ_gl(x), respectively. Finally, the SNR at coordinate x was calculated as S_l (x) = ḡl (x)/σ_gl (x). We averaged S_l(x) over 100 pixels inside the phantom region on the projection image, and denoted the averaged SNR as . It was our objective to find proper thresholds such that approximately equal among all energy channels l = 1, 2, 3. As such, we considered the metric γ = max_l / min_l . We enumerated all possible combinations of thresholds T² and T³, and found out the combination that minimized γ as the optimal thresholds for subsequent studies.

2.3

Detector calibration

To calibrate the non-uniform pixel responses of PCD, we acquired air scan under the threshold-sweeping mode as in the actual CBCT scan using the energy thresholds determined from Sec. 2.2. We repeated this with x-ray tube current I in the range of [0, 3.4] mA with 0.2 mA increment. Let us denote the count image at energy threshold T^l, l = 1, …, 3 with the jth tube current , j = 1, 2, …, as . Note that we repeatedly acquired 200 frames and took the average of these frames to reduce noise. Count images at the three energy channels were computed by taking the difference between data with subsequent energy thresholds.

With the count images acquired for each energy channel l = 1, 2, 3 and tube current level j = 0, 1, 2, …, for each pixel x and energy channel, we fit the tube current as a function of counts using and in a polynomial form. The coefficient of determination R² was used to evaluate the goodness of the fitting model. In the fitting procedure, we gradually increased the number of polynomial order, until the resulting R² reached . If this was not achievable at 4^th order polynomial, the pixel was considered as a ‘dead pixel’.

2.4

Image processing, reconstruction and material decomposition

After data acquisition of a CBCT scan, we converted raw count images C_l(x, θ) of energy channel l to the corresponding images of current I_l (x, θ) using the established calibration model of each pixel. Here we explicitly write the projection angle θ. In this step, nearest-neighbor interpolation was applied to handle all ‘dead pixels’. For each flat-fielded projection image at angle θ, we identified a region that was not blocked by the phantom and computed the averaged pixel value in this region as the air norm . Compared to using the nominal tube current in the scan as the air norm, this approach deriving considered the tube current fluctuation among projections during the CBCT data acquisition.

After these stpes, we computed the projection image . Finally, CBCT image f = (f₁, f₂, f₃)^T were reconstructed using the projection data and the GPU-based Feldkamp-Davis-Kress (FDK) reconstruction code,^{10, 11} where f_l, l = 1, 2, 3 represents reconstructed image of energy channel l.

In this study, we considered a two-material decomposition model that expressed each image voxel as a linear combination of two basis materials: water and bone. However, with three energy channels, in principle, it is possible to perform three-material decomposition. With our model, f = Am, where m = (m_w, m_b)^T is a two-component image representing mass density image of the two basis materials, and A ∈ ℛ^(3×2) is the system matrix. The two columns of A are x-ray mass attenuation coefficients of corresponding materials in the three energy channels, which was obtained in the reconstructed CT calibration phantom images. To solve for the material images without amplifying noise, we performed material decomposition by solving an optimization problem

where the first term of the objective function ensured fidelity of the solution to the CBCT images, and the second one was a Total variation regularization term to reduce noise while preserving edges in the solution. λ is a weighting factor controlling the relative importance of the two terms. Note that the operator ∇ was applied to the spatial direction only. This model was solved with the Alternating Direction Method of Multipliers.

2.5

Evaluation

Two phantoms (SmART Scientific Solutions B.V., Maastricht, Netherlands) were utilized to calibrate the ME-CBCT system and evaluate the imaging procedure and image quality. The first one was a preclinical CT calibration insert phantom that contained 10 inserts of different materials plugged into a 30 mm diameter background base representing water. After the CBCT images were reconstructed, we picked two regions of interest with water and bone materials and used the average voxel values normalized by the known material density to construct the material decomposition system matrix A. We then decomposed the CBCT images into the two material images using the model in Eq. (1). The second phantom was a plastinated mouse specimen. Due to the small size of the PCDs (5.0 cm along the rotation axis), only the head was scanned.

3.1

Energy thresholds and detector calibration

With the method described in Sec. 2.2, we computed the SNR ratios γ for all threshold combinations of (T², T³). The SNR ratio that was closest to 1 was γ = 1.05, when the two energy thresholds were T² = 56 keV and T³ = 90 keV. The corresponding SNRs of projection images were 7.07, 7.32 and 6.98 for the three energy channels, respectively.

Fig. 3 presents the impacts of detector calibration. Without calibration, we can clearly see different detector sub-panels to assemble the entire detector. The boundaries among panels would cause ring artifacts in the reconstructed images. We selected two regions of interest (ROIs) that were not blocked by the phantom. Due to flat x-ray beam fluence, the projection data in these two ROIs are expected to be the same. However, without calibration, we observed a large difference between the raw counts in the two ROIs across all projection angles. The calibration step was able to reduce this discrepancy. The fluctuation among projections was ascribed to x-ray tube output fluctuation, which was taken care of before image reconstruction by the choice of projection-specific air norm.

Figure 3.

Left: projection images at θ = 90 degree in energy channel E₁ with and without detector calibrations; blue boxes show the location of two ROIs. Right-top: average counts of the two ROIs. Right-bottom: averaged mA of the two ROIs after detector calibration.

3.2

CBCT images and material decomposition

The reconstructed CBCT images of the CT calibration insert phantom and the mouse phantom with three energy channels are shown in Fig. 4. The object contrast reduced with increasing photon energies. For the calibration phantom, we compared measured x-ray attenuation coefficients and those computed by averaging the known x-ray attenuation data weighted by x-ray beam spectrum in the three energy channels. The difference was 10% on average. The difference can be ascribed to multiple factors including the reconstruction model, the discrepancy between nominal and actual energy thresholds, the ignored detector response variation over energy, etc.

Figure 4.

CBCT images of the CT calibration insert phantom (top) and the plastinated mouse specimen (bottom) in the three energy channel. Display window: [0,0.5] cm^–1. The blue and red circles indicate the solid water and the bone basis materials.

To perform image-based material decomposition, we selected solid water and bone as two basis materials, and obtained the attenuation coefficients by averaging the voxel values over the corresponding ROIs indicated by the blue and red circles in Fig. 4. The mass attenuation coefficients in the three energy channels were [0.240, 0.198, 0.186] cm²/g and [0.489, 0.296, 0.261] cm²/g for water and bone materials, respectively.

Fig. 5 presents the material images of water and bone of CT calibration insert phantom and the plastinated mouse phantom in the axial plane. The two basis materials were differentiated clearly in the material images. Note that the animal skin also appeared in the bone image. The manufacture process of the plastinated mouse caused accumulation of polymer materials at the skin, which appeared as bone-like materials under x-ray imaging.

Figure 5.

Water (left) and bone (right) maps of the CT calibration insert phantom (top) and the plastinated mouse specimen in the axial plane (bottom).