Tensor framelet based iterative image reconstruction algorithm for low-dose multislice helical CT

In this study, we investigate the feasibility of improving the imaging quality for low-dose multislice helical computed tomography (CT) via iterative reconstruction with tensor framelet (TF) regularization. TF based algorithm is a high-order generalization of isotropic total variation regularization. It is implemented on a GPU platform for a fast parallel algorithm of X-ray forward band backward projections, with the flying focal spot into account. The solution algorithm for image reconstruction is based on the alternating direction method of multipliers or the so-called split Bregman method. The proposed method is validated using the experimental data from a Siemens SOMATOM Definition 64-slice helical CT scanner, in comparison with FDK, the Katsevich and the total variation (TV) algorithm. To test the algorithm performance with low-dose data, ACR and Rando phantoms were scanned with different dosages and the data was equally undersampled with various factors. The proposed method is robust for the low-dose data with 25% undersampling factor. Quantitative metrics have demonstrated that the proposed algorithm achieves superior results over other existing methods.


Introduction
X-ray computed tomography (CT) has been one of the most widely used medical imaging techniques since Hounsfield invented the first commercial medical X-ray machine in 1972 [1]. The Helical CT was first invented by I. Mori [2] in the late 1980s and was developed by W. Kalender [3] in the 1990s. The number of detector rows has been increased to achieve larger volume coverage with a reduced scan time and improved z-resolution. The 8-slice CT system was first introduced in 2000, Siemens SOMATOM Definition scanner has 64-slice rows for up PLOS  FDK, the Katsevich and TF algorithms, with low-dose and sparse-view data. Section IV summarizes this work.

Minimization problem
The mathematical formulation of an iterative CT reconstruction can be expressed by a leastsquare minimization problem as x ¼ arg min where x is the three-dimensional image to be reconstructed with given projection data y and the projection matrix A. The first term indicates the data fidelity in the L 2 -norm. The second term consists of R(x) as a regularization function with regularization parameter λ. For example, the TV norm is a popular regularization choice for sparsity-based CT image reconstruction [31,32].
In this paper, we solve Eq (1) with the given data y from the multislice helical CT system. The projection matrix A contains the helical geometry with the flying focal spot [38]. For the forward projection A and its adjoint A T , parallelized algorithms with an infinitely narrow beam are used with GPU implementation [39].
Tensor framelet regularization. Consider a 3D image x as a tensor, where x ijk is the (i, j, k)-th voxel in three-dimensional image space, N x , N y , and N z are the number of voxels in the x, y and z-axis respectively. We define x x , x y , and x z as 1D unfolded matrices of x along the x, y, and z-axes, respectively. The TF transform is constructed using the standard 1D framelet transform [40], e.g., the 1D piecewise linear tight frame with the following refinement masks.
The operator ω 0 is an averaging operator, and the two other operators ω 1 and ω 2 are the first and second differential operators, respectively. Note that ω 0 smoothes the image, while ω 1 and ω 2 enhance the edges. Define where � denotes the convolution operator. The TF regularization function W and its adjoint W T are respectively defined as below. and The TF norm is defined as λkWxk 1 = λ 0 kM 0 xk 1 + λ 1 kM 1 xk 1 + λ 2 kM 2 xk 1 , where ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi , for all j = 0, 1, and 2. TF transform W is left invertible and W T W = I, by the simple calculation [27]. If λ 0 = 0, λ 1 6 ¼ 0, and λ 2 = 0, kWxk 1 corresponds to the isotropic TV norm of x. In other words, TF regularization is a high-order generalization of TV. The TF transform W can be extended to the multilevel by diluting the masks to o l i such that : and its adjoint W T is defined as Similarly, Eqs (4) and (5) are a generalization of TV to multilevel, and it keeps the framelet features such as W T (Wx) = x. With the TF regularization, Eq (1) becomes The TF regularization term is defined as the isotropic shrinkage TF norm [27]: Optimization algorithm. The TF regularization (7) is the summation of L 1 -norm. To solve the non-differentiable L 1 minimization problem (6), we choose the alternating direction method of multipliers (ADMM) [41] or the so-called Split Bregman method [42]. In general it is difficult to solve the L 1 -regularized minimization problem because it has non-differentiable L 1 term. The basic idea of ADMM is to split L 1 and L 2 components by introducing auxiliary variables d, and v. Eq (6) becomes which can be split into three steps: Because of the decoupled form, step 1 is the sum of two differentiable L 2 -norm terms. Thus, we can efficiently solve it from its optimal condition by the conjugate gradient method. Note that TF is more computationally efficient than TV due to W T W = I.
Step 2 can be solved efficiently using the TF shrinkage formula [28].
Step 3 is in its explicit form, thus it is easy to implement.

