A scale-based forward-and-backward diffusion process for adaptive image enhancement and denoising

This work presents a scale-based forward-and-backward diffusion (SFABD) scheme. The main idea of this scheme is to perform local adaptive diffusion using local scale information. To this end, we propose a diffusivity function based on the Minimum Reliable Scale (MRS) of Elder and Zucker (IEEE Trans. Pattern Anal. Mach. Intell. 20(7), 699-716, 1998) to detect the details of local structures. The magnitude of the diffusion coefficient at each pixel is determined by taking into account the local property of the image through the scales. A scale-based variable weight is incorporated into the diffusivity function for balancing the forward and backward diffusion. Furthermore, as numerical scheme, we propose a modification of the Perona-Malik scheme (IEEE Trans. Pattern Anal. Mach. Intell. 12(7), 629-639, 1990) by incorporating edge orientations. The article describes the main principles of our method and illustrates image enhancement results on a set of standard images as well as simulated medical images, together with qualitative and quantitative comparisons with a variety of anisotropic diffusion schemes.

Different attributes such as noise, due to image acquisition, quantization, compression and transmission, blur or artefacts can influence the perceived quality of digital images [1], and requires post-processing such as image smoothing and sharpening steps for further image analysis including image segmentation, feature extraction, classification and recognition. In order to reduce noise while preserving spatial resolution, recent approaches use an adaptive approach by applying a combination of smoothing and enhancing filter to the image: image areas with little edges or sharpness are selectively smoothed while sharper image areas are instead processed with edge enhancing filters. Thus, the optimal strategy for noisy image enhancement is the combination of smoothing and sharpening that is adaptive to local structure of the image [2] with the aim of improving signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) [3–8] of the image.

Scale-space methods in image processing have proven to be fundamental tools for image diffusion and enhancement. The scale-space concept was first presented by Iijima [9–11] and became popular later on by the works of Witkin [12] and Koenderink [13]. The theory of linear scale-space supports edge detection and localization, while suppressing noise by tracking features across multiple scales [12–17]. In fact, the linear scale-space is equivalent to a linear heat diffusion equation [13, 14]. However, this equation was found to be problematic as edge features are smeared and distorted after a few iterations. In order to overcome this problem, Perona and Malik [18] proposed an anisotropic diffusion partial differential equation (PDE), with a spatially constant diffusion coefficient. In their work, the term "anisotropic" refers to the case where the diffusivity is a scalar function varying with the location, it limits the smoothing of an image near pixels with a large gradient magnitude, which are essentially the edge pixels. Perona and Malik's work paved the way for a variety of anisotropic diffusive filtering methods [19–49], which have attempted to overcome the drawbacks and limitations of the original model and has produced significant advances. The main motivation for anisotropic diffusion is to reduce noise while minimizing image blurring across boundaries, but this process has several drawbacks, among them the disappearance of fine structures in low SNR or CNR regions and increased blurring in fuzzy boundaries. This is mainly due to the fact that the image gradient magnitude, on which the diffusion coefficient is estimated, is noisy and makes it difficult to distinguish between significant features and noise due to overlocalization, hence deciding edginess using the diffusion coefficient could be unreliable. In addition, traditional nonlinear diffusion filtering process does not offer any image-dependent guidance for selecting the optimum gradient magnitude at which the diffusion flow must have a maximum value [50]. Moreover, as it was expressed by Black et al. [29], the choice of the diffusion coefficients could greatly affect the level of preservation of the edges.

In this article, based on early works on forward-and-backward (FAB) diffusion schemes [38, 50], where the smoothing and sharpening actions are combined in the same diffusion process system, we propose a scale-based forward-and-backward diffusion (SFABD) scheme. The main idea of the proposed scheme is that the magnitude of the diffusion coefficient at each pixel is determined by taking into account the local property of the image through the scales. This is performed using the notion of the Minimum Reliable Scale (MRS) as proposed by Elder and Zucker [18]. This technique is based on statistical reliability of the edge detection operator responses at different scales [51]. The reliable scale as defined by Elder and Zucker, means that at the selected scale and larger ones, the likelihood of error due to sensor noise is below a standard tolerance. By choosing the MRS, for edge detection at each pixel in the image, we prevent errors due to sensor noise while simultaneously minimizing errors due to interference from nearby structure. Such a multiscale measure along with the gradient can capture more accurately edges over a wide range of blur and contrasts. Using this concept, a MRS-based diffusivity function is proposed. As a result, the proposed scheme can adaptively encourage strong smoothing in homogeneous regions, while suitable sharpening in detail and edge regions. Furthermore, we modify the Perona-Malik [50] discrete scheme by taking edge orientations into account.

The remainder of this article is organized as follows: Section 2 gives an overview of the state-of-the-art anisotropic diffusion filtering; Sect. 3 presents the proposed SFABD algorithm; In Sect. 4, we illustrate image enhancement results on a set of well known test images as well as artificial medical images, and perform a qualitative and quantitative comparison of our method with a variety of anisotropic diffusion schemes. Finally, Sect. 5 states our concluding remarks.

