# A Numerical Method for Preserving Curve Edges in Nonlinear Anisotropic Smoothing

Authors: Shaoxiang Hu; Zhiwu Liao; Dan Sun; Wufan Chen
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

We focus on nonlinearity for images and propose a new method which can preserve curve edges in image smoothing using nonlinear anisotropic diffusion
(NAD). Unlike existing methods which diffuse only among the spatial variants, the new method suggests that the diffusion should be performed both among the time variants and spatial variants, named time and space nonlinear anisotropic diffusion (TSNAD). That is, not only the differences of the spatial variants should be estimated by the nearby spatial points but also the differences of the time variants should be approximated by the weighted time differences of nearby points, according to the differences of gray levels between them and the consideration point. Since the time differences of nearby points using NAD can find more points with similar gray levels which form a curve belt for the center pixel on a curve edge, TSNAD can provide satisfied smoothing results while preserving curve edges. The experiments for digital images also show us the ability of TSNAD to preserve curve edges.


---

## Body

## 1. Introduction

Nonlinear systems are systems that cannot be mathematically described as the sum of their components, and their mathematical modeling is often very difficult or impossible. As a result, nonlinear systems are often studied through use of simulations [1]. Nonlinearity for many real applications has been widely discussed [2–11], and some methods based on NAD have been proposed in the early 90s [12–15]. Recently, some efforts for improving the performance of NAS have been proposed related to numerical method for partial different function (PDF) theory, adaptive smoothing, and multiresolution analysis [16–22].Although NAD provides prefect mathematical theory for preserving edges, how to preserve curve edges in image smoothing using NAD is still an unsolved problem [23, 24].In a digital image, a curve edge is composed by discrete points instead of a continuous curve line. Moreover, the discrete points on a curve edge are composed by two types of points: one is edge and the other is corner. Note that most of curve edge points are corners instead of edge points since the curve or declining parts of a curve edge are composed by corners. Note that unlike the definition of corners in mathematics, a corner point defined in a digital image is a point with both large modules of gradients and tangents [12–14] or two large eigenvalues for the tensor matrix [15].In order to preserve curve edges, some schemes are proposed for different types of points. For edge points, two famous diffusion schemes, which are diffusion only along tangent lines or diffusion along backward directions of the gradients, are proposed in [14, 15]. For corners, Alvarez et al. suggest to stop diffusion at [13].However, these schemes lead to noise amplification, oscillation effects, and blurred curve edges [23, 24]. Analyzing behavior of these schemes, we can find that the diffusion for most of discrete points of a curve edge, corners, is stopped or only is carried on a very small scale. Therefore, the noise is amplified near edges, and corners are blurred even with a very small-scale diffusion.Motivated by diffusion along tangents, we propose that the diffusion should be performed along curve edges to preserve curve edges. However, until now, we can not find a method to perform the diffusion along the curve edges. It is the result of one-step diffusion along tangent! Note that diffusion along the tangent means that the diffusion is performed on a straight line passing through the consideration point and with the direction of the tangent. Therefore, even with an elaborately selected scale, a straight line can not approximate a curve edge.Fortunately, it is well known that a curve edge can be approximated well by multistraight lines. Although one-step method only gets one straight line, we have to design a multistep scheme to obtain multistraight lines. This multistep scheme can be realized when we try to approximate the differences of the time variant of the consideration point using the time differences of nearby points! Since each time difference relates to one straight diffusion line, the time differences of nearby points can find multistraight diffusion lines along the curve edge. Thus, it forms a curve belt along the curve edge!Based on this theory, a new scheme, named time and space NAD (TSNAD), is proposed that approximating the difference of time variant uses the time difference of nearby points with NAD manner. That is, the time differences of nearby points are weighted by the absolute value of the gray differences between the nearby points and the consideration point. Then, all time differences are computed by the space differences according to the diffusion equation.Since more time differences allow us to find more straight lines, TSNAD can find enough points with similar gray levels along the curve edges. Therefore, it can preserve curve edges well.The rest of this paper is as follows. Section2 is the method of TSNAD. In Section 3, the experimental results are given and discussed. We also give conclusions, future works, and finally acknowledgments.


## 2. The Method

Just as the above discussion, TSNAD provides a more NAD for time variant to preserve curve edges. In this section, we will introduce the theory of TSNAD in detail.


### 2.1. Motivation

