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

Abstract

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.

1. 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𝑥,𝑦,𝑡∇𝑢𝑥,𝑦,𝑡 (1) 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: 𝑔∇𝑢𝑥,𝑦,𝑡=𝑒^{−(‖∇𝑢(𝑥,𝑦,𝑡)‖/𝜎)2},𝑔=1∇𝑢𝑥,𝑦,𝑡1+∇𝑢𝑥,𝑦,𝑡/𝜎², (2)(3) 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–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].

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–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–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–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–24]. These methods include some state-of-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.

2. Backgrounds

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

2.1. 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, 𝑁₀(𝑖,𝛼) 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 𝑟_𝑖𝑁()=−ln𝑖,𝛼𝑁₀(),𝑖,𝛼 (4) where 𝑁₀(𝑖,𝛼) is a constant and 𝑁(𝑖,𝛼) is Poisson distribution with mean 𝑁.

Thus, we have 𝑁(𝑖,𝛼)=𝑁₀(𝑖,𝛼)exp−𝑟_𝑖. (5)

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 𝜎²_𝑖𝜇_𝑖=𝑓_𝑖𝜇exp^𝑖𝛾, (6) where 𝜇_𝑖 is the mean and 𝜎²_𝑖 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.

2.2. 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].

The discretization for Laplacian operator is 𝑢+1𝑖,𝑗,𝑡+1=𝑢𝑖,𝑗,𝑡4𝑐_𝑁⋅∇²_𝑁𝑢𝑖,𝑗,𝑡+𝑐_𝑆⋅∇²_𝑆𝑢𝑖,𝑗,𝑡+𝑐_𝐸⋅∇²_𝐸𝑢𝑖,𝑗,𝑡+𝑐_𝑊⋅∇²_𝑊𝑢𝑖,𝑗,𝑡, (7) where ∇²_𝑁𝑢,∇𝑖,𝑗,𝑡=𝑢𝑖−1,𝑗,𝑡−𝑢𝑖,𝑗,𝑡_2𝑆𝑢,∇𝑖,𝑗,𝑡=𝑢𝑖+1,𝑗,𝑡−𝑢𝑖,𝑗,𝑡_2𝐸𝑢,∇𝑖,𝑗,𝑡=𝑢𝑖,𝑗+1,𝑡−𝑢𝑖,𝑗,𝑡_2𝑊𝑢𝑖,𝑗,𝑡=𝑢𝑖,𝑗−1,𝑡−𝑢𝑖,𝑗,𝑡. (8) According to (1)–(3), the diffusion coefficient is defined as a function of module of the gradient. However, computing a gradient accurately in discrete data is very complex and the module of the gradient is simplified as the absolute values of four directions and diffusion coefficients are 𝑐_𝑁𝑖,𝑗,𝑡=𝑔∇²_𝑁𝑢,𝑐𝑖,𝑗,𝑡_𝑆𝑖,𝑗,𝑡=𝑔∇²_𝑆𝑢,𝑐𝑖,𝑗,𝑡_𝐸𝑖,𝑗,𝑡=𝑔∇²_𝐸𝑢,𝑐𝑖,𝑗,𝑡_𝑊𝑖,𝑗,𝑡=𝑔∇²_𝑊𝑢𝑖,𝑗,𝑡, (9) 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 𝜕𝑢𝑥,𝑦,𝑡𝑔𝐺𝜕𝑡=div_𝜎1⋅∇𝑢𝑥,𝑦,𝑡∇𝑢𝑥,𝑦,𝑡. (10)