Perona and Malik [50] formulated anisotropic diffusion filtering as a process that encourages intraregional smoothing, while inhibiting interregional denoising. The Perona-Malik (P-M) nonlinear diffusion equation is of the form:

(1)

where ∇ is the gradient operator, div is the divergence operator and c(·) is the diffusion coefficient, which is a non-negative monotonically decreasing function of local spatial gradient. If c(·) is constant, then isotropic diffusion is enacted. In this case, all locations in the image, including the edges, are equally smoothed. This is an undesirable effect because the process cannot maintain the natural boundaries of objects. The P-M discrete version of Equation 1 is given by:

(2)

where (x, y) denotes the coordinates of a pixel to be smoothed in the 2-D image domain, t denotes the discrete time step (iterations). The constant λ is a scalar that determines the rate of diffusion, η(x, y) represents the neighbourhood of pixel (x, y) and |η(x, y)| is the number of neighbours of pixel (x, y). indicates the image intensity difference between two pixels at (x, y) and (p, q) to approximate the image gradient. For the 4-connected neighbourhood's case, the directional gradients are calculated in the following way:

(3)

In Perona-Malik's work [50], the diffusivity function has been defined as: , where α > 0

or

(4)

where is the gradient magnitude and the parameter k serves as a gradient threshold: a smaller gradient is diffused and positions with larger gradient are treated as edges. The P-M equation has several practical and theoretical drawbacks. As mentioned by Alvarez et al. [20], the continuous P-M equation is ill posed with the diffusion coefficients in (4); very close pictures can produce divergent solutions and therefore very different edges. This is caused by the fact that the diffusion coefficient c used in [50] leads to flux decreasing for some gradient magnitudes and the scheme may work locally as the inverse diffusion that is known to be ill posed, and can develop singularities of any order in arbitrarily small time. As a result, a large gradient magnitude no longer represents true edges and the diffusion coefficients are not reliable, resulting in unsatisfactory enhancement performance.

So far, much research has been devoted for improving the Perona-Malik's anisotropic diffusion method. For example, Catte et al. [19] showed that the P-M equation can be made well posed by smoothing isotropically the image with a scaling parameter σ, before computing the image gradient used by the diffusion coefficient:

(5)

Equation 5 corresponds to the regularized version of the P-M PDE, and I _σ = G _σ (I)*I is a smoothed version of I obtained by convolving the image with a zero-mean Gaussian kernel G _σ of variance σ². Similarly, Torkamani-Azar et al. [52] replaced the Gaussian filter with a symmetric exponential filter and the diffusion coefficient is calculated from the convolved image. Although these improvements can convert the ill-posed problem [53] in the Perona-Malik's anisotropic diffusion method into a well-posed one, their reported enhancement and denoising performance has been further improved. Weickert [54] proposed later a nonlinear diffusion coefficient that produces segmentation-like results given by:

(6)

where segmentation-like results are obtained using m = 4 and C₄ = 3.31488.

Black et al. [29] proposed a more robust "edge-stopping" function derived from Tukey's biweight:

(7)

where σ _e is the "robust scale". Anisotropic smoothing with such edge stopping function can preserve sharper boundaries than previous schemes. Another diffusivity function, based on sigmoid function, has been proposed by Monteil and Beghdadi [33]:

(8)

where γ controls the steepness of the min-max transition region of anisotropic diffusion, and k controls the extent of the diffusion region in terms of gradient gray-level.

Notice that all of anisotropic diffusion filters presented above, utilize a scalar-valued diffusion coefficient (diffusivity function) c that is adapted to the underlying image structure, Weickert [26, 30, 55] defined this "pseudoanisotropy" as "isotropic nonlinear", and pointed out that the consequence of isotropic nonlinear diffusion is that only the magnitude, but not the direction of the diffusion flux can be controlled at each image location. Noise close to edge features remains unfiltered due to the small flux in the vicinity of edges. To enable smoothing parallel to edges, Weickert [30] proposed edge enhancing diffusion by constructing the diffusion tensor based on an orientation estimate obtained from observing the image at a larger scale:

(9)

where D is a 2 × 2 diffusion tensor and is specially designed from the spectral elements of the structure tensor to anisotropically smooth the image, while taking into account its intrinsic local geometry, preserving its global discontinuities.

For simultaneously enhance, sharpen and denoise images, Gilboa et al. [38] proposed a FAB adaptive diffusion process, denoted here as GSZFABD, where a negative diffusion coefficient is incorporated into image-sharpening and enhancement processes to deblur and enhance the extremes of the initial signal:

(10)

where: k _f has similar role as the k parameter in the P-M diffusion equation; k _b and w define the range of backward diffusion, and are determined by the value of the gradient that is emphasized; α controls the ratio between the forward and backward diffusion; and the exponent parameters (n, m) are chosen as (n = 4, m = 1). Equation 10 is locally adjusted according to image features, such as edges, textures and moments. The GSZFABD diffusion process can therefore enhance features while locally denoising the smoother segments of images. Following the same ideas, and in order to avoid that the transition length between the maximum and minimum coefficient values varies with the gradient threshold, which makes controlling diffusion difficult, we proposed in [44] the local variance controlled forward-and-backward diffusion (LVCFABD) coefficient:

(11)

where β₁ and β₂ control the steepness for the min-max transition region of forward diffusion and backward diffusion, respectively. These two parameters are vital to the FAB diffusion behaviour and the transition width from isotropic to oriented flux can be altered by modulating them. In addition, the obtained diffusion process can preserve the transition length from isotropic to oriented flux, and thus it is better at controlling the diffusion behaviour than the FAB diffusion of Gilboa et al. [38].

In this article, we propose a SFABD scheme combining the forward-backward scheme given by Equation 10 and the notion of MRS as proposed by Elder and Zucker [18]. The MRS allows defining a classification map R(x, y), where each pixel (x, y) is classified into homogenous, edge or detail pixel. R(x, y) is then used in the coefficient α of Equation 10 to locally adapt the anisotropic diffusion. Finally, for implementing the SFABD scheme, we propose a modification of the P-M scheme by taking edge orientation into account.

3.1 Local scale-based classification map

In anisotropic diffusion scheme the rate of diffusion at each pixel is determined by diffusion coefficients that are monotonically decreasing functions of the gradient, thereby mainly ensuring strong smoothing in flat areas and weak diffusing near edge features. Thus, the strategy of identifying homogeneous and edge regions is very significant. Gradient is widely used to detect variable boundary in image processing, however, it is difficult for this measure to distinguish significant discontinuities from noise due to overlocalization. In addition, during anisotropic diffusion process, fine structures often disappear and increasing blurring occurs in fuzzy boundaries. To overcome this problem, we follow the idea of Elder and Zucker [18] of multiscale approach for edge detection, and explore the selection of proper scales for local estimation that depends upon the local structure of edges. The estimated scale is then used as a critical value and represents the MRS for each pixel in an image. The MRS proposed by Elder and Zucker [18] is based on two assumptions: (1) the noise comes from a stationary, zero-mean white noise process; (2) the standard deviation of the noise can be estimated from the image itself or by calibration. For the sake of clarity, the MRS is briefly described.

In edge detection, it is very important to estimate the nonzero gradient at each pixel in an image. The gradient computation from discrete data is an ill-posed problem. Smoothing the data with a Gaussian filter is the well-known regularization method. Hence, the gradient can be estimated using steerable Gaussian first derivative basis filters:

(12)

(13)

where σ₁ denotes the scale of the first derivative Gaussian kernel g(x, y, σ₁). The magnitude and direction of the gradient in an image I(x, y) are given by:

(14)

where:

(15)

θ is the gradient vector direction at (x, y). I _x (x, y, σ₁) and I _y (x, y, σ₁) are defined as:

(16)

(17)

In gradient-based edge detection, the local gradients in a homogeneous region due to noise will have a nonzero value. Thus, the likelihood that the response of the gradient operator induced by noise should be respected. Considering the derivative operation as a random process transformation, the probability distribution function (PDF) of a noise gradient can be represented as [56, 57]:

(18)

(19)

where the L² norm of the first derivative operator is , σ _n is the standard deviation of the sensor noise, and σ₁ is the scale of the Gaussian kernel. Elder and Zucker [18] defined the MRS as the unique scale at which the events can be reliably detected. By reliable, they mean that at this and larger scales, the likelihood of error due to sensor noise is equal to or below a predetermined false positive rate. Reliability is defined in terms of a Type I (false positive) error, α _I , for the entire image and a point-wise Type I error, α _p . Under the assumption of i.i.d. noise, the point-wise Type I error α _p can be computed from the probability of having no false positive edges as follows [18]:

(20)

where N is the total number of pixels in the image. By selecting the MRS, the error due to sensor noise is limited while the interference of neighbourhood structures is minimized. Given the probability distribution function (pdf) of gradient of the noise in equation (18), point-wise Type I error α _p is defined when using a gradient threshold c₁ to detect an edge:

(21)

Using the above equation, and considering a fixed type I error, we can define a critical value function:

(22)

Giving a point-wise Type I error α _p , c₁(σ₁) represents the statistically reliable minimum gradient response based on the sensor noise and operator scale. Figure 1 depicts the critical value function c₁(σ₁) for different noise levels and different Type I error rates. It is easy to observe that c₁(·) is a non-negative monotonically decreasing function of σ₁, which is helpful in detecting blurred boundaries. Comparing Figure 1a and 1b, we notice that c₁(σ₁) is more sensitive to the standard deviation of sensor noise σ _n than to the Type I error α _I . Furthermore, c₁(σ₁) grows as σ _n increases, for eliminating spurious edges in the presence of highly corrupted images. In this article, a thin-plate smoothing spline model is used to smooth a given image. It is assumed that the model whose generalized cross-validation score is minimum can provide the variance of the sensor noise σ _n [58].

Plots of the critical value function for different parameters settings. (a) The critical value function with respect to different noise levels (α _I = 0.05). (b) The critical value function with respect to different Type I error rates (σ _n = 20).

Full size image

