Ring-Artifact Correction With Total-Variation Regularization for Material Images in Photon-Counting CT

We propose a ring-artifact correction method with a compressed sensing for material images obtained with a photon-counting computed tomography (CT) system. The ring-artifacts are caused by nonuniformity of detector properties. Conventional ring-artifact correction methods tend to degrade the quality of images. In contrast, compressed sensing methods can correct ring-artifacts with less degradation of the image quality owing to a priori knowledge that ring-artifacts appeared as stripes in sinograms. In this study, we extend the compressed sensing methods for material sinograms obtained with a photon-counting CT system. This is because material sinograms tend to be simpler and sparser, for which a compressed-sensing method can be more effective. We introduced a cost function with a total variation-regularization term and positivity constraint, and optimized it with a prime-dual splitting method. The feasibility of this method was confirmed by simulations and an experiment. In both the simulations and experiment, the proposed method better corrected the ring artifacts than those on attenuation domain and without a priori knowledge. The comparison with previous methods in literature also showed the same results. These results suggest that our method is effective for correcting ring-artifacts in material images.

reconstruction methods can be applied. On the other hand, the post-processing approaches operate in reconstructed images. Usually, they apply a polar-coordinate transformation and remove artifact components in the transformed images [18,19]. Although these methods are effective in ring-artifact removals, they lead to other problems such as blurring.
In contrast, a compressed sensing method can remove ringartifacts with less blurring [20][21][22][23]. A compressed sensing is a technique making use of a priori knowledge on the images. In medical images, it is a smoothness of the pixel values so that a total variation (TV) is usually applied in the regularization [24]. The TV regularization is often adopted to reduce noise and artifacts in CT images [25,26]. Recently, some studies focus on a machine-learning method to remove artifacts [27,28]. Nonetheless, a compressed-sensing based method has a clearer principle and still worth devising.
Previous studies adopted TV regularization on the reconstructed images, and removed ring-artifacts during the image reconstruction process. Paleo et al. defined a data-fidelity term in which the sinogram was split into two components: a "genuine" sinogram and spurious straight lines producing rings in a reconstructed image [21]. In addition, they adopted an L1 norm of the ring components and dictionary learning. They applied their method in the conventional CT images. Schmidt et al. attributed ring-artifacts to the non-uniformity of the gain in their photon-counting detector, and applied TV regularization [23]. They introduced a gain term in their objective function, and optimized it to remove the ring-artifacts. These studies showed that the ring-artifacts can clearly be removed with less degradation of the images.
In this study, we extend the above methods and apply TV regularization on material-domain sinograms. Since the material images and sinograms tend to be simpler than those of attenuation owing to the decomposition process, their gradients also tend to be sparser. Hence, the TV regularization could be more effective. In addition, this method is a kind of preprocessing approach, which can easily be combined with other pre-processing techniques. The feasibility of this method is investigated via simulations and an experiment. This paper is structured as follows. In section II, we define a cost function and explain an optimization method. In sections III and IV, we investigate the feasibility of our method by simulations and an experiment, respectively. We discuss the results in section V and conclude our paper in section VI.

II. THEORY
In this section, we describe our ring-artifact correction method. We define a cost function in section II-A and describe an optimization method in section II-B.

A. Definition of a cost function
Sinograms without artifacts should be smooth, and so we can assume that a total variation is small. Hence, an ideal sinogram reduces the following cost function: The first term is a data-fidelity term, where y and x are observed and material-domain sinograms, respectively. A is a transform matrix. While many studies adopt A as a system matrix, ours corresponds to attenuation coefficient. The subscripts k and m indicate energy-bin and material indices, and i and j are the projection-view and detector-position indices, respectively. The second term indicates TV regularization. As indicated in the last term of Eq.1, the TV term is an L1 norm of the L2 norm of the differences among adjacent pixel values, so that this term is a mixed norm of L1 and L2 (hereafter L 1,2 norm). The parameter λ regularizes weights of data-fidelity and TV terms in the cost function. A ring-artifact results from non-uniformity of the detector properties such as sensitivity and energy thresholds. We assume that these properties are constant during the data acquisition, and the ring component is independent of a projection view. Hence, we adopt the following equation: where x m,0 and z m are an initial estimated density sinogram and its stripe component of mth material, respectively. The initial density can be estimated with a singular value decomposition (SVD) method. Since the ring components can be assumed to be constant during the data acquisition, z m does not have a projection dimension subscribed by i. In addition, we applied an L1 norm of the ring component ( z L1 ) following Paleo et al. [21]. This is effective to avoid an over-smoothing due to the TV regularization. We also adopted a positivity constraint. Then, the ring-artifact removal method becomes the following problem: where λ 1 and λ 2 are regularization parameters of a datafidelity term and sparsity of the ring component, respectively. We note that x 0 and z do not have the same dimensionality as indicated in Eq. 2: x 0 + z means that the ring-components z are added each row of the sinogram x 0 as suggested by Paleo et al. [21].