Here 𝐺_𝜎1 defined as 𝐺_𝜎1=1𝐶𝑒^{−(𝑥2+𝑦2)/𝜎21} (11) 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 𝑐_𝑁𝑖,𝑗,𝑡=𝑔𝐺_𝜎1⋅∇²_𝑁𝑢,𝑐𝑖,𝑗,𝑡_𝑆𝑖,𝑗,𝑡=𝑔𝐺_𝜎1⋅∇²_𝑆𝑢,𝑐𝑖,𝑗,𝑡_𝐸𝑖,𝑗,𝑡=𝑔𝐺_𝜎1⋅∇²_𝐸𝑢,𝑐𝑖,𝑗,𝑡_𝑊𝑖,𝑗,𝑡=𝑔𝐺_𝜎1⋅∇²_𝑊𝑢𝑖,𝑗,𝑡, (12) where ‖⋅‖ is the absolute value of the number, 𝐺_𝜎1⋅ is defined in (11), and 𝑔(⋅) is defined in (2) or (3).

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.

3.1. 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)).

(a) The first-scale second-order enclosed difference of 𝑢(𝑖,𝑗)

(b) The second-scale second-order enclosed difference of 𝑢(𝑖,𝑗)

Figure 1 

The second-order enclosed difference of 𝑢(𝑖,𝑗).

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 ∇𝑢𝑖,𝑗=∇₀𝑢𝑖,𝑗,∇₁𝑢,∇𝑖,𝑗₂𝑢𝑖,𝑗,∇₃𝑢𝑖,𝑗,∇₄𝑢,∇𝑖,𝑗₅𝑢𝑖,𝑗,∇₆𝑢𝑖,𝑗,∇₇𝑢𝑖,𝑗^𝑇, (13) where 𝑇 represents the transpose of the vector and ∇𝑢_𝑘(𝑖,𝑗), 𝑘=0,…,7 are defined as ∇₀𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗+0.5−𝑢𝑖,𝑗₁𝑢,∇𝑖,𝑗=𝑢𝑖−0.5,𝑗+0.5−𝑢𝑖,𝑗₂𝑢,∇𝑖,𝑗=𝑢𝑖−0.5,𝑗−𝑢𝑖,𝑗₃𝑢,∇𝑖,𝑗=𝑢𝑖−0.5,𝑗−0.5−𝑢𝑖,𝑗₄𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗−0.5−𝑢𝑖,𝑗₅𝑢,∇𝑖,𝑗=𝑢𝑖+0.5,𝑗−0.5−𝑢𝑖,𝑗₆𝑢,∇𝑖,𝑗=𝑢𝑖+0.5,𝑗−𝑢𝑖,𝑗₇𝑢𝑖,𝑗=𝑢𝑖+0.5,𝑗+0.5−𝑢𝑖,𝑗. (14)

Thus the first-scale second-order difference of 𝑢(𝑖,𝑗) is ∇²𝑢𝑖,𝑗=∇²₀𝑢𝑖,𝑗,∇²₁𝑢,∇𝑖,𝑗₂₂𝑢𝑖,𝑗,∇²₃𝑢𝑖,𝑗,∇²₄𝑢,∇𝑖,𝑗₂₅𝑢𝑖,𝑗,∇²₆𝑢𝑖,𝑗,∇²₇𝑢𝑖,𝑗^𝑇, (15) where 𝑇 represents the transpose of the vector. From (13) we have ∇²₀𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗+1−𝑢𝑖,𝑗₂₁𝑢,∇𝑖,𝑗=𝑢𝑖−1,𝑗+1−𝑢𝑖,𝑗₂₂𝑢,∇𝑖,𝑗=𝑢𝑖−1,𝑗−𝑢𝑖,𝑗₂₃𝑢,∇𝑖,𝑗=𝑢𝑖−1,𝑗−1−𝑢𝑖,𝑗₂₄𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗−1−𝑢𝑖,𝑗₂₅𝑢,∇𝑖,𝑗=𝑢𝑖+1,𝑗−1−𝑢𝑖,𝑗₂₆𝑢,∇𝑖,𝑗=𝑢𝑖+1,𝑗−𝑢𝑖,𝑗₂₇𝑢𝑖,𝑗=𝑢𝑖+1,𝑗+1−𝑢𝑖,𝑗. (16)

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 eight-nearest-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 ⁷_𝑘=0∇²_𝑘𝑢𝑖,𝑗 (17) is called the first-scale enclosed Laplacian operator (ELO) of 𝑈(𝑖,𝑗).

