A spectral interior CT by a framelet-based reconstruction algorithm

Abstract

Reducing radiation dose is an important goal in medical computed tomography (CT), for which interior tomography is an effective approach. There have been interior reconstruction algorithms for monochromatic CT, but in reality, X-ray sources are polychromatic. Using a polychromatic acquisition model and motivated by framelet-based image processing algorithms, in this paper, we propose an interior reconstruction algorithm to obtain an image with spectral information assuming only one scan with a current energy-integrating detector. This algorithm is a new nonlinear iterative method by minimizing a special functional under a polychromatic acquisition model for X-ray CT, where the attenuation coefficients are energy-dependent. Experimental results validate that our algorithm can effectively reduce the beam-hardening artifacts and metal artifacts. It also produces color overlays which are useful in tumor identification and quantification.

Keywords

Spectral interior CT beam-hardening artifacts maximum likelihood soft thresholding iterative reconstruction

1 Introduction

Since X-ray computed tomography (CT) was introduced in 1972, its most popular reconstruction method is the filtered backprojection (FBP) algorithm. But for a good reconstruction the FBP algorithm needs measurements through the whole cross-section, which results in a high radiation dose. As more and more CT scans are performed every year, controlling the radiation dose becomes urgent and necessary to reduce the risk of the radiation-induced cancer. Interior tomography uses the projection data through the region of interest (ROI) to reconstruct the ROI image. Hence interior tomography is an effective way to reduce the radiation dose in the medical CT field. However conventional wisdom is that interior problem does not have a unique solution [23]. Nevertheless, many efforts have been done in this direction for reducing the radiation dose and other important reasons.

Recently, it has been reported that the interior problem is solvable under some conditions. In this context, Noo et al. [24] proved that exact and stable reconstruction of a small region of the ROI can be realized if two ends of the chord are outside the object using the link between the local projections and the Hilbert transform of the ROI along a chord. Then, Defrise et al. in [9] relaxed the conditions of Noo et al. and proved that there can be exact reconstruction of a larger region of the ROI on the chord if at least one end of the chord is outside the object. After that, it was proved that if a subregion in the ROI is known, there can be exact and stable reconstruction of the whole ROI in [5 , 46].

Inspired by the compressed sensing (CS) theory [4, 12], if the ROI is piecewise constant or piecewise polynomial, the ROI can be reconstructed exactly and stably from under-sampled data by the total variation minimization [15 , 45]. Motivated by the convergence analysis of the soft-thresholding [6], Wang et al. developed the simultaneous algebraic reconstruction technique (SART) combined with soft thresholding [33 , 44].

All methods mentioned above do not take the statistical nature of the projection data into account. Hence they are relatively sensitive to the noise and may not work well for the low-dose CT. As for the statistical properties of the projection data for low-dose CT, there are two main models [16, 21], i.e. Poisson model and Non-stationary Gaussian model. In this context, there have been several statistical iterative algorithms [8 , 32]. In [8], De Man et al. proposed an iterative polychromatic algorithm for CT to reduce the beam-hardening artifacts under the Poisson model, which gives a relatively satisfactory result.

After some efforts of several years, there have been some great progresses in monochromatic interior CT [5 , 39–46]. But in reality, the X-ray is polychromatic. If we do not take the spectral information of X-rays into account, it will cause beam-hardening artifacts. Spectral CT (dual-energy CT or multi-energy CT) can reduce this kind of artifacts and produce the spectral information of an object. Now there are mainly two methods to realize dual-energy CT either by using two scans, one with low-energy X-rays, the other with high-energy X-rays, or using one scan by energy-discriminative detectors. One such CT equipment is GE Discovery CT750 HD [1].

Here we give a new reconstruction algorithm for spectral CT assuming only one scan with a current energy-integrating detector, which can provide more information for diagnosis. That is to say, using our algorithm, no hardware of CT scanners needs to be changed or adjusted. Using the polychromatic acquisition model [8, 31] and motivated by framelet-based image processing algorithms [3], in our previous work [34], we proposed a nonlinear CT reconstruction scheme using one standard X-ray scan. This is a three-step iterative reconstruction algorithm by minimizing some special functional with some regularization terms. In [34], we do not report interior reconstruction results. In this paper, we apply our previous algorithm in interior reconstruction for spectral CT.

In Section 2 we will introduce the data acquisition model we use in this paper, framelet-based image processing methods, and our proposed algorithm. Then in Section 3 we will give our experimental results on the synthesized data and the real-world data. Conclusions are given in Section 4.

2 Model and framelets

In this section, we introduce the data acquisition model and framelets as preliminaries for our proposed algorithm, and then present our method.

2.1 Data acquisition model

First, let us go through the monochromatic data acquisition model for X-ray CT.

Suppose the X-ray spectrum is monochromatic. Let L presents an X-ray beam passing through an object Ω, μ (r) be the X-ray attenuation coefficient at the point r ∈ Ω, and I₀ the initial intensity of the beam L. Then the intensity I of the beam L measured by the detector is [23] $I = I_{0} exp (- \int_{L} μ (r) dr) .$ (1)

After taking logarithm on both sides of Equation (1), it is obvious that the reconstruction of the attenuation coefficients is linear.

However, in reality, the X-ray spectrum is polychromatic. Now we introduce the polychromatic acquisition model we use in our proposed method and the corresponding CT reconstruction procedure, which has been introduced in our previous work [34]. Consider an X-ray tube with a spectral curve I₀ (E) as in Fig. 1, which is composed of the continuous part (due to the effect of bremsstrahlung) and the spikes part (from characteristic κ lines for certain atoms).

Let μ (r, E) be the X-ray attenuation coefficient at the point r ∈ Ω at the energy E of the X-rays. The attenuation coefficient can be approximated by [2 , 31]

$μ (r, E) = φ (r) \cdot Φ (E) + θ (r) \cdot Θ (E)$ (2) for the photoelectric component φ (r) and the Compton scatter component θ (r) which are independent of the energy E, where Θ (E) is the Klein-Nishina function $Θ (E) = \frac{1 + ɛ}{ɛ^{2}} (\frac{2 (1 + ɛ)}{1 + 2 ɛ} - \frac{1}{ɛ} ln (1 + 2 ɛ)) + \frac{1}{2 ɛ} ln (1 + 2 ɛ) - \frac{1 + 3 ɛ}{{(1 + 2 ɛ)}^{2}},$ (3) where ɛ = E/511 keV.

Let I₀ (E) be the initial intensity of the beam L at photon energy E. Then the intensity I of the beam L measured by the detector is [23] $I = \int_{0}^{\infty} I_{0} (E) exp (- \int_{L} μ (r, E) dr) dE .$ (4)

From Equation (4), we can see that the reconstruction of μ from I is not a linear problem any more.

In the classical theory, to linearize the reconstruction problem, by mean value theorem of integrals, Equation (4) can be reduced to $I = \int_{0}^{\infty} I_{0} (E) dE \cdot exp (- \int_{L} μ (r, \bar{E}) dr),$ (5) where $\bar{E}$ is the equivalent energy. And denoting $\int_{0}^{\infty} I_{0} (E) dE$ in Equation (5) as I₀, it is actually Equation (1).

Using Equation (5), the FBP reconstruction of $μ (r, \bar{E})$ for r ∈ Ω becomes a linear problem again. However, the linear attenuation coefficients reconstructed by Equation (5) are actually the linear coefficients at the equivalent energy $\bar{E},$ which is regarded to be independent of the energy E. Hence the approximation $μ (r, E) = μ (r, \bar{E})$ does not provide spectral information. As we all know, if we do not take into account the continuous spectrum of X-rays, it will cause beam-hardening artifacts.

As described in [8, 34], in the discrete setting, we denote the photon energy E as E_k, k = 1, 2, ..., K, the observed measurements as y_i, i = 1, 2, ..., M, the linear attenuation coefficients μ_jk, j = 1, 2, ..., J, and k = 1, 2, ..., K, where i, j, k are indexes of the detector, pixel, and energy respectively, and M, J, K are total numbers of the detector, pixel and energy.

Discretizing Equation (2), we express the attenuation coefficients μ_jk as $μ_{jk} = φ_{j} \cdot Φ (E_{k}) + θ_{j} \cdot Θ (E_{k}),$ (6) where $Φ (E_{k}) = \frac{1 / E_{k}^{3}}{1 / E_{0}^{3 \cdot}}$ (7)

Here E₀ is a reference energy (e.g., E₀ = 70 keV), $Θ (E_{k}) = \frac{Θ (E_{k})}{Θ (E_{0})},$ and φ_j and θ_j actually represent the photoeletric component and Compton scatter part at the reference energy E₀. With l_ij the length of intersection between the X-ray i and the pixel j, the ideal measurements ${\hat{y}}_{i}$ should be ${\hat{y}}_{i} = \sum_{k = 1}^{K} b_{ik} \cdot exp (- \sum_{j = 1}^{J} l_{ij} (φ_{j} \cdot Φ (E_{k}) + θ_{j} \cdot Θ (E_{k}))),$ (8) where b_ik is the total intensity detected by the detector i without the scanned object for incident photons of the energy E_k_, and is given by $b_{ik} = I_{ik} \cdot S_{k} \cdot E_{k},$ (9) where I_ik is the incident number of photons, S_k is the detector sensitivity factor, and E_k is the photon energy (keV). From Equation (9), we can see that b_ik is of the same unit as the energy. That is because of the energy-integrating detectors for CT.

It was reported that the observed measurement y_i approximately follows a Poisson distribution [16] with the probability density function $p (y_{i}) = \frac{{\hat{y}}_{i}^{y_{i}} exp (- {\hat{y}}_{i})}{y_{i}!},$ with ${\hat{y}}_{i}$ as the expectation. Since for any i, y_i is an independent Poisson random variable with expectation ${\hat{y}}_{i}$ , the log-likelihood function is derived as $L = \sum_{i = 1}^{M} (y_{i} \cdot log {\hat{y}}_{i} - {\hat{y}}_{i}) .$ (10)

Now the number of unknowns is 2J. To further simplify the problem (to reduce the number of unknowns), De Man et al. [8] applied the least squares fit to φ_j and θ_j, j = 1, 2, ..., J, and obtained an analytic expression of φ and θ with respect to μ at E = 70 keV. Consequently, they reduced the number of unknowns to J.

In our study, we do not further reduce the number of unknowns in order to reduce the approximation, we just substitute the μ_jk in Equation (8) by Equation (6) to get the ideal measurements ${\hat{y}}_{i} = \sum_{k = 1}^{K} b_{ik} \cdot exp (- \sum_{j = 1}^{J} l_{ij} (φ_{j} \cdot Φ (E_{k}) + θ_{j} \cdot Θ (E_{k}))),$ (11)

Substituting Equation (11) in Equation (10), we get our new log-likelihood L (φ, θ) as $\begin{matrix} L (φ, θ) = \sum_{i = 1}^{M} y_{i} \cdot log (\sum_{k = 1}^{K} b_{ik} \cdot exp (- \sum_{j = 1}^{J} l_{ij} (φ_{j} \cdot Φ (E_{k}) + θ_{j} \cdot Θ (E_{k})))) \\ - \sum_{k = 1}^{K} b_{ik} \cdot exp (- \sum_{j = 1}^{J} l_{ij} (φ_{j} \cdot Φ (E_{k}) + θ_{j} \cdot Θ (E_{k}))), \end{matrix}$ (12) where φ and θ are column vectors of φ_j and θ_j for j = 1, 2, ... , J, respectively.

2.2 Framelets

In this section, some preliminaries of tight framelets are presented as follows.

Let’s first introduce some notation. $L_{2} (ℝ)$ denotes the space of square integrable functions on $ℝ$ . $L_{2} (ℝ)$ is the space of square summable sequences on $ℝ$ . <· , · > is the inner product in $L_{2} (ℝ)$ , i.e. $< f, g > = \int_{ℝ} f (x) \bar{g (x)} dx,$

where $\bar{g (x)}$ denotes the complex conjugate of g (x). ∥ · ∥ ₂ = < · , · > ^1/2 is the norm in $L_{2} (ℝ)$ .

We start with the introduction of the well-known wavelet system. Let ψ be a finite subset of $L_{2} (ℝ)$ . A wavelet system X (Ψ) is defined by $X (Ψ) = {ψ_{j, k} (t) = 2^{j / 2} ψ (2^{j} t - k) | ψ \in Ψ, j, k \in ℤ},$ where ψ ∈ Ψ is called mother wavelet. Any function in $L_{2} (ℝ)$ can be represented by the linear combination of the elements in X (Ψ). In this context, the analysis operator (decomposition operator) $D : L_{2} (ℝ) \mapsto l_{2} (ℝ)$ can be defined by Df = { < f, g > } _g∈X(Ψ) for any $f \in L_{2} (ℝ)$ . The synthesis operator $R : l_{2} (ℝ) \mapsto L_{2} (ℝ)$ is defined by Rw = ∑_g∈X(Ψ)w_gg, for any $w_{g} \in l_{2} (ℝ)$ .

A system $X (Ψ) \subset L_{2} (ℝ)$ is called a frame if the analysis operator is bounded above and below, i.e. $C_{1} {∥ f ∥}_{2}^{2} \leq \sum_{g \in X (Ψ)} {| < f, g > |}^{2} \leq C_{2} {∥ f ∥}_{2}^{2} .$ (13)

If C₁ = C₂ = 1, the system X (Ψ) is called a (Parseval) tight frame, i.e., ${∥ f ∥}^{2} = \sum_{g \in X (Ψ)} {| < f, g > |}^{2}, for all f \in L_{2} (ℝ) .$ (14)

This is equivalent to $f = \sum_{g \in X (Ψ)} < f, g > g, for all f \in L_{2} (ℝ) .$ (15)

When the system X (Ψ) is a tight frame, the corresponding elements g ∈ X (Ψ) is called tight framelets [7]. Obviously orthonormal basis is a tight frame, because conditions (14) and (15) hold for any orthonormal basis in $L_{2} (ℝ)$ . Hence orthonormal basis is a special case of tight frames. From Equation (15), we can see that all of the tight framelets g ∈ X (Ψ) may not be a basis for $L_{2} (ℝ)$ , but the closed closure of X (Ψ) is. Consequently the linear independency of tight framelets is not necessary, but the tight frame still keep the unitary property, that is, the synthesis operator R (reconstruction operator) is the adjoint (transpose) of the analysis operator D.

In [10], to construct a framelet system with compact support, one starts with a compactly supported refinable function $φ \in L_{2} (ℝ)$ with a refinement mask ζ_φ, which follows the refinement equation: $\hat{φ} (2 \cdot) = ζ_{φ} \hat{φ}$ , where $\hat{φ}$ is the Fourier transform of φ, ζ_φ is trigonometric polynomial with ζ_φ (0) = 1. And then the compactly supported framelets ψ are defined by $\hat{ψ} (2 \cdot) = ζ_{ψ} \hat{φ}$ in the Fourier domain, where ζ_ψ is a trignometric polynomial. By the unitary extension principle in [29], if ζ_φ and {ζ_ψ } _ψ∈Ψ satisfy $ζ_{φ} (w) \bar{ζ_{φ} (w + γ π)} + \sum_{ψ \in Ψ} ζ_{ψ} (w) \bar{ζ_{ψ} (w + γ π)} = δ_{γ, 0}, γ = 0, 1,$ (16) for almost all $w \in ℝ$ , then a system X (Ψ) forms a tight frame in $L_{2} (ℝ)$ . Here the sequence ζ_ψ are called tight framelet mask (high-pass filters). In short, given a refinement mask ζ_φ, the construction of tight framelets Ψ actually is to compute the tight framelet masks {ζ_ψ } _ψ∈Ψ by Equation (16). When ζ_φ satisfies ${| ζ_{φ} |}^{2} + {| ζ_{φ} (\cdot + π) |}^{2} \leq 1,$ it is easy to compute {ζ_ψ } _ψ∈Ψ [7, 10].

Here we give one example we used in our numerical experiments. We choose the simplest Haar framelets in our experiments. The refinement mask is $h_{0} = \frac{1}{2} [1, 1],$ and the corresponding tight framelet mask is $h_{1} = \frac{1}{2} [1, - 1] .$ This is the one-dimensional form. In order to get the two-dimensional form, we just need to do the tensor product of the masks as follows: $h^{(00)} = h_{0} \otimes h_{0}, h^{(01)} = h_{0} \otimes h_{1}, h^{(10)} = h_{1} \otimes h_{0}, h^{(11)} = h_{1} \otimes h_{1},$ where ⊗ denotes the tensor product defined by a ⊗ b = ab^T, where a and b are column vectors. That is $h^{(00)} = \frac{1}{4} [\begin{matrix} 1 & 1 \\ 1 & 1 \end{matrix}], h^{(01)} = \frac{1}{4} [\begin{matrix} 1 & - 1 \\ 1 & - 1 \end{matrix}], h^{(10)} = \frac{1}{4} [\begin{matrix} 1 & 1 \\ - 1 & - 1 \end{matrix}], h^{(11)} = \frac{1}{4} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] .$

2.3 Our method based on framelets

In this section, we introduce our spectral interior CT by a framelet-based iterative reconstruction algorithm. Using model (2), all we should compute are the photoelectronic coefficients φ and the Compton coefficients θ. We assume that φ and θ are piecewise smooth. This assumption is reasonable, because from Equation (6), we can see that the φ and θ actually represent the attenuation coefficients at the reference energy E₀ and as we all know the linear attenuation coefficients of human body are usually modeled as piecewise smooth.

Using the polychromatic data acquisition model, under the piecewise smooth assumption, and motivated by the framelet-based image inpainting algorithm by Cai et al. [3], we developed an framelet-based reconstruction method for polychromatic CT, which first transforms the data into the framelet domain, and then perform soft thresholding, and finally transform it back to the original domain. This algorithm is firstly developed for whole reconstruction in our previous work [34]. Now we apply it in the interior reconstruction problem. We will introduce this algorithm as follows.

We aim at minimizing the negative log-likelihood function L (φ, θ) [17 , 30], with sparsity constraints. Let Ω be the whole cross-section of the object to be reconstructed and Ω_in the ROI of Ω. The minimization problem becomes $min_{φ, θ \in Ω_{in}} F (φ, θ) s . t . W φ and W θ are sparse,$ (17) where F (φ, θ) = - L (φ, θ) , W denotes the tight frame analysis transform, α = Wφ, β = Wθ are the tight frame coefficients of φ and θ respectively. In (17), the objective function is composed of the smooth part F (φ, θ) and sparsity constraints on Wφ and Wθ. In order to minimize F (φ, θ), we can apply the popular scaled gradient descent method [17 , 30]. And as for the sparsity constraints, we use the soft thresholding method.

All steps are as follows:

First initialize φ, θ as φ⁰, θ⁰ and then do the following two steps.

1. To minimize the smooth part F (φ, θ), we take the scaled gradient descent method $φ^{n + \frac{1}{3}} = φ^{n} + δ_{1} Δ φ^{n}, θ^{n + \frac{1}{3}} = θ^{n} + δ_{2} Δ θ^{n},$ where δ₁ and δ₂ are step sizes, and $\begin{matrix} Δ φ^{n} = - \frac{\partial F}{\partial φ} (φ^{n}, θ^{n}) {(\sum_{h = 1}^{J} (\frac{\partial^{2} F}{\partial φ \partial φ_{h}} + \frac{\partial^{2} F}{\partial θ \partial φ_{h}}) (φ^{n}, θ^{n}))}^{- 1}, \\ Δ θ^{n} = - \frac{\partial F}{\partial θ} (φ^{n}, θ^{n}) {(\sum_{h = 1}^{J} (\frac{\partial^{2} F}{\partial φ \partial θ_{h}} + \frac{\partial^{2} F}{\partial θ \partial θ_{h}}) (φ^{n}, θ^{n}))}^{- 1} . \end{matrix}$

Here we use the superscript $n + \frac{1}{3}$ to show that this is just a transitional solution that needs to be regularized later, and the derivation of vartriangleφⁿand vartriangleθⁿ are given in our previous work [34] in detail. In order to keep the nonnegativity of the photoelectric component and the Compton scatter component of the attenuation coefficients, we project $φ^{n + \frac{1}{3}}$ and $θ^{n + \frac{1}{3}}$ onto the nonnegativity cone, i.e. $φ^{n + \frac{2}{3}} = max {φ^{n + \frac{1}{3}}, 0}, θ^{n + \frac{2}{3}} = max {θ^{n + \frac{1}{3}}, 0} .$

2. For the sparsity constraints on Wφ and Wθ, we use soft-thresholding $α^{n + 1} = T_{λ_{1}} (W φ^{n + \frac{2}{3}}), β^{n + 1} = T_{λ_{2}} (W θ^{n + \frac{2}{3}}),$ (18) where T_λ (·) denotes the soft thresholding operator [11] defined by $T_{λ} (x) = \frac{| x | - 1}{| x |} x$ if |x| > λ, and zero otherwise. Note that (18) is also the solution of ${\begin{matrix} α^{n + 1} = arg min_{α} \frac{1}{2} {∥ α - W φ^{n + \frac{2}{3}} ∥}_{2}^{2} + λ_{1} {∥ α ∥}_{1}, \\ β^{n + 1} = arg min_{α} \frac{1}{2} {∥ β - W θ^{n + \frac{2}{3}} ∥}_{2}^{2} + λ_{1} {∥ β ∥}_{1} \end{matrix}$

After we get αⁿ⁺¹ and βⁿ⁺¹, we then obtain φⁿ⁺¹ and θⁿ⁺¹ by $φ^{n + 1} = W^{T} T_{λ_{1}} (W φ^{n + \frac{2}{3}}), θ^{n + 1} = W^{T} T_{λ_{2}} (W θ^{n + \frac{2}{3}}) .$

In summary, we present our developed algorithm as follows:

Algorithm.

(1) Set initial guess φ⁰, θ⁰.

(2) Iterate on n ${\begin{matrix} \begin{matrix} φ^{n + \frac{1}{3}} = φ^{n} - δ_{1} \frac{\partial F}{\partial φ} (φ^{n}, θ^{n}) {(\sum_{h = 1}^{J} (\frac{\partial^{2} F}{\partial φ \partial φ_{h}} + \frac{\partial^{2} F}{\partial θ \partial φ_{h}}) (φ^{n}, θ^{n}))}^{- 1} \\ φ^{n + \frac{2}{3}} = max {φ^{n + \frac{1}{3}}, 0} \\ θ^{n + \frac{1}{3}} = θ^{n} - δ_{2} \frac{\partial F}{\partial θ} (φ^{n}, θ^{n}) {(\sum_{h = 1}^{J} (\frac{\partial^{2} F}{\partial φ \partial θ_{h}}) (φ^{n}, θ^{n}))}^{- 1} \end{matrix} \\ \begin{matrix} θ^{n + \frac{2}{3}} = max {θ^{n + \frac{1}{3}}, 0} \\ φ^{n + 1} = W^{T} T_{λ_{1}} (W θ^{n + \frac{2}{3}}) \\ θ^{n + 1} = W^{T} T_{λ_{2}} (W θ^{n + \frac{2}{3}}) \end{matrix} \end{matrix}$

(3) The outcome of Step (2), φ^* and θ^*, is a solution of (17), which represents the attenuation coefficients at the reference energyE₀. And then the energy-dependent attenuation coefficients μ (r, E) at any energy can be computed by Equation (2).

3 Experimental results

In this section, we give our results on the Shepp-Logan Phantom and real CT data.

3.1 Experiments with Shepp-Logan Phantom

In this subsection, we present our results using the 2D Shepp-Logan phantom, a widely-used synthesized phantom for testing algorithms in CT field. In our simulations, there are 10 ellipses in the phantoms. To generate our test Phantoms 1 and 2, we place different substances in the ellipses: the most outer ellipsoidal shell is filled by bone, and the part next to it by plexiglas; the two slant ellipses in the middle are air, and the remaining six ellipses are made of aluminium (Phantom 1) and iron (Phantom 2) respectively. The intensities [8] of these different substances filled in the Shepp-Logan phantom are listed in Table 1. Using these data, the underlying photoelectric coefficients φ and Compton scatter coefficients θ are generated for Phantoms 1 and 2 as shown in Figs. 2(a)(b) and 5(a)(b), respectively.

A fan-beam imaging geometry is used to simulate the imaging system, which is similar to the one used by Yu et al. [42]. In this system, the radius of the circular scanning locus is 57 cm, and for interior reconstruction the diameter of the FOV is 12 cm corresponding to 360 detector bins, with each element 0.033 cm. The reconstructed images are of size 256 × 256, with each pixel of a size of 0.078125 cm by 0.078125 cm. In our simulations, we use 180 projection views for a 360° scan for our test Phantoms. That is to say, the system matrix is composed of 672 × 180 rows and 256 × 256 columns. This system matrix is determined by the CT system and is corresponding to l_ij in Equation (11). We simulate the X-ray spectrum by 25 discrete different energies in the range [20 140] keV. All simulated measurements follow Poisson distribution with the mean in Equation (8). Using our proposed algorithm, we can reconstruct the photoelectric coefficients φ and Compton scatter coefficients θ. The step size δ₁ is 0.48 for φ and δ₂ equals 0.9 for θ in the experiment for Phantom 1, and δ₁ = 1 for φ and δ₂ = 0.25 for Phantom 2. We choose the threshold parameters λ₁ and λ₂ empirically and take 5000 and 3000 iterations for Phantoms 1 and 2 when it reaches the stop criterion.

We show the reconstructed images by our algorithm in Figs. 2–7. From Figs. 2(c)–(d), we can see that our method based on one scan can reconstruct the attenuation coefficients of Phantom 1 very well. This shows that not very dense substances, such as soft tissue, bone, and aluminium can be reconstructed well by our algorithm. In Figs. 5(c)–(d), Phantom 2 are also restored very well. Obviously most of beam-hardening artifacts are removed, and just a little noise in the iron area is remained. And we can keep all reconstructed images of our phantom in good shape.

To further show the effectiveness of our method, we plot profiles of the 64th row of the reconstructed φ and θ of Phantoms 1 and 2 in Figs. 4(a)(b) and 7(a)(b). From these figures, we can observe that all the profiles of the reconstructed φ and θ are very close to those of the corresponding simulated phantoms, which can further validate that our method performs well in removing beam-hardening artifacts and metal artifacts. In Figs. 4(a)(b) for Phantom 1, the profile is very steady and close to that of the real data. In Figs. 7(a)(b), there is just a little deviation in the Fe area, which shows that using our method, most of the artifacts are removed.

To do quantitative validation of our method, we calculate the mean square error (MSE) $MSE = \frac{1}{N} \sum_{i = 1}^{m} \sum_{j = 1}^{n} {(f_{ij} - {\hat{f}}_{i j})}^{2},$ where N are total pixel number, f and $\hat{f}$ are pixel value of the reconstructed image and the true image respectively. The MSE against the number of iterations are shown in Figs. 3 and 6, from which one can observe that the MSE drops down very quickly to a steady value, and then it decreases very slightly. This implies that our algorithm converges to a stable solution.

Finally, we do some comparison between our algorithm with the algorithm where the sparsity transform in (17) replaced by the gradient operator. Specifically it is $min_{φ, θ \in Ω_{in}} F (φ, θ) s . t . D φ and D θ are sparse,$ (19)

where D is the gradient transform. In the following we will call this algorithm as TV method. Using the TV method, we do the reconstruction using the same amount of projection views as those in our method. The reconstruction result of the TV method is shown in Figs. 2(f) and 5(f). From Figs. 2(f) and 5(f), we can see there are still many obvious artifacts in the reconstructed images, especially in the Fe region of Phantom 2 (c.f. Fig. 5(f)). The profile figure in Figs. 4(c) and 7(c) can validate this point further. In these two figures, the profile of the TV method is not very steady, and is still with many oscillations, especially in the Fe regions. In contrast to the TV method, the profile of our method is much steadier and smoothier and is much closer to the real data. In addition, in the experiments, it has been found that our method can use much fewer views to obtain a relatively good image than the TV method. Hence it is preferred for reduction of radiation dose. Furthermore, the tight frame in our proposed method captures variations in a more sophisticated fashion than the TV method. Hence, our proposed tight frame method is generally better than the TV method in terms of details and edges in reconstructed images. Consequently, we can draw a conclusion that compared with the TV method, our method performs better in artifacts reduction and details preservation.

3.2 Experiments with real-world measurements

In this section, interior reconstructions of real data will be presented. A real sheep head was scanned using a cone-beam imaging geometry, which is the same data as what we used in our previous work [34]. The sheep head was bought dead at a local meat market. The rotation center was located with the method proposed by Yang et al. [37, 38]. The distances from the X-ray source to the rotation center and the detector are 109.76 cm and 139.83 cm, respectively. 1195 projections are collected under a fast-continuous-rotation scanning mode. Then the central slice data corresponding to 1565 detector bins and 1195 projection views were extracted with a structure similarity (SSIM) based method proposed by Yang et al. [36]. In order to get measurements through the ROI, we delete the left and right parts of the measurements for all projection views, leaving the middle measurements with 513 detector bins. For interior reconstructions using our algorithm, two 1024 × 1024 images were reconstructed using 513 detector bins and 200 equi-angular projection views, and used to obtain any spectral images μ at any energy according to Equation (2).

In our reconstruction, we use 21 spectra from 20 keV to 120 keV of X-rays at 110 kV in Fig. 1. The reconstructed images are shown in Fig. 8. The reconstructed μ at monochromatic 70 keV is shown in Fig. 8(a) and Fig. 8(b) is the corresponding ROI of Fig. 8(a). And Fig. 8(d) is the color μ with color overlay by reconstructed μ at monochromatic 40, 70, 100 keV, which is a real color CT image. The color overlay is performed as follows: Set the reconstructed μ at monochromatic energy levels 40, 70, 100 as Red, Green, and Blue respectively, and then linearly merge them and empirically adjusted the weights. The advantage of the color overlay is the improved visibility of the reconstructed image. This is the main innovative point of our work.

We also compare our method with the TV method and the reconstructed image is given in Fig. 8(c). Comparing the Figs. 8(b) and 8(c), we can see that our method performs better at artifacts reduction and details preservation.

4 Conclusion

A spectral interior CT by a framelet-based reconstruction algorithm is developed in this study. Excellent results on simulations and real data show that our algorithm can prevent beam-hardening artifacts effectively. In addition to reducing beam-hardening artifacts, our algorithm can reduce metal streak artifacts a lot. Our main work is that we can reconstruct any spectral image using just one scan and use the attenuation coefficients at different energies to merge a real-color CT image.

Footnotes

Acknowledgments

The first author was supported in part by the State Scholarship Fund of China (File No. 201406220189). She is grateful to the University of Iowa for hosting a visit during which this paper was prepared.

References

GE Healthcare, Computed Tomography–Revolution HD, http://www3.gehealthcare.com/en/products/categories/computed_tomography/revolution_hd.

Avrin

D.E.

, Macovski

, Zatz

L.M.

Clinical application of compton and photo-electric reconstruction in computed tomography: Preliminary results. Investigative Radiology 13(3) (1978)217–222.

Cai

J.-F.

, Chan

R.H.

, Shen

A framelet-based image inpainting algorithm. Applied and Computational Harmonic Analysis 24(2) (2008)131–149.

Candès

E.J.

, Romberg

, Tao

Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory 52(2) (2006)489–509.

Courdurier

, Noo

, Defrise

, Kudo

Solving the interior problem of computed tomography using a priori knowledge. Inverse Problems 24(6) (2008)065001.

Daubechies

, Defrise

, De

Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. arXiv preprint math/0307152 2003.

Daubechies

, Han

, Ron

, Shen

Framelets: MRA-based constructions of wavelet frames. Applied and Computational Harmonic Analysis 14(1) (2003)1–46.

De Man

, Nuyts

, Dupont

, Marchal

and Suetens

, An iterative maximum-likelihood polychromatic algorithm for CT. IEEE Transactions on Medical Imaging 20(10) (2001)999–1008.

Defrise

, Noo

, Clackdoyle

, Kudo

Truncated Hilbert transform and image reconstruction from limited tomographic data. Inverse Problems 22(3) (2006)1037–1053.

10.

Dong

, Shen

Pseudo-splines, wavelets and framelets. Applied and Computational Harmonic Analysis 22(1) (2007)78–104.

11.

Donoho

D.L.

De-noising by soft-thresholding. IEEE Transactions on Information Theory 41(3) (1995)613–627.

12.

Donoho

D.L.

Compressed sensing. IEEE Transactions on Information Theory 52(4) (2006)1289–1306.

13.

Elbakri

I.A.

, Fessler

J.A.

Statistical image reconstruction for polyenergetic X-ray computed tomography. IEEE Transactions on Medical Imaging 21(2) (2002)89–99.

14.

Elbakri

I.A.

, Fessler

J.A.

Segmentation-free statistical image reconstruction for polyenergetic x-ray computed tomography with experimental validation. Physics in Medicine and Biology 48(15) (2003)2453–2477.

15.

Han

, Yu

, Wang

A general total variation minimization theorem for compressed sensing based interior tomography. Journal of Biomedical Imaging 2009 (2009)1–3.

16.

Hsieh

Adaptive streak artifact reduction in computed tomography resulting from excessive x-ray photon noise. Medical Physics 25(11) (1998)2139–2147.

17.

Hudson

H.M.

, Larkin

R.S.

Accelerated image reconstruction using ordered subsets of projection data. IEEE Transactions on Medical Imaging 13(4) (1994)601–609.

18.

Kole

Statistical image reconstruction for transmission tomography using relaxed ordered subset algorithms. Physics in Medicine and Biology 50(7) (2005)1533–1545.

19.

Kudo

, Courdurier

, Noo

, Defrise

Tiny a priori knowledge solves the interior problem in computed tomography. Physics in Medicine and Biology 53(9) (2008)2207–2231.

20.

Lange

, Carson

et al., EM reconstruction algorithms for emission and transmission tomography. Journal of Computer Assisted Tomography 8(2) (1984)306–316.

21.

, Li

, Wang

, Wen

, Lu

, Hsieh

, Liang

Nonlinear sonogram smoothing for low-dose X-ray CT. IEEE Transactions on Nuclear Science 51(5) (2004)2505–2513.

22.

Mumcuoglu

E.U.

, Leahy

, Cherry

S.R.

, Zhou

Fast gradient-based methods for Bayesian reconstruction of transmission and emission PET images. IEEE Transactions on Medical Imaging 13(4) (1994)687–701.

23.

Natterer

The mathematics of computerized tomography, , volume 32, SIAM. 1986.

24.

Noo

, Clackdoyle

, Pack

J.D.

A two-step Hilbert transform method for 2D image reconstruction. Physics in Medicine and Biology 49(17) (2004)3903–3923.

25.

Nuyts

, DeMan

, Dupont

, Defrise

, Suetens

, Mortelmans

Iterative reconstruction for helical CT: A simulation study. Physics in Medicine and Biology 43(4) (1998)729–737.

26.

O’Sullivan

J.A.

, Benac

Alternating minimization algorithms for transmission tomography. IEEE Transactions on Medical Imaging 26(3) (2007)283–297.

27.

Phelps

M.E.

, Hoffman

E.J.

, Ter-Pogossian

M.M.

Attenuation coefficients of various body tissues, fluids, and lesions at photon energies of 18 to 136keV 1. Radiology 117(3) (1975)573–583.

28.

Rockmore

A.J.

, Macovski

A maximum likelihood approach to emission image reconstruction from projections. IEEE Transactions on Nuclear Science 23(4) (1976)1428–1432.

29.

Ron

, Shen

Affine systems in

L_{2} (ℝ^{d})

: The analysis of the analysis operator. Journal of Functional Analysis 148(2) (1997)408–447.

30.

Schirra

C.O.

, Roessl

, Koehler

, Brendel

, Thran

, Pan

D.Z.

, Anastasio

, Proksa

et al., Statistical reconstruction of material decomposed data in spectral CT. IEEE Transactions on Medical Imaging 32(7) (2013)1249–1257.

31.

Stonestrom

J.P.

, Alvarez

R.E.

, Macovski

A framework for spectral artifact corrections in X-ray CT. IEEE Transactions on Biomedical Engineering 2 (1981)128–141.

32.

Sukovic

, Clinthorne

N.H.

Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography. IEEE Transactions on Medical Imaging 19(11) (2000)1075–1081.

33.

Wang

, Jiang

Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART). Journal of X-ray Science and Technology 12(3) (2004)169–178.

34.

Wang

, Wang

, Mao

, Cong

, Cai

J.-F.

, Ye

A framelet-based iterative maximum-likelihood reconstruction algorithm for spectral CT, preprint.

35.

Yang

, Yu

, Jiang

, Wang

High-order total variation minimization for interior tomography. Inverse Problems 26(3) (2010)035013–3501329.

36.

Yang

, Duan

, Wang

, Li

, Meng

, Zhang

Extra projection data identification method for fast-continuous-rotation industrial cone-beam CT. Journal of X-ray Science and Technology 21(4) (2013)467–479.

37.

Yang

, Gao

, Li

, Meng

, Wei

A new method to determine the center of rotation shift in 2D-CT scanning system using image cross correlation. NDT & E International 46 (2012)48–54.

38.

Yang

, Pan

, Zhang

, Song

S.-J.

, Meng

, Li

, Wei

Center of rotation automatic measurement for fan-beam CT system based on sinogram image features. Neurocomputing 120 (2013)250–257.

39.

, Yu

, Wang

Exact interior reconstruction with cone-beam CT, International Journal of. Biomedical Imaging 2007 (2007)1–5.

40.

, Yu

, Wang

Exact interior reconstruction from truncated limited-angle projection data. International Journal of Biomedical Imaging 2008 (2008)1–5.

41.

, Yu

, Wei

, Wang

A general local reconstruction approach based on a truncated Hilbert transform. International Journal of Biomedical Imaging 2007(1) (2007)1–8.

42.

, Wang

Compressed sensing based interior tomography. Physics in Medicine and Biology 54(9) (2009)2791–2805.

43.

, Wang

A soft-threshold filtering approach for reconstruction from a limited number of projections. Physics in Medicine and Biology 55(13) (2010)3905–3916.

44.

, Wang

SART-type image reconstruction from a limited number of projections with the sparsity constraint. Journal of Biomedical Imaging 2010 (2010)1–9.

45.

, Yang

, Jiang

, Wang

Supplemental analysis on compressed sensing based interior tomography. Physics in Medicine and Biology 54(18) (2009)N425–N432.

46.

, Ye

, Wang

Interior reconstruction using the truncated Hilbert transform via singular value decomposition. Journal of X-ray Science and Technology 16(4) (2008)243–251.