Optimized energy thresholds in a spectral computed tomography scan for contrast agent imaging

ACCELERATOR, RAY TECHNOLOGY AND APPLICATIONS

Optimized energy thresholds in a spectral computed tomography scan for contrast agent imaging

Kai-Xin Huang，

Zhi Deng，

Xiao-Fei Xu，

Yu-Xiang Xing

Nuclear Science and Techniques

Vol.30, No.3

Article number 38

Published in print 01 Mar 2019

Available online 12 Feb 2019

DOI：10.1007/s41365-019-0563-9

32800

Spectral computed tomography (CT) based on photon counting detectors (PCDs) is a well-researched topic in the field of X-ray imaging. When PCD is applied in a spectral CT system, the PCD energy thresholds must be carefully selected, especially for K-edge imaging, which is an important spectral CT application. This paper presents a threshold selection method that yields better quality images in K-edge imaging. The main idea is to optimize the energy thresholds ray-by-ray according to the targeted component coefficients, followed by obtaining an overall optimal energy threshold by frequency voting. A low-dose pre-scan is used in practical implementations to estimate the line integrals of the component coefficients for the basis functions. The variance of the decomposed component coefficients is then minimized using the Cramer–Rao lower bound (CRLB) method with respect to the energy thresholds. The optimal energy thresholds are then used to take a full scan and gain better image reconstruction with less noise than would be given by a full scan using the nonoptimal energy thresholds. Simulations and practical experiments on imaging iodine (I) and gadolinium (Gd) solutions, which are commonly used as contrast agents in medical applications, were used to validate the method. The noise was significantly reduced with the same dose relative to the nonoptimal energy thresholds in both simulations and in practical experiments.

Spectral CTContrast agent imagingCramer–Rao lower boundThresholds optimizationK-edge

1 Introduction

X-ray CT imaging is an important tool in clinical diagnoses, non-destructive testing, security inspections, and other fields. X-ray attenuation is a function of the photon energy and the material properties, but traditional CT detectors cannot identify the photon energy [1]. Therefore, traditional CT systems can only provide an average linear attenuation coefficient in the material for a polychromatic source. Generally, the higher-energy spectral components are less strongly attenuated, while the lower-energy spectral components are more strongly attenuated. This results in the beam hardening effect, which produces shading and cupping artifacts in the CT reconstructions and is inevitable in traditional CT systems.

Spectral CT imaging can acquire signals for multiple energy spectra, so the energy-dependent attenuation can be fully utilized to identify different material properties to discriminate different materials. A current popular realization involves dual-energy imaging [2], which uses kVp-switching techniques with two sources running at different energies [3] and dual-layer detector approaches [4]. Dual-energy imaging can effectively suppress beam-hardening artifacts and yield specific material information, such as the effective atomic number that could not be obtained from traditional CT imaging.

Spectral CT based on PCD has generated much interest in recent years. PCD uses different energy thresholds to identify the energy of each incident photon so that signals collected by PCD carry energy information. K-edge imaging [5, 6], is one of the most important PCD applications. The K-edge energies are used as thresholds, where K-edge imaging is then used to identify and quantify high Z elements within a specimen, an animal, or a patient. Some high Z elements serve as contrast agents for diagnostic X-ray imaging, such as I, Gd, and Au solutions. The raw CT data is then analyzed with different energy bins to enable spectral CT to have a high signal-to-noise ratio (SNR) and spatial resolution. The multiple spectra are then used with material decomposition methods to discriminate between materials. Furthermore, it is also possible to reduce the radiation dose that is used in each scan [7]. The common spectral CT reconstruction methods are pre-reconstruction, post-reconstruction, and "one-step" reconstruction [8, 9].

The pre-reconstruction method decomposes the material information in the projection domain. Dual effect and material basis function techniques are commonly used for the decomposition. The analysis also requires K-edge material attenuation coefficients. The linear combination coefficients of the basis functions are calculated according to the raw projection data collected in different energy bins to reconstruct the spatial distribution of the basis functions. The physical model given by the pre-reconstruction algorithm is fairly accurate. Since a pair of ray paths is simultaneously measured using low and high energy bins, this method can be implemented using dual-source or high-low switching modes. The post-reconstruction method decomposes materials in the spatial domain and is mostly used in medical dual-energy CT [10]. First, the CT images are reconstructed using the low and high energy bin information. Then, the material basis function distributions are obtained based on the reconstructions. This is easily implemented and is not limited by the system scanning method. However, beam-hardening artifacts can affect the image quality. The "one-step" method directly estimates the material reconstruction using the raw measured CT data in the low and high energy bins. Iterative techniques are then normally used for the "one-step" method, which makes the method more complex than the other two methods.

This work focuses on the pre-reconstruction method due to its robustness and usefulness using pairs of matched ray paths, which can be easily obtained in spectral CT.

