Image denoising using non linear diffusion tensors

Image denoising is an important pre-processing step for many image analysis and computer vision system. It refers to the task of recovering a good estimate of the true image from a degraded observation without altering and changing useful structure in the image such as discontinuities and edges. In this paper, we propose a new approach for image denoising based on the combination of two non linear diffusion tensors. One allows diffusion along the orientation of greatest coherences, while the other allows diffusion along orthogonal directions. The idea is to track perfectly the local geometry of the degraded image and applying anisotropie diffusion mainly along the preferred structure direction. To illustrate the effective performance of our model, we present some experimental results on a test and real photographic color images.


INTRODUCTION
Image denoising has been one of the most important and widely studied problems in image processing and computer vision. The need to have a very good image quality is increasingly required with the advent of the new technologies in a various areas such as multimedia, medical image analysis, aerospace, video systems and others. Indeed, the acquired image is often marred by noise which may have a multiple origins such as: thermal fluctuations; quantify effects and properties of communication channels. It affects the perceptual quality of the image, decreasing not only the appreciation of the image but also the performance of the task for which the image has been intended. The challenge is to design methods, which can selectively smooth a degraded image without altering edges, losing significant features and producing reliable results. Traditionally, linear models have been commonly used to reduce noise. It is shown that these methods perform well in the flat regions of images but do not preserve edges and discontinuities which are often smeared out. In contrast, nonlinear models can handle edges in a much better way than those linear models. Many approaches have been proposed to remove the noise effectively while preserving the original image details and features as much as possible. In the past few years, the use of non linear PDEs methods involving anisotropic diffusion has significantly grown and becomes an important tool in contemporary image processing. The key idea behind the anisotropic diffusion is to incorporate an adaptative smoothness constraint in the denoising process. That is, the smooth is encouraged in a homogeneous region and discourage across boundaries, in order to preserve the natural edge of the image. One of the most successful tools for image denoising is the Total Variation (TV) model [9] [8] [10] and the anisotropic smoothing model [1] which has since been expanded and improved upon [3] [5] [20]. Over the years, other very interesting denoising methods have been emerged such as: Bilateral filter and its derivatives [6][11] [12]. In our work, we address image denoising problem by using the so-called structure tensors [14] which have proven their effectiveness in several areas such as: texture segmentation [17], motion analysis [18] and corner detection [16][19] [21]. The structure tensor provides a more powerful description of local pattern images better than a simple gradient. Based on its eigenvalues and the corresponding eigenvectors, the tensor summarizes the predominant directions of the gradient in a specified neighborhood of a point and the degree to which those directions are coherent. Our contribution in this work lies in the use of two coupled diffusion tensors which allow a complete and coherence regularization process which significantly improves the quality image. This paper is organized as follows. In Section 2, we introduce the non linear diffusion PDEs and discuss the various options that have been implemented for the anisotropic diffusion. In Section 3, we present the non linear structure tensor formalism and its mathematical concept. Section 4 focuses on our proposed denoising approach. Numerical experiments and results on test and real photographic images are shown in Section 5.