B. Minimizing the cost function
In this study, we applied a prime-dual splitting method [29] to solve the problem in Eq. (3). Compared with a similar method like Alternating Direction Method of Multipliers (ADMM, [30]), the convergence is slow but requires less memory, because it does not need to calculate an inverse matrix. The prime-dual splitting method updates sinograms in the following way: where X is the dual variable of x, prox a proximal operator, f the data-fidelity term, and L the TV transformation matrix.
x is a function of z as indicated in Eq.2. The subscript g indicates the L1 norm of the ring-component, a projectionindependency constraint and a positivity constraint. Also, h indicates the L 1,2 norm for the TV regularization while h * is the convex conjugate of h. The proximal parameters γ 1 and γ 2 are related to only the update rates. As can be seen, the primedual method repeats calculations with proximal operators to solve the optimization problems with various constraints.
For convergence of Eq. 4, the followings should hold, where β is a Lipschitz constant and Eq. 5 indicates the Lipschitz continuous. The convergence speed of Eq. 4 is basically faster with larger γ 1 unless it is greater than the above limit.

III. SIMULATION
We investigated the feasibility of our method with simulations. We projected phantoms onto a detector and adopted our ring-artifact correction method on the sinograms. In the following subsections, we describe our simulations and present the results.

A. Phantoms
We used circular and CT-image phantoms, as shown in Fig.1 (b) and (c). The circular phantom is simple and regarded as a contrast enhanced tumor in a brain. The CT image phantom is useful for evaluating the ring-removal capability for a realistic complicated phantom.
In the circular phantom, we assumed a gold solution region surrounded by a water region. The diameters of the water (gray) and gold-solution (white) regions were 2.0 and 0.8 cm, respectively. The gold-solution region was located at a position 0.4 cm shifted from the center. The concentration of the gold solution was 3.0wt% (3.09 g/cm 3 ).
For the CT-image phantom, we retrieved a pelvis image from the Cancer Image Archive 1 . The image was divided into soft-tissue and bone regions. We approximated the energy dependence of the attenuation coefficient for these regions as those with H 2 O and Ca, and calculated the corresponding densities. The horizontal size of the image was about 50 cm.

B. Data acquisition and image reconstruction
We performed simulations with a two-dimensional fan-beam geometry. The simulation geometry is shown in Fig. 1 (a). For the circular (CT-image) phantom, the distance between the Xray tube and target was 72 (120) cm, and the target-detector distance was 8 (40) cm. The detector length was 4.0 (80) cm and its pixel size was 0.2 mm × 0.2 mm (2.0 mm × 2.0 mm). We set an X-ray tube voltage to 100 kV, and used three energy windows: 60-70, 70-80, and 80-90 keV. The geometry and detector specifications were similar to those of the experiment in section IV.
We assumed that the detector has 0.1, 1.0, 5.0 and 10.0% non-uniformity of sensitivity in each energy bin. The large non-uniformity reproduced strong artifacts. Instead, we assumed no uncertainty in the energy calibration for the sake of simplicity. In the following, the non-uniformities are temporarily fixed to 0.1% and 1.0% for the circular and CT-image phantoms, respectively. The results with other non-uniformity levels are discussed in section V.
CT-image phantom with 1.0% non-uniformity  Figure 2 shows example sinograms of the Ca density of the CT-image phantom. The densities were estimated by material decomposition with an SVD method. The 1.0% nonuniformity was amplified during the material decomposition. As can be seen, the sinogram (b) shows significant stripe artifacts.
We reconstructed the density sinograms with a filtered backprojection (FBP) method. Their matrix sizes are 200×200 and 256×256 for circular and CT-image phantoms, respectively. The example images are shown in Fig. 3 [(a), (g)] and Fig. 4 [(a), (g)]. They are the water and gold-solution regions in the circular phantom, and the water and calcium regions in the CTimage phantom, respectively. As can be seen, they have strong ring-artifacts due to the stripe artifacts in the sinograms.