The PCD signals are used to choose suitable energy thresholds and obtain high quality images. Two empirical methods can be used to determine the energy thresholds. One method simply distributes the thresholds evenly over the entire spectrum. The other method averages the number of incident photons over all energy bins. This study uses Cramer–Rao lower bound (CRLB) theory to calculate the energy thresholds. CRLB theory, originally proposed by Cramer in 1945 [11], provides a theoretical noise bound for unbiased estimates of random variables. In 1976, Alvarez and Macovski used CRLB with a dual-energy CT to analyze the noise in the basis function decomposition [2]. In 2009, Rossel et al. utilized CRLB with spectral CT data to analyze the noise in the material decomposition results in the projection domain and in reconstructions in the spatial domain [12, 13]. They used CRLB to obtain an analytical expression for the lower bound of the decomposition noise. They then optimized the thresholds based on CRLB analysis in the projection domain using line integrals over the basis function coefficients in the special situation for a vertical ray path [14]. Meng et al. optimized the energy thresholds in the spatial domain using image subtraction [15]. Since the line integrals of the decomposition coefficients for the scanned objects are unknown and could vary greatly in practice, this study proposes an energy threshold optimization method for K-edge imaging. The main idea is to first use a low-dose pre-scan to obtain the line integrals of the basis function coefficients. The line integral along a single ray path is then used to reduce the local variance of the basis function using CRLB theory. This yields a minimum variance threshold. Frequency voting is then used to calculate the overall optimal thresholds. Finally, a full scan is taken using the optimal energy thresholds in order to obtain a high-quality image.

2 Theory and methods

2.1 System model of a spectral CT

A key step in spectral CT imaging is to obtain an energy-independent component coefficient distribution for the scanned object. In a general model[16], the linear attenuation coefficient curve for a certain material is decomposed into a linear combination of several basis functions:

μ (\vec{r}, E) = \sum_{τ = 1}^{Γ} a_{τ} (\vec{r}) ϕ_{τ} (E),

(1)

where ϕτ (E) represents the energy-dependent component in $μ (\vec{r}, E)$ and $a_{τ} (\vec{r})$ represents the spatial distribution of the component coefficient at location $\vec{r}$ corresponding to the τ^th basis function.

Therefore, the projection of linear attenuation coefficients in CT imaging can be formulated as:

\int μ (\vec{r}, E) d \vec{r} = \int \sum_{τ = 1}^{Γ} a_{τ} (\vec{r}) ϕ_{τ} (E) d \vec{r} .

(2)

By defining the line integral of the τ^th basis function as Aτ:

\int a_{τ} (\vec{r}) d \vec{r} = A_{τ} .

(3)

Equation (2) becomes:

\int μ (\vec{r}, E)) d \vec{r} = \sum_{τ = 1}^{Γ} A_{τ} ϕ_{τ} (E) .

(4)

Consider a spectral CT image obtained using a photon counting detector with N_E independent energy windows. Let ni (i=1,2,...,N_E) denote the number of photons detected in the ith energy channel with ${\bar{n}}_{i}$ as its mean value. Assuming Poisson statistics for the counting measurements, the probability is:

Pr (n_{1}, \dots, n_{N_{E}}) = \prod_{i = 1}^{N_{E}} \frac{{\bar{n}}_{i}^{n_{i}}}{n_{i}!} e^{- {\bar{n}}_{i}} .

(5)

According to Beer’s law:

\begin{array}{l} {\bar{n}}_{i} & = \int_{E_{i}}^{E_{i + 1}} d E^{'} \int_{0}^{+ \infty} I_{0} (E) \exp (- \sum_{τ = 1}^{n} A_{τ} ϕ_{τ} (E)) h (E^{'}, E) d E \\ = \int_{0}^{+ \infty} I_{0} (E) \exp (- \sum_{τ = 1}^{n} A_{τ} ϕ_{τ} (E)) d E \int_{E_{i}}^{E_{i + 1}} h (E^{'}, E) d E^{'} \\ = \int_{0}^{+ \infty} I_{0} (E) \exp (- \sum_{τ = 1}^{n} A_{τ} ϕ_{τ} (E)) H (i, E) d E, \end{array}

(6)

where I₀(E) denotes the spectrum of the X-ray source and Ei is the energy threshold. h(E’,E) denotes the detector response function. We also define $H (i, E) = \int t_{E_{i}}^{E_{i + 1}} h (E^{'}, E) d E^{'}$ as the integrated detector response for the ith energy channel. We define the vector A=(A₁,...,AΓ)^T. Using the negative log-likelihood and ignoring the constant term, the objective function to be minimized is:

\begin{array}{l} Φ (A_{1}, \dots, A_{Γ}) = \sum_{i = 1}^{N_{E}} ({\bar{n}}_{i} - n_{i} ln {\bar{n}}_{i}) \\ = \sum_{i = 1}^{N_{E}} (\int_{0}^{+ \infty} I_{0} (E) e^{- \sum_{τ = 1}^{Γ} A_{τ} Φ_{τ} (E)} H (i, E) d E \\ - n_{i} ln \int_{0}^{+ \infty} I_{0} (E) e^{- \sum_{τ = 1}^{Γ} A_{τ} Φ_{τ} (E)} H (i, E) d E) \end{array}

(7)

The function Φ(A₁,...,AΓ) is nonlinear. The ni measurements are used to calculate the line integrals A by minimizing Eq.(7) in the case of Γ NE:

({\hat{A}}_{1}, \dots, {\hat{A}}_{Γ}) = \underset{A_{τ}}{\arg \min} Φ (A_{1}, \dots, A_{Γ}) .

(8)

Newton’s method can be used to iteratively solve Eq.(8). According to Eqs.(5-7):

T_{τ} = \frac{\partial Φ}{\partial A_{τ}} = \sum_{i = 1}^{N_{E}} \frac{\partial {\bar{n}}_{i}}{\partial A_{τ}} - \frac{n_{i}}{{\bar{n}}_{i}} \frac{\partial {\bar{n}}_{i}}{\partial A_{τ}} .

(9)

Therefore, the updating formula is:

| \begin{matrix} A_{1}^{k + 1} \\ ⋮ \\ A_{Γ}^{k + 1} \end{matrix} | = {| \begin{matrix} A_{1}^{k} \\ ⋮ \\ A_{Γ}^{k} \end{matrix} | - | \begin{matrix} \frac{\partial T_{1}}{\partial A_{1}} & \dots & \frac{\partial T_{1}}{\partial A_{Γ}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial T_{Γ}}{\partial A_{1}} & \dots & \frac{\partial T_{Γ}}{\partial A_{Γ}} \end{matrix} | |}_{A_{τ}^{k}}^{- 1} {| \begin{matrix} \frac{\partial Φ}{\partial A_{1}} \\ ⋮ \\ \frac{\partial Φ}{\partial A_{Γ}} \end{matrix} | |}_{A_{τ}^{k}} .

(10)

2.2 Optimizing energy thresholds for spectral CT imaging

The line integral of decomposition coefficient Aτ is further reconstructed to determine its spatial distribution $a_{τ} (\vec{r})$ via commonly used CT reconstruction methods. Hence, the accuracy and noise of Aτ will propagate to $a_{τ} (\vec{r})$ , which will affect the output image quality. The key factors in Eq. (7) affecting the accuracy and noise of Aτ are the energy thresholds Ei when I₀(E) is fixed for the system. Hence, this study focuses on finding the optimal energy thresholds.

The variance is minimized to optimize the energy thresholds. According to references [12, 17, 18], there is a lower bound for the variance of the estimate Âτ:

σ_{τ}^{2} \geq {[\begin{matrix} F^{- 1} \end{matrix}]}_{τ, τ},

(11)

where 𝑭 is the Fisher information matrix of the function Φ(A₁,...,AΓ) and [·]_τ,τ represents an element on the diagonal.

The Fisher information matrix with respect to Aτ is calculated as:

F_{τ_{1}, τ_{2}} = {〈 \frac{\partial^{2} Φ}{\partial A_{τ_{1}} \partial A_{τ_{2}}} 〉 |}_{{\hat{A}}_{τ_{1}}, {\hat{A}}_{τ_{2}}},

(12)

\begin{array}{l} \frac{\partial^{2} Φ}{\partial A_{τ_{1}} \partial A_{τ_{2}}} = \\ \sum_{i = 1}^{N_{E}} (\frac{\partial^{2} {\bar{n}}_{i}}{\partial A_{τ_{1}} \partial A_{τ_{2}}} + \frac{n_{i}}{{\bar{n}}_{i}^{2}} \frac{\partial {\bar{n}}_{i}}{\partial A_{τ_{1}}} \frac{\partial {\bar{n}}_{i}}{\partial A_{τ_{2}}} - \frac{n_{i}}{{\bar{n}}_{i}} \frac{\partial^{2} {\bar{n}}_{i}}{\partial A_{τ_{1}} \partial A_{τ_{2}}}), \end{array}

(13)

where 〈·〉 denotes the expectation value. With

〈 n_{i} 〉 = {\bar{n}}_{i}

(14)

Eqs.(12-13) can be simplified to:

{F_{τ_{1}, τ_{2}} = \sum_{i = 1}^{N_{E}} \frac{1}{{\bar{n}}_{i}} \frac{\partial {\bar{n}}_{i}}{\partial A_{τ_{1}}} \frac{\partial {\bar{n}}_{i}}{\partial A_{τ_{2}}} |}_{{\hat{A}}_{τ_{1}}, {\hat{A}}_{τ_{2}}} .

(15)

The optimal energy thresholds are then:

E_{i}^{*} = \underset{{E_{i}}, i \in [1, N_{E}]}{\arg \min} {[F^{- 1}]}_{τ, τ} .

(16)

However, there are three problems:

1) The line integral Aτ is not known for practical scans.

2) There are a huge number of ray paths in CT images, and the optimal energy thresholds can differ for different paths.

3) The optimal energy thresholds can also differ for different basis components.

In light of these problems, the following method was developed to optimize contrast agent imaging in spectral CT images:

1) Pre-scan an object with an ultralow dose.

2) Estimate Âτ from the pre-scan data using Eqs. (8-10).

3) Find the optimal energy thresholds for each ray for the basis function of interest, which are the contrast fractions in this task according to Eqs. (15-16) (the basis function of interest is the component of contrast agent in this paper).

4) Set the optimal energy thresholds as the thresholds with maximum frequency for the ray paths of interest.

5) Scan the object with a normal dose.

6) Decompose the material components and reconstruct the spatial distribution of the coefficients.

Since there is no closed form solution for the noise as a function of the energy thresholds, we can only optimize by threshold scan.

The pseudo-code for the decomposition procedure and the procedure used to find the optimal thresholds for every single ray path is as follows (Table 1).

