Spectral CT imaging method based on blind separation of polychromatic projections with Poisson prior

Xiaojie Zhao; Ping Chen; Jiaotong Wei; Zhaoyan Qu

doi:10.1364/OE.392675

1. Introduction

X-ray computed tomography (CT) plays an irreplaceable role in medical diagnosis, material analysis and characterization, and nondestructive testing [1]. An incident X-ray beam in a conventional CT system is composed of photons with a wide energy range, and the linear attenuation coefficient (LAC) of the material is related to the photon energy [2]. Generally, the higher the energy is, the smaller the LAC. As X-rays traverse an object, the center of the spectrum shifts to higher energy, known as beam hardening. The beam hardening effect causes artifacts in reconstructed images from conventional methods, e.g., cupping and streaking artifacts, which limits the quantitative analysis of object components [3]. Spectral CT obtains additional information about the energy dependence of the LACs by using multiple different photon energy spectra. Considering the polychromatism of X-ray energy spectra, this technique has been applied to hardening artifact correction, quantitative imaging and contrast-to-noise ratio improvement [4].

Spectral CT reconstruction approaches are mainly divided into image-based, direct inversion and projection-based approaches [5,6]. Although relatively computationally simple, the image-based approach cannot eliminate hardening artifacts intrinsically [5]. The direct inversion approach can simultaneously realize measurement decomposition and image reconstruction. Long et al., [2] proposed a dual-energy CT (DECT) multimaterial decomposition algorithm capable of preserving edges. Barber et al., [6] developed a primal-dual algorithm to implement basis map decomposition. Chen et al., [7] developed a nonconvex optimization statistical reconstruction method that is also applicable to nonstandard scan configurations. Mechlem et al., [4] proposed a basis material statistical reconstruction algorithm based on a semiempirical forward model. Wu et al., [8] proposed a material simultaneous algebraic reconstruction technique and introduced a block matching framework to enhance the anti-noising ability. The direct inversion approach is mathematically elegant, and there is no risk of losing information without intermediate processes [9]. Some corresponding regularization methods are also mature, but the direct inversion approach is computationally expensive and slower to converge.

The projection-based approach separates decomposition and reconstruction processes. It first decomposes the line integrals of the base maps (material-specific line integrals) from measurements and then linearly reconstructs the line integrals [6,10]. The key to this approach is determining how to decompose measurements quickly and accurately. Schlomka et al., [11] used the Nelder-Mead (NM) method to solve the line integrals of basic effect and high-Z material coefficients. Subsequently, Schirra et al., [10] improved the quality of a reconstructed image by estimating the variance of each material-specific line integral using Fisher information. Ducros et al., [5] proposed a regularized weighted least squares Gauss-Newton (GN) algorithm to quickly solve the line integrals of basic materials. Recently, Abascal et al., [12] developed a regularized iterative method based on the Bregman distance to improve the global convergence of multimaterial decomposition. The projection-based approach can achieve parallelized decomposition and accelerate the calculation speed based on making full use of polychromatic data, but the decomposition error is easily introduced into the reconstruction process.

The above spectral CT reconstruction approaches decompose polychromatic attenuation measurements from the point of solving line integrals of basis maps, and most of these approaches need to know X-ray energy spectra in advance. Since actual X-ray energy spectra are difficult to measure [3] and the accuracy and stability of measured energy spectra are difficult to guarantee, these approaches are limited. In the case of unknown spectra, Chang et al., [13] proposed a DECT iterative reconstruction algorithm that simultaneously recovers the spectra and the material fraction images.

Note that each element response of the detector in a conventional CT system contains attenuation information of photons at all transmission energies of each continuous X-ray spectrum. Inspired by blind source separation methods [14], to achieve spectral separation of polychromatic measurements, Wei et al., [15,16] decomposed measurements to solve narrow-energy-width projections in the case of unknown X-ray energy spectra. They proposed weighted least squares and minimization for the sum of local variance of residual error with nonnegative constraints and used the first-order Karush-Kuhn-Tucker (KKT) condition to derive multiplicative iterative formulas for alternating updating; both algorithms improved the contrast ratio of materials with approximately linear attenuation coefficients to some extent. Some decomposed projections have the characteristics of narrow-energy-width projections.

However, although the two methods have certain effects for reducing hardening artifacts and improving material contrast, many decomposed projections are not narrow-energy-width projections and have higher noise than the original projections. Moreover, the energy interval corresponding to each decomposed projection is not clear, and an effective initialization approach is not given for the nonconvexity of decomposed problems. Recently, to determine the energy orientation of the narrow-energy-width projection, Zhao et al., [17] proposed a blind separation algorithm of multi-voltage projection sequences based on basic effect decomposition, but the algorithm still relies on good initial values.

Thus, a new X-ray spectral CT blind separation algorithm for solving narrow-energy-width projections is presented. It uses an X-ray multispectral forward model with the material prior. Moreover, measurements are considered statistically independent and follow a Poisson distribution. An efficient alternating optimization algorithm is employed, and an initialization method is presented in numerical and real experiments. The experimental results show that the new algorithm improves the accuracy of decomposed narrow-energy-width projections.

The rest of this paper is organized as follows: Section 2 derives the X-ray multispectral forward model and the proposed optimization algorithm. Section 3 describes the results of the simulation and real experiments. Section 4 discusses relevant issues, and the last section concludes the paper.

2. Method

2.1 X-ray multispectral forward model

In a conventional X-ray CT system, when an X-ray beam traverses an object, according to the Beer-Lambert theorem, the mean of the projection measurements at the pixel position ${\textbf u} \in {\mathbf \Phi }$ of the detector can be expressed as

(1)$$\bar{I}({\textbf u} )= {I^0}({\textbf u} )\int_E {S(E )\exp \left( { - \int_{L({\textbf u} )} {\mu ({{\textbf x},E} )dl} } \right)dE} ,$$

where L(u) denotes a ray path depending on the imaging geometry (parallel beam, fan beam, etc.), I⁰(u) denotes the incident X-ray intensity at the ray path L(u), μ(x,E) denotes the LAC for the object at energy E and spatial position ${\textbf x} \in {\mathbf \Omega }$ and S(E) denotes the normalized spectrum function. The notations ${\mathbf \Phi }$ and ${\mathbf \Omega }$ represent the detector plane and the object space, respectively.

Depending on the materials, μ(x,E) is decomposed into [18]

(2)$$\mu ({{\textbf x},E} )= \sum\limits_{k = 1}^K {{\mu _k}(E ){b_k}({\textbf x} )} ,$$

