A spectral interior CT by a framelet-based reconstruction algorithm

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 ( - ∫ 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 ) · Φ ( E ) + θ ( r ) · Θ ( 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 ) = 1 + ɛ ɛ 2 ( 2 ( 1 + ɛ ) 1 + 2 ɛ - 1 ɛ ln ( 1 + 2 ɛ ) ) + 1 2 ɛ ln ( 1 + 2 ɛ ) - 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 = ∫ 0 ∞ I 0 ( E ) exp ( - ∫ 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 = ∫ 0 ∞ I 0 ( E ) dE · exp ( - ∫ L μ ( r , E ¯ ) dr ) ,(5) where E ¯ is the equivalent energy. And denoting ∫ 0 ∞ I 0 ( E ) dE in Equation (5) as I₀, it is actually Equation (1).

Using Equation (5), the FBP reconstruction of μ ( r , 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 E ¯ , which is regarded to be independent of the energy E. Hence the approximation μ ( r , E ) = μ ( r , 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 · Φ ( E k ) + θ j · Θ ( E k ) ,(6) where Φ ( E k ) = 1 / E k 3 1 / E 0 3 ·(7)

Here E₀ is a reference energy (e.g., E₀ = 70 keV), Θ ( E k ) = Θ ( 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 y ˆ i should be y ˆ i = ∑ k = 1 K b ik · exp ( - ∑ j = 1 J l ij ( φ j · Φ ( E k ) + θ j · Θ ( 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 · S k · 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 ) = y ˆ i y i exp ( - y ˆ i ) y i ! , with y ˆ i as the expectation. Since for any i,  y _i is an independent Poisson random variable with expectation y ˆ i, the log-likelihood function is derived as L = ∑ i = 1 M ( y i · log y ˆ i - 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 y ˆ i = ∑ k = 1 K b ik · exp ( - ∑ j = 1 J l ij ( φ j · Φ ( E k ) + θ j · Θ ( E k ) ) ) ,(11)

Substituting Equation (11) in Equation (10), we get our new log-likelihood L (φ, θ) as L ( φ , θ ) = ∑ i = 1 M y i · log ( ∑ k = 1 K b ik · exp ( - ∑ j = 1 J l ij ( φ j · Φ ( E k ) + θ j · Θ ( E k ) ) ) ) - ∑ k = 1 K b ik · exp ( - ∑ j = 1 J l ij ( φ j · Φ ( E k ) + θ j · Θ ( E k ) ) ) ,(12) where φ and θ are column vectors of φ _j and θ _j for j = 1, 2, ... , J, respectively.

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 > = ∫ ℝ f ( x ) g ( x ) ¯ dx ,

where 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 ) | ψ ∈ Ψ , j , k ∈ ℤ } , 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 ( ℝ ) ↦ l 2 ( ℝ ) can be defined by Df = { <  f, g > } _g∈X(Ψ) for any f ∈ L 2 ( ℝ ). The synthesis operator R : l 2 ( ℝ ) ↦ L 2 ( ℝ ) is defined by Rw = ∑_g∈X(Ψ)w _g g, for any w g ∈ l 2 ( ℝ ).

A system X ( Ψ ) ⊂ L 2 ( ℝ ) is called a frame if the analysis operator is bounded above and below, i.e. C 1 ∥ f ∥ 2 2 ≤ ∑ g ∈ X ( Ψ ) | < f , g > | 2 ≤ C 2 ∥ f ∥ 2 2 .(13)

If C₁ = C₂ = 1, the system X (Ψ) is called a (Parseval) tight frame, i.e., ∥ f ∥ 2 = ∑ g ∈ X ( Ψ ) | < f , g > | 2 , for all f ∈ L 2 ( ℝ ) .(14)

This is equivalent to f = ∑ g ∈ X ( Ψ ) < f , g > g , for all f ∈ 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 φ ∈ L 2 ( ℝ ) with a refinement mask ζ _φ , which follows the refinement equation: φ ˆ ( 2 · ) = ζ φ φ ˆ, where φ ˆ is the Fourier transform of φ,  ζ _φ is trigonometric polynomial with ζ _φ  (0) = 1. And then the compactly supported framelets ψ are defined by ψ ˆ ( 2 · ) = ζ ψ φ ˆ in the Fourier domain, where ζ _ψ is a trignometric polynomial. By the unitary extension principle in [29], if ζ _φ and {ζ _ψ  } _ψ∈Ψ satisfy ζ φ ( w ) ζ φ ( w + γ π ) ¯ + ∑ ψ ∈ Ψ ζ ψ ( w ) ζ ψ ( w + γ π ) ¯ = δ γ , 0 , γ = 0 , 1 ,(16) for almost all w ∈ ℝ, 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 + | ζ φ ( · + π ) | 2 ≤ 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 = 1 2 [ 1 , 1 ] , and the corresponding tight framelet mask is h 1 = 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 ⊗ h 0 , h ( 01 ) = h 0 ⊗ h 1 , h ( 10 ) = h 1 ⊗ h 0 , h ( 11 ) = h 1 ⊗ h 1 , where ⊗ denotes the tensor product defined by a ⊗ b = ab ^T , where a and b are column vectors. That is h ( 00 ) = 1 4 [ 1 1 1 1 ] , h ( 01 ) = 1 4 [ 1 - 1 1 - 1 ] , h ( 10 ) = 1 4 [ 1 1 - 1 - 1 ] , h ( 11 ) = 1 4 [ 1 - 1 - 1 1 ] .

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 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 + 1 3 = φ n + δ 1 Δ φ n , θ n + 1 3 = θ n + δ 2 Δ θ n , where δ₁ and δ₂ are step sizes, and Δ φ n = - ∂ F ∂ φ ( φ n , θ n ) ( ∑ h = 1 J ( ∂ 2 F ∂ φ ∂ φ h + ∂ 2 F ∂ θ ∂ φ h ) ( φ n , θ n ) ) - 1 , Δ θ n = - ∂ F ∂ θ ( φ n , θ n ) ( ∑ h = 1 J ( ∂ 2 F ∂ φ ∂ θ h + ∂ 2 F ∂ θ ∂ θ h ) ( φ n , θ n ) ) - 1 .

Here we use the superscript n + 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 + 1 3 and θ n + 1 3 onto the nonnegativity cone, i.e. φ n + 2 3 = max { φ n + 1 3 , 0 } , θ n + 2 3 = max { θ n + 1 3 , 0 } .

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

After we get αⁿ⁺¹ and βⁿ⁺¹, we then obtain φⁿ⁺¹ and θⁿ⁺¹ by φ n + 1 = W T T λ 1 ( W φ n + 2 3 ) , θ n + 1 = W T T λ 2 ( W θ n + 2 3 ) .

In summary, we present our developed algorithm as follows:

Algorithm.

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

(2) Iterate on n { φ n + 1 3 = φ n - δ 1 ∂ F ∂ φ ( φ n , θ n ) ( ∑ h = 1 J ( ∂ 2 F ∂ φ ∂ φ h + ∂ 2 F ∂ θ ∂ φ h ) ( φ n , θ n ) ) - 1 φ n + 2 3 = max { φ n + 1 3 , 0 } θ n + 1 3 = θ n - δ 2 ∂ F ∂ θ ( φ n , θ n ) ( ∑ h = 1 J ( ∂ 2 F ∂ φ ∂ θ h ) ( φ n , θ n ) ) - 1 θ n + 2 3 = max { θ n + 1 3 , 0 } φ n + 1 = W T T λ 1 ( W θ n + 2 3 ) θ n + 1 = W T T λ 2 ( W θ n + 2 3 )

(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.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 = 1 N ∑ i = 1 m ∑ j = 1 n ( f ij - f ^ i j ) 2 , where N are total pixel number, f and 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 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.

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.