C. Removing ring artifacts
To remove these ring-artifacts, we performed the method described in section II. We set the regularization parameter λ 1 =30 and 15 for the circular and CT-image phantoms, respectively. In general, lower λ 1 is optimized for stronger artifacts and simpler phantoms. The second parameter λ 2 was fixed to (a) Before corr. shows the images without a ring-artifact correction. The second, third and fourth columns shows those with our proposed, correction in attenuation domain, and that without a z term, respectively. The fifth and sixth columns show the reconstructed images with the Eldib17 and Münch09 methods. The non-uniformity level is 0.1%.
No error Before corr. Proposed 0.02. The first parameter of the prime-dual splitting method, γ 1 , was initially set to half of the limit value (Eq. 6), and then gradually decreased in the iterations for fine tuning. The γ 2 was fixed to 1.0, and the number of iterations was fixed to 200. The dependence of these parameters was investigated in the section III-D. In the optimization, regions without any materials were masked out.  Table I, where reconstructed images without non-uniformity were regarded as the ground truth. The first and second rows indicate the results before and after the correction. They show that the RMSEs were significantly reduced due to the correction. These results indicate that our method is effective to correct the ring artifacts in both simple and complicated phantoms.
To investigate the effectiveness of our method (Eq. 3), we performed additional simulations. In our method, the ring correction is conducted after the material decomposition, and a ring component z is adopted. Hence, we firstly performed a ring correction before material decomposition, that is, a correction in attenuation sinograms (corr. in att.). In this simulation, A in the Eq.3 became a unit matrix, x 0 the attenuation sinograms. The parameters were λ 1 =10 and 1.0 for the circular and CT-image phantoms, respectively. Secondly, we performed a simulation without the ring component (corr. w/o z). In this simulation, the parameters were λ 1 =30 and 500 for the circular and CT-image phantoms, respectively. These parameters were determined to minimize the RMSEs for each simulation.
Third and fourth columns of Figs. 3 and 4 show the reconstructed images of these simulations. As can be seen, ring artifacts are significantly remained, especially in the CTimage phantoms [ Fig. 4 (c), (d)]. The third and fourth rows of Table I indicate the RMSEs of these images. They show that the RMSEs are basically larger than those with our proposed method (second rows). The only exception is the H 2 O region in the circular phantom. In this case, the ring-artifact correction in attenuation images show the lowest RMSE (0.011 g/cm 3 ; third row), although our proposed method show a similar value (0.013 g/cm 3 ). This could be because the material decomposition does not make the region sparse, and our proposed method could not demonstrate an advantage for this phantom. Further comparison and evaluation are performed in the discussion section.

D. Parameter dependence
The resultant images can strongly depend on parameters in the cost function. Here we investigate this dependency. The most important parameter is λ 1 , which controls the balance between the data-fidelity and TV regularization. Higher λ 1 makes a larger weight in the data-fidelity term. Higher λ 2 leads to sparser ring-components and suppresses an over-smoothing effect due to the TV regularization. The proximal parameters γ 1 and γ 2 relate to only the convergence speed and not to the image quality.
We performed our ring-artifact correction method with various λ 1 where λ 2 was fixed to 0.02. We evaluated the quality of reconstructed images with RMSEs. Figure 6 show the results for the circular phantom [(a), (b)] and CT-image phantom [(c), (d)]. In Fig.6 (a) and (b), lower λ 1 leads to larger RMSEs, indicating that larger weights in the TV term leads to too blurred images. This is also seen in Fig.6 (c) and (d) for the correction without a z term (red lines). This is because TV regularization leads to blurred images if we do not adopted a ring-component in the cost function. On the other hand, larger λ 1 also leads to higher RMSEs. This is because it reduces the effect of TV regularization.
These figures also show that the best parameters are slightly different in different material. Although optimizing the parameters for each material leads to better reconstructed images, it doubles the time. Considering the best parameter is similar in each material, we adopted the same parameters for each material. We confirmed that our results less depend on the parameters. Fig. 7: Geometry of the experiment. The target was a gold solution with 3.0wt% in a 2 cm-diameter vial. The distance between the X-ray tube and target on the rotation table was 72 cm, and the target-detector distance was 8 cm.

IV. EXPERIMENT
We also performed an experiment to confirm the effectiveness of our method. In the following subsections, we describe the experiment and data analysis, and present the results.