Experiments
Data acquisition. The multislice helical CT reconstruction quality was evaluated using the American College of Radiology (ACR) CT accreditation phantom (Data Spectrum Corporation. Model: ECT/DLX/P) and the Rando phantom. Siemens SOMATOM Definition 64-slice helical CT scanner was used to generate the helical CT projection data. Details of the scan parameters for ACR phantom were as follows: Various voltage parameters with effective mAs, CTDI vol , and DLP are described in Table 1. For every voltage level, there was a 3.05 s scan time, 0.5 s gantry rotation time, and 64 � 0.6 mm collimation with z-flying focal spot. The helical pitch is set to be p = 1, with 2304 projections per rotation. Image volume resolution is: 2 mm slice thickness and 0.9766 × 0.9766 mm 2 axial resolution. The whole image volume has 512 × 512 × 88 voxels. A 21.6 cm inside diameter cylindrical ACR phantom is used. Parameter details for the Rando phantom scan were as follows: 120kV with 350 effective mAs are used. There was a 17 s scan time and 20 � 0.6 mm collimation with z-flying focal spot. The helical pitch is set to be p = 1, with 4608 projections per rotation. Image volume resolution is: 4 mm slice thickness and 0.9766 × 0.9766 mm 2 axial resolution. The whole image volume has 512 × 512 × 53 voxels.
Quantitative metrics. To evaluate the performance of the proposed algorithm quantitatively in comparison to FDK and the Katsevich algorithm, four different quantitative metrics are selected. The Universal Quality Index (QUI) measures the intensity similarity between the reconstructed and true images. Image noise is measured by Signal-to-Noise Ratio (SNR) and Contrast-to-Noise Ratio (CNR). These two metrics quantify the noise level of the reconstructed images. The Modulation Transfer Function (MTF) is used to evaluate the resolution of the reconstructed images.
Image similarity-Universal Quality Index (UQI). The Universal Quality Index (UQI) [43] was measured to evaluate the similarity between the reconstructed and true images. We considered the image from the scanner to be the true image. Given the ROI within the reconstructed and true images, the associative mean of the image μ, the variance and covariance of μ with the true image μ true over the ROI are denoted as � m, σ 2 , and Cov(μ, μ true ), respectively. The definition of UQI is given as The UQI measures the intensity similarity between two images, and its value ranges [0, 1]. A UQI value close to 1 indicates a better level of similarity between the reconstructed and true images. We chose two ROIs: The whole ACR phantom body on slices 10 and 50. We calculated the UQI scores for all three methods under comparison. Image noise-SNR and CNR. To evaluate the quantitative noise level of the reconstructed images, we chose two different metrics, SNR and CNR. The definitions are as follows. where σ ROI and s ROI air refer to the standard deviations and � m ROI and � m ROI air refer to the mean pixel value in a ROI inside and the background of the phantom, respectively. We chose four Regions Of Interest (ROI) to compare the reconstructed images from all three methods with that from the scanner. For convenience, the CT numbers are normalized with 1 as the maximum.
Image resolution-MTF. The Modulation Transfer Function (MTF) [43,44] is calculated to measure resolution of the reconstructed images. An Edge Spread Function (ESF) was obtained along the profile of the red line on Fig 1. The Line Spread Function (LSF) was achieved by differentiating the ESF. The MTF was obtained from the Fourier transformation of the LSF. Normalization was performed as MTF(0) = 1.

