Iterative image-domain ring artifact removal in cone-beam CT

Ring artifacts in cone beam computed tomography (CBCT) images are caused by pixel gain variations using flat-panel detectors, and may lead to structured non-uniformities and deterioration of image quality. The purpose of this study is to propose a method of general ring artifact removal in CBCT images. This method is based on the polar coordinate system, where the ring artifacts manifest as stripe artifacts. Using relative total variation, the CBCT images are first smoothed to generate template images with fewer image details and ring artifacts. By subtracting the template images from the CBCT images, residual images with image details and ring artifacts are generated. As the ring artifact manifests as a stripe artifact in a polar coordinate system, the artifact image can be extracted by mean value from the residual image; the image details are generated by subtracting the artifact image from the residual image. Finally, the image details are compensated to the template image to generate the corrected images. The proposed framework is iterated until the differences in the extracted ring artifacts are minimized. We use a 3D Shepp–Logan phantom, Catphan©504 phantom, uniform acrylic cylinder, and images from a head patient to evaluate the proposed method. In the experiments using simulated data, the spatial uniformity is increased by 1.68 times and the structural similarity index is increased from 87.12% to 95.50% using the proposed method. In the experiment using clinical data, our method shows high efficiency in ring artifact removal while preserving the image structure and detail. The iterative approach we propose for ring artifact removal in cone-beam CT is practical and attractive for CBCT guided radiation therapy.