NON LINEAR DIFFUSION PDE: OVERVIEW
In this section, we review some basic mathematical concepts of the nonlinear diffusion PDEs proposed by Perona and Malik [1]. Let ‫ݑ‬ ሺ‫ݔ‬ǡ ‫ݕ‬ǡ ‫ݐ‬ሻǣ ߗ ՜ ܴ be the grayscaled intensity image with a diffusion time t, for the where ߲ ௧ ‫ݑ‬ǣ denotes the first derivative regarding the diffusion time t; ȁ‫ݑ‬ȁ: denotes the gradient modulus and g( ) is a non-increasing function, known as the diffusivity function which allow isotropic diffusion in flat regions and no diffusion near edges. By developing the divergence term of (1), we obtain: where: ‫ݑ‬ ఎఎ ൌ ߟ ୄ ‫ܪ‬ ߟ and ‫ݑ‬ కక ൌ ߦ ୄ ‫ܪ‬ ߦ are respectively the second spatial derivatives of ‫ݑ‬ in the directions of the gradient ߟ ൌ ‫,‪ȁ‬ݑ‪Ȁȁ‬ݑ‬ and its orthogonal ߦ ൌ ߟ ୄ ; H denotes the Hessian of u. According to these definitions, on the image discontinuities, we have the diffusion along ߟ (normal to the edge) weighted with ܿ ఎ ൌ ݃ ᇱᇱ ሺȁߘ‫ݑ‬ȁሻ and a diffusion along ߦ (tangential to the edge) weighted with ܿ క ൌ ݃ ᇱ ሺȁߘ‫ݑ‬ȁሻȀȁ‫ݑ‬ȁ. To understand the principle of the anisotropic diffusion, let represent a contour C (figure 1) separating two homogeneous regions of the image, the isophote lines (e.g. level curves of equal graylevels) correspond to u(x,y) = c. In this case, the vector ߟ is normal to the contour C, the set ሺߦǡ ߟሻ is then a moving orthonormal basis whose configuration depends on the current coordinate point (x, y). In the neighborhood of a contour C, the image presents a strong gradient. To better preserve these discontinuities, it is preferable to diffuse only in the direction parallel to C (i.e. in the ߦ -direction).
In this case, we have to inhibit the coefficient of ‫ݑ‬ ఎఎ (e.g. ܿ ఎ ൌ Ͳሻ, and to suppose that the coefficient of ‫ݑ‬ కక does not vanish.

Figure 1 Image contour and its moving orthonormal basis ሺߦǡ ߟሻ
So it appears the following conditions for the choice of the g( ) functions to be edge preserving : allows anisotropic diffusion to preserve discontinuities for the high gradient.
An extension of the nonlinear diffusion filtering to vector-valued image (e.g. color image) has been proposed [20] [3]. It evolves: ߗ ՜ ܴ under the diffusion equations: Where ‫ݑ‬ :denotes the i th component channels of ‫ݑ‬ Note that σ ୀଵ ‫ݑ‪ȁ‬‬ ȁ ଶ , represents the luminance function, which coupled all vectors channels taking the strong correlations among channels. However, this luminance function is not being able to detect isoluminance contours.

NON LINEAR DIFFUSION TENSOR
The non linear diffusion PDE seen previously does not give reliable information in the presence of flow-like structures (e.g. fingerprints). It would be desirable to rotate the flow towards the orientation of interesting features. This can be easily achieved by using the structure tensor, also referred to the second moment matrix. For a multivalued image, the structure tensor has the following form: With ‫ݑ‬ ఙ ൌ ‫ܭ‬ ఙ ‫כ‬ ‫ݑ‬ ൌ ‫ܭ‬ ఙ ‫כ‬ ሺ ǡ ሻ: the smoothed version of the gradient which is obtained by a convolution with a Gaussian kernel ‫ܭ‬ ఙ . The structure scale ߪ determines the size of the resulting flow-like patterns. Increasing ߪ gives an increased distance between the resulting flow lines. These new gradient features allow a more precise description of the local gradient characteristics. However, it is more convenient to use a smoothed version of ܵ ఙ , that is: Where ‫ܭ‬ ఘ : a Gaussian kernel with standard deviation ߩ. The integration scale ߩ averages orientation information. Therefore, it helps to stabilize the directional behavior of the filter. In particular, it is possible to close interrupted lines if ߩ is equal or larger than the gap size. In order to enhance coherent structures, the integration scale ߩ should be larger than the structure scale ߪ (݁‫ݔ‬ǣ ߩ ൌ ͵ߪሻ.In summary, the convolution with the Gaussian kernels ‫ܭ‬ ఙ ǡ ‫ܭ‬ ఘ , make the structure tensor measure more coherent. To go further into the formalism, the structure tensor ‫ܬ‬ ఘ can be written over its eigenvalues (ߣ ା , ߣ ି ) and eigenvectors ሺߠ ା ǡ ߠ ି ሻ , that is: The eigenvectors of ‫ܬ‬ ఘ give the preferred local orientations, and the corresponding eigenvalues denote the local contrast along these directions.
The eigenvector ߠ ା which is associated to the larger eigenvalue ߣ ା defines the direction of largest spatial change (i.e. the "gradient" direction). There are several ways to express the norm of the vector gradient to detect edges and corners; the most used is ܰ ൌ ඥߣ ା ߣ ି [20]. The eigenvalues ሺߣ ା ǡ ߣ ି ሻ are indeed well adapted to discriminate different geometric cases: -If ߣ ା ؆ ߣ ି ؆ Ͳ , the region doesn't contain any edges or corners. For this configuration, the variation norm ܰ should be low. -If ߣ ା ‫ب‬ ߣ ି there are a lot of vector variations.
The current point may be located on a vector edge. For this configuration, the variation norm ܰ should be high. -If ߣ ା ؆ ߣ ି ‫ب‬ Ͳ, There is a saddle point of the vector surface, which can possibly be a vector corner in the image In this case ܰ should be even higher than the case above. We note that for the case of the scalar image (e.g. n=1), ߣ ା ൌ ȁ‫ݑ‬ȁ ଶ ,ߠ ା ൌ ߟ ൌ ‫‪ȁ‬ݑ‪Ȁȁ‬ݑ‬ and ܰ ൌ ȁ‫ݑ‬ȁǤ Moreover, Weickert [3] proposed a non linear diffusion tensor by replacing the diffusivity function g() in (1) with a structure tensor, to create a truly anisotropic scheme, that is: Where ሺ ሻ is the diffusion tensor which is positive definite symmetric 2x2 matrix. This tensor posses the same eigenvectors ߠ ି , ߠ ା as the structure tensor ‫ܬ‬ ఘ and uses Ȝ 1 and Ȝ 2 to control the diffusion speeds in these two directions, that is : The Diffusion tensor D takes the form of an ellipsoid as is represented in figure 2. We note that for the scalar image from (1), the diffusion tensor is reduced to D ൌ ሺȁ‫‬ȁሻ ୢ ; where I d : Identity matrix.

Figure 2 Diffusion tensor 2D representation
In the coherence enhancing diffusion (CED) proposed by Weiker [4], the eigenvalues are assembled via: -In flat regions, we should have Ȝ ା = Ȝ ି ൌ Ͳ, and then Ȝ ଵ ൌ Ȝ ଶ ൌ ଵ ; D = ଵ I d where Id is the identity matrix. The tensor D is defined to be isotropic in these regions and takes the form of a circle of radius ଵ .
-Along image contours, we have ߣ ା ‫ب‬ ߣ ି ‫ب‬ Ͳ, and then ɉ ଶ ɉ ଵ Ͳ. The diffusion tensor D is then anisotropic, mainly directed by the smoothed direction Ʌ ି of the image isophotes͘ The idea from the CED approach is that broken boundaries of a single structure could be reconnected by allowing diffusion along the orientation of greatest coherence

PROPOSED APPROACH
Our idea is to combine two types of tensors: one allows diffusion along the orientation of greatest coherence, while the other allows diffusion along orthogonal directions. It is viewed as a regularization process. The proposed equation is as follows: The parameter Į can be viewed as a parameter of regularization which ensures the compromise between the two tensors. The tensor ଵ is derived from (11), that is ଵ = ൫‫ܬ‬ ఘ ൯ with the following diffusion weight functions: The tensor D 2 is constructed from the local coordinate system ሺߦǡ ߟሻ . For a scalar image, the tensor D 2 takes the form: (14) (ߦ ‫כ‬ ǡ ߟ ‫כ‬ ሻ : denote the local coordinate system which are the smoothed version of ሺߦǡ ߟሻ; with: Where ‫ܩ‬ ఙ : denotes a Gaussian Kernel; and ȁ‫ݑ‬ȁ ൌ ඥ‫ݑ‬ ‫ݔߪ‬ By developing (14), the tensor ‫ܦ‬ ଶ can be expressed by: ሺ݃ ଵ ǡ ݃ ଶ ሻ: denote the conductivity in the direction of the gradient and along the isophotes respectively. There are several choices for these conductivities, the wise choice is:݃ ଶ ሺȁ‫ݑ‬ȁሻ ൌ ݁ ିሺȁ‫‬௨ȁȀሻ మ and ݃ ଵ ሺȁ‫ݑ‬ȁሻ ൌ ߚ݃ ଶ ሺȁ‫ݑ‬ȁሻ, with Ͳ ൏ ߚ ൏ ͳ. This allows preserving and enhancing edges while smoothing within flat regions. Indeed, when the gradient modulus ȁ‫ݑ‬ȁ is high (e.g. edges, corners region), both ሺ݃ ଶǡ ݃ ଵ ሻ tends to 0, inhibiting the effect of diffusion . In contrast, when ȁ‫ݑ‬ȁ is low (e.g. flat region), ሺ݃ ଶǡ ݃ ଵ ሻ tends to two constants: 1 and ߚ respectively, the diffusion is isotropic in (ߦ ‫כ‬ ǡ ߟ ‫כ‬ ሻ directions. The extension of (16) to multivalued (e.g. color) images require the integration of the spectral components, that is: ; i=1..n Furthermore, equation (12) can be solved numerically using finite differences [2]. The time derivative μ ୲ at (i,j,t n ) is approximated by the forward difference μ ୲ ൌ ሺ ୧ ୬ାଵ െ ୧ ୬ ሻȀο, which leads to the iterative scheme: The main steps of the algorithm are: -Step 1: Parameters initialization: N, ߙ, ߪ, ߩ,k -Step 2: while n ܰ (iterations number) -Make discretization of the tensor components : ݆ ଵଵ ǡ ݆ ଵଶ ǡ ݆ ଶଵ ǡ ݆ ଶଶ of ‫ܬ‬ ఘ from eq.(5), by finite differences -Reconstruct the tensor ଵ using eq.(11) -Reconstruct the tensor ଶ using eq.(17) -Apply the iterative scheme of (18)  In the following section, we will give some experimental results.

EXPERIMENTAL RESULTS
We test the performance of the proposed algorithm with the famous "Lena" image sized 256x256 pixels figure 3-(a). A white gaussian noise with SNR=7.27 dB is added to the original image to obtain the noisy version showed in figure 3(b).The model's parameters are fixed to: The restored image in figure 3-f shows a significantly improvement: edges and discontinuities have been recovered and preserved with a good suppression of noise. Compared to the other methods, we note that the TV one (figure 3(c)) has similar performance. Nevertheless, our method seems to better preserve and enhance discontinuities. The CED method (figure 3(e)) is efficient preserving details, but introduces artifacts and streamlines in flat regions. The BF method is efficient removing noise, but loses details. To evaluate and quantify the quality image, we use two measures: the classic peak signal-to-noise-ratio (PSNR) and the mean structural similarity (MSSIM) index [15], which compares the structure of two images after subtracting luminance, and normalizing variance. The MSSIM approximates the perceived visual quality of an image better than PSNR. It takes values in [0,1] and increases as the quality increases. Table 1 confirms the effectiveness of our model with the highest score of MSSIM.  It can be observed that the restored image seems to be sharper, less noisy, and having a good edges preservation.

CONCLUSION
In this paper, we have proposed a new approach for image denoising based on the diffusion tensors. The idea is to combine two types of tensors: one allows diffusion along the orientation of greatest coherence, while the other allows diffusion along orthogonal directions. This offers a flexible and effective control on the diffusion process. Experimental results on test and real digital pictures are very promising and provide very good quality images in terms of noise reduction and discontinuities preservation. Future work will include automatic parameters estimation and computational models that can automatically predict perceptual image quality.