Ring Artifacts Reduction in CBCT:Pixels Detection and Patch Based Correction

Haitong Zhao,Yi Li and Shouhua Luo,＊

1School of Biological Science and Medical Engineering,Southeast University,Nanjing,China.

Abstract:The ring artifacts introduced by the defective pixels with non-linear responses in the high-resolution detector,have a great impact on subsequent processing and quantitative analysis of the reconstructed images.In this paper,a multistep method is proposed to suppress the ring artifacts of micro CT images,which firstly locates the positions of the defective pixels in the sinogram,and then corrects the corresponding value in the projections.Since the defective pixels always appear as vertical stripes in the sinogram,a horizontal curve is derived by summing the pixel values along vertical direction,thus the abrupt segments related to the defective stripes are enhanced notably,and a proportion coefficient based on the second derivative of the curve is taken as the indicator for the position and the severity of the defective pixels.Then,the detected defective pixels in the sinogram are transferred and relocated in the projections,an improved 3D block matching filtering (BM3D) algorithm is applied to restore the defective pixels in corresponding projection images.In the end,the tomographic images are reconstructed from the corrected projections.In the experiment,a small piece of the motherwort’s rhizome and a part of a mouse’s lung are imaged by micro-CT,and the result shows that,compared with the other four state-of-art methods,the proposed method has a great reduction on the ring artifacts of the reconstructed images,and makes less impact in spatial resolution and contrast in the same time.

Keywords:Ring artifacts; defective pixels; micro CT; BM3D; image inpainting

1 Introduction

Micro-CT is a non-destructive 3D imaging technique for small objects,which plays an important role in basic medical research,such as bone structure observation [1],advanced material exploring etc.The defective pixels are the intrinsic characteristics of the detector [2],which always result in some ring artifacts in the reconstructed images.Existing ring artifacts reduction methods generally fall into two categories:hardware based and software based.

Hardware based methods,which always shift the positions of the detector or object movement and calculate the average response of a pixel in different positions [3-5],require special hardware arrangement with high precision.Altunbas et al.[6] use a Cu beam filter in flat field correction,which needs acquisition twice and fails to remove the rings completely if the response function of different detector element is not same.

In contrast with the hardware based ones,the software based approaches are easier in practice,which often transform the reconstructed images into polar coordinates to remove ring artifacts [7-8].However,it is hard for keeping same resolution due to interpolation during the transformation via different coordinate systems.Xiao et al.proposed a method [9] for extracting ring artifacts in polar coordinates,though the image resolution is kept well,but it may not be suitable for the complicated scenario as it supposes that the ring artifacts in the θ direction have the same gray value.The dictionary representation [10] and the deep learning based methods [11] are also used in reconstructed images.These methods achieve better results than before,though the problem of image blurring still exist.Removing ring artifacts is more feasible in sinogram domain,since the defective pixels always appears as single or multiple stripe artifacts.Various filters,such as center-weighted median filter [12],wavelet-Fourier filter [13] and Gaussian filter [14] have been utilized in removing ring artifacts.Podgorsak et al.proposed a median filtering-based method [19]based on the assumption that the values measured by the same defective pixel have same attenuation against the ideal values.These methods have achieved some effects in ring artifacts removal,but some detailed features of the images may be lost during the filtering process [15].A multistep method based on the detection and correction of the defective pixels can remove ring artifacts more effectively without imparting noticeable distortion in the image [16-17].However,these methods tend to introduce radioactive artifacts with wrong correction of defective pixels.

In views of these problems,we propose a novel approach,which utilizing a derivative based algorithm firstly to locate the position of defective pixels in sinogram,then applying an improved 3D block matching filtering (BM3D) [18] algorithm to restore the defective pixels in corresponding projection domain.The block-matching method is employed to search for the patches similar to the ones with the defective pixels.The 3D hard-threshold filter and wiener filter act on the 3D patches to estimate the value of defective pixels,of which the parameters depend on the severity of the defective pixels.The specific procedures including algorithm details,data experiments will be elucidated in the following chapters.

2 Method

The proposed method here is divided into two steps:defective pixel detection step and corresponding correction step.

2.1 Defective Pixel Detection

The sinogram is shown in Fig.1(a).Sub-sinograms [16] with a width of m and a height equals to the height of the sinogram are cropped from the sinogram.Figs.1(b) and 1(c) are both the sub-sinogram with m = 9.As are shown in Fig.1(b) and Fig.1(c) respectively,the responses of each defective pixel in every θ are much lower or higher than those of the neighboring uncorrupted ones.

Figure1:(a) The sinogram of shepp-logan phantom.the red arrows indicate the stripe artifacts.(b) the sub-sinogram ① is zoomed in.(c) the sub-sinogram ② is zoomed in

Suppose that the actual sinogram is P(i,θ),i is the one of pixels in a row of detector and θ is the one of views.The sum of all the values in each column can magnify the difference between the defective pixels and the normal ones.The summation is as follows:

S(i) is the sum of all the gray values of each column.The curve of S(i) is shown in the second row in Fig.2.If there has a defective pixel in the sub-sinogram,a peak or valley will appear in S(i).In order to highlight the stripe artifacts,the second derivative D(i) of S(i) is taken as an indicator to distinguish the defective from the normal pixels,as shown in Eq.(2):

The D(i) is enough for a single pixel stripe identification.There may exist some continuous defective pixels in the detector,thus the width radius r is introduced to detect the stripe with a width of 2∗r or 2∗r+1,and the second derivative D(ic) at the center of sub-sinogram iccan be rewritten as Dr1(ic) or