The main objective for image smoothing is reducing image details or slightly noises without removing significant parts of the image content. Nonlinear anisotropic diffusion (NAD) is a well-known method in this task [12].Formally, NAD is defined as(2.1)∂u(x,y,t)∂t=div⁡(g(x,y,t)∇u(x,y,t)),
where u(x,y,0) is the initial gray-scale image, u(x,y,t) is the smooth gray-scale image at time t, ∇ denotes the gradient, div(·) is the divergence operator, and g(x,y,t) is the diffusion coefficient. g(x,y,t) controls the rate of diffusion and is usually chosen as a monotonically decreasing function of the module of the image gradient. Two functions proposed in [12] are(2.2)g(‖∇u(x,y,t)‖)=e-(‖∇u(x,y,t)‖/σ)2,(2.3)g(‖∇u(x,y,t)‖)=11+(‖∇u(x,y,t)‖/σ)2,
where ∥·∥ is the module of the vector, and the constant σ controls the sensitivity to edges.Perona and Malik also propose a simple method to approach the modules of gradients which is called PM method as  follows from the paper [12]. Its discretization for Laplacian operator is(2.4)u(i,j,t+1)=u(i,j,t)+14[cN⋅∇N2u(i,j,t)+cS⋅∇S2u(i,j,t)+cE⋅∇E2u(i,j,t)+cW⋅∇W2u(i,j,t)],
where(2.5)∇N2u(i,j,t)=u(i-1,j,t)-u(i,j,t),∇S2u(i,j,t)=u(i+1,j,t)-u(i,j,t),∇E2u(i,j,t)=u(i,j+1,t)-u(i,j,t),∇W2u(i,j,t)=u(i,j-1,t)-u(i,j,t).
According to (2.2) and (2.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(2.6)cN(i,j,t)=g(|∇Nu(i,j,t)|),cS(i,j,t)=g(|∇Su(i,j,t)|),cE(i,j,t)=g(|∇Eu(i,j,t)|),cW(i,j,t)=g(|∇Wu(i,j,t)|),
where |·| is the absolute value of the number and g(·) is defined in (2.2) or (2.3).Although this scheme is not the exact discretization of (2.1) from the view of mathematics, it is the best way for preserving singularities in existing NAD methods since it can find neighbors with similar gray levels correctly without computing the gradients. Thus, it reduces the discretization errors greatly.In addition, unlike other methods which have to distinguish corners form the edges, PM method allows to handle corners and edges using a unified scheme diffusion only among neighbors with similar gray levels. This advantage also decreases the influence of incorrect classification for corners and edges.However, even for this best scheme in existing methods, how to preserve curve edges in image smoothing is an unsolved problem. The main difficulty for PM method and other existing methods is that approximation of a curve edge using a straight line segment (a line passing through the consideration point and along the tangent with a certain scale) is impossible.Thus, we need a new scheme to find multistraight lines to approximate a curve edge. Intuitively, these multistraight lines should be related to some tangent lines passing through some of neighbors of the consideration points. Thus, these multistraight lines can be sought out using a quasimultisteps. That is, differences of the time variant are computed by the time differences of nearby pixels. It is obvious that each time difference relates to one tangent line. Thus, multitime differences relate to multitangent lines, which form a curve belt along the curve edge.By this way, the new scheme, TSNAD, can provide satisfied smoothing results for preserving curve edges and deleting details.


### 2.2. The Space Difference

Existing methods approximate the time difference for the left hand of (2.1) by(2.7)∂u(x,y,t)∂t=u(i,j,t+1)-u(i,j,t)Δt,
where u(i,j) is the consideration point, and t is the time variant.The space difference for a two-dimension digital imageu(i,j) is defined as half-point differences between the center point u(i,j) and the half points of its eight nearest neighbors(2.8)∇u(i,j)=[∇0u(i,j),∇1u(i,j),∇2u(i,j),∇3u(i,j),∇4u(i,j),∇5u(i,j),∇6u(i,j),∇7u(i,j)]T,
where T represents the transpose of the vector, and ∇uk(i,j), k=0,…,7 are defined as(2.9)∇0u(i,j)=u(i,j+0.5)-u(i,j),∇1u(i,j)=u(i-0.5,j+0.5)-u(i,j),∇2u(i,j)=u(i-0.5,j)-u(i,j),∇3u(i,j)=u(i-0.5,j-0.5)-u(i,j),∇4u(i,j)=u(i,j-0.5)-u(i,j),∇5u(i,j)=u(i+0.5,j-0.5)-u(i,j),∇6u(i,j)=u(i+0.5,j)-u(i,j),∇7u(i,j)=u(i+0.5,j+0.5)-u(i,j).Thus, the first-scale second-order difference ofu(i,j) is(2.10)∇2u(i,j)=[∇02u(i,j),∇12u(i,j),∇22u(i,j),∇32u(i,j),∇42u(i,j),∇52u(i,j),∇62u(i,j),∇72u(i,j)]T,
where T represents the transpose of the vector. From (2.7), we have(2.11)∇02u(i,j)=u(i,j+1)-u(i,j),∇12u(i,j)=u(i-1,j+1)-u(i,j),∇22u(i,j)=u(i-1,j)-u(i,j),∇32u(i,j)=u(i-1,j-1)-u(i,j),∇42u(i,j)=u(i,j-1)-u(i,j),∇52u(i,j)=u(i+1,j-1)-u(i,j),∇62u(i,j)=u(i+1,j)-u(i,j),∇72u(i,j)=u(i+1,j+1)-u(i,j).Let(2.12)g(i,j)=[g0(i,j),g1(i,j),g2(i,j),g3(i,j),g4(i,j),g5(i,j),g6(i,j),g7(i,j)]T,
where T represents the transpose of the vector, and gk(i,j), k=0,…,7 are defined as(2.13)gk(i,j)=g(|∇ku(i,j)|),k=0,1,…,7,
where ∇ku(i,j), k=0,…,7 defined in (2.8) are the components of vector ∇u(i,j), and g is the decreasing function of absolute value of ∇ku(i,j), k=0,…,7. Following (2.2) and (2.3), g(|∇uk(x,y,t)|) can be defined as(2.14)g(|∇uk(i,j,t)|)=e-(|∇uk(i,j,t)|/σ)2,k=0,…,7,
or(2.15)g(|∇uk(i,j,t)|)=11+(|∇uk(i,j,t)|/σ)2,k=0,…,7,
where |·| is the absolute value of the number, and the constant σ controls the sensitivity to edges. However, half-point difference for ∇u(i,j) defined in (2.8) cannot be computed directly. Thus, the second-order integral-point difference can be used to approximate the first-order half-point difference. The (2.14) and (2.15) become(2.16)g(|∇uk(i,j,t)|)=e-(|∇uk(i,j,t)|/σ)2≈e-(|∇2uk(i,j,t)|/σ)2,k=0,…,7,(2.17)g(|∇uk(i,j,t)|)=11+(|∇uk(i,j,t)|/σ)2≈11+(|∇2uk(i,j,t)|/σ)2,k=0,…,7.The NAD is converted to(2.18)∂u(i,j,t)∂t=λdiv⁡(g0(i,j)∇0u(i,j,t)g1(i,j)∇1u(i,j,t)g2(i,j)∇2u(i,j,t)g3(i,j)∇3u(i,j,t)g4(i,j)∇4u(i,j,t)g5(i,j)∇5u(i,j,t)g6(i,j)∇6u(i,j,t)g7(i,j)∇7u(i,j,t)),
where λ is a constant to ensure the convergence of the iteration, the ∇ku(i,j,t), k=0,…,7 are the components of vector ∇u(i,j,t) in (2.10), and gk(i,j), k=0,…,7 defined in (2.13) are the components of g(i,j) in (2.12). Moreover, gk(i,j) in (2.13) can be represented by (2.16) or (2.17).The above equation also is(2.19)∂u(i,j,t)∂t=λ∑k=07gk(i,j)∇k2u(i,j,t),
where λ is a constant to ensure the convergence of the iteration. ∇k2u(i,j,t) is the second-order difference of the kth components of u(i,j) which can be computed according to (2.10).Thus, we have(2.20)∂u(i,j,t)∂t=u(i,j,t+1)-u(i,j,t)=λ∑k=07gk(i,j)∇k2u(i,j,t),
where u(i,j,t+1) is the gray level of (i,j) at time t+1, and gk(i,j), ∇k2u(i,j,t), k=0,…,7 are defined in the same way as in (2.19).


### 2.3. The New Time Difference

In this paper, we propose that not only the space differences but also the time differences should be estimated using NAD. That is, the time difference ofu(i,j,t) should be estimated using the time differences of nearby points. In this section, we will discuss the simplest form for new time differences approximated by the weighted time differences of its eight nearest neighbors.Representing this approximation using the terms of high-dimensional vector analysis, a vector is defined as(2.21)∂u(i,j,t)∂t=[∂u0(i,j,t)∂t,∂u1(i,j,t)∂t,∂u2(i,j,t)∂t,∂u3(i,j,t)∂t,∂u4(i,j,t)∂t,∂u5(i,j,t)∂t,∂u6(i,j,t)∂t,∂u7(i,j,t)∂t]T,
where T represents the transpose of the vector, and ∂uk(i,j,t)/∂t, k=0,…,7 are defined as(2.22)∂u0(i,j,t)∂t=∂u(i,j+1,t)∂t,∂u1(i,j,t)∂t=∂u(i-1,j+1,t)∂t,∂u2(i,j,t)∂t=∂u(i-1,j,t)∂t,∂u3(i,j,t)∂t=∂u(i-1,j-1,t)∂t,∂u4(i,j,t)∂t=∂u(i,j-1,t)∂t,∂u5(i,j,t)∂t=∂u(i+1,j-1,t)∂t,∂u6(i,j,t)∂t=∂u(i+1,j,t)∂t,∂u7(i,j,t)∂t=∂u(i+1,j+1,t)∂t.The diffusion coefficients for∂uk(i,j,t)/∂t, k=0,…,7 are represented by a vector g(i,j) defined in (2.12), and its components are defined in (2.16). Thus, the time difference of u(i,j) is(2.23)∂u(i,j,t)∂t=λ1[∂u(i,j,t)∂t+div⁡(g(i,j)T∂u(i,j,t)∂t)],
where T represents the transpose of the vector, λ1 is a constant to ensure the convergence of the iteration, and ∂u(i,j,t)/∂t is defined in (2.21). Thus, (2.23) is(2.24)∂u(i,j,t)∂t=λ1[∂u(i,j,t)∂t+∑k=07gk(i,j)∂uk(i,j,t)∂t],
where ∂uk(i,j,t)/∂t, k=0,…,7 are the components of vector ∂u(i,j,t)/∂t defined in (2.21), and gk(i,j), k=0,…,7 are the components of vector g(i,j) defined in (2.12).According to (2.24), the time difference of u(i,j,t) is estimated by the weighed time difference of its neighbors according to the difference of the gray levels between them and u(i,j,t). That is, the neighbors whose gray levels are more similar to u(i,j,t) have larger weights.


### 2.4. Time and Space NAD (TSNAD)

Substituting (2.19) and (2.21) into (2.24), we have the equation of TSNAD,(2.25)∂u(i,j,t)∂t=λ1[∂u(i,j,t)∂t+∑k=07gk(i,j)∂uk(i,j,t)∂t]=λ1[λ∑k=07gk(i,j)∇k2u(i,j,t)+λg0(i,j)∑k=07gk(i,j+1,t)∇k2u(i,j+1,t)+λg1(i,j)∑k=07gk(i-1,j+1,t)∇k2u(i-1,j+1,t)+λg2(i,j)∑k=07gk(i-1,j,t)∇k2u(i-1,j,t)+λg3(i,j)∑k=07gk(i-1,j-1,t)∇k2u(i-1,j-1,t)+λg4(i,j)∑k=07gk(i,j-1,t)∇k2u(i,j-1,t)+λg5(i,j)∑k=07gk(i+1,j-1,t)∇k2u(i+1,j-1,t)+λg6(i,j)∑k=07gk(i+1,j,t)∇k2u(i+1,j,t)+λg7(i,j)∑k=07gk(i+1,j+1,t)∇k2u(i+1,j+1,t)].Letλ1×λ=λ be a constant to ensure the convergence of the iteration, (2.25) becomes
(2.26)∂u(i,j,t)∂t=λ[∑k=07gk(i,j)∇k2u(i,j,t)+g0(i,j)∑k=07gk(i,j+1,t)∇k2u(i,j+1,t)+g1(i,j)∑k=07gk(i-1,j+1,t)∇k2u(i-1,j+1,t)+g2(i,j)∑k=07gk(i-1,j,t)∇k2u(i-1,j,t)+g3(i,j)∑k=07gk(i-1,j-1,t)∇k2u(i-1,j-1,t)+g4(i,j)∑k=07gk(i,j-1,t)∇k2u(i,j-1,t)+g5(i,j)∑k=07gk(i+1,j-1,t)∇k2u(i+1,j-1,t)+g6(i,j)∑k=07gk(i+1,j,t)∇k2u(i+1,j,t)+g7(i,j)∑k=07gk(i+1,j+1,t)∇k2u(i+1,j+1,t)],
where gk(i,j), k=0,…,7 defined in (2.13) are the components of vector g(i,j) defined in (2.12), and ∇k2u(i,j), k=0,…,7 are the components of ∇2u(i,j) defined in (2.10).


## 2.1. Motivation

The main objective for image smoothing is reducing image details or slightly noises without removing significant parts of the image content. Nonlinear anisotropic diffusion (NAD) is a well-known method in this task [12].Formally, NAD is defined as(2.1)∂u(x,y,t)∂t=div⁡(g(x,y,t)∇u(x,y,t)),
where u(x,y,0) is the initial gray-scale image, u(x,y,t) is the smooth gray-scale image at time t, ∇ denotes the gradient, div(·) is the divergence operator, and g(x,y,t) is the diffusion coefficient. g(x,y,t) controls the rate of diffusion and is usually chosen as a monotonically decreasing function of the module of the image gradient. Two functions proposed in [12] are(2.2)g(‖∇u(x,y,t)‖)=e-(‖∇u(x,y,t)‖/σ)2,(2.3)g(‖∇u(x,y,t)‖)=11+(‖∇u(x,y,t)‖/σ)2,
where ∥·∥ is the module of the vector, and the constant σ controls the sensitivity to edges.Perona and Malik also propose a simple method to approach the modules of gradients which is called PM method as  follows from the paper [12]. Its discretization for Laplacian operator is(2.4)u(i,j,t+1)=u(i,j,t)+14[cN⋅∇N2u(i,j,t)+cS⋅∇S2u(i,j,t)+cE⋅∇E2u(i,j,t)+cW⋅∇W2u(i,j,t)],
where(2.5)∇N2u(i,j,t)=u(i-1,j,t)-u(i,j,t),∇S2u(i,j,t)=u(i+1,j,t)-u(i,j,t),∇E2u(i,j,t)=u(i,j+1,t)-u(i,j,t),∇W2u(i,j,t)=u(i,j-1,t)-u(i,j,t).
According to (2.2) and (2.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(2.6)cN(i,j,t)=g(|∇Nu(i,j,t)|),cS(i,j,t)=g(|∇Su(i,j,t)|),cE(i,j,t)=g(|∇Eu(i,j,t)|),cW(i,j,t)=g(|∇Wu(i,j,t)|),
where |·| is the absolute value of the number and g(·) is defined in (2.2) or (2.3).Although this scheme is not the exact discretization of (2.1) from the view of mathematics, it is the best way for preserving singularities in existing NAD methods since it can find neighbors with similar gray levels correctly without computing the gradients. Thus, it reduces the discretization errors greatly.In addition, unlike other methods which have to distinguish corners form the edges, PM method allows to handle corners and edges using a unified scheme diffusion only among neighbors with similar gray levels. This advantage also decreases the influence of incorrect classification for corners and edges.However, even for this best scheme in existing methods, how to preserve curve edges in image smoothing is an unsolved problem. The main difficulty for PM method and other existing methods is that approximation of a curve edge using a straight line segment (a line passing through the consideration point and along the tangent with a certain scale) is impossible.Thus, we need a new scheme to find multistraight lines to approximate a curve edge. Intuitively, these multistraight lines should be related to some tangent lines passing through some of neighbors of the consideration points. Thus, these multistraight lines can be sought out using a quasimultisteps. That is, differences of the time variant are computed by the time differences of nearby pixels. It is obvious that each time difference relates to one tangent line. Thus, multitime differences relate to multitangent lines, which form a curve belt along the curve edge.By this way, the new scheme, TSNAD, can provide satisfied smoothing results for preserving curve edges and deleting details.


## 2.2. The Space Difference

Existing methods approximate the time difference for the left hand of (2.1) by(2.7)∂u(x,y,t)∂t=u(i,j,t+1)-u(i,j,t)Δt,
where u(i,j) is the consideration point, and t is the time variant.The space difference for a two-dimension digital imageu(i,j) is defined as half-point differences between the center point u(i,j) and the half points of its eight nearest neighbors(2.8)∇u(i,j)=[∇0u(i,j),∇1u(i,j),∇2u(i,j),∇3u(i,j),∇4u(i,j),∇5u(i,j),∇6u(i,j),∇7u(i,j)]T,
where T represents the transpose of the vector, and ∇uk(i,j), k=0,…,7 are defined as(2.9)∇0u(i,j)=u(i,j+0.5)-u(i,j),∇1u(i,j)=u(i-0.5,j+0.5)-u(i,j),∇2u(i,j)=u(i-0.5,j)-u(i,j),∇3u(i,j)=u(i-0.5,j-0.5)-u(i,j),∇4u(i,j)=u(i,j-0.5)-u(i,j),∇5u(i,j)=u(i+0.5,j-0.5)-u(i,j),∇6u(i,j)=u(i+0.5,j)-u(i,j),∇7u(i,j)=u(i+0.5,j+0.5)-u(i,j).Thus, the first-scale second-order difference ofu(i,j) is(2.10)∇2u(i,j)=[∇02u(i,j),∇12u(i,j),∇22u(i,j),∇32u(i,j),∇42u(i,j),∇52u(i,j),∇62u(i,j),∇72u(i,j)]T,
where T represents the transpose of the vector. From (2.7), we have(2.11)∇02u(i,j)=u(i,j+1)-u(i,j),∇12u(i,j)=u(i-1,j+1)-u(i,j),∇22u(i,j)=u(i-1,j)-u(i,j),∇32u(i,j)=u(i-1,j-1)-u(i,j),∇42u(i,j)=u(i,j-1)-u(i,j),∇52u(i,j)=u(i+1,j-1)-u(i,j),∇62u(i,j)=u(i+1,j)-u(i,j),∇72u(i,j)=u(i+1,j+1)-u(i,j).Let(2.12)g(i,j)=[g0(i,j),g1(i,j),g2(i,j),g3(i,j),g4(i,j),g5(i,j),g6(i,j),g7(i,j)]T,
where T represents the transpose of the vector, and gk(i,j), k=0,…,7 are defined as(2.13)gk(i,j)=g(|∇ku(i,j)|),k=0,1,…,7,
where ∇ku(i,j), k=0,…,7 defined in (2.8) are the components of vector ∇u(i,j), and g is the decreasing function of absolute value of ∇ku(i,j), k=0,…,7. Following (2.2) and (2.3), g(|∇uk(x,y,t)|) can be defined as(2.14)g(|∇uk(i,j,t)|)=e-(|∇uk(i,j,t)|/σ)2,k=0,…,7,
or(2.15)g(|∇uk(i,j,t)|)=11+(|∇uk(i,j,t)|/σ)2,k=0,…,7,
where |·| is the absolute value of the number, and the constant σ controls the sensitivity to edges. However, half-point difference for ∇u(i,j) defined in (2.8) cannot be computed directly. Thus, the second-order integral-point difference can be used to approximate the first-order half-point difference. The (2.14) and (2.15) become(2.16)g(|∇uk(i,j,t)|)=e-(|∇uk(i,j,t)|/σ)2≈e-(|∇2uk(i,j,t)|/σ)2,k=0,…,7,(2.17)g(|∇uk(i,j,t)|)=11+(|∇uk(i,j,t)|/σ)2≈11+(|∇2uk(i,j,t)|/σ)2,k=0,…,7.The NAD is converted to(2.18)∂u(i,j,t)∂t=λdiv⁡(g0(i,j)∇0u(i,j,t)g1(i,j)∇1u(i,j,t)g2(i,j)∇2u(i,j,t)g3(i,j)∇3u(i,j,t)g4(i,j)∇4u(i,j,t)g5(i,j)∇5u(i,j,t)g6(i,j)∇6u(i,j,t)g7(i,j)∇7u(i,j,t)),
where λ is a constant to ensure the convergence of the iteration, the ∇ku(i,j,t), k=0,…,7 are the components of vector ∇u(i,j,t) in (2.10), and gk(i,j), k=0,…,7 defined in (2.13) are the components of g(i,j) in (2.12). Moreover, gk(i,j) in (2.13) can be represented by (2.16) or (2.17).The above equation also is(2.19)∂u(i,j,t)∂t=λ∑k=07gk(i,j)∇k2u(i,j,t),
where λ is a constant to ensure the convergence of the iteration. ∇k2u(i,j,t) is the second-order difference of the kth components of u(i,j) which can be computed according to (2.10).Thus, we have(2.20)∂u(i,j,t)∂t=u(i,j,t+1)-u(i,j,t)=λ∑k=07gk(i,j)∇k2u(i,j,t),
where u(i,j,t+1) is the gray level of (i,j) at time t+1, and gk(i,j), ∇k2u(i,j,t), k=0,…,7 are defined in the same way as in (2.19).


## 2.3. The New Time Difference

In this paper, we propose that not only the space differences but also the time differences should be estimated using NAD. That is, the time difference ofu(i,j,t) should be estimated using the time differences of nearby points. In this section, we will discuss the simplest form for new time differences approximated by the weighted time differences of its eight nearest neighbors.Representing this approximation using the terms of high-dimensional vector analysis, a vector is defined as(2.21)∂u(i,j,t)∂t=[∂u0(i,j,t)∂t,∂u1(i,j,t)∂t,∂u2(i,j,t)∂t,∂u3(i,j,t)∂t,∂u4(i,j,t)∂t,∂u5(i,j,t)∂t,∂u6(i,j,t)∂t,∂u7(i,j,t)∂t]T,
where T represents the transpose of the vector, and ∂uk(i,j,t)/∂t, k=0,…,7 are defined as(2.22)∂u0(i,j,t)∂t=∂u(i,j+1,t)∂t,∂u1(i,j,t)∂t=∂u(i-1,j+1,t)∂t,∂u2(i,j,t)∂t=∂u(i-1,j,t)∂t,∂u3(i,j,t)∂t=∂u(i-1,j-1,t)∂t,∂u4(i,j,t)∂t=∂u(i,j-1,t)∂t,∂u5(i,j,t)∂t=∂u(i+1,j-1,t)∂t,∂u6(i,j,t)∂t=∂u(i+1,j,t)∂t,∂u7(i,j,t)∂t=∂u(i+1,j+1,t)∂t.The diffusion coefficients for∂uk(i,j,t)/∂t, k=0,…,7 are represented by a vector g(i,j) defined in (2.12), and its components are defined in (2.16). Thus, the time difference of u(i,j) is(2.23)∂u(i,j,t)∂t=λ1[∂u(i,j,t)∂t+div⁡(g(i,j)T∂u(i,j,t)∂t)],
where T represents the transpose of the vector, λ1 is a constant to ensure the convergence of the iteration, and ∂u(i,j,t)/∂t is defined in (2.21). Thus, (2.23) is(2.24)∂u(i,j,t)∂t=λ1[∂u(i,j,t)∂t+∑k=07gk(i,j)∂uk(i,j,t)∂t],
where ∂uk(i,j,t)/∂t, k=0,…,7 are the components of vector ∂u(i,j,t)/∂t defined in (2.21), and gk(i,j), k=0,…,7 are the components of vector g(i,j) defined in (2.12).According to (2.24), the time difference of u(i,j,t) is estimated by the weighed time difference of its neighbors according to the difference of the gray levels between them and u(i,j,t). That is, the neighbors whose gray levels are more similar to u(i,j,t) have larger weights.


## 2.4. Time and Space NAD (TSNAD)

Substituting (2.19) and (2.21) into (2.24), we have the equation of TSNAD,(2.25)∂u(i,j,t)∂t=λ1[∂u(i,j,t)∂t+∑k=07gk(i,j)∂uk(i,j,t)∂t]=λ1[λ∑k=07gk(i,j)∇k2u(i,j,t)+λg0(i,j)∑k=07gk(i,j+1,t)∇k2u(i,j+1,t)+λg1(i,j)∑k=07gk(i-1,j+1,t)∇k2u(i-1,j+1,t)+λg2(i,j)∑k=07gk(i-1,j,t)∇k2u(i-1,j,t)+λg3(i,j)∑k=07gk(i-1,j-1,t)∇k2u(i-1,j-1,t)+λg4(i,j)∑k=07gk(i,j-1,t)∇k2u(i,j-1,t)+λg5(i,j)∑k=07gk(i+1,j-1,t)∇k2u(i+1,j-1,t)+λg6(i,j)∑k=07gk(i+1,j,t)∇k2u(i+1,j,t)+λg7(i,j)∑k=07gk(i+1,j+1,t)∇k2u(i+1,j+1,t)].Letλ1×λ=λ be a constant to ensure the convergence of the iteration, (2.25) becomes
(2.26)∂u(i,j,t)∂t=λ[∑k=07gk(i,j)∇k2u(i,j,t)+g0(i,j)∑k=07gk(i,j+1,t)∇k2u(i,j+1,t)+g1(i,j)∑k=07gk(i-1,j+1,t)∇k2u(i-1,j+1,t)+g2(i,j)∑k=07gk(i-1,j,t)∇k2u(i-1,j,t)+g3(i,j)∑k=07gk(i-1,j-1,t)∇k2u(i-1,j-1,t)+g4(i,j)∑k=07gk(i,j-1,t)∇k2u(i,j-1,t)+g5(i,j)∑k=07gk(i+1,j-1,t)∇k2u(i+1,j-1,t)+g6(i,j)∑k=07gk(i+1,j,t)∇k2u(i+1,j,t)+g7(i,j)∑k=07gk(i+1,j+1,t)∇k2u(i+1,j+1,t)],
where gk(i,j), k=0,…,7 defined in (2.13) are the components of vector g(i,j) defined in (2.12), and ∇k2u(i,j), k=0,…,7 are the components of ∇2u(i,j) defined in (2.10).


## 3. The Experiments and Discussion

In order to analyse the performance of TSNAD and compare it with existing methods, two images are selected in our experiments: one is a gray level image of an autotire, most of whose edges are curve edges (see Figure1(d)); the other is a binary image with some test patterns including circles with different width, filled circles and some shapes with curve edges, and so forth (see Figure 2(d)). Both images are selected from the standard test images of MatLab.Smoothing results for the noisy image (g) and the recoded difference images between the relative smoothing images at their left and the original image (d) using different schemes.
(a)(b)
PM method results(c)
Recoded difference image(d)
Original image(e)
Backward filter results(f)
Recoded difference image(g)
Noisy image(h)
Coherence filter results(i)
Recoded difference image(j)(k)
Proposed method results(l)
Recoded difference imageSmoothing results for the image (d) and the difference images between the relative smoothing images at their left and the original image (d) using different schemes.
(a)(b)
PM method results(c)
Difference image(d)
The image(e)
Backward filter results(f)
Difference image(g)(h)
Coherence filter results(i)
Difference image(j)(k)
Proposed method results(l)
Difference imageEssentially, there are only three schemes for NAD until now: PM method [12, 13], backward filter [14], and coherence filter [15].PM provides a simple discretization scheme for (2.1) by diffusion among points in four directions with similar gray levels to the consideration point. Although PM method is not the exact discretization of (2.1) from the view of mathematics, it has the best performance in existing methods because of the direct computing of diffusion coefficients using the differences of gray levels without computing the gradients [12].Backward filter suggests that the diffusion near edges should be performed along the backward directions of the gradients [14], while coherence filter suggests that the diffusion should be performed along the tangents [15]. These two methods have to approximate the gradients and tangents using few discrete directions which leads to discretization errors. Moreover, since these two methods assume that the curve edges are composed by both corners and edges which should be handled separately, the diffusion has to stop at the corners or perform in a very small scale. The former will amplify noises at the corners while the latter will blur the corners. Therefore, both schemes can not provide satisfied smoothing results for curve edges.TSNAD has all advantages for the above schemes. Firstly, its diffusion coefficients are computed according to PM method which can reduce the discretization errors greatly. Secondly, it adopts the scheme of diffusion along the curve edges to reduce the influence of the noises. Thirdly, it provides a more NAD for time variant to track the curve edges. Thus TSNAD has good performance both in image smoothing and curve edge preserving.In the first image (see Figure1(d)) is added a Gaussian white noise (GWN) with standard deviation 5 (see Figure 1(g)) to test the smoothing performance at slight noises. The NAD is performed on the noisy image, while the difference image is absolute difference value images between the relative smoothing images and the original image (see Figures 1(c), 1(f), 1(i), and 1(l)). Note that the difference images are  recoded to show the relation in difference images more clearly. Thus, although Figure 1(l), which is the difference image between the smoothing image using TSNAD and the original image, has the much more white points, it does not mean that TSNAD has the biggest absolute difference to the original image because of “recorded" image. However, the difference images at least can provide what parts of the original image represented by white points are corrupted more seriously. For example, the white points in Figures 1(f) and 1(i) are all on the curve edges means that backward filter and coherence filter will blur the curve edges, while the white points in Figure 1(l) are distributed very randomly means that the curve edges will be kept well.The second image is a binary image with many test patterns (see Figure2(d)). Observing the difference images in Figures 2(c), 2(f), 2(i), and 2(l), we can find that the difference image of the TSNAD has no white points means that it is very similar to the original image (see Figure 2(l)). Moreover, the difference image for the PM method has big differences at thinning edges (see Figure 2(c)), while the backward filter only lost curve edges (see Figure 2(f)). The difference image for the coherence filter shows that all edges for test patterns are lost (see Figure 2(i)).From the above experiments and discussion, we can conclude that TANAD can preserve curve edge well in image smoothing.


## 4. Conclusions

In this paper, we propose TSNAD to preserve curve edges in image smoothing. Unlike existing methods which only perform NAD for the space variant, TSNAD also allows us to approximate the differences of time variant using NAD. The more time differences for NAD form curve belts along the curve edges, which can provide better fitness for curve edges. Thus, TSNAD can get satisfied results both in curve edge preserving and image smoothing.


## 5. Future Works

Recently, fractional-order transforms become a hot topic in many fields both in theory and in application [5–7, 25–31] for their attractive natures in image and signal processing. Thus, our future works will be devoted to the combination of NAD and the fractal; it includes the following.(1)
It is a work trying the kernel of fractional order in (2.2) instead of power of 2 [32].(2)
The future work will consider the fractional Gaussian noise (fGn) following [25–27, 32, 33]. In addition to fGn, other classes of Gaussian noise, such as the generalized Cauchy (GC) process [29–31], are also worth experiments. Further, multiscaled fGn or the GC noise [28–31] are worth noting in experiments.



---

*Source: 186507-2011-05-10.xml*