1.

INTRODUCTION

Photon-counting detector (PCD) based spectral CT has attracted an increasing attention [1]. However, a single energy bin contains only a fraction of the total photon, and most PCDs can only accommodate a limited count rate. The multiple projection datasets obtained by PCD usually contain very strong Poisson noise [2]. This makes it difficult to meet the challenges of practical applications.

Inspired by the image domain non-local feature similarity, Zhang et al. extended the traditional vectorized dictionary learning to tensor dictionary learning (TDL) for spectral CT reconstruction. The TDL based reconstruction algorithms for spectral CT have some limitations [3]. When the noise is too large, it is impossible to distinguish the noise from the organizational structures, and block artifacts appear in the reconstructed images. Recently, tensor decomposition methods are widely used in image denoising. A high dimensional tensor can be represented approximately by the weighted sum of a series of low-rank tensor data, leading to effective noise suppression and artifact reduction.

The correlation of reconstructed image in each spectral channel has attracted more and more attention. Rigie and Patrick developed a spectral CT reconstruction method based on vector total variation (VTV) [4], which combines the nuclear norm to promote the sparsity of the multi-channel gradient vector field. Knoll et al. used the second-order total generalized variation as a special multi-channel regularization function (MTGV), and the structural information was shared in the process of reconstruction while the unique differences were retained [5].

In order to make better use of sparsity in the image domain and the correlation of the reconstructed image for multi spectral channels, we propose a spectral CT reconstruction algorithm based on joint multi-channel total generalized variation and tensor decomposition, and we call it MTGV-TD.

2.

MODELS AND METHOD

2.1

Joint multi-channel TGV regularization based on nuclear norm

In 2014, Rigie and Patrick extended the single-channel gradient vector to the vector field and defined the discrete Jacobian matrix in the following form:

where

When two images have the same curve, the two images have the same direction gradient and the converse is also true [6]. If all the gradient vectors of each channel image are parallel or antiparallel, then the rank of the Jacobian matrix will be 1. Hence, there will be only one non-zero singular value. Based on those facts, if the image gradients of each channel are parallel, the nuclear norm will be minimized. Our algorithm will extend the total generalized variation to the vector:

where

A,B ∈ U, σ is the singular value vector of the (∇u − w)(i, j) matrix, and ||·||_Frob is the sum of squares of all elements in the matrix (the Frobenius norm of the matrix). The discrete gradient operator of the vector field is ∇: U^L → U^L×2, The symmetric gradient operator is ε: U^{L×2 →}U^L×4, and ∇u(i, j)= J_u (i, j), w ∈ U^L×2, ε(w) ∈ U ^L×4.

if

then

where , represents the first order backward difference between horizontal and vertical directions, w_ij ∈ U.

2.2

Tensor Decomposition

Tensor decomposition has been widely used in image reconstruction and image processing, etc. Usually, there are two main tensor decomposition methods: the Tucker decomposition and the Candecomp/Parafac (CP) decomposition. Here, we only address the CP decomposition.

CP decomposition can decompose a N^th tensor χ into the weighted sum of rank-one tensors, which can be expressed as

where λ_h is the weight function, a_nh ∈ R¹ⁿ, (n = 1,2,···,N) is a normalized vector and the symbol “°” represents the outer product. The sparsity level of representation can be controlled by adjusting the parameter H. We can apply tensor decomposition for image restoration. The mathematical model can then be expressed as:

where χ₀ and χ are the corrupted and restored images, respectively. χ can be obtained using the ALS method by solving each A_n alternatively [7].

2.3

Mathematical Model for MTGV-TD

In this algorithm, tensor decomposition is used to improve the image blocks quality, and joint multi-channel TGV function is sharing information between channels. Combining multi-channel TGV function and tensor decomposition a spectral CT reconstruction algorithm is proposed. Its objective function can be expressed as the following convex minimization problem:

where u = (u₁, u₂, ··, u_s) , g = (g₁, g₂, ···, g_s), S indicates the number of channels, u_i (i = 1,2,···S) represents the reconstructed image in i^th energy channel, g_i (i = 1,2,···S) represents the projection data in i^th energy channel, C and Z_c denote the group number and group extraction operator.

2.4

Solution

Refer to (7), there are two variables that need to optimization. We divides it into the following two sub-problems and then adopts the method of alternating optimization to solve them:

By introducing variables, equation (7a) is transformed into a constrained problem, which can be summarized as:

Equation (8) can be converted into another unconstrained optimization function:

Equation (9) contains two variables, which are optimized by using the method of alternating iteration. Equation (9) is divided into two sub-problems:

Gradient descent method is adopted to solve equation (10), then:

A first-order Primal-Dual algorithm is used to approximate the global optimal solution of equation (11) [8]. Because each group is independent regarding the optimization process, equation (7b) can be rewritten into

Equation (13) is solved by using ALS method [7].

3.

EXPERIMENTAL RESULTS

The major goal of this paper is to evaluate the performance of the MTGV-TD for spectral CT. The algorithms of SART, TV, TDL and multi-channel TGV(MTGV) are implemented for comparison [9,10]. All the algorithms are implemented in a hybrid mode of Matlab and C++. While the interface is implemented in Matlab, all the extensive computational parts are implemented in C++ and complied via MEX function.

A numerical mouse thorax phantom generated by the MOBY software is used for simulation experiments and 1.2% (by volume) iodine contrast agent is introduced into the blood circulation. The model size is 20 × 20 mm², and the resolution is 512 × 512. The scanning radius is 100 mm, and the length of virtual detector located at the center of the object is 20 mm. There are 320 detector units in total, each of which is 0.0625mm. The voltage is 50kVp and energy spectrum is divided into five energy channels: WE1 = {11keV-26keV}, WE2 = {27keV-30keV}, WE3 = {31keV-34keV}, WE4= {35keV-39keV} and WE5 = {40keV-50keV}, as shown in figure 1. In the simulation experiments, an isometric fan-beam scanning is adopted. For each X-ray path, 50000 photons are assumed emitted from the X-ray source. The projection data with Poisson noises are generated with expectations being the number of photons received in the corresponding noise-free case. The parameters used in this algorithm are λ=50, β=0.8, μ=0.2, C=120.

Figure 1.

Spectrum used for numerical simulation

In this experiment, table 1 shows the comparison results of PSNR, NRMSE and SSIM values between reconstructed images and comparison images in representative channels with full scan 360 projections. It can be seen from table 1 that the reconstruction accuracy of MTGV-TD is superior to the comparison algorithms for all the channels. Figure 2 gives the corresponding reconstructed images of channels 1, 3 and 5.

Figure 2

The reconstruction results of the mouse model form 360 projections. From left to right, the rows are SART, TV, TDL, MTGV and MTGV_TD. From top to bottom, the columns are 1st, 3rd and 5th energy channel.

Table. 1

Quantitative evaluation results of different algorithms

From figure 2, one can see that SART algorithm has strong noise in image reconstruction. Image blurring and detail missing appear in TV algorithm, and it is difficult to distinguish noise from details. The overall denoising ability of TDL is weaker than that of MTGV. The reconstructed image of the MTGV algorithm shows that it can guarantee clearer edges while denoising. The proposed MTGV-TD algorithm in this paper reconstructs the image with clear edges and obvious details, and it achieves better noise reduction effect and detail preservation ability. figure 3 shows the comparison results of PSNR, NRMSE and SSIM values between reconstructed images and comparison images with 360, 180, 120, and 90 uniformly sampled full scan projections. It can be seen that the reconstruction effect will become better with the increase of projection number.

Figure 3

Convergence curves of (a)NRMSE, (b) PSNR and (c) SSIM as a function of iterations under different angle

4.

CONCLUSION

In this paper, we propose a spectral CT reconstruction algorithm based on tensor decomposition and joint multi-channel total generalized variation. The algorithm improves the reconstructed image quality by using the image domain sparse condition and the information correlation among channels. The experiment results show this method can not only effectively suppress the noise, but also protect the edge and detail features of the images.