The pseudo-code for the optimization method

Define N_iter and initialize A1
for iteration k=1 to N_iter
	for i=1 to NE, Calculate ${\bar{n}}_{i}^{k}$ by Eq. (6)
	end
	Update A^k+1 by Eq. (8-10)
end
$\begin{matrix} A = A^{N_{iter} + 1} \end{matrix}$
To optimize Ei, fix $E_{1}, \dots, E_{i - 1}, E_{i + 1}, \dots, E_{N_{E}}$
Calculate ${\bar{n}}_{1}, \dots, {\bar{n}}_{i - 1}, {\bar{n}}_{i + 1}, \dots, {\bar{n}}_{N_{E}}$ by Eq. (6)
for Ei=E_i-1 to E_i+1:
	Calculate ${\bar{n}}_{i}$ by Eq. (6)
	Calculate F by Eq.(15)
	Noise(Ei)=[F^-1]_τ,τ
end
Search $E_{i}^{*} = {\arg \min}_{{E_{i}}} Noise {(E_{i})}_{τ, τ}$

3 Experimental studies

The method was validated using simulations and practical experiments on commonly used I and Gd solutions. These experiments used the photoelectric effect function, Compton effect function, and attenuation function for the contrast agent as the basis functions for component decomposition.

3.1 Simulation experiments

The simulations used the phantoms shown in Fig. 1. Data was collected in 360 views over 2π using the energy spectrum shown in Fig. 2. This method can be extended to cases with other numbers of energy bins. We would like to focus on 4 energy bins because the advantage of using more than 4 energy bins is limited in our studies. The detector response function was assumed to be ideal, i.e., H(i,E)=1. The distance from the source to the center of rotation is 48.5 cm, and the distance from the source to the detectors is 73.5 cm. The detector is 15.36 cm in width with 100 um bin size. The pixel size of the reconstruction image is 100 um. These geometric parameters are suitable for both simulation and practical experiments.

Fig. 1.

(a) Numerical Phantom used for imaging Gd; (b) Numerical Phantom used for imaging I.

Fig. 2.

Spectrum and energy thresholds for the Gd phantom. The dashed lines are the non-optimal thresholds. The solid lines are the optimal thresholds.

3.1.1 Optimize energy thresholds for imaging the Gd contrast agent

As mentioned in Section 2.2, a pre-scan of the phantom was performed using the non-optimized energy thresholds. Since the K-edge of Gd is 50 keV, one threshold was set to 50 keV. In the normal dose case, I₀ was 210⁶. To yield reasonable comparisons, the normal dose was divided into two parts with 10% of the normal dose used for the pre-scan so that the pre-scan dose was 0.1I_0=210⁵. The other 90% of the normal dose was used for imaging with the optimal thresholds. The pre-scan thresholds and the spectrum are shown in Fig. 2. The simulations used E1=25 keV, E3=50 keV, and E5=100 keV, and the thresholds for E2 and E4 were determined for each ray.

Eqs.(8-10) were used to compute the line integrals of the decomposition coefficients for each basis function using the pre-scan data. The distribution of the line integral A_Gd is shown in Fig. 3.

Fig. 3.

The distribution for the line integral of Gd component

We could see that there are some line integrals of Gd smaller than 0 because of noise. In theory, zero values of A_Gd are expected for rays that avoid the contrast agent. These negative values are considered to result from noise. We set a threshold to pick A_Gd related to the contrast agent as:

Thr = - min (A_{Gd}) \forall A_{Gd} .

(17)

Therefore, only ray paths satisfying the condition A_Gd gt;Thr are taken into account for energy threshold optimization. This strategy was applied to all our simulations and practical experiments.

Then, the CRLB of the basis functions were calculated for Gd using Eqs.(11-15). We could then obtain the thresholds of E2 and E4, which responds to the minimum CRLB value for each ray path taken into account. The frequencies for the optimal thresholds of E2 and E4 for each path were taken into account, as shown in Fig. 4. The maximum frequencies for E2 and E4 were 39 keV and 60 keV. The optimized thresholds are shown in Fig. 2 for this imaging case. With these thresholds, the number of emitted photons in each energy bin from low to high were 4.361910⁵, 4.059310⁵, 4.021610⁵, and 5.259510⁵. These photon counts in each energy bin differ by about 30%. The phantom was then scanned with the optimal thresholds and an incident photon flux of 0.9I_0=1.810⁶. The results are compared with a normal scan with a non-optimal energy threshold in Fig. 5. Newton’s method was used for decomposition with the FBP reconstruction for each basis function. The results show that the image quality is significantly improved in all the coefficient images. The results were validated quantitatively by calculating the standard deviation (std) over the regions of interest (ROI) indicated in Fig. 5. These ROIs were all circular regions containing a uniform material with the edge pixels excluded. The standard deviations are listed in Table 2, where "PE" denotes the basis component of the photoelectric effect function, and "Cmp" denotes the basis component of the Compton effect function. The results show that the decomposed images of the Gd contrast agent obtained with this method contain less noise, even though only they used only 90% of the dose of the non-optimal energy thresholds. The noise with the photoelectric effect and Compton Effect functions are also lower than the non-optimal threshold case.

