The small animal material discrimination study based on equivalent monochromatic energy projection decomposition method with dual-energy CT system

Abstract

Material discrimination is an important application of dual-energy computed tomography (CT) techniques. Projection decomposition is a key problem for pre-reconstruction material discrimination. In this study, we focused on the pre-reconstruction space based on the photoelectric and Compton effect decomposition model to characterize different material components, and proposed an efficient method to calculate the projection decomposition coefficient. We converted the complex projection integral into a linear equation by calculating the equivalent monochromatic energy from the high and low energy spectrum. Meanwhile, we constructed a dual-energy CT system based on a photon-counting detector to take small animal scan and material discrimination analysis. Finally, the results of simulation and experimental study demonstrated the feasibility of our proposed new method, and explained the characteristics of photoelectric absorption and Compton scattering reconstruction images.

Keywords

Keywords: Dual-energy CT photon-counting detector material discrimination projection decomposition

1 Introduction

X-ray computed tomography (CT) have been widely applied in medical and industrial applications [1]. Conventional CT systems employ the digital integrating detector whose output is proportional to the energy integrated over the entire incidence spectrum [2]. With the development of energy-solved photon-counting detector [3 –5], absorption features in the multiple ranges of photon energy can be identified [6, 7], which has facilitated dual-energy CT to distinguish different materials. Dual-energy data acquisition by photon-counting detector has two main approaches: switching different tube voltages of the x-ray source and setting different energy thresholds of the photon-counting detector [8, 9]. The two sets of spectral measurements provide sufficient absorption features to decompose materials and characterize components. Thus, two basis materials such as bone and soft tissue can be discriminated and represented by their physical characteristics in decomposition images.

Current material decomposition methods based on different x-ray spectrum measurements include projection-based decomposition of pre-reconstruction and image-based decomposition of post-reconstruction [10 –14]. Alvarez and Macovski first proposed the attenuation coefficient function based on dual-energy measurements to separate the photoelectric and Compton contributions in the projection domain [10]. E Roessl and R Proksa introduced K-edge characteristics of the contrast agent and the maximum-likelihood estimator to decompose sinogram data, and represented a nonlinear optimization problem [15]. Zhang L. et al. set up a pre-created lookup table to match the basis material coefficients via a nonlinear integral model with x-ray high-energy spectrum (more than 1 MeV) [16]. Heismann et al. developed a material decomposition method based on post-reconstruction space, which used density and atomic number as functions of attenuation values [13]. Liu X. et al. introduced the principle of mass conservation as a third criteria to achieve three-material decomposition with basis material method in image domain [17], and the accuracy of this method is limited on the elements with low mass fractions. Mendonca et al. established a material library containing basis material triplets to derive an image-based triple material decomposition [18], which obtained linear-attenuation dual-material-density pair from dual-energy measurements and applied volume conservation assumption to decompose materials pixel by pixel.

The accuracy of material decomposition for post-reconstruction method depends on raw data and reconstruction algorithm [19], which is easily influenced by noises and artifacts in image domain. Pre-reconstruction method is performed in projection domain, which generally accompanies some complicated spectrum characteristics analysis and calculation. Thus, this paper proposed an efficient method to calculate the projection decomposition coefficient. The proposed method introduced an equivalent monochromatic energy to simplify projection decomposition calculation, which could generate an effective and accurate decomposition for different material components with two spectrum measurements. In the next section, the material decomposition algorithm is described. In the third section, corresponding experimental results are presented. In the last section, we will discuss the related issues and conclude the paper.

2 Materials and methods

In medical CT imaging, the diagnostic energy range is within 140 keV. Photoelectric absorption and Compton scattering are the two dominant x-ray attenuation processes [20]. Hence the linear attenuation coefficient can be decomposed into a linear combination of photoelectric absorption and Compton scattering: $μ (E) = (σ_{ph} + σ_{co}) ρ N_{A} / A_{r}$ (1)

where ρ, N_A and A_r are mass density, Avogadro’s number and atomic mass respectively. The photoelectric absorption and Compton scattering cross section σ_ph, σ_co represent the attenuation probability of incident photons and objects [21], which are formulated as: $σ_{ph} = k_{p} Z^{4} / E^{3}$ (2) $σ_{co} = k_{c} Zf (ɛ)$ (3)