3.2. Multidirections of the 𝑛th Scale

In this section, multidirections of the second scale are introduced firstly to enhance the intuition of multidirections for the 𝑛th scale; and then, the related definitions of multidirections of the 𝑛th scale are given.

Just as shown in Figure 1(b), there are 16-nearest-neighbors in the second scale to form enclosed directions. Thus, the second-scale enclosed gradient vector is defined as =∇∇𝑢𝑖,𝑗₀𝑢𝑖,𝑗,∇₁𝑢𝑖,𝑗,…,∇₁₅𝑢𝑖,𝑗^𝑇, (18) where 𝑇 represents the transpose of the vector and ∇𝑢_𝑘(𝑖,𝑗), 𝑘=0,…,15, are defined as ∇₀𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗+1−𝑢𝑖,𝑗₁𝑢,∇𝑖,𝑗=𝑢𝑖−0.5,𝑗+1−𝑢𝑖,𝑗₂𝑢,∇𝑖,𝑗=𝑢𝑖−1,𝑗+1−𝑢𝑖,𝑗₃𝑢,⋮∇𝑖,𝑗=𝑢𝑖−1,𝑗+0.5−𝑢𝑖,𝑗₁₄𝑢,∇𝑖,𝑗=𝑢𝑖+1,𝑗+1−𝑢𝑖,𝑗₁₅𝑢𝑖,𝑗=𝑢𝑖+0.5,𝑗+1−𝑢𝑖,𝑗. (19)

The second-scale second-order difference of 𝑢(𝑖,𝑗) is ∇²𝑢𝑖,𝑗=∇²₀𝑢𝑖,𝑗,∇²₁𝑢𝑖,𝑗,…,∇²₁₅𝑢𝑖,𝑗^𝑇, (20) where 𝑇 represents the transpose of the vector. We have ∇²₀𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗+2−𝑢𝑖,𝑗₂₁𝑢,∇𝑖,𝑗=𝑢𝑖−1,𝑗+2−𝑢𝑖,𝑗₂₂𝑢,∇𝑖,𝑗=𝑢𝑖−2,𝑗+2−𝑢𝑖,𝑗₂₃𝑢,⋮∇𝑖,𝑗=𝑢𝑖−2,𝑗+1−𝑢𝑖,𝑗₂₁₄𝑢,∇𝑖,𝑗=𝑢𝑖+2,𝑗+2−𝑢𝑖,𝑗₂₁₅𝑢𝑖,𝑗=𝑢𝑖+1,𝑗+2−𝑢𝑖,𝑗. (21)

Similar to previous discussion, the second-scale second-order 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 ¹⁵_𝑘=0∇²_𝑘𝑢𝑖,𝑗 (22) 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: (2𝑛+1)²−(2𝑛−1)²=8𝑛. (23)

Therefore, the 𝑛th scale enclosed gradient vector is defined as ∇𝑢𝑖,𝑗=∇₀𝑢𝑖,𝑗,∇₁𝑢𝑖,𝑗,…,∇_8𝑛−1𝑢𝑖,𝑗^𝑇, (24) where 𝑇 represents the transpose of the vector and ∇𝑢_𝑘(𝑖,𝑗), 𝑘=0,…,8𝑛−1, are defined as ∇₀𝑢𝑛𝑖,𝑗=𝑢𝑖,𝑗+2,∇−𝑢𝑖,𝑗₁𝑢𝑛𝑖,𝑗=𝑢𝑖−0.5,𝑗+2,∇−𝑢𝑖,𝑗₂𝑢𝑛𝑖,𝑗=𝑢𝑖−1,𝑗+2,∇−𝑢𝑖,𝑗₃𝑢𝑛𝑖,𝑗=𝑢𝑖−1.5,𝑗+2,⋮∇−𝑢𝑖,𝑗_8𝑛−2𝑢𝑛𝑖,𝑗=𝑢𝑖+1,𝑗+2,∇−𝑢𝑖,𝑗_8𝑛−1𝑢𝑛𝑖,𝑗=𝑢𝑖+0.5,𝑗+2−𝑢𝑖,𝑗. (25)

