MISD-IR: material-image subspace decomposition-based iterative reconstruction with spectrum estimation for dual-energy computed tomography

Introduction

Dual-energy computed tomography (DECT) is a promising technique, which has been increasingly applied in medical imaging (1), security inspection (2), and industrial production (3). Compared with traditional computed tomography (CT), DECT can obtain quantitative properties of different materials, which can be exploited by using the attenuation information under different spectra called material decomposition (MD) (4). However, it is difficult to obtain accurate material maps due to noise (5), scatter (6), beam-hardening (7), and so on (8).

Research on MD can be divided into three categories including image-domain-based reconstruction, projection-domain-based, and one-step iterative reconstruction (9,10). Since the image-domain-based reconstruction carries out the decomposition on reconstructed CT images, it is convenient in clinical application (11-14). The main problems of image-domain methods are the magnification of both noise and beam-hardening in material maps from reconstructed images (10). A common way to deal with the noise is to introduce prior information (5) of material maps. However, the imaging model is assumed to be linear in these methods, which does not fit with reality. Therefore, in theory, the image-domain decomposition methods cannot completely eliminate hardening artifacts.

Compared with image-domain-based methods, projection domain-based methods and one-step iterative reconstruction directly work on the transmission data. Projection domain-based methods first decompose the materials in the projection domain and then obtain the material maps by classical reconstruction methods such as filter back-projection algorithm (FBP) (15). These kinds of methods (16) can process the projection of each view independently, which can not only help to solve the nonlinear effect but also improve the calculation efficiency. The projection domain-based approach was proposed by Alvarez and Macovski who laid down the basic principles in their pioneering paper published in 1976 (9). This approach was further explored by the Hasegawa group, which developed several DECT prototypes in the 80s and 90s. In particular, the Hasegawa group was concerned also with the systematic bias of DECT arising from the polychromatic X-ray source (17), and how to correct for the systematic bias, for example, using the basis materials transformation method (18), so that the errors propagation when used in conjunction with other image modalities, such as single-photon X-ray CT, will be minimized (19). However, the geometric consistency of different transmission data and the requirement of accurate spectral distribution hinder the application of projection-domain-based methods.

Different from projection domain-based methods, one-step iterative reconstruction methods establish the imaging model based on transmission data and obtain material maps directly. These methods are not only flexible for scanning conditions but also have the potential to eliminate nonlinear effects, which has received increasing attention. However, a primary problem for this kind of method is the estimation of accurate spectral distribution information. Common approaches to estimating the spectrum are based on the physical phantom with known attenuation information and geometry, which is called indirect spectrum estimation. These methods transform spectrum estimation into finding the solutions of linear equations. Research on indirect spectrum estimation mainly focuses on tackling the pathology of linear equations. Designing and optimizing the regularization terms is the main research direction (20). The adoption of expectation-maximization (EM) is a widely applied classical method (21), solving the equation iteratively, and often the early stop is utilized. Recently, Ha et al. (22) proposed to transform energy estimation into an optimization problem with Kullback-Leibler divergence. It is convenient to incorporate prior information in this method and the effectiveness can be comparable to that of the EM method. Another research direction is to use a few parameters to fit the spectrum with a set of basis models. Zhao et al. (23) proposed to use weighted known spectra to fit the estimated spectrum. Since the known spectra are simulated according to the conditions of actual scanning, these spectra are close to the actual spectra. Thus, it is highly stable with the potential for broader application (23,24). However, since energy spectrum estimation and MD are independent of each other, spectral mismatch may seriously affect the decomposition accuracy (25).

Referring to the reconstruction of material maps, many methods have been developed in the past decades. Extended algebraic reconstruction technique (E-ART) (26) and statistical iterative reconstruction are two classical methods (27). To suppress the noise and artifacts, statistical dependencies in channels (28,29), and a semi-empirical forward model (30) have been proposed and incorporated into material image reconstructions. The construction of regularization terms with sparsity is another way to suppress noise and artifacts. For example, Cai et al. (31) proposed a full-spectral Bayesian reconstruction approach with the Huber function as a regularization term. Liu et al. (32) developed a total image-constrained method with non-local total variation (TV) regularization. It utilized the high signal-to-noise ratio (high-SNR) reconstructed CT images as a reference to construct non-local TV as a regularization term achieving impressive results. The similarity between the reconstructed CT images and material maps has also been applied in ring-artifact correction (33) and obtained promising results. Additionally, considering clinical application, convergence speed is an important research direction. In related works, E-ART was modified with the simultaneous algebraic reconstruction technique (ESART) (34), oblique projection modification technique (35), and the monochromatic images guided iteration method (36) in succession and obtained improved results.

Recently, research on the joint estimation of spectrum distribution and material maps has attracted much attention. In 2020, Chang et al. (37) proposed to simultaneously recover the spectrum and reconstruct the image. In their work, a muti-variable optimization was established and solved by the block-coordinate-descent (BCD) algorithm. Further, Zhao et al. (38) followed this framework and developed a spectral CT imaging method based on blind separation of polychromatic projections with Poisson prior, which improved the accuracy of narrow-energy-width projections. Similar frameworks have been developed by various modifications and applied to broader application scenarios (39). These methods eliminate the extra workload of spectral estimation and show emerging potential in alleviating the negative impact caused by errors from spectral mismatch. Despite the progress achieved in previous works, the capabilities of noise depression and convergence rate cannot fully meet the requirements in various medical imaging scenarios.

Deep learning is a promising technology, which has received increasing attention from imaging communities. For image domain-based methods, it is an effective method to compensate for the nonlinear effect. For example, in 2019, Zhang et al. (10) designed a butterfly network to realize end-to-end MD and performed well in clinical datasets. Subsequently, various improvements for network architecture, and loss function have been proposed (40,41). Recently, to overcome the large demand for data, the semi-supervised (42) method has also been used in this problem. For projection domain-based and one-step iterative reconstruction methods, it cannot only realize end-to-end mapping (43) but can also combine with the model-based iterative reconstruction (44). Despite the promising results, deep learning methods are data-driven with a strong dependency on data scale, meaning that they are difficult to apply when lacking datasets.

In this paper, we propose a material-image subspace decomposition-based iterative reconstruction (MISD-IR) method with spectrum estimation for DECT. Despite the promising results of previous works, there are more factors to be considered in clinical application. In clinical practice, convergence speed is a key criterion in considering whether a method can be applied. Besides, the noise is complex and needs to be taken seriously. In this work, model spectrum and ESART were introduced to improve computational efficiency. To avoid noise accumulation, we explored the similarity between reconstructed images and material images based on subspace decomposition. To further enhance the stability, the proximal block coordinate algorithm was introduced and weighted by adaptive weight. The proposed method was validated by numerical simulation and practical experiments. The experimental results showed that the proposed method can obtain accurate energy spectrum and material images while improving the convergence speed.

Methods

Forward transmission model for reconstruction

Assuming P represents the transmission data of DECT, the forward transmission model can be modeled as

Pk,L=∫sk(E)exp(−∑iμi(E)∫LFi(x)dl)dE[1]

where L∈Ωk,k=High,Low.L denotes the transmission path and Ωk denotes the collection of x-ray transmission path indexes corresponding to k-th spectrum. sk(E) represents the k-th normalized spectrum. It is the product of the emission spectrum of the X-ray tube, the material and thickness of the detector scintillator, and the sensitivity of the detector (23). e stands for X-ray energy. μi(⋅) and Fi(x) stands for the mass attenuation coefficient and the material maps in spatial x for material i. sk is the product of the emission spectrum, the material and thickness of the detector scintillator, and the sensitivity of the detector. Bones and water are the two basis materials we studied. Eq. [1] can be transferred as

Pk,L=∫sk(E)exp(−∫L(μb(E)Fb(x)+μw(E)Fw(x))dl)dE[2]