where μ_k(E) denotes the LAC function of material k independent of the object space and b_k(x) is the spatial distribution of material k. Each pixel position of the detector contains the attenuation information of photons at all transmitted energies of the continuous X-ray spectrum, which provides a basis for extracting narrow-energy-width projections from polychromatic measurements. The X-ray spectrum is divided into R narrow-energy-width spectra, each of which has a narrow energy range of [E_r_-1, E_r], where r = 1,…,R, E₀= 0, E_R= E_max and E_max denotes the maximum energy of the X-ray photons. Generally, each narrow energy range can be limited to 10 keV. According to the generalized first integral mean value theorem, Eq. (1) can therefore be expressed as

(3)$$\bar{I}({\textbf u} )= {I^0}({\textbf u} )\sum\limits_{r = 1}^R {{S_r}\exp \left( { - \sum\limits_{k = 1}^K {{\mu_k}({{\xi_r}} ){d_k}({\textbf u} )} } \right)} ,$$

where ${d_k}({\textbf u} )= \int_{L({\textbf u} )} {{b_k}({\textbf x} )dl} $, ${S_r} = \int_{{E_{r - 1}}}^{{E_r}} {{\mathop{\rm S}\nolimits} (E )dE}$ and ${\xi _r} \in [{{E_{r - 1}},{E_r}} ]$. Note that the energy range of each narrow-energy-width spectrum is small compared to the complete spectrum, within which the LACs of each material vary little. Thus, assume that the LACs of each material in the object are constant in each narrow-energy-width spectrum. Let ${\xi _r}\textrm{ = }{{({{E_{r - 1}}\textrm{ + }{E_r}} )} \mathord{\left/ {\vphantom {{({{E_{r - 1}}\textrm{ + }{E_r}} )} 2}} \right.} 2}$. Since S(E) is the normalized spectrum function, $\sum\nolimits_{r = 1}^R {{S_r}} = 1$ is obtained and called the sum-to-one constraint. In short, d_k(u) is the thickness, S_r is the spectral decomposition weight, ${p_r}({\mathbf \Phi } )= \sum\nolimits_{k = 1}^K {{\mu _k}({{\xi_r}} ){d_k}({\mathbf \Phi } )}$ is the r-th projected component, and p_r(Φ)'s reconstructed image is the r-th reconstructed component. Due to reflecting the attenuation information of the object on the narrow-energy-width spectrum, p_r(Φ) is the narrow-energy-width projection. An unknown X-ray spectrum is considered, that is, the spectral decomposition weight S_r is blind.

Assume that the number of spectra is N, and E_max is still the maximum energy of the X-ray photons under all spectra. The same narrow-energy-width spectrum division [E_r-₁, E_r] and energy value ξ_r are taken for each projection measurement of each spectrum. The index n is introduced to represent the variables Ī_n(u), $I_{n}^{0}(\mathbf{u})$ and S_n,r of different spectra. In addition, the spatial distribution of the detector pixels is discretized. Suppose that the detector is composed of a pixel array P located at u(p), p=1,…P, the number of angles during CT scanning is Y, and M = PY is the total number of rays.

Let $\bar{{\textbf I}} \in {{\textbf R}^{NM}}$ and ${\textbf I} \in {{\textbf R}^{NM}}$ be the mean and measured projection measurement vectors, which are defined by

(4)$$\bar{{\textbf I}} = {[{{{\bar{I}}_{1,1}}\; \cdots \;{{\bar{I}}_{N,1}}\; \cdots \; \cdots \;{{\bar{I}}_{1,M}}\; \cdots \;{{\bar{I}}_{N,M}}} ]^T},$$

(5)$${\textbf I} = {[{{I_{\textrm{1}, \textrm{1}}}\; \cdots \;{I_{N,1}}\; \cdots \; \cdots \;{I_{1,M}}\; \cdots \;{I_{N,M}}} ]^T}.$$

Similarly, ${\textbf s} \in {{\textbf R}^{NR}}$ denotes the spectral decomposition weight vector, and ${\textbf d} \in {{\textbf R}^{KM}}$ denotes the thickness vector, which are defined by

(6)$${\textbf s} = {[{{S_{\textrm{1,1}}} \cdots \;{S_{N,\textrm{1}}}\; \cdots \; \cdots \;{S_{\textrm{1},R}}\; \cdots \;{S_{N,R}}} ]^T},$$

(7)$${\textbf d} = {[{{d_{\textrm{1},\textrm{1}}} \cdots \;{d_{K,\textrm{1}}}\; \cdots \; \cdots \;{d_{\textrm{1},M}}\; \cdots \;{d_{K,M}}} ]^T}.$$

Thus, Eq. (3) can be expressed as

(8)$$\bar{{\textbf I}} = {\bf {\cal F}}({{\textbf s},{\textbf d}} ),$$

where ${\bf {\cal F}}$ denotes the nonlinear mapping induced by Eq. (3). Equation (8) is the X-ray multispectral forward model.

2.2 Regularized maximum likelihood estimation with constraints

The projection measurements are assumed to be subject to a Poisson distribution, i.e.,

(9)$${I_{n,m}} \sim P({\lambda = {{\bar{I}}_{n,m}}} ),$$

where m${\in} ${1,…,M}, and P($\lambda $) denotes a Poisson distribution with a mean of $\lambda $.

To decompose narrow-energy-width projections from polychromatic attenuation measurements, it is necessary to solve the spectral decomposition weight vector and thickness vector. To this end, according to the maximum likelihood principle, the constraint optimization problem is established by

(10)$$\begin{array}{l} \mathop {\min }\limits_{{\textbf s},{\textbf d}} \frac{1}{2}||{{\textbf I} - {\bf {\cal F}}({{\textbf s},{\textbf d}} )} ||_{\textbf W}^2 + \alpha {{\cal G}}({\textbf d} ),\\ s.\,t.\;{\textbf s} \ge {\textbf 0},{\textbf d} \ge {\textbf 0},\sum\nolimits_{r = 1}^R {{S_{n,r}}} = 1,\forall n, \end{array}$$

where $||{\textbf x} ||_{\textbf W}^2 = {{\textbf x}^T}{{\textbf W}^T}{\textbf Wx}$. The matrix ${\textbf W}$ is the weight matrix, which is chosen as

(11)$${\textbf W} = \textrm{diag}\left( {{1 \mathord{\left/ {\vphantom {1 {\sqrt {\textbf I} }}} \right.} {\sqrt {\textbf I} }}} \right)$$