For the MRS algorithm, how to sample the scale space is an open question. In scale space theory and for natural images, it is known that logarithmic scale is sufficient to represent the scale space completely [13]. For example, Elder and Zucker [18] used six logarithmic scales σ₁ = {0.5, 1, 2, 4, 8, 16} in their experiments. Table 1 summarises the alternative sampling schemes for scale space, both the Logarithmic and Limited-Log methods are logarithmic scales, while the latter has a shorter coverage. The Linear method samples the scale uniformly, and the Linear-Log one is a combination of Linear and Logarithmic. In this work, we empirically found the following linear sampling gives good results: σ₁ = {0.6, 0.9, 1.2, 1.5,..., 2.4}. In our implementation, we select the MRS at each pixel as the smallest scale at which the gradient estimate exceeds the critical value function:

(23)

Strictly speaking, if , the pixel usually resides in homogeneous regions and the MRS can be defined as , while , the pixel may correspond to edge or detailed feature and the MRS is chosen as .

Finally, we define the scale-based classification map R(x, y) as follows:

(24)

where R(x, y) denotes the region type of pixel (x, y). It has to be noted that the proper modulation of the thresholds σ_smooth and σ_edge is required for a robust classification map. As an example, the classification map of the Cameraman image and its noisy version (σ² = 225) are illustrated in Figures 2b and 2d, respectively. In the map, black regions are homogeneous, gray regions represent detail regions, while white regions manifest inhomogeneities that indicate most of the important features, e.g., the camera and tripod. This example clearly illustrates that the scale-based classification map readily indicate locations of highly homogeneous, detail and edge regions.

Local scale-based classification map of the Cameraman image. (a) Original image. (b) Classification map of (a) (σ_smooth = 2, σ_edge = 1). (c) Noisy image (σ² = 225). (d) Classification map of (c) (σ_smooth = 2, σ_edge = 1).

Full size image

3.2 Scale-based forward-and-backward diffusion coefficient

As shown in [59, 60], if the weight α in (10-11) is constant for all pixels in an image, this diffusion coefficient (10) is positive for small image gradients, while it becomes negative for large ones, and finally becomes positive again. Different nonlinear diffusion behaviours can be obtained by varying the value of α. For example, when α is large, the backward diffusion force is dominant. The larger α is, the more sharpening occurs. However, a too large α easily results in oscillations. Conversely, when α is small, the diffusion process performs image smoothing. Small values of α increase the noise attenuation at the price of a possible increase of detail blur. Thus, the optimal choice depends on the amount of noise variance and is typically a trade-off between smoothing and sharpening. Nevertheless, Gilboa et al. [38] proposed that the same α is used for the entire image, regardless of local structures of an image, which leads to over-smoothing in edge or detail regions and under-smoothing in homogeneous regions. In this article, we propose the balancing weight α, with different values δ_smooth, δ_edge and δ_detail , selectively applied at each pixel following the value of the local scale-based classification map R(x, y). Indeed, in the edge areas, the image requires more sharpening to highlight important features embedded in it, while in the detail regions, the diffusion process should compromise the effects of smoothing and sharpening, ensuring that the backward force can emphasize the fine structures in the image and the stabilizing forward force is strong enough to avoid oscillations. Whereas, in homogeneous regions, the gradient magnitude is very slow after the Gaussian pre-smoothing is applied to reduce the background noise. Thus, the approximate isotropic diffusion is performed to uniformly smooth the gentle and flat areas. In this way, each pixel is adaptively assigned a different parameter by evaluating the local structures. This is made possible using the MRS-based diffusivity function:

(25)

with

(26)

where δ_smooth, δ_edge and δ_detail are the scale-based weights, selected empirically such that δ_edge ≥ δ_detail ≥ δ_smooth ≥ 0. K_f and k_b control the MRS-based diffusivity function for forward and backward diffusion, respectively. As mentioned above, the parameter k_f has the same role as the gradient threshold in the P-M diffusion equation. Thus, the mean of local intensity differences in homogeneous regions is effective for controlling the forward diffusion; while k_b is determined by the value of the gradient that is emphasized. Previous works [38, 59] demonstrated that k_b is several times larger than k_f, in our case, we empirically defined the two gradient thresholds in (25) as [k_f, k_b] = [1, 2]*k. This strategy is indeed required in cases of natural signals or images because of their nonstationary structure. Usually, a minimal value of forward diffusion should be kept, so that large smooth areas do not become noisy. For the estimation of k, we use the assumption of i.i.d. noise, indeed, since the noise is assumed to be randomly distributed in the image space, a practical way of estimating its variance is to consider homogeneous regions where small variations or textures are mainly due to noise. Thus, k is estimated as the mean of the local intensity differences on the homogeneity map, i.e.,

(27)

where and B _xy is the neighbourhood set of pixel (x, y), and N _h is the total number of pixels in the homogeneous regions as defined by the classification map R(x, y). When Ω _h is empty, the simplest idea might be to setup k as a user defined constant, or using a "noise estimator": a histogram of the absolute values of the gradient throughout the image is computed, and k is set greater than or equal to e.g. 90% value of its integral at each iteration.

3.3 Edge orientation driven discretization scheme (EODDS)