A. Experiment with a photon-counting CT system
The experiment was performed with a photon-counting detector developed by Ogawa et al. [1]. It is a CdTe detector with four energy windows, eighteen modules with 40×40 pixels. The pixel size is 0.2 mm × 0.2 mm and length of the module gap is 0.4 mm, corresponding to two pixels. The energy resolution (FWHM) is 4-5 keV at 60-120 keV range.
As a target, we used a gold solution (size 5-10 nm; Tanaka holdings 2 ) in a 2-cm diameter vial. The gold concentration was 3.0wt%. Gold has a large attenuation coefficient and is chemically stable, so that it is used as a contrast agent [31]. It has a K-edge absorption feature at 80.7 keV, which is used for an energy calibration of a photon-counting detector [32].
The conditions of the experiment were as the follows. We used an X-ray tube TRIX-150S (Toreck, Japan) with a tube voltage of 100 kV and a current of 1.2 mA. An aluminum filter with a thickness of 2 mm was used. The energy thresholds were set at 60, 70, 80, and 90 keV, and we used three bins between the thresholds. We did not use the highest energy bin (> 90 keV) to minimize the pulse-pile-up effects. The distance between the X-ray tube and target was 72 cm and that between the target and detector was 8 cm. The experimental set-up is shown in Fig. 7. The data acquisition time was 15 sec. in each projection angle and 180 projection views were obtained.

B. Data analysis
Data processing was performed in the following way. A flat-field correction was done with a PMMA image for a target region and with an air image for the other region. We masked pixels whose values deviated more than three times normalized median-absolute-deviation from the median values in each module image. Module gaps were interpolated with a cubic spline function.
We calculated median values along the height direction of the obtained images to produce a sinogram. This is because producing the three dimensional sinogram was difficult due to the presence of many bad pixels. We assumed a fan-beam projection geometry. We calculated the attenuation in each energy bin in each pixel, and performed material decomposition to estimate an initial material density x 0 . Material decomposition was performed with a pseudo inverse matrix of the attenuation coefficient derived from an SVD technique. In this calculation, attenuation coefficients in each energy bin were derived from weighted average with a 100 kV spectrum, where the energy resolution was also considered.
We performed the ring-artifact correction method described in section II. The regularization parameters were λ 1 =0.2 and λ 2 =0.02. The optimization method was the same as the simulation in the previous section. The resultant sinograms were reconstructed with an FBP method. We can see that the stripes disappeared after the correction for both H 2 O and Au sinograms. The remaining weak stripe pattern is caused by module-gap residuals. Since the residuals depend on projections, they are smoothed out in a reconstruction process. Figure 9 shows comparisons of the reconstructed images. In both H 2 O and Au images, ring-artifacts are significantly reduced after the correction [(b), (f)]. A region with larger values in H 2 O images corresponds to an edge of the vial. Since we adopted the cross-section of H 2 O rather than PMMA for simplicity, it is difficult to estimate the density in this region.

C. Results
Hence, we evaluated the reconstructed images in the inner regions, which are within 0.9 cm (45 pixels) from the center. In Table II, we show the averages and standard deviations (σ H2O , σ Au ) of the densities in this region. Without the ring-artifact correction, σ H2O and σ Au were 10.2 and 15.5%, respectively. In contrast, they were reduced to 8.2 and 11.2% in the corrected images. On the other hand, the average densities were changed only slightly during the correction, indicating that our correction removed only the ring components.

V. DISCUSSION
In the previous sections, we showed that our method significantly reduced the ring-artifacts and the density variances. In this section, we compare our method with others and discuss the advantages of our method.  Before corr. Eldib17 Münch09 Fig. 10: RMSEs against non-uniformity levels for each material image in each phantom for each ring-artifact correction method. Blue, green, and red solid lines indicate the proposed, attenuation domain, without z term methods, respectively. Black, magenta, and cyan broken lines indicate those before correction, with Eldib17 and Münch09, respectively.

