Combining Acceleration Techniques for Low-Dose X-Ray Cone Beam Computed Tomography Image Reconstruction

Background and Objective Over the past decade, image quality in low-dose computed tomography has been greatly improved by various compressive sensing- (CS-) based reconstruction methods. However, these methods have some disadvantages including high computational cost and slow convergence rate. Many different speed-up techniques for CS-based reconstruction algorithms have been developed. The purpose of this paper is to propose a fast reconstruction framework that combines a CS-based reconstruction algorithm with several speed-up techniques. Methods First, total difference minimization (TDM) was implemented using the soft-threshold filtering (STF). Second, we combined TDM-STF with the ordered subsets transmission (OSTR) algorithm for accelerating the convergence. To further speed up the convergence of the proposed method, we applied the power factor and the fast iterative shrinkage thresholding algorithm to OSTR and TDM-STF, respectively. Results Results obtained from simulation and phantom studies showed that many speed-up techniques could be combined to greatly improve the convergence speed of a CS-based reconstruction algorithm. More importantly, the increased computation time (≤10%) was minor as compared to the acceleration provided by the proposed method. Conclusions In this paper, we have presented a CS-based reconstruction framework that combines several acceleration techniques. Both simulation and phantom studies provide evidence that the proposed method has the potential to satisfy the requirement of fast image reconstruction in practical CT.


Introduction
Based on the theory of compressive sensing (CS) [1,2], near optimal computed tomography (CT) image can be reconstructed from very few projection data. This new methodology indicates a potential for substantially reducing radiation dose. Over the past decade, many CS-based image reconstruction methods have been shown to improve image quality in low-dose CT [3][4][5][6][7][8]. Such methods are often referred to as total variation (TV). To solve the TV problem in the field of CT reconstruction, a two-step alternating minimization framework is commonly used. In the first step, a general iterative reconstruction algorithm is performed to reduce data discrepancy. In the second step, the TV of the reconstructed image is minimized by a gradient descent method. Despite a great improvement of image quality with CS-based reconstruction methods, both high computational load and slow convergence rate limit their practical use.
Recently, the one-step minimization scheme such as firstorder gradient-projection backtracking-line search method [9] and gradient-projection Barzilai-Borwein method was found to converge faster than two-step alternating minimization scheme [9][10][11]. However, there are various ways to improve the convergence rate of the two-step alternating minimization scheme. For example, image reconstruction based on ordered subsets (OS) of projection data is a common acceleration technique used in emission tomography [12] and CT [13][14][15]. Using OS acceleration [12], data discrepancy 2 BioMed Research International can be reduced rapidly compared with non-OS methods. In addition to OS-type simultaneous algebraic reconstruction technique (SART) [15], other faster methods such as ordered subsets transmission (OSTR) algorithm [13] and ordered subsets convex (OSC) algorithm [14] were derived previously. It was also reported that OS-type reconstruction methods can be further accelerated using a bigger step size or a power factor [16][17][18]. Previous studies showed that accelerated OStype algorithms using a power factor can converge three or even four times faster than conventional OS-type algorithms [16][17][18]. These methods, although originally used in emission tomography, can be directly applied for the CT reconstruction. Besides reducing data discrepancy, many different minimization techniques, including fast iterative shrinkage thresholding algorithm (FISTA) [19], TV minimization with dual dictionaries [20], anisotropic TV minimization [21], total difference (TD) with soft-threshold filtering (STF) [6,8], weighted TD (WTD) minimization with STF [22], and TV minimization with half-threshold filtering [23], can be introduced to improve the convergence of the two-step alternating minimization scheme.
Although many acceleration techniques have been developed previously, a combination of acceleration techniques has not been completely studied. To investigate this, we present a reconstruction framework that applies some abovementioned acceleration techniques to the two-step alternating minimization. Specifically, we implemented the TD minimization with STF (TDM-STF), which is one type of TV [6,8,24]. In the TDM-STF, the OSTR algorithm was chosen to accelerate the convergence. To further speed up the convergence of the proposed method, we applied the power factor [16][17][18] and the FISTA algorithm [19] to OSTR and TDM-STF, respectively. This study is different from our recent work [25] that investigated the feasibility of using the power factor [16][17][18] to accelerate the TV-based reconstruction [5]. The purpose of this paper is to study whether combining these techniques can further accelerate the convergence of the two-step alternating minimization. We used simulation and phantom data to evaluate the performance of the proposed algorithm.

CT Image Reconstruction Problem.
According to the idea of CS [1,2] and the TDM algorithm proposed by Yu and Wang [6], the CT image reconstruction problem is to solve the constrained convex optimization problem of the following form: where is the image estimate, is the system matrix, is the measured projection data, is a regularization factor, and TD( ) denotes the total difference of the image estimate defined as [6] TD ( ) = ∑ , , , , − +1, , , , in (2) is called a discrete difference transform (DDT) [6,8]. Similar to the TV problem [3][4][5], the TD problem in (1) can be solved iteratively using a two-step alternating minimization scheme [6,8]. The TDM-STF algorithm [6] was a two-step alternating minimization algorithm. In this study, we used the TDM-STF algorithm [6] to minimize (1). Our aim is to improve the convergence rate of the TDM-STF algorithm [6]. The TDM-STF algorithm and its accelerated algorithm are summarized in the following three sections, respectively. Finally, the implementation of the proposed reconstruction algorithm is summarized in a pseudocode.

TDM-STF and TDM-STF-FISTA.
The STF algorithm proposed by Daubechies et al. [26] was originally developed to solve the linear inverse problems regularized by a sparsity constraint. Yu and Wang [6] adapted the STF method for CT image reconstruction and developed the TDM-STF algorithm to minimize (1). In brief, the TDM-STF method involves three steps. In the first step, the data-fidelity term (i.e., ‖ − ‖ 2 2 ) was minimized via a typical iterative reconstruction algorithm. In the second step, a soft-threshold filtration was performed on the DDT of the current image estimate (i.e.,̂, +1 ). In the third step, the inverse of DDT was computed to obtain image estimate. However, DDT is not invertible [6]. Instead, the second and third steps of the TDM-STF algorithm can be performed based on a pseudoinverse of DDT [6]: where , , and denote the three-dimensional location of the voxel and where is a threshold value. As pointed out by Liu et al. [8], the TDM-STF method can be further accelerated by using a portion of FISTA algorithm [19], which is performed by the following steps: , This accelerated technique is also called Nesterov's momentum algorithms [27,28]. Note that 0 = 1 and old = 0 when = 0.

Acceleration of OSTR Using a Power Factor ℎ.
In the first step of the TDM-STF algorithm, the data-fidelity term was minimized via a typical iterative reconstruction algorithm. To speed up the data-fidelity minimization, we used the OSTR algorithm [13]. The OSTR algorithm is chosen because it is a fast, efficient, and easily implemented algorithm [13]. Moreover, the OSTR algorithm [13] can be accelerated by a power factor ℎ [16][17][18]. The accelerated OSTR (AOSTR) algorithm can be expressed as follows: where , is the estimated attenuation coefficient at voxel and at the th iteration and th subiteration, is the measured projection data at detector bin , is the blank scan at detector bin , = ∑ , and is the number of subsets. Note that the AOSTR algorithm becomes the OSTR algorithm when power factor ℎ = 1. Using the Taylor series expansion, (9) can be approximated as follows: More details can be found in [16][17][18]. Interestingly, the AOSTR algorithm is the same as the OSTR algorithm with a fixed step size. However, in order to preserve the total counts of the forward projections [16][17][18], the reconstructed image updated by (10) is rescaled by the following equation: Note that the solution for (1) from the OSTR algorithm and its accelerated version is not exact, but approximate [29]. An initial condition of uniform image ( 0,0 ) with a value of 0.0002 was used for all reconstructions.

AOSTR-TDM-STF-FISTA.
In summary, we applied the AOSTR reconstruction rather than the conventional iterative reconstruction method in the first step of TDM-STF with FISTA, and the proposed method was called AOSTR-TDM-STF-FISTA. In fact, other combinations of OSTR, AOSTR, and FISTA into the TDM-STF method are possible. For example, the accelerations of OSTR-TDM-STF using a power factor on OSTR and FISTA on TDM-STF are called AOSTR-TDM-STF and OSTR-TDM-STF-FISTA, respectively. In addition to the proposed AOSTR-TDM-STF-FISTA, difference combinations of the acceleration methods are also explored in this study.

Implementation of AOSTR-TDM-STF-FISTA.
Note that the implementation of AOSTR-TDM-STF-FISTA requires considerably more computation per iteration as compared with OSTR. This is due to the computation of the rescaling factor (i.e., (11)) and the TDM step (i.e., (3)) at each subiteration. To reduce the computation time, the TDM step was applied after the last subiteration of each iteration, indicating that STF is performed only once at each iteration of the AOSTR algorithm. Because of this implementation, the rescaling step can be combined with the forward projection of the next subiteration except for the last subiteration [16][17][18]. This means that the rescaling step is computed only at the end of each iteration (i.e., the last subiteration). Similarly, the FISTA algorithm is performed once per iteration rather than once per subiteration. Such modifications make the proposed method an efficient approach for CT reconstruction. The final practical and efficient implementation of AOSTR-TDM-STF-FISTA can be summarized in Pseudocode 1.
As reported by Kim et al. [28], OS-type reconstruction algorithms combined with FISTA become unstable when a large number of subsets were used. To prevent an unstable convergence behavior, FISTA used in the proposed algorithm was run for the first ten iterations only. Moreover, to make the proposed algorithm more stable and computationally efficient, the bit-reversal order of subsets [28,30] was applied instead of the traditional (sequential) order of subsets.

Simulation Study.
To investigate whether OSTR can be accelerated by a power factor ℎ, Figure 2 shows RRMSE values obtained using OSTR and AOSTR (ℎ = 1.5, 2.0, and 2.9) at different iterations. We used 30 subsets for all algorithms. As seen in the comparison, AOSTR can reach lower RRMSE values faster than OSTR. This indicates that the present AOSTR algorithm can converge faster than the OSTR algorithm. The result also shows that the power ℎ of the AOSTR algorithm can be up to 2.9. Furthermore, as illustrated in Figure 3, the present AOSTR algorithm (ℎ = 2.9) provides faster reconstruction than the OSTR algorithm, regardless of numbers of subsets. However, due to the illposed nature of the reconstruction process, the image noise increases with the number of iterations. As a result, the RRMSE values quickly increase as the number of iterations is increased. In particular, the RRMSE value of the AOSTR algorithm with appropriate parameters increases earlier than others because it converges faster than the OSTR algorithm.   reconstructions, the images from the proposed AOSTR-TDM-STF-FISTA algorithm are visually much closer to the true image, reflecting faster convergence. In addition, the AOSTR-TDM-STF algorithm is slightly slower than the AOSTR-TDM-STF-FISTA algorithm. Reconstruction images of chest and abdomen were also shown in Figures 6 and 7.
In Table 1, we compared the computational time required by each algorithm in C program on a Linux workstation with 64 GB memory and an eight-core Intel i7-5960 3.0 GHz CPU. Note that the computation time per iteration of AOSTR-TDM-STF-FISTA is slightly higher (i.e., 6%∼10%) than that of OSTR, but it is minor as compared to the higher acceleration provided by the proposed method.

Summary and Discussion
Since the introduction of CT in 1970s, a wide variety of techniques have been developed to reduce radiation dose, scan time, and image reconstruction time while providing sufficient image quality for diagnosis. However, images reconstructed from low-dose CT data are degraded by noise and artifacts. To address this issue, many iterative reconstruction algorithms were proposed [13][14][15]. Over the past few years, TV has rapidly become a popular and powerful tool for low-dose CT image reconstruction [3][4][5][6][7][8]. In general, TV-based reconstruction problems are solved by a two-step alternating minimization scheme [3][4][5][6][7][8], including an iterative reconstruction algorithm for minimizing data discrepancy and a gradient-based search method for solving the TV regularization problem. Despite a substantial improvement in image quality, TV-based reconstruction methods suffer from high computational cost and slow convergence rate.
Motivated by this issue, we presented a CS-based reconstruction framework that combines several techniques to speed up the convergence rate of the two-step alternating minimization algorithm. First, we implemented the  TDM-STF algorithm which was shown to be superior to conventional TV-based algorithms [6,8]. Second, we combined TDM-STF with OSTR [13] instead of OS-SART [15]. Third, we applied the power acceleration scheme [16][17][18] and the FISTA algorithm [19] to OSTR and TDM-STF, respectively. To investigate whether combining the abovementioned techniques can provide powerful acceleration for the two-step alternating minimization, we compared the performance of the proposed AOSTR-TDM-STF-FISTA algorithm and other algorithms including OSTR, AOSTR, In the simulation study, we showed the feasibility of combining the above-described acceleration techniques to improve the convergence rate of the two-step alternating minimization algorithm. The proposed AOSTR-TDM-STF-FISTA algorithm requires less iterations for convergence than all other algorithms. As also shown in Figures 5-7, the proposed AOSTR-TDM-STF-FISTA algorithm can provide better image quality compared with all other algorithms. We also evaluated the performance of the proposed method using experimental phantom data [32], which support the results observed in the simulation study. More importantly, the above-described acceleration techniques can greatly improve reconstruction speed with minor additional computational time (6%∼10%). However, as the size of reconstructed images and projection data increases, the additional computational time can increase accordingly. Fortunately, image reconstruction using graphics processing unit can drastically reduce the reconstruction time [33][34][35].
Although the convergence proof for the proposed AOSTR-TDM-STF-FISTA algorithm is not available, it could rapidly achieve stable and lower RRMSE as compared with other presented algorithms. This suggests that the proposed AOSTR-TDM-STF-FISTA algorithm has the potential to satisfy the requirement of fast image reconstruction in practical CT applications. We also observed that the proposed AOSTR-TDM-STF-FISTA algorithm did not show any unstable convergence behavior. This may be due to the fact that the FISTA algorithm is terminated after the first ten iterations. In fact, we found out that it is not necessary to perform the FISTA algorithm per iteration. This is simply because the proposed AOSTR-TDM-STF-FISTA algorithm was observed to approach a stable solution after several iterations. Alternatively, the FISTA algorithm that combines with a diminishing step size [28] can be performed in each iteration while maintaining a stable convergence behavior. However, the adaptive step-size method requires additional parameters which may make image reconstruction more complicated.
In this study, we applied the power acceleration scheme to OS-type reconstruction algorithms. The highest acceleration can be achieved using the power factor of ℎ = 2.9. This value is close to the upper limit (=3) found in previous results for emission tomography [16,17] and might imply that the upper limit of the power ℎ was not sensitive to imaging systems and other factors including the number of subsets, image objects, and noise levels [16][17][18]. In addition to the power factor ℎ, the threshold value determines the reconstruction result and the convergence rate of the proposed algorithm. For simulated data, we adjusted the threshold value ( ) to generate the best performance in terms of RRMSE. For phantom data, the value of was adjusted in terms of lower image noise while maintaining acceptable spatial resolution similar to or better than that without TV. Unlike the previous studies [6,8] that used the projected gradient method [36] to automatically determine , our study used fixed values of . This is because the automatic determination of the threshold would increase computation cost [8]. More importantly, dynamic adaptation of may not guarantee an optimal image quality [23]. An efficient optimization of the step size is necessary. This is beyond the goal of this study but will be studied in our future work.
In addition to the acceleration techniques used in this study, many techniques for improving the convergence rate of the two-step alternating minimization algorithm are available. For example, one can use improved TV-based reconstruction algorithms such as TV minimization with dual dictionaries [20], anisotropic TV minimization [21], WTD minimization [22], and TV minimization with half-threshold filtering [23]. However, combining these algorithms does not guarantee optimal results. For example, we had applied WTD [22] to our AOSTR-TDM-STF-FISTA algorithm, and there was a slight improvement in terms of image quality (data not shown). However, weighting values calculated using the neighboring voxels would increase computation time considerably, especially in three-dimensional cases. Despite the previous study [22] that showed the advantage of WTD over TD, combining our AOSTR-TDM-STF-FISTA algorithm with WTD may not be beneficial. More importantly, acceleration techniques that add extra parameters to the existing algorithm may lead to unstable and inaccurate results. Further study is possible on exploring this.

Conclusion
We have presented a CS-based reconstruction framework that combines several acceleration techniques. Both simulation and phantom studies show that, by using the proposed method, the convergence rate of the two-step alternating minimization algorithm can be improved substantially.