Evaluations with low-dose data
Four evaluations metrics were compared on the different dosage levels of 80, 100, 120, and 140 kVs. A different x-ray source has a different effective dosage (see Table 1). We chose two slices for the evaluation process, slices 10 and 50.

Evaluations with sparse-view data
To evaluate with sparse-view performance, we fixed the dose level at 100kV. Images were reconstructed at four different sampling steps, 1, 4, 8, and 16. The full view data has 2304 views per 360˚. Sampling step 4 was achieved by taking 576 data uniformly per 360˚. Similarly, sampling steps 8 and 16 were achieved with 288 and 144 views per 360˚, respectively. For sampling  algorithm. As shown in the first row, reconstruction images at sampling step 1 are streak-free for all reconstruction algorithms. However, streaks appeared on the images with FDK and Katsevich for sparse-view data. The last column of the Figs 6 and 7 showed that visually TV and TF reconstruction outperformed other two reconstruction methods. On Figs 6(a 0 ) and 7(a 0 ), ROI's are defined as in the previous section. Visual comparison between TV and TF is given in the next subsection.
For the quantitative evaluation of similarity between the reconstructed image and the scanner image, we computed the UQI for each slices 10 and 50. The ROI for the UQI is set as the whole phantom area on a given slice. Fig 8 shows the result of UQI with various sampling step sizes. Both plots (a) and (b) show that the TF algorithm achieved the highest value except one case, which means that the image reconstructed using the TF algorithm was the most similar to the scanner image.
For the quantitative evaluation of the noise level of the reconstructed images, we computed the SNR and CNR on ROIs 1-4. algorithm achieved the highest MTF, especially when the fewest sample generated the highest MTF difference among other reconstructed methods.

Comparison with TV
As shown in Figs 1~7, image qualities of TF and TV are hard to compare. Each quantitative metric shows a slight superiority of TF. To show some good points of the proposed algorithm, we have tested Rando phantom data, which has more realistic and complicated structure than ACR phantom. Rando phantom is scanned and reconstructed with sparse-view as done in the previous subsection. Fig 11 shows the results by TV(top rows) and TF(bottom rows)  Tensor framelet based iterative image reconstruction algorithm algorithms. The sampling size 1 (which is the right column) results shows similar to each other. But as the step size increases, the TV results are more blurry but clean, while that of TF maintains sharpen edges even with a large step size. Same results can be shown in Fig 6.  Streaking artifacts due to partial projection data are shown less in the TF results. As indicated in red box, TF image has more sharpen edges than that of TV. Overall, we can conclude that TV and TF image qualities are similarly good, but TF has more sharpened edge and less artifacts.
One of the key factor to evaluate iterative algorithms is the reconstruction time. TV and TF elapsed times are summarized in Table 2. TF algorithm requires about 25% less time than TV algorithm.

Discussion and conclusion
To summarize, we have successfully developed a GPU-based TF iterative image reconstruction algorithm for low-dose multislice helical CT, and have shown that the TF method provided improved image quality over the FDK, the Katsevich and TF algorithms when dealing with low-dose and sparse-view data, using UQI, SNR, CNR, and MTF measurements as evaluation metrics. High quality images are reconstructed by the proposed algorithm even with partial view data. TF algorithm is more computationally efficient than that of TV, because of the leftinvertibility of the TF transform property [27]. Moreover, TV reconstructed images show more blurry and flattened than TF. The computational complexity of the TF algorithm is O(1), which is the cost of the x-ray transform and its adjoint per parallel thread [27].