As mentioned in Sect. 3.1, three different regions are classified before diffusion evolution. However, edge orientation is not taken care in the discrete scheme of P-M anisotropic diffusion. As a result, they are always considered to be displaced vertically or horizontally [61]. Moreover, one cannot recognize whether a slight intensity variation is mainly due to a slow varying edge or noise, so it is unreasonable that both situations are treated in the same way. The anisotropic diffusion discrete scheme should be modified to take edge orientations into account in the detail and edge regions, i.e. filtering action should be rather stronger on the direction parallel to the edge, and weaker on the perpendicular direction. Hence, we discretized the original anisotropic diffusion equation as follows:

(28)

where the mnemonic subscripts N, S, E, W, NE, SW, NW and SE denote the eight directions North, South, East, West, North-East, South-West, North-West and South-East, and the symbol ∇ stands for nearest-neighbour differences. λ is the time step for the numerical scheme; θ is the edge direction at pixel (x, y), W_V(θ), W_H(θ), and are weights for different edge directions.

For a nonlinear diffusion scheme, stability is an important issue that concerns possible unbounded growth or boundness of the final result of the diffusion scheme. The essential criterion defining stability is that the numerical process must restrict the amplification of all components from the initial conditions. In the following, we describe how to find the maximum value of λ assuring the stability. Assuming N _d the dimension of the neighborhood in direction d (in the vertical or horizontal direction for 4-connected neighbourhood, N _d = 1), the stability condition is given by [30]:

where D is the dimension of a given image. For our case (2-D images and 8-connected neighborhood), the condition becomes:

In this article, the step of keypoint orientation in scale-invariant feature transform (or SIFT) [62] algorithm is used for estimating the edge direction. The image is subdivided into nonoverlapping blocks of the same size, typically between 8 × 8 and 32 × 32 pixels. The gradient-based edge orientation histogram is then calculated in each block. If we let N be the total number of pixels in the image and n be the total number of bins, the histogram H _i meets the following conditions:

(29)

In the histogram, 360 degree is grouped in 36 groups, each of which is π/18 degree, and we obtain n = 36. Thus, the main orientation in each block is defined as follows:

(30)

and

(31)

where ϑ is the main gradient direction, by calculating the histogram of the gradient direction for each pixel (x, y) in the block, and "arctan" is the inverse tangent function. We assume that if an intensity variation between two zones is present, the edge has to be located along the perpendicular direction. The calculation of orientation histogram can be performed in real time. Furthermore, the comparison of orientation histograms can be performed using Euclidian distance that is very fast to compute for vectors whose dimensions are 36.

Once the estimation of the edge direction has been performed, the weights W_v(θ), W_H(θ), and have to be defined, in such a way that they satisfy the following constraint, with the aim of maintaining the numerical stability of the process:

(32)

In order to illustrate the way the weights are estimated, we divide the x - y plane into five domains as follows (see Figure 3):

Relating edge direction to direction in an image.

Full size image

(33)

Taking the constraint (32) and the trigonometric relation into account, the weights W_v(θ), W_H(θ), and are estimated as:

(34)

(35)

(36)

(37)

For instance, if θ ∈ Ω₀, substituting these weights in the modified anisotropic diffusion Equation 29 leads to the following:

(38)

In this case, the edge orientation should approximate the vertical direction according to the fact that the edge direction is always perpendicular to the gradient direction. During the diffusion process, a relatively large weight cos²θ is assigned in the vertical direction to guarantee that the diffusion should mainly occur in the direction parallel to the edge, while a relatively small weight sin²θ is assigned in the horizontal direction to ultimately avoid diffusion across the edge.

3.4 SFABD algorithm

The algorithm for the proposed SFABD scheme is summarised in Algorithm 1.

Chen [63] classified the existing performance evaluation methods into three categories; i.e. subjective, objective and application-based methodologies. By the subjective methodology, a noisy image and its enhanced images are illustrated. Thus, the evaluation on the performance of an algorithm is dependent on human's common sense gained from very much sophisticated visual perception experience. By the objective methodology, an evaluation is performed by comparing the enhanced image and its original uncorrupted version to see how much noise has been removed from a noisy image. By the application-based methodology, images in a certain application field are used for test and the enhancing results are assessed by a specialist who has expertise in the field or a comparison with an anticipated result set up prior to the test.

To assess the proposed approach, we follow the above-described methodology and demonstrate the effectiveness of SFABD in enhancing fine edge structures, i.e. we applied it to a variety of blurred and noisy images by comparing its results to five counterparts, namely, the Catte's anisotropic diffusion (CAD) [19], the robust anisotropic diffusion (RAD) [29], the Monteil's anisotropic diffusion (MAD) [33], the Weickert's anisotropic diffusion (WAD) [54], and the edge-enhancing diffusion (EED) [30]. The gradient threshold k should be chosen according to the noise level and the edge strength. In our experiments, we set k in different diffusion algorithms by referring to the original papers. The ultimate goal of image enhancement is to facilitate the subsequent processing for early vision. To demonstrate the usefulness of our algorithm in an early vision task, we apply our algorithm for performing edge-enhancing filtering on medical images, for an application-based evaluation.

