Regularized Multidirections and Multiscales Anisotropic Diffusion for Sinogram Restoration of Low-Dosed Computed Tomography

Although most of existing anisotropic diffusion (AD) methods are supported by prefect mathematical theories, they still lead to smoothed edges and anatomy details (EADs). They are caused by not considering the discrete nature of digital signal. In order to improve the performance of AD in sinogram restoration of low-dosed computed tomography (LDCT), we propose a new AD method, named regularized multidirections and multiscales anisotropic diffusion (RMDMS-AD), by extending AD to regularized AD (RAD) in multidirections and multiscales. Since the multidirections can reduce the discrete errors to the maximum extent, meanwhile multiscales and RAD make searching neighborhood of solution be as large as possible which can get more optimal solution to AD, the new proposed method can improve the performance of AD both in denoising and in stability of solution. Moreover, the discrete errors and ill-posed solutions occur mostly near the EADs; the RMDMS-AD will also preserve EADs well. Comparing the proposed new method to existing AD methods using real sinogram, the new method shows good performance in EADs preserving while denoising and suppressing artifacts.


Introduction
Anisotropic diffusion (AD) proposed in 1990, also called Perona-Malik diffusion (PMD), is a technique aiming at reducing image details without removing significant parts of the image contents, typically edges, lines, or textures which are important for the image [1]. Formally, AD is defined as ( , , ) = div ( ( , , ) ∇ ( , , )) , where ( , , 0) is the initial gray scale image, ( , , ) is the smoothed gray scale image at time , ∇ denotes the gradient, div(⋅) is the divergence operator, and ( , , ) is the diffusion coefficient. ( , , ) controls the rate of diffusion and is usually chosen as a monotonically decreasing function of the module of the image gradient. The coefficients usually take the following two forms: ( ∇ ( , , ) ) = 1 where ‖ ⋅ ‖ is the module of the vector and the constant controls the sensitivity to edges.
The key idea of PMD is to smooth the homogenous regions with small ‖∇ ( , , )‖ while, near singularities with big ‖∇ ( , , )‖, PMD is only smootheed along the perpendicular direction of the gradient.
So far, much research has been devoted to improving the performance of PMD from the view of mathematical theory [2][3][4][5][6][7][8][9][10][11]. In 1992, Catté et al. indicate that PMD is ill-posed, and they also propose a new well-posed method named regularized anisotropic diffusion (RAD) [2]. 2 Computational and Mathematical Methods in Medicine Some recent efforts also focus on the illness of PMD and propose that posing some conditions on diffusion coefficients will make the PMD become well posed [3][4][5][6]. These conditions include posing edge influence functions near edges [3] and substituting the module of the gradient operator to a nearby differentiable function [6].
Since the main idea for AD is that diffusion should be along the normal near singularity, how to approximate to directions of gradient and normal (DsGN) in discrete data turns very important in designing AD filtering [1,2]. From the coherence filtering proposed by Weickert [8] in 1999, vector tensor becomes a preferable alternative for DsGN since the tensor not only can be extended to high dimensional data directly but also can avoid computing DsGN [9][10][11].
Essentially, tensor is also a 2-dimension method for data analysis, which can not provide enough directions for diffusing in discrete situation. Thus, diffusion according to tensor will also make serious artifacts, which hampers AD to be used for low-dose CT imaging.
Recently, how to minimize the radiation exposure to patients has been one of the major efforts in modern clinical X-ray CT radiology because of the associated risk of cancer for patients receiving CT examination [12][13][14][15][16]. But LDCT results in the presentation of strong noise and artifacts which will degrade the quality of LDCT images dramatically and decrease the accuracy of diagnosis.
Filtering noise from clinical scans is a challenging task since these scans contain artifacts and noise and consist of many structures with different shapes, sizes, and contrast, which should be preserved for making correct diagnosis. Many strategies have been proposed to improve the imaging quality of LDCT [17][18][19][20][21][22][23][24]. These methods include some stateof-the-art techniques, such as nonlocal means [18].
Very recently, a new AD, named time and space AD (TSAD), which extends classical space AD to time and space AD by adding an extra diffusion dimension of time, is proposed and shows promising potential in preserving singularities of LDCT scans while denoising [25]. Because much more neighbor pixels will be used in searching solution of ( , ), TSAD can obtain more optimal solution than traditional AD and RAD.
However, the searching neighborhood of TSAD is still not big enough especially near singularities, which will blur the edges. Therefore, how to expand searching neighborhood for resolution of AD is a key problem to improve performance of AD. In this paper, we argue that diffusion should be carried on multiscales in multidirections in order to search a solution in larger neighborhood and to reduce discrete errors.
It should be indicated that multiscales in this paper are different to multiscales in multiscale analysis, which used one scale in a time while, in this paper, multiscales are used simultaneously to extend one scale searching neighborhood, whose details will be introduced in Section 3. The proposed method also is different from RPM since RPM only regularizes diffusion coefficients while our method provides more pixels for diffusion itself.
The arrangement of this paper is as follows. In Section 2, the backgrounds about noise model of LDCT and AD are introduced; and then the new diffusion method is given in Section 3; the experiment results are shown and discussed in Section 4; the final part is the conclusions and acknowledgment.