Noise of decomposed component image for every basis function.

Material	Component	SD
		Non-optimal	Optimal
50 mg/mL Gd	Gd	3.14e-3	1.44e-3
100 mg/mL Gd		3.22e-3	1.60e-3
150 mg/mL Gd		3.34e-3	1.67e-3
Bone	PE	1.60e3	1.42e3
	Cmp	2.20e-2	1.61e-2

Fig. 4.

Frequency of the optimal values for E2 and E4 for the Gd phantom

Fig. 5.

Decomposed component images of the three basis functions. The left panel shows the non-optimal energy threshold case with 100% dose, and the right panel shows the optimal energy threshold image with 90% dose. (a) Component images of Gd (display window: [-0.015,0.015]). (b) Component image of the photoelectric effect (display window: [-6000,12000]). (c) Component image of the Compton effect (display window: [-0.1 0.12]).

We further tested our method on a more complex FORBILD head phantom[19]. Four Gd solutions were inserted, as shown in Fig 6. All the percentages refer to mass fractions. The settings were the same as before. The pre-scan used a dose of 0.1I_0=110⁵. The source spectrum with pre-scan thresholds are shown in Fig. 7. E2 and E4 are optimized for each ray and the distributions of the resulting optimal values are shown in Fig. 8. E2 and E4 were set to 39 keV and 60 keV (with maximum frequencies) by frequency voting. The decomposed component images for the non-optimal energy thresholds cases were reconstructed using data with a 100% dose while the optimal energy thresholds cases used a 90% dose. The reconstruction results are shown in Fig 9. Standard deviations within the ROIs were again calculated to evaluate the reconstruction image quality (listed in Table 3). The results show that this method is much less noisy, even though only 90% of the dose was used to determine the non-optimal energy thresholds.

Noise in the decomposed component images for each basis function with the FORBILD phantom

Material	Component	SD
		Non-optimal	Optimal
1% Gd	Gd	8.13e-3	2.74e-3
2% Gd		6.33e-3	2.35e-3
3% Gd		5.88e-3	2.45e-3
4% Gd		5.67e-3	2.35e-3
Bone	PE	6.79e3	6.53e3
	Cmp	2.28e-2	2.25e-2

Fig. 6.

FORBILD head phantom with contrast agent added

Fig. 7.

Spectrum and energy thresholds for the FORBILD phantom. The dashed lines are the non-optimal thresholds. The solid lines are the optimal thresholds.

Fig. 8.

The frequency of the optimal values for E2 and E4.

Fig. 9.

Decomposed component images of the three basis functions for the FORBILD phantom. The left panel shows the non-optimal energy threshold case with 100% dose, and the right panel shows the optimal energy threshold result with 90% dose. (a) Component images of Gd (display window: [-0.03,0.1]). (b) Component image of the photoelectric effect (display window: [-20000,90000]). (c) Component image of the Compton effect (display window: [-0.08 0.38]).

3.1.2 Optimized energy thresholds for imaging contrast agent I

In this section, contrast agent I was the target material; hence, I was used as a basis function. The phantom is shown in Fig.1 (b). Since the K-edge of I is 33 keV and the minimum threshold is 25 keV, the analysis did not use another threshold between those two limits, as was the case with Gd. The analysis fixed E1=25 keV, E2=33 keV, and E5=100 keV, and only E3 and E4 were optimized. The pre-scan thresholds are shown in Fig. 10. The pre-scan dose was again 0.1I_0=210⁵. The optimal thresholds were then computed for each ray. The results are shown in Fig. 11. The optimal thresholds are 39 keV and 52 keV for E3 and E4, respectively. The number of emitted photons in the four energy bins are 2.156210⁵, 4.854810⁵, 2.073010⁵, and 8.609410⁵, which is a rather unbalanced division of photon counts. A spectral CT scan was then performed with a dose of 0.9I_0=1.810⁶ using the optimized threshold shown in Fig. 10. The result is compared with the non-optimal threshold case with a dose of I_0=210⁶. The reconstructed material component images shown in Fig. 12 show visually lower noise levels in the reconstructions using these methods with 90% dose. The standard deviations are listed in Table 4. Similar to the Gd contrast agent imaging, the image quality with the decomposed coefficients is significantly improved, though not as much as in the Gd case.

Standard deviation of the decomposed component images for each basis function

Material	Component	SD
		Non-optimal	Optimal
10 mg/mL I	I	1.73e-3	1.57e-3
20 mg/mL I		1.84e-3	1.68e-3
30 mg/mL I		2.03e-3	1.88e-3
Bone	PE	1.85e3	1.75e3
	Cmp	1.96e-2	1.77e-2

Fig. 10.

Spectrum and energy thresholds for the I phantom. The dashed lines indicate the non-optimal energy thresholds. The solid lines are the optimized energy thresholds.

Fig. 11.

The frequency of the optimal value for E3 and E4 for the I phantom

Fig. 12.

Decomposed component images of each basis function for different energy thresholds. The left side shows the non-optimal energy thresholds with 100% (normal) dose, while the right side shows the optimal energy thresholds with 90% dose. (a) Component images for I (display window: [-0.008 0.011]). (b) Component images for the photoelectric effect (display window: [-5000 12000]). (c) Component images for the Compton effect (display window: [-0.06 0.11]).