Dr2(ic) respectively in Eq.(3) and Eq.(4):

When D(ic) and D(ic±r) are of opposite sign and the amplitude |D(ic)| is greater than |D(ic±r)|,the pixel icmay be the center of a defective band with the width of 2∗r or 2∗r+1.Fig.2 shows the S(i)and D(i) curve of stripes with different widths.

Figure2:Characteristics of defective pixels with different widths.The second row is S(i),and the third row is D(i).(a-c) Normal pixel feature.(d-f) Feature of single defective pixel.(g-i) Feature of two consecutive defective pixels.(j-l) Feature of three consecutive defective pixels

In order to determine the defective pixels,a proportion coefficient c is denoted as the severity level in Eq.(5):

Thus,the pixel i is as a defective pixel when c＞c0,where c0is a fixed threshold.In the experiment of this paper,c0is 1000.

2.2 Defective Pixel Correction

As there is more correlative geometrical information among the defective pixels in the projection images than that in sinogram,it is feasible to recover the value of defective pixels while mapping them in the projection domain.The specific procedures are presented based on the BM3D algorithm [18],which consist of two steps:basic estimate and final estimate,as shown in Fig.3.

Figure3:Flowchart of the projection image inpainting algorithm

Supposing the projection image is P,the location and the severity information of the defective pixels are set A ({(xi,yi)∈A | 0＜i≤n}) and set C({(ci)∈C | 0＜i≤n}),respectively.For every defective pixel Aiin A,we construct a reference block whose center is Ai.For each such block:

· Find blocks in Ns×Nsneighborhood that are similar to the reference one (block-matching,as shown in Fig.4) and stack them together to form a 3D array (group).

· Perform collaborative filtering of the group and return the estimate of Ai.

Euclidean distance is used to evaluate the similarity between two blocks.That is only the distance is smaller than a fixed threshold τmatcℎ,it can be considered similar to the reference block,and be grouped together.The block size is (N1×N1) and the number of grouped blocks is restricted by setting an upper bound N2.The distance could be calculated as:

where ‖·‖2denotes the l2-norm and PAiand Pxare represent two blocks.One constrain is introduced to guarantee the accuracy of similar blocks we found:If a block contains defective pixels,it cannot be a candidate of group.

Figure4:A simple example of grouping in a simulated image,where for each reference block (with thick borders) there exist perfectly similar ones (with dotting borders)

In the basic estimate,block-matching is carried out for each point in A to obtain the group S,and the value of the defective pixel is estimated by hard threshold filtering:

where c∈C is the proportional coefficient calculated from previous,representing the extent of a defective pixels,and λ3Dis a constant.The basic estimate can be expressed as:

The value of the defective pixels is estimated by the hard threshold,and the other points remain unchanged.

In the final estimate,we define the empirical Wiener shrinkage coefficients from the energy of the 3D transform coefficients of the basic estimate as:

The final estimate of the whole projection can be obtained as follows:

The final estimate is the restored projection image which can be directly used for CT reconstruction.

3 Results

To verify the effect of our algorithm,a Catphan Phantom,a mouse and small piece of the motherwort’s rhizome were scanned by the high resolution micro CT in our laboratory.After acquiring the high-quality projection data,the value of some random selected pixels is increased or decreased proportionally to simulate non-linear responses.

The parameter in these experiments is presented in Tab.1.The effect of our method is illustrated in Fig.5.Among them,the first column is the original images,the second column is the corrected images,the third column is the residual images,and the fourth column is the zoomed image details.

Table1:Parameter sets for experiment

Figure5:Comparison between the original images and the results of ring correction.(a) cross section of the phantom.(b) lung consolidation of the mouse.(c) the normal lung of the mouse.(d) the rhizome of the motherwort

Figure6:The ring removal results of the motherwort’s rhizome made by the different methods:(a) the original image.(b) the polar coordinate based method.(c) the wavelet-Fourier method.(d) the GAN-based method.(e) the median filtering-based method.(f) the polynomial interpolation method.(g) the proposed method

We compare the proposed algorithm with several state-of-the-art ring removal methods via sinogram domain,including polar coordinate based method [9],wavelet-Fourier correction [13],GAN-based method[11],median filtering-based method [19],and polynomial interpolation method [20],and the correction images are shown in Fig.6,of which all images are set with the same window width and level,and the first column in each row is in original size,the 2ndand 3rdcolumn are the enlarged part of the original one.In Fig.6,the first row (a) shows original images of the motherwort’s rhizome,and the 2ndto 6throws (b)-(g)are the results of different ring removal methods.Note that the result by polar coordinate based method (b)cannot completely remove ring artifacts,as the ring artifacts voxels have the different gray values.The results by wavelet-Fourier filtering method (c) and GAN-based method (d) is accompanied with blurry effect,and the result by median filtering method (e) is accompanied with sharpening effects.The residuals show that the loss of original information is worse in these two.The result of the polynomial interpolation method (f),which has heavy radioactive artifacts around the ring artifacts.The result by the proposed method (g) shows that the high frequency information is kept well,which outperforms the other five in ring artifact correction.

4 Conclusion

Ring artifact is one of the most common artifacts in CT image.This paper develops a new method to suppress ring artifacts by locating the positions of the defective pixels by inconsistent pixel value of detector and correcting the error response of them.The detecting algorithm is suitable for both single and continuous stripes.Compared with other ring removal methods,the proposed method has an advantage on removal of ring artifacts in CT images without less impact in spatial resolution and contrast.

Acknowledgement:This work was supported in part by National Key R&D Program of China(2017YFA0104302) and National Natural Science Foundation of China (61871126).