A. Comparison with other methods
To evaluate the feasibility of our method, we performed ring-correction methods introduced by Münch et al. [15, here after Münch09] and Eldib et al. [16,here after Eldib17]. Both these methods are pre-processing and can be directly compared with our method. The Münch09 method utilizes a wavelet-Fourier filtering to remove stripe artifacts in sinograms. It performs a wavelet decomposition and the vertical component is blurred with a Gaussian damping function on the Fourier domain. In this study, we applied a Daubechies wavelet (db15). Standard deviations in the damping function were set to 1.0-10.0 to minimize the RMSEs, depending on the nonuniformity levels.
The Eldib17 method is divided into three steps. First, defective pixels are identified and corrected. Second, stripe artifacts are corrected. In this step, sinograms are averaged over the projection direction to produce the 1D vector, by which sinograms are divided by row by row over all view angles. Although this process corrects stripe artifacts caused by non-uniformity, it also leads to a contrast change. Third, the contrast is compensated. In this step, both original and contrast changed sinograms are reconstructed and the difference between them is blurred with a Gaussian filter. The blurred difference is added into the contrast changed reconstructed images. This method successfully corrects ring artifacts and less affects the contrast in the reconstructed images. Furthermore, this method is independent of image-reconstruction method like a pre-processing method.
In this work, we adopted their methods on both the simulations and experiment. We slightly modified the Eldib17 method for our study. We did not apply their first step because this can independently be adopted for any pre-processing method. In the second step, we did not take a division but a difference because our sinograms were material domain, which were derived from a logarithm of the X-ray intensity. Finally, we optimized the standard deviation of the Gaussian filter, which is 10-pixel in most cases according to Eldib et al.. We changed them to be 0.1∼10 and found that it is usually around one pixel to minimize the RMSE. In the following comparison, we adopted the optimized standard deviation for each case.
In Figs. 3 and 4, we compare reconstructed images of the circular and CT-image phantoms, respectively. The second, fifth and sixth rows show those with our proposed, Eldib17 and Münch09 methods. The third and the fourth show, again, results of our method on attenuation sinograms and without a z term, respectively. As can be seen, our proposed method most clearly corrects the ring artifacts, especially for the CTimage phantom. Quantitative comparison in Table I  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. the same results. The main principle of Eldib17 and Münch09 methods is to smooth out only stripe artifacts in the sinogram. However, strong artifacts can be remained after the smoothing. In contrast, our method directly compensates stripe artifacts and can be adopted even for the strong artifacts. We also compare our proposed method with the others for different non-uniformity levels. In Fig. 10, we show RM-SEs of reconstructed images against several non-uniformity levels for different ring-artifact-correction methods. We did not show the results of the circular phantom with 10% nonuniformity because we could not reconstruct the images due to large uncertainty. Also, those of the CT-image phantom with 0.1% non-uniformity were not shown because this low non-uniformity did not produce ring-artifacts. As expected, all methods show larger RMSEs for larger non-uniformity. The best method, which shows the minimum RMSE, depends on the phantom and non-uniformity level. Nonetheless, our proposed method (blue) basically shows the minimum RMSE.
Finally, we compare ours with Eldib17 and Münch09 methods using the reconstructed images from the experiment (Fig.  9). Similar to the simulation results, our method better corrects the ring artifacts. However, the RMSEs of our proposed, El-dib17 and Münch09 methods are different only slightly (Table   II). This is due to the original data have a large statistical noise which dominates the RMSEs. Thus, making use of the simplicity of our experimental phantom, we produced radial profiles of these images (Fig. 11). This is because a radial profile is less affected by statistical noise and emphasizes the ring artifacts. We do not show those with the Münch09 method because it is very similar to those with the Eldib17 method and showing multiple lines make the graph very complicated. Whereas the profiles without correction (gray) and with the Eldib17 method (black) show a significant fluctuation, those with our method (pink) are much flatter. These results indicate again that our method better corrects ring artifacts.

B. Other advantages of our method
In this study, we applied a simple cost function compared with previous studies. Paleo et al. used an L1 norm of the ring component and dictionary learning in order to avoid a Cartoon effect resulting from TV regularization [21]. On the other hand, we did not apply dictionary learning despite which no artifact due to an over-smoothing effect in the sinogram was observed. This could be due to the fact that the material sinograms tend to be simpler than attenuation ones. Hence, our results suggest that our cost function is appropriate for correcting ring-artifacts in material images.
We applied TV regularization on projection domain. However, this is usually applied in reconstructed images because they tend to be flatter than the sinograms. Although the TV regularization to non-flat images could lead to an oversmoothing effect, we did not observe such effects (e.g. Fig.  3). This result suggests that the TV regularization is also appropriate in sinograms.
Thus, different from the previous studies [19,21], we applied the ring-artifact correction on the projection domain. That is, our method is a pre-processing technique. Hence, our method can be easier to combine with other pre-processing techniques. For example, corrections of spectral distortion caused by pulse-pile-up and charge-sharing effects are easier to perform on the projection domain. This applicability is one advantage of our method.

VI. CONCLUSION
In this study, we proposed a ring-artifact correction method with a compressed sensing for photon-counting CT images. We used a priori knowledge that ring-artifacts appeared as stripes in sinograms, and produced a cost-function with TV regularization. The cost function was optimized with a prime-dual splitting method. To evaluate the feasibility of our method, we performed simulations and an experiment with a photon-counting detector. We used a circular phantom with water and gold-solution regions and CT-image phantom with water and calcium regions. Applying our method on a material domain, ring-artifacts were significantly reduced. We evaluated the reconstructed images, and confirmed that the RMSEs and density variances were significantly reduced. We also confirmed that both the correction on material domain and use of a ring-component z term are crucial for ring-artifact correction. These results indicate that the proposed method is appropriate for correcting ring-artifacts in photon-counting CT images.