These three simulations show that energy threshold optimization is an important issue in spectral CT imaging. Moreover, the optimal energy thresholds vary greatly for different imaging tasks. Thus, real applications require case-by-case optimization of the energy thresholds.

3.2 Practical experiments

The energy threshold optimization method was further validated by imaging Gd and I solutions. All the solution percentages indicated here refer to mass fractions.

The measurements used an XCounter Flite X1 photon counting detector. Its energy response function h(E’,E) was calibrated with fluorescence signals using a hybrid model combining Monte Carlo predictions with a linear correction[20, 21]. The X-ray source was 100 kvp, and the current was 1 mA. The source spectrum was then measured using the pre-calibrated detector response function by a blank scan with threshold scanning [22]:

n_{i}^{Blank} = \int_{0}^{+ \infty} I_{0} (E) \int_{E_{i}}^{E_{i + 1}} h (E^{'}, E) d E^{'} d E,

(18)

where Ei are the discretized thresholds with small steps. This equation was calculated numerically.

Since the detector response is calibrated to be a discretized function with step-size 1.8 keV, the resolution of the optimized energy threshold was no better than 1.8 keV.

The exposure time per view for the pre-scan was 0.33 second with data collected for 72 views over 2π during the pre-scan. The normal dose imaging used a 1 s exposure time per view, with projections using 360 views over 2π. The optimal energy threshold scans with the reduced dose used projections with 360 views over 2π and 0.93 second exposure time per view. This ensures that the total pre-scan dose plus the optimal energy threshold scan were nearly equal to the dose used in normal imaging. The analyses used the same Newton’s method decomposition and FBP reconstruction as in the simulations.

3.2.1 Optimized energy thresholds for the Gd contrast agent imaging

The measurements used the phantom shown in Fig. 13. The materials that were placed in the phantom are listed in Table 5.

Materials in the phantom

Region ID	Material
1	1.5% Gd
2	1% Gd
3	0.5% Gd
4	0.4% Gd
5	Water
6	10% CaCl₂
7	Air
8	PMMA

Fig. 13.

Phantom used in the measurements

The measured data was analyzed using the photoelectric effect function, Compton effect function, and attenuation function of 50μm Gd as basis functions. The non-optimal pre-scan energy thresholds shown in Fig. 14 were used to optimize the energy thresholds for each ray. The optimal E2 and E4 frequencies are shown in Fig. 15.

Fig. 14.

Spectrum and thresholds for the Gd phantom. The dashed lines are the non-optimal energy thresholds used in the pre-scan. The solid lines are the optimal energy thresholds.

Fig. 15.

Frequency of the optimal values for E2 and E4

The optimal energy thresholds shown in Fig. 14 were then used to record reduced dose scans. The ring-artifact reduction (RAR) method was used during reconstruction due to strong non-uniformities in the detector response. The decomposition coefficient images are shown in Fig. 16. The noise is significantly reduced, both with and without RAR. The improvement is further confirmed by the ROI standard deviations listed in Table 6.

Standard deviations of the decomposed component images for the various basis functions

Material	Component	SD w/o correction		SD w/o correction
		Non-optimal	Optimal	Non-optimal	Optimal
0.4% Gd	Gd	7.61e-3	7.16e-3	7.59e-3	7.14e-3
0.5% Gd		7.78e-3	7.33e-3	7.82e-3	7.31e-3
1% Gd		8.18e-3	7.62e-3	8.19e-3	7.59e-3
1.5% Gd		8.22e-3	7.72e-3	8.22e-3	7.65e-3
10% CaCl₂	PE	107.80	107.36	99.32	98.80
Water		95.63	95.16	83.74	82.95
10% CaCl₂	Cmp	2.78e-3	2.60e-3	2.62e-3	2.59e-3
Water		2.72e-3	2.50e-3	2.59e-3	2.44e-3

Fig. 16.

Reconstruction results of the decomposed components. The images in the left two columns are without RAR for the non-optimal and optimal thresholds. The images in the right two columns are shown with the RAR corrections. First row: Images of the Gd contrast agent after threshold optimization. Second Row: Images of the photoelectric effect component. Third row: Images of the Compton effect component. Images from the optimal-threshold scan have lower noise than images obtained from the non-optimal-threshold scan.

As in the simulation study, the optimized energy thresholds significantly improve the image quality of the Gd-decomposed image, which is the optimization target. The image qualities with the photoelectric effect components and the Compton effect components are also slightly better than the images obtained with the non-optimized thresholds and with a normal dose. The gains due to energy threshold optimization with and without RAR are similar. Therefore, only the results with RAR from the other experiment will be shown.

3.2.2 Optimize energy thresholds for contrast agent I imaging

These measurements used the phantom shown in Fig. 17 with the materials listed in Table 7.

The materials in the phantom

Region ID	Material
1	0.5% I
2	1% I
3	1.5% I
4	2% I
5	10% NaCl
6	10% CaCl₂
7	Water
8	PMMA

Fig. 17.

Phantom used in the measurements

The analyses used the photoelectric effect, Compton effect, and I mass attenuation functions as basis functions. The pre-scan thresholds and spectrum are shown in Fig. 18. The energy thresholds E1, E2, and E5 were fixed, while only E3 and E4 were optimized as in the simulations. The optimized energy threshold frequencies are plotted in Fig. 19 for the optimal threshold settings shown in Fig. 18. The reconstructions from the optimal energy and non-optimal energy threshold scans are shown in Fig. 20. The ROI standard deviations for the component images are listed in Table 8. The results show reduced component noise for the I contract agent. The results show that using the photoelectric and Compton effect basis functions also improves the images.

Std of decomposed component images

Material	Component	Std w/ RAR
		Non-optimal	Optimal
2% I	I	2.74e-4	2.34e-4
1.5% I		2.63e-4	2.28e-4
1% I		2.58e-4	2.26e-4
0.5% I		2.54e-4	2.22e-4
10% CaCl₂	PE	250.24	203.87
10% NaCl		239.00	197.04
10% CaCl₂	Cmp	1.35e-3	1.29e-3
10% NaCl		1.37e-3	1.30e-3

Fig. 18.

Spectrum and energy thresholds for the I phantom. The dashed lines indicate the non-optimal energy thresholds, And the solid lines indicate the optimal energy thresholds.

Fig. 19.

Frequency of the optimal value for E3 and E4

Fig. 20.

Decomposed component images of the phantom measurements. The left column shows the non-optimal energy thresholds, while the right column shows the optimal energy thresholds. (a) Component images of I with display window [-0.0025 0.002]. (b) Component images of the photoelectric effect with display window [-1000 3000]. (c) Component images of the Compton effect with display window [-0.005 0.025].

Thus, there is a clear upgrade of the image quality of decomposed I component. The image quality of the photoelectric effect and Compton effect components are also better than the non-optimized results.

3.3 Influence from the accuracy of coefficient projections to optimal energy thresholds

The line integrals of the decomposition coefficients for the material basis functions could be noisy in the low dose pre-scan so that we would obtain a noise-affected distribution for the ray-by-ray optimal energy thresholds. Frequency voting will find the location of the maximum frequency value to remove this noise effect. For demonstration, we compared the optimal energy thresholds results from a low-dose scan and a normal dose scan for the simulation study of I imaging. The dose is 2e5 for the normal-dose situation and 2e4 for the low-dose situation.

The ray-by-ray optimal energy thresholds for E3 and E4 are shown in Figs. 21 and 22, respectively. We can see that the optimal energy thresholds may be different in some ray paths due to the larger amount of noise in low-dose pre-scan. The location of the frequency maximum is the same for both E3 and E4. Hence, the frequency voting strategy can mostly remove the influence of noise from a low-dose scan. This tells us that a low-dose pre-scan is sufficient to find the comprehensive optimal energy thresholds.

Fig. 21.

The normalized frequency of optimal energy threshold E3.

Fig. 22.

The normalized frequency of optimal energy threshold E4.

4 Discussion and Conclusion

This work presents a method for optimizing the energy thresholds in a spectral CT for contrast agent imaging. The optimization problem minimizes noise in the contrast agent components. The method is based on CRLB and basis function decomposition in the projection domain. The basic FBP reconstruction method is used to avoid bias in the spatial CT reconstruction. The optimization is customized for each scan using an ultra-low dose pre-scan. Scans using the optimal thresholds then use slightly reduced doses so that the total dose is no higher than for a normal scan. The method was tested in simulations and practical experiments.

We notice that there would be an unavoidable error in the estimations of h(E’,E) and I₀(E), which are required in our CRLB computation. This could affect our optimal thresholds and reconstructions. Hence, improvements in practical experiments could be affected by accuracy of the estimated h(E’,E) and I₀(E) functions. The reconstructed image qualities are much better. The level of improvement is related to the distance between the optimal and non-optimal energy thresholds. Non-target component images can also be improved by optimizing the energy threshold, though not as much as the targeted components. Thus, energy threshold optimization can improve the image quality, with task-oriented optimization yielding the best results. Moreover, we would like to point out that the optimal energy thresholds in our study do not form uniform photon counts among different energy windows, which is a common choice in the field. We observe a difference of more than 300% between photon counts in the optimal setting. We think this is because noise in the reconstructions is mainly affected by the detected signals, which are affected by the source spectrum, energy bin separation, and scanned object. Uniform incident photon counts for different energy bins does not necessarily yield an optimal noise distribution of the detected signals. Therefore, optimal thresholds are task-dependent and could be very different from an even incident photon distribution in energy bins. We also performed tests on freeing the thresholds at K-edge points and optimizing these thresholds. Without exception, setting thresholds at K-edge points of the contrast agents gives us optimal results. Therefore, we can conclude that optimizing the two energy thresholds in our four energy windows case is sufficient to produce images with improved quality.

This study presents a method for optimizing the energy thresholds in multi-energy spectral CT. This method can be computationally expensive because it calculates the probability distribution of the optimal thresholds for all rays. A lookup table could be developed for the decomposition coefficients of various material compositions for the ray paths based on a large set of clinical data. Further work will consider practical implementations of this method.

References

[1]

X. Wang, D. Meier, K. Taguchi, et al.

Material separation in x-ray CT with energy resolved photon-counting detectors

. Med Phys. 38:1534-1546(2011). doi: 10.1118/1.3553401

Baidu Scholar

Google Scholar

[2]

R.E. Alvarez, A. Macovski.

Energy-selective reconstructions in X-ray computerized tomography

. Phys Med Biol. 21:733-744(1976). doi: 10.1088/0031-9155/21/5/002

Baidu Scholar

Google Scholar

[3]

A. So, J. Hsieh, S. Narayanan, et al.

Dual-energy CT and its potential use for quantitative myocardial CT perfusion

. J Cardiovasc Comput Tomogr. 6:308-317(2012). doi: 10.1016/j.jcct.2012.07.002

Baidu Scholar

Google Scholar

[4]

D. Panetta.

Advances in X-ray detectors for clinical and preclinical Computed Tomography

. Nucl Instrum. Meth. A. 809:2-12(2016). doi: 10.1016/j.nima.2015.10.034

Baidu Scholar

Google Scholar

[5]

E. Roessl, R. Proksa.

K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors

. Phys Med Biol. 52:4679-4696(2007). doi: 10.1088/0031-9155/52/15/020

Baidu Scholar

Google Scholar

[6]

J.P. Schlomka, E. Roessl, R. Dorscheid, et al.

Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography

. Phys Med Biol. 53:4031-4047(2008). doi: 10.1088/0031-9155/53/15/002

Baidu Scholar

Google Scholar

[7]

K. Taguchi, J.S. Iwanczyk.

Single photon counting x-ray detectors in medical imaging

. Med Phys. 40:100901(2013). doi: 10.1118/1.4820371

Baidu Scholar

Google Scholar

[8]

L. Shen, Y. Xing, L. Zhang,et al.

Hybrid decomposition method for dual energy CT

. In: 2014 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC). 2014:1-4. doi: 10.1109/NSSMIC.2014.7430897