The 𝑛th scale second-order difference of 𝑢(𝑖,𝑗) is ∇²𝑢𝑖,𝑗=∇²₀𝑢𝑖,𝑗,∇²₁𝑢𝑖,𝑗,…,∇²_8𝑛−1𝑢𝑖,𝑗^𝑇, (26) where 𝑇 represents the transpose of the vector. We have ∇²₀𝑢,∇𝑖,𝑗=𝑢𝑖,𝑗+𝑛−𝑢𝑖,𝑗₂₁𝑢,∇𝑖,𝑗=𝑢𝑖−1,𝑗+𝑛−𝑢𝑖,𝑗₂₂𝑢,∇𝑖,𝑗=𝑢𝑖−2,𝑗+𝑛−𝑢𝑖,𝑗₂₃𝑢,⋮∇𝑖,𝑗=𝑢𝑖−3,𝑗+𝑛−𝑢𝑖,𝑗_28𝑛−2𝑢,∇𝑖,𝑗=𝑢𝑖+2,𝑗+𝑛−𝑢𝑖,𝑗_28𝑛−1𝑢𝑖,𝑗=𝑢𝑖+1,𝑗+𝑛−𝑢𝑖,𝑗. (27)

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 ^8𝑛−1_𝑘=0∇²_𝑘𝑢𝑖,𝑗 (28) is called the 𝑛th scale enclosed Laplacian operator (ELO) of 𝑈(𝑖,𝑗).

3.3. 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 𝑔=𝑔₀,𝑔₁,…,𝑔_(2𝑛+1)²₋₁^𝑇, (29) where 𝑛>2 is an integer, 𝑇 represents the transpose of the vector and 𝑔_𝑘, 𝑘=0,…,(2𝑛+1)²−1, is defined as 𝑔_𝑘=𝑔𝐺_𝜎1⋅∇_𝑘𝑢𝑖,𝑗∑^{(2𝑛+1)2−1}_𝑡=0𝑔𝐺_𝜎1⋅∇_𝑡𝑢,𝑖,𝑗 (30) where 𝐺_𝜎1 is defined in (11), ∇_𝑘𝑢(𝑖,𝑗), 𝑘=0,…,(2𝑛+1)²−1, defined in (24) are the components of vector ∇𝑢(𝑖,𝑗), ∑^{(2𝑛+1)2−1}_𝑡=0𝑔(|𝐺_𝜎1⋅∇_𝑡𝑢(𝑖,𝑗)|) is the normalized constant, and 𝑔 is the decreasing function of absolute value of ∇_𝑘𝑢(𝑖,𝑗), 𝑘=0,…,(2𝑛+1)²−1. Following (2) and (3), 𝑔(|𝐺_𝜎1⋅∇𝑢_𝑘(𝑥,𝑦,𝑡)|) can be defined as 𝑔𝐺_𝜎1⋅∇𝑢_𝑘𝑥,𝑦,𝑡=𝑒^{−(|𝐺𝜎1⋅∇𝑢𝑘(𝑥,𝑦,𝑡)|/𝜎)2} (31) or 𝑔𝐺_𝜎1⋅∇𝑢_𝑘=1𝑥,𝑦,𝑡1+𝐺_𝜎1⋅∇𝑢_𝑘𝑥,𝑦,𝑡/𝜎², (32) 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 𝑔𝐺_𝜎1⋅∇𝑢_𝑘𝑥,𝑦,𝑡=𝑒^{−(|𝐺𝜎1⋅∇𝑢𝑘(𝑥,𝑦,𝑡)|/𝜎)2}≈𝑒^{−(|𝐺𝜎1⋅∇2𝑢𝑘(𝑥,𝑦,𝑡)|/𝜎)2},𝑔𝐺_𝜎1⋅∇𝑢_𝑘=1𝑥,𝑦,𝑡1+𝐺_𝜎1⋅∇𝑢_𝑘𝑥,𝑦,𝑡/𝜎²≈11+𝐺_𝜎1⋅∇²𝑢_𝑘𝑥,𝑦,𝑡/𝜎². (33)(34)

