Low-Dose Cone-Beam CT Iterative Reconstruction via Total Variation and Gradient Total Variation

: The compressed sensing (CS) technique has been utilized to reconstruct Cone-beam computed tomography (CBCT) images via limited projection from undersampled measurements. However, the condition of limited projection is an ill-posed problem. Since the CBCT image itself doesn’t have sparse features, the total variation (TV) transform has been widely adopted in CBCT reconstruction. This method, which penalizes the weight of each voxel at a constant rate regardless of different spatial gradient, may not recover qualified CBCT images from ill-posed projection data. This work presents a new strategy to deal with the deficits stated above by utilizing non-uniform weighting penalization in CBCT reconstruction. The proposed new strategy combines TV and gradient total variation (GTV) for reconstruction in a hybrid weighting penalization way, where the total variation is penalized by the gradient total variation in advance. The proposed penalty not only retains the benefits of TV, including artifact and noise suppression, but also maintains the structures in regions with gradual gradient intensity transition more effectively. This study tested the proposed method by under-sampled projections of 2 objects and 2 experiments (2 digital phantom). We assessed its performance against the OS-SART method, FDK method, conventional TV method and TV+GTV method in the tissue contrast, reconstruction accuracy, and imaging resolution by comparing the root mean squared error (RMSE), the correlation coefficient (CC), the structural similarity (SSIM), and profiles intensity of the reconstructed images. The proposed method produced the reconstructed image with the lowest RMSEs and the highest CCs and SSIMs for each experiment.