In order to objectively evaluate the performance of the different diffusion algorithms, we adopt two noise-reduction measures: peak signal-to-noise ratio (PSNR) and the universal image quality index (UIQI). The measure of PSNR has been widely used in evaluating performance of a smoothing algorithm in the objective methodology. For a given noisy image I, I(i, j, T) denotes the intensity of pixel (x, y) ∈ I at iteration T while an anisotropic diffusion algorithm is applied to the noisy image. G(i, j) is its uncorrupted ground-truth. As a result, the PSNR is defined as follows:

(39)

Here, MAX _I is the maximum gray value of the image. When the pixels are represented using 8 bits per sample, MAX _I = 225. Typical values for the PSNR in lossy image and video compression are between 30 and 50 dB, where higher is better. Acceptable values for wireless transmission quality loss are considered to be about 20 to 25 dB [64, 65]. Recently, the UIQI has been used to better evaluate image quality due to its strong ability in measuring structural distortion occurred during the image degradation processes [66]:

(40)

where M is the total step number and Q _j denotes the local quality index computed within a sliding window. In this article, a sliding window of size 8 × 8 is applied to estimate an entire image. The dynamic range of Q is [-1,1], the value 1 is only achieved if the compared images are identical and the value of -1 means lowest quality of the distorted image.

4.1 General images

The performance of the proposed algorithm is evaluated using four standard images of size 512 × 512 and 256 gray-scale values. The image of Peppers is employed as an example of piecewise-constant image. The Lena and Cameraman images are two examples with both textures and smooth regions. The Boat image is an example with different edge features. For performance evaluation, the images have been corrupted with additive Gaussian white noise with different noise levels. The PSNR and UIQI values of the four noisy images with respect to different noise variance are listed in Table 2. The Lena and Boat images and their noisy versions with noise variance 225 are displayed in Figures 4 and 5, respectively. For clarity, only selected regions of the images are displayed.

Lena image. (a) Original image. (b) Noisy image with a noise variance of 225.

Full size image

Boat image. (a) Original image. (b) Noisy image with a noise variance of 225.

Full size image

Figures 6 and 7 depict the restored images using the six algorithms, for visual quality assessment. The results yielded by CAD and WAD schemes are depicted in Figures 6a, b and 7a, b, respectively. Both methods can well clean noise but blur the details of the restored results, such as the hat, its decoration and the hair in the image of Lena (see Figure 4a)), and the ground texture at the end of the Boat image of (see Figure 5a). This conforms our analysis that using the gradient, as only local discontinuity measure, would yield difficulties in distinguishing between edge details and noise and detecting fine structure. For RAD, a lot of noise still survives in the restored images. The restored results indicate that this method is very sensitive to noise. In Figures 6d and 7d, very large oscillations of gradient introduced by noise cannot be fully attenuated by MAD. The two resultant images present insufficient diffusion for restoration, in which the homogeneous background, such as Lena's face and bare shoulder (see Figure 4a) and the sky in the Boat image (see Figure 5a), cannot be completely eliminated because the diffusion process is terminated in early iterations. A better edge-preserving filtering is yielded by the EED process and the corresponding results are shown in Figures 6e and 7e, respectively. Finally, the images produced by the proposed SFABD scheme are represented in Figures 6f and 7f, respectively. The noise is removed and this is due to the forward diffusion. Meanwhile, edge features, including most of the fine details, are sharply reproduced. By comparing the resultant images of SFABD with the other five classical algorithms, we can notice that the SFABD algorithm achieves better visual quality. The reason for this is twofold: First, the multiscale discontinuity measure of the MRS-based diffusivity function is more effective than the gradient in detecting edge features and fine structure under a noisy environment, which is helpful for correctly classifying regions and estimating the gradient thresholds. Second, the proposed diffusion method incorporates both of the two discontinuity measures in the FAB diffusion coefficient by adopting a scale-based weight for balancing the forward diffusion and backward one. This strategy can ensure the elegant property of effectively smoothing noise while simultaneously sharpening edges and fine details of a noisy image. Table 3 lists the PSNR and UIQI values that are reported by the different algorithms, applied on the test images with different noise levels. For clarify, a noise variance of 400 is used for comparison. The experimental results demonstrate that the SFABD scheme can efficiently improve the PSNR value by around 8.6 dB better than the other algorithms. Additionally, the proposed diffusion scheme can produce an image with around 22% less structural distortion according to the UIQI values, which is the best among the six algorithms. Thus, we can say that the SFABD scheme outperforms the state-of-the-art diffusion methods. In addition, the performance of the EODDS has also been revealed in Table 4. It is evident that our algorithm using EODDS has achieved better statistical results than that of our algorithm without it, which confirms the validity of the EODDS.

Enhanced Lena image. (a) CAD. (b) WAD. (c) RAD. (d) MAD (γ = 0.1). (e) EED. (f) SFABD (σ = 0.1, σ_smooth = 2, σ_edge = 1, δ_smooth = 0.3, δ_edge = 0.6, δ_edge = 0.9) (10 iterations).

Full size image

Enhanced Boat image. (a) CAD. (b) WAD. (c) RAD. (d) MAD (γ = 0.1). (e) EED. (f) SFABD (σ = 0.1, σ_smooth = 2, σ_edge = 1, δ_smooth = 0.3, δ_detail = 0.6, δ_edge = 0.9) (10 iterations).