The new AD based on the enclosed second-order difference is defined as 𝜕𝑢𝑖,𝑗,𝑡𝑔𝜕𝑡=div⁰∇₀𝑢𝑔𝑖,𝑗,𝑡₁∇₁𝑢𝑔𝑖,𝑗,𝑡₂∇₂𝑢⋮𝑔𝑖,𝑗,𝑡_{(2𝑛+1)2−3}∇_(2𝑛+1)²₋₃𝑢𝑔𝑖,𝑗,𝑡_{(2𝑛+1)2−2}∇_(2𝑛+1)²₋₂𝑢𝑔𝑖,𝑗,𝑡_{(2𝑛+1)2−1}∇_(2𝑛+1)²₋₁𝑢,𝑖,𝑗,𝑡 (35) where the ∇_𝑘𝑢(𝑖,𝑗,𝑡), 𝑘=0,…,(2𝑛+1)²−1, are the components of vector ∇𝑢(𝑖,𝑗,𝑡) in (24) and 𝑔_𝑘, 𝑘=0,…,(2𝑛+1)²−1, defined in (30) are the components of 𝑔 in (29). Moreover, the numerator of 𝑔_𝑘 in (30) can be approximated by (33) or (34).

The above equation can be represented as 𝜕𝑢𝑖,𝑗,𝑡=𝜕𝑡^{(2𝑛+1)2−1}_𝑘=0𝑔_𝑘∇²_𝑘𝑢𝑖,𝑗,𝑡, (36) where ∑^{(2𝑛+1)2−1}_𝑘=0𝑔_𝑘=1 and the right-hand side of (36) is the normalized weight ELO. ∇²_𝑘𝑢(𝑖,𝑗,𝑡) is the second-order difference of the 𝑘th components of 𝑢(𝑖,𝑗) which can be computed according to (26).

Thus the explicit form for solving (36) is 𝑢𝑖,𝑗,𝑡+1=𝑢𝑖,𝑗,𝑡+𝜆^{(2𝑛+1)2−1}_𝑘=0𝑔_𝑘∇²_𝑘𝑢𝑖,𝑗,𝑡, (37) where 𝑢(𝑖,𝑗,𝑡+1) is the gray level of (𝑖,𝑗) at time 𝑡+1 and 𝑔_𝑘, ∇²_𝑘𝑢(𝑖,𝑗,𝑡) are the same as in (36).

4. 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.

4.1. 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].

(a) LDCT image with tube current time product 30 mAs

(b) LDCT image with tube current time product 60 mAs

(c) LDCT image with tube current time product 30 mAs

(d) LDCT image with tube current time product 30 mAs

Figure 2 

LDCT images.

(a) Related SDCT image of Figure 2(a) with tube current time product 150 mAs

(b) Related SDCT image of Figures 2(b) and 2(c) with tube current time product 150 mAs

(c) Related SDCT image of Figure 2(d) with tube current time product 150 mAs

Figure 3 

Related SDCT images.

4.2. 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 function; finally, 1 is subtracted by the fuzzy edge function to get 𝛼(𝑥,𝑦,𝑡). Therefore, the diffusion equation becomes 𝜕𝑢𝑥,𝑦,𝑡𝛼𝑔.𝜕𝑡=div𝑥,𝑦,𝑡𝑥,𝑦,𝑡∇𝑢𝑥,𝑦,𝑡 (38)

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.