Introduction
Cone-beam computed tomography (CBCT) is widely utilized in tumor visualization and localization which is critical for clinical image guided radiation therapy [1], since it is closely related to tumor control and radiation toxicity. However, the cumulative imaging dose of repeated CBCT scans is clinically significant and superfluous ionization radiation exposure can increase the venture of secondary cancer induction [2,3]. Lowering the dose itself not only decreases the amount of x-ray exposure but also reduces the image quality [4,5]. Therefore, more and more researches focus on decreasing the venture of exposed radiation while keeping qualified CBCT image.
Generally, the dose from CBCT scans can be reduced by two ways: one way is to lower the X-ray tube current per projection view and another way is to decrease the number of projection data. Unfortunately, reducing one of the factors or both not only reduces the amount of patient dose but also reduce the CBCT image quality. For example, when using the conventional analytical methods, such as FDK algorithm [6], reducing the number of x-ray projections will cause serious streaking artifacts and lowering exposure level (mAs) per projection will amplify the noise level. To deal with the drawback of the analytical methods, iterative algorithms such as the ART [7], SART [8], and OS-SART [9] have been developed. Although the reconstructed image quality has been improved, it is still very challenging for low-dose CBCT reconstruction.
The compressed sensing (CS) algorithms [10,11] have been developed for lowdose CBCT reconstruction under the assumption that the solution is sparse. As the clinical CBCT image is made up of numerous non-zero elements, which means CBCT image would not be sparse by itself, therefore, a variety of transforms were utilized to make image sparse, such as wavelet transform, curvelet transform, and total variation (TV) transform. For instance, statistical iterative reconstruction (SIR) methods [12,13] based on the TV [14], the Huber [15], the isotropic quadratic [16], was widely used to low-dose CBCT image reconstruction. Generally, the TV-based methods have been widely utilized to CBCT reconstruction for low-dose projections data. Adaptivesteepest-descent (ASD) projection-onto-convex-sets (POCS), based on TV, was designed to reconstruct CBCT image from under-sampled projections data [17]. As a matter of fact, TV's variants have been investigated in the last several decades, them were originally utilized to remove out noise and artifacts as a regularization function [14,[18][19][20][21][22][23][24][25][26][27][28][29][30]. For example, the penalized weighted least-squares (PWLS) based on TV was developed for noise-suppressed and edge-preserved [18,31].
This study has shown that TV-based method was more robust than conventional filtered-backprojection (FBP) method in low-dose CBCT reconstruction [27]. The selection of TV minimization parameters has been investigated with visualization and quantitative assessments [24]. Even though TV removes the overall noise effectively, it attempts to penalize gradient uniformly across the entire CBCT image, which would result in anatomical edge structures inevitably over-smoothed. Therefore, edgepreserving total variation (EPTV) was designed to identify the edges and penalty them with lower weight [32]. In order to improve limitations of TV method, more and more TV-based methods were proposed. A reweighted anisotropic TV method was designed to solve the limited-angle CT image reconstruction [33]. A simultaneous deblurring and iterative reconstruction method based on TV was indicated that CBCT image unsharpness was caused by several factors [20]. Meanwhile, this work combined the TV minimization and the ordered-subsets accelerated with a power factor into low-dose CBCT reconstruction [21]. Despite reducing the imaging dose, low-dose CBCT has not gained widespread application in the clinical. Therefore, reducing the dose while achieving high-quality images for radiotherapy is still urgently demanded in the clinical.
In order to avoid edge structures over-smoothed, edge-preserved methods based on TV have been developed. Edge guided total variation minimization (EGTVM) was used to reconstruct CT image from under-sampled projections data [34]. Adaptiveweighted total variation (AwTV)-POCS was utilized to preserve edges, which assumes that the edges have anisotropic property [35]. For the sake of edges information preserved, adaptive-weighted TV was used in the steepest descent part by adaptive-weighted projection-controlled steepest descent (AwPCSD) method [36]. Furthermore, an optimization algorithm based on TV and alternating direction method (ADM) was designed for sparse reconstruction [37]. A method based on TV and gradient total variation (GTV) was proposed for low-dose CBCT reconstruction [38]. Because TV penalizes weight to voxel at the same rate regardless of different spatial gradient, it may not recover qualified CBCT images from ill-posed projection data, a non-local operator was proposed by [39]. Non-local TV imposes non-uniform weight on a more global area centered on each voxel [40][41][42]. With those TV-based methods proposed, the imaging dose was reduced. However, qualified image is still urgently demanded in the clinical low-dose CBCT application.
In recent years, high-order differential operators were usually used to penalize image, aiming to suppress over-smoothing effects for various inverse problems in image processing [43][44][45]. High-order penalties can essentially preserve the oversharpening of areas with smooth intensity transitions because of their piecewisevanishing property [46,47], which can achieve piecewise-linear solutions that better fit smooth intensity changes. The total generalized variation (TGV) model with symmetric tensors developed by Bredies [48], which involved and balanced image second-order derivatives. The high-order penalties were proposed by Chen [46], which was used to image restoration for image processing. Furthermore, the high-order TV penalties were redesigned by Hu for image restoration [49]. Based on the high-order TV minimization, CT reconstruction algorithm has been proposed [50]. However, high-order penalties can effectively reduce the piecewise-constant effect, but might also introduce additional edge blurry. Besides those methods, other methods are further investigated for low-dose CBCT reconstruction, such as the dictionary learning [51], and Hessian Schatten penalties [52,53]. In addition, based on the deep learning methods [54,55] are also used for the CBCT reconstruction, and good results have been achieved. However, the edge of image still could not be preserved accurately via sparse projections reconstruction and would cause the wrong judgement in the clinical.
Due to the drawback of TV regularization term, which penalizes the edges uniformly, it would result in the edge structures over-smoothed. This work presents a new strategy to deal with the deficits stated above by utilizing non-uniform weighting penalization in CBCT reconstruction. The proposed new strategy combines TV and GTV for reconstruction in a hybrid weighting penalization way, where the total variation is penalized by the gradient total variation in advance. The proposed penalty not only retains the benefits of TV, including artifact and noise suppression, but also maintains the structures in regions with gradual gradient intensity transition more effectively. This study tested the proposed method by under-sampled projections of 2 objects and 2 experiments (2 digital phantom). We assessed its performance against the OS-SART method, FDK method, conventional TV method and TV+GTV method in the tissue contrast, reconstruction accuracy, and imaging resolution by comparing the root mean squared error (RMSE), the correlation coefficient (CC), the structural similarity (SSIM), and profiles intensity of the reconstructed images. The proposed method produced the reconstructed image with the lowest RMSEs and the highest CCs and SSIMs for each experiment.

