Statistical iterative spectral CT imaging method based on blind separation of polychromatic projections

Xiaojie Zhao; Yihong Li; Yan Han; Ping Chen; Jiaotong Wei

doi:10.1364/OE.456184

1. Introduction

X-ray computed tomography (CT) can non-destructively acquire the internal features of objects. Currently, it has become an irreplaceable non-destructive testing method for medical and industrial applications [1–3]. However, beam hardening artifacts (BHA) are widely present in CT images obtained using traditional reconstruction algorithms (e.g., filtered backprojection and algebraic reconstruction technique) because X-ray beams generated by the conventional CT system are continuous and polychromatic, whereas the traditional CT reconstruction algorithm is based on the monochromatic assumption [4,5]. When such algorithms are used to reconstruct polychromatic projections, CT images will appear as cupping and banding artifacts because the beam hardening process is neglected [5]. BHA causes inconsistencies in the linear attenuation coefficient (LAC) of the material, affecting the non-destructive evaluation and quantitative analysis of imaged objects [6,7]. Although numerous BHA correction methods (e.g., linearization) are used for conventional CT, the artifacts are usually not fully corrected [6,8]. They may also destroy the energy characteristics of LAC. Spectral CT has been employed to reduce BHA by exploiting multiple different spectra to obtain additional information regarding the energy dependence of LACs [9]. Implementing virtual monochromatic imaging (VMI) based on spectral CT can fundamentally eliminate BHA while providing rich energy attenuation information for component identification and characterization [8,10,11].

VMI based on spectral CT is indirectly achieved by acquiring material-specific images. The sum of the attenuations of all material-specific images at a specific energy constitutes a virtual monoenergetic image at that energy. According to the acquisition method of material-specific images, spectral CT can be divided into projection-based methods, image-based methods, and direct decomposition methods [12]. Projection-based methods first obtain material-specific projections from polychromatic projections and then reconstruct them using linear reconstruction methods. Alvarez et al. introduced the Compton effect and photoelectric effect basis functions and indirectly estimated the projections of two basic effects [13]. Abascal et al. developed a regularized iterative algorithm based on the Bregman distance to obtain basis material projections [14].

Image-based methods first perform traditional reconstruction of polychromatic projections to obtain CT images at different voltages. Then, the CT images are decomposed to obtain material-specific images [15,16]. Projection-domain or image-domain correction is usually performed to suppress artifacts in acquiring CT images. Zhang et al. developed a butterfly network based on the image-domain decomposition model for bi-material decomposition, which significantly reduces the noise of basis material images [17]. Li et al. combined penalized weighted least squares with regularization based on a mixed union of learned transform models to reduce the noise and artifacts of material images in dual-energy CT [18].

Direct decomposition methods attempt to solve material-specific images directly from polychromatic projections through statistical iteration, and simultaneously realizes material decomposition and image reconstruction [19]. Tilley et al. proposed a penalized weighted least squares material decomposition method based on the forward model and used a separable surrogate function to directly solve material images [9]. Zhang et al. introduced material decomposition into coded spectrum CT reconstruction, which can obtain higher-quality material images [20].

Complex system calibration or image correction processes [8,18,21,22] are usually required to suppress artifacts and achieve an accurate solution of material images, such as X-ray spectrum estimation, nonlinear model parameter estimation, and projection-domain hardening correction. These three types of spectral CT methods are limited in the practical implementation of VMI.

Recently, scholars have attempted to achieve narrow-energy-width imaging in blind scenarios in which the system parameters of the forward model are unknown, avoiding the energy spectrum calibration in the early stages of scanning. In cases in which the spectrum and material are unknown, Wei et al. proposed a narrow-energy-width imaging method based on the scattering low-frequency prior [23]. The method uses the local variance sum of the residuals as the objective function and adopts nonnegative matrix factorization (NMF) to solve the problem. Chen et al. proposed another narrow-energy-width imaging method based on basis effect decomposition to optimize the sum of the squared error [24]. Zhao et al. proposed a narrow-energy-width imaging algorithm using a Poisson prior. The method introduces quadratic potential function regularization to the projection-domain decomposition to suppress noise, and an efficient initialization method was proposed [25].

To further improve the accuracy of narrow-energy-width images, a statistical iterative spectral CT imaging method based on the blind separation of polychromatic projection is presented in this study. To realize direct inversion of material images in blind scenes and reduce the decomposition error, the system matrix is introduced into the X-ray multispectral forward model. The measured data are assumed to be statistically independent Poisson random variables. Based on the maximum likelihood principle, a constrained optimization problem for the spectral decomposition weight vector and volume fraction vector is established. An efficient alternating optimization algorithm is adopted to solve the problem under the block coordinate descent framework. Simulation and actual experiments indicate that the proposed algorithm improves the accuracy of acquiring narrow-energy-width images.

The remainder of this paper is organized as follows. Section 2 introduces the X-ray multispectral forward model with the system matrix, establishes the constrained optimization problem, and solves it. The simulation and actual experimental results are described in Section 3. Section 4 provides a discussion based on the experimental results. The final section concludes the paper.

2. Methods

2.1 X-ray multispectral forward model

Conventional X-ray CT systems produce continuous polychromatic X-rays. When a beam of polychromatic X-rays passes through an object, the average energy of the X-ray beam gradually increases with the transmitted thickness because low-energy photons are more easily absorbed than high-energy photons. This phenomenon is known as the beam hardening effect. According to Beer’s law, the ideal X-ray intensity at pixel $\mathbf {u}\in \Phi$ of the flat detector after beam hardening is

(1)$$\bar{I}_{n}(\mathbf{u}) = I_{n}^{(0)}(\mathbf{u})\int_{E_{min}}^{E_{max}}\mathrm{S}_{n}(E)exp\left ( -\int_{L(\mathbf{u})}\mu (\mathbf{x},E)\mathrm{d}l \right )\mathrm{d}E,$$

where subscript $n$ denotes the $n$th scan voltage, $I_{n}^{(0)}(\mathbf {u})$ denotes the initial X-ray intensity, $\mathrm {S}_n(E)$ denotes the X-ray normalized energy spectrum function which includes the X-ray source spectrum and detector energy dependence, $E_{min}$ and $E_{max}$ denote the minimum and maximum X-ray photon energies, $\mu (\mathbf {x},E)$ denotes the LAC at spatial position $\mathbf {x}\in \Omega$ and energy $E$, and $L(\mathbf {u})$ is the X-ray traversing path, which depends on the imaging geometry, such as the parallel beam, fan beam, and cone beam. The notations $\int _{E_{min}}^{E_{max}}$ and $\int _{L(\mathbf {u})}$ denote the integration over energy and ray paths, respectively. The notations $\Phi$ and $\Omega$ are the detector plane and object space, respectively. The diagram of the calculation system is shown in Fig. 1.

Fig. 1. Diagram of the calculation system. Cone-beam geometry is chosen here for illustration.

Download Full Size | PDF

Equation (1) shows that the detector collects attenuation information of photons at all transmitted energies of the continuous X-ray spectrum, which provides a basis for obtaining narrow-energy-width images from polychromatic measurement data [25]. We divide the continuous X-ray spectrum into $R$ narrow-energy-width spectra. The energy range of the $r$th narrow-energy-width spectrum is $[E_{r-1}, E_{r}]$, where $r=1,2,\ldots,R$, $E_0=E_{min}$, and $E_R=E_{max}$. Using the basis material decomposition model and generalized first integral mean value theorem, Eq. (1) is transformed into [13]:

(2)$$\bar{I}_{n}(\mathbf{u}) = I_{n}^{(0)}(\mathbf{u})\sum_{r=1}^{R}s_{n,r}exp\left ( -\sum_{k=1}^{K}\mu_{k}(\xi_r)\int_{L(\mathbf{u})}\mathrm{b}_{k}(\mathbf{x})\mathrm{d}l \right ),$$