Baidu Scholar

Google Scholar

[9]

T. Gilat Schmidt,E.Y. Sidky.

Investigation of a one-step spectral CT reconstruction algorithm for direct inversion into basis material images

. In: SPIE Medical Imaging.; 2015:94124Y. doi: 10.1117/12.2082057

Baidu Scholar

Google Scholar

[10]

K. Taguchi,M. Zhang, E.C. Frey, et al.

Image-domain material decomposition using photon-counting CT

. Proc SPIE - Int Soc Opt Eng.:651008(2007). doi: 10.1117/12.713508

Baidu Scholar

Google Scholar

[11]

H. Cramer. Mathematical Methods of Statistics.(Princeton: Princeton university press, 1945)

[12]

E. Roessl,C. Herrmann.

Cramer-Rao lower bound of basis image noise in multiple-energy x-ray imaging

. Phys. Med. Biol. 54:1307-1318(2009). doi: 10.1088/0031-9155/54/5/014

Baidu Scholar

Google Scholar

[13]

E. Roessl, B. Brendel,K.J. Engel, et al.

Sensitivity of photon-counting based K-edge imaging in X-ray computed tomography

. IEEE Trans Med Imaging. 30(9):1678-1690(2011).

Baidu Scholar

Google Scholar

[14]

E. Roessl, R. Proksa.