Here, k_p, k_c are the constants of photoelectric absorption and Compton scattering, Z is the atomic number, and f (ɛ) is the Klein-Nishina function: $f (ɛ) = \frac{1 + ɛ}{ɛ^{2}} [\frac{2 (1 + ɛ)}{1 + 2 ɛ} - \frac{ln (1 + 2 ɛ)}{ɛ}] + \frac{ln (1 + 2 ɛ)}{2 ɛ} - \frac{1 + 3 ɛ}{{(1 + 2 ɛ)}^{2}}$ (4) $ɛ = \frac{hv}{m_{e} c^{2}} \approx \frac{E}{511 keV}$ (5)

When x-ray source emits a continuous polychromatic spectrum, the x-ray intensity I depends on the object material composition and the photon energy after x-ray beam passes through an object. At a given energy threshold, the x-ray intensity I measured by the detector is represented as: $I (E) = \int_{E_{\min}}^{E_{\max}} S (E) e^{- \int μ (E, r) d r} d E$ (6)

where S (E) is the detector response spectrum related to energy distribution of x-ray source.

Assuming f_p (E) = 1/E³, f_c (E) = f (ɛ) and substituting Equations (1–3) to Equation (6), the above expression can be inverted into:

$I (E) = \int_{E_{\min}}^{E_{\max}} S (E) e^{- \int [n_{p} f_{p} (E) + n_{c} f_{c} (E)] d r} d E$ (7)

The two new variables n_p, n_c are the material dependent coefficients of photoelectric absorption and Compton scattering respectively, which are related to the physical characteristics of materials such as mass density and atomic number. The Equation (7) denotes the integral of energy collected over the whole given x-ray spectrum, and the integrand is only related to energy which is expressed as: $i (E) = S (E) e^{- N_{p} f_{p} (E) - N_{c} f_{c} (E)}$ (8)

where $N_{p} = \int n_{p} d r$ , $N_{c} = \int n_{c} d r$ . According to mean value theorem for integrals, the integrating measurement of x-rays intensity can be transformed to: $I (\tilde{E}) = i (\tilde{E}) (E_{\max} - E_{\min})$ (9)

Here, we introduce the equivalent monochromatic energy $\tilde{E}$ to indicate the polychromatic spectrum between E_min and E_max. In terms of this mathematical processing, the complex projection integral can be simplified into the product of the equivalent integrand function value $i (\tilde{E})$ and the length of integral interval.