where $\xi _r\in [E_{r-1}, E_{r}]$, $s_{n,r}=\int _{E_{r-1}}^{E_{r}}\mathrm {S}_{n}(E)\mathrm {d}E$ is the spectral decomposition weight satisfying $s_{n,r}\geqslant 0$ and $\sum _{r=1}^{R} s_{n,r}=1$ (called the sum-to-one constraint), $\mathrm {b}_{k}(\mathbf {x})$ is the spatial distribution of material $k$, $\int _{L(\mathbf {u})}\mathrm {b}_{k}(\mathbf {x})\mathrm {d}l$ denotes the thickness through the material $k$ along the ray path $L(\mathbf {u})$, $\mu _{k}(\xi _r)$ is the LAC of material $k$ at energy $\xi _{r}$, and $K$ is the number of basis materials. When the energy width of the narrow-energy-width spectrum is relatively small, such as less than 10 keV, $\xi _r=(E_{r-1}+E_r) / 2$ is desirable [25].

We discretize the object space $\Omega$ and assume that the object space consists of $J$ pixels located at $\mathbf {x}_j$. Likewise, we suppose that the detector plane $\Phi$ consists of an array of $P$ pixels located at $\mathbf {u}_p$. If the number of rotation angles is $Q$, then the total number of ray paths $M$ is $PQ$. The ideal X-ray intensity of the $m$th ray path is given as follows:

(3)$$\bar{I}_{n,m}= I_{n,m}^{(0)}\sum_{r=1}^{R}s_{n,r}exp\left ( -\sum_{k=1}^{K}\sum_{j=1}^{J}\mu_{k}(\xi_r)b_{k,j}a_{j,m} \right ),$$

where $a_{j,m}$ is the intersection length of the $m$th ray path and pixel $j$, and $b_{k,j}$ is the volume fraction of material $k$ at pixel $j$. The variable $b_{k,j}$ satisfy $b_{k,j}\geqslant 0$ and $\sum _{k=1}^{K} b_{k,j}=1$ when the basis materials are all object components. For simplicity, $\mathbf {q}_{r}=\left ( q_{r,1},\ldots,q_{r,J} \right )^{\mathrm {T}}$ is called the $r$th reconstructed image, where $q_{r,j}=\sum _{k=1}^{K}\mu _{k}(\xi _{r})b_{k,j}$, $\mathbf {p}_{r}=\left ( p_{r,1},\ldots,p_{r,M} \right )^{\mathrm {T}}$ is called the $r$th decomposed projection, where $p_{r,m}=\sum _{j=1}^{J}q_{r,j}a_{j,m}$. Because $\mathbf {p}_r$ can reflect the attenuation information of the imaging object in the $r$th narrow-energy-width spectrum, it is a narrow-energy-width projection. Correspondingly, $\mathbf {q}_r$ is a narrow-energy-width reconstructed image. The intersection lengths of all rays and pixels form the system matrix $\mathbf {A}=\left [ a_{j,m} \right ]_{J\times M}$. The introduction of $\mathbf {A}$ can avoid the indirect calculation of thickness and reduce the decomposition error.

In practice, the X-ray spectrum is widely obtained through indirect measurements [26]. However, the process is cumbersome, and even if the measurement result is obtained, it is difficult to maintain its accuracy and stability for a long time. When the imaging parameters such as voltage change and use time increase, the spectrum must be re-measured [27]. Therefore, we assume that the spectrum is unknown to avoid repeated measurements, i.e., the spectral decomposition weight $s_{n,r}$ is "blind". Let $\mathbf {b}=\left ( b_{1,1},\ldots,b_{K,1},\ldots,b_{1,J},\ldots,b_{K,J} \right )^{\mathrm {T}}$ be the volume fraction vector, and $\mathbf {s}=\left ( s_{1,1},\ldots,s_{N,1},\ldots,s_{1,R},\ldots,s_{N,R} \right )^{\mathrm {T}}$ be the spectral decomposition weight vector. Let $\mathbf {\bar {I}}=\left ( \bar {I}_{1,1},\ldots,\bar {I}_{N,1},\ldots,\bar {I}_{1,M},\ldots,\bar {I}_{N,M} \right )^{\mathrm {T}}$ be the ideal X-ray intensity vector. Thus, Eq. (3) can be expressed as

(4)$$\mathbf{\bar{I}}=\mathcal{F}\left ( \mathbf{s},\mathbf{b} \right ),$$

where $\mathcal {F}$ denotes the nonlinear mapping induced by the right-hand term in Eq. (3). Equation (4) is called the X-ray multispectral forward model. Here, "multispectral" refers to the spectra produced by multiple scanning voltages.

2.2 Maximum penalty likelihood estimation

The actual X-ray intensity $I_{n,m}$ acquired by the detector is assumed to obey a Poisson distribution, that is, $I_{n,m}\sim \mathrm {P}\left ( \lambda =\bar {I}_{n,m} \right )$, where $\lambda$ is the expectation of the Poisson distribution. Let $\mathbf {I}=\left ( I_{1,1},\ldots,I_{N,1},\ldots,I_{1,M},\ldots,I_{N,M} \right )^{\mathrm {T}}$ be the actual X-ray intensity vector. To estimate narrow-energy-width images from polychromatic measurements, the spectral decomposition weight vector $\mathbf {s}$ and the volume fraction vector $\mathbf {b}$ in Eq. (4) must be simultaneously solved. Thus, according to the maximum likelihood principle, the following negative log-likelihood function can be minimized:

(5)$$-\mathrm{L}(\mathbf{s}, \mathbf{b})=\sum_{n=1}^{N}\sum_{m=1}^{M}\mathcal{F}_{n,m}(\mathbf{s}, \mathbf{b})-I_{n,m}\ln\left(\mathcal{F}_{n,m}(\mathbf{s}, \mathbf{b})\right),$$

where $\mathcal {F}_{n,m}(\mathbf {s}, \mathbf {b})=\bar {I}_{n,m}$. Meanwhile, to improve the ill-posed nature of nonlinear decomposition problems, the following regularization function with edge preservation is obtained:

(6)$$\begin{aligned}[c] \mathrm{R}(\mathbf{b}) & =\sum_{k=1}^{K} \sum_{j=1}^{J} \sum_{j^{\prime} \in N_{j}} \beta_{k} \psi\left(b_{k, j}-b_{k, j^{\prime}}\right), \\ \psi\left(\Delta\right) & =\left\{\begin{matrix} \Delta^{2} / 2 & \left|\Delta\right|\leq \gamma \\ \gamma\left|\Delta\right|-\Delta^2/2 & \left|\Delta\right|> \gamma, \\ \end{matrix}\right. \end{aligned}$$

where $\psi$ denotes the Huber penalty functional, $\gamma$ is a tuning parameter, $N_j$ is the geometrical neighborhood of voxel $j$, and $\beta _k$ is the regularization parameter of material $k$ used to balance the regularization strengths of different materials.

Therefore, by referring to the maximum penalized likelihood estimation [28], we can establish the following constrained optimization problem:

(7)$$\begin{array}{ll} \min_{\mathbf{s}, \mathbf{b}} \quad & \theta(\mathbf{s}, \mathbf{b})={-}\mathrm{L}(\mathbf{s}, \mathbf{b})+\mathrm{R}(\mathbf{b}) \\ s.t. & \mathbf{s} \geq \mathbf{0}, \sum_{r=1}^{R} s_{n, r}=1, \forall n \\ & \mathbf{b} \geq \mathbf{0}, \sum_{k=1}^{K} b_{k, j}=1, \sum_{k=1}^{K} I_{\left\{b_{k, j} \neq 0\right\}} \leq 2, \forall j. \end{array}$$

where $\sum _{k=1}^{K} b_{k, j}=1$ is the volume conservation [29], and $I_{\left \{b_{k, j} \neq 0\right \}}$ is the indicator function. The constraint $I_{\left \{b_{k, j} \neq 0\right \}}\leq 2$ indicates that pixel $j$ comprises two types of materials. Similar to the literature [30,31], the constraint is used here to reduce the ill-posedness of the spectral CT multi-material decomposition. Notably, this constraint does not mean that the object contains only two components. Because the objective function $\theta \left (\mathbf {s}, \mathbf {b}\right )$ is a nonconvex function, by solving Eq. (7), we aim to solve the optimal local solution. The feasible region composed of the constraints in Eq. (7) can effectively reduce the solution space and improve the solution efficiency. In addition, to ensure that the problem is non-underdetermined, the scan voltage number $N$ is generally set equal to the material number $K$ (including air).

2.3 Optimization algorithm

Owing to the differences in the quantity scale, position, and constraints, the solution variables in Eq. (7) can be divided into two parts: $\mathbf {s}$ and $\mathbf {b}$. Inspired by the block coordinate descent approach [32], Eq. (7) is split into two subproblems: Eqs. (8) and (10). Eq. (7) can be optimized by alternately solving the two subproblems. The following are the two subproblems and their solution processes.

1) The volume fraction vector $\mathbf {b}$ is fixed, and the subproblem regarding the spectral decomposition weight vector $\mathbf {s}$ is