where b,w denote bone and water, respectively. To facilitate numerical implementation, the discretized representation of the Eq. [2] can be written as

Pk,L=∑E=1NkSk,Eδexp(−μb,EALFb−μw,EALFw)[3]

where Nk is the discretization number and δ is the step size of discretization for each spectrum. Sk,E is the short form of Sk(E) that means the value of k-th spectrum in energy E. μb,E,μw,E are also short forms for the mass attenuation in energy of bone and water, respectively. AL is the projection operator of path L, aL,x denotes a component of AL in position x. Fb,Fw are material maps of bone and water, respectively.

Spectrum estimation with model spectra

Inspired by Zhao’s work (23), the normalized spectrum can be represented by weighted known spectra. Assuming Sk,mdenotes known spectra, the k-th spectrum can be represented as

sk(E)=∑m=1Mkck,msk,m(E)[4]

where ck,mdenotes the weight coefficient and m denotes the tag of different model spectra. Thus, the transmission data Pk,L can be written as

Pk,L=∑E=1Nk∑m=1Mkck,msk,m,Eδexp(−μb,EALFb−μw,EALFw)[5]

where sk,m,E is a short form of sk,m(E). The model spectrum can be obtained by simulating the spectrum of different filter thicknesses under the same voltage and current conditions. Given material maps, when the estimated spectrum and material maps are consistent with the real spectrum and material maps, the estimated data is closest to the transmission data Pk,L. Thus, the solution of spectrum and material maps can be expressed as:

argminc,Fb,Fw∑k12‖FL(ck,Fb,Fw)−P˜k,L‖22,s.t.ck=(ck,1,⋯,ck,m),ck,m≥0,∑m=1Mkck,m=1[6]

where Pk,L=FL(ck,Fb,Fw), F(⋅) denotes the forward transmission function along the path L and ck denotes the set of weight coefficients of the different model spectrum. Considering the physical meaning, ck,m is non-negative and the sum of different coefficients is 1. More details can be found in the Appendix 1.

Extended simultaneous algebraic reconstruction technique

ESART (34) is an MD algorithm with fast convergence and certain noise suppression capability. For calculation convenience, the forward transmission model is expressed as:

Qk,L=−ln∑E=1Nksk,Eδexp(−μb,EALFb−μw,EALFw)[7]

According to ESART (34), material maps can be expressed as:

Fb,jn+1=Fb,jn+λb∑L∈ΩkaL,x∑L∈ΩkaL,x(eFb,Ln∑j=0JaL,j)[8]

Fw,jn+1=Fw,jn+λw∑L∈ΩkaL,x∑L∈ΩkaL,x(eFw,Ln∑j=0JaL,j)[9]

where x is the pixel in the material maps. a is the projection matrix. λb,λw are relaxation parameters. e represents the estimation error, and eb,Ln,ew,Ln can be expressed as:

eb,Ln=Φk,LnPk,Ln(Φk,Ln)2+(Θk,Ln)2(Q˜k,L−Qk,Ln),ew,Ln=Θk,LnPk,Ln(Φk,Ln)2+(Θk,Ln)2(Q˜k,L−Qk,Ln),Φk,Ln=∑E=1Nksk,Eδμb,Eexp(−μb,EALFbn−μw,EALFwn),Θk,Ln=∑E=1Nksk,Eδμw,Eexp(−μb,EALFbn−μw,EALFwn)[10]

where Q˜k,L is the negative logarithmic transmission data and Qk,Ln is the data calculated with n−th estimated material maps in Eq. [7].

MISD-IR with spectrum estimation

In this paper, the spectrum and material maps needed to be estimated simultaneously. For stability, the model spectral-based method was used. Considering the noise in the transmission data and the errors in estimated spectra, prior information was introduced into model , which can be written as:

argminc,Fb,Fw∑k12‖FL(ck,Fb,Fw)−P˜k,L‖22+αϕ(Fb,Fw),s.t.ck=(ck,1,⋯,ck,m),ck,m≥0,∑m=1Mkck,m=1[11]

where ϕ(⋅) denotes the regularization term and α is the parameter to balance the fidelity term and regularization term. In previous work, Chang et al. (37) introduced TV into the iteration and achieved certain effects. However, TV is suitable for piece-wise constant images, and medical images cannot simply be considered as piece-wise constant images, which can easily lead to blur or loss of image edges or texture details. It is necessary to design a new regularization term to suppress the noise.

The regularization term is designed based on the similarity between the reconstructed images corresponding to different spectra and the estimated material maps (5). In image domain-based decomposition methods, the reconstructed images can be considered as a linear combination of material maps, which can be written as:

(IHighILow)=(μbHighμwHighμbLowμwLow)(FbFw)[12]

where IHigh,ILow denote the reconstructed images corresponding to high and low energy spectra. μik denotes the mass attenuation coefficient of material i, namely, bone or water corresponding to k-th spectrum. Fb,Fw are bone and water material maps, respectively. Thus, it is a natural idea that there exists similarity between material maps and reconstructed images in structure and texture. According to the above idea, the regularization term was constructed based on the correlation between material maps and reconstructed images which are called prior images. In addition, the self-similarity of material maps was also introduced. In this paper, the exploration was conducted by subspace decomposition. First, a matrix was constructed based on material maps and reconstructed images, which can be expressed as:

y=[βbFb,βwFw,IHigh,ILow]T∈ℝ4×ns.t.βb≠0,βw≠0[13]

where βb,βw are constant scalars to distinguish different materials from each other, and T denotes matrix transpose. In y, for computational convenience, the combination of material maps and reconstructed images were transformed into a one-dimensional vector. n denotes the total number of pixels for material maps. According to Eq. [13], material maps can be expressed as:

Fb=y(1)βb,Fw=y(2)βw[14]

Further, according to the similarity, y has sparse and low-rank properties (45). Therefore, it is reasonable that y can be represented in subspace with a lower dimension. Assuming a 4×a dimensional subspace sa⊂ℝ4×awith a≤4, thus y can be expressed as:

y=HaZ,s.t.HaTHa=In[15]

where Ha:=[h1,⋯,sa]∈ℝ4×a contains the basis of subspace sa with orthogonality. In∈ℝn×n denotes identity matrix. Z∈ℝa×n contains the eigen-image of y. As a low-dimensional representation of y, Z provides an opportunity to exploit the non-local similarity of y. Thus, the regularization term can be expressed as:

argminZ,Ha12‖HaZ−y‖F2+αZαφ(Z),s.t.HaTHa=In[16]

where 12‖HaZ−y‖F2 is the fidelity term to ensure that the subspace is learned from the transmission data. φ(⋅) is the module to exploit self-similarity of eigen-image such as BM3D (46), EPLL (47), or self-supervised methods (48,49). αZα is the parameter to balance the two terms. In this work, the BM3D framework was utilized as the regularization term in the above problem, thus, a stable solution for Eq. [16] could be directly used in the subspace decomposition algorithm proposed by Zhuang (45). Combined with Eq. [11] and Eq. [14], the objective function can be written as:

argminc,Fb,Fw,y,Ha,Z∑k12‖FL(ck,Fb,Fw)−P˜k,L‖22+λc∑k(‖ck‖1−1)2+α2‖Z−HaTy‖F2+αZφ(Z)+γb2‖Fb−y(1)βb‖22+γw2‖Fb−y(2)βw‖22+γH2‖IH−y(3)‖22+γL2‖IL−y(4)‖22s.t.ck,m≥0,βb≠0,βw≠0,∀m,k[17]