so that ${{||{{\textbf I} - {\bf {\cal F}}({{\textbf s},{\textbf d}} )} ||_{\textbf W}^2} \mathord{\left/ {\vphantom {{||{{\textbf I} - {\bf {\cal F}}({{\textbf s},{\textbf d}} )} ||_{\textbf W}^2} 2}} \right.} 2}$ is the weighted least squares fidelity term applicable to the Poisson distribution. Due to the ill-posed nature of CT decomposition problems, a regularization function ${\cal G}$ is introduced, where α denotes the regularization parameter.

Given the differences in the priori information of materials in the object, the following regularization function is used:

(12)$${\cal G}({\textbf d} )= \sum\limits_{k = 1}^K {{\beta _k}{{\cal G}_k}({{{\textbf L}_k}{{\textbf d}_k}} )} ,$$

where ${\beta _k}$, ${{\cal G}_k}$, ${{\textbf d}_k} = {[{{d_{k,1}}\; \cdots \;{d_{k,M}}} ]^T}$ and ${{\textbf L}_k} \in {{\textbf R}^{{Q_k} \times M}}$ are the regularization parameter, regularization function, thickness vector and linear transformation of the k-th material, respectively. The linear transformation L_k is chosen as a discrete Laplace, where Q_k = 2M.

The internal compositions of most objects in nature are generally continuous distributions (components with uniform mixing of several materials can be regarded as a material), so the difference in the thickness d_k,m and d_k,m+₁ of each material on the adjacent ray path is small. Consequently, given the prior of smoothness, a quadratic potential function is used as the regularization function; that is, ${{\cal G}_k}({{{\textbf L}_k}{{\textbf d}_k}} )= ||{{{\textbf L}_k}{{\textbf d}_k}} ||_2^2$.

In addition, to ensure that Eq. (10) is not undetermined, N, R, K, and M must meet the following inequality:

(13)$$NM \ge NR + KM.$$

2.3 Optimization algorithm

In Eq. (10), although the feasible domain determined by the constraints is a convex set, the objective function is a nonconvex function. The optimization algorithm below aims to find an interesting (local) solution to the problem, and thus, an effective initialization approach is necessary.

Due to the differences in number scale and constraints, the solved variables are naturally divided into two blocks, s and d. Inspired by the block coordinate descent approach [19], s and d are alternately iterated using different algorithms to seek the solution of Eq. (10), where nonnegative matrix factorization (NMF) [20] and the GN [21] algorithm are used to update s and d, respectively.

NMF is a novel and effective blind separation method with inherent nonnegative constraints, and the sum-to-one constraint is introduced to NMF by employing a method in the literature [22]. The number scale of the thickness vector d is much larger than that of s, and the objective function for d is nonlinear in the framework of alternating optimization. Using a first-order solver together with s will slow the convergence speed of the algorithm as a whole. Given that the GN algorithm is a classic nonlinear optimization method with second-order convergence speed, it is adopted as the solver of d.

This paper minimizes Eq. (10) by alternately updating s and d using steps 1) and 2) below, where the i-th iteration proceeds as follows:

1) (NMF) The thickness vector is fixed as dⁱ^-1, and NMF is adopted to deal with the following subproblem of Eq. (10): $(14)$$\mathop {\min }\limits_{\textbf s} \frac{1}{2}||{{\textbf I} - {\bf {\cal F}}({{\textbf s},{{\textbf d}^{i - 1}}} )} ||_{\textbf W}^2,s.\,t.\;{\textbf s} \ge {\textbf 0},\sum\nolimits_{r = 1}^R {{S_{n,r}}} = 1,\forall n.$$$

Regardless of the sum-to-one constraint, the following formula can be obtained from the literature [20]:

(15)$$S_{n,r}^i = S_{n,r}^{i - 1}{{\left( {\sum\limits_{m = 1}^M {I_{n,m}^0t_{r,m}^{i - 1}} } \right)} \mathord{\left/ {\vphantom {{\left( {\sum\limits_{m = 1}^M {I_{n,m}^0t_{r,m}^{i - 1}} } \right)} {\left( {\sum\limits_{m = 1}^M {\frac{{{{({I_{n,m}^0} )}^2}}}{{{I_{n,m}}}}t_{r,m}^{i - 1}\left( {\sum\limits_{r = 1}^R {S_{n,r}^{i - 1}t_{r,m}^{i - 1}} } \right)} } \right)}}} \right.} {\left( {\sum\limits_{m = 1}^M {\frac{{{{({I_{n,m}^0} )}^2}}}{{{I_{n,m}}}}t_{r,m}^{i - 1}\left( {\sum\limits_{r = 1}^R {S_{n,r}^{i - 1}t_{r,m}^{i - 1}} } \right)} } \right)}},\forall n,r,$$

where $t_{r,m}^{i - 1} = \exp \left( { - \sum\nolimits_{k = 1}^K {{\mu_k}({{\xi_r}} )d_{k,m}^{i - 1}} } \right)$. In addition, to satisfy the sum-to-one constraint, I is extended to

(16)$${{\textbf I}_A} = {[{{{\textbf I}^T}\;{I_{1,M + 1}}\; \cdots \;{I_{N,M + 1}}} ]^T},\,{I_{n,M + 1}} = I_{n,M}^0\delta ,\,\forall n,$$

and let $t_{r,M + 1}^{i - 1} = \delta ,\,\forall r$, $I_{n,M + 1}^0 = I_{n,M}^0,\forall n$. The iterative formula of Eq. (14) is as follows by expanding M in Eq. (16) to M+1:

(17)$$S_{n,r}^i = \frac{{S_{n,r}^{i - 1}\left( {\sum\limits_{m = 1}^M {I_{n,m}^0t_{r,m}^{i - 1}} + I_{n,M}^0\delta } \right)}}{{\sum\limits_{m = 1}^M {\frac{{{{({I_{n,m}^0} )}^2}}}{{{I_{n,m}}}}t_{r,m}^{i - 1}\left( {\sum\limits_{r = 1}^R {S_{n,r}^{i - 1}t_{r,m}^{i - 1}} } \right) + I_{n,M}^0\delta \sum\limits_{r = 1}^R {S_{n,r}^{i - 1}} } }},\forall n,r,$$

where $\delta $ controls the impact of the sum-to-one constraint. The larger that $\delta $ is, the closer that $\sum\nolimits_{r = 1}^R {{S_{n,r}}}$ is to 1.

2) (GN) The spectral decomposition weight vector is fixed as sⁱ, and the GN algorithm is employed to solve the following subproblem of Eq. (10):

(18)$$\mathop {\min }\limits_{\textbf d} \frac{1}{2}||{{\textbf I} - {\bf {\cal F}}({{{\textbf s}^i},{\textbf d}} )} ||_{\textbf W}^2 + \alpha {{\cal G}}({\textbf d} ),s.\,t.\,{\textbf d} \ge {\textbf 0}.$$