Full size image

Second, the proposed SFABD algorithm has been also compared to three existing FAB diffusion schemes, namely the GSZFABD [38], LVCFABD [44] and tunable FAB diffusion (TFABD) [47], using visual quality and the PSNR and UIQI values. Figure 8 depicts the obtained results of the considered FAB diffusion schemes. One can notice that all the four FAB processes can achieve a good compromise between sharpening and denoising. However, as illustrated in Figure 8a, the GSZFABD process blurs edges and detail features. From Figures 8b, c, it can be seen that the LVCFABD and TFABD schemes are sensitive to noise: the LVCFABD results in developing singularities in homogeneous regions, such as the inner parts of peppers, while the TFABD causes oscillations in the vicinity of edges, e.g. the outer contour of peppers. However, the proposed SFABD scheme exhibits the best edge-enhancing diffusion behaviour. The quantitative results of the four schemes are given in Table 4. It is evident from Table 4 that the SFABD scheme is much more efficient than the other three schemes for the four images. Hence, we can say that SFABD outperforms the existing FAB enhancement techniques.

Peppers image. (a) Original image. (b) Noisy image with a noise variance of 225. (c) Enhanced with GSZFABD. (d) LVCFABD. (e) TFABD. (f) SFABD. (σ = 0.1, σ_smooth = 2, σ_edge = 1, δ_smooth = 0.3, δ_detail = 0.6, δ_edge = 0.9) (10 iterations).

Full size image

In order to appraise the effectiveness of the adaptive gradient threshold, the gradient threshold k _f curves for four noisy images (σ² = 400) are depicted in Figure 9. It can be seen that all the curves, representing the evolution of this parameter, share the same decreasing behaviour as already demonstrated in other works, allowing lesser and lesser gradients to take part in the diffusing process. Moreover, after 20 iterations, k _f decreases slower and slower and the scheme converges to a steady state where for t → ∞, we get c(|∇I|) → 0, which means that almost no diffusion is performing. Note that, the estimation of an optimum threshold value k has been addressed by several authors [29, 50, 67, 68]. However, to our knowledge, these authors do not explain how to determine the homogenous regions during the process. In this work, an appropriate solution for automatically adapting the gradient threshold at each iteration has been proposed.

Gradient threshold evolution curves for the noisy test images Peppers, Lena, Cameraman and Boat (noise variance 400).

Full size image

4.2 Medical images

In medical images, low SNR and CNR often degrade the information and affect several image processing tasks, such as segmentation, classification and registration. Therefore, it is of considerable interest to improve SNR and CNR to reduce the deterioration of image information. In this section, we report the results of the proposed SFABD scheme on two three-dimensional MR images [69, 70], both of which have been simulated using two sequences (T1- and T2-weighted) with 1 mm of slice thickness, 9% noise level and 20% of intensity non-uniformity downloaded from Brainweb [71] using default acquisition parameters for each modality. These simulations are based on an anatomical model of normal brain, which can serve as ground truth for any analysis procedure.

Figures 10 and 11 show two examples of enhanced MR image using different diffusion schemes. The original noise-free images and their corrupted versions are illustrated in Figures 10a, b and 11a, b, respectively. As expected, the six algorithms remove noise and simultaneously smooth the homogeneous regions, such as white matter. However, for RAD, noise is still remaining in the resulting images. Some structure details are not visible in the images restored by the CAD, WAD and MAD algorithms, though they can greatly attenuate the effect of noise. According to the visual analyses of the image quality, the results given by the EED diffusion and the proposed SFABD are comparable, because the two processes perform edge-enhancing diffusion. Nevertheless, the SFABD scheme achieves better contrast and produces more reliable edges, which is especially useful for segmentation and classification purposes necessary in medical image applications.

Enhanced images for the 3-dimensional data volume of a T1-weighted MR simulated image. (a) Original MR image (slice 80). (b) Corrupted MR image. (c) Enhanced MR image with CAD. (d) WAD. (e) RAD. (f) MAD (γ = 0.1). (g) EED. (h) SFABD (σ = 0.1, σ_smooth = 2, σ_edge = 1, δ_smooth = 0.2, δ_detail = 0.4, δ_edge = 0.6) (10 iterations).

Full size image

Enhanced images for the 3-dimensional data volume of a T2-weighted MR simulated image. (a) Original MR image (slice 80). (b) Corrupted MR image. (c) Enhanced MR image with CAD. (d) WAD. (e) RAD. (f) MAD (γ = 0.1). (g) EED. (h) SFABD (σ = 0.1, σ_smooth = 2, σ_edge = 1, δ_smooth = 0.2, δ_detail = 0.4, δ_edge = 0.6) (10 iterations).

Full size image

In order to objectively evaluate the performances of the different diffusion algorithms on medical images, we adopt the PSNR and the Structural Similarity (SSIM) measure [72]. SSIM is a quality metric that measures the presence of the image structure details in the restored images and the value one is only achieved if the compared images are identical. The lowest value is zero if the images show no similarity at all. Since both the considered MR simulated images are three-dimensional data volume, we compare the PSNR and SSIM values at each slice for objective evaluation. As shown in Figure 12, the PSNR values of the restored images achieved by the proposed SFABD scheme are comparable or higher than the other diffusion algorithms, and the SSIM values of SFABD are significantly higher. Finally, the proposed scheme enhances boundary sharpness and fine structures better than other considered diffusion methods.