(8)$$\min_{\mathbf{s}} -\mathrm{L}\left(\mathbf{s},\mathbf{b}^{\left(i\right)}\right), s.t. \ \mathbf{s} \geq \mathbf{0}, \sum_{r=1}^{R} s_{n, r}=1, \forall n,$$

where $\mathbf {b}^{\left (i\right )}$ is a constant vector representing the volume fraction vector after the $i$th iteration. When the volume fraction vector is fixed, $\bar {I}_{n,m}$ is linear with respect to $\mathbf {s}$. Therefore, the objective function of Eq. (8) is a convex function, and the solution to Eq. (8) is a convex optimization. However, the problem is ill-conditioned, and the solution is susceptible to noise [33]. NMF [34] was adopted to derive the solution formula of Eq. (8). NMF is an effective method for blind source separation and is a first-order solver. It is less susceptible to noise and has non-negative constraints. Furthermore, the sum-to-one constraint can be easily introduced into the NMF using the method proposed by [35]. The final iteration formula is given as follows:

(9)$$s_{n,r}^{\left(i+1\right)}=s_{n,r}^{\left(i\right)}\frac{\sum_{m=1}^{M}I_{n,m}\exp\left({-}p_{r,m}^{\left(i\right)}\right)\bigg/\left(\sum_{{r}'=1}^{R}s_{n,{r}'}^{\left(i\right)}\exp\left({-}p_{r,m}^{\left(i\right)}\right)\right) + \delta I_{n,M}^{\left(0\right)}\bigg/\left(\sum_{{r}'=1}^{R}s_{n,{r}'}^{\left(i\right)}\right)}{\sum_{m=1}^{M}I_{n,m}^{\left(0\right)}\exp\left({-}p_{r,m}^{\left(i\right)}\right)+\delta I_{n,M}^{\left(0\right)}},$$

where $p_{r,m}^{(i)}=\sum _{j=1}^{J}\sum _{k=1}^{K}\mu _k\left (\xi _r\right )b_{k,j}^{(i)}a_{j,m}$ and $\delta$ controls the strength of the sum-to-one constraint. The larger the value of $\delta$, the closer $\sum _{r=1}^{R}s_{n,r}$ is to 1. We choose $\delta =\exp \left (-p_{r,M}^{(0)}\right )/10$ to ensure that the sum-to-one constraint is satisfied.

2) The spectral decomposition weight vector $\mathbf {s}$ is fixed, and the subproblem regarding the volume fraction vector $\mathbf {b}$ is

(10)$$\begin{aligned} \min_{\mathbf{b}} \ & \theta\left(\mathbf{s}^{(i)}, \mathbf{b}\right)={-}\mathrm{L}(\mathbf{s}^{(i)}, \mathbf{b})+\mathrm{R}(\mathbf{b}) \\ s.t. \ & \mathbf{b} \geq \mathbf{0}, \sum_{k=1}^{K} b_{k, j}=1, \sum_{k=1}^{K} I_{\left\{b_{k, j} \neq 0\right\}} \leq 2, \forall j. \end{aligned}$$

The objective function of Eq. (10) is nonconvex, solving which directly is difficult [30]. Therefore, the optimization transfer principle [36] is used to construct the pixel-wise separable quadratic surrogate function $\varphi \left (\mathbf {b};\mathbf {b}^{(i)}\right )$ (see Section 6 for details). The minimization of the original objective function $\theta \left (\mathbf {s}^{(i)}, \mathbf {b}\right )$ can be achieved by minimizing the surrogate function at each iteration step ($i$). Meanwhile, the minimization problem consisting of $\varphi \left (\mathbf {b};\mathbf {b}^{(i)}\right )$ and the constraints in Eq. (10) is easier to solve and is suitable for parallel computation.

To solve $\varphi \left (\mathbf {b};\mathbf {b}^{(i)}\right )$ under the constraints of Eq. (10), we have

(11)$$\mathbf{b}^{(i+1)}=\mathop{\arg\min}_{\mathbf{b}\geq \mathbf{0}} \left\{ \varphi\left(\mathbf{b};\mathbf{b}^{(i)}\right), s.t. \ \sum_{k=1}^{K} b_{k, j}=1, \sum_{k=1}^{K} I_{\left\{b_{k, j} \neq 0\right\}} \leq 2, \forall j\right\}.$$

Owing to the separable nature of the pixels, solving Eq. (11) is equivalent to solving the following problem for all $j$ one by one:

(12)$$\mathbf{b}_{j}^{(i+1)}=\mathop{\arg\min}_{\mathbf{b}_j\geq \mathbf{0}} \left\{ \varphi_{j}\left(\mathbf{b}_j;\mathbf{b}^{(i)}\right), s.t. \ \sum_{k=1}^{K} b_{k, j}=1, \sum_{k=1}^{K} I_{\left\{b_{k, j} \neq 0\right\}} \leq 2 \right\}.$$

The solution process in Eq. (12) is divided into the following two steps: (1) For $\forall k_1 \in \left \{1,\ldots,K-1\right \}$ and $\forall k_2 \in \left \{k_1+1,\ldots,K\right \}$, solve

(13)$$\mathbf{b}_j^{(k_1, k_2)} = \arg\min_{\mathbf{b}_j\geq\mathbf{0}}\left\{\varphi_{j}\left(\mathbf{b}_j;\mathbf{b}^{(i)}\right), s.t. \ \sum_{k=1}^{K} b_{k, j}=1, b_{k,j}=0, \forall k\notin \left\{k_1,k_2\right\} \right\}.$$

Equation (13) requires solving a quadratic programming problem of scale $K$, which can be efficiently implemented by the interior point method [37]. (2) We solve

(14)$$\left(k_1^{*}, k_2^{*}\right) = \mathop{\arg\min}_{(k_1,k_2)}\left\{\varphi_j\left(\mathbf{b}_j^{(k_1,k_2)};\mathbf{b}^{(i)}\right)\right\}$$

to obtain the final solution $\mathbf {b}_j^{(i+1)}\equiv \mathbf {b}_j^{\left (k_1^{*}, k_2^{*}\right )}$ of Eq. (12).

Under the block coordinate descent framework, the iterative formulas of Eq. (9) and Eq. (11) alternately iterate. To simplify the description, the solution process is called the NMF-SQS algorithm. The convergence conditions are given as follows: $i>N_{\mathrm {max}}$ or $1-\theta \left (\mathbf {s}^{(i+1)}, \mathbf {b}^{(i+1)}\right ) /\theta \left (\mathbf {s}^{(i)}, \mathbf {b}^{(i)}\right )<\varepsilon$, where $N_{\mathrm {max}}$ is the maximum number of iterations and $\varepsilon$ is a small convergence threshold. We empirically chose $N_{\mathrm {max}}=300$ and $\varepsilon =10^{-5}$. The objective function value and solving result change slightly when one of the two conditions is satisfied. We consider the algorithm to have converged at this point.

3. Result

In this section, numerical simulations and actual experiments were performed to evaluate the performance of the proposed algorithm. NMF-GN and joint statistical iterative reconstruction (JI-SIR) were used as comparison algorithms. NMF-GN is a projection-domain narrow-energy-width CT imaging method. It first decomposes polychromatic measurements to obtain narrow-energy-width projections; these projections then undergo linear reconstruction to provide narrow-energy-width images. The reconstruction algorithm generally adopts the filtered backprojection. JI-SIR is a direct material decomposition algorithm. The spectral parameters must be estimated in advance. Here, the method is used to generate a virtual monochromatic image of specific energy for comparison with the proposed algorithm. The spectral parameters directly use the known X-ray spectra in numerical simulations. Simultaneously, the initial state of the spectra (used to set the initial value of Eq. (9)) is substituted into JI-SIR for experimental comparison, which is called JI-SIR-S0.