where ∑k(‖ck‖1−1)2 is normalized constraint. ‖⋅‖1 is the L1-norm and λc is set to be ∞. Considering y is constructed by material maps and reconstructed images, γb2‖Fb−y(1)βb‖22, γw2‖Fb−y(2)βw‖22, γHigh2‖IHigh−y(3)‖22 and γLow2‖ILow−y(4)‖22 are added. γb,γw,γHigh,γLow are parameters to balance different regularization and fidelity terms. To fully utilize the prior information in reconstructed images, γHigh,γLow are set to be ∞ to make sure the effective information in reconstructed images not lost.

Optimization algorithm

Considering that the solution of the Eq. [17] is nonconvex for all the variables in the equation, the proximal regularized block coordinate descent (P-BCD) (50) was introduced and modified. Due to the influence of noise and errors in spectra and material maps, the iteration became unstable as the convergence speed became faster. Compared with BCD, P-BCD was more stable for iteration. Despite multiple variables that needed to be solved in Eq. [17], it could be divided into four blocks, namely, spectral weight c, material maps Fb,Fw, regularization term Z,Ha and auxiliary variable y.

BFGS-based estimation of spectra

Assuming when estimating the spectrum, material maps are fixed as Fbn,Fwn,yn,Zn,Han, removing the terms not related to c, Eq. [17] can be written as:

argminc∑k12‖FL(ck,Fbn,Fwn)−P˜k,L‖22+λc∑k(‖ck‖1−1)2s.t.ck,m≥0,∀m,k[18]

At nth iteration, the objective function is a convex subproblem. On account of the noise and errors in material maps, a variant of BCD, namely, proximal regularized BCD (P-BCD), was introduced and modified to improve the stability of the algorithm. When a proximal operator is introduced, Eq. [18] can be expressed as:

cn+1=argminck∑k12‖FL(ck,Fbn,Fwn)−P˜k,L‖22+λc∑k(‖ck‖1−1)2+12ηc,kn‖ck−ckn‖22s.t.ck,m≥0,∀m,k[19]

where {ηcn} are positive scalar quantities. When the spectra is more accurate, it will be more sensitive to errors in material maps and noise in transmission data. Thus, we modified P-BCD based on the value of {ηcn} to control the step size. In this paper, the value of {ηcn} was determined by the difference between the estimated data and transmission data, which can be expressed as:

ηc,kn=χk‖FL(ckn,Fbn,Fwn)−P˜k,L‖22[20]

where χk is a positive scalar quantity. Referring to the solution of Eq. [19], the main consideration in the design for the solution of c was the rate of convergence. The quasi-Newton method has the characteristics of fast convergence and low computational complexity. BFGS-B (51) is one of the most popular quasi-Newtonian methods. Its validity in spectrum estimation problem has been verified in (37).

The solution of material maps based on ESART and proximal iteration

When other variables are fixed, the solution of material maps can be expressed as:

argminFb,Fw∑k12‖FL(ckn,Fb,Fw)−P˜k,L‖22+γb2‖Fb−yn(1)βb‖22+γw2‖Fw−yn(2)βw‖22s.t.λb≠0,λw≠0[21]

Note that the theoretical exact analytical solution for the above problem is difficult to obtain since the forward transmission model is not linear. Moreover, it is more valuable and useful to obtain an approximate but computationally efficient solution in an iterative algorithm. To get an approximate solution, a two-step method was proposed by dividing the Eq. [21] into two subproblems. In this method, ESART and proximal iteration are utilized to efficiently solve each problem. First, the material maps are updated by the fidelity term, which can be expressed as:

argminFb,Fw∑k12‖FL(ckn,Fb,Fw)−P˜k,L‖22[22]

In this paper, the scanning sets of different energy spectra are consistent. Under such conditions, ESART has a high convergence efficiency and accuracy (34), which has been widely used in medical diagnosis and industrial testing (52). Additionally, ESART has good noise suppression characteristics, which is conducive to dealing with noise during MD (34). Thus, ESART was chosen to solve Eq. [22], which can be expressed as:

Fb,jn+12=Fb,jn+λb∑L∈TaL,x∑L∈ΩkaL,x(eFb,Ln∑j=0JaL,j)[23]

Fw,jn+12=Fw,jn+λw∑L∈TaL,x∑L∈ΩkaL,x(eFw,Ln∑j=0JaL,j)[24]

where x is the pixel in the material maps, a is the projection matrix, λb,λw are relaxation parameters, e represents the estimation error and eFb,Ln,eFw,Ln can be calculated by Eq. [10], and n+12 denotes the transition state between n and n+1 iteration.

Second, the material maps are updated according to the proximal terms of:

argminFb,Fw12‖Fb−Fbn+12‖22+γb2‖Fb−yn(1)βb‖22+12‖Fw−Fwn+12‖22+γw2‖Fw−yn(2)βw‖22[25]

and the final solution can be expressed as:

(Fbn+1Fwn+1)=(ρbFbn+12ρwFwn+12)+(1−ρbβbyn(1)1−ρwβwyn(2))[26]

where ρb=11+γb，ρw=11+γw. As γb,γw increase, (n+1)−th material maps would be closer to n−th material maps.

The solution of regularization term based on subspace decomposition

In the solution of the regularization term, Eq. [17] can be rewritten as:

argminHa,Zα2‖Z−HaTyn‖F2+αZφ(Z)[27]

For the solution of Eq. [27], the main consideration was to fully utilize the information in reconstructed images and improve the computational efficiency. Considering the similarity between reconstructed images and material images, subspace decomposition is a useful tool. In this paper, it was solved by the subspace decomposition method proposed by Zhuang and Bioucas-Dias (45). It is a fast and stable solution for eigen-image Z. However, we believe that with the development of deep learning, it has the potential to be improved by self-supervised methods. The application of the self-supervised method is described in the discussion section.

The solution of y based on the stagnation point

Having updated Fb,Fw,c,Ha, and Z, the solution of y can be expressed as:

argminyγb2‖Fb−y(1)βb‖22+γw2‖Fw−y(2)βw‖22+α2‖Z−HaTy‖F2+γHigh2‖IHigh−y(3)‖22+γLow2‖ILow−y(4)‖22[28]

y(3),y(4) are forced to be close to reconstructed images according to the definition of y. Thus, Eq. [28] can be expressed as:

argminy(1),y(2)γb2‖Fbn+1−y(1)βb‖22+γw2‖Fwn+1−y(2)βw‖22+α2‖Zn+1−(HaT)n+1y‖F2[29]

Assuming that function G is defined as:

G(y(1),y(2))=γb2‖Fbn+1−y(1)βb‖22+γw2‖Fwn+1−y(2)βw‖22+α2‖Zn+1−(HaT)n+1y‖F2[30]

Considering Eq. [30] is convex, the minima of this function can be taken when its derivative is an all-zero vector, which can be expressed as:

∂G∂y(1)=γbβb(y(1)βb−Fbn+1)+∂α2‖Zn+1−(HaT)n+1y‖F2∂y(1),∂G∂y(2)=γwβw(y(2)βw−Fwn+1)+∂α2‖Zn+1−(HaT)n+1y‖F2∂y(2),[31]

where

∂α2‖Zn+1−(HaT)n+1y‖F2∂y(1)=α(y(1)−y^(1)),∂α2‖Zn+1−(HaT)n+1y‖F2∂y(2)=α(y(2)−y^(2)),y^=Han+1Zn+1[32]

Y^ is obtained by estimated eigen-image. Thus, at n−th iteration, Y can be expressed as:

yn+1(1)=βb2αβb2+γb(γbβb2Fbn+1+αy^(1))，yn+1(2)=βw2αβw2+γw(γwβw2Fwn+1+αy^(2)),yn+1(3)=IHigh,yn+1(4)=ILow[33]

Considering that Eq. [29] is the sum of quadratic terms, it is easy to prove that is a convex function. Thus, Eq. [33] is the solution of Eq. [28].