TV-based reconstruction methods
Then the general mathematical reconstruction model can be expressed as follows: Where ‖ − ‖ 2 2 is the data fidelity term indicating the difference between the generate projection and observe projection data . In the expression above, is an

Gradient TV minimization
Recently, a technique based on gradient total variation (GTV) has been utilized in CBCT reconstruction [38]. With more gradient prior constraints in model, the quantified image can be recovered from under-determined measurements. The CBCT image can be de estimated by solving the following optimization problem.
Where the minimization process consists of the data fidelity term (‖ − ‖ 2 2 ) and regularization function ( ‖ ‖ + ‖ ‖ ). The ‖ ‖ is the gradient total variation regularization term. It can be represented as: Where ‖ * ‖ is derived form ‖ * ‖ , TV achieves the assumption that gradient is approximately sparse in CBCT image, GTV also has sparsity property. The Eq. (3) shows that the gradient magnitude will not be penalized by ‖ ‖ term. However, regardless of different amounts of image gradient, TV regularization term always penalizes the weight to voxel uniformly. Due to the noise, artifacts, and low contrast tissue, TV property might critically exterminate the reconstructed image quality.

The proposed method
In order to solve the problem of the TV penalty gradient uniformly, we propose an improvement model based on GTV for constrained TV regularization reconstruction algorithm. We redesign CBCT reconstruction models as follows: Where ( ) is the penalized TV regularization term, and is the parameter of ( ) regularization term, ‖ ‖ is an object term in Eq. (6), which is penalized by ‖ ‖ , and is the parameter of ‖ ‖ regularization term.
Solving Eq. (5) and (6) by a gradient decent method, we need to choose quantified regularization parameter and . The step-size of ( ) regularizing is controlled by , which can only be positive. The step-size of ‖ ‖ is controlled by , unlike , which can be positive or negative.
In equation (6), as the prior information, the ‖ ‖ is incorporated into the TV regularizing function. Therefore, the penalized TV term ( * ( )) would not any more weight each voxel at the same rate, that means the edge structures information can be preserved. The only drawback is the time cost of each TV iterative will increase during the minimization process.
In Table 1, Algorithm 1 summarizes the entire algorithmic procedure of our proposed process. In this section, the GTV process is carried out two times, and we set up the number of iterations for OS-SART loop and TV loop in Algorithm 1 to be 30.
The developed algorithm is implemented by redesigning the ASD-POCS reconstruction approach.

Image quality evaluation criteria
In this work, to quantify of the reconstructed images, we selected the root mean squared error (RMSE), increase in the correlation coefficient (CC) and structural similarity (SSIM). As shown in (7), (8) and (9). RMSE is defined as: Where the is the reconstructed image, and the ̂ is the reference image. CC is defined as: In (8), the is the standard deviation of the reconstruction image, and ̂ is the standard deviations of the reference image. The value of CC is between -1 and 1, where value=1 is the total positive linear correlation, value=0 is no linear correlation and value=-1 is the total negative linear correlation. SSIM is defined as: In (9), the and ̂ are the average of the reconstruction image and reference image, 2 is the variance of the reconstruction image, ̂2 is the variance of the reference image. The 1 and 2 are constants guaranteed to be non-zero.

Simulation experiments
To evaluate the reconstruction performance of the proposed algorithm, we tested our proposed method on two cases (the Sheep-Logan phantom, the XCAT phantom) and compared it with OS-SART, FDK, TV, and TV+GTV algorithm. In the phantom experiments, Siddon's ray tracing algorithm was used for the projection. The parameters of both phantoms acquisition are listed in the following Table 2. All the experiments were implemented on a personal computer (16 GB memories, 3.6 GHz CPU, Windows 10 64bit system environments).

The shepp-logan phantom reconstruction
In  Figure 1, where the original ( Figure 1(A)), the reconstructed image using OS-SART (Figure 1(B)), the analytical FDK reconstructed image (Figure 1(C)), the reconstructed image using TV ( Figure   1(D)), the reconstructed image using TV+GTV (Figure 1(E)) and the reconstructed image using proposed method (Figure 1(F)) are shown.
OS-SART slice and FDK slice are seriously attacked by the streak artifacts. TV slice and TV+GTV slice can prevent the artifacts very well, however, some edge structures are smoothed out. While the proposed method performs well both in terms of artifact removal and edges protection. In general, the proposed method is much better than the OS-SART and FDK in terms of artifacts and noise suppression and performs better than the TV and TV+GTV for edge structures protection.
In order to analyze the edge-preserving property of the experimental results, Figure 2 shows the voxel curves along the coronal axis at 159th row, 128th slice.
Comparing the several methods, the proposed method curves are the most consistent with the phantom curves, that means more edge structures information preserved. The consequences of the quality evaluation are shown in the Table 3. Consistent with the visual quality, the proposed method achieves the best scores.     Table 4.

The XCAT phantom reconstruction
In this section, a digital XCAT phantom with a matrix size of 512*512*512 voxels was used to evaluate the proposed method. The tests were conducted with 121 uniformly under-sampled projections for the XCAT phantom. The sampling angle is evenly distributed in the range of (0, 2*pi), and the projected signal at each angle contains 512*512 measurements. The reconstructed image has a size of 512*512*512.
The iteration number for the reconstruction algorithms was set to 30. Reconstruction methods of OS-SART, FDK, TV, and TV+GTV, are used for comparison.
Reconstructed slice of the XCAT phantom obtained by several reconstruction methods is shown in Figure 4, where the original (Figure 4(A)), the reconstructed image using OS-SART (Figure 4(B)), the analytical FDK reconstructed image (Figure 4(C)), the reconstructed image using TV (Figure 4(D)), the reconstructed image using TV+GTV (Figure 4(E)) and the reconstructed image using proposed method ( Figure   4(F)) are shown. We can see that the OS-SART method (Figure 4(B)) and the FDK method (Figure 4(C)) are not able to reconstruct the CBCT images with few projections and obvious streaking artifacts are observed. In contrast, both the TV method ( Figure   4(D)) and the TV+GTV method (Figure 4(E)) can still achieve most of the edge structures, leading to visually much better reconstruction results. However, one drawback of the TV method, including the TV+GTV method, is potentially missing small details in the reconstructed images due to the edge structures over-smoothed. By comparison, the proposed method (Figure 4(F)) has a good performance in terms of artifacts reduction and preservation of edge structures. Generally, in terms of artifacts suppression, the proposed method is much better than the OS-SART and FDK. As far as edges information protection is concerned, the proposed method is better than the TV and TV+GTV.
To visualize the difference in detail, Figure 5 presents the voxel curves along the coronal axis (253th row, 256th slice). The images obtained by use of the proposed method are reasonably accurate with only small distortions. Our proposed method curves are the most closet to the phantom curves, which means more edge structures information preserved.    For a clearer detail comparison, the edge structures (indicated by the red rectangular) area is zoomed in and shown in Figure 6. As for comparison between the OS-SART method, the FDK method, the TV method, the TV+GTV method, and the proposed method, the spinal bone details of low-contrast indicated by the red arrow in  Table 6, consistent with the visual quality, the proposed method achieves the best scores.  Another reconstructed slice (73th slice) of the XCAT phantom obtained by different reconstruction algorithms is demonstrated in Figure 7. OS-SART slice ( Figure   7(B)) and FDK slice (Figure 7(C)) both introduces serious artifacts. TV slice ( Figure   7(D)) and TV+GTV slice (Figure 7(E)) have a non-uniform intensity distribution in the edge structures, therefore, some details are over-smoothed. But the proposed method ( Figure 7(F)) performs better in the preservation of edges structure. In order to indicate the difference between voxels intuitively, Figure 8 indicates the voxel curves along the coronal axis at 250th row, 73th slice. The proposed method curves are the most consistent with the phantom curves, that means more edge structures information achieved. The results of the quality evaluation are shown in the Table 7, the proposed method achieves the best scores.   For a visualized detail comparison, the edge structures in the yellow rectangular area are zoomed in and shown in Figure 9. As indicated by the green arrow, OS-SART edges (see Figure 9(B)) and FDK edges (see Figure 9(C)) are attacked by the noise effects. TV edges (see Figure 9(D)) and TV+GTV edges (see Figure 9(E)) are slightly over-smoothed. As shown in Figure 9(F), the proposed method edges are cleaner and sharpener than the others. The results of the quality evaluation are shown in the Table   8, the proposed method achieves the best scores.  The edge structures in the green rectangular area are zoomed in and shown in As indicated by the yellow arrow, OS-SART edges (see Figure 10(B)) and FDK edges (see Figure 10(C)) are attacked by the noise effects. TV edges (see Figure 10(D)) and TV+GTV edges (see Figure 10(E)) are slightly over-smoothed. As shown in Figure   10(F), the proposed method edges are cleaner and sharpener than the others. The results of the quality evaluation are shown in the Table 9. Consistent with the visual quality, the proposed method achieves the best scores.

Discussion
CBCT imaging dose can be reduced by both a small number of projections and a low tube current (mAs) per projection. However, the conventional reconstruction methods result in streaking artifacts and obvious image noise in these two imaging circumstances, which can be clearly observed in the reconstructed results.
The CS technique has been utilized to reconstruct CBCT images with limited projection from under-sampled measurements. But, the condition of limited projection is an ill-posed problem. Since the reconstructed image itself doesn't have sparse property, TV transform has been widely adopted in CBCT reconstruction. The iterative method with TV norm as a regularization function demonstrated its effective in CBCT reconstruction with under-sampled projections. This TV minimization process, which penalizes the weight of each voxel at a same rate regardless of different spatial gradient, may not recover qualified CBCT images from ill-posed projection data and would result in edge structures over-smoothed.
In this work, we proposed a new strategy to deal with the deficits stated above by Although the proposed method achieves a better reconstruction quality than the others, it shows a slightly inferior noise suppression as compared to the conventional TV method, attributed to the GTV penalizes TV in the calculate of image total variation minimization process. Technically, compared with the other methods, the proposed method might increase the computational expenses due to: (1) additional iterates in the entire reconstruction with the concept of the GTV minimization, and (2) choosing optimal parameters. In fact, the extra iterates for the GTV process were confirmed to be not significant in our study since it converged much faster than the other methods.
One of the crucial challenges in optimizing the model of the proposed method is selecting several optimal regularization parameters (Equation 5) and (Equation 6). The selection of the parameters is very importance to the success of this proposed reconstruction method. and should to be carefully selected by considering the number of projections, and the characteristics of the anatomical features of the scanning object. In this study, we manually select and parameters to adjust the relative ratio between the data fidelity term and the regularization term to obtain the best CBCT image quality. It is found that the optimal value of these parameters is case dependent, the setting of or parameters depend on the experiences of simulation experiments. In future, we would study the parameter setting and try to find an automatic approach to guarantee a qualified selection of parameters. Note that GPU implementation and the parallel process of GTV minimization is currently an active research area. In next studies, we would work on optimizing our method to improve the speed further and implement GTV minimization using the GPU.
The value of our proposed method lies in the development of an algorithm that can effectively reconstruct CBCT images with the help of GTV penalizes TV. The proposed method could be used in various applications, such as 4D CT/CBCT reconstruction, where projections reduction in the scanning would be important. In the further we would also concentrate on the mathematics analysis of the proposed method.

Conclusions
In this study, we propose a CBCT image reconstruction method based on TV and GTV and redesigned ASD-POCS algorithm to minimize the objective function. The