Because the objective function in Eq. (7) is nonconvex, the initial value must be set carefully. The initial value of spectra was obtained by SpectrumGUI_1.03 software [38]. The material-related initial value was obtained by segmenting the directly reconstructed images of polychromatic projections. In particular, two materials with similar attenuation characteristics overlap in grayscale when BHA exists. At this time, they are treated as the "same" material for segmentation; after segmentation, the initial value of each material became half that of the segmentation result. For fair comparison, material regularization parameters in NMF-SQS and JI-SIR were set such that the final material images have similar noise levels. To compare the accuracy of narrow-energy-width images, each reconstructed image decomposed by NMF-GN and NMF-SQS was assumed to correspond to a certain energy in the experiments. We approximated each decomposed image obtained by NMF-GN and NMF-SQS as a virtual monoenergetic image at the median energy of the narrow-energy-width spectrum.

3.1 Simulation example

A fan-beam CT scan was performed on a numerical phantom, as shown in Fig. 2(a). The diameter of the phantom was 10.78 mm, and the diameter of each inner circle was 2.33 mm. The phantom was composed of four materials: air, polymethyl methacrylate (PMMA), magnesium (Mg), and aluminum (Al). The image size was $1024\times 1024$ (0.011 mm/pixel). For simplicity, a point X-ray source was simulated. The X-ray spectra generated by Spekpy v2.0 [39] are shown in Fig. 2(b). The anode material was tungsten. Each spectrum was filtered with 1.5 mm Al. Tube voltages of 60, 70, 80, and 90 kV correspond to the following number of photons: $8.29\times 10^{4}$, $1.17\times 10^{5}$, $1.57\times 10^{5}$, and $2.00\times 10^{5}$. Figure 2(c) depicts the LAC curve of each material obtained from the website of the National Institute of Standards and Technology (NIST) [40]. Suppose the number of detector elements of a flat panel detector is 832, and the pixel size is $0.127 mm\times 0.127 mm$. The distances from the source to the rotation center of the object (SOD) and detector center (SDD) were 140.00 mm and 780.57 mm, respectively. A total of 360 projections were sampled at equal intervals between $0^{\circ }$ and $360^{\circ }$. Poisson noise was added to each projection, and the influence of scattering was not considered.

Fig. 2. (a) Numerical phantom. (b) Spectra generated by the software Spekpy v2.0. (c) Attenuation coefficient curves.

Download Full Size | PDF

The directly reconstructed images at the four voltages are shown in Fig. 3. The image size was $512\times 512$ (0.022 mm/pixel). Severe banding artifacts exist in directly reconstructed images, and low-density information is missing. The outer grayscale of the Al region was larger than that of the central grayscale. Hence, cupping artifacts were observed in the region.

Fig. 3. Directly reconstructed images (a)–(d) at 60, 70, 80, and 90 kV. The display window is [-0.05, 0.3].

Download Full Size | PDF

The proposed algorithm was employed to decompose the simulated measurement data. Air was considered for decomposition, that is, $K = 4$. The energy range $\left [6, 90\right ]$ was divided equally into 17 intervals. Therefore, the number of narrow-energy-width spectra was $R = 17$, and the energy $E_{r}$ was $(r + 1) \times 5$. The parameter $\gamma = 0.01$ was set empirically, and the values of the regularization parameter $\beta _{k}$ were 80, 50, 200, and 200 for air, PMMA, Mg, and Al, respectively. The 5th, 9th, and 13th reconstructed images are shown in Fig. 4, corresponding to energies of 28, 48, and 68 keV. The first column is the ideal monoenergetic image at each corresponding energy. The remaining columns are obtained from JI-SIR, JI-SIR-S0, NMF-GN, and NMF-SQS. JI-SIR-S0 adopted the simulated spectra generated by SpectrumGUI_1.03, which is different from the spectra shown in Fig. 2(b). The peak signal-to-noise ratio (PSNR) [41] of each reconstructed image is shown in the upper left corner. JI-SIR-S0 still has severe banding artifacts and overly bright grayscale values in the Mg and Al regions. For the JI-SIR, NMF-GN, and NMF-SQS algorithms, the banding artifacts largely disappeared, and the grayscale in the Mg and Al regions was similar to that of the ideal monoenergetic image. In terms of the PSNR metric, NMF-SQS significantly outperformed the other algorithms. The PSNR of NMF-GN approaches or exceeds that of JI-SIR, whereas that of JI-SIR-S0 is below 20 dB with severe image degradation.

Fig. 4. Reconstructed images obtained by different algorithms. The 5th, 9th, and 13th reconstructed images are shown from the top to bottom, corresponding to the energies of 28, 48, and 68 keV. The first column is the ideal monoenergetic image. The remaining columns are obtained from JI-SIR, JI-SIR-S0, NMF-GN, and NMF-SQS. The display windows for rows 1 to 3 are [-0.05, 0.45], [0, 0.15], [0, 0.12].

Download Full Size | PDF

Figure 5 shows profiles of the 5th and 13th reconstructed images at the marked lines in Fig. 3(a). The graphs from top to bottom show the profiles of PMMA, Mg, and Al. Columns 1 and 2 correspond to the 5th and 13th reconstructed images, respectively. The profile errors of JI-SIR and NMF-GN are close in the Mg and Al regions, and the former is smaller in the PMMA region. The NMF-SQS profiles showed the lowest error.

Fig. 5. The profiles of the 5th and 13th reconstructed images at the marked lines in Fig. 3(a). The graphs from the top to bottom show the profiles of PMMA, Mg, and Al. Columns 1 and 2 correspond to the 5th and 13th reconstructed images, respectively.

Download Full Size | PDF

To compare the reduction degree of the hardening artifacts of NMF-GN and NMF-SQS, we selected two rectangular regions near and far from the center of the Al region ( Fig. 3(a)). The relative difference in the average LAC in the two regions was calculated. The smaller the relative difference, the stronger the consistency of the LAC in the region and the smaller the hardening artifacts. The relative differences of the four directly reconstructed images were 7.31%, 6.76%, 7.21%, and 6.72%, respectively. The results of NMF-GN and NMF-SQS are shown in Fig. 6(a). The relative difference for NMF-GN was below 3.80%, which is approximately one fold lower than that of directly reconstructed images. The relative difference for NMF-SQS was approximately 0.38%. In addition, although the beam hardening decreases with the increase of energy in actual imaging, due to the influence of algorithm error, the LAC difference of NMF-GN at higher energies was more than that at lower energies. This may be because projection domain regularization overly suppresses projection edges. The problem was improved in NMF-SQS.

Fig. 6. (a) The relative differences between the average LACs of each decomposed image in the two rectangular regions of Fig. 3(a). (b) NRMSE between each decomposed image and the corresponding ideal monoenergetic image.

Download Full Size | PDF

We calculated the normalized root mean square error (NRMSE) of the decomposed images. The NRMSEs of JI-SIR-S0 were all greater than 20%, and the NRMSEs of the other algorithms are shown in Fig. 6(b). Those of JI-SIR and NMF-GN are between 7% and 9%, and that of NMF-SQS is between 2% and 4%. The NRMSE of NMF-SQS was improved by 4% compared with that of JI-SIR and NMF-GN.

3.2 Real data examples

A full-angle scan of an actual phantom was performed using the FF20CT device with angular intervals of $1^{\circ }$. The phantom was a cylinder made of PMMA, Mg, and Al. Its diameter and height were 30 and 40 mm, respectively. The diameter of each inner cylinder was 6 mm. The focal spot size of the FF20CT device is 2 microns. Its target material is tungsten. The scan voltages were set to 70, 80, 90, and 100 kV, the current to 90, 80, 70, and 60 $\mathrm{\mu}$A, and the integration time to 769, 667, 588, and 556 ms. Since the size of the actual phantom is larger than the numerical phantom, a higher voltage setting was used here to improve the image quality. The SOD and SDD were 234.375 mm and 750 mm, respectively. The directly reconstructed images and LAC profiles at 80, 90, and 100 kV are shown in Fig. 7. Similar to the numerical simulations, directly reconstructed images have severe banding artifacts in the PMMA region, resulting in loss of low-density information. Cupping artifacts were observed in the Mg and Al regions.

Fig. 7. Directly reconstructed images (a)–(c) at 80, 90, and 100 kV. The display window is [0, 0.15]. (d)–(f) Attenuation coefficient profiles of the three reconstructed images at the marked line in (a).

Download Full Size | PDF

The proposed algorithm was employed to decompose the measurements at the four voltages. The energy range $\left [6, 100\right ]$ was divided equally into 19 intervals. The number of narrow-energy-width spectra was $R = 19$, and $E_{r} = (r + 1) \times 5$. We empirically set $\gamma = 0.01$, and the regularization parameters $\beta _{k}$ were 32, 40, 80, and 80 for air, PMMA, Mg, and Al, respectively. The 5th, 10th, 12th, and 15th reconstructed images are shown in Fig. 8, corresponding to energies of 28, 53, 63, and 78 keV. The reconstructed images from the top to bottom were obtained using the JI-SIR-S0, NMF-GN, and NMF-SQS algorithms. JI-SIR-S0 adopted the simulated spectra generated by SpectrumGUI_1.03. The spectra were also used to set the initial spectral decomposition weight $s_{n,r}^{(0)}$ of NMF-SQS. For JI-SIR-S0, the banding artifacts were reduced in the 10th reconstructed image but were still evident in the 5th, 12th, and 15th reconstructed images. The banding artifacts of both NMF-GN and NMF-SQS were significantly attenuated. Compared with NMF-GN, less noise exists in the decomposed images obtained by NMF-SQS, and the noise changes of the different decomposed images are relatively consistent.

Fig. 8. Reconstructed images decomposed by different algorithms. The 5th, 10th, 12th, and 15th reconstructed images are shown from the left to right, corresponding to energies of 28, 53, 63, and 78 keV. The reconstructed images from the top to bottom were obtained using the JI-SIR-S0, NMF-GN, and NMF-SQS algorithms. The display windows of columns 1 to 4 are [-0.05, 0.45], [0, 0.09], [0, 0.08], and [0, 0.06], respectively.

Download Full Size | PDF

Figure 9 shows the profiles of the 5th and 15th reconstructed images at the marked line in Fig. 7(a). The upper and lower horizontal dashed lines represent the theoretical LACs of Mg and PMMA at the corresponding energies. The JI-SIR-S0 profiles shift far from the reference values, especially at the Mg position. The NMF-GN and NMF-SQS profiles are close to the reference values. The NMF-GN profiles fluctuate more than those of NMF-SQS, especially in the middle position of PMMA. NMF-GN also has some errors in the Mg position of the 15th reconstructed image. The overall fluctuation of the NMF-SQS profiles is smaller. Notably, NMF-SQS has a slight error at the Mg edge.

Fig. 9. Profiles (a)–(b) of the 5th and 15th reconstructed images at the marked line in Fig. 7(a). The upper and lower horizontal dashed lines represent the theoretical LACs of Mg and PMMA at the corresponding energies, respectively.

Download Full Size | PDF

To evaluate the accuracy of the decomposed images, the regions of interest (ROI1, ROI2, and ROI3 in Fig. 7(a)) of PMMA, Mg, and Al were selected to calculate the mean, absolute percent error (APE) [42], root mean square relative error (RMSRE) [31], and standard deviation (STD). The calculated results of the 5th and 15th reconstructed images are listed in Table 1. In the PMMA region of the 5th reconstructed image, the APE and RMSRE of NMF-GN are 2.5 and 7.3 times higher than those of NMF-SQS, respectively. The difference between the APE and RMSRE of the other two materials was within 0.7% in the reconstructed image. In the 15th reconstructed image, the APE and RMSRE of NMF-SQS were both smaller than those of NMF-GN. Meanwhile, the STD of NMF-SQS was smaller than that of NMF-GN, indicating that the decomposed image of NMF-SQS is less noisy.

Table 1. Mean, APE, RMSRE, and STD for each material region in the 5th and 15th reconstructed images.

View Table

Figure 10 shows the APE and RMSRE curves for each material with respect to the decomposed image number. The NMF-GN error significantly varies with the decomposed image number. The APE fluctuation ranges of the three materials are 0.20%–11.5%, 0.05%–5.15%, 1.08%–2.90%, and the RMSRE fluctuation ranges are 7.02%–52.8%, 3.55%–6.33%, 2.37%–3.90%. The error distribution of NMF-SQS is relatively stable, and the APE and RMSRE are approximately 2.00% and 3.00%, respectively. Overall, the error of the NMF-SQS algorithm was smaller than those of the other algorithms.

Fig. 10. The APE curves (a)–(c) and RMSRE curves (d)–(f) for each material with respect to the decomposed image number.

Download Full Size | PDF

To verify the effectiveness of the proposed algorithm on natural objects, a full-angle scan of a 5 mm diameter granite cylinder was performed using the FF20CT device. The angle interval is $1^{\circ }$. The scan voltages were set to 60, 70, 80, and 90 kV. The SOD and SDD were 78 mm and 780 mm, respectively. A directly reconstructed image at 60 kV is shown in Fig. 11(a). The profiles marked by the blue and red lines are shown in Figs. 11(d) and 11(g), respectively. The outer edge of the granite is brighter in the directly reconstructed image. The LAC profiles exhibited cupping artifacts, and BHA was evident. The LACs of quartz and potassium feldspar overlapped, making the two materials indistinguishable.

Fig. 11. (a) The directly reconstructed image at 60 kV. The display window is [0, 0.77]. (b)–(c) 5th and 10th reconstructed images. The display windows are [0.00, 1.01] and [0, 0.20], respectively. (d)–(f) Profiles of reconstructed images (a)–(c) at the blue line. (g)–(i) Profiles of reconstructed images (a)–(c) at the red line.

Download Full Size | PDF

The measurement data were decomposed using the proposed algorithm. The photon energy range $\left [6, 90\right ]$ was equally divided into 17 intervals. The number of narrow-energy-width spectra was $R = 17$, and $E_{r} = (r + 1) \times 5$. The basis materials were set as air, quartz, potassium feldspar, and biotite. We empirically set $\gamma = 0.01$, and the regularization parameters $\beta _{k}$ were 40, 20, 50, and 100 for air, quartz, potassium feldspar, and biotite respectively. The 5th and 10th reconstructed images are shown in Figs. 11(b)–(c). Figures 11(e)–(f) show the profiles of the 5th and 10th reconstructed images at the blue marked line in Fig. 11(a). Figures 11(h)–(i) show the profiles of the 5th and 10th reconstructed images at the red marked line in Fig. 11(a). The grayscale decreased at the outer edges of the granite in the decomposed images. Cupping artifacts were significantly reduced, especially in the biotite region. Quartz and potassium feldspar were distinguished because of reduced hardening and increased material contrast in the decomposed images. Notably, the decomposed image was lighter in the middle region of the granite. Likewise, it was faintly visible in the directly reconstructed image but was not apparent. This phenomenon occurred because the middle position of the current slice is at the boundary of the two materials, which also verifies that the proposed algorithm improves the material contrast.

4. Discussion

Beam hardening causes banding, striping artifacts in the reconstructed images, and cupping artifacts for homogeneous objects. Narrow-energy-width imaging or VMI can effectively suppress these artifacts. Simultaneously, the narrow-energy-width or virtual monochromatic images at different energy intervals (energy) also have the specific LACs. In the above simulation and actual experiments, the directly reconstructed images have noticeable banding artifacts in the PMMA region, and the low-density information is largely missing. The outer values of the Mg and Al regions are higher than the central values, and cupping artifacts were observed (Figs. 3 and 7). JI-SIR-S0 led to bright or dark banding artifacts in the PMMA region (the 3rd column in Fig. 4 and the 1st row in Fig. 8), making the LACs poorly consistent. The NMRSEs in the numerical simulations were greater than 20%. The LACs in the Al and Mg regions also significantly differed from the reference values (Figs. 4 and 9), indicating that the virtual monochromatic images obtained by JI-SIR-S0 are inaccurate. This is because JI-SIR-S0 adopted the X-ray spectra that mismatch the forward model. Compared with JI-SIR-S0, owing to the adoption of accurate X-ray spectra, the images obtained by JI-SIR have accurate LACs, and the banding artifacts disappear (the 2nd column in Fig. 4 and Fig. 5).