Optimal energy threshold arrangement in photon-counting spectral x-ray imaging

. In: IEEE Nucl. Sci. Conf. R 3; 2007:1950-1954. doi: 10.1109/NSSMIC.2006.354276

Baidu Scholar

Google Scholar

[15]

B. Meng, W. Cong, Y. Xi, et al.

Energy window optimization for X-Ray K-edge tomographic imaging

. IEEE Trans Biomed Eng. 63:1623-1630(2016). doi: 10.1109/TBME.2015.2413816

Baidu Scholar

Google Scholar

[16]

Y. Xing, Y. Li, L. Shen.

A General Adaptive Decomposition method for Multi-Energy Spectral CT

, Proc.2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC) (2013).

Baidu Scholar

Google Scholar

[17]

S. K. Sengupta, Fundamentals of Statistical Signal Processing: Estimation Theory. PTR Prentice hall, (1993).

[18]

W. T. Eadie, L. D. Kirkpatrick, Statistical Methods in Experimental Physics. North-Holland, (1982).

[19]

Arbeitsgruppe Phantome

[EB/OL] 2013-05-02. http://www.imp.uni-erlangen.de/phantoms/index.html

Baidu Scholar

Google Scholar

[20]

D. Wu,X. Xu, L. Zhang, et al.,

A hybrid Monte Carlo model for the energy response functions of X-ray photon counting detectors

. Nucl. Instrum. Methods. 830:397-406(2016) doi: 10.1016/j.nima.2016.05.097

Baidu Scholar

Google Scholar

[21]

X. Xu, L. Zhang,D. Wu,et al.,

Response function estimation for the XCounter Flite X1 photon counting detector using Monte Carlo method

. In: Nuclear Science Symposium, Medical Imaging Conference and Room-Temperature Semiconductor Detector Workshop.(2017).

Baidu Scholar

Google Scholar

[22]

S. Wang, H. Gao, L. Zhang,et al..

Quasi-monochromatic imaging in x-ray CT via spectral deconvolution using photon-counting detectors

. Phys Med Biol. 62:2208-2223(2017). doi: 10.1088/1361-6560/aa5a47

Baidu Scholar

Google Scholar