MISD-IR aims to reconstruct material maps and recover spectra simultaneously with fast convergence speed and obtain accurate results. Three techniques were introduced to improve the convergence speed and noise suppression capability. Firstly, to improve the convergence speed, ESART and model spectral were incorporated into this method. Secondly, subspace decomposition was adopted to explore the similarity between material maps and reconstructed CT images to suppress the noise. Finally, the idea of P-BCD was introduced to stabilize the solution and suppress the noise propagation. The numerical data obtained by Siemens SOMATOM Definition Flash CT scanner (Siemens, Erlangen, Germany) of our university and the real scanning data obtained by Thales Hawkeye-130 (Thales, La Defense, France) and Varian 4030E (Varian, Palo Alto, CA, USA) were used to verify the effectiveness of the proposed MISD-IR. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study has been granted exemption from ethical approval and patients’ consent by the Ethics Committee of PLA Strategic Support Force Information Engineering University due to not involving human trials, clinical diagnosis, treatment information, sensitive personal information, commercial interests, or causing harm to the human body as well as the retrospective nature of this study.

Results

In this paper, numerical simulations, real scanning, and ablation experiments were conducted to evaluate the proposed method. Direct inverse, dual-energy methods based on iterative reconstruction (DE-IR) with FBP images (5) and spectrum estimation-guided iterative reconstruction algorithm (SG-IR) (37) were considered for comparison. All the calculations in the experiments were conducted in a workstation containing an Intel^® Xeon^® Gold 6226R Processor (Intel, Santa Clara, CA, USA).

Numerical simulations setting

In this section, preclinical data of head and abdomen were used to evaluate the performance of the proposed method. These data were obtained with the help of radiologists. They corrected the material maps with Syngo.Via equipped in Siemens SOMATOM Definition Flash CT scanner. Under the constraint of confidentiality, patient information was treated confidentially. For every dataset, a representative slice was chosen to apply in the experiments, and shown in Figure 1.

Figure 1 The slices used in simulation experiments, the normalized spectral used in numerical simulation, and the linear attenuation coefficient of bone and water. (A) The slice of head and it is displayed in [−300, 1,100] HU; (B) the slice of abdomen and it is displayed in [−400, 400] HU; (C) the normalized spectrum; (D) the linear attenuation coefficient. HU, Hounsfield unit.

For simplicity, a fan beam geometry was applied and the geometry was kept consistent in dual-energy scanning. The polychromatic X-ray spectra were simulated with ImaSim (53). As shown in Figure 1, the X-ray spectra were simulated with a tube voltage of 80 and 140 kV, which was filtered with 5 mm aluminum and 1,000 mm air. The model spectra were simulated with different thicknesses of aluminum, which can be seen in Figure S1. The sampling interval was 1 keV for all the spectra. Bone and water were used as basis materials. The mass attenuation coefficients were from the National Institute of Standard Technology (NIST) tables of X-ray mass attenuation coefficients (54). The scanning configuration was set as follows. The distance from the X-ray source to the rotation center (SOD) and the detection center (SDD) were 300 and 600 mm, respectively. The detector was composed of 512 probes and the size of the probe was 0.124 mm. For each spectrum, the scan angle was 360 degrees with 720 projections. The size of reconstructed images and material maps was 512×512. To verify the noise suppression capability of the proposed methods, three different levels of Poisson noise were added and the photon number was set to be 1×10⁶, 1×10⁵ and 5×10⁴.

To verify the generalization performance, two more different sets of kVps were set under the photon number of 1×10⁶. The tube voltage of the two sets was 70/120 kVp and 60/100 kVp. The rest of the experimental settings remained the same as before.

To further quantitatively evaluate the proposed method and estimate whether the CT number [Hounsfield unit (HU)] of the CT image is preserved by the sum of the images obtained by MD, virtual monoenergetic images (VMIs) were introduced. First, the ground truth and the generated material images corresponding to different methods were used to construct VMIs of 40 kVp. Second, the value of these VMIs was converted to CT number. Then, the mean value and standard deviations were calculated based on VMIs.

Additionally, to verify the application range of the proposed method, we further tested the proposed method in experiments with different sets of materials. In this experiment, considering the practical application in medical detection, contrast media-5% iodine and water materials were used with a newly designed phantom. The spectra and scanning parameters used in this experiment were the same as those used in the previous numerical experiments. Figure 2 shows the VMI in 40 keV of the phantom used in this experiment.

Figure 2 The virtual monoenergetic image in 40 keV of the phantom used in this experiment. The display window is [0, 1,000] HU. HU, Hounsfield unit.

In numerical simulation, structural similarity (SSIM), and peak signal-to-noise ratio (PSNR) are used to evaluate the performance of different methods in material images. PSNR can evaluate the distortion degree of the generated material images compared with the ground truth, and SSIM can evaluate the structural similarity between the generated material image and the ground truth. Root mean square error (RMSE) is used to evaluate the accuracy of the estimated spectra. To verify the convergence speed, the running time and iterative rounds for SG-IR and MISD-IR were recorded and shown.

Real scanning setting

Practical experiments were carried out based on real scanning data with laboratory dual-energy CT. An anthropomorphic head phantom (55) was used to be scanned under 80 and 140 kVp. The X-ray source and detector were Thales Hawkeye-130 and Varian 4030E, respectively. The distance between the source and the object was 623.09 mm and the distance between the source and detector was 865.05 mm. The detector contained 512 bins and the size of the bin was 0.83 mm. The sampling interval was 1 degree, and 360 angles of data were collected. The size of reconstruction images was 512×512.

Ablation experiments setting

To evaluate the performance of different parts in the proposed method, we modified the model based on these parts and carried out experiments. In this part, the head dataset with 1×10⁶ Poisson noise was used to validate various parts of the proposed method. The results were evaluated by PSNR. The proposed algorithm was divided into three parts, namely, spectral estimation, regularization term, and ESART. The regularization term and spectrum estimation were eliminated in turn to verify the performance of these two terms. When spectrum estimation was eliminated, the spectra were given as the initial value.

Numerical simulations results

Figure 3 shows the reconstruction results by different methods. The results indicate that direct inversion is much more seriously degraded by the noise compared to the latter three methods. Besides, the edge of DE-IR was significantly darker than the center of the image, which resulted from beam-hardening artifacts. Compared with the image-domain-based methods, SG-IR and the proposed methods had a remarkable effect on noise suppression and beam-hardening artifacts. However, the proposed method showed advantages in detail preservation over direct inversion, DE-IR, and SG-IR. To further clarify this point, two regions of interest (ROIs) were selected and marked with yellow squares. They were enlarged in the lower left corner of every single image. It is obvious that the edges became blurred and tiny structures were missed in DE-IR and SG-IR. Using the proposed method, the edges were sharper and the structures were preserved to a larger extent.

Figure 3 Reconstruction results of preclinical data by different methods. (A1-A5) Bone images of the head, their displayed grayscale is [0, 0.9]; (B1-B5) water images of the head, their displayed grayscale is [0, 1.0]; (C1-C5) bone images of abdomen, their displayed grayscale is [0, 1.0]; (D1-D5) water images of abdomen, their displayed grayscale is [0, 0.9]. (A-D) The results of ground truth, Direct Inverse, DE-IR, SG-IR, and the proposed method, respectively. Regions of interest are marked by yellow boxes and enlarged in the lower left corner. DI, direct inverse; DE-IR, dual-energy methods based on iterative reconstruction; SG-IR, spectrum estimation-guided iterative reconstruction algorithm; MISD-IR, material-image subspace decomposition-based iterative reconstruction with spectrum estimation.

Figure 4 shows the difference images between the ground truth and the results for different methods. To better show the performance of different methods, the difference images were displayed in the JET color space. Compared with image-domain-based methods, SG-IR and the proposed method were closer to ground truth as a whole and their difference maps were more uniform. Compared with SG-IR, the difference between the proposed method and the ground truth on edges was obviously smaller.

Figure 4 Difference images between the results of different methods and ground truth. (A1-D1) Difference images of the head’s bone material, their displayed grayscale is [0, 0.2]. (A2-D2) Difference images of the head’s water material, their displayed grayscale is [0, 0.4]. (A3-D3) Difference images of abdomen’s bone material, their displayed grayscale is [0, 0.2]. (A4-D4) Difference images of abdomen’s water material, their displayed grayscale is [0, 0.2]. DI, direct inverse; DE-IR, dual-energy methods based on iterative reconstruction; SG-IR, spectrum estimation-guided iterative reconstruction algorithm; MISD-IR, material-image subspace decomposition-based iterative reconstruction with spectrum estimation.

To quantitatively evaluate the decomposition accuracy for different methods, SSIM and PSNR were calculated and are listed in Table 1. There were two different levels of noise added in the experiments. In different data sets with different noise levels, the quantitative results of the proposed method indicated obvious advantages over the others. Compared with the image-domain-based methods, SG-IR and the proposed method had higher precision, confirming the findings verified by Figures 3,4. Compared with SG-IR, the proposed method also showed better performance. In the head dataset, the SSIM and PSNR of the proposed method were 0.01 and 3.02 dB higher, respectively, than those of SG-IR in bone material with photon number of 1×10⁶. For water material, the proposed method was 6.23 dB higher than SG-IR in PSNR. In the abdomen dataset, the proposed method was 0.03 and 2.95 dB higher, respectively, than SG-IR in bone material. In water material, the proposed method was 0.03 and 5.54 dB higher, respectively, than SG-IR. When the photon number became less, specifically, the noise level became higher, the proposed method also performed better than other methods, which demonstrated that the proposed method is robust to noise. Additionally, to further test the performance of the proposed method in the case of severe noise, the experiment under photon number of 5×10⁴ was carried out with the abdomen dataset. According to the results, the PSNR of the bone material map generated by the proposed method was 1.92 dB better than others at least. The SSIM of the proposed method was 0.26 higher than others at least. Regarding water material maps, the PSNR of the proposed method was 2.99 dB higher than others at least. The SSIM of the proposed method was 0.22 higher than others at least. Even with severe noise, the proposed method still had a better performance than other methods.

To demonstrate the convergence of the proposed method, the behavior curves of RMSEs based on the estimated water material of the experiment under photon number of 1×10⁶ are plotted in Figure 5.

Figure 5 The convergence curve of SG-IR and the proposed method. The blue line denotes the root mean square error of SG-IR and the red line denotes the root mean square error of MISD-IR in the abdomen dataset. RMSE, root mean square error; MISD-IR, material-image subspace decomposition-based iterative reconstruction with spectrum estimation; SG-IR, spectrum estimation-guided iterative reconstruction algorithm.

Additionally, the converged iteration number of the proposed method was around 50 and that of SG-IR was around 200. The average time for the proposed method was 1,425 seconds whereas that for SG-IR was 5,700 seconds, demonstrating the fast convergence of the proposed method.

Table 2 shows the RMSEs of the estimated spectra with respect to the reference spectra. Compared with SG-IR, the estimated spectra of the proposed method was closer to the reference spectra. When the photon number was 1×10⁶, the estimated spectra of the proposed method was 0.005 and 0.01 lower than that of SG-IR in the head dataset. When the photon number was 1×10⁵, the estimated high and low spectra of the proposed method were 0.005 and 0.007 lower than those of SG-IR. In the abdomen dataset, when the photon number was 1×10⁶, the estimated high and low spectra were 0.023 and 0.018 lower than those of SG-IR. When the photon number was 1×10⁵, the estimated high and low spectra were 0.021 and 0.022 lower than those of SG-IR. When the photon number was 5×10⁴, the estimated high and low spectra were 0.022 and 0.028, lower than those of SG-IR. The numerical results showed that the proposed method is robust to different noise level.

Figure 6 shows the VMIs generated by the material maps of different methods. The VMIs of two image domain-based methods were obviously different from the reference VMI. In the VMI of DI, there were serious noise and artifacts. In the VMI of DE-IR, the noise was reduced to some extent while the artifacts were still serious. The VMIs of SG-IR and MISD-IR were close to the ground truth. There were little artifacts in the two VMIs. However, there was still some noise in SG-IR, which influenced the image quality. In the VMI of MISD-IR, there were no visual noise and artifacts. Generally, the VMI of the proposed MISD-IR is close to the ground truth with high image quality. To quantitatively evaluate the proposed method and estimate whether the CT number was preserved by the sum of the images obtained by MD, the ROIs in the same position of different VMIs were chosen to calculate mean values and SDs of different methods. The ROI of the ground truth was marked with yellow boxes. According to the mean values of different methods, the image domain-based methods were further away from the ground truth compared with SG-IR and MISD-IR. The difference between the mean value of the proposed method and the mean value of ground truth was about 26.1927. According to the SDs of different methods, the proposed MISD-IR had a better noise suppression ability. Despite that the SD of MISD-IR was lower that of ground truth, the SD of MISD-IR was closer to that of ground truth compared with other methods.

Figure 6 Virtual monoenergetic images of 40 keV generated by the material maps of different methods. (A1) The reference VMI generated by reference material maps; (A2) the VMI generated by the material maps of DI; (A3) the VMI generated by the material maps of DE-IR; (A4) the VMI generated by the material maps of SG-IR; (A5) the VMI generated by the material maps of MISD-IR. A ROI is marked by a yellow box. The ROI in each VMI is used to calculate their mean values and SDs that are shown below the VMI of various methods. The unit of the mean value is HU. The display window is [−400, 400] HU. DI, direct inversion; DE-IR, dual-energy methods based on iterative reconstruction; SG-IR, spectrum estimation-guided iterative reconstruction algorithm; MISD-IR, material-image subspace decomposition-based iterative reconstruction with spectrum estimation; ROI, region of interest; SD, standard deviation; VMI, virtual monoenergetic image; HU, HU, Hounsfield unit.

Table 3 shows the numerical results of different sets of different kVps. According to the results, the SSIMs of the three sets were similar. The difference of PSNR was limited in 1.14 dB. The proposed method had similar effectiveness for different kVps in MD. In the results of spectra, the difference for the RMSEs of high kVp spectra was limited to 0.007 and the difference for the RMSEs of low kVp spectra is limited to 0.006. The numerical results showed that the proposed method has generalization performance to some extent. To our knowledge, the difference in the experiments under different kVps is acceptable because the information corresponding to different kVps is different such as image contrast and attenuation information (56).

Figure 7 shows the results of different methods. According to Figure 2, in the results of the two image-domain-based methods, there are obvious artifacts. Despite the noise suppression ability, the image quality of DE-IR was unsatisfactory. In the results of SG-IR and MISD-IR, artifacts were almost invisible and the image quality had been improved. According to the difference maps, the results of MISD-IR were closer to the ground truth compared with SG-IR. Generally, the results generated by MISD-IR were better than those of other methods, which preliminarily demonstrates the effectiveness of MISD-IR for different material. To further test the effectiveness of the proposed method, we evaluated these results quantitatively. Table 4 shows the quantitative results of different materials. The SSIM and PSNR of MISD-IR were obviously higher than those of other methods which confirms the effectiveness of the proposed method in different sets of materials. Additionally, the RMSEs of spectral estimation of SG-IR for the spectra of 140 and 80 kVp were 0.0369 and 0.0470, respectively. Those for MISD-IR were 0.0311 and 0.0511, respectively. Generally, MISD-IR is still effective for spectral estimation and can be applied in different sets of materials effectively.

Figure 7 The results generated by different methods. (A1-A5) The results of iodine material maps, the display grayscale is [0.2, 1.2]; (B1-B5) the difference maps between the results of iodine material maps generated by different methods and the ground truth, the display grayscale is [0, 0.2]; (C1-C5) the results of water material maps, the display grayscale is [0.2, 1.2]; (D1-D5) the difference maps between the results of water material maps generated by different methods and the ground truth, the display grayscale is [0, 0.2]. DI, direct inverse; DE-IR, dual-energy methods based on iterative reconstruction; SG-IR, spectrum estimation-guided iterative reconstruction algorithm; MISD-IR, material-image subspace decomposition-based iterative reconstruction with spectrum estimation.

Real scanning results

Figure 8 shows the reconstruction results of different materials. Three iterative methods had a significant effect on noise suppression. In the results of the selected slices, DE-IR, and the proposed method showed better results on noise suppression. Different from the ideal dataset, the real scanning data was greatly affected by noise and beam-hardening artifacts. Therefore, the results of SG-IR using total variation were not so effective in slice 2.

Figure 8 Reconstruction results of real scanning data. The first row shows the reconstruction results of bone material for slice 1 and the display scale is [0.3, 0.9]. The second row shows the reconstruction results of water material for slice 1 and the display scale is [0.6, 1.2]. The third row shows the reconstruction results of bone material for slice 2 and the display scale is [0.2, 1.0]. The fourth row shows the reconstruction results of water material for slice 2 and the display scale is [0.4, 1.4]. Columns 1–4 show the results of different methods, which are direct inverse, DE-IR, SG-IR and, MISD-IR from left to right. The region of interest is signed by squares and enlarged in the lower left corner. Two ROIs are marked with yellow boxes to calculated mean values and SDs. The red arrows in the figure are set to point out typical detail and structure to show the characteristics of different methods. DI, direct inverse; DE-IR, dual-energy methods based on iterative reconstruction; SG-IR, spectrum estimation-guided iterative reconstruction algorithm; MISD-IR, material-image subspace decomposition-based iterative reconstruction with spectrum estimation; ROI, region of interest; SD, standard deviation.

As shown by the enlarged ROI and the edge pointed by the arrow in slice 1, there is obviously visible improvement for edge preservation in the results of the proposed method. The edges in the results of DE-IR and SG-IR are seriously blurred whereas those of the proposed method are sharper and clearer. In the results in slice 2, an area of teeth is enlarged and shown in the third row. In the proposed method, the boundaries between different teeth are clearer and the edges of each tooth are also sharper. In the fourth row, a yellow arrow points at the region of soft tissue around the teeth. In DE-IR and SG-IR, some details are missing in these areas yet they are well preserved in the results of direct inversion and the proposed method, which demonstrates the promising detail retention capability of the proposed methods. Thus, the merits of the proposed method were also further verified in real scanning.

We further evaluated the proposed method by the mean and standard deviation of ROIs, as shown in Figure 8 (marked with yellow boxes). Table 5 shows the quantitative results of decomposition images. Generally, the proposed method had lower SDs, which stands for better noise suppression ability. In the results of direct inversion, the SDs of ROIA and ROIB reached 0.1423 and 0.1545, respectively. Compared with direct inversion, the proposed MISD-IR reduced SDs by 86.09% and 78.25%, respectively. Compared with image-domain iterative methods DE-IR, MISD-IR reduced SDs by 22.05% and 52.20%, respectively. Compared with SG-IR, MISD-IR reduced SDs by 47.34% and 49.85%, respectively. The quantitative results indicated that the proposed MISD-IR has the potential to achieve high quality decomposition.

Ablation experiments results

Table 6 shows the PSNRs of the proposed method with various modifications. As shown in the table, when the material images were estimated by initial spectrum without regularization term, the accuracy was the worst. When the spectral estimation term was added, the material images could be much better estimated. When the spectra could not be updated, the accuracy of the material images could not be well guaranteed despite the regularization term. When the spectra were updated, the regularization term could significantly improve the accuracy of the material images.

Discussion

In this paper, we have investigated the problem of simultaneous estimation of the spectral distribution and material maps. We have proposed a MISD-IR method. To improve the convergence speed, ESART and model spectral are introduced and combined in a united optimized model. Besides, we explored the similarity between reconstructed CT images and material maps to suppress the noise. Finally, to stabilize the iteration and better integrate material reconstruction and spectrum estimation, the idea of P-BCD is introduced and modified in this paper. Numerical and real scanning results demonstrated the improvement in convergence speed and the accuracy of the proposed method. According to the results under different levels of noise, the application potential of MISD-IR in low-dose DECT imaging could be seen. In future work, we plan to extend MISD-IR to different application scenarios. We have noted the role of DECT in security screening (57,58), and we believe that MISD-IR has the potential to be applied in security screening as well. However, compared with noise, metal artifacts may cause information loss and involve information completion problems, which we plan to study in future work.

In the ablation experiment, spectrum estimation and regularization terms were eliminated in turn, which verified the performance of different parts of the method. However, it is worth noting that when initial spectra are fixed in iteration, the regularization term is of limited use. The main reason is when initial spectra are fixed, there is too much error in the estimated material maps to explore the similarity between material maps and reconstructed images. Actually, the accuracy of spectral estimation is greatly affected by the initial value. Although the spectrum can be updated in MISD-IR, the convergence speed will be affected by the initial value. Image domain methods can obtain material attenuation information without spectra. Thus, material maps obtained by image-domain methods can be used to estimate the initial value of MISD-IR in practice.

MISD-IR is designed under the assumption of a lack of data. Despite the promising results, the determination of iteration parameters requires careful consideration. Besides, the application of deep learning methods in noise suppression is promising. Self-supervised methods can be applied with the noise suppression of eigen-images. In the future, these ideas will be the direction of our work. Under the assumption of access to large-scale data becoming possible, we have a preliminary design about the above ideas. In our opinion, the main challenge is the construction of dataset for training. Considering that the noise in eigen-images is different from that in the material images and the reconstructed images, it is not feasible to simply apply the pre-trained model. In our design, first, according to the idea in reference (59-61), we plan to construct the proposed iterative algorithm into an unrolling network, which can optimize the process of parameter selection. Second, to optimize the effect of noise suppression, the deep learning network was designed to be embedded in the unrolling network. Considering that compared with unpaired data, the acquisition of paired data sets is more difficult, the application of a self-supervised method is easier and cheaper. To our knowledge, the methods mentioned in blind spot scheme such as in (62) and deep generative model as in (63-66) are feasible for the suppression of intrinsic image noise, of which the noise suppression ability has been verified. Third, to avoid information loss in the process of noise suppression, in the loss function, the difference term between the different types of images is required. Compared with the method proposed in this paper, the above design requires a large amount of data, but the parameter selection and noise suppression can be further optimized.

Conclusions

In conclusion, we have proposed an iterative reconstruction based on material-image subspace decomposition with spectrum estimation for DECT. Numerical and real scanning experimental results show that the proposed method has better image edge and detail retention ability compared with other methods in this paper. Compared with the other method of the same type, the proposed method has a faster convergence rate, which added its clinical practicability. Our future work will focus on combining the proposed approach with unsupervised learning to further optimize its performance and promote its practical application.

Acknowledgments

Funding: This work was supported by the National Natural Science Foundation of China (Grant Nos. 62101596, 62201616, and 62271504). This work was also supported by Technology Innovation Leading Talent Project of Zhongyuan (Grant No. 244200510015) and China Postdoctoral Science Foundation (Grant No. 2023T160792).

Footnote

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://qims.amegroups.com/article/view/10.21037/qims-23-1681/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work, including ensuring that any questions related to the accuracy or integrity of any part of the work have been appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The study has been granted exemption from ethical approval and patient consent by the Ethics Committee of PLA Strategic Support Force Information Engineering University due to not involving human trials, clinical diagnosis, treatment information, sensitive personal information, commercial interests, or causing harm to the human body as well as the retrospective nature of this study.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Zheng Y, Han X, Jia X, Ding C, Zhang K, Li H, Cao X, Zhang X, Zhang X, Shi H. Dual-energy CT-based radiomics for predicting invasiveness of lung adenocarcinoma appearing as ground-glass nodules. Front Oncol 2023;13:1208758. [Crossref] [PubMed]
2. Manerikar A, Li F, Kak AC. DEBISim: A simulation pipeline for dual energy CT-based baggage inspection systems1. J Xray Sci Technol 2021;29:259-85. [Crossref] [PubMed]
3. Bauer C, Wagner R, Orberger B, Firsching M, Ennen A, Garcia Pina C, Wagner C, Honarmand M, Nabatian G, Monsef I. Potential of Dual and Multi Energy XRT and CT Analyses on Iron Formations. Sensors (Basel) 2021;21:2455. [Crossref] [PubMed]
4. McCollough CH, Boedeker K, Cody D, Duan X, Flohr T, Halliburton SS, Hsieh J, Layman RR, Pelc NJ. Principles and applications of multienergy CT: Report of AAPM Task Group 291. Med Phys 2020;47:e881-912. [Crossref] [PubMed]
5. Niu T, Dong X, Petrongolo M, Zhu L. Iterative image-domain decomposition for dual-energy CT. Med Phys 2014;41:041901. [Crossref] [PubMed]
6. Altunbas C. Feasibility of dual-energy CBCT material decomposition in the human torso with 2D anti-scatter grids and grid-based scatter sampling. Med Phys 2024;51:334-47. [Crossref] [PubMed]
7. Winklhofer S, Lambert JW, Sun Y, Wang ZJ, Sun DS, Yeh BM. Pelvic Beam-Hardening Artifacts in Dual-Energy CT Image Reconstructions: Occurrence and Impact on Image Quality. AJR Am J Roentgenol 2017;208:114-23. [Crossref] [PubMed]
8. Wang Y, Cai A, Liang N, Yu X, Zhong X, Li L, Yan B. One half-scan dual-energy CT imaging using the Dual-domain Dual-way Estimated Network (DoDa-Net) model. Quant Imaging Med Surg 2022;12:653-74. [Crossref] [PubMed]
9. Alvarez RE, Macovski A. Energy-selective reconstructions in X-ray computerized tomography. Phys Med Biol 1976;21:733-44. [Crossref] [PubMed]
10. Zhang W, Zhang H, Wang L, Wang X, Hu X, Cai A, Li L, Niu T, Yan B. Image domain dual material decomposition for dual-energy CT using butterfly network. Med Phys 2019;46:2037-51. [Crossref] [PubMed]
11. Zhang T, Yu H, Xi Y, Wang S, Liu F, Spectral CT. Image-Domain Material Decomposition via Sparsity Residual Prior and Dictionary Learning. IEEE Transactions on Instrumentation and Measurement 2023;72:1-13.
12. Wang S, Wu W, Cai A, Xu Y, Vardhanabhuti V, Liu F, Yu H. Image-spectral decomposition extended-learning assisted by sparsity for multi-energy computed tomography reconstruction. Quant Imaging Med Surg 2023;13:610-30. [Crossref] [PubMed]
13. Wang S, Yu H, Xi Y, Gong C, Wu W, Liu F. Spectral-Image Decomposition With Energy-Fusion Sensing for Spectral CT Reconstruction. IEEE Transactions on Instrumentation and Measurement 2021;70:1-11.
14. Wang S, Cai A, Wu W, Tao Z, Liu F, Yu H. IMD-MTFC: Image-Domain Material Decomposition via Material-Image Tensor Factorization and Clustering for Spectral CT. IEEE Transactions on Radiation and Plasma Medical Sciences 2023;7:382-93.
15. Kalender WA. X-ray computed tomography. Phys Med Biol 2006;51:R29-43. [Crossref] [PubMed]
16. Ducros N, Abascal JFP, Sixou B, Rit S, Peyrin F. Regularization of nonlinear decomposition of spectral x-ray projection images. Med Phys 2017;44:e174-87. [Crossref] [PubMed]
17. Goh KL, Liew SC, Hasegawa BH, editors. Energy dependent systematic errors in dual-energy X-ray CT. 995 IEEE Nuclear Science Symposium and Medical Imaging Conference Record, San Francisco, CA, USA, 1995;3:1587-91.
18. Goh KL, Liew SC, Hasegawa BH, editors. Correction of energy-dependent systematic errors in dual-energy X-ray CT using a basis material coefficients transformation method. EEE Transactions on Nuclear Science, 1997;44:2419-24.
19. Goh KL, Liew SC. Dual-energy x-ray approach for object/energy-specific attenuation coefficient correction in single-photon emission computed tomography: effects of contrast agent. J Med Imaging (Bellingham) 2021;8:052106. [Crossref] [PubMed]
20. Ruth C, Joseph PM. Estimation of a photon energy spectrum for a computed tomography scanner. Med Phys 1997;24:695-702. [Crossref] [PubMed]
21. Sidky EY, Yu L, Pan X, Zou Y, Vannier M. A robust method of x-ray source spectrum estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data. J Appl Phys 2005;97:124701.
22. Ha W, Sidky EY, Barber RF, Schmidt TG, Pan X. Estimating the spectrum in computed tomography via Kullback-Leibler divergence constrained optimization. Med Phys 2019;46:81-92. [Crossref] [PubMed]
23. Zhao W, Niu K, Schafer S, Royalty K. An indirect transmission measurement-based spectrum estimation method for computed tomography. Phys Med Biol 2015;60:339-57. [Crossref] [PubMed]
24. Chang S, Mou X. A statistical iterative reconstruction framework for dual energy computed tomography without knowing tube spectrum. Proc SPIE 2016;99671L.
25. Wu W, Wang Q, Liu F, Zhu Y, Yu H. Block matching frame based material reconstruction for spectral CT. Phys Med Biol 2019;64:235011. [Crossref] [PubMed]
26. Zhao Y, Zhao X, Zhang P. An extended algebraic reconstruction technique (E-ART) for dual spectral CT. IEEE Trans Med Imaging 2015;34:761-8. [Crossref] [PubMed]
27. Xu Q, Mou X, Tang S, Hong W, Zhang Y, Luo T, editors. Implementation of penalized-likelihood statistical reconstruction for polychromatic dual-energy CT. Medical Imaging 2009;7258:72585I.
28. Zhang R, Thibault JB, Bouman CA, Sauer KD, Hsieh J. Model-Based Iterative Reconstruction for Dual-Energy X-Ray CT Using a Joint Quadratic Likelihood Model. IEEE Trans Med Imaging 2014;33:117-34. [Crossref] [PubMed]
29. Yu X, Cai A, Li L, Jiao Z, Yan B. Low-dose spectral reconstruction with global, local, and nonlocal priors based on subspace decomposition. Quant Imaging Med Surg 2023;13:889-911. [Crossref] [PubMed]
30. Mechlem K, Ehn S, Sellerer T, Braig E, Munzel D, Pfeiffer F, Noel PB. Joint Statistical Iterative Material Image Reconstruction for Spectral Computed Tomography Using a Semi-Empirical Forward Model. IEEE Trans Med Imaging 2018;37:68-80. [Crossref] [PubMed]
31. Cai C, Rodet T, Legoupil S, Mohammad-Djafari A. A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography. Med Phys 2013;40:111916. [Crossref] [PubMed]
32. Liu J, Ding H, Molloi S, Zhang X, Gao H. TICMR: Total Image Constrained Material Reconstruction via Nonlocal Total Variation Regularization for Spectral CT. IEEE Trans Med Imaging 2016;35:2578-86. [Crossref] [PubMed]
33. Murata K, Ogawa K. Ring-Artifact Correction With Total-Variation Regularization for Material Images in Photon-Counting CT. EEE Transactions on Radiation and Plasma Medical Sciences 2021;5:568-77.
34. Hu J, Zhao X, Wang F. An extended simultaneous algebraic reconstruction technique (E-SART) for X-ray dual spectral computed tomography. Scanning 2016;38:599-611. [Crossref] [PubMed]
35. Zhao S, Pan H, Zhang W, Xia D, Zhao X. An oblique projection modification technique (OPMT) for fast multispectral CT reconstruction. Phys Med Biol 2021;66:065003. [Crossref] [PubMed]
36. Zhang W, Zhao S, Pan H, Zhao Y, Zhao X. An iterative reconstruction method based on monochromatic images for dual energy CT. Med Phys 2021;48:6437-52. [Crossref] [PubMed]
37. Chang S, Li M, Yu H, Chen X, Deng S, Zhang P, Mou X. Spectrum Estimation-Guided Iterative Reconstruction Algorithm for Dual Energy CT. IEEE Trans Med Imaging 2020;39:246-58. [Crossref] [PubMed]
38. Zhao X, Chen P, Wei J, Qu Z. Spectral CT imaging method based on blind separation of polychromatic projections with Poisson prior. Opt Express 2020;28:12780-94. [Crossref] [PubMed]
39. Shen L, Xing Y, Zhang L. Joint Reconstruction and Spectrum Refinement for Photon-Counting-Detector Spectral CT. IEEE Trans Med Imaging 2023;42:2653-65. [Crossref] [PubMed]
40. Shi Z, Li H, Cao Q, Wang Z, Cheng M. A material decomposition method for dual-energy CT via dual interactive Wasserstein generative adversarial networks. Med Phys 2021;48:2891-905. [Crossref] [PubMed]
41. Wu W, Guo X, Chen Y, Wang S, Chen J. Deep Embedding-Attention-Refinement for Sparse-View CT Reconstruction. IEEE Transactions on Instrumentation and Measurement 2023;72:1-11.
42. Duan Z, Li D, Zeng D, Bian Z, Ma J. A semi-supervised material quantitative intelligent imaging algorithm for spectral CT based on prior information perception learning. Nan Fang Yi Ke Da Xue Xue Bao 2023;43:620-30. [Crossref] [PubMed]
43. Wang W, Xia XG, He C, Ren Z, Lu J, Wang T, Lei B. An End-to-End Deep Network for Reconstructing CT Images Directly From Sparse Sinograms. IEEE Transactions on Computational Imaging 2020;6:1548-60.
44. Li D, Bian Z, Li S, He J, Zeng D, Ma J. Noise Characteristics Modeled Unsupervised Network for Robust CT Image Reconstruction. IEEE Trans Med Imaging 2022;41:3849-61. [Crossref] [PubMed]
45. Zhuang L, Bioucas-Dias JM. Fast Hyperspectral Image Denoising and Inpainting Based on Low-Rank and Sparse Representations. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 2018;11:730-42.
46. Dabov K, Foi A, Katkovnik V, Egiazarian K. Image denoising by sparse 3-D transform-domain collaborative filtering. IEEE Trans Image Process 2007;16:2080-95. [Crossref] [PubMed]
47. Zoran D, Weiss Y. From learning models of natural image patches to whole image restoration. 2011 International Conference on Computer Vision, Barcelona, Spain, 2011:479-86,
48. Zainulina E, Chernyavskiy A, Dylov D. No-reference denoising of low-dose CT projections. 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI), 2021:77-81.
49. Hasan AM, Mohebbian MR, Wahid KA, Babyn P. Hybrid-Collaborative Noise2Noise Denoiser for Low-Dose CT Images. IEEE Transactions on Radiation and Plasma Medical Sciences 2021;5:235-44.
50. Xu Y, Yin W. A Block Coordinate Descent Method for Regularized Multiconvex Optimization with Applications to Nonnegative Tensor Factorization and Completion. SIAM J Imaging Sci 2013;6:1758-89.
51. Byrd RH, Lu P, Nocedal J, Zhu C. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM J Sci Comput 1995;16:1190-208.
52. Zhang S, Geng G, Li Z, Zhang Y. Fast implementation of area integral model SART algorithm based on look-up table. Cluster Computing 2018;22:15195-203.
53. Landry G, Deblois F, Verhaegen F. ImaSim, a software tool for basic education of medical x-ray imaging in radiotherapy and radiology. Front Phys 2013;1:22.
54. Hubbell J, Stephen S. Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients from 1 keV to 20 MeV for Elements Z = 1 to 92 and 48 Additional Substances of Dosimetric Interest*. Available online: http://physics.nist.gov/PhysRefData/XrayMassCoef/cover.html. 1995.
55. Jung B. Book Review: Phantoms and Computational Models in Therapy, Diagnosis and Protection. ICRU Report 48. Acta Radiologica 1994;35:202.
56. Feldle P, Grunz JP, Kunz AS, Pannenbecker P, Patzer TS, Pichlmeier S, Sauer ST, Hendel R, Ergün S, Bley TA, Huflage H. Influence of spectral shaping and tube voltage modulation in ultralow-dose computed tomography of the abdomen. BMC Med Imaging 2024;24:49. [Crossref] [PubMed]
57. Li F, Manerikar AV, Prakash T, Kak AC. A Splitting-Based Iterative Algorithm for GPU-Accelerated Statistical Dual-Energy X-Ray CT Reconstruction. CoRR 2019:abs/1905.00934.
58. Martin L, Karl WC, Ishwar P, editors. Structure-preserving dual-energy CT for luggage screening. 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Florence, Italy, 2014:1175-9.
59. Chen X, Xia W, Yang Z, Chen H, Liu Y, Zhou J, Wang Z, Chen Y, Wen B, Zhang Y. SOUL-Net: A Sparse and Low-Rank Unrolling Network for Spectral CT Image Reconstruction. IEEE Trans Neural Netw Learn Syst 2023; Epub ahead of print. [Crossref]
60. Wu W, Hu D, Niu C, Yu H, Vardhanabhuti V, Wang G. DRONE: Dual-Domain Residual-based Optimization NEtwork for Sparse-View CT Reconstruction. IEEE Trans Med Imaging 2021;40:3002-14. [Crossref] [PubMed]
61. Wang Y, Ren J, Cai A, Wang S, Liang N, Li L, et al. Hybrid-Domain Integrative Transformer Iterative Network for Spectral CT Imaging. IEEE Transactions on Instrumentation and Measurement 2024;73:1-13.
62. Lee W, Son S, Lee KM. AP-BSN: Self-Supervised Denoising for Real-World Images via Asymmetric PD and Blind-Spot Network. 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2022.
63. Song Y, Shen L, Xing L, Ermon S. Solving Inverse Problems in Medical Imaging with Score-Based Generative Models. arXiv 2021:abs/2111.08005.
64. Wu W, Pan J, Wang Y, Wang S, Zhang J. Multi-channel Optimization Generative Model for Stable Ultra-Sparse-View CT Reconstruction. IEEE Trans Med Imaging 2024; Epub ahead of print. [Crossref]
65. Pan J, Yu H, Gao Z, Wang S, Zhang H, Wu W. Iterative Residual Optimization Network for Limited-Angle Tomographic Reconstruction. IEEE Trans Image Process 2024;33:910-25. [Crossref] [PubMed]
66. Zhang J, Mao H, Wang X, Guo Y, Wu W. Wavelet-Inspired Multi-channel Score-based Model for Limited-angle CT Reconstruction. IEEE Trans Med Imaging 2024; Epub ahead of print. [Crossref]