When the spectra are unknown, the banding and cupping artifacts of the reconstructed images obtained by NMF-GN and NMF-SQS were also significantly reduced (columns 4–5 in Fig. 4 and rows 2–3 in Fig. 8). The LACs were distributed around the reference values (Figs. 5 and 9). The results show that the two algorithms can achieve effective narrow-energy-width imaging under unknown energy spectra. Different from JI-SIR-S0, although these two algorithms also adopted the mismatched spectral decomposition weight vector, the model error was reduced due to the continuous update of the spectral decomposition weight in the iteration. Compared with the directly reconstructed images, the grayscale consistency of NMF-GN and NMF-SQS in the Al region was improved, and the hardening artifacts were significantly reduced. At the same time, the grayscale consistency of NMF-SQS was higher than that of NMF-GN, and the hardening artifacts were lower (Fig. 6(a)). In the numerical simulations, the NRMSEs of NMF-SQS were improved by 4% compared with that of NMF-GN and JI-SIR (Fig. 6(b)). In the actual phantom experiment, the APE and RMSRE of NMF-GN considerably fluctuated (Fig. 9). The RMSREs of NMF-SQS are smaller than those of NMF-GN. In the partially reconstructed images, the APE of NMF-GN was also several times higher than that of NMF-SQS, especially in the PMMA region. Overall, the proposed algorithm can obtain more accurate narrow-energy-width images. And from the experimental results, the proposed algorithm does not lose the spatial resolution of the reconstructed image. In the granite experiment, the resulting reconstructed images show reduced cupping artifacts and enhanced material contrast; furthermore, quartz and potassium feldspar could be distinguished (Fig. 11), verifying that NMF-SQS can also achieve good results for natural objects.

Why does NMF-SQS obtain more accurate narrow-energy-width images than NMF-GN? First, NMF-SQS introduces the system matrix, which replaces the projection-domain material decomposition of NMF-GN with the direct inversion of material images. The direct decomposition method avoids the risk of missing information owing to the introduction of intermediate processes. Second, image-domain regularization with edge-preserving properties can better suppress noise and is more effective than projection-domain regularization. Finally, because of the ill-posed nature of the nonlinear inversion problem, constraining each pixel with two materials avoids crosstalk between materials with similar attenuation properties and has a specific regularization effect. Another problem introduced by direct inversion is the enormous computational cost, for which NMF-SQS uses the optimization transfer principle to solve the nonlinear constrained optimization problem. Each iteration only needs to solve multiple independent small-scale problems, facilitating parallelization and reducing computational cost and time.

However, notably, the LACs at the material edge were slightly lower in the actual phantom experiment (Fig. 9). This phenomenon was not observed in the numerical simulations. A possible reason for this is that factors such as scattering cannot be ignored in actual imaging. Because the scattering of conventional CT is unavoidable, for the next step, we can introduce scatter correction [43,44] to the proposed algorithm. In addition, the material regularization parameters $\beta _{k}$ affect each other, and to achieve the best denoising effect, the parameter combination must be set carefully. The method of setting the regularization parameters requires further study. In the above experiments, the materials PMMA, Mg, and Al were selected as the experimental objects to facilitate the algorithm verification. Subsequently, we will conduct more experimental studies based on different application scenarios (medicine or industry) to further refine the algorithm.

5. Conclusion

In this study, we proposed a statistical iterative spectral CT imaging method based on blind separation of polychromatic projections to improve the accuracy of narrow-energy-width image decomposition. To achieve direct inversion in blind scenarios, the system matrix is introduced into the X-ray multispectral forward model to directly solve material images, avoiding the indirect introduction of material thickness, which reduces the decomposition error. By treating the measurement data as statistically independent Poisson random variables, based on the maximum likelihood principle, a constrained optimization problem with edge-preserving regularization of the spectral decomposition weight vector and volume fraction vector is established. Based on the block coordinate descent framework, the problem is decomposed into two sub-problems which are solved alternately, and the optimization transformation principle is used to simplify the problem. Experiments show that, compared with existing algorithms, decomposed images of the proposed algorithm not only have fewer beam hardening artifacts, but also have more accurate LACs. The proposed algorithm improves the accuracy of obtaining narrow-energy-width images. Simultaneously, relevant results can also be applied to material characterization. The next steps will focus on scatter correction, regularization parameter adjustment, experimental extension, and multispectral image material characterization.

6. Derivation of the surrogate function

Here we use the optimization transfer principle [36] to derive the pixel-wise separable quadratic surrogate function $\varphi \left (\mathbf {b};\mathbf {b}^{(i)}\right )$ of $\theta \left (\mathbf {s}^{(i)}, \mathbf {b}\right )$. If $\varphi \left (\mathbf {b};\mathbf {b}^{(i)}\right )$ satisfies the condition [36]

(15)\begin{align}\varphi\left(\mathbf{b}^{(i)};\mathbf{b}^{(i)}\right) & =\theta\left(\mathbf{s}^{(i)},\mathbf{b}^{(i)}\right), \end{align}

(16)\begin{align}\frac{\partial \varphi\left(\mathbf{b};\mathbf{b}^{(i)}\right)}{\partial b_{k,j}}\bigg| _{\mathbf{b} =\mathbf{b}^{(i)}} & =\frac{\partial \theta\left(\mathbf{s}^{(i)},\mathbf{b}\right)}{\partial b_{k,j}}\bigg| _{\mathbf{b}=\mathbf{b}^{(i)}},\forall k,j , \end{align}

(17)\begin{align}\varphi\left(\mathbf{b};\mathbf{b}^{(i)}\right) & \geq\theta\left(\mathbf{s}^{(i)},\mathbf{b}\right), \end{align}

minimizing the surrogate function $\varphi \left (\mathbf {b};\mathbf {b}^{(i)}\right )$ of $\theta \left (\mathbf {s}^{(i)}, \mathbf {b}\right )$ at each iteration step ($i$) also minimizes $\theta \left (\mathbf {s}^{(i)}, \mathbf {b}\right )$.

To simplify the derivation of the surrogate function, we define that

(18)$$\begin{array}{c} h_{n,m}\left(\bar{I}_{n,m}\right)=\bar{I}_{n,m}-I_{n,m}\ln\left(\bar{I}_{n,m}\right),\ l_{k,r,m}\left(\mathbf{b}\right)=\sum_{j=1}^{J}\mu_{k}\left(\xi_r\right)b_{k,j}a_{j,m} , \\ \beta_{n,r,m}^{(i)}=\left[\sum_{{r}'=1}^{R}s_{n,{r}'}^{(i)}\exp\left(-\sum_{k=1}^{K}l_{k,r,m}\left(\mathbf{b}^{(i)}\right)\right)\right]\big/ \exp\left(l_{k,r,m}\left(\mathbf{b}^{(i)}\right)\right). \end{array}$$

Because $h_{n,m}\left (\bar {I}_{n,m}\right )$ is a convex function and $\sum _{r=1}^{R}s_{n,r}^{(i)} / \beta _{n,r,m}^{(i)}=1$, it follows from Jensen’s inequality that

(19)$$-\mathrm{L}\left(\mathbf{s}^{(i)}, \mathbf{b}\right) \leq \sum_{n=1}^{N}\sum_{m=1}^{R}\sum_{r=1}^{R}\frac{s_{n,r}^{(i)}}{\beta_{n,r,m}^{(i)}}h_{n,m}\left(I_{n,m}^{(0)}\beta_{n,r,m}^{(i)}\exp\left(-\sum_{k=1}^{K}l_{k,r,m}\left(\mathbf{b}\right)\right)\right) \equiv Q_{1}\left(\mathbf{b} ; \mathbf{b}^{(i)}\right).$$

Let

(20)$$g_{n,r,m}^{(i)}\left(\boldsymbol{l}_{r,m}\right) \equiv h_{n,m}\left(I_{n,m}^{(0)}\beta_{n,r,m}^{(i)}\exp\left(-\sum_{k=1}^{K}l_{k,r,m}\left(\mathbf{b}\right)\right)\right),$$

where $\boldsymbol {l}_{r,m}=\left (l_{1,r,m},\ldots,l_{K,r,m}\right )^{\mathrm {T}}$. The second-order expansion of $g_{n,r,m}^{(i)}\left (\boldsymbol {l}_{r,m}\right )$ is used to define the quadratic function:

(21)$$\begin{aligned} q_{n, r, m}^{(i)}\left(\boldsymbol{l}_{r, m}\right) = & g_{n, r, m}^{(i)}\left(\boldsymbol{l}_{r, m}^{(i)}\right)+\left.\sum_{k=1}^{K} \frac{\partial g_{n, r, m}^{(i)}}{\partial l_{k, r, m}}\right|_{\boldsymbol{l}_{r, m}=\boldsymbol{l}_{r, m}^{(i)}}\left(l_{k, r, m}-l_{k, r, m}^{(i)}\right) \\ & +\frac{1}{2} \sum_{k=1}^{K} \sum_{{k}'=1}^{K} T_{n, r, m}^{(i)}\left(l_{k, r, m}-l_{k, r, m}^{(i)}\right)\left(l_{k^{\prime}, r, m}-l_{k^{\prime}, r, m}^{(i)}\right), \end{aligned}$$

where $l_{k, r, m}^{(i)}=l_{k, r, m}\left (\mathbf {b}^{(i)}\right )$, $\boldsymbol {l}_{r, m}^{(i)}=\left (l_{1, r, m}^{(i)}, \ldots, l_{K, r, m}^{(i)}\right )^{\mathrm {T}}$, $T_{n,r,m}^{(i)}$ is the curvature. When choosing a suitable curvature to satisfy $g_{n, r, m}^{(i)}\left (\boldsymbol {l}_{r, m}\right ) \leq q_{n, r, m}^{(i)}\left (\boldsymbol {l}_{r, m}\right )$, we obtain

(22)$$Q_{1}\left(\mathbf{b} ; \mathbf{b}^{(i)}\right) \leq \sum_{n=1}^{N} \sum_{m=1}^{M} \sum_{r=1}^{R} \frac{s_{n, r}^{(i)}}{\beta_{n, r, m}^{(i)}} q_{n, r, m}^{(i)}\left(\boldsymbol{l}_{r, m}\right) \equiv Q_{2}\left(\mathbf{b} ; \mathbf{b}^{(i)}\right).$$

The curvature is approximated by selecting $I_{n,m}$[36]. Transform $l_{k,r,m}(\mathbf {b})$ as follows:

(23)$$l_{k, r, m}(\mathbf{b})=\sum_{j=1}^{J} \omega_{j, m}\left(\frac{a_{j, m} \mu_{k}\left(\xi_{r}\right)}{\omega_{j, m}}\left(b_{k, j}-b_{k, j}^{(i)}\right)+l_{k, r, m}^{(i)}\right) \equiv \sum_{j=1}^{J} \omega_{j, m} \lambda_{j, k, r, m}\left(\mathbf{b}_{j}\right),$$

where $\omega _{j, m}=a_{j, m} \big /\left (\sum _{j^{\prime }=1}^{J} a_{j^{\prime }, m}\right )$. Because Eq. (21) is a convex function, and $\sum _{j=1}^{J} \omega _{j, m}=1$, it follows from Jensen’s inequality that

(24)$$Q_{2}\left(\mathbf{b} ; \mathbf{b}^{(i)}\right) \leq \sum_{n=1}^{N} \sum_{m=1}^{M} \sum_{r=1}^{R} \frac{s_{n, r}^{(i)}}{\beta_{n, r, m}^{(i)}} \sum_{j=1}^{J} \omega_{j, m} q_{n, r, m}^{(i)}\left(\lambda_{j, k, r, m}\left(\mathbf{b}_{j}\right)\right)\equiv \sum_{j=1}^{J} \gamma_{n, r, j, m}^{(i)}\left(\boldsymbol{b}_{j}\right),$$

where $\mathbf {b}_{j}=\left (b_{1, j}, \ldots, b_{K, j}\right )^{\mathrm {T}}$.

For the regularization function in Eq. (10), because $\psi$ is a symmetric convex function, we get

(25)$$\begin{aligned} R(\mathbf{b}) & \leq \sum_{k=1}^{K} \sum_{j=1}^{J} \sum_{j^{\prime} \in N_{j}} \beta_{k}\left[\frac{1}{2} \psi\left(2 b_{k, j}-b_{k, j}^{(i)}-b_{k, j^{\prime}}^{(i)}\right)+\frac{1}{2} \psi\left({-}2 b_{k, j^{\prime}}+b_{k, j^{\prime}}^{(i)}+b_{k, j}^{(i)}\right)\right] \\ & =\sum_{k=1}^{K} \sum_{j=1}^{J} \sum_{j^{\prime} \in N_{j}} \beta_{k} \psi\left(2 b_{k, j}-b_{k, j}^{(i)}-b_{k, j^{\prime}}^{(i)}\right). \end{aligned}$$

Therefore, the pixel-wise separable quadratic surrogate function that satisfies conditions (15)–(17) is defined as follows:

(26)$$\begin{array}{cc} \varphi \left(\mathbf{b};\mathbf{b}^{(i)}\right) = \sum_{j=1}^J\varphi_{j}\left(\mathbf{b}_j;\mathbf{b}^{(i)}\right), \\ \varphi_{j}\left(\mathbf{b}_j;\mathbf{b}^{(i)}\right) = \gamma_{n,r,j,m}^{(i)}\left(\mathbf{b}_j\right) + \sum_{k=1}^{K}\sum_{{j}'\in N_j}\beta_k\psi \left(2b_{k,j}-b_{k,j}^{(i)}-b_{k,{j}'}^{(i)}\right). \end{array}$$

Funding

National Natural Science Foundation of China (61871351, 61971381, 62122070).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. H. Villarraga-Gómez, E. L. Herazo, and S. T. Smith, “X-ray computed tomography: from medical imaging to dimensional metrology,” Precision Engineering 60, 544–569 (2019). [CrossRef]

2. Y. Li, K. Li, C. Zhang, J. Montoya, and G. H. Chen, “Learning to reconstruct computed tomography images directly from sinogram data under a variety of data acquisition conditions,” IEEE Trans. on Med. Imaging 38, 2469–2481 (2019). [CrossRef]

3. Q. Zhao, X. Ma, A. Cuadros, T. Mao, and G. R. Arce, “Single-snapshot X-ray imaging for nonlinear compressive tomosynthesis,” Opt. Express 28, 29390–29407 (2020). [CrossRef]

4. C. Cai, T. Rodet, S. Legoupil, and A. Mohammad-Djafari, “A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography,” Med. Phys. 40, 111916 (2013). [CrossRef]

5. W. Zhao, D. Li, K. Niu, W. Qin, H. Peng, and T. Niu, “Robust beam hardening artifacts reduction for computed tomography using spectrum modeling,” IEEE Trans. on Comput. Imaging 5, 333–342 (2019). [CrossRef]

6. Q. Yang, W. K. Fullagar, G. R. Myers, S. J. Latham, T. Varslot, A. P. Sheppard, and A. M. Kingston, “X-ray attenuation models to account for beam hardening in computed tomography,” Appl. Opt. 59, 9126–9136 (2020). [CrossRef]

7. C. R. Romano, J. M. Minto, Z. K. Shipton, and R. J. Lunn, “Automated high accuracy, rapid beam hardening correction in X-Ray computed tomography of multi-mineral, heterogeneous core samples,” Comput. Geosci. 131, 144–157 (2019). [CrossRef]

8. L. Yu, S. Leng, and C. H. McCollough, “Dual-energy CT-based monochromatic imaging,” American Journal of Roentgenology 199, S9–S15 (2012). [CrossRef]

9. S. Tilley II, W. Zbijewski, and J. W. Stayman, “Model-based material decomposition with a penalized nonlinear least-squares CT reconstruction algorithm,” Phys. Med. Biol. 64, 035005 (2019). [CrossRef]

10. M. C. Jacobsen, S. L. Thrower, R. B. Ger, S. Leng, L. E. Court, K. K. Brock, E. P. Tamm, E. N. Cressman, D. D. Cody, and R. R. Layman, “Multi-energy computed tomography and material quantification: current barriers and opportunities for advancement,” Med. Phys. 47, 3752–3771 (2020). [CrossRef]

11. Z. Zhou, X. Zhang, R. Xin, L. Mao, J. Jia, S. Wei, T. Sheng, and J. Zheng, “Direct iterative basis image reconstruction based on MAP-EM algorithm for spectral CT,” Journal of Nondestructive Evaluation 40, 1–10 (2021). [CrossRef]

12. W. Fang, D. Wu, K. Kim, M. K. Kalra, R. Singh, L. Li, and Q. Li, “Iterative material decomposition for spectral CT using self-supervised Noise2Noise prior,” Phys. Med. Biol. 66, 15 (2021). [CrossRef]

13. R. E. Alvarez and A. Macovski, “Energy-selective reconstructions in X-ray computerised tomography,” Phys. Med. Biol. 21, 733 (1976). [CrossRef]