Backgrounds
In this section some related backgrounds will be introduced. These backgrounds include noise models and numerical schemes for AD and RAD.

Noise Model.
Poisson noise (PN) and Gaussian dependent noise (GDN) after some transform of original sinogram are generally accepted as two noise models of LDCT [26,27]. These two models are built both on repeated phantom experiments and theory analysis.
The photon noise is due to the limited number of photons collected by the detector [27]. For a given attenuating path in the imaged subject, 0 ( , ) and ( , ) denote the incident and the penetrated photon numbers, respectively. Here, denotes the index of detector channel or bin and is the index of projection angle. In the presence of noises, the sinogram should be considered as a random process and the attenuating path is given by where 0 ( , ) is a constant and ( , ) is Poisson distribution with mean . Thus, we have Both its mean value and variance are . Gaussian distributions of ployenergetic systems were assumed based on limited theorem for high-flux levels, and followed by many repeated experiments in [26], we have where is the mean and 2 is the variance of the projection data at detector channel or bin , is a scaling parameter, and is a parameter adaptive to different detector bins. The most common conclusion for the relation between Poisson distribution and Gaussian distribution is that the photon count will obey Gaussian distribution for the case with large incident intensity and Poisson distribution with feeble intensity [27]. In addition, in [26], the authors deduce the equivalency between Poisson model and Gaussian model. Therefore, both theories indicate that these two noises have similar statistical properties and can be unified into a whole framework.

The Numerical Schemes for PM and RPM.
In Section 1, AD is defined as in (1)-(3), which encourage diffusion (hence smoothing) within regions and stop it near strong edges. Hence the edges can be preserved while smoothing from the image [1]. Medicine   3 The discretization for Laplacian operator is

Computational and Mathematical Methods in
where where ‖ ⋅ ‖ is the absolute value of the number and (⋅) is defined in (2) or (3). However, PMD is an ill-posed equation and it can be well posed by RAD [2], which is defined as Here 1 defined as is a Gaussian function and is a constant. The diffusion coefficients (⋅) are defined in (2) or (3). The discretization for Laplacian operator is defined in (7) and (8). But the module of the gradient is simplified as the absolute values of smoothed four directions and diffusion coefficients are where ‖ ⋅ ‖ is the absolute value of the number, 1 ⋅ is defined in (11), and (⋅) is defined in (2) or (3).

The New Method
In this section, related definitions of multidirections for the first and second scale are given firstly; then the general definitions for the th scale ( > 2) are provided. Finally, the new numerical scheme of AD is proposed by mixing the multi scale information together.