4.3. 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 𝜎₁ used for RPM, WT-RPM, and RMDMS-AD is also set to 1.5 since, in [3], the authors suggest that 𝜎₁ 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.

4.4. Results and Discussion

From 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.

(a) Recovered LDCT image with scale 𝑠=1 and iteration number 𝑡=3

(b) Recovered LDCT image with scale 𝑠=2 and iteration number 𝑡=3

(c) Recovered LDCT image with scale 𝑠=3 and iteration number 𝑡=3

(d) Recovered LDCT image with scale 𝑠=1 and iteration number 𝑡=5

(e) Recovered LDCT image with scale 𝑠=2 and iteration number 𝑡=5

(f) Recovered LDCT image with scale 𝑠=3 and iteration number 𝑡=5

(g) Recovered LDCT image with scale 𝑠=1 and iteration number 𝑡=10

(h) Recovered LDCT image with scale 𝑠=2 and iteration number 𝑡=10

(i) Recovered LDCT image with scale 𝑠=3 and iteration number 𝑡=10

(j) Recovered LDCT image with scale 𝑠=1 and iteration number 𝑡=20

(k) Recovered LDCT image with scale 𝑠=2 and iteration number 𝑡=20

(l) Recovered LDCT image with scale 𝑠=3 and iteration number 𝑡=30

Figure 4 

Recovered LDCT images of Figure 2(a) using RMDMS-AD with different scales 𝑠 and different iteration numbers 𝑡. The first column: 𝑠=1; the second column: 𝑠=2; the three column: 𝑠=3. The first row: 𝑡=3; the second row: 𝑡=5; the third row: 𝑡=10; the fourth row: 𝑡=30.

(a) Recovered LDCT image using RMDMS-AD with 𝑡=5

(b) Recovered LDCT image using RMDMS-AD with iteration number 𝑡=10

(c) Recovered LDCT image using RMDMS-AD with iteration number 𝑡=20

(d) Recovered LDCT image using RPM [2] with iteration number 𝑡=5

(e) Recovered LDCT image using RPM [2] with iteration number 𝑡=10

(f) Recovered LDCT image using RPM [2] with iteration number 𝑡=20

(g) Recovered LDCT image using PM [1] with iteration number 𝑡=5

(h) Recovered LDCT image using PM [1] with iteration number 𝑡=10

(i) Recovered LDCT image using PM [1] with iteration number 𝑡=20

(j) Recovered LDCT image using WE-PM [3] with iteration number 𝑡=5

(k) Recovered LDCT image using WE-PM [3] with iteration number 𝑡=10

(l) Recovered LDCT image using WE-PM [3] with iteration number 𝑡=20

(m) Recovered LDCT image using WT-PM [11] with iteration number 𝑡=5

(n) Recovered LDCT image using WT-PM [11] with iteration number 𝑡=10

(o) Recovered LDCT image using WT-PM [11] with iteration number 𝑡=20

Figure 5 

Recovered LDCT images of Figure 2(c) using RMDMS-AD (the first row), RPM (the second row [2]), PM (the third row [1]), WE-PM (the fourth row [3]), and WT-PM (the fifth row [11]) with different iteration numbers 𝑡=5 (the first column), 𝑡=10 (the second column), and 𝑡=20 (the third column). The other parameters include 𝜎1=1.5, 𝜎=2, and 𝜆=1/7 and scale 𝑠=3 for the proposed method.

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 greatly which remind us that “regularization” is an essential technology in AD.

5. 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.

Acknowledgments

This paper is supported by the National Natural Science Foundation of China (no. 60873102), Major State Basic Research Development Program (no. 2010CB732501), and Open Foundation of Visual Computing and Virtual Reality Key Laboratory of Sichuan Province (no. J2010N03). This work was supported by a Grant from the National High Technology Research and Development Program of China (no. 2009AA12Z140).