14. J. F. Abascal, N. Ducros, and F. Peyrin, “Nonlinear material decomposition using a regularized iterative scheme based on the Bregman distance,” Inverse Problems 34, 124003 (2018). [CrossRef]

15. W. Wu, P. Chen, V. V. Vardhanabhuti, W. Wu, and H. Yu, “Improved material decomposition with a two-step regularization for spectral CT,” IEEE Access 7, 158770–158781 (2019). [CrossRef]

16. J. Feng, H. Yu, S. Wang, and F. Liu, “Image-domain based material decomposition by multi-constraint optimization for spectral CT,” IEEE Access 8, 155450–155458 (2020). [CrossRef]

17. W. Zhang, H. Zhang, L. Wang, X. Wang, X. Hu, A. Cai, L. Li, T. Niu, and B. Yan, “Image domain dual material decomposition for dual-energy CT using butterfly network,” Med. Phys. 46, 2037–2051 (2019). [CrossRef]

18. Z. Li, S. Ravishankar, Y. Long, and J. A. Fessler, “DECT-MULTRA: Dual-energy CT image decomposition with learned mixed material models and efficient clustering,” IEEE Trans. on Med. Imaging 39, 1223–1234 (2020). [CrossRef]

19. K. Mechlem, S. Ehn, T. Sellerer, E. Braig, D. Münzel, F. Pfeiffer, and P. B. Noël, “Joint statistical iterative material image reconstruction for spectral computed tomography using a semi-empirical forward model,” IEEE Trans. on Med. Imaging 37, 68–80 (2017). [CrossRef]

20. T. Zhang, S. Zhao, X. Ma, A. P. Cuadros, Q. Zhao, and G. R. Arce, “Nonlinear reconstruction of coded spectral X-ray CT based on material decomposition,” Opt. Express 29, 19319–19339 (2021). [CrossRef]

21. M. M. Goodsitt, E. G. Christodoulou, and S. C. Larson, “Accuracies of the synthesized monochromatic CT numbers and effective atomic numbers obtained with a rapid kVp switching dual energy CT scanner,” Med. Phys. 38, 2222–2232 (2011). [CrossRef]

22. N. Ducros, J. F. P. J. Abascal, B. Sixou, S. Rit, and F. Peyrin, “Regularization of nonlinear decomposition of spectral X-ray projection images,” Med. Phys. 44, e174–e187 (2017). [CrossRef]

23. J. Wei, Y. Han, and P. Chen, “Narrow-energy-width CT based on multivoltage X-ray image decomposition,” International Journal of Biomedical Imaging 2017, 8126019 (2017). [CrossRef]

24. P. Chen, Y. Han, and Y. Li, “X-Ray multispectrum CT imaging by projection sequences blind separation based on basis-effect decomposition,” IEEE Trans. on Instrument. Meas. 70, 4502208 (2021). [CrossRef]

25. X. Zhao, P. Chen, J. Wei, and Z. Qu, “Spectral CT imaging method based on blind separation of polychromatic projections with Poisson prior,” Opt. Express 28, 12780–12794 (2020). [CrossRef]

26. B. Liu, H. Yang, H. Lv, L. Li, X. Gao, J. Zhu, and F. Jing, “A method of X-ray source spectrum estimation from transmission measurements based on compressed sensing,” Nucl. Eng. Technol. 52, 1495–1502 (2020). [CrossRef]

27. R. Gu and A. Dogandzic, “Blind x-ray ct image reconstruction from polychromatic poisson measurements,” IEEE Trans. on Comput. Imaging 2, 150–165 (2016). [CrossRef]

28. I. A. Elbakri and J. A. Fessler, “Statistical image reconstruction for polyenergetic X-ray computed tomography,” IEEE Trans. on Med. Imaging 21, 89–99 (2002). [CrossRef]

29. P. R. S. Mendonca, P. Lamb, and D. V. Sahani, “A flexible method for multi-material decomposition of dual-energy CT images,” IEEE Trans. on Med. Imaging 33, 99–116 (2014). [CrossRef]

30. Y. Long and J. A. Fessler, “Multi-material decomposition using statistical image reconstruction for spectral CT,” IEEE Trans. on Med. Imaging 33, 1614–1626 (2014). [CrossRef]

31. Y. Xue, W. Qin, C. Luo, P. Yang, Y. Jiang, T. Tsui, H. He, L. Wang, J. Qin, Y. Xie, and T. Niu, “Multi-material decomposition for single energy CT using material sparsity constraint,” IEEE Trans. on Med. Imaging 40, 1303–1318 (2021). [CrossRef]

32. Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” SIAM J. Imaging Sci. 6, 1758–1789 (2013). [CrossRef]

33. W. Ha, E. Y. Sidky, R. F. Barber, T. G. Schmidt, and X. Pan, “Estimating the spectrum in computed tomography via Kullback–Leibler divergence constrained optimization,” Med. Phys. 46, 81–92 (2019). [CrossRef]

34. D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in Proceedings of the 13th International Conference on Neural Information Processing Systems (2000), pp. 535–541.

35. D. C. Heinz and C. I. Chang, “Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE Trans. on Geosci. Remote Sens. 39, 529–545 (2001). [CrossRef]

36. H. Erdogan and J. A. Fessler, “Monotonic algorithms for transmission tomography,” IEEE Trans. on Med. Imaging 18, 801–814 (1999). [CrossRef]

37. A. Forsgren, P. E. Gill, and M. H. Wright, “Interior methods for nonlinear optimization,” SIAM Rev. 44, 525–597 (2002). [CrossRef]

38. “SpectrumGUI 1.03,” https://sourceforge.net/projects/spectrumgui/ (2007).

39. G. Poludniowski, A. Omar, R. Bujila, and P. Andreo, “Technical Note: SpekPy v2.0—a software toolkit for modeling X-ray tube spectra,” Med. Phys. 48, 3630–3637 (2021). [CrossRef]

40. M. Berger, J. Hubbell, S. Seltzer, J. Chang, J. Coursey, R. Sukumar, D. Zucker, and K. Olsen, “XCOM: photon cross sections database,” https://www.nist.gov/pml/xcom-photon-cross-sections-database (2010).

41. A. Horé and D. Ziou, “Image quality metrics: PSNR vs. SSIM,” in 2010 20th International Conference on Pattern Recognition (2010), pp. 2366–2369.

42. S. Kim and H. Kim, “A new metric of absolute percentage error for intermittent demand forecasts,” Int. J. Forecasting 32, 669–679 (2016). [CrossRef]

43. J. S. Lee and J. C. Chen, “A single scatter model for X-ray CT energy spectrum estimation and polychromatic reconstruction,” IEEE Trans. on Med. Imaging 34, 1403–1413 (2015). [CrossRef]

44. W. Zhao, D. Vernekohl, J. Zhu, L. Wang, and L. Xing, “A model-based scatter artifacts correction for cone beam CT,” Med. Phys. 43, 1736–1753 (2016). [CrossRef]

Material		The 5th reconstructed image			The 15th reconstructed image
		PMMA	Mg	Al	PMMA	Mg	Al
Reference: ${d m}^{- 1}$		3.92	19.2	36.4	2.10	3.45	5.56
Mean: ${d m}^{- 1}$	NMF-GN	4.08	19.6	36.8	2.05	3.62	5.71
	NMF-SQS	3.86	18.8	35.7	2.06	3.38	5.46
APE: %	NMF-GN	4.07	1.96	1.22	2.45	4.77	2.67
	NMF-SQS	1.65	2.05	1.89	1.90	2.11	1.87
RMSRE: %	NMF-GN	22.6	3.72	2.37	7.79	5.84	3.56
	NMF-SQS	3.08	3.13	2.38	2.79	3.11	2.37
STD: %	NMF-GN	0.872	0.608	0.741	0.155	0.117	0.131
	NMF-SQS	0.102	0.454	0.529	0.0429	0.0790	0.0812

Statistical iterative spectral CT imaging method based on blind separation of polychromatic projections

Abstract

1. Introduction

2. Methods

2.1 X-ray multispectral forward model

2.2 Maximum penalty likelihood estimation

2.3 Optimization algorithm

3. Result

3.1 Simulation example

3.2 Real data examples

4. Discussion

5. Conclusion

6. Derivation of the surrogate function

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (11)

Tables (1)

Equations (26)

Optics Express