Dual-energy imaging CT system with energy-resolved photon-counting detector can acquire two spectrally distinct measurements I_l, I_h by setting two different energy thresholds. The x-ray intensities I_l, I_h with low- and high-energy spectrum are described by the following equations: ${\begin{matrix} I_{l} = [S ({\tilde{E}}_{l}) e^{- N_{p} f_{p} ({\tilde{E}}_{l}) - N_{c} f_{c} ({\tilde{E}}_{l})}] (E_{maxl} - E_{minl}) \\ I_{h} = [S ({\tilde{E}}_{h}) e^{- N_{p} f_{p} ({\tilde{E}}_{h}) - N_{c} f_{c} ({\tilde{E}}_{h})}] (E_{maxh} - E_{minh}) \end{matrix}$ (10)

In Equation (10), E_max, E_min are the high and low energy thresholds, and N_p, N_c are the unknowns needed to be solved. $S (\tilde{E})$ can be directly acquired by detector measurements. Then the equivalent monochromatic energy ${\tilde{E}}_{h}$ and ${\tilde{E}}_{l}$ can be calculated, thereby obtaining $f_{p} (\tilde{E})$ and $f_{c} (\tilde{E})$ . Next by solving equations, we are able to complete projection decomposition and manage to separate materials. The key process of calculating ${\tilde{E}}_{h}$ , ${\tilde{E}}_{l}$ is as follows.

Step 1. Dual-energy CT image reconstruction

Split-Bregman algorithm is used to reconstruct the dual-energy CT images f_l, f_h. CT reconstruction problem based on compressive sensing theory can be described by a mathematical model: $f = \underset{f}{argmin} {∥ Φ (f) ∥}_{1} + λ ∥ mf - b ∥_{2}^{2}$ (11)

where Φ (f) is the gradient transform of reconstructed image f, m is system matrix, b is projection data, and λ is penalty weight. By using an intermediate variable d = Φ (f), Equation (11) can be converted into: $f = \underset{f, d}{argmin} {∥ d ∥}_{1}^{1} + λ {∥ mf - b ∥}_{2}^{2} + γ {∥ d - Φ (f) ∥}_{2}^{2}$ (12)

where γ is the convergence coefficient. Then Equation (12) is converted into an unconstrained optimal problem: $(f^{k + 1}, d^{k + 1}) = \underset{f, d}{argmin} {∥ d ∥}_{1}^{1} + λ {∥ mf - b ∥}_{2}^{2} + γ {∥ d - Φ (f) - b^{k} ∥}_{2}^{2}$ (13) $b^{k + 1} = b^{k} + (Φ (f^{k + 1}) - d^{k + 1})$ (14)

Equation (13) can be solved by splitting into two equations: $f^{k + 1} = \underset{f}{argmin} λ {∥ mf - b ∥}_{2}^{2} + γ {∥ d^{k} - Φ (f) - b^{k} ∥}_{2}^{2}$ (15) $d^{k + 1} = arg min_{d} {∥ d ∥}_{1}^{1} + γ {∥ d - Φ (f^{k + 1}) - b^{k} ∥}_{2}^{2}$ (16)

Spilt-Bregman method can effectively suppress noise and artifacts, and accelerate iterative convergence in CT image reconstruction [22], which facilitates the following material decomposition.

Step 2. Attenuation coefficient calculation

Select the same material component regions from the reconstructed images in low- and high-energy spectrum, and calculate the average gray value of the same component regions. The gray value of material component in reconstructed image represents the x-ray attenuation coefficient, so the x-ray attenuation coefficients μ_l and μ_h in low- and high-energy spectrum are obtained.

Step 3. The equivalent monochromatic energy analysis

Fit the theoretical attenuation curve μ (E) using x-ray attenuation characteristic data reported by the US National Institute of Standards and Technology (NIST) [23] and calculate the equivalent energy ${\tilde{E}}_{l}$ and ${\tilde{E}}_{h}$ according to x-ray attenuation coefficients μ_l and μ_h in low- and high-energy spectrum. Then by solving Equation (10), the photoelectric absorption and Compton scattering effect decomposition coefficients integrals N_p, N_c can be figured out, and n_p, n_c are thereby obtained after reconstruction process in step 1.

In this section, the equivalent monochromatic energy in the two spectral measurements are acquired according to relationship between the material gray value in reconstructed images and the theoretically fitting attenuation curve. We use the equivalent monochromatic energy to simplify projection integrals equations and achieve material decomposition.

3 Results and analysis

3.1 Numerical simulation

To evaluate the feasibility of our proposed method for dual-energy projection decomposition, we designed a multi-material phantom as shown in Fig. 1. The circular phantom has a diameter of 4 cm, and contains different size of regions that are filled with soft tissue and bone materials. It is discretized into a 400×400 matrix (the pixel size: 0.01 cm×0.01 cm). The goal of the multi-material phantom design is to mimic a simple small biological sample (i.e., a mouse), so the background region represents soft tissue and the eight sub-regions represent bone. In Fig. 2. we plotted the linear attenuation coefficient curve of bone and soft tissue according to the x-ray attenuation characteristic database reported by NIST.

Fig.1

Numerical dual-energy phantoms. (a) Image at energy of 30–40 keV, (b) Image at energy of 60–70 keV. Display window: [0 2].

Fig.2

The linear attenuation coefficient of two materials.

In the simulation study, the x-ray spectrum was simulated from 15 keV to 100 keV. We selected two energy bins (30–40 keV and 60–70 keV) to analyze the attenuation characteristics of bone and tissue, and calculated the average linear attenuation coefficients in the two energy bins. According to the attenuation coefficient curve in Fig. 2, we could obtain the equivalent monochromatic energy of the two energy spectrum: ${\tilde{E}}_{l} = 32.8 keV, {\tilde{E}}_{h} = 64.7 keV$ . The phantom projection data could be acquired by Radon transformation. We used the previously mentioned method to decompose projection data, and reconstructed the photoelectric absorption and Compton scattering images, which are shown in Fig. 3.

Fig.3

Reconstructed images through projection decomposition. (a) The reconstructed photoelectric absorption image, display window: [0 0.5] (b) The reconstructed Compton scattering image, display window: [0 1.2].

3.2 Experimental study

Meanwhile, we applied our proposed method to study a mouse specimen. In the experimental study, we constructed a dual-energy CT imaging system with a photon-counting detector to distinguish x-ray attenuation characteristics in different energy bins, as shown in Fig. 4. The system uses a transmitted micro-focus x-ray source, the voltage of which ranges from 20 kVp to 225 kVp and the current ranges from 0.05 mA to 1.0 mA. The detector is the CdTe-CMOS photon counting detector manufactured by X-Counter Corporate of Sweden, and utilizes two energy thresholds to select and record incoming photons. It has 2048×64 pixels with an effective pixel pitch of 100μm, and the x-ray energy detection range is from 15 keV to 250 keV.

Fig.4

The dual-energy CT imaging system based on x-ray photon-counting detector.

The mouse was injected anesthetic and put in the bottom of a plastic bottle which was placed on the scanning rotation table, as shown in Fig. 5. Before the specimen scan, we took the system geometric calibration to determine the detector to axis distance and x-ray source to axis distance [24]. Then we used the dual-energy CT imaging system to scan the chest of mouse. The energy thresholds of detector were set to 15 keV and 60 keV, and we scanned specimen twice to acquire two sets of projection data in low- and high-energy spectrum respectively. The experimental parameters were summarized in Table 1.

Fig.5

The scanned mouse specimen.

Table 1

The experimental parameters

Parameter	Value
Tube voltage	100 kV
Tube current	70μA
Source to detector distance	480 mm
Source to object distance	160 mm
Total number of angles	500
Angle increment	0.72°
Length of detector	2048

Figure 6 shows the projections of the mouse specimen at one angle with three scan positions, in which the bone joint of the mouse could be clearly seen. The detector had 64 layers for one position and we selected the 32nd layer to extract the sinogram for one image slice along the marked line in Fig. 6. Then, we used Split-Bregman method to reconstruct mouse specimen dual-energy CT images at three positions, as shown in Fig. 7. The bone regions are clearly visible in reconstructed images at energy of 15–60 keV, and the images at energy of 60–100 keV contain some noises and lose some details due to the decrease of photon number.

Fig.6

Three projections of the mouse specimen. (a) Projection at position 1, (b) Projection at position 2, (c) Projection at position 3.

Fig.7

Reconstructed CT images at three positions. From left to right is position 1 to positon 3. (a) Images at energy of 15–60 keV, (b) Images at energy of 60–100 keV. Display window: [0 0.5].

Then, we took the bone regions of the mouse as reference to calculate the attenuation coefficients followed by the step 2 described previously. The measurement results of each position with dual-energy CT are shown in Table 2. There are some differences in attenuation coefficients of bone with different energy ranges. From the linear attenuation coefficient curve of bone, we could approximately obtain equivalent monochromatic energy value in low and high energy spectrum at three positions.

Table 2

The measurement results of three positions with dual-energy CT

Scan position	μ _l	E _l	μ _h	E _h
Position 1	0.3707	55 keV	0.2998	64 keV
Position 2	0.3729	54.7 keV	0.3015	63.7 keV
Position 3	0.3665	55.3 keV	0.2954	64.5 keV

After we obtained equivalent monochromatic energy value in two spectrum, the photoelectric and Compton decomposition coefficient linear integrals were calculated using Equation (10). Following the same reconstruction step, we could finally obtain the photoelectric absorption and Compton scattering images respectively. The reconstructed decomposition images are shown in Fig. 8.

From Fig. 8, we can see that two sets of decomposition images have their own distinct features. The Compton scattering image highlights the soft tissue of the mouse, whereas in the photoelectric absorption image, the skeletal portion is clearly visible. The experimental results demonstrated the performance of our proposed material discrimination algorithm, which distinguished the bone and soft tissue of the mouse.

Fig.8

Reconstructed images after projection decomposition. (a) The reconstructed Compton scattering images, display window: [0.01 0.15] (b) The reconstructed photoelectric absorption images, display window: [0.02 0.1].

4 Discussions and conclusions

In this paper, the proposed projection decomposition algorithm based on equivalent monochromatic energy analysis with dual-energy pre-reconstruction is efficient and accurate. We used mean value theorem for integral transform to simplify projection integrals equations, and obtained the equivalent monochromatic energy, which avoided complex calculations of energy-related parameters. In CT image reconstruction, we used Spilt-Bregman method to reduce noises and accelerate iterative convergence, since high-quality reconstruction images are the important prerequisite for material decomposition.

There are some key factors on material discrimination worthy of further discussion. The bone and soft tissue discrimination depend on photoelectric absorption effect and Compton scattering. Due to different material atomic number and incident photon energy, each effects have different physical performance [25]. In general, photoelectric absorption effect easily occurs when the material atomic number is higher and the photon energy is lower. Moreover, the Compton scattering tends to dominate as the incident energy of photon increases [13, 26]. The main component of bones is calcium carbonate, and the muscular tissues are mainly composed of water and protein. Therefore, the overall atomic number of bones is higher and it principally reflects the photoelectric absorption effect. In muscles and other tissues, it is ruled by the Compton scattering. Although the two main components of the specimen are well separated from the reconstructed images, different types of soft tissues could not be further distinguished due to the insufficiency spectrum data provided by two spectrum measurements.

In further research, we consider to take contrast agents and increase the number of energy spectrum to achieve multi-material decomposition and acquire more image information that we are interested in. Moreover, the method described in the study is not only applied to the discrimination of biological specimen, but also to the mixtures with different atomic number materials. We focus on solving the equivalent monochromatic energy of different spectrum to simplify the projection decomposition, which is not limited to the material types but the number of spectrum. The basis material model can also be combined with this method to achieve material discrimination.

In conclusion, we have proposed a projection decomposition algorithm based on pre-reconstruction method, and demonstrated the feasibility of the algorithm with simulation and experimental study. Finally, we analyzed the characteristics of photoelectric absorption effect and Compton scattering to distinguish different components of materials. We will next collect multi-energy spectrum data to facilitate further improvements in material discrimination.

Footnotes

Acknowledgments

This work was supported by National Key R&D Program of China (No. 2016YFC0104609), National Natural Science Foundation of China (No. 61401049), the Strategic Industry Key Generic Technology Innovation Project of Chongqing (No. cstc2015zdcy-ztzxX0002), the Chongqing Foundation and Frontier Research Project (No. cstc2016jcyjA0473), the Fundamental Research Funds for the Central Universities (No. 10611CDJXZ238826), Visiting Scholar Foundation of Key Laboratory of Optoelectronic Technology & Systems (Chongqing University), Ministry of Education.

References

Wang

, Yu

and De Man

, An outlook on X-ray CT research and development, Med Phys 35(3) (2008), 1051–1064.

, Wei

, Feng

, Chen

and Mi

, Material discrimination based on k-edge characteristics, Comput Math Methods Med 2013(3) (2013), 308520.

Tanguay

, Kim

H.K.

and Cunningham

I.A.

, The role of x-ray Swank factor in energy-resolving photon-counting imaging, Med Phys 37(3) (2010), 6205–6211.

Shikhaliev

P.M.

, Photon counting spectral CT: Improved material decomposition with K-edge-filtered x-rays, Phys Med Biol 57(6) (2012), 1595.

Ren

, Zheng

and Liu

, Tutorial on X-ray photon counting detector characterization, J Xray Sci Technol 26(1) (2018), 1–28.

, Yu

, Bennett

, Ronaldson

, Zainon

, Butler

, Wei

and Wang

, Energy-discriminative performance of a spectral micro-CT system, J Xray Sci Technol 21(3) (2013), 335–345.

Wang

, Meier

, Taguchi

, Wagenaar

D.J.

, Patt

B.E.

and Frey

E.C.

, Material separation in x-ray CT with energy resolved photon-counting detectors, Med Phys 38(3) (2011), 1534–1546.

Hao

, Kang

, Zhang

and Chen

, A novel image optimization method for dual-energy computed tomography, Nucl Instrum Methods Phys Res 722(722) (2013), 34–42.

, Chen

, Cong

and Wang

, Spectral ct modeling and reconstruction with hybrid detectors in dynamic-threshold-based counting and integrating modes, IEEE Trans Med Imag 34(3) (2015), 716–728.

10.

Alvarez

R.E.

and Macovski

, Energy-Selective Reconstructions in X-Ray Computerized Tomography, Phys Med Biol 21(5) (1976), 733–744.

11.

Lehmann

L.A.

, Alvarez

R.E.

, Macovski

, Brody

W.R.

, Pelc

N.J.

and Riederer

S.J.

, Generalized image combinations in dual kvp digital radiography, Med Phys 8(5) (1981), 659–667.

12.

Kalender

W.A.

, Perman

W.H.

, Vetter

J.R.

, Mazess

R.B.

and Holden

J.E.

, Evaluation of a prototype dual-energy computed tomographic apparatus. I. Phantom studies, Med Phys 13(3) (1986), 334–339.

13.

Heismann

B.J.

, Leppert

and Stierstorfer

, Density and atomic number measurements with spectral x-ray attenuation method, J Appl Phys 94(3) (2003), 2073–2079.

14.

Tremblay

J.É.

, Bedwani

and Bouchard

, A theoretical comparison of tissue parameter extraction methods for dual energy computed tomography, Med Phys 41(8) (2014), 081905.

15.

Roessl

and Proksa

, K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors, Phys Med Biol 52(15) (2007), 4679–4696.

16.

Zhang

, Chen

, Kang

, Hu

, Xing

, Duan

, Li

, Liu

and Huang

, Image reconstruction method for high-energy, dual-energy CT system, U.S. Patent, (2012).

17.

Liu

, Yu

, Primak

A.N.

and Mccollough

C.H.

, Quantitative imaging of element composition and mass fraction using dual-energy ct: Three-material decomposition, Med Phys 36(5) (2009), 1602–1609.

18.

Mendonça

P.R.S.

, Bhotika

, Maddah

, Thomsen

, Dutta

, Licato

P.E.

, Joshi

M.C.

, Multi-material decomposition of spectral CT images, in, SPIE 2010 Med Imag 7622(3) (2010), 640–643.

19.

Lee

, Choi

Y.N.

and Kim

H.J.

, Quantitative material decomposition using spectral computed tomography, Phys Med Biol 59(18) (2014), 5457–5482.

20.

Schirra

C.O.

, Roessl

, Koehler

, Brendel

, Thran

and Pan

, Statistical reconstruction of material decomposed data in spectral CT, IEEE Trans Med Imag 32(7) (2013), 1249–1257.

21.

Podgoršak

E.B.

, Radiation Physics for Medical Physicists, Springer Berlin Heidelberg, 2006.

22.

Chen

, Mi

, He

, Deng

and Wei

, A CT reconstruction algorithm based on L1/2 regularization. Comput Math Methods Med 2014(3) (2014), 862910.

23.

http://www.nist.gov/pml/data/xraycoef/index.cfm.

24.

, Yu

, Thayer

, Jin

, Xu

, Bennett

, Tappenden

, Wei

, Goldstein

, Renaud

, Butler

and Wang

, Preliminary experimental results from a MARS micro-CT system, J Xray Sci Technol 20(2012), 199–211.

25.

Cullen

D.E.

, A simple model of photon transport, Nucl Instrum Methods Phys Res 101(4) (1995), 499–510.

26.

Brown

J.M.C.

, Dimmock

M.R.

, Gillam

J.E.

and Paganin

D.M.

, A low energy bound atomic electron compton scattering model for geant4, Nucl Instrum Methods Phys Res 338(1) (2014), 77–88.