Introduction
Flat-panel detectors can acquire 2D projections using x-rays (Cho et al 1995, Scarfe et al 2006, and have been widely used in cone-beam CT (CBCT) systems. In addition, flat-panel detectors can be used in many fields such as biomedical imaging, image-guided radiotherapy (IGRT) (Jaffray et al 2002, Yan et al 2014, Liang et al 2016, industrial x-ray flaw detection (Hofmann et al 2004, Zhang et al 2005, and materials research. However, some ring artifacts are always located at the center of the rotation axis in the reconstructed images, mainly due to the pixel gain variations in the detector (Vidal et al 2005). They appear in the form of concentric circles in the reconstructed image. The ring artifacts severely impair the quantification of the diagnosis. In addition, these artifacts are considered an integral part of the image (e.g. tumor) by the computer, which may lead to a wrong diagnosis (Davis andElliott 1997, Sadi et al 2010). In radiotherapy, severe ring artifacts may be identified as bone, resulting in incorrect positioning. Since ring artifacts make post-processing and quantitative evaluation complicated, the development of an effective and convincing method of ring artifact removal is crucial.
In recent years, many researchers have been developing CBCT methods for ring artifact removal (Raven 1998, Tang et al 1999, Boin and Haibel 2006, Ketcham 2006, Sadi et al 2010, Wu et al 2015a. These methods can be roughly divided into two categories: pre-processing in the projection domain and post-processing in the image domain. One of these pre-processing methods is a sinogram processing method (Abu Anas et al 2010, Ashrafuzzaman et al 2011, Guo et al 2015, in which the main principle is that the ring artifacts in projections are manifestd as stripe artifacts in the vertical direction. Compared to ring artifacts, stripe artifacts are relatively easy to detect and eliminate. As early as 1978 (Kowalski 1978), pre-processing in the projection domain has been proposed. A simple low-pass filtering was proposed to eliminate inconsistent information in projections caused by ring artifacts. However, this disturbs the high-frequency information in the projections and deteriorates the quality of the reconstructed images. In recent years, many methods have been proposed to overcome this deficiency. In 1998, Raven (1998) eliminated ring artifacts in sinograms because the artifacts manifest as stripe artifacts in the parallel vertical direction and there is a severe change of the gray value in the horizontal direction. As a result, stripe artifacts in the vertical direction are located in the center of image after Fourier transform. Thus, low-pass filters are applied in the horizontal direction to remove the artifacts while preserving high-frequency image details. This approach was further developed in 2009 (Munch et al 2009). Combining wavelet decomposition and a Fourier low-pass filter, a more rigorous distinction of artifacts and image details was achieved. More recently, Altunbasa et al (2014) have removed ring artifacts by identifying clusters of isolated pixels; however, this method of pixel gain correction needs a unique data set corresponding to each tube current setting, which makes the correction process incompatible.
Pre-processing in the projection domain needs to access the raw projection data, which restricts the applicability of the algorithm. Thus, a post-processing method which was independent of the raw data was proposed in order to significantly extend the range of application. A most popular post-processing method for removing ring artifacts is transforming the reconstructed image from a Cartesian to polar coordinate system (Chen et al 2010, Wei et al 2013. Here, the ring artifacts are manifest as stripe artifacts in polar coordinate systems and, after image processing in polar coordinates, the image is transformed back into a Cartesian coordinate system. In 2004, Jan Sijbers et al (Sijbers and Postnovz 2004) proposed a ring artifact correction method based on morphological operators. First, a region of interest (ROI) to be processed is selected. Then a sliding window is applied to generate the artifact template and acquire the corrected image by subtraction. The parameter selection is crucial to the accuracy of the artifact correction in this method. Subsequently, in 2009, Chen et al (2009 presented an effective method using independent component analysis (ICA), but this process deteriorates the image details. Kyriakou et al (2009) andPrell et al (2009) have used a median filter to remove stripe artifacts in polar coordinates with the assumption that the gray value feature in the ring artifacts are local extrema. Recently, Yan et al (2016) proposed a method of variationbased ring artifact removal with sparse constraint; this method showed effective correction, but the process of coordinate transformation deteriorated the spatial resolution of the images.
In this paper, a general iterative image-domain ring artifact removal method is proposed in CBCT to tackle the disadvantage of the aforementioned algorithm. The proposed algorithm removes the ring artifacts from pixel gain variations using flat-panel detectors. The relative total variance (RTV) is applied in a polar coordinate system. To improve performance, an iterative framework is applied in the processing of ring artifact removal (Wu et al 2015b). The framework of the ring artifact removal can be widely used in clincal CT images without deteriorating the image details.
Both numerical and real data are applied in the evaluation of the proposed method. In the numerical simulations, a 3D Shepp-Logan phantom is used for quantitative evaluation of image quality. With the real data, a uniform acrylic cylinder and Catphan©504 phantom are used with the Varian onboard imaging (OBI) system, whose detector has pixel gain variations. Finally, we use the clinical data of a head patient, in which ring artifact removal is more challenging.

Workflow
As is known, ring artifacts come from pixel gain variation during data acquisition, which deteriorates the structure and texture of images. In order to isolate the ring artifacts from the texture, we transform the raw CBCT images into polar coordinate images. After conversion, the ring artifacts manifest as stripe artifacts and an edge-preserving smoothing process is applied to acquire an ideal template image without texture, so that we can extract the stripe artifacts from the texture while maintaining structure and detail. Figure 1 shows the framework of iterative image-domain ring artifact removal. To start the framework, we input the reconstructed CBCT images with ring artifacts, which are transformed into the polar coordinate system. Using the RTV algorithm, the template image is generated; then the image is subtracted from the CBCT image and the residual image, including image details and ring artifacts, is generated.
Since the ring artifacts manifest as stripe artifacts in the polar coordinate system, they can be extracted by mean value from the residual image. The image details are generated by subtracting the artifact template from residual images. The dashed rectangle in figure 1 shows how the image details are separated from the residual image. Then we can compensate the image detail to the template image and update the uncorrected images. To reduce the error in the generation of the template image, an iteration process is implemented until the difference in the extracted ring artifact image is minimized. Finally, the extracted ring artifacts are transformed into Cartesian coordinates and subtracted from the raw CBCT images to obtain corrected images without loss of spatial resolution.
The following sections will list the detailed operations, including structure extraction, ring artifact extraction, and stopping criteria as shown in figure 1.

Structure extraction using RTV
The uncorrected CBCT images are transformed into the polar coordinate system using spline interpolation so that the ring artifacts can be manifested as stripe artifacts . Then a template image is generated using a smoothing technique to remove ring artifacts and image details. Since ring artifacts and image details are considered to be texture in the image, a smoothing algorithm with good edge preserving is crucial to removing the texture while preserving the structure of the CBCT images. In this work, we use an RTV algorithm as a smoothing algorithm with edge preserving. The objective function (Xu et al 2012) is expressed as (1) where I p is the input CBCT image and p indexes pixels in images. S is the template image without ring artifacts and image detail. The data term (S p − I p ) 2 is the smoothing item. Due to the complexity of human structure, two pixel-wise windowed TV q∈R(p) g p,q · (∂ x S) q o and q∈R(p) g p,q · (∂ y S) q are applied to the count of the absolute spatial difference within the window R(p) in the x and y directions, respectively. Pixel q belongs to R(p), where pixel p is located at the center of the window. ∂ x and ∂ y are the partial derivatives in the x and y directions. g p,q is a weighting function. Two windowed inherent variations q∈R(p) g p,q · (∂ x S) q and q∈R(p) g p,q · (∂ y S) q capturing the overall spatial variation help to distinguish prominent structures (e.g. bony structures) from the image detail. λ is a weight which is crucial to the intensity of smoothing. In order to avoid division by zero, ε is set to a small positive number. λ (•) is called RTV, which can separate the image detail and ring artifact due to inherent variation and RTV measures. The objective function of template images in equation (1) is nonconvex and can be solved with the penalty of a quadratic measure. The optimization process of the objective function was proposed by Xu et al (2012) and is shown in the appendix.

Extraction of ring artifacts
By subtracting the smoothing images from raw CBCT images, residual images are generated which contain ring artifacts and image details. The ring artifacts manifest as stripe artifacts in the polar coordinate system, and thus the stripe artifacts in the θ direction have the same gray value.   Based on the fact that stripe artifacts voxels have the same gray values in the θ direction caused by the same pixels of detectors in polar coordinates, the mean value calculation in the θ direction can be estimated as the artifact shown in figure 2. Since the detector's dead pixels have higher and lower gray values in the image domain, we take the range of 80% in the middle size of the gray value points for mean value calculation. This can rule out the effect of dead pixels and results in a more accurate extraction. A better template image of stripe artifacts can be efficiently generated, as shown in figure 1. Subsequently, a detailed image without the stripe artifacts is obtained by subtracting the stripe artifact template from the residual image. Finally, the image details are compensated to the smoothing image for stripe artifact removal while maintaining image details.

Stopping criteria
In order to handle structure extraction error when using RTV with serious ring artifacts, the RTV process is updated in every iteration. If the relative difference between two sequential extracted ring artifacts reaches a pre-set value, the iteration is terminated. We can quantify this difference using L2 norm, and the stopping criteria can be defined as where r k+1 and r k represent the extracted ring artifacts of the current and previous iterations, respectively. In this study, s d was set to 0.002.

Evaluation
First, the proposed method was applied to simulated data of a 3D Shepp-Logan phantom to validate the performance of this study. The source trajectory had a radius of R = 1000.0 mm and the distance from source to detector was S = 1500.0 mm. The projection consisted of 1024 × 768 elements and its total area was 400 × 300 mm 2 , exactly matching that of the CBCT system on the Varian Trilogy. We acquired projection data of the 3D Shepp-Logan phantom at 660 views uniformly distributed over 2π. We generated the ring artifacts in projection space. Since the ring artifacts are caused by pixel gain variations using flat-panel detectors, we can amplify some of the detector elements randomly, which shows as parallel noise-free stripes in the sinograms. Images with ring artifacts were then reconstructed using the sinograms with pixel gain variations. The simulation of the uncorrected image is shown in figure 3(a). A reconstructed image using Feldkamp-Davis-Kress (FDK) algorithm without any artifacts was generated for comparison. After achieving good image performance in the simulated data of the 3D Shepp-Logan phantom, the proposed method was further evaluated using real data, with a uniform acrylic cylinder and Catphan©504 on the Varian Trilogy OBI system. Side-by-side image compariso ns were carried out between the uncorrected and corrected images.
Furthermore, a more challenging evaluation using clinical images of a head patient was carried out due to the complex structure of the patient image. The patient image was acquired  on a Varian Trilogy CBCT system. For comparison, a diagnostic CT image in the same patient was acquired to register to the CBCT image. Table 1 lists the scanning parameters for the 3D Shepp-Logan phantom, uniform acrylic cylinder phantom, Catphan©504 phantom, and patient head. In studies presented in this paper, we implemented the entire data processing procedure in Matlab, including RTV, coordinate transform, and FDK reconstruction. The evaluation of image detail maintenance is calculated using the structural similarity index (SSIM), which can be defined as where c 1 = (K 1 L) 2 and c 2 = (K 2 L) 2 . L is the dynamic range of the pixel values. Two small constants c 1 and c 2 are calculated by setting K 1 = 0.01 and K 2 = 0.03, respectively. Note that, if the two input images are the same, the value of SSIM is one. Furthermore, we quantify the global quality of the image in the selected uniform ROI using spatial uniformity (SU), which can be defined as where I tissue is the selected uniform ROI, num (•) represents the pixel number in ROI, and hist is the histogram of the images. SU = 1 means a completely uniform ROI while SU between 0 and 1 means a variation of CT number in ROI.   Figure 5 shows the reconstructed CBCT images of a uniform acrylic cylinder phantom. Ring artifacts are visible in the axial slice without ring artifact removal in figure 5(a). Using the proposed method, the intensity of the ring was   Figure 7 shows the reconstructed CBCT images using Cat-phan©504, where the first row shows the uncorrected images, the second row shows the corrected image, and the third row shows the extracted ring artifacts using the proposed method. The three columns represent images of different slices of the phantom. The ring artifact is obvious in the first row, while no obvious ring artifact can be detected in the second row. The coronal views of the reconstructed images are shown in figure 8, where the ring artifacts manifest as stripe artifacts. The ring artifacts are clearly shown, as indicated by the white arrow in figure 8(a), while no ring artifacts can be observed in the center of figure 8(b).

Patient head studies.
For greater challenge, a patient head image obtained by clinical CBCT was evaluated using the proposed framework. Four rows of axial images of the patient's head are shown in figure 9. Column (a) shows the uncorrected images and column (b) the corrected images. Column (c) shows the registered MDCT image as a ground truth for comparison. The two small images to the left are enlarged regions indicated by the red square. As shown in figure 9(b), the rings were suppressed significantly after applying the proposed method, and the corrected images match well with the reference image. The performance is further evaluated using SSIM in the red square of figure 9(a1). Compared with the MDCT image in column (c1), the SSIM is increased from 70.23% to 85.34%.

Error in RTV smoothing
Accurate structure extraction is not easily performed due to inevitable smoothing errors in the RTV smoothing image, especially for CBCT images with serious ring artifacts. A worstcase scenario is applied to validate the performance of the proposed framework, where two inaccurate RTV smoothing images are generated as shown in figure 10. Rows 1 and row 2 are the head patient CBCT images at λ = 0.008 and λ = 0.05, respectively. The ring artifacts are removed, indicating that the proposed framework can tolerate the RTV smoothing error. This suggests that the framework is robust and the result of the corrected images is not overdependent on the parameter setting.

Iterative scheme and stopping criteria
As image details are difficult to generate on the uncorrected image, an iterative scheme can improve the effectiveness of the ring artifact removal, which is important to the improvement of smoothing as well as the separation of texture and structure. The iteration efficiency is shown in figure 11, where each column shows the images generated in the successive iteration steps. Row (a) in figure 11 shows the RTV smoothing images; row (b) shows the residual images; row (c) shows the extracted ring artifacts; and row (d) shows the images with corresponding correction. It is clear that the smoothing using RTV becomes increasingly accurate after the iterations-especially in the region with serious ring artifacts, as shown by the dashed box in the region of figure 11(a). The proposed scheme corrects the ring artifacts effectively, while maintaining structure and detail.
If the algorithm reaches the pre-set stopping criteria defined in equation (2), the iterative scheme will stop. The convergence is stable and fast, as shown in figure 12. These four evaluations with stable convergence suggest that the framework of this method is robust.
In our experiments, we set the stopping criteria at 0.005, and the whole procedure stops after 10, 9, 9, and 11 iterations for the four data sets (Shepp-Logan phantom, acrylic cylindrical phantom, Catphan©504, and head patient), respectively.

Advantages of the study
Unlike the previous study, this study uses an RTV algorithm to separate 'structure + texture' images, so that the image detail can be well preserved. The proposed method shows good performance in a clinical CBCT system. The advantages of the proposed method are listed as follows. First of all, the proposed algorithm can preserve image details in CBCT images, which avoids the possibility of bias towards the template images of the corrected CBCT image. Furthermore, the post-processing method can remove the ring artifacts without the raw projections. The process is implemented in the reconstructed images and it is more compatible with the IGRT, without changing the software and hardware configuration. Finally, the algorithm has low time complexity, which can be accelerated by a graphics processing unit. We will implement this acceleration in the near future.

Conclusions
In this study, we propose a robust iterative scheme for ring artifact removal. The corrected images are acquired using RTV without losing the spatial resolution of the images. Our method is evaluated on a 3D Shepp-Logan phantom, Catphan©504 phantom, uniform acrylic cylinder, and images from a head patient. In the simulated data, the proposed method increases the SU by a factor of 1.68 and SSIM from 87.12% to 95.50%. The ring artifacts are suppressed in the CBCT image, while the image structure and detail are maintained. The proposed method is practical and attractive for CBCT-guided radiation therapy.
The expression (A.2) indicates that the u x of each pixel combines the information of the adjacent gradient in an isotropic spatial filter manner. With the standard deviation σ, G σ is a Gaussian filter. w x and the pixel-wise gradient are related. The element-wise gradient is shown in the division in (A.2) and * is the operator of convolution. The y-directional penalty can be expressed in a similar manner to the x-direction. After this operation, equation (1) is written in matrix form: The vector representation of S and I are denoted by v s and v I , respectively. C x and C y , from the operators of discrete gradient, are Toeplitz matrices. A special iterative optimization process can be applied to (A.4). We can obtain a numerically stable approximation naturally because of the quadratic parts and nonlinear decomposition. The approximation, which is stable, is found to be very effective for estimating the texture and structure images. The optimization steps are as follows: 1. From the estimated structure images S in the early iteration, the values of u and w are acquired according to equations (A.2) and (A.3). They can construct the matrices of U and W . 2. Minimization using the values of U x , W x , C x , and W y comes down to obtaining the solution of linear system in the iteration; this can be defined as where L t = C T x U t x W t x C x + C T y U t y W t y C y is the weight matrix and 1 is a unit matrix.
(1 + λL t ) is the symmetric positive definite Laplace matrix. Discrete gradients were approximated from the forward difference. This leads to a sparse five-point Laplace matrix. An effective solver can be used for this.