Multidirections-Enclosed Difference of the First Scale.
Intuitively, more directions will obtain better performance in image smoothing. In order to extend four directions to the most directions, enclosed directions for ( , ) which are defined as the differences between ( , ) and its enclosed nearest neighbors in [28] are adopt. For example, the differences between ( , ) and its eight-nearest-neighbors are their first-scale enclosed directions (see Figure 1(a)).
Although the motivation of multidirections is very simple, its theory is based on the high dimensional vector analysis for AD. The high dimensional vector used in it is half-point differences between the center point ( , ) and its nearest eight or more half points. The first-scale enclosed gradient vector is defined as Thus the first-scale second-order difference of ( , ) is Computational and Mathematical Methods in Medicine The first-scale second-order enclosed difference of ( , ) The second-scale second-order enclosed difference of ( , ) Figure 1: The second-order enclosed difference of ( , ).
where represents the transpose of the vector. From (13) we have It can be easily concluded that the first scale second-order difference of ( , ) is composed by eight components which are the differences between ( , ) and its enclosed eightnearest-neighbors. Except for the points on border of an image, each point is enclosed by its eight-nearest-neighbors, so the first-scale second-order difference of ( , ) is called the enclosed first-scale second-order difference of ( , ). Definition 1. The first-scale second-order-difference of ( , ) defined in (15) is called the first-scale second-order enclosed difference of ( , ).
Based on Definition 1, we have the following.

Definition 2.
Let the first-scale second-order difference of ( , ) defined in (15). The following equation is called the first-scale enclosed Laplacian operator (ELO) of ( , ).

Multidirections of the th
The second-scale second-order difference of ( , ) is Similar to previous discussion, the second-scale secondorder difference is defined as follows.
Definition 3. The second-scale second-order difference of ( , ) defined in (20) is called the second-scale second-order enclosed difference of ( , ).

Based on Definition 3, we have the following.
Definition 4. Let the second-scale second-order difference of ( , ) defined in (20). The following equation is called the second-scale enclosed Laplacian operator (ELO) of ( , ).
In order to discuss the general definition of multidirections of the th scale ( > 2), we have to determine the number of multidirections for the th scale firstly: Therefore, the th scale enclosed gradient vector is defined as The th scale second-order difference of ( , ) is where represents the transpose of the vector. We have The th scale second-order difference is defined as follows.
Definition 5. The th scale second-order difference of ( , ) defined in (26) is called the th scale second-order enclosed difference of ( , ).
Based on Definition 5, we have the following. Definition 6. Let the th scale second-order difference of ( , ) defined in (26). The following equation is called the th scale enclosed Laplacian operator (ELO) of ( , ).

The New Numerical Scheme.
Although the form of the PDE is the same as the RAD model equations (10)-(11), the numerical scheme in this paper is quite different from the RAD. Let where > 2 is an integer, represents the transpose of the vector and , = 0, . . . , (2 + 1) 2 − 1, is defined as where 1 is defined in (11), ∇ ( , ), = 0, . . . , (2 + 1) 2 − 1, defined in (24) are the components of vector ⃗ ∇ ( , ), is the normalized constant, and is the decreasing function of absolute value of ∇ ( , ), = 0, . . . , (2 +1) 2 −1. Following (2) and (3) where |⋅| is the absolute value of the number and the constant controls the sensitivity to edges. However, half-point difference for ⃗ ∇ ( , ) defined in (24) cannot be computed directly. Thus second-order integral point difference can be used to approximate the one-order half-point difference. Equations (31)-(32) become The new AD based on the enclosed second-order difference is defined as . . .

Experiments and Discussion
The main objective for smoothing LDCT images is to delete the noise while preserving anatomy details for the images. Experiments in this section are designed to check performance of different AD methods in LDCT imaging which used both real LDCT sinogram and standard-dose CT (SDCT) sinogram, where sinogram of SDCT will provide visual reference for recovered LDCT images.

Data. Three abdominal CT images with different doses
were scanned from a 16-multidetector row CT unit (Somatom Sensation 16; Siemens Medical Solutions) using 120 kVp and 5 mm slice thickness ( Figure 2). Other remaining scanning parameters are gantry rotation time, 0.5 second; detector configuration (number of detector rows section thickness), 16 × 1.5 mm; table feed per gantry rotation, 24 mm; pitch, 1 : 1; and reconstruction method, Filtered Back Projection (FBP) algorithm with the soft-tissue convolution kernel "B30f. " Different CT doses were controlled by using two different fixed tube current 30 mAs and 150 mAs (or 60 mA and 300 mAs) for LDCT and standard-dose CT (SDCT) protocols, respectively (see Figures 2 and 3). The CT dose index volume (CTDIvol) for LDCT images and SDCT images is in positive linear correlation to the tube current and is calculated to be approximately ranged from 15.32 mGy to 3.16 mGy [18].

Compared
Methods. PM and RPM have been discussed in detail in Section 2. In this subsection, we will briefly introduce two new AD methods in [3,11].
In [3], the authors propose a new method to make AD well posed by posing edge constraint. That is, edge is detected by canny operator to form an edge map; then a Gaussian smoothing is carried on the edge map to design a fuzzy edge The authors of [11] propose a hybrid scheme that combines forward and backward differences with AD by tensor. The tensors obtained by our approach depend on four directional derivatives of the intensity of an image, and hence they are adaptively determined by local image structure. The proposed diffusion filter is isotropic in the interior of a region, whereas it is anisotropic at edges.

The Parameter
Choice. There are five methods will be compared in this paper: PM, RPM, the method proposed in [3] which is named as well-posed PM (WE-PM), the method proposed in [11] which is named as new tensor PM (WT-PM), and our method (RMDMS-AD). In order to ensure that comparison is put on a fair level, the common used parameters are set to the same value.
The common used parameters for five methods include gradient modulus threshold that controls the conduction, integration constant (0 ≤ ≤ 1/7), and iteration number .
Due to numerical stability, is set to its maximum value 1/7.
Intuitively, the bigger the is, the smoother the recovered image is. Therefore, should be set to a small number. In this paper, it is set to 2 to preserve EADs.
The iteration number is very important in AD. That is, too big will make over-smooth image while too small will still leave many noises. In order to study the performance of five compared methods with different iteration numbers and fixed other parameters, is set to 3, 5, 10, and 20, respectively.
The standard deviation of smoothed Gaussian kernel for the image 1 used for RPM, WT-RPM, and RMDMS-AD is also set to 1.5 since, in [3], the authors suggest that 1 should be a small number. Another parameter is the number of scale , which is very important to support the start point of our method where the big scale will have better performance than the small scale. Thus, the is set to 1, 2, and 3 to study the different behavior in different scales. Figure 4, noises in all reconstructed images are suppressed in different extent and EADs in these methods are preserved, which show that the proposed method is a promising LDCT imaging way. It also can be seen that the multiscale has better performance in preserving EADs. That is, from left to right of each row much more EADs are preserved. In addition, we also can observe that big leads to over-smooth image. Thus, from top to down of each column in Figure 4, the smoother and smoother images can be obtained. Like most of existing AD methods, iteration numbers for the proposed scheme should be chosen carefully. Therefore, we will compare performance of comparison methods on Section 4.2 with = 5, 10, 20 (see Figure 5). The reconstructed image of WT-PM has many odd artifacts from = 5 to = 20. These artifacts will lead to error diagnosis. Thus WT-PM does not suit for LDCT imaging.

Results and Discussion. From
In Figure 5 the reconstructed LDCT images for PM (the third row of Figure 5 [1]) and WE-PM (the fourth row of Figure 5 [3]) left so many noises even after 20 times of iteration (the last column of related row). The other reconstructed LDCT images also have similar situation, which is caused by the illness of AD equation and serious noise for LDCT.
RPM can be considered as a special case of the proposed method. That is, when scale = 1, RMDMS-AD becomes RPM. Just as discussed in this section, RMDMS-AD has better performance than RPM. Thus RMDMS-AD can preserve more EADs for LDCT, which is very important in clinical diagnosis. RPM and RMDMS-AD can suppress noise

WT-RPM
(o) Recovered LDCT image using WT-PM [11] with iteration number = 20 greatly which remind us that "regularization" is an essential technology in AD.

Conclusions
In this paper, we propose a new AD for LDCT sinogram imaging using multiscales and multidirections, named RMDMS-AD. Since not only multidirections can reduce discrete error greatly but also regularization and multiscales can provide big searching regions for solutions of AD equation, the new proposed method has good performance both in noise suppressing and DsGN preserving. Moreover, RMDMS-AD does not produce new artifacts in denoising which is very important in clinical diagnosis.