A sequential regularization based image reconstruction method for limited-angle spectral CT

In spectral computed tomography (CT), the object is respectively scanned under different x-ray spectra. Multiple projection data can be collectively used for reconstructing basis images and virtual monochromatic images, which have been used in material decomposition, beam-hardening correction, bone removal, and so on. In practice, projection data may be obtained in a limited scanning angular range. Images reconstructed from limited-angle data by conventional spectral CT reconstruction methods will be deteriorated by limited-angle related artifacts and basis image decomposition errors. Motivated by observations of limited-angle spectral CT, we propose a sequential regularization-based limited-angle spectral CT reconstruction model and its numerical solver. Both simulated and real data experiments validate that our method is capable of suppressing artifacts, preserving edges and reducing decomposition errors.


Introduction
X-ray computed tomography (CT) is a non-destructive imaging technique, which has been widely used in medical diagnosis (Garvey and Hanlon 2002), industrial detection (Chiffre et al 2014) and security inspection (McCollough et al 2015). The reconstructed image is related to the materials of scanned object and the x-ray spectrum. In standard CT, spectral information is ignored, which leads to beam hardening (Brooks and Chiro 1976) and the difficulty in distinguishing some different materials (Riederer and Mistretta 1977). While in spectral CT, the object is respectively scanned under different x-ray spectrum, and more information can be employed to reconstruct material-selective images and energy-selective images Alvarez and Macovski (1976), Vetter et al (1986), Chuang and Huang (1988). In medical applications, spectral CT can be applied to non-destructive imaging with the capability for material decomposition , Niu et al 2019, beam-hardening correction (Coleman andSinclair 1985, Gao et al 2006), bone removal (Morhard et al 2009, Johnson et al 2011, and so on. According to the type of detector, the scanning mode of spectral CT system can be mainly classified into two categories: energy-integrating detector-based scanning, including dual-source CT scanning (Primak et al 2009) and fast kVp switching scanning (Xu et al 2009); and energy-discriminating detector-based scanning, including multilayer detector scanning (Carmi et al 2005) and photon-counting detector scanning (Schlomka et al 2008). Each scanning mode is almost equal to respectively scan the object under different spectra.
Just like standard CT, spectral CT scanning may also encounter the problem of limited-angle imaging. Objectively, some projections in specific views may be unavailable owing to the scanning geometry, or X-rays may be hard to pass through the scanned object due to its high density. Subjectively, the scanning angular range can be narrowed for non-standard scanning configurations, which have the advantages of reducing scanning time and radiation dose, lowering hardware cost, and so on (Chen et al 2017). In spectral CT scanning, when the projection data measured under certain spectrum is angle limited and cannot satisfy complete sampling condition, this scanning faces the problem of limited-angle imaging. Taking circular fan-beam CT as an example, the complete sampling condition is that the effective scanning range covers 180 • plus fan angle (Noo et al 2002). In limited-angle spectral CT, as a result of lacking essential projections to reconstruct the whole unknowns, reconstructed images are degenerated by limited-angle related artifacts and basis material decomposition errors.
Spectral CT reconstruction methods can be basically divided into image space (IS) based methods and projection space (PS) based methods. IS based methods separately reconstruct image from each polychromatic projection data and then linearly combine the reconstructed images to generate basis material images (Brooks 1997, Maaß et al 2009. In the case of limited-angle imaging, limited-angle related artifacts deteriorate the reconstructed images and are hard to be suppressed in the subsequent combination. PS based methods include directly calibration methods and iterative methods. In PS based directly calibration methods, multiple polychromatic projection data are combined to generate the monochromatic projection data of basis material images, then one can directly reconstructed the basis material images from their corresponding projection data (Alvarez andMacovski 1976, Stenner et al 2007). To improve the reconstruction quality, these methods usually require product terms of polychromatic projections measured under different spectra. Thereby the projection data are required to be geometrically consistent, which is difficult to be satisfied in limited-angle spectral CT scanning. PS based iterative methods reconstruct basis material images from polychromatic projection data by iteratively solving a non-linear imaging model (Long and Fessler 2014, Zhao et al 2015, Barber et al 2016. Even though PS based iterative methods are computationally intensive, they have the following advantages: polychromatic projection is naturally taken into account to achieve a more accurate reconstruction; both geometrically consistent and geometrically inconsistent projection data can be used for reconstruction; they are easy to incorporate regularization term, which is a generally used technique to constrain the process of reconstruction. Therefore, our proposed method belongs to PS based iterative methods. The limited-angle spectral CT reconstruction is an ill-posed inverse problem (Natterer 1986), and a widely used method to solve this kind of problem is regularization method: where f i (i = 1, · · · , I) is ith basis image to be reconstructed. D(f 1 ; · · · ; f I ) is data fidelity term, which denotes the discrepancy between the measured polychromatic projection data and the polychromatic projection data of reconstructed basis images. C(f i ) is regularization term, which represents prior information about basis image f i and is used for forcing f i to possess the desired characteristic. The penalty parameter τ is used for adjusting relative strength between data fidelity term and regularization terms. The sparsity of image in specific transformed space is a kind of frequently used prior knowledge (Candes andRomberg 2006, Donoho 2006). The sparsity-exploiting CT reconstruction methods can base on discrete gradient transform (Sidky and Pan 2008), prior image similarity (Chen et al 2008), wavelet tight frame (Zhou et al 2013), dictionary learning (Chen et al 2014), deep learning (Chen et al 2018), and so on. In the application of CT reconstruction, the image to be reconstructed usually has sparse representation on discrete gradient space (Sidky and Pan 2008). Besides, the sparsity based on discrete gradient transform has wide application scope and has no need of any training set. Standard total variation (TV), formulated as the l 1 norm of gradient image, has been applied to CT reconstruction of limited-angle (Song et al 2007), circular cone-beam (Sidky and Pan 2008), low dose (Tian et al 2011), sparse angle (Bian et al 2013), and so on. Nevertheless, standard TV has the following two limitations: it is isotropic in nature, and is incapable of suppressing artifacts in appointed direction; it bases on l 1 norm, and is prone to blur the edges of object. Some modified regularization terms have been developed to improve CT reconstruction methods. On the one hand, anisotropic TV constrains the sparsity of reconstructed image in appointed directions, and leads to higher image quality than its TV counterpart in limited-angle (Chen et al 2013) and exterior problem (Guo et al 2017) CT reconstruction. On the other hand, l 0 gradient sparsity can suppress artifacts and preserve edges in limited-angle (Yu and Zeng 2015, Yu et al 2017a, Xu et al 2019, few view (Yu et al 2017b) and swinging multi-source (Yu et al 2019) CT reconstruction. Even though the above regularization terms are developed for standard CT, they can be generalized to spectral CT as well. For instance, TV term (Chen et al 2015, Chen et al 2017, Yao et al 2019, l 0 gradient sparsity term (Wu et al 2018a, Wu et al 2018b have been applied to limited-angle, sparse angle and low dose spectral CT reconstruction. Differing from standard CT, spectral CT reconstruction is an illness process in nature, which is sensitive to disturbance. Specially, under the circumstance of limited-angle scanning, this illness nature will boost limited-angle related artifacts. Besides, to the best of my knowledge, the existing methods for limited-angle spectral CT have not taken scanning angular range into account. Therefore, it is still a challenging task to be discussed.
In this work, we observe images reconstructed from limited-angle spectral CT, and summarize three observations: first, the artifacts induced by limited-angle scanning are related to scanning angular range and their scales are substantial; second, the artifacts induced by modelling inaccuracy and computational error are independent on scanning angular range and their scales are small; third, the artifacts induced by different limited-angle projection data may have different geometric distributions. These observations will be specified with examples in next section. Motivated by these observations, we present the idea of sequential regularization, namely we sequentially employ multiple regularization terms into iterative process to constrain characteristics of reconstructed images. Specifically, one kind of term is weighted directional gradient sparsity term, which is devised to suppress artifacts in observation 1. Another kind of term is standard TV term, which is employed to suppress artifacts in observation 2. It is unsuited to separately use these terms. Because the first one used alone is prone to generate streak artifacts, and the second one used alone tend to obscure edges of object. Therefore, multiple terms should be cooperatively employed to suppress corresponding artifacts and alleviate the opposing adverse impacts. In observation 3, we discover that the appointed sparsity directions should be separately selected according to the scanning range of different projection data. Based on this idea, we formulate a sequential regularization based limited-angle spectral CT reconstruction model. Then, we develop a numerical solver, namely sequential regularization based PS iterative reconstruction algorithm, to solve this model. Both simulated and real data experiments validate that our method can suppress limited-angle related artifacts, preserve inherent edges and reduce basis material decomposition errors in limited-angle spectral CT.
The remainder of this paper is organized as follows. In section 2, we propose sequential regularization based limited-angle spectral CT reconstruction method. First, we introduce the general spectral CT reconstruction model. Second, limited-angle scanning situations of spectral CT are described. Then, we observe images reconstructed from limited-angle projection data by conventional spectral CT reconstruction method, and present the idea of sequential regularization. Last, we propose a sequential regularization based limited-angle spectral CT reconstruction model and its numerical solver. In section 3, simulated and real data experiments are taken to validate our method. In section 4, we discuss some details of this work. In section 5, we summarize this work.

Spectral CT reconstruction model
In clinical applications, energy-integrating detector based spectral CT scanning systems are widely used (McCollough et al 2015). In this paper, we take these scanning systems (hereinafter referred to as spectral CT systems) for instance to present our method.
Omitting scattered photons and noises and combining with basis material decomposition method , discrete spectral CT reconstruction model can be derived as follows: where I is the number of basis materials, K is the number of energy spectra, L k is the set of x-ray paths corresponding to kth spectrum, p k,l k is the projection measured along path l k . Let M k is the sampling size of kth spectrum, and δ is the sampling interval, then {S k,m } M k m=1 are sampling values of kth normalized effective spectrum, and {ψ i,m } M k m=1 are sampling values of mass attenuation coefficient of the ith selected basis material. r l k = (r l k ,1 , r l k ,2 , . . . , r l k ,J ) denotes projection vector measured along path l k . f i = (f i,1 , . . . , f i,J ) T denotes ith basis material image, namely the density image of the ith basis material. J is total pixel number of each image. And r l k ,j represents the contribution of the jth pixel of f i to the projection along path l k .
Under ideal conditions, spectral CT reconstruction can be boiled down to solve basis images (f 1 ; · · · ; f I ) from the complete non-linear system equation (2). Then, virtual monochromatic image can be approximately estimated by linear combination of basis images . However, inverse problem equation (2) is ill-posed in nature, which will magnify disturbances, such as data incomplete, modelling inaccuracy and computational error. Besides, when projection data are incompletely measured, equation (2) is ill-conditioned and its solution is non-unique. Under this circumstance, it is difficult to correctly reconstruct basis images without any constraint to decrease unknowns. Specially, if the projection data are scanning angle limited, the mistaken reconstruction is characterized by inducing limited-angle related artifacts.

Sequential regularization (SR)
In this subsection, we first briefly introduce two kinds of spectral CT systems, and exemplify their limited-angle scanning situations. Then, we observe images reconstructed from limited-angle spectral CT. Last, we specify the idea of SR to deal with limited-angle spectral CT reconstruction. The mathematical formula of SR is presented with reconstruction model in the next subsection.

Limited-angle scanning of spectral CT
As shown in figure1(a) and (b), fast kVp switching and dual-source CT scanning mode are sketched, and their limited-angle scanning situations are exemplified. In this coordinate system, we define the unit vector ⃗ e α as the reverse direction from x-ray scource to detector, and scanning angle α ∈ [0, 360 • ) is the angle between ⃗ e α and the positive direction of x-axis. As introduced in section 1, if the effective scanning angular range of circular fan-beam CT unable to cover 180 • plus fan angle, this scanning encounter a limited-angle problem.
In these examples, we respectively fix the scanning angular range size and interval of each projection data to 140 • and 0.5 • .
Fast kVp switching CT system changes voltages between each view so that the odd (or even) projections correspond to the low (or high) tube voltage, and the high-and low-energy projection data have one view mismatching (Xu et al 2009). In figure 1(a), a pair of x-ray source and data acquisition system is employed. The kVp alternates between each view, and high-and low-energy projection data have 0.25 • mismatching. The interior and exterior circular ring are respectively represent the scanning angular range of high-and low-energy projection data. The solid arcs represent their limited-angle scanning ranges. In this scanning mode, the global scanning angular range size reaches 140.25 • .
Dual-source CT is a CT system where two pairs of x-ray source and data acquisition system are mounted on the same gantry, positioned orthogonally to each other, and the projection data are approximately 90 • out of phase (Primak et al 2009). In figure 1(b), two pairs of x-ray source and data acquisition system are orthogonally positioned. The interior and exterior circular ring are respectively represent the scanning angular range of low-and high-energy projection data. The solid arcs represent their limited-angle scanning ranges. In this scanning mode, the global scanning angular range size reaches 230 • .

Observations of limited-angle spectral CT reconstruction
To observe the results of limited-angle spectral CT reconstruction, we simulate the limited-angle scanning situations in figures 1(a) and (b) and scan the phantom in figure 1(c). The reconstructed images are respectively shown in figures 1(d) and (e). Experimental configurations are omitted here, and please refer to section 3 for more details. To analyze the influences of each limited-angle projection data, we respectively reconstruct images from them by conventional Algebraic Reconstruction Techniques (ART) (Gordon et al 1970). In this work, Extended Algebraic Reconstruction Technique (E-ART) (Zhao et al 2015), which is a PS based spectral CT iterative reconstruction method and proposed by our team, is taken as an example to solve equation (2).
As shown in figure 1, observations of limited-angle spectral CT reconstruction can be summarized as follows: • observation 1: The artifacts induced by limited-angle scanning are related to scanning angular range and their scales are substantial. In CT, limited-angle artifacts have directional relativity, which is related to scanning angular range. In standard CT, Quinto has introduced two principles that describe features of limited-angle artifacts, and has presented theoretical explanation (Quinto 2017). The first principle is that edges and details tangent to measured projections are more likely to be recovered, otherwise they are difficult to be recovered. The second principle is that streak artifacts can occur around the edges along directions of the end of scanning angular range. As shown in figure 1(d) and (e), the images reconstructed from each limited-angle projection data by ART are deteriorated by limted-angle artifacts, which conform to the aforementioned principles. We take the image reconstructed from high-energy projection data in figure 1(d) as an instance to elaborate on the directional relativity of limited-angle artifacts. In this case, high-energy projection data are measured in scanning angular range [110 • , 250 • ), which means that vertical projections are missing and the scanning angles at the end of range are 110 • and 250 • . In reconstructed image, the vertical edges are blurry, and the streak artifacts occur around the edges along directions of 110 • and 250 • . These phenomena confirm that limited-angle artifacts are related to scanning angular range. In spectral CT, multiple projection data are collectively employed to reconstruct basis images, and the limited-angle artifacts have the similar directional relativity. In figure 1(d), because the scanning angular ranges of high-and low-energy projection data are almost the same, the limited-angle artifacts in basis images and monochromatic image are similar to the ones in standard CT images. In figure 1(e), even though the scanning angular ranges of each projection data have 90 • difference, the artifacts have directional relativity as well. In basis image 1, the diagonal edges and back-diagonal edges are difficult to be recovered due to the absence of diagonal high-energy projections and the absence of back-diagonal low-energy projections. In basis image 2, the streak artifacts exist along directions of 65 • , 205 • , 155 • and 295 • . In addition to directional relativity, the limited-angle artifacts are too severe to be distinguished with inherent image edges, which means that their scales are substantial.
As shown in figure 1(d) and (e), the basis images and monochromatic images are also deteriorated by slight artifats, which are not obviously ralated to scanning angular range. These artifacts mainly come from modelling inaccuracy and computational error, and are magnified by the illness nature of equation (2). Other than limited-angle induced artifacts, this kind of artifacts also damage reconstruction quality, and should be suppressed by proper method. Figure 1(e) also illustrates that the limited-angle artifacts induced by each projection data may have different geometric distributions, and these artifacts impact the spectral CT reconstruction together. In limited-angle dual-source CT scanning, the projection data measured under different spectra have 90 • difference. On the one hand, they can mutually supplement missing information in the reconstruction process. Thus the streak artifacts of basis images and monochromatic image in figure 1(e) does not spread too long. On the other hand, the scanning angular ranges are doubled. Thus the streak artifacts occur along more directions.

The idea of SR
As stated in section 2.1, limited-angle spectral CT reconstruction is an ill-conditioned inverse problem, and regularization term should be incorporated into reconstruction model to decrease unknowns. Motivated by the aforementioned observations and previous works, we present the idea of SR, namely we sequentially employ multiple regularization terms into iterative process to constrain characteristics of reconstructed images.
The first kind of regularization term is weighted directional gradient sparsity term, which is devised to suppress artifacts in observation 1. In this term, scanning angular range is taken as an additional prior information, and we accordingly assign multiple sparsity directions to suppress scanning angular range related artifacts. Besides, this term can suppress substantial artifacts by increasing weights or employing an appropriate norm to achieve sparser gradient image. When suppressing this kind of substantial artifacts, the inherent image edges should be preserved as well.
The second kind of regularization term is standard TV term, which is applied to suppress artifacts in observation 2. This term is isotropic, which means that it has no obvious directivity. Besides, TV is formulated as l 1 norm of gradient image. Comparing to l p (p < 1) norm based regularization term, its effect of sparsity is relatively mild.
In order to explain the effects of these two kinds of regularization terms, we take the image reconstructed from high-energy projection data in figure 1(d) as an example to implement image smoothing, and the results are shown in figure 2. It should be noted that these are simply image post-processing, which are not incorporated into iterative reconstruction. Thus the artifacts can not be completely suppressed. It is illustrated from figure 2 that the single directional gradient minimization is capable of suppressing artifacts in observation 1. But it is prone to generate streak artifacts along sparsity direction. The standard TV minimization is good at suppressing artifacts in observation 2. But it tends to obscure edges of object. To suppress corresponding artifacts and alleviate the opposing adverse impacts, we sequentially employ these two kinds of terms in the whole iterative process. Specifically, the standard TV term should be employed after directional gradient sparsity term to weaken the streak artifacts induced by it. And directional gradient sparsity term should be employed after standard TV term to sharpen the edges. Therefore, we present the sequential regularization in this work.
Besides, as revealed in observation 3, the limited-angle artifacts are related to the scanning angular range of each projection data. Thus, we should separately select sparsity directions for each data.

SR based reconstruction model 2.3.1. Reconstruction model
Based on the idea of SR, we formulate a limited-angle spectral CT reconstruction model: where D(f 1 ; · · · ; f I ) is the data fidelity term, which denotes the discrepancy between the measured polychromatic projection data and the polychromatic projection data of reconstructed basis images. C which represents a prior information of f i . This regularization term bases on sparsity of gradient image, because we assume that the reconstructed images are sparse in their discrete gradient space. The penalty parameter τ is used for adjusting relative strength between data fidelity term and regularization terms.

Data fidelity term
The data fidelity term is as follows: Under kth spectrum, p k is the measured polychromatic projection data, and R k represents the system matrix of scanning, where L k denotes the number of x-ray paths: The data fidelity term is derived from equation (2), thus our model omits scattered photons and noises as it does.

Regularization term
The SR term C (n) SR (f i ) is composed of the following terms: where the first one is directional gradient sparsity term based on direction ⃗ e α d and l p1 (p 1 ∈ [0, 1]) norm. In this work, we select N directions ⃗ e α d , d = 1, · · · , N to constrain the sparsity. ∇ α d denotes discrete directional difference operator along direction ⃗ e α d , and | · | denotes a spatial magnitude operator with regard to single variable. The second one is standard gradient sparsity term based on l p2 (p 2 ∈ [0, 1]) norm. ∇ x and ∇ y respectively denotes discrete directional difference operator along horizontal direction x and vertical direction y, and √ (·) 2 + (·) 2 denotes a spatial magnitude operator with regard to two variables. Weights ω i,d and ω i are respectively used for controlling strength of two terms. A typical choice for l p1 and l p2 are l 0 and l 1 respectively. Where ∥x∥ 0 = ∑ p #{p | |x p | ̸ = 0}, and #{·} is the counting operator, outputting the number of pixel that satisfies |x p | ̸ = 0. And ∥x∥ 1 = ∑ p |x p |. Under the assumption that gradient image is sparse, l 0 norm of gradient magnitude can achieve sparser gradient image and is more powerful on edge preservation than its l 1 norm counterpart (Xu et al 2011). Thus, we employ l 0 norm for C p1,d (f i ) to suppress substantial artifacts in observation 1 and preserve inherent image edges. And we employ l 1 norm for C p2 (f i ) to suppress the slight artifacts in observation 2 and weaken the streak artifacts induced by C p1,d (f i ).
As explained in observation 3, the limited-angle artifacts are related to the scanning angular range of each projection data. Thus, we separately select sparsity directions for each projection data. We take the high-energy projection data, of which the missing projections range from 70 • to 110 • , as an example to explain how to select sparsity directions. To recover the blurry edges, one should constrain sparsity along the directions of angles in [70 • , 110 • ]. To suppress streak artifacts, one should constrain sparsity perpendicular to the directions of 70 • and 110 • . In practice, in order to simplify the implementation and fully utilize information of the 8-neighborhoods of pixels, we preset four designative directions ⃗ e 0 , ⃗ e π/4 , ⃗ e π/2 , ⃗ e 3π/4 . Then, we select N directions ⃗ e α d , d = 1, · · · , N from them to approximate the above sparsity directions. In this example, we select directions ⃗ e π/4 , ⃗ e π/2 , ⃗ e 3π/4 to cover the missing range [70 • , 110 • ], and select directions ⃗ e 0 , ⃗ e π/4 , ⃗ e 3π/4 to approximate directions perpendicular to the directions of 70 • and 110 • . Thus, the selected sparsity directions are ⃗ e 0 , ⃗ e π/4 , ⃗ e π/2 , ⃗ e 3π/4 . In most cases, the sparsity directions are ⃗ e 0 , ⃗ e π/4 , ⃗ e π/2 , ⃗ e 3π/4 , but the effect of them may be different. Some of them are used for recovering blurry edges, and some of them are used for suppressing streak artifacts. When sparsity directions are determined, the weights ω i,d , d = 1, · · · , N are used for controlling strengths of corresponding directional gradient sparsity terms.
As shown in figure 3, discrete directional difference operators employed in this work are expressed as follows: As explained in subsection 2.2.3, regularization terms C 0,d (f i ) and C 1 (f i ) should be sequentially employed in the whole iterative process to suppress corresponding artifacts and alleviate the opposing adverse impacts. The principle to select regularization term in iteration is as follows: we first employ one kind of term to suppress its corresponding artifacts (such as we employ C 0,d (f i ) to suppress artifacts in observation 1). In the next iteration, we change another term to suppress the other artifacts (such as we employ C 1 (f i ) to suppress artifacts in observation 2) and mitigate the damage induced by previous regularization term. Generally, we select one term in each iteration, and sequentially employ different terms in iterative process. Given the scanning angular range, we can select N directions ⃗ e α d , d = 1, · · · , N. Then we equally allocate C 0,d (f i ) into iterations. In adjacent iterations, we employ C 0,d (f i ) and C 1 (f i ) respectively. For example, the regularization term sequence could be The difference between SR and regularization with total terms is that SR can flexibly employ different regularization terms in iterative process. These regularization terms can cooperatively suppress artifacts and mutually restrict each other. Besides, if we incorporate too much regularization terms into each iteration, the data fidelity term is easier to be neglected. The SR can help us to keep a balance between data fidelity term and regularization terms.

The E-ART+SR algorithm
To solve reconstruction model equation (3), we present a SR based PS iterative reconstruction algorithm, referred to as E-ART+SR. First, we substitute the non-linear data fidelity term equation (4) with a linear one by E-ART (Zhao et al 2015). Then, we employ an alternating optimization strategy to approximately solve this optimization problem.

The linearization of data fidelity term
In this work, we take E-ART (Zhao et al 2015) for example to solve the non-linear spectral CT reconstruction model equation (2), which forms the data fidelity term equation (4). In E-ART, the number of basis materials I and the number of energy spectra K are both fixed to 2. E-ART iteratively solves equation (2) by first-order Taylor approximation at n th iteration point (f where q (n) Then, the projection data p k , can be linearly approximated as where Thus, the non-linear data fidelity term equation (4) can be linearised in iterative process: The optimization model equation (3) can be rewritten as follows:

Iterative algorithm
The original optimization problem equation (14) is hard to be solved due to the presence of l 0 norm in SR term. Similarly to (Yu and Zeng 2015), we employ alternating optimization strategy to approximately solve it. In nth iteration, equation (14) is split into two subproblems: Equation (15) is equal to spectral CT reconstruction by gradient descent update with step size λ (n) , and this subproblem is solved by E-ART in this work.
Equation (16) is SR minimization with similarity constraint. In its objective function, basis images are mutually independent, thus it can be divided into two identical problems: where λ * = 2λ (n) · τ . Because λ * will multiply with weight ω i,d or ω i , we set it to 1. Equation (17)  i . We employ the algorithm framework in reference (Xu et al 2011) to solve it. This algorithm framework is applied to l 0 gradient minimization, so two modifications must be made. First, we extend it to directional l 0 gradient minimization by replacing the difference operator ∇ = (∇ x , ∇ y ) T with directional difference operator ∇ α d , and the subsequent operations are taken on direction α d . Second, we extend it to l 1 gradient minimization by replacing the element-wise hard threshold shrinkage operator with element-wise soft threshold shrinkage operator (Elad 2010). Except for weights ω i,d and ω i , the parameters of algorithm framework stay the same. The critical parameters ω i,d and ω i are used for controlling the level of smoothing.
The pseudo-code for the proposed E-ART+SR algorithm is listed in algorithm 1: Scanning views permutation (Mueller et al 1997): get a list of views (v0, v1, . . . , v views−1 ) 3 for Spectral CT reconstruction (Zhao et al 2015) from projections measured under view v 5 Images updating: g (n) 1 , g 2 6 end 7 SR minimization with g (n) 1 , g 2 by framework in reference (Xu et al 2011) 8 Images updating: f In the input stage of this algorithm, we specify main iteration times to N ite . Parameters 2 are used for spectral CT reconstruction. Where limited-angle projection data p k are measured under kth spectrum. Sampling values of kth normalized effective spectrum S k,m (m = 1, · · · , M k ) are simulated with software (SpectrumGUI 2014) or estimated by algorithm . Sampling values of mass attenuation coefficient of the ith basis material ψ i,m (m = 1, · · · , M k ) are retrieved from NIST (Hubbell and Seltzer 2004). System matrix R k are calculated by scanning configuration. The initial images f 2 are both set to 0. Parameters α d , ω i,d , ω i are used for SR minimization. The directions α d , d = 1, · · · , N are determined by scanning angular range. The weights ω i,d , d = 1, · · · , N and ω i are problem dependent, which are discussed in the next section.
This algorithm is implemented in iterative process, and is divided into three steps. First, we reorder the scanning views of projection data for speeding up (Mueller et al 1997). Second, we carry out once PS based spectral CT iterative reconstruction (Zhao et al 2015) with sorted views, and update the intermediate images 2 ). Third, we employ algorithm framework in (Xu et al 2011) to solve SR minimization with similarity constraint. Then the basis images (f ) are updated for the next iteration until iteration termination condition is reached.

Experimental results
In this section, we perform numerical experiments with both simulated and real data to investigate the performance of our proposed E-ART+SR method for limited-angle spectral CT reconstruction. First, we employ simulated circle phantom to validate the effectiveness of E-ART+SR. Second, we introduce principles to select parameters. Third, we employ simulated thorax phantom to compare E-ART+SR with two other methods. Fourth, we employ a bone submerged into water as an authentic specimen to test the feasibility of E-ART+SR for practical applications.
Unless otherwise specified, two limited-angle scanning situations described in figure 1(a) and (b) are simulated. In fast kVp switching scanning mode, high-and low-energy scanning angular range are [110 • , 250 • ) and [110.25 • , 250.25 • ) respectively. In dual-source CT scanning mode, high-and low-energy scanning angular range are [65 • , 205 • ) and [155 • , 295 • ) respectively. The scanning angular range size and interval of each projection data are 140 • and 0.5 • respectively. To simulate limited-angle scanning of spectral CT, we fully scan the object with high-and low-energy spectrum, and extract the projections in designated angle.
In this work, we employ visual comparison or combination of visual comparison and quantitative analysis to verify the performance of E-ART+SR on artifacts suppression, edges preservation and basis material decomposition. To quantitatively analyze the performance, three metrics are used: peak signal to noise ratio (PSNR), root mean square error (RMSE), and feature similarity (FSIM) index (Zhang et al 2011). There is no authentic reference images in real data experiments, thus quantitative analysis does not apply to them.
In comparing experiments, E-ART+SR is compared with E-ART and E-ART+TV methods. E-ART is employed to illustrate the deteriorated results without any regularization item, and its reconstructed images can be regarded as lower bounds. TV is a widely used regularization item in CT reconstruction, thus we combine it with E-ART to make a comparison with our proposed SR item. The parameters of regularization based methods are determined by trial and error. The principles of parameters selection are introduced in subsection 3.2.
All experiments are implemented on a 2.40 GHz Intel® Xeon® CPU processor and a NVIDIA Quadro K2200 graphics card with 4 GB global memory and 640 CUDA cores. The reconstruction codes are written in C language and Compute Unified Device Architecture (CUDA 2020). And the quality analysis codes are written in Matlab language.

Validation of E-ART+SR
In this subsection, we validate the effectiveness of E-ART+SR by mechanism verification, convergence test and multiple scanning congifurations test.
As shown in figure 4(a), we decompose the circle phantom into bone density image and water density image, and calculate the virtual monochromatic image at 70 keV. The resolutions of reference images and reconstructed images are 512 × 512. We simulate the x-ray spectra of GE Maxiray 125 x-ray tube voltages 80 and 140 kV, where the latter is filtered with 1 mm copper. The detector consists of 512 bins with length 1.345 mm. The source-object distance is 1000 mm and the object-detector distance is 300 mm.
Polychromatic projection data are calculated by equation (2). Where the mass attenuation coefficients of basis materials are retrieved from NIST (Hubbell and Seltzer 2004). The polychromatic x-ray spectra are simulated with the software (SpectrumGUI 2014). In this subsection, the polychromatic projection data are noise-free, which means that there is no inconsistency in the projection data. And the phantom is piecewise constant, which means that it has sparse gradient image and the regularization term can exactly describe its sparsity. In these ideally simulated experiments, we can expect the data residuals to decrease to zero. Comparing to standard CT, spectral CT reconstruction model is solved for double or more unknowns, and it is ill-posed in nature, which will magnify disturbances. Thus, sufficient main iteration times is necessary to investigate the convergence behavior of E-ART+SR, and ensure it obtain stable solution. We fix main iteration times to 5000.

Mechanism verification and convergence test
We observe the effect of E-ART+SR in iterative process of limited-angle spectral CT reconstruction, and test its convergence.
As shown in figure 4(b) and (c), E-ART+SR gradually suppresses limited-angle related artifacts in iterative process, and it is capable of preserving inherent edges. At the end of iteration, basis images are exactly decomposed without visible artifacts. Reconstructed basis images and virtual monochromatic images can recover to reference images.
As shown in figure 4(d) and (e), RMSEs of reference images and reconstructed intermediate images are discribed to test the convergence of E-ART+SR. It is illustrated that E-ART+SR can obtain stable solution, which approaches the authentic images if proper parameters are chosen. Besides, focuses are worthy to be given to figure 4(e). First, as a result of increasing global scanning range, dual-source CT scanning is easier to reconstruct high quality images from limited-angle projection data than fast kVp switching scanning. Second, because we sequentially employ different regularization terms in whole iterative process, the basis images reconstructed in different iteration are smoothed by different minimization methods. Thus, the convergence curves of basis images vibrate with iteration. The virtual monochromatic image is a linear combination of basis images, and its convergence curve vibrates as well. This phenomenon is obvious in early iterations, because of the incomplete materials decomposition. In some extent, these vibrations can verify the effectiveness of employing regularization terms in sequence.

Multiple scanning configurations test
In this subsection, we test E-ART+SR with more scanning configurations. First, we supplement experiments with 90 • and 120 • scanning angular range sizes to observe the effects of E-ART+SR in inferior scanning situations. Then, we fix the scanning angular range size to 100 • and rotate the scanning directions to illustrate the necessity to select sparsity directions for each projection data. As shown in figure 5, the reconstructed images are compared. To better observe the details in regions of interest (ROIs), we zoom in on them and display them in the upper left of reconstructed images.
As shown in figure 5(a), both fast kVp switching and dual-source CT scanning mode are tested to investigate the ability of E-ART+SR in inferior scanning situations.
In configuration 1, the scanning range size is 90 • , and the projections in range [45 • , 135 • ) are missing. Comparing to 140 • scanning range size, the edges tangent to the missing projections are increased. Therefore, the images reconstructed by E-ART are seriously damaged. Even though E-ART+SR compensates the missing information by image smoothing along veritical, diagonal and back-diagonal directions, it is incapable of recovering the vertical edges. Becase the projections near to vertical direction are all missing, the ability of information compensation is limited in this case.
In configuration 2, the scanning range size is 90 • , and high-and low-energy projection data have no intersection. What makes the matter worse is that the end of scanning angular ranges of each projection data are coincided. As stated in subsection 2.2.2, the streak artifacts appear along vertical and horizontal directions. In this case, the streak artifacts induced by each limited-angle projection data are mutually reinforced. Thus, the inner circles reconstructed by E-ART are seriously damaged. E-ART+SR can suppress the artifacts in some extent, but it is unable to completely recover the edges of inner circles due to the serious deteriorations.
In configuration 3 and 4, the scanning range sizes are 120 • . Comparing to configuration 1 and 2, the missing ranges are decreased. Thus the limited-angle induced artifacts in images reconstructed by E-ART are attenuated. In configuration 3, the bone density image can be well reconstructed. And the water density image can also be well reconstructed except for the vertical edges in the right of outer circle. The projections passing through the right of outer circle are scarce. Thus the information of blurry edges in the right is harder to be recovered. In configuration 4, E-ART+SR can suppress artifacts and perfectly recover the image structures.
As shown in figure 5(b), the scanning angular range sizes are fixed to 100 • . Then, we rotate the scanning directions to observe the effect of different sparsity direction in SR, and testify the statement in subsection 2.3.3 that the sparsity directions should be selected for each projection data. E-ART is employed to observe the deteriorated reconstruction results. Then, we merely select sparsity directions for high-energy projection data to implement the incomplete E-ART+SR. This can illustrate that if sparsity direction for low-energy projection data is absence, the artifacts induced by it are lack of constraint. To better display the necessity of sparsity directions for low-energy projection data, we decrease all strengths of regularization terms. Last, we carry out E-ART+SR with fine-tuned parameters to explain that it can perfectively reconstruct images in these two cases.
In configuration 5, the high-energy projections in [5 • , 85 • ) are missing, and the low-energy projections in [95 • , 175 • ) are missing. In images reconstructed by E-ART, the diagonal edges are damaged by high-energy scanning, and the back-diagonal edges are damaged by low-energy scanning. We select ⃗ e 0 , ⃗ e π/4 , ⃗ e π/2 for high-energy projection data, and reconstruct images by incomplete E-ART+SR. As shown in figure 5(b), the diagonal edges can be recovered by sparsity direction ⃗ e π/4 . Whereas the back-diagonal edges are hard to be recovered due to the absence of sparsity direction ⃗ e 3π/4 , which should has been selected for low-energy projection data.
In configuration 6, the high-energy projections in [320 • , 40 • ) are missing, and the low-energy projections in [50 • , 130 • ) are missing. In images reconstructed by E-ART, the horizontal edges are damaged by high-energy scanning, and the vertical edges are damaged by low-energy scanning. We select ⃗ e 0 , ⃗ e π/4 , ⃗ e 3π/4 for high-energy projection data, and reconstruct images by incomplete E-ART+SR. As shown in figure 5(b), the horizontal edges can be recovered by sparsity direction ⃗ e 0 . Whereas the vertical edges are hard to be recovered due to the absence of sparsity direction ⃗ e π/2 , which should has been selected for low-energy projection data. The experiments in configuration 5 and 6 can verify the necessity for selecting sparsity directions for low-energy projection data. Thus we need to select sparsity directions for each projection data.

Parameters selection
In this subsection, we introduce principles to select parameters for E-ART+SR. First, we present the possible ranges for parameters, and explain how to narrow them. Then, we introduce the method for fine-tuning. Last, we recommend how to select parameters in practical application.
In E-ART+SR, two kinds of parameters should be predetermined, namely weights ω i,d , d = 1, · · · , N and ω i . Where i corresponding to the ith basis image. Because the basis images have different image features, such as image structure and image contrast, these two kinds of parameters should be selected for each image. In our experiments, we unify the directional weights for each image. Then the parameters are decreased to 1D weight for basis image 1 (refered to as ω 1,1d ), 1D weight for basis image 2 (refered to as ω 2,1d ), 2D weight for basis image 1 (refered to as ω 1,2d ) and 2D weight for basis image 2 (refered to as ω 2,2d ).
In our experiments, basis material 1 is bone, and basis material 2 is water. Comparing to water density image, bone density image usually has sparser gradient image. Thus, the weights for bone density image are no less than the ones for water density image. The possible range for ω 1,1d is [10 −4 , 10 −2 ], and the possible range for ω 2,1d is [10 −5 , 10 −3 ]. Besides, the effect of TV minimization is weaker than the effect of directional l 0 gradient minimization. Thus the 2D weights are no less than their corresponding 1D weights. The possible ranges for 2D weights ω i,2d (i = 1, 2) are [ω i,1d , 2ω i,1d ]. We first narrow the search ranges, and then fine-tune the parameters in the narrowed ranges.
To narrow the search ranges, we divide the possible ranges for 1D weights into multiple test values. The range for ω 1,1d is divided into 10 −4 , 5 × 10 −4 , 10 −3 , 5 × 10 −3 , 10 −2 . The range for ω 2,1d is divided into 10 −5 , 5 × 10 −5 , 10 −4 , 5 × 10 −4 , 10 −3 . We temporarily set ω i,2d = 2ω i,1d (i = 1, 2) for convenience. Then, we can search the parameters in finite cases. As shown in figure 6, we take the 140 • fast kVp switching scanning as an instance to explain how to narrow the search ranges. This process is divided in two steps. In step 1, we fix ω 2,1d to its middle value 10 −4 , and successively test the values of ω 1,1d . The test value numbers 1, 2, 3, 4, 5 respectively corresponding to test values 10 −4 , 5 × 10 −4 , 10 −3 , 5 × 10 −3 , 10 −2 . We assess these values by reconstructed images and metrics. It is illustrated from metrics that the basis images reconstructed with the third test value have lowest RMSEs and highest FSIMs, which means that third test value 10 −3 is the best one among five test values. As shown in reconstructed images, when ω 1,1d is smaller than its proper value, the deteriorated vertical edges still be blurry. When ω 1,1d is larger than its proper value, not only the blurry edges can not be recovered, but also several feak structures are generated. These feak structures have high-contrast and their intensities may exceed the normal range. Thus, the strength of ω 1,1d can be distinguished by visual observation. After determining ω 1,1d = 10 −3 , we carry out step 2 to select ω 2,1d from its test values. The test value numbers 1, 2, 3, 4, 5 respectively corresponding to test values 10 −5 , 5 × 10 −5 , 10 −4 , 5 × 10 −4 , 10 −3 . It is illustrated from metrics that the basis images reconstructed with the second test value 5 × 10 −5 have lowest RMSEs and highest FSIMs. As shown in the lower-left water density image, when ω 2,1d is smaller than its proper value, the vertical edges in the right are blurry. Visually, the basis images reconstructed with second and third test values all have the accurate image structures. But the basis images reconstructed with second value have lower RMSEs. When ω 2,1d is larger than its proper value, the edges in ROIs are distorted. The larger ω 2,1d is, the worse distortion is. Thus, the strength of ω 2,1d can be observed by vision. The parameters ω 1,1d = 10 −3 and ω 2,1d = 5 × 10 −5 are the best ones in this example. If the final images and metrics can satisfy the need for application, the paramters are all determined. Otherwise, we continue to select paramters by fine-tuning. If there are several test values have approximate effects in step 1, we carry out step 2 for each Figure 6. An example to illustrate the process of narrowing search ranges of parameters. In step 1, we fix ω 2,1d and select ω 1,1d . From left to right, the ω 1,1d is increased. In step 2, we fix ω 1,1d and select ω 2,1d . From left to right, the ω 2,1d is increased. Reconstructed images and their ROIs are compared, and the metrics are described. Display windows of bone density image, water density image are [2.0,3.0], [0.4,1.0] respectively. of them. If there are several test values have approximate effects in step 2, we select the range between them. Thus, the search ranges are narrowed by excluding wrong values.
If the parameters are not yet determined, we employ fine-tuning to further select them. We first select 1D weights and fix corresponding 2D weights as above. The ranges for ω i,1d , (i = 1, 2) are averaged into several values. Just as the previous steps, we fix ω 2,1d to its middle value, and test the possible ω 1,1d by metrics and visual observation. Then, ω 1,1d can be determined or its range can be narrowed. We fix ω 1,1d to the best value or the middle value, and select ω 2,1d . These steps are repeated until the final parameters can satisfy the practical need or the ranges are narrow enough. In the latter case, we further select 2D weights. Given ω i,1d , we employ the aforementioned method to select ω i,2d . In the end, the parameters are all determined.
In practice, the reference images are absence, and we can not select parameters by metrics. Furthermore, parameter fine-tuning is a time-consuming process. Thus, parameter selection is an open problem. The parameters are related to several factors, such as scanning angular range, image structure, ROIs. We can simulate the scanning configuration and imaging subject to select parameters (Sidky and Pan 2008).

Experimental study on simulated data
In this subsection, we compare our proposed E-ART+SR with two other methods, and assess the reconstruction qualities by visual comparison and quantitative analysis.
As shown in figure 7, we decompose a thorax phantom into bone and water density images, and calculate the virtual monochromatic image at 70 keV. The resolutions of reference images and reconstructed images are 512 × 512. In reference images, some ROIs and profile lines are selected and marked for reconstruction quality assessment. The scanning configurations are the same as ones in subsection 3.1. Both noise-free and noisy data are tested, and the noisy projection data is generated by adding Poisson noise with 10 7 quanta number. Even though noises are added into projection data, the inconsistencies are smaller than real data experiments. Thus, the main iteration times still be fixed to 5000 to ensure all algorithms obtain their stable solutions.

Fast kVp switching cases
In 140 • fast kVp switching scanning mode, reconstructed images, their quality metrics and their profiles are given in figure 8 and 9.
It is illustrated from reconstructed images that E-ART is incapable of suppressing limited-angle related artifacts and reducing material decomposition errors. The deterioration becomes serious when noises are added into projections. Even though E-ART+TV can improve reconstruction quality in some extent, it is not good at keeping balance of artifacts removal and edges preservation. In water density images reconstructed by E-ART+TV, the limited-angle related artifacts still exist but the black circles in W-3 and W-4 undergo excessive smoothing. The images reconstructed by E-ART+SR are visually better than ones reconstructed by E-ART+TV on both artifacts suppression and edges preservation. Besides, we have noticed that the details of object in the right of images are hard to be recovered due to deficiency of too much projections. On the contrary, details passed through by sufficient projections, such as left part of object, are easier to be reconstructed.
PSNRs, RMSEs and FSIMs in selected ROIs are compared to quantitatively analyze the performances of different methods. Images reconstructed by E-ART+SR have the lowest RMSEs, the highest PSNRs and FSIMs in ROIs, which illustrate that E-ART+SR outperforms other two methods on quantitative analysis.
Image profiles are depicted to make further comparison. In noise-free case, E-ART+SR can exactly recover information of image profiles. In noisy case, even though our method unable to exactly recover all information, it outperforms other methods on noises reduction and edges preservation. When noises are added into projection data, some inherent edges are attenuated together with the magnified aitifacts. Because SR alternatingly employs TV term, and TV term is prone to blur edges. If we aim at suppressing the serious artifacts, the regularization terms may be slightly overused.

Dual-source CT cases
In 140 • dual-source CT scanning mode, reconstructed images, their quality metrics and their profiles are respectively given in figure 10 and 11.
Comparing to images reconstructed by E-ART in fast kVp switching scanning mode, the reconstructed images in dual-source CT scanning mode have different features. On the one hand, the streak artifacts not spread too long. On the other hand, the edges in W-1 suffer from more serious distortions, and material decomposition errors in W-4 are strengthened. The reconstruction qualities are further reduced due to introduction of noises. E-ART+TV are capable of suppressing artifacts in some extent. But it is difficult to keep balance of artifacts suppression, edges preservation and decomposition errors reduction. The distortions of edges still exist in W-1, and the material decomposition errors still exist in W-4. These denote that the strength of regularization term is inadequate. Whereas blurring effects appear in W-3, which means that the regularization term is overused. E-ART+SR can suppress artifacts, reduce material decomposition errors and preserve edges at the same time. Different with fast kVp switching scanning mode, global scanning range are naturally increased in this mode, thus the details of object in the right of images are easier to be reconstructed.
PSNRs, RMSEs and FSIMs in selected ROIs are compared to quantitatively analyze the performance of different methods. In noise-free case, images reconstructed by E-ART+SR have lowest RMSEs, highest PSNRs and FSIMs in ROIs. In noisy case, even though PSNR and RMSE of images reconstructed by E-ART+SR in W-2 are inferior to E-ART+TV counterpart, the differences are quite slight. And from the overall perspective, E-ART+SR outperforms E-ART+TV in the majority.
Image profiles are depicted to make further comparison. In noise-free case, each method can approximately recover information of image profiles. Among them, E-ART+SR recovers edges better. In noisy case, E-ART+SR can suppress artifacts and basically recover edges. But the inherent edges may be attenuated together with the magnified aitifacts due to introduction of TV term in SR.

Experimental study on real projection data
In this subsection, we test the feasibility of E-ART+SR for practical applications. As shown in figure 12, the specimen is a real bone submerged into a cup of water. An x-ray source (YXLON Y.TU450 D09 tube) is operated at the tube voltage of 90 kV and 140 kV for low-and high-energy spectra scanning, and the tube current is 5 mA and 3 mA respectively. The source spectrum estimation please refer to , and the mass attenuation coefficient please refer to (Hubbell and Seltzer 2004). The detector consists of 800 bins with length 0.083 mm. The source-object distance is 677.107 mm, and the object-detector distance is 222.563 mm.
We choose bone and water as basis materials, and reconstruct reference images from full angle projection data by E-ART. The reconstruction resolution is 1024 × 1024. In real data experiments, several factors may lead to inconsistency in spectral CT reconstruction model, such as noises, scattered photons, inaccurate spectrum estimation, inaccurate geometric calibration, non-sparse gradient images, and so on. In these cases, a simplified reconstruction model may reconstruct inauthentic images. Thus, too many iteration numbers are unsuitable. In these experiments, we first reconstruct reference images from full angle projection data by E-ART, and find that 160 times iterations are enough to properly decompose basis materials. Then, we fix the iteration times to 160, and try to recover reference images from limited-angle imaging.
As shown in figure 12, the images reconstructed in fast kVp switching and dual-source CT scanning mode are respectively compared to test our method in practical application. The images reconstructed by E-ART are damaged by obvious limited-angle related artifacts and material decomposition errors. E-ART+TV has not yet suppress artifacts well but obscure edges of object. This confirms that E-ART+TV is unable to keep a good balance between artifacts suppression and edges preservation. In dual-source CT scanning mode, the global scanning range is enlarged, thus E-ART and E-ART+TV can reconstruct fine virtual monochromatic images. But their basis material decompositions are mistaken. E-ART+SR can suppress artifacts, preserve edges and reduce the decomposition errors. It means that our method surpasses other two methods on improving reconstruction quality in practical application.

Discussions
In this section, we disscuss several details about this work.

Convergence behavior of E-ART+SR
The convergence behavior of l p , p ∈ [0, 1) norm involved regularization method is not an unique problem induced by our proposed SR. In this work, we take l 0 norm as an instance to achieve sparser gradient images. l 0 norm based gradient sparsity has been widely used in limited-angle (  this work is identical to abovementioned works. First, the motivations for employing l 0 gradient sparsity regularization term are all suppressing artifacts or denoising, as well as preserving inherent image edges. Second, the reconstruction models are all divided into several subproblems, including l 0 gradient sparsity based subproblem, which has the form as where g (n) is an intermediate image calculated by other step or steps, ∇ denotes discrete difference operator or discrete directional difference operator, mag(·) calculates the magnitude. As all known, l 0 gradient minimization is a non-convex and non-deterministic polynomial hard (NP-hard) problem. Thus, we and all the above authors employ framework in reference (Xu et al 2011) to approximately solve the equation (18)-like subproblems. This framework adopts a special alternating optimization strategy to expand the original l 0 involved terms and update them iteratively. It makes the problem easier to tackle and upholds the property to maintain and enhance salient structures.
Just as abovementioned authors did, we carry out numerical experiments to demonstrate the convergence behavior of our proposed algorithm in this work. In fact, it is difficult to provide a theoretical  characterization of the behavior of our algorithm. Because l 0 norm involved regularization lacks proper regularity properties, e.g. smoothness and convexity, which are usually required for carrying out theoretical analysis (Xu et al 2019). Even the proof of local minimization of l 0 norm involved reconstruction algorithm is an important problem. In reference , authors analyze the convergence of their l 0 norm involved algorithm and prove that the sequence is convergent to a local minimization under certain conditions. If necessary, the theoretical analysis of our algorithm can be regarded as a further work to be studied later.

Computational time
In our numerical implementation, E-ART+SR algorithm is mainly divided into three modules: polychromatic forward projection and residuals calculation; backprojection and residuals allocation; basis images smoothing. If we omit the time for memory copy and scanning views permutation, which is relative little, the computational time of this algorithm can be calculated by summing the computational time of three modules.
In each main iteration, we carry out first two modules for each projection data, and then implement image smoothing for each basis image. The average computational time in each main iteration is derived from experiments in this work, and is presented in table 2. The computational time of each module is positively correlated to scanning angular range size and reconstruction resolution. It should be noticed that we implement this algorithm in Compute Unified Device Architecture (CUDA), thus these positively correlations are non-linear.
Given the main iteration times, the total computational time is calculated by multiplying main iteration times and average computational time together. In simulated data experiments, the main iteration times is 5000. When the scanning angular range size of each projection data is 140 • and the reconstruction resolution is 512 × 512, the computational time approaches 2750 s. In real data experiments, the main iteration times is 160. The scanning angular range size of each projection data is 140 • and the reconstruction resolution is 1024 × 1024. Thus, the computational time approaches 227 s.
In simulated data experiments, there is few inconsistency in reconstruction model. Thus, we fix the main iteration times to a large number to investigate the convergence behavior of our proposed method, and ensure algorithm obtain stable solution. It is a common problem that regularization based iterative reconstruction methods usually have many iteration times. For example, the work (Xu et al 2019) tests the convergence of method with 1000th iterations, and the work (Guo et al 2017) fixes the main iteration times to 3000. These works are based on standard CT rather than spectral CT. In spectral CT, the reconstruction model is solved for double or more unknowns, and it is ill-posed in nature, which will magnify disturbances. Therefore, we increase the main iteration times and fix it to 5000. It should be noticed that the main iteration times to ensure stable solution is problem dependent. In our simulated data experiments, if the scanning angular range size of each projection data is no less than 100 • , the dual-source scanning mode can achieve well decomposed artifacts-free basis images within 1000th iterations. Given 5000th main iterations, each simulated data experiment usually needs tens of minutes, which is really time-consuming. But the key point of this work is to devise regularization term. And this work can be regarded as an extension of spectral CT reconstruction algorithm on limited-angle imaging. Thus, the computational time of our method is dependent on basic spectral CT reconstruction algorithm. As shown in table 2, the time for image smoothing is about 8%−14% in total computational time. As stated in section 1, intensive computation is one of the major disadvantages of projection space (PS) based iterative reconstruction methods. But other kinds of spectral CT reconstruction methods are unsuitable for limited-angle imaging. In order to reduce computational time of this work, efforts should be made on researching basic PS based spectral CT iterative reconstruction methods. This is a valuable reseach direction, and will benefit for many CT applications.
Besides, in real data experiments, their main iteration times are 160, which reduces the computational time to several minutes. Even though improvements must be made for clinical applications, the feasibility of our method is tested by experiments and the computational time will not be a tricky question. Or at least reducing computational time can be regarded as another reseach direction, which can be seperated with this work.

Conclusions
In this work, we focus on limited-angle spectral CT reconstruction. Motivated by three observations, we first present the idea of sequential regularization. Then, a SR based limited-angle spectral CT reconstruction model is proposed and its numerical solver is developed. Last, we investigate the performance of our proposed E-ART+SR method with simulated and real data experiments. The experimental results illustrate that E-ART+SR surpasses E-ART and E-ART+TV on artifacts suppression, edges preservation and decomposition errors reduction.
The validity of our method mainly from two factors. First, limited-angle related artifacts are classified into two types, which are respectively suppressed by different regularization terms. Second, we sequentially employ different regularization terms in iterative process to suppress corresponding artifacts and alleviate the opposing adverse impacts. The key point of this work is devising regularization term. On the one hand, our work can be regarded as an extension of spectral CT reconstruction algorithm on limited-angle imaging. On the other hand, the idea of SR is also suitable for other CT systems, such as energy discriminating detector based spectral CT, standard CT, cone-beam CT and so on.
In future work, three improvements could be taken. First, our reconstruction model is based on some assumptions, such as sparsity of gradient image, absence of scattered photons and noises. In some practical applications, these assumptions may not be satisfied. Even though real data experiments have tested the feasibility of our method, the assumptions should be adjusted or expanded according to practical application scenarios. Second, our method is a regularization method, which relies on some parameters. In this work, we select parameters by some principles. Some algorithms to automatically select parameters can be combined with our method. Third, our work is dependent on PS spectral CT iterative reconstruction method, some reconstruction method to accelerate convergence is of great interests.