For nonnegative constraints, an interior point (IP) method [23] is used. The logarithmic barrier function is written as

(19)$$P({{\textbf d},u} )= \frac{1}{2}||{{\textbf I} - {\bf {\cal F}}({{{\textbf s}^i},{\textbf d}} )} ||_{\textbf W}^2 + \alpha {{\cal G}}({\textbf d} )- u\sum\limits_{m = 1}^M {\sum\limits_{k = 1}^K {\ln ({{d_{k,m}}} )} } ,$$

where u denotes the barrier parameter. For a given barrier parameter ${u^j}$, the iteration formula of Eq. (19) is as follows:

(20)$${{\textbf d}^i} = {{\textbf d}^{i - 1}} + {\gamma ^i}{{\mathbf \delta }^i},$$

where ${{\mathbf \delta }^i} \in {{\textbf R}^{KM}}$ denotes the GN direction and ${\gamma ^i} \in {\mathop{\rm R}\nolimits} $ denotes the search step.

The GN direction is obtained by solving the following equation:

(21)$$({{{({{{\textbf J}^{i - 1}}} )}^T}{{\textbf W}^T}{\textbf W}{{\textbf J}^{i - 1}} + \alpha {\nabla^2}{\cal G}({{{\textbf d}^{i - 1}}} )+ {u^j}{{\textbf D}^{i - 1}}{{\textbf D}^{i - 1}}} ){{\mathbf \delta }^i} ={-} \nabla P({{{\textbf d}^{i - 1}},{u^j}} ),\,$$

where ${{\textbf D}^{i - 1}} = {\mathop{\rm diag}\nolimits} ({{\textrm{1} \mathord{\left/ {\vphantom {\textrm{1} {{{\textbf d}^{i - 1}}}}} \right.} {{{\textbf d}^{i - 1}}}}} )$, ${{\textbf J}^{j - 1}}$ denotes the Jacobi matrix of ${\bf {\cal F}}({{{\textbf s}^i},{\textbf d}} )$ at ${{\textbf d}^{i - 1}}$, ${\nabla ^2}{\cal G}({{{\textbf d}^{i - 1}}} )$ denotes the Hessian matrix of ${\cal G}$ at ${{\textbf d}^{i - 1}}$ and $\nabla P({{{\textbf d}^{i - 1}},{u^j}} )$ denotes the gradient of $P({{\textbf d},{u^j}} )$ at ${{\textbf d}^{i - 1}}$. $\nabla P({{{\textbf d}^{i - 1}},{u^j}} )$ is given by

(22)$$\nabla P({{{\textbf d}^{i - 1}},{u^j}} )={-} {({{{\textbf J}^{i - 1}}} )^T}{{\textbf W}^T}{\textbf W}({{\textbf I} - {\bf {\cal F}}({{{\textbf s}^i},{{\textbf d}^{i - 1}}} )} )+ \alpha \nabla {\cal G}({{{\textbf d}^{i - 1}}} )- {u^j}\frac{1}{{{{\textbf d}^{i - 1}}}},$$

where $\nabla {\cal G}({{{\textbf d}^{i - 1}}} )$ denotes the gradient of ${\cal G}$ at ${{\textbf d}^{i - 1}}$.

The Jacobi matrix ${\textbf J} \in {{\textbf R}^{NM \times KM}}$ introduced in Eq. (21) is defined as

(23)$${\textbf J} = \frac{{\partial {\bf {\cal F}}({{{\textbf s}^i},{\textbf d}} )}}{{\partial {{\textbf d}^T}}}.$$

Let ${\bar{{\textbf I}}^i} = { \mathbf{\cal F}}({{{\textbf s}^i},{\textbf d}} )= {[{\bar{I}_{1,1}^i\; \cdots \;\bar{I}_{N,1}^i\; \cdots \; \cdots \;\bar{I}_{1,M}^i\; \cdots \;\bar{I}_{N,M}^i} ]^T}$; then,

(24)$$\frac{{\partial \bar{I}_{n,m}^i}}{{\partial {d_{k,m^{\prime}}}}} = \left\{ {\begin{array}{{cc}} {0,}&{m \ne m^{\prime}}\\ {l_{n,k,m}^i,}&{m = m^{\prime}} \end{array}} \right.$$

with

(25)$$l_{n,k,m}^i ={-} I_{n,m}^0\sum\limits_{r = 1}^R {S_{n,r}^i{\mu _k}({{\xi_r}} )\exp \left( { - \sum\limits_{k = 1}^K {{\mu_k}({{\xi_r}} ){d_{k,m}}} } \right)} .$$

Thus, ${\textbf J}$ can be represented as a block diagonal matrix

(26)$${\textbf J} = {\mathop{\rm diag}\nolimits} ({[{{{\textbf J}_1}\; \cdots \;{{\textbf J}_M}} ]} ),$$

where J_m ${\in} {{\textbf R}^{N \times K}}$ and (J_m)_n,k = $l_{n,k,m}^i$. Obviously, J is a sparse matrix with only NKM nonzero elements.

The gradient of the regularization function ${\cal G}$ is calculated as follows:

(27)$$\nabla {\cal G}({\textbf d} )= \sum\limits_{k = 1}^K {{\beta _k}{\nabla _k}{{\cal G}_k}({{{\textbf L}_k}{{\textbf d}_k}} )\otimes {{\textbf e}_k}} ,$$

where ${\nabla _k}$ denotes the gradient with respect to ${{\textbf d}_k}$ and is defined as

(28)$${\nabla _k} = {\left[ {\frac{\partial }{{\partial {d_{k,1}}}}\; \cdots \;\frac{\partial }{{\partial {d_{k,M}}}}} \right]^T}, $$

${\nabla _k}{{\cal G}_k}({{{\textbf L}_k}{{\textbf d}_k}} )= 2{\textbf L}_k^T{{\textbf L}_k}{{\textbf d}_k}$, ${{\textbf e}_k}$ denotes a K-dimensional natural basis vector, and ${\otimes} $ denotes the Kronecker product.

Similarly, the Hessian matrix of the regularization function ${\cal G}$ is given as

(29)$${\nabla ^2}{\cal G}({\textbf d} )= \sum\limits_{k = 1}^K {{\beta _k}\nabla _k^2{{\cal G}_k}({{{\textbf L}_k}{{\textbf d}_k}} )\otimes {{\textbf e}_k} \otimes {\textbf e}_k^T} ,$$

where $\nabla _k^2{{\cal G}_k}({{{\textbf L}_k}{{\textbf d}_k}} )= 2{\textbf L}_k^T{{\textbf L}_k}$.

The selection of the search step ${\gamma ^i}$ will determine whether dⁱ in Eq. (20) can be sustained as a feasible solution. Backtracking linear search [24] is used to determine the step size as follows:

(1) Let ${\bar{\gamma }^i} = \max \left\{ {{1 \mathord{\left/ {\vphantom {1 {\left( {1 + \sqrt { - \nabla P{{({{{\textbf d}^{i - 1}},{u^j}} )}^T}{{\mathbf \delta }^i}} } \right)}}} \right.} {\left( {1 + \sqrt { - \nabla P{{({{{\textbf d}^{i - 1}},{u^j}} )}^T}{{\mathbf \delta }^i}} } \right)}},{\gamma_0}} \right\}$, where ${\gamma _0}$ denotes a constant.
(2) Select ${\gamma ^i} = {\tau ^c}{\bar{\gamma }^i}$, where $c \ge 0$ is the smallest integer that enables ${\gamma ^i}$ to meet conditions Eq. (30) and Eq. (31), $\tau \in ({0,1} )$ denotes the step size adjustment parameter and $\tau ^{\prime}$ denotes the backtracking parameter. $(30)$${{\textbf d}^{i - 1}} + {\gamma ^i}{{\mathbf \delta }^i} \ge {\textbf 0}.$$$ $(31)$$P({{{\textbf d}^{i - 1}} + {\gamma^i}{{\mathbf \delta }^i},{u^j}} )\le P({{{\textbf d}^{i - 1}},{u^j}} )+ {\gamma ^i}\tau ^{\prime}{({\nabla P({{{\textbf d}^{i - 1}},{u^j}} )} )^T}{{\mathbf \delta }^i}.$$$

Iterative formulas of Eq. (17) and Eq. (20) iterate alternately. When ${||{{{\textbf d}^i} - {{\textbf d}^{i - 1}}} ||_2} < {\varepsilon _1}{||{{{\textbf d}^i}} ||_2}$ and $M \times {u^j} \ge {\varepsilon _2}$ meet, the barrier parameter is updated to ${u^{j + 1}} = \theta {u^j}$; then, proceed to the (i+1)-th iteration, where ${\varepsilon _1} > 0$ denotes the GN stop threshold, ${\varepsilon _2} > 0$ denotes the IP stop threshold and $\theta $ denotes the IP adjustment factor. The iteration continues until $M \times {u^j} < {\varepsilon _2}$ or $i > T$ is satisfied, where T denotes the maximum number of iterations.

The proposed algorithm is referred to as the NMF-GN algorithm for simplicity.

3. Results

To evaluate the proposed NMF-GN algorithm, simulation and real experiments are presented. The reconstruction algorithm adopts the filtered backprojection (FBP). Set ${\gamma _0} = \textrm{1}$, $\tau = 0.\textrm{8}$, $\tau ^{\prime} = 0.01$, ${\varepsilon _1} = {10^{ - 4}}$, ${u^0} = 0.01$, $\theta = 0.1$, ${\varepsilon _2} = 1$ and $T = 500$ in the GN step. The values of these parameters are relatively broad because they are not sensitive to the solved results. The regularization parameter $\alpha $ is selected by the L-curve method [25].

An effective initialization approach is provided for the nonconvexity of Eq. (10). The initial value d⁰ is obtained by projecting forward the digitized model based on the imaging parameters, in which a slight disturbance is added to each projection value. The digital model is determined by segmenting an FBP reconstructed image of polychromatic projections by thresholds under a certain spectrum. In particular, when the gray values of two materials in the image overlap (such as silicon and aluminum), these two materials are divided into “the same material”, and the thickness vector of both materials is one-half of the calculated result. The initial value s⁰ is calculated in advance by Eq. (17). In the calculation of s⁰, the thickness vector is fixed to d⁰ and the initial value of s satisfies

(32)$$\begin{array}{l} \quad \quad \quad \quad \quad \sum\limits_{r = 1}^R {{S_{n,r}}} = 1,\\ {S_{n,r}} = 0,\forall ({n,r} )\in \{{({n,r} )|{{E_{n,\max }} \le {E_{r - 1}}} } \}, \end{array}$$

where E_n,_max is the maximum energy of the n-th energy spectrum. The initial value is determined by using the spectra of SpectrumGUI_1.03 software in the experiments. In particular, the spectra are filtered by 2 mm aluminum in the simulation experiment. Stop the iteration and output s⁰ when the maximum number of iterations, 1000, is reached or ${||{{{\textbf s}^t} - {{\textbf s}^{t - 1}}} ||_2} < {||{{{\textbf s}^t}} ||_2} \times {10^{ - 5}}$ is satisfied. Since Eq. (17) is multiplicative, S_n,r set to 0 will always be constant. The parameter $\delta = \sum\nolimits_{m = 1}^M {t_{1,m}^0} $ is employed in this paper.

3.1 Simulation example

The cylindrical phantom described in Fig. 1(a) consists of eight circles with the same diameter, each representing silicon or aluminum, with an image size of 256 × 256. A fan-beam CT scan of the phantom was simulated. The tube voltages were chosen as 70, 80 and 90 kV, and the tube currents were set to 0.7 mA. Spectrum curves generated by SpectrumGUI_1.03 are shown in Fig. 1(b). The flat panel detector width and pixel size were set to 832 and 0.127 mm × 0.127 mm, respectively. The distance from the source to the detector (SDD) is 780.57 mm, and the distance from the source to the object (SOD) is 140 mm. Each voltage was evenly spaced from 0 to 360 by 360 angles and Poisson noise was added to each polychromatic measurement. The influence of scattering and other factors was ignored in the simulation.

Fig. 1. The cylindrical phantom (a) and spectra (b) obtained from SpectrumGUI_1.03.

Download Full Size | PDF

The directly reconstructed image at 70 kV is shown in Fig. 2(a). The energy range [0, 90] was equally divided into 18 intervals; that is, the number of narrow-energy-width spectra is R = 18, and let E_r = r × 5. The proposed algorithm was employed to decompose polychromatic simulation measurements. Since the distributions of the two materials in the simulated phantom are not very different, ${\beta _k} = 1$ was set as in the experiments below. The regularization parameter was selected as $\alpha = 1$. The 6th, 12th, and 18th projection components were selected from the decomposed results, and their reconstructed images are shown in Figs. 2(b)–2(d). Figures 2(e)–2(h) show the grayscale curves of each reconstructed image in Figs. 2(a)–2(d) at the blue line position shown in Fig. 2(a). The green and red lines on these curves represent the grayscale mean levels of silicon and aluminum, respectively. In the simulation experiment, due to the large noise, each reconstructed image was denoised by using 3 × 3 median filter.

Fig. 2. The directly reconstructed image (a) at 70 kV and the 6th, 12th and 18th reconstructed components (b)-(d). The grayscale curves (e)-(h) of the directly reconstructed image and the 6th, 12th, and 18th reconstructed components at the blue line position shown in (a). The green and red lines on these curves represent the grayscale mean levels of silicon and aluminum, respectively.

Download Full Size | PDF

Obvious cupping artifacts exist in the directly reconstructed image, resulting in overlap of the grayscale levels of aluminum and silicon (Figs. 2(a) and 2(e)). In contrast, the cupping artifacts nearly disappear in the reconstructed components, and the grayscale levels of silicon and aluminum in the 6th and 18th reconstructed components are different (Figs. 2(b)–2(d), 2(f)–2(h)).

To evaluate the relationship between the reconstructed components and the narrow-energy-width spectra, the average LAC ${\bar{g}_{r,k}}$ (third and fifth columns in Table 1) of each material in each reconstructed component was calculated. We searched for the closest energy ${E_{r,k}}$ (fourth and sixth columns in Table 1) to ${\bar{g}_{r,k}}$ from the attenuation coefficient mapping table. The attenuation coefficient mapping table was obtained by cubic spline interpolation of attenuation data from the National Institute of Standards and Technology (NIST). In Table 1, k = 1 represents silicon, and k = 2 represents aluminum. Except for the last two reconstructed components with slightly larger energy differences (1.08 and 1.33 keV), the energy difference between E_r_,1 and E_r_,2 is within the range of -0.65 to 0.65 keV (Table 1). In addition, there are ${E_{r,k}} \in [{{E_{r - 1}},{E_r}} ]$ for any r and k in Table 1.

Table 1. The average attenuation coefficient of each material and its corresponding energy for partially reconstructed components in the simulation experiment.

View Table | View all tables in this article

3.2 Real data examples

Three samples were scanned using a YXLON FF20CT device and labeled as samples 1 to 3. The flat panel detector used by the device has 1122 × 1122 pixels, and each pixel is 0.127 mm × 0.127 mm in size. The three samples were scanned at 360, 360, and 1080 angles from 0 to 360 at equal intervals. For the sake of simplicity, only the fan-beam data at the center of the source were processed.

Sample 1 is composed of 4 irregular small pieces (2 silicon, 2 aluminum). The scan voltages are 40, 50 and 60 kV, and the current is 100 uA. To reduce the noise at low voltage, the exposure time was chosen as 2 s, 1 s and 0.5 s, respectively. The SDD is 600 mm, and the SOD is 80 mm. The directly reconstructed image at 50 kV is shown in Fig. 3(a). Polychromatic attenuation measurements obtained by CT scanning were decomposed by the proposed NMF-GN algorithm. The energy range [0, 60] was equally divided into 12 intervals; that is, the number of narrow-energy-width spectra is R = 12, and let E_r = r × 5. The regularization parameters were chosen as α = 10 and β_k = 1 during the decomposition process.

Fig. 3. The directly reconstructed image (a) at 50 kV and the 6th and 12th reconstructed components (b)-(c) for sample 1 decomposed by the NMF-GN algorithm. (d)-(f) show the grayscale curves of the reconstructed images in (a)-(c) at the blue line position in (a). The green and red lines on these curves represent the grayscale mean levels of silicon and aluminum, respectively.

Download Full Size | PDF

The 6th and 12th projection components were selected from the decomposition results, and their reconstructed images are shown in Figs. 3(b)–3(c). Figures 3(d)–3(f) show the grayscale curves of the reconstructed images in Figs. 3(a)–3(c) at the blue line position shown in Fig. 3(a). The green and red lines on these curves represent the grayscale mean levels of silicon and aluminum, respectively. Compared with direct reconstruction in Fig. 3(a), the reconstructed images in Figs. 3(b)–3(c) after decomposition not only reduce the cupping artifacts but also show that the LAC of silicon is higher than that of aluminum in one component and close to that of aluminum in another component (Figs. 3(e)–3(f)). Thus, the proposed algorithm avoids overlap of the gray values of the two materials and enables the two materials to be distinguished.

Sample 2 is also composed of four irregular small pieces (2 magnesium, 2 aluminum), and the attenuation characteristics of the two materials constituting the sample differ greatly. The scan voltages are 60, 80 and 100 kV, and the currents and exposure times are 100, 80, and 60 uA and 833, 500, and 333 ms, respectively. The SDD is 750 mm, and the SOD is 120 mm. The directly reconstructed image at 60 kV is shown in Fig. 4(a). The energy range [0, 100] was equally divided into 10 intervals; that is, the number of narrow-energy-width spectra is R = 20, and let E_r = r × 5. The regularization parameters were set to α = 10 and β_k = 1 during the decomposition process. The 8th and 16h projection components were selected from the decomposition results, and their reconstructed images are shown in Figs. 4(b)–4(c). Figures 4(d)–4(f) show the grayscale curves of the reconstructed images in Figs. 4(a)–4(c) at the blue line position shown in Fig. 4(a). The green and red lines on these curves represent the grayscale mean levels of aluminum and magnesium, respectively.

Fig. 4. The directly reconstructed image (a) at 60 kV and the 8th and 16th reconstructed components (b)-(c) for sample 2 decomposed by the NMF-GN algorithm. (d)-(f) show the grayscale curves of the reconstructed images in (a)-(c) at the blue line position in (a). The green and red lines on these curves represent the grayscale mean levels of aluminum and magnesium, respectively.

Download Full Size | PDF

Sample 3 is a cylinder of aluminum foam. The scan voltages are 110 and 120 kV, and the currents and exposure times are 55, 50 uA and 454, 400 ms, respectively. The SDD is 780.58 mm, and the SOD is 260 mm. The energy range [0, 120] was equally divided into 24 intervals; that is, the number of narrow-energy-width spectra is R = 24, and let E_r = r × 5. The regularization parameters were set to α = 1 and β_k = 1 during the decomposition process. The directly reconstructed image at 110 kV and the 12th reconstructed component are shown in Figs. 5(a)–5(b). Figures 5(c)–5(d) show the grayscale curves of the two reconstructed images at the blue line position shown in Fig. 5(a).

Fig. 5. The directly reconstructed image (a) at 110 kV and the 12th reconstructed component (b) for sample 3 decomposed by the NMF -GN algorithm. (c)-(d) show the grayscale curves of the two reconstructed images at the blue line position in (a).

Download Full Size | PDF

Similar to the results for sample 1, the cupping artifacts in the reconstructed components of samples 2 and 3 are weakened, and the gray values of the same material change relatively steadily (Figs. 4 and 5).

Bean hardening artifacts cause large fluctuations in the gray values of the same material area, and the standard deviation (STD) can measure the fluctuation and then quantify the degree of hardening artifacts. Because there are differences in the gray levels of different reconstructed images, the coefficient of variation (CV) was selected instead of the STD. In the case of similar noise, the smaller the CV, the lower the degree of hardening artifacts. Tables 2 and 3 show the mean values, STDs, and CVs of each material region in the partially reconstructed images of the three samples. The symbol 50kV represents the directly reconstructed image at 50kV. The last two columns in Table 3 show the numerical results of sample 3. The CVs after decomposition are significantly smaller than those in the directly reconstructed images (Tables 2–3).

Table 2. The means, STDs and CVs of each material region in the reconstructed images of sample 1.

View Table | View all tables in this article

Table 3. The means, STDs and CVs of each material region in the reconstructed images of samples 2 and 3. The last two columns show the numerical results of sample 3.

View Table | View all tables in this article

To evaluate the relationship between the reconstructed components and the narrow energy spectra, Fig. 6 shows the curves of the LAC of each material as a function of energy. Figures 6(a) and 6(b) show the LAC curves of the materials Si and Al in sample 1. Figures 6(c) and 6(d) show the LAC curves of the materials Mg and Al in sample 2. The blue lines represent the true values of the LAC at the middle energy of narrow-energy-width spectra. The red lines represent the mean LACs of reconstruction components corresponding to narrow-energy-width spectra. There is a small error between the true value and the calculated value of the LAC, and the change trend of the two is basically the same (Fig. 6).

Fig. 6. The LAC curves of each material as a function of energy. (a)-(b) show the LAC curves of the materials Si and Al in sample 1. (c)-(d) show the LAC curves of the materials Mg and Al in sample 2. The blue lines represent the true values of the LAC at the middle energy of narrow-energy-width spectra. The red lines represent the mean LACs of reconstruction components corresponding to narrow-energy-width spectra.

Download Full Size | PDF

4. Discussion

The scale of hardening artifacts in reconstructed images is an important feature for distinguishing polychromatic projections from narrow-energy-width projections. In the above experiments, the lower levels of hardening artifacts in reconstructed components conform to the characteristics of the reconstructed image of the narrow-energy-width projection (Figs. 2–5, Tables 2–3). The contrast of materials with approximately LACs was improved due to the reduction in hardening artifacts (Figs. 2(f), 2(h) and 3(e)). The formula E_r,k${\in} $[E_r-₁, E_r] holds for all reconstructed components in simulation. In fact, E_r,k${\in} $[E_r-₁, E_r] indicates that the r-th reconstructed component is energy-dependent on the r-th narrow-energy-width spectrum. The error between the true value and the calculated value of the LAC in each narrow-energy-width spectrum is relatively small, and the trends of the true value and the calculated value are basically the same. This fact verifies that each projection component obtained by decomposition has a clear orientation of narrow-energy-width spectrum, which conforms to the characteristics of narrow-energy-width projection.

The experiments on samples 1 and 2 showed that the proposed algorithm is effective regardless of whether the LACs of the materials are very different.

The above experiments verified that compared with the methods in [15] and [16], each decomposition projection by the proposed algorithm has a clear narrow-energy-width spectrum orientation. This is mainly attributed to the material prior in the X-ray multispectral forward model. A more accurate optimization function is established by using the Poisson statistical properties of the measured data. The addition of the regularization term ${\cal G}({\textbf d} )$ in the optimization function also suppresses the noise amplification of the reconstructed components (Figs. 2(f)–2(h), 3(e)–3(f), 4(e)–4(f) and 5(d)). In addition, employing the GN algorithm to solve d is the key to obtaining effective decomposition results by the proposed algorithm. In the framework of the block coordinate descent approach, the objective function with respect to d is nonlinear, while the GN algorithm is suitable for solving nonlinear optimization problems and has second-order convergence speed. Nonconvex optimization is a common problem in spectral CT, especially in the case of unknown X-ray spectra. Therefore, a good initialization approach is necessary. Compared with random initialization in [15], the initialization approach adopted in this paper is more efficient.

However, for some reconstructed components in real experiments, although their hardening artifacts have been effectively suppressed, the LACs of the same material change slightly. This may be due to the interference of scattering, electronic read-out noise and other factors in the actual CT imaging, which were not estimated in this paper. To further improve the accuracy of decomposed projections, these factors can be estimated by referring to the literature [26]. In addition, in order to decompose complex objects more effectively, we intend to introduce basic effect decomposition into the framework proposed in this paper.

5. Conclusion

We developed an X-ray spectral CT blind separation algorithm for obtaining narrow-energy-width projections from measured polychromatic attenuation data. Since the spectrum is updated as an unknown with other variables, detailed information about the spectrum is not required. Furthermore, the material prior is introduced into the X-ray multispectral forward model and is critical to determining the energy orientation of the narrow-energy-width projection. Using the statistical properties of the measured data, a regularized weighted least squares optimization problem with constraints is established. An alternating algorithm called NMF -GN is combined with the proposed initialization approach to solve this nonconvex optimization problem. Simulations and real experiments show that the reconstructed components obtained by decomposition have smaller hardening artifacts and that the projected components are energy-dependent on the narrow-energy-width spectrum. This indicates that the projected components have the characteristics of a true narrow-energy-width projection. The proposed algorithm suppresses the noise amplification of the projected components, determines their energy orientation and improves the accuracy of obtaining the narrow-energy-width projections compared with previous methods. However, in some reconstructed components in actual experiments, the LACs of the same material vary slightly, which indicates that the accuracy of narrow-energy-width projections needs to continue to improve. In addition, to decompose complex objects more effectively, the basic effect decomposition can be introduced into the decomposition framework of this paper.

Funding

National Natural Science Foundation of China (61801437, 61871351, 61971381); Natural Science Foundation of Shanxi Province (201801D221206, 201801D221207).

Disclosures

The authors declare no conflicts of interest.

References

1. 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. Med. Imaging 38(10), 2469–2481 (2019). [CrossRef]

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

3. R. Gu and A. Dogandzic, “Blind X-Ray CT Image Reconstruction From Polychromatic Poisson Measurements,” IEEE Trans. Comput. Imaging 2(2), 150–165 (2016). [CrossRef]

4. 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. Med. Imaging 37(1), 68–80 (2018). [CrossRef]

5. 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(9), e174–e187 (2017). [CrossRef]

6. R. Foygel Barber, E. Y. Sidky, T. Gilat Schmidt, and X. Pan, “An algorithm for constrained one-step inversion of spectral CT data,” Phys. Med. Biol. 61(10), 3784–3818 (2016). [CrossRef]

7. B. Chen, Z. Zhang, E. Y. Sidky, D. Xia, and X. Pan, “Image reconstruction and scan configurations enabled by optimization-based algorithms in multispectral CT,” Phys. Med. Biol. 62(22), 8763–8793 (2017). [CrossRef]

8. W. Wu, Q. Wang, F. Liu, Y. Zhu, and H. Yu, “Block matching frame based material reconstruction for spectral CT,” Phys. Med. Biol. 64(23), 235011 (2019). [CrossRef]

9. C. Mory, B. Sixou, S. Si-Mohamed, L. Boussel, and S. Rit, “Comparison of five one-step reconstruction algorithms for spectral CT,” Phys. Med. Biol. 63(23), 235001 (2018). [CrossRef]

10. C. O. Schirra, E. Roessl, T. Koehler, B. Brendel, A. Thran, D. Pan, M. A. Anastasio, and R. Proksa, “Statistical reconstruction of material decomposed data in spectral CT,” IEEE Trans. Med. Imaging 32(7), 1249–1257 (2013). [CrossRef]

11. J. P. Schlomka, E. Roessl, R. Dorscheid, S. Dill, G. Martens, T. Istel, C. Bäumer, C. Herrmann, R. Steadman, G. Zeitler, A. Livne, and R. Proksa, “Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography,” Phys. Med. Biol. 53(15), 4031–4047 (2008). [CrossRef]

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

13. S. Chang, M. Li, H. Yu, X. Chen, S. Deng, P. Zhang, and X. Mou, “Spectrum Estimation-Guided Iterative Reconstruction Algorithm for Dual Energy CT,” IEEE Trans. Med. Imaging 39(1), 246–258 (2020). [CrossRef]

14. S. Haykin and Z. Chen, “The cocktail party problem,” Neural Comput. 17(9), 1875–1902 (2005). [CrossRef]

15. J. Wei, Y. Han, and P. Chen, “Improved contrast of materials based on multi-voltage images decomposition in X-ray CT,” Meas. Sci. Technol. 27(2), 025402 (2016). [CrossRef]

16. J. Wei, Y. Han, and P. Chen, “Narrow-Energy-Width CT Based on Multivoltage X-Ray Image Decomposition,” Int. J. Biomed. Imaging 2017, 1–9 (2017). [CrossRef]

17. Y. X. Zhao, Y. Han, and P. Chen, “Blind Separation of Multi-Voltage Projection Sequence Based on Fundamental Effect Decomposition,” Guang Pu Xue Yu Guang Pu Fen Xi/Spectroscopy Spectr. Anal. 39, 2763–2768 (2019). [CrossRef]

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

19. 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(3), 1758–1789 (2013). [CrossRef]

20. D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” Advances in Neural Information Processing Systems 13, 556–562 (2001).

21. M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol. 50(10), 2365–2386 (2005). [CrossRef]

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

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

24. N. Andrei, “An acceleration of gradient descent algorithm with backtracking for unconstrained optimization,” Numer. Algorithms 42(1), 63–73 (2006). [CrossRef]

25. P. C. Hansen, “Analysis of Discrete Ill-Posed Problems by Means of the L-Curve,” SIAM Rev. 34(4), 561–580 (1992). [CrossRef]

26. A. P. Colijn and F. J. Beekman, “Accelerated simulation of cone beam X-Ray scatter projections,” IEEE Trans. Med. Imaging 23(5), 584–590 (2004). [CrossRef]

r	[E_r_-1, E_r]	${\bar{g}}_{r, 1}$	E_r _,1	${\bar{g}}_{r, 2}$	E_r _,2
1	[0,5]	224	3.01	210	3.00
3	[10,15]	3.46	13.05	3.12	13.02
5	[20,25]	0.671	23.08	0.608	23.03
7	[30,35]	0.256	33.09	0.237	33.04
9	[40,45]	0.137	43.05	0.131	43.05
11	[50,55]	0.0912	52.93	0.0899	53.05
13	[60,65]	0.0696	62.70	0.0705	63.03
15	[70,75]	0.0578	72.36	0.0598	73.01
17	[80,85]	0.0507	81.90	0.0533	82.98
18	[85,90]	0.0481	86.63	0.0509	87.96

component material	50 kV		6		12
component material	Si	Al	Si	Al	Si	Al
mean	0.3339	0.3225	0.3621	0.3427	0.0754	0.0738
STD	0.0208	0.0147	0.0087	0.0083	0.0023	0.0019
CV	0 . 0671	0.0509	0.0241	0.0241	0.0304	0.0259

component material	60 kV		8		16		110 kV	12
component material	Mg	Al	Mg	Al	Mg	Al	Al Foam	Al Foam
mean	0.1512	0.2439	0.0916	0.1551	0.0326	0.0512	0.0810	0.0616
STD	0.0106	0.0181	0.0039	0.0053	0.0015	0.0020	0.0093	0.0045
CV	0 . 0702	0.0740	0.0428	0.0339	0.0460	0.0395	0.1146	0.0726

r	[E_r_-1, E_r]	${\bar{g}}_{r, 1}$	E_r _,1	${\bar{g}}_{r, 2}$	E_r _,2
1	[0,5]	224	3.01	210	3.00
3	[10,15]	3.46	13.05	3.12	13.02
5	[20,25]	0.671	23.08	0.608	23.03
7	[30,35]	0.256	33.09	0.237	33.04
9	[40,45]	0.137	43.05	0.131	43.05
11	[50,55]	0.0912	52.93	0.0899	53.05
13	[60,65]	0.0696	62.70	0.0705	63.03
15	[70,75]	0.0578	72.36	0.0598	73.01
17	[80,85]	0.0507	81.90	0.0533	82.98
18	[85,90]	0.0481	86.63	0.0509	87.96

component material	50 kV		6		12
component material	Si	Al	Si	Al	Si	Al
mean	0.3339	0.3225	0.3621	0.3427	0.0754	0.0738
STD	0.0208	0.0147	0.0087	0.0083	0.0023	0.0019
CV	0 . 0671	0.0509	0.0241	0.0241	0.0304	0.0259

Spectral CT imaging method based on blind separation of polychromatic projections with Poisson prior

Abstract

1. Introduction

2. Method

2.1 X-ray multispectral forward model

2.2 Regularized maximum likelihood estimation with constraints

2.3 Optimization algorithm

3. Results

3.1 Simulation example

3.2 Real data examples

4. Discussion

5. Conclusion

Funding

Disclosures

References

Cited By

Figures (6)

Tables (3)

Equations (32)

Optics Express