The PSNR and SSIM measures for the different diffusion algorithms at each slice of the T1- and T2-weighted MR images.

Full size image

We have presented a novel SFABD scheme for image restoration and enhancement. In the proposed scheme, the magnitude of the diffusion coefficients at each pixel is determined by taking into account the property of the image through scale-space, using a classification map obtained via the MRS. According to the type of the considered pixel (belonging to a homogenous, detail or edge region), a variable weight is incorporated into the anisotropic diffusion PDE to adaptively encourage strong smoothing in homogeneous regions and suitable sharpening in detail and edge regions. Moreover, we propose a method to estimate the parameter k of MRS-based diffusivity function, as the mean of the local intensity differences on homogeneous regions as determined by the MRS-based classification map. Finally, a numerical scheme, taking into account the edge orientation has been proposed. Furthermore, extensive qualitative and quantitative comparisons with a variety of existing diffusion schemes demonstrate the effectiveness of the proposed scheme, along with its potential use for medical image applications. Future work will involve two main aspect of the proposed approach, namely an adaptive approach for the estimation of the parameters, as well as establishing an automatic stopping criterion to replace the prefixed numbers of iteration for anisotropic diffusion.

1. 1.

   Initialize the image data I. I (x, y, 0) denotes the original intensity of pixel (x, y).

2. 2.

   Initialize the diffusion parameters. Set the values of the of the noise scale σ, the maximum number of iterations T, the classification map thresholds σ _smooth and σ _edge, and the scale-based weights δ _smooth, δ _edge and δ _detail.

3. 3.

   Calculate the critical value for each pixel and determine its region type.

4. a.

   Obtain the regularized image I _σ .

5. b.

   Compute the gradient of the smoothed image, ∇I _σ = (d _x ,d _y ) ^T .

6. c.

   Calculate the critical value for each pixel.

7. d.

   Determine the minimum reliable scale of each pixel by using the relationship between the spatial gradient and the critical value (23).

8. e.

   Estimate the classification map R(x, y) in (24)

9. 4.

   Iterate the diffusion filtering until t = T.

10. a.

    The gradient thresholds k _f and k _b are estimated as discussed in Section 3.2.

11. b.

    For each pixel (x, y), the diffusion coefficient c(∇) is computed using Eq. (25). In homogeneous and detail regions, the traditional 4-connected neighbourhood diffusion discretization equation is performed to update I(x, y, t); while in edge regions, the 8-connected neighbourhood diffusion discretization equation (28) is performed to update I(x, y, t).

CAD:

Catte's anisotropic diffusion

CNR:

contrast-to-noise ratio

EED:

edge-enhancing diffusion

EODDS:

edge orientation driven discretization scheme

FAB:

forward-and-backward

LVCFABD:

local variance controlled forward-and-backward diffusion

MAD:

Monteil's anisotropic diffusion

MRS:

Minimum Reliable Scale

PDE:

partial differential equation

pdf:

probability distribution function

P-M:

Perona-Malik

PSNR:

peak signal-to-noise ratio

RAD:

robust anisotropic diffusion

SIFT:

scale-invariant feature transform

SFABD:

scale-based forward-and-backward diffusion

SNR:

signal-to-noise ratio

SSIM:

Structural Similarity

TFABD:

tunable FAB diffusion

UIQI:

universal image quality index

WAD:

Weickert's anisotropic diffusion.

This work was supported in part by the National Natural Science Foundation of China under Grant 40901205, in part by the National Basic Research Program of China (973) under Grant 2009CB723905, in part by the Special Fund for Basic Scientific Research of Central Colleges, China University of Geosciences, Wuhan, under Grant CUGL090210, in part by the Foundation of Key Laboratory of Geo-informatics of State Bureau of Surveying and Mapping under Grant 201022, in part by the Foundation of Key Laboratory of Resources Remote Sensing & Digital Agriculture, Ministry of Agriculture under Grant RDA1005, in part by the Foundation of Key Laboratory of Education Ministry for Image Processing and Intelligent Control under Grant 200908, in part by the Foundation of Digital Land Key Laboratory of Jiangxi Province under Grant DLLJ201004. The authors would also like to thank the anonymous reviewers for their valuable comments and suggestions which significantly improved the quality of this article.

Authors and Affiliations

1. Institute of Geophysics and Geomatics, China University of Geosciences, People's Republic of China

   Yi Wang, Ruiqing Niu & Ke Wu

2. State Key Laboratory of Information Engineering in Surveying, Mapping, and Remote Sensing, Wuhan University, People's Republic of China

   Liangpei Zhang

3. Department of Electronics & Informatics (ETRO), Vrije Universiteit Brussel (VUB), Belgium

   Hichem Sahli

Corresponding author

Correspondence to Yi Wang.

Competing interests

The authors declare that they have no competing interests.

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions