Image Processing for Low-Dose CT via Novel Anisotropic Fourth-Order Diffusion Model

Low-dose CT images contain severe mottle noise and streak artifacts, which seriously affect the physician’s diagnosis of the disease. Hence, in this paper, we propose a novel anisotropic fourth-order diffusion model for low-dose CT image processing. The proposed diffusion model uses both image gradient magnitude and weighted residual local energy to determine the diffusion coefficient. Gradient magnitude is used to detect the image edges, while the weighted residual local energy preserves textures and details in the image. In addition, the fidelity term is introduced into the diffusion model to avoid excessive smoothing and weaken the blocky effects. Experimental results show that when compared with the anisotropic fourth-order diffusion model, the proposed algorithm protects the texture details and suppresses the blocky effects. In comparison with other state-of-the-art algorithms, the proposed model effectively suppresses mottle noise and streak artifacts while simultaneously improving the low-dose CT image quality.


I. INTRODUCTION
In recent years, computed tomography (CT) has been used as an effective auxiliary tool in clinical diagnosis and disease detection, as it can provide clear images without tissue overlap. However, the large amount of radiation applied during the CT scan is likely to induce cancer, leukemia, and other diseases, causing serious harm to patients [1], [2]. Therefore, numerous researchers have attempted to address this concern, and the as low as reasonably achievable (ALARA) principle [3] was proposed and implemented in clinical practice. Naidich et al. [4] first proposed the concept of low-dose CT (LDCT); subsequently, the LDCT technology has undergone rapid development. The general methods to lower the radiation dose are as follows: using a sparse angle to obtain the projection data; reducing the tube voltage and current of the X-ray tube; improving the hardware conditions of the CT system, etc. Among them, reducing the tube current is most commonly used in clinical practice because it is easy to implement and has limited influence on the reconstruction results. Presently, low-dose spiral CT is an effective method for early screening of lung cancer in middle and old aged smokers. Moreover, the technology has been extended to other detection fields. However, the low-dose CT technology adds noise to the projection data, causing severe noise in the reconstructed image, which affect the physician's diagnosis of the disease.
Denoising in the projection domain is the most direct method to improve the low-dose CT image quality. After removing the noise from the projection data, an approximate standard dose CT (SDCT) image can be reconstructed by the filtered back projection algorithm. The iterative reconstruction algorithms mainly reduce noise by optimizing the objective function with prior terms, and the design of the prior terms becomes the key to constructing the objective function. However, both these methods require projected data, which restricts their development.
The post-processing algorithm removes noise from the reconstructed images; therefore, these do not require projection data and are not affected by hardware devices. This has significantly improved their applicability and attracted the attention of numerous researchers. Based on the features of low-dose CT images, researchers have modified the classical denoising algorithm and proposed methods such as dictionary learning [5], [6], guided filtering [7], non-local mean [8], weighted kernel norm minimization [9], and BM3D [10], [11]. Meanwhile, deep learning-based algorithms [12]- [16] have attracted much research attention due to their excellent noise reduction performance. However, the noise reduction results are heavily dependent on the training dataset, which has a problem of insufficient robustness; moreover, training the network also requires a significant amount of time.
In recent years, algorithms based on partial differential equation (PDE) have achieved good results in image noise reduction. Perona and Malik proposed the famous anisotropic diffusion model, namely the PM model [17], and Rudin et al. proposed the total variation (TV) model [18]. They pioneered the PDE in the field of image edge preservation and noise reduction, but there are problems such as the generation of blocky effects and inaccurate identification of weak edges. Since then, researchers have proposed several improved algorithms. Wang et al. [19] modified the PM model combine directional Laplacian operator, which effectively alleviated the blocky effects. Chao et al. [20] fused the image local variance with the diffusion coefficient to better retain the image edges and fine details. Rafsanjani et al. [21] incorporated the residual local energy into the PM model to better preserve the texture and details in the diffusion process. To suppress the blocky effects caused by traditional secondorder PDE, many researchers have utilized fractional PDE and fourth-order PDE. You and Kaveh presented a fourth-order isotropic diffusion model, namely the YK model [22], which avoided the blocky effect but produced speckle noise. Hajiaboli [23] first replaced the Laplacian operator with the gradient magnitude in the diffusion coefficient to attenuate the speck noise in the YK model. In a subsequent study, Hajiaboli [24] proposed the anisotropic fourth-order diffusion (AFOD) model, in which the diffusion strength varies according to the direction of the gradient and edge to achieve a better edge protection effect. However, when the noise intensity increases, blocky effects are generated. The researchers also applied the PDE methods to achieve noise reduction in low-dose CT images. Mendrik et al. [25] proposed a denoising algorithm for low-dose CT images by combining edge enhancement and coherent enhancement diffusion. Wang et al. [26] combined the fractional order PM model with fractional order TV model and presented a novel fractional order diffusion model (FPMTV) to improve the quality of low-dose CT images. Liu et al. [27] improved the PM model and fused the image variance and residual image variance with the diffusion coefficient, named ASNDF model, which effectively suppressed the streak artifact and mottle noise in low-dose CT images. Based on the above analysis and inspired by the idea presented in [21], we propose a novel fourth-order diffusion model of weighted residual local energy applied to low-dose CT image noise reduction. The model is named the novel anisotropic fourth-order diffusion (NAFOD) model, which overcomes the shortcomings of the AFOD [24] model, i.e., generating blocky effects and neglecting the protection of texture and details. The NAFOD model effectively suppresses the noise in low-dose CT images.
This paper is organized as follows: Section II briefly introduces the relevant fourth-order PDE model. Section III introduces the algorithm and explains its implementation in this study. Section IV presents the parameter setting and experimental analysis, and section V summarizes the study.

II. RELATED FOURTH-ORDER PDE MODELS
Although the PM model and the subsequent improved secondorder PDE algorithm have achieved good results in the field of image noise reduction, their denoising results often contain blocky effects, which deteriorate the visual effect and blur the image details. Researchers have been exploring methods to remove the blocky effect, and the higher order PDE model is one of the main techniques.

A. YK MODEL
Taking the Laplace operator of the image as the energy function, You and Kaveh [22] presented a fourth-order isotropic diffusion model, the energy function is where 2  represents the Laplacian operator and ( ) f  is a non-negative increasing function. The Laplace operator is zero in the smooth region, therefore, minimizing the above equation is equivalent to smoothing the image and removing noise.
The gradient descent flow of the above equation can be obtained as follows Experimental results revealed that the YK model converges to a piecewise planar image, and no blocky effects are generated. However, there is severe speckle noise in the denoising results.

B. ANISOTROPIC FOURTH-ORDER DIFFUSION MODEL
Hajiaboli pointed out that the YK model is an isotropic diffusion model and proposed the AFOD [24] model. In the diffusion coefficient function, the gradient magnitude | | u  is used instead of the Laplace operator 2 | | u  , obtaining faster convergence rate and reducing speckle noise. The diffusion model is expressed as where 2  is the Laplace operator, u  and u  are the second order derivatives of the image in the directions of the gradient and edge respectively, and 0 (| |) 1 c u    is the diffusion coefficient function, which decreases monotonically with the image gradient magnitude. In the diffusion process, the local properties of the image are considered, and the diffusion strength in the image gradient direction is weaker than in the edge direction, therefore, this model can obtain stronger edge protection. However, the AFOD model has two major flaws: firstly, the blocky effects are caused by the inconsistency of the diffusion strength in the directions of the gradient and edge, the blocky effects are more serious when the noise is severe. Secondly, the gradient magnitude is used as the edge detector, which ignores the protection of the texture and fine details.

III. PROPOSED MODEL
Inspired by the concept presented in [21], the proposed NAFOD model uses the residual local energy as the texture and detail detector, and image gradient magnitude as the edge detector in the diffusion coefficient. Because the residual local energy and image gradient magnitude are insensitive to mottle noise and streak artifact in low-dose CT images, the diffusion strengths in the noise and artifacts region is large, which achieves the effect of suppressing the mottle noise and streak artifact while preserving the valid information. Different diffusion strengths in the directions of the edge and gradient lead to better edge preservation capability. At the same time, the fidelity item is added to avoid excessive smoothness and reduce the blocky effect.

A. RESIDUAL LOCAL ENERGY
Some denoising algorithms such as the PM model [17], YK model [22], and AFOD model [24] remove some valid information from the images, such as weak edges and fine details. Based on this, the local residual energy is used to define the texture and detail detector in the proposed algorithm to obtain valid information that was incorrectly removed from the residual image. Unlike [21], we found that the algorithmic implementation of the proposed model, achieved better results using the AFOD model to obtain the residual images than the PM model, hence, the AFOD model was used instead of the PM model. The residual image is defined as: where o u denotes the distorted image, s u denotes the valid information such as texture details, n u includes some noise components. The residual energy can be expressed as where R P , s P , and n P are the local energies of R u , s u , n u , respectively. Assuming that the valid information and noise in the residual image can be separated accurately, the valid information can be returned to the denoising result. The iteration stop condition of the AFOD model is described as follows: where 1   is constant, and 2 n k denotes the noise variance. In this case, the residual image has rich texture and fine details. To remove noise and artifacts from the residual image, a median filter and Gaussian filter are added to the residual local energy, which are described as follows: In the diffusion model, we normalize the residual local energy and propose a new diffusion coefficient function as (10) Fig. 1 shows an example of the residual local energy of the low-dose CT image. We can see the texture and details of the tracheas and blood vessels in the lungs in Fig. 1(a), as well as the mottle noise and streak artifacts caused by the insufficient X-ray dose. After processing with the AFOD model, the residual image is obtained, and processed by the median filter and local Gaussian filter. The residual local energy image mainly contains edges and texture details wrongly removed as shown in Fig. 1(b), and there is almost no streak artifact, indicating that the residual local energy is insensitive to the artifacts in the low-dose CT image. Therefore, as the texture detail detector, the residual local energy can protect the image texture details and restore some of the valid information to the denoising results wrongly removed by the AFOD model.

B. NOVEL ANISOTROPIC FOURTH-ORDER DIFFUSION MODEL
According to the literature [24], we propose the following objective function with a fidelity term: where ( ) f  is a non-negative increasing function. Minimizing the above equation, the gradient descent flow is obtained by the Euler equation as follows: The approach of the literature [24] is adopted, using the gradient magnitude in the diffusion function and by dividing the diffusion process into the gradient direction and the edge direction. Thereafter, combined with the residual local energy, the model is proposed as follows:  Compared to the AFOD algorithm, the NAFOD algorithm has two main improvements. First, regarding the diffusion coefficient, the weighted texture detail detection operator is used in the model, which offers better protection for low-dose CT images rich in texture and detail. Second, the evolving image u and distorted image 0 u are added as fidelity items to avoid excessive smoothing and reduce the blocky effects.

C. NUMERICAL SOLUTION OF THE MODEL
The finite difference scheme is used to solve the proposed model. Assuming a time step size t  and space grid size =1 h , the discrete time and space coordinates can be expressed as , 0 ,1 ,2 , 0 ,1 ,2 , 0 ,1 ,2 the difference operator is , 1 , with symmetric boundary conditions , 1 , the image gradient magnitude is u  and u  can be expressed as x xx x y xy y yy x y y xx x y xy x yy the discrete expressions for the residual local energy and diffusion coefficient are as follows: , , Finally, the discrete result is obtained as follows: Applying the above discrete scheme, the specific steps of the model are given below.

IV. EXPERIMENTAL RESULTS
To evaluate the performance of the proposed NAFOD algorithm, we tested it on the Mayo public dataset [28], the low-dose CT images were reconstructed using standard-dose projection data added with Poisson noise, simulating a 1/4 SDCT image. The CT images of the abdominal cavity, chest cavity, and pelvic cavity were used in the following experiments.
In the quantitative assessment of images, the peak signalto-noise ratio (PSNR) is an effective criterion for noise evaluation, the larger the value, the smaller is the image noise Mean structural similarity (MSSIM) [29] is built on the human eye perception system, it combines the image brightness, contrast, and structural information. The larger the value, the greater is the similarity between the denoising image and the original image in terms of structure. Feature similarity (FSIM) [30] and gradient magnitude similarity deviation (GMSD) [31] were also used as quality indices of the noise reduction performance in low-dose CT images. FSIM is an improvement over the MSSIM, and the larger the value, the better protected are the feature parts. The GMSD is proposed due to the feature that the image gradient is more sensitive to image degradation, and the smaller the value, the better the denoising image preserves the gradient information.
The proposed algorithm was compared with the YK model [22], AFOD model [24], FPMTV model [26], ASNDF model [27], and REDCNN [32]. The YK and AFOD models are fourth-order PDE models; FPMTV and ASNDF models are the PDE algorithms for processing low-dose CT images, and the REDCNN was proposed by Chen et al. as an excellent deep learning based method for low-dose CT image processing. We used the above mentioned PSNR, MSSIM, FSIM, and GMSD as the quantitative evaluation indexes. The parameter setting and analysis of the experimental results are presented below.

A. PARAMETER SETTING
The proposed NAFOD model mainly involves two parameters: the residual weight coefficient  in Eq. (10) and the weight coefficient  of the fidelity term in Eq. (11). In the AFOD model, the time step t  proves that when 0.0313 t   , the model has a stable solution. Therefore, in the proposed model, the time step is set 0.03 t   . We initially used the method applied in literature [33] to estimate the noise variance 2 n k in Eq. (7), due to the complicated noise of the low-dose CT image, the noise estimation was not ideal; therefore, we used the image processed by Gaussian filtering as the reference standard dose image to estimate the noise variance and achieve better results, the size of the Gaussian operator was 5 5  , and 5   . In Eq. (7), the noise intensity coefficient 1  , and 5   . Because the dataset contains SDCT images, we used the MSSIM as the algorithm iteration stopping criterion to retain the best structural properties of the denoising results. In the comparison experiments, we determined the parameters by MSSIM metric and fine-tuned them according to visual effects. MSSIM was also used as the iteration stopping criterion. The REDCNN used 760 sets of images from the Mayo dataset for network training, with the learning rate set to -4 10 and batch-size of 16. Table 1 shows the model parameter settings.  In the next section, the abdomen CT image in the folder L014 will be used as an example to introduce the parameter settings of the proposed model and explain the effect of the parameter settings on the denoising performance. We determined the weight coefficient  and  according to the evaluation index and visual effect of the processed result.

1) Residual local energy weight coefficient α
Literature [21] proposed the residual local energy, and the ASNDF model added the residual local energy to the diffusion coefficient, but the weight was simply set as 1. In this study, different weights were selected according to the image properties to enhance the protection of the texture details and improve the image quality. We can see that the MSSIM and PSNR curves increase first and then decrease in Fig. 3(a), which indicates that the quality of the denoising image improves first and then deteriorates. When =2  , MSSIM has the maximum value and PSNR has a large value, hence, it is considered the optimal weight coefficient. In this case, the details and edges are clear, as shown in the zoomed regions of interest (ROI) in Fig. 4(b) (indicated by the green arrow). When =0  , we can see that some image edges are blurred, i.e., the interior tissue has noise in Fig. 4(a). When  is large, the noise is enhanced, as indicated by the red arrow in Fig. 4(c). It can also be seen from Fig. 5 that the best metrics values are obtained when =2  .

2) Fidelity term weight coefficient λ
Added the fidelity term can avoid excessive smoothing and reduce the blocky effects. The MSSIM and PSNR curves increase first and then decrease in Fig. 3(b). When =0.6  , both MSSIM and PSNR have large values; therefore, it is regarded as the optimal parameter value. In such a case, the processed image has good noise reduction effect, as shown in Fig. 6(b). When =0  , from Fig. 6(a), a small number of blocky effects is generated (indicated by the yellow arrow).
When  is large, noise is introduced in the processing result, as shown by the red arrow in Fig. 6(c). It can also be seen from Fig. 7 that the best metrics values are obtained when =0.6  . In general, an image with a high level of texture and detail should have a large residual local energy weight coefficient  . The more severe the noise, the smaller is the fidelity term weight coefficient  . Based on numerous experimental results, we give the reference ranges for two parameters [1,6 ), (24) We used Matlab 2016a to implement the proposed model, using a CPU with Intel i7-9700K @ 3.6GHz, 32GB of RAM, and an NVIDIA GeForce RTX 2080s GPU. From Fig. 2, it is evident that the time complexity and space complexity of the proposed model are consistent with the literature [24]. To be precise, let the number of image pixels be n . The proposed model contains two iterative processes, the first iteration determines the residual image, the second iteration is to solve the result image, and the time and space complexity are ( ) O n for both iterations. Therefore, the time and space complexity of the proposed model are ( ) O n . For the Mayo dataset images used in the experiment, the time to process a low-dose CT image is approximately 10 seconds using the GPU parallel acceleration toolbox in Matlab, since the time step t  is a small value and the number of two iteration processes add up to approximately 2000 times.

B. VISUAL ANALYSIS
The results of the image processing are depicted below, and the display window of all the images is [-160,240] HU. Fig. 8(a1) shows a low-dose abdominal cavity CT image selected from the folder L014. The image contains mottle noise and a small amount of streak artifacts, which affect the detection of lesions and classification of tissues and organs. Fig. 8(c1-h1) show the processed results from different models. Fig. 8(a2-h4) show the corresponding zoomed ROIs of Fig. 8(a1-h1), respectively. The three ROIs contain the liver region with lesion, the marginal region of the spleen, and the smooth area of the liver with veins. The YK model overcomes the blocky effect, there is severe speckle noise in the image, indicated by the red arrow in Fig. 8(c2), moreover, the YK model blurs the edge of the image. Blocky effects exist in the AFOD model at the tissue edge (indicated by the yellow arrow), as shown in Fig. 8(d3). The FPMTV model has the problem of incomplete noise reduction, as shown in Fig. 8(e4) (indicated by the red arrow). The ASNDF model has a high similarity with the SDCT image, but it also has the problem of a small amount of residual noise. The overall visual effect of REDCNN is good, however, there is the problem of blurred edges (indicated by the blue arrow in Fig. 8(g3)). In this study, the proposed algorithm demonstrated the best visual performance. There is almost no blocky artifact, and the edge and small details are retained. In Fig. 8(h2), the lesion is clearly visible; in Fig. 8(h3), the spleen edge is better protected (indicated by the green arrow). Denoising performance is evident in the smooth area in Fig. 8(h4), moreover, the liver's blood vessels are clearly visible. 2) Low-dose chest CT Fig. 9(a1) depicts a low-dose chest cavity CT image selected from the folder L019, the ROIs in the red box contain the artery edge region and the smooth area on the right side of the heart. From the processed results, the YK model has severe speckle noise (indicated by the red arrow in Fig. 9(c3)). There are a few blocky effects in the AFOD model (indicated by the yellow arrows in Fig. 9(d2) and (d3)). Residual noise exists in the FPMTV model, which is evident in the interior vascular and smooth areas of the heart (indicated by the red arrows in Fig. 9(e2) and (e3)). The ASNDF model had a good protective effect on the artery edge (indicated by the green arrow in Fig.  9(f2)), whereas it produced more serious blocky effects than the AFOD model on the cardiac edge (indicated by the yellow arrow in Fig. 9(f3)). Despite the high similarity between the REDCNN and SDCT images, there is still the problem of blurred edges as shown in Fig. 9(g3) (indicated by the yellow arrow). The NAFOD model achieved a good balance between image noise suppression and edge protection; the inner area of the blood vessel is smooth (indicated by the direction of the green arrow in Fig. 9(h2)), and the heart edge is clear (indicated by the green arrow in Fig. 9(h3)), in addition, almost no blocky effect is generated in the entire image.

1) Low-dose abdominal cavity CT
3) Low-dose pelvic CT Fig. 10(a1) depicts a low-dose pelvic cavity CT image selected from the folder L006. The ROIs in the red box contain the hip bone region and pelvic tissue region. We see that a large number of speckles exists in the YK model (indicated by the red arrow in Fig. 10(c3)). The AFOD model has good edge protection effect, but a small number of blocky effects are generated in the pelvic tissue (indicated by the yellow arrow in Fig. 10(d3)). The FPMTV model and ASDNF model have the problem of residual noise, as indicated by the red arrows in Fig. 10(e3) and (f3). The REDCNN has strong noise suppression and no blocky effect; however, its edges are blurred, as indicated by the yellow arrow on the edge of the tissue in Fig. 10(g3). The proposed algorithm achieved the best results for the pelvic cavity image, the noise in the smooth area of the muscle was completely reduced (indicated by the green arrow in Fig. 10(h2)), Fig. 10 (h3) shows that there are clear tissue edges and the blocky effect is not produced. The quantitative result of the proposed algorithm and competing algorithms are presented in Table 2. It is evident that only the REDCNN has some higher metrics than the proposed algorithm. The REDCNN is effective at suppressing noise and has metrics that closely resemble the proposed model. The proposed NAFOD model achieves high metrics and it exhibits the best performance among all the PDE methods.

C. QUANTITATIVE ASSESSMENT
Considering the ROIs contain more valuable information, we selected three ROIs in the red box from the abdominal cavity low-dose CT image in Fig. 8(a1). The calculated values of MSSIM and PSNR are shown in Fig. 11. In the three ROIs, the MSSIM and PSNR of the proposed model achieved the highest values among all PDE methods, which confirms that the NAFOD model is effective in reducing noise and maintaining the structures. While the REDCNN still achieves good results, the proposed model has better metrics. Thus, the visual effects and quantitative assessment show that the proposed NAFOD model demonstrates the best result in noise suppression and structure maintain of low-dose CT images.

V. CONCLUSIONS
This paper proposes a novel anisotropic fourth-order diffusion model for suppressing mottle noise and streak artifacts in lowdose CT images. The edges of the low-dose CT image are detected using the image gradient magnitude, whereas the texture details are detected by the residual local energy. The parameter  is introduced to detect texture details more accurately to achieve the effect of preserving the valid information and suppressing the noise. In addition, a fidelity term is introduced into the diffusion model to avoid excessive smoothing and weaken the blocky effects. The experiment of low-dose CT images of the abdominal cavity, thoracic cavity, and pelvic cavity from the Mayo open dataset verifies the effectiveness of the proposed algorithm. The proposed model has certain flaws: the criterion for stopping the iterations utilizes the MSSIM standard to make better preservation of image structure, but this criterion relies on reference images. Moreover, the MSSIM standards introduce uncertainty into the number of iterations, therefore, the model running time is not only related to the image size, but also to the image structure. Besides, as exposures are reduced further, low-dose CT images contain more severe noise, the proposed model will gradually show the problem of noise residue, in which case, the parameter k does not conform to the adaptive formula, and needs to be manually adjusted. Our future work will focus on the iterative stopping condition and adaptive parameters selection. We have made the code publicly available to facilitate future studies [35].