# A Novel Multiscale Edge Detection Approach Based on Nonsubsampled Contourlet Transform and Edge Tracking

Authors: Enqing Chen; Jianbo Wang; Lin Qi; Weijun Lv
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Edge detection is a fundamental task in many computer vision applications. In this paper, we propose a novel multiscale edge detection approach based on the nonsubsampled contourlet transform (NSCT): a fully shift-invariant, multiscale, and multidirection transform. Indeed, unlike traditional wavelets, contourlets have the ability to fully capture directional and other geometrical features for images with edges. Firstly, compute the NSCT of the input image. Secondly, theK-means clustering algorithm is applied to each level of the NSCT for distinguishing noises from edges. Thirdly, we select the edge point candidates of the input image by identifying the NSCT modulus maximum at each scale. Finally, the edge tracking algorithm from coarser to finer is proposed to improve robustness against spurious responses and accuracy in the location of the edges. Experimental results show that the proposed method achieves better edge detection performance compared with the typical methods. Furthermore, the proposed method also works well for noisy images.


---

## Body

## 1. Introduction

Edges are prominent features in images and their detection and analysis are an essential task in computer vision and image processing. However, edge detection is a very difficult task, mainly because of the abstract nature of the edge detection problem. For instance, there is no clear definition of what an edge is. Typically, edges are defined as those points of an image where the gradient is noticeably large [1]. When viewing an image, humans can easily determine the boundaries within that image without needing to do so consciously. However, no single edge detection approach, up to the present, has been devised which will successfully determine every different type of edges.The conventional edge detectors include different differential operators, such as Sobel operator [2], Prewitt operator [3], and Robert operator [4]. These operators are easy to implement and, however, sensitive to noise due to lack of smoothing stage. They are inappropriate for the real images. Canny [5] puts forward well known three criteria of edge detectors. The criteria expressed by Canny are as follows [5, 6].(i)
Good detection (C1): there should be a low probability of failing to mark real edge points and low probability of falsely marking nonedge points.(ii)
Good localization (C2): the points marked as edge points by the operator should be as close as possible to the center of the true edge.(iii)
Only one response to a single edge (C3): each edge in the image should produce a single response.Besides, Canny also proposed a scheme for combining the outputs from different scales. His strategy is fine-to-coarse and this method is called feature synthesis. It was observed in [7, 8] that, at a single scale, the Canny edge detector is equivalent to the detection of the local maxima of the wavelet transform of an image, for some particular choices of the analyzing wavelet. However, the edges are still obtained by globally threshold processing. In this kind of approaches, many edges are extracted accompanying noises when the threshold value is too low; fewer significant edges remain when the threshold value is too high [9]. A richer description can be obtained by considering the response of the image to a family of filters of different scales and orientations.In the past, many multiscale methods based on wavelets have been successfully applied to the analysis and detection of edges [7–12]. For instance, Mallat et al. [7, 8] have implemented multiscale Canny edge detection with the discrete wavelet transform and proved that the wavelet transform modulus maxima detect all the singularities of a function. Shih and Tseng [9] proposed a two-stage edge extraction approach with contextual filter edge detector and multiscale edge tracker for further analysis of computer vision. Heric and Zazula [10] presented an edge detection algorithm using Haar wavelet transform and signal registration. Brannock and Weeks [11] have proposed an edge detection method based on the discrete wavelet transform (DWT), which combines DWT with other methods to achieve an optimal solution to the edge detection problem. Jiang et al. [12] propose a set of simplified version of Gabor wavelets and an efficient algorithm for extracting the features for edge detection. For one-dimensional (1D) piecewise smooth signals, wavelets have provided an optimal representation for these signals in a certain sense. However, natural images contain intrinsic geometrical structures that are key features in visual information. As a result of a separable extension from 1D bases, wavelets in two-dimension (2D) are not very effective in representing images containing distributed discontinuities such as edges. In addition, separable wavelets have a limited capability in dealing with directional information: an important and unique feature of multidimensional signals. These disappointing behaviors indicate that edge detection based on wavelet methods cannot get satisfactory results.In recent years, Do and Vetterli [13] proposed an efficient directional multiresolution image representation method, namely, contourlet transform. Indeed, unlike wavelet, the contourlet transform has good direction sensitivity, the anisotropy, and can catch accurately the image edge information. Owing to the geometric information, the contourlet transform achieves better results than wavelet transform in image analysis applications such as enhancement and denoising. Due to downsampling and upsampling, the contourlet transform lacks shift-invariance, which is desirable in many image applications such as edge detection, contour characterization, and image enhancement [14]. da Cunha et al. [15] propose an overcomplete transform called the nonsubsampled contourlet transform, which is a fully shift-invariant, multiscale, and multidirection expansion that has better directional frequency localization and a fast implementation.In this paper, a novel multiscale edge detection method based on the NSCT is proposed. Our edge detection approach has similarities with a number of other methods based on curvelet, complex wavelet, and shearlet. By contrast, the proposed approach is very competitive for detecting both the location and orientation of edges and can effectively reduce the noise influence and produce proper edges for images.The rest of this paper is organized as follows. Section2 introduces NSCT and K-means clustering algorithm. The proposed algorithm is described in Section 3. Section 4 shows the experimental results to demonstrate the effectiveness of the proposed method. Section 5 concludes this paper.


## 2. NSCT andK-Means Clustering

In this section, we will briefly review the theory of the NSCT and theK-means clustering, presented in [14, 15, 17] and in [18, 19], respectively. They will be applied in the subsequent sections.


### 2.1. NSCT

The NSCT is a shift-invariant version of the contourlet transform. The contourlet transform employs Laplacian pyramids (LPs) for multiscale decomposition and directional filter banks (DFBs) for directional decomposition [20]. Due to downsamplers and upsamplers present in both the LPs and DFBs, the contourlet transform is not shift-invariant. To achieve the shift-invariance, the NSCT is built upon nonsubsampled pyramids and nonsubsampled DFBs as shown in Figure 1.(a) NSFB structure that implements the NSCT. (b) Idealised frequency partitioning obtained with NSFB structure.
(a)
(b)(1) Nonsubsampled Pyramid (NSP). The nonsubsampled pyramid structure ensures the multiscale property of the NSCT. The building block of the nonsubsampled pyramid is two-channel nonsubsampled 2D filter banks. Figure 2(a) illustrates the nonsubsampled pyramid decomposition with J=3 stages [15]. The perfect reconstruction condition is given as
(1)H0zG0z+H1zG1z=1,
where H0(z) is low pass decomposition filter, H1(z) is high pass decomposition filter, G0(z) is low pass reconstruction filter, and G1(z) is high pass reconstruction filter.Nonsubsampled pyramid. (a) Three-stage pyramid decomposition. (b) Frequency divisions of a nonsubsampled pyramid.
(a)
(b)The ideal frequency support of the low pass filter at thejth stage is the region [-π/2j,π/2j]2. Accordingly, the support of the high pass filter is the complement of the low pass support region on the [-π/2j-1,π/2j-1]2∖[-π/2j,π/2j]2. The equivalent filters of a J level cascading NSPFB are given by
(2)Hnz=H1z2n-1∏j=0n-2H0z2j1≤n<2k∏j=0n-1H0z2jn=2k,
where H0(z) and H1(z) denote the low pass filter and the corresponding high pass filter at the first stage, respectively. The ideal frequency response of the nonsubsampled pyramid is given in Figure 2(b).(
2) Nonsubsampled Directional Filter Bank (NSDFB). The nonsubsampled DFB structure gives varying directions. The building block of the nonsubsampled DFB is also two-channel nonsubsampled 2D filter banks. The NSDFB is constructed by eliminating the downsamplers and upsamplers in the DFB [15]. To achieve multidirection decomposition, the NSDFB is iteratively computed. Figure 3 illustrates a four-channel directional decomposition with two-channel fan filter banks. The equivalent filter in each channel is given by
(3)Ukz=UizUjzQ.Four-channel NSDFB constructed with two-channel fan filter bank. (a) Filtering structure. (b) Corresponding frequency partitioning.
(a)
(b)(
3) Combining the Nonsubsampled Pyramid and Nonsubsampled Directional Filter Bank in the NSCT. The NSCT is constructed by combining the NSP and the NSDFB as shown in Figure 1(a). The NSP provides multiscale decomposition and captures the point discontinuities. The NSDFB provides directional decomposition and links point discontinuities into linear structures [21]. This scheme can be iterated repeatedly on the low pass subband outputs of nonsubsampled pyramids. The NSCT is very suitable for image edge detection because it has such important properties as multiresolution, localization, shift-invariance, and multidirection. In this paper, the “maxflat” filters and the “dmaxflat7” filters are, respectively, selected for NSPs and NSDFBs. The details of design for filters are referred to [15].


### 2.2.K-Means Clustering

TheK-means is a classical unsupervised clustering algorithm [22]. As an unsupervised technique, K-means clustering is effective in producing clusters for many practical applications such as image segmentation and image compression and determining the boundaries of the objects in an image.The aim ofK-means algorithm is to divide N points in D-dimensional space into K clusters so that the within-cluster sum of squares is minimized and the between-cluster sum of squares is maximized. Suppose we have a set of N data points {x1,x2,…,xN∣xn∈RD};  {ck∣k=1,2,…,K} is cluster centers. The measure of the distance is commonly chosen as Euclidian distance to minimize the following mean square error (MSE) cost function:
(4)CMSE=∑n=1N∑k=1Kon,kxn-ck2,
where on,k is a binary variable that indicates the cluster assignment of the point; if xn is assigned to cluster i then on,i=1 and on,j=0, for j≠i, j=1,2,…,K. In order to minimize (4) we need to find the set of cluster ownerships {on,k} and the set of cluster centers {ck}. Starting from a set of initial values of ck, this can be accomplished through an iterative procedure that interleaves between optimizing with respect to the cluster ownerships, keeping the cluster centers fixed, and optimizing with respect to the cluster centers while keeping the ownerships fixed [18]. The K-means algorithm can be summarized as follows.Step 1.
Generate the starting condition by defining the number of clusters and randomly select the initial cluster centers. For instance, select a set ofK candidate cluster centers {c1,c2,…,cK}.Step 2.
Assign eachxn to its nearest cluster center ck.Step 3.
Recalculate the cluster centers to the mean value of the points in each cluster.Step 4.
Repeat Steps2 and 3 until a distance convergence criterion is met.


## 2.1. NSCT

The NSCT is a shift-invariant version of the contourlet transform. The contourlet transform employs Laplacian pyramids (LPs) for multiscale decomposition and directional filter banks (DFBs) for directional decomposition [20]. Due to downsamplers and upsamplers present in both the LPs and DFBs, the contourlet transform is not shift-invariant. To achieve the shift-invariance, the NSCT is built upon nonsubsampled pyramids and nonsubsampled DFBs as shown in Figure 1.(a) NSFB structure that implements the NSCT. (b) Idealised frequency partitioning obtained with NSFB structure.
(a)
(b)(1) Nonsubsampled Pyramid (NSP). The nonsubsampled pyramid structure ensures the multiscale property of the NSCT. The building block of the nonsubsampled pyramid is two-channel nonsubsampled 2D filter banks. Figure 2(a) illustrates the nonsubsampled pyramid decomposition with J=3 stages [15]. The perfect reconstruction condition is given as
(1)H0zG0z+H1zG1z=1,
where H0(z) is low pass decomposition filter, H1(z) is high pass decomposition filter, G0(z) is low pass reconstruction filter, and G1(z) is high pass reconstruction filter.Nonsubsampled pyramid. (a) Three-stage pyramid decomposition. (b) Frequency divisions of a nonsubsampled pyramid.
(a)
(b)The ideal frequency support of the low pass filter at thejth stage is the region [-π/2j,π/2j]2. Accordingly, the support of the high pass filter is the complement of the low pass support region on the [-π/2j-1,π/2j-1]2∖[-π/2j,π/2j]2. The equivalent filters of a J level cascading NSPFB are given by
(2)Hnz=H1z2n-1∏j=0n-2H0z2j1≤n<2k∏j=0n-1H0z2jn=2k,
where H0(z) and H1(z) denote the low pass filter and the corresponding high pass filter at the first stage, respectively. The ideal frequency response of the nonsubsampled pyramid is given in Figure 2(b).(
2) Nonsubsampled Directional Filter Bank (NSDFB). The nonsubsampled DFB structure gives varying directions. The building block of the nonsubsampled DFB is also two-channel nonsubsampled 2D filter banks. The NSDFB is constructed by eliminating the downsamplers and upsamplers in the DFB [15]. To achieve multidirection decomposition, the NSDFB is iteratively computed. Figure 3 illustrates a four-channel directional decomposition with two-channel fan filter banks. The equivalent filter in each channel is given by
(3)Ukz=UizUjzQ.Four-channel NSDFB constructed with two-channel fan filter bank. (a) Filtering structure. (b) Corresponding frequency partitioning.
(a)
(b)(
3) Combining the Nonsubsampled Pyramid and Nonsubsampled Directional Filter Bank in the NSCT. The NSCT is constructed by combining the NSP and the NSDFB as shown in Figure 1(a). The NSP provides multiscale decomposition and captures the point discontinuities. The NSDFB provides directional decomposition and links point discontinuities into linear structures [21]. This scheme can be iterated repeatedly on the low pass subband outputs of nonsubsampled pyramids. The NSCT is very suitable for image edge detection because it has such important properties as multiresolution, localization, shift-invariance, and multidirection. In this paper, the “maxflat” filters and the “dmaxflat7” filters are, respectively, selected for NSPs and NSDFBs. The details of design for filters are referred to [15].


## 2.2.K-Means Clustering

TheK-means is a classical unsupervised clustering algorithm [22]. As an unsupervised technique, K-means clustering is effective in producing clusters for many practical applications such as image segmentation and image compression and determining the boundaries of the objects in an image.The aim ofK-means algorithm is to divide N points in D-dimensional space into K clusters so that the within-cluster sum of squares is minimized and the between-cluster sum of squares is maximized. Suppose we have a set of N data points {x1,x2,…,xN∣xn∈RD};  {ck∣k=1,2,…,K} is cluster centers. The measure of the distance is commonly chosen as Euclidian distance to minimize the following mean square error (MSE) cost function:
(4)CMSE=∑n=1N∑k=1Kon,kxn-ck2,
where on,k is a binary variable that indicates the cluster assignment of the point; if xn is assigned to cluster i then on,i=1 and on,j=0, for j≠i, j=1,2,…,K. In order to minimize (4) we need to find the set of cluster ownerships {on,k} and the set of cluster centers {ck}. Starting from a set of initial values of ck, this can be accomplished through an iterative procedure that interleaves between optimizing with respect to the cluster ownerships, keeping the cluster centers fixed, and optimizing with respect to the cluster centers while keeping the ownerships fixed [18]. The K-means algorithm can be summarized as follows.Step 1.
Generate the starting condition by defining the number of clusters and randomly select the initial cluster centers. For instance, select a set ofK candidate cluster centers {c1,c2,…,cK}.Step 2.
Assign eachxn to its nearest cluster center ck.Step 3.
Recalculate the cluster centers to the mean value of the points in each cluster.Step 4.
Repeat Steps2 and 3 until a distance convergence criterion is met.


## 3. Research Methodology

In this section, the proposed multiscale edge detection algorithm will be described in detail. This includes the feature classification usingK-means clustering algorithm, the NSCT modulus maxima processing, and the edge tracking for coarse-to-fine. A schematic representation of the proposed method is shown in Figure 4.Figure 4
Schematic representation of the proposed method.


### 3.1. Feature Classification UsingK-Means Clustering

In this section, we will first analyze and discuss the coefficients distribution of the NSCT in different subbands at a fixed level, for some typical locations. Then, we describe the feature classification usingK-means clustering. As described in Section 2.1, the NSCT is a fully shift-invariant, multiscale, and multidirection transform. It provides not only multiresolution analysis, but also geometric and directional representation [14]; thus the NSCT can efficiently capture the geometrical structures in natural images. Since edges are geometric structures, while noises are not, we can use the NSCT to distinguish them.The NSCT is shift-invariant, which leads to each pixel of the transform subbands corresponding to that of the original image in the same spatial location [15]. Therefore, we gather the geometric information from the NSCT coefficients for some typical spatial locations. Consider a simple image μ containing edges and smooth regions, like the one illustrated in Figure 5(a), and let NSCTμ[j,k,n] be the NSCT of μ, where j is the scale variable, k is the direction variable, and n is the location variable. We examine the behavior of the NSCTμ[j,k,n], at a fixed scale j0, for some typical locations n0. As Figure 5 shows, we can classify pixels into three categories by analyzing the plot of dn0(k)=NSCTμ(j0,k,n0), namely, strong edge points, weak edge (near edge) points, and smooth (noise) points. For instance, at the point n0=A on smooth edges and the junction point n0=C, dn0(k) have the large magnitude coefficients in all subbands (strong edge points); at the point n0=B, dn0(k) has large magnitude coefficients in some directional subbands and small magnitude coefficients in other directional subbands within the same scale (weak edge points); at the point n0=D,  dn0(k) has small magnitude coefficients in all subbands (smooth points). The same behavior holds for more general images, even in the presence of additive white Gaussian noise.The coefficients pattern for some typical locations across directional subbands at a fixed scale. (a) Representative points on the test image. (b)–(e) NSCT patterns(k), as a function of k for the points indicated in (a). (f) Representative points on the test image.
(a)
(b)
(c)
(d)
(e)
(f)This suggests devising a strategy for classifying strong edge points, weak edge points, and smooth points based on the plot patterns of Figure5(f). Let j0 be a fixed (fine) scale. In this paper, we choose to use two statistical feature variables (the energy and the maximum denoted by E and Max, resp.) of the coefficients for each pixel across directional subbands at this level. The main reason is the magnitude of the NSCT coefficients at the points along the boundaries (including strong edges and weak edges) which is significantly larger than at the other points. On the one hand, the E value has large difference to some boundary points and smooth (noise) points; they can be easily separated by looking at the energy function at n:
(5)En=∑kNSCTμj0,k,n2.On the other hand, theE value is similar to other weak edge points and smooth points but the Max with large difference; they can also be easily separated by looking at the maximum function at n:
(6)Max⁡n=argmax⁡kNSCTμj0,k,n.Therefore, we construct a two-dimensional feature vector (E, Max) based on the energy value E and the maximum Max⁡ for each pixel as shown in Figure 6. Next, using the (E, Max) as the feature function, we apply the K-means clustering algorithm, as described in Section 2.2, to cluster the image points into three sets X1, X2, and X3. Let S, W, and N denote the set of points Ximax⁡, Ximiddle, and Ximin⁡, respectively. This can be defined by
(7)imax⁡=argmax⁡i1Xi∑n∈XiEn,imin⁡=argmin⁡i1Xi∑n∈XiEn,
where S is the set of strong edge points and N is the set of smooth points. The remaining set W contains the weak edge points. In order to obtain a complete edge, we reserve the sets S and W, and the set N is set to 0. We denote the original edge map of an image μ as Mj,n={En,n∈S∪W} for j=1,2,…,J, where ∪ represents union. Figure 7 illustrates the application of the feature classification using K-means clustering algorithm to a test image. The results show that the feature classification using K-means clustering can effectively distinguish nonedge points from edge points.(a) The distribution of two-dimensional feature vectors corresponds to some typical locations. (b) The distribution of two-dimensional feature vectors corresponds to weak edge points and smooth points.
(a)
(b)Results of feature classification usingK-means clustering. (a) Original image. (b) Low pass subband image of the NSCT. (c)–(f) The results of feature classification of high pass subbands from coarser to finer.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.2. Identifying of the NSCT Modulus Maxima

In this section, we will describe the modulus maxima processing in NSCT domain. The theory of modulus maxima to edge detection lies in the fact that modulus maxima carry a significant degree of information about the position of edges. Mallat and Zhong [7] proved that the wavelet transform modulus maxima detect all the singularities of a function. They calculated the local maxima of discrete wavelet transform at each scale and formed a multiscale edge representation of an image. In what follows, we apply the modulus maxima in the NSCT domain.The NSCT provides a multidirectional structure, which can be used to efficiently obtain the edge direction for each pixel point. Figure8 illustrates an eight directions decomposition of the NSCT. As in the wavelet method, we will select the edge point candidates of an image μ by identifying the local maxima of the function Mj,n at fine scales. The NSCT modulus maxima algorithm consists of the following steps.Figure 8
The gradient map with directional decomposition number is 23.Step 5.
Input the raw edge subbandM(j,n) at each scale. The M(j,n) is obtained in Section 3.1.Step 6.
Estimate the edge directionϕ(j,n) by
(8)ϕj,n=argmax⁡kNSCTμj,k,n.In order to compare modulus values in any direction, we use the interpolation method (in our case, linear interpolation).Step 7.
At each edge point candidate, the magnitude ofM(j,n) is compared with the values of its neighbors along the gradient direction (this is obtained from (8)). If the magnitude is the largest, the point is kept. Otherwise, it is discarded and M(j,n) is set to 0.These produce a set of possible edge points denoted by{B(j,n),j=1,2,…,J}, where J is the decomposition level of the NSCT. The results of the NSCT modulus maxima processing are shown in Figure 9. We can find that the NSCT modulus maxima processing can effectively remove pseudo-edge responses.(a)–(d) Results of the NSCT modulus maxima processing from coarser to finer.
(a)
(b)
(c)
(d)


### 3.3. Edge Tracking

After modulus maxima processing, it is still a challenging problem to extract accurately edge image from the multiscale edge images because very few works propose techniques for combining multiscale information. To this end, we propose to combine the edges appearing at different scales by performing coarse-to-fine edge tracking. This algorithm is used to extract the complete and true edges and improve the robustness against noise. According to the multiscale property of the NSCT, for an imageμ, it provides scale space decomposition in the frequency domain. In general, the fine scales tend to produce results with a higher spatial accuracy, but also to be particularly sensitive to noise. Alternatively, the coarse scales are more robust against noise but tend to suffer from displacements of the edges from their true positions [23]. As shown in Figure 9, we observe that the nonedges disappear and the true edges progressively move away from their true positions and eventually merge or vanish when scale j increases. Besides, we also found that the displacements of the edges are along its gradient direction and the distance d=1 commonly between two consecutive scales. We chose Chebyshev distance as a measure of distance. We can see that the finest scale edge image has the problem of double edges and broken edges as shown in Figure 9(d). Therefore, we replace the finest scale edge image with canny edge image to obtain a complete and continuous edge image. The following sections describe whole edge tracking algorithm process in detail.Let the decomposition level of the NSCT beJ, and B={B1,B2,…,BJ} and Φ={Φ1,Φ2,…,ΦJ} represent the set of edge images and gradient orientation maps, respectively. In our implementation, the tracking starts from the coarsest scale, that is, the edge image BJ. In the first step, we remove isolated pixels by morphological processing of the Bj, where j is currently tracking scale variable. Then for each edge pixel (n1,n2)∈Bj, where Bj(n1,n2)≠0, we obtain the corresponding gradient direction Φj(n1,n2). Next, find the pixel (n1,n2)∈Bj-1, and search nonzero pixels Bj-1(x,y)≠0 in Bj-1, where (x,y) are these pixels satisfying that the maximum displacement d((x,y),(n1,n2))≤Tdist⁡ in direction Φj(n1,n2). If there are nonzero pixels in (x,y), then Bj(n1,n2) is set to 0 and these pixels satisfying Bj-1(x,y)≠0 are set to 255; namely, Bj(n1,n2)=0 and Bj-1(x,y)=255, respectively. We get a strong edge map S satisfying Bj-1(n1,n2)=255 and a weak edge map W satisfying Bj-1(n1,n2)=1. Then, we put S as a high threshold map and put W as a low threshold map, and perform hysteresis algorithm [5]. Therefore, we obtain an edge tracking map Bj-1, and the tracking starts from Bj-1 until the finest scale.We summarize our edge tracking algorithm as shown in Algorithm1.Algorithm 1
Input: the set of edge images B and the set of gradient orientation maps Φ.
Output:the final edge tracking map I.
Begin
Setj=J, where j is currently tracking scale variable and J is decomposition level of the NSCT.
While   j>1
(1) Remove isolated pixels for edge imageBj.
(2) Get gradient directionΦj(n1,n2) for each pixel (n1,n2)∈Bj, where Bj(n1,n2)≠0.
(3) Form the set of edge pixelsP={Bj-1(x,y)} for each pixel (n1,n2)∈Bj-1 corresponds to Bj(n1,n2)≠0,
whereP are these pixels satisfying: the maximum displacement d((x,y),(n1,n2))≤Tdist⁡ in direction
Φj(n1,n2) and Bj-1(x,y)≠0. If P∉⌀, then Bj(n1,n2)=0 and Bj-1(x,y)=255.
(4) Generate a strong edge mapS and a weak edge map W for Bj-1(n1,n2)=255 and Bj-1(n1,n2)=1,
respectively. Then perform edge connector forS and W, as is described in Section 3.3, the
result is assigned toBj-1.
(5) Perform the operationj=j-1 and go back to (1).
End
I
=
B.
EndThe edge tracking results at each scale are shown in Figures10(a)–10(d). We can observe that spurious responses are removed and produce edge results with a higher spatial accuracy. Meanwhile, we can also find that some significant edges are removed at finest scale but it is retained at subfinest scale. Therefore, we integrate next the edges at finest scale and the significant edges by removing small fragments at subfinest. To this end, the final edge map is obtained by the edge subsequent processing (e.g., edge connector and thinning processing). The final edge image using the proposed approach is shown in Figure 10(e).(a)–(d) The results of the edge tracking from coarser to finer. (e) The final edge map on original in Figure7(a).
(a)
(b)
(c)
(d)
(e)


## 3.1. Feature Classification UsingK-Means Clustering

In this section, we will first analyze and discuss the coefficients distribution of the NSCT in different subbands at a fixed level, for some typical locations. Then, we describe the feature classification usingK-means clustering. As described in Section 2.1, the NSCT is a fully shift-invariant, multiscale, and multidirection transform. It provides not only multiresolution analysis, but also geometric and directional representation [14]; thus the NSCT can efficiently capture the geometrical structures in natural images. Since edges are geometric structures, while noises are not, we can use the NSCT to distinguish them.The NSCT is shift-invariant, which leads to each pixel of the transform subbands corresponding to that of the original image in the same spatial location [15]. Therefore, we gather the geometric information from the NSCT coefficients for some typical spatial locations. Consider a simple image μ containing edges and smooth regions, like the one illustrated in Figure 5(a), and let NSCTμ[j,k,n] be the NSCT of μ, where j is the scale variable, k is the direction variable, and n is the location variable. We examine the behavior of the NSCTμ[j,k,n], at a fixed scale j0, for some typical locations n0. As Figure 5 shows, we can classify pixels into three categories by analyzing the plot of dn0(k)=NSCTμ(j0,k,n0), namely, strong edge points, weak edge (near edge) points, and smooth (noise) points. For instance, at the point n0=A on smooth edges and the junction point n0=C, dn0(k) have the large magnitude coefficients in all subbands (strong edge points); at the point n0=B, dn0(k) has large magnitude coefficients in some directional subbands and small magnitude coefficients in other directional subbands within the same scale (weak edge points); at the point n0=D,  dn0(k) has small magnitude coefficients in all subbands (smooth points). The same behavior holds for more general images, even in the presence of additive white Gaussian noise.The coefficients pattern for some typical locations across directional subbands at a fixed scale. (a) Representative points on the test image. (b)–(e) NSCT patterns(k), as a function of k for the points indicated in (a). (f) Representative points on the test image.
(a)
(b)
(c)
(d)
(e)
(f)This suggests devising a strategy for classifying strong edge points, weak edge points, and smooth points based on the plot patterns of Figure5(f). Let j0 be a fixed (fine) scale. In this paper, we choose to use two statistical feature variables (the energy and the maximum denoted by E and Max, resp.) of the coefficients for each pixel across directional subbands at this level. The main reason is the magnitude of the NSCT coefficients at the points along the boundaries (including strong edges and weak edges) which is significantly larger than at the other points. On the one hand, the E value has large difference to some boundary points and smooth (noise) points; they can be easily separated by looking at the energy function at n:
(5)En=∑kNSCTμj0,k,n2.On the other hand, theE value is similar to other weak edge points and smooth points but the Max with large difference; they can also be easily separated by looking at the maximum function at n:
(6)Max⁡n=argmax⁡kNSCTμj0,k,n.Therefore, we construct a two-dimensional feature vector (E, Max) based on the energy value E and the maximum Max⁡ for each pixel as shown in Figure 6. Next, using the (E, Max) as the feature function, we apply the K-means clustering algorithm, as described in Section 2.2, to cluster the image points into three sets X1, X2, and X3. Let S, W, and N denote the set of points Ximax⁡, Ximiddle, and Ximin⁡, respectively. This can be defined by
(7)imax⁡=argmax⁡i1Xi∑n∈XiEn,imin⁡=argmin⁡i1Xi∑n∈XiEn,
where S is the set of strong edge points and N is the set of smooth points. The remaining set W contains the weak edge points. In order to obtain a complete edge, we reserve the sets S and W, and the set N is set to 0. We denote the original edge map of an image μ as Mj,n={En,n∈S∪W} for j=1,2,…,J, where ∪ represents union. Figure 7 illustrates the application of the feature classification using K-means clustering algorithm to a test image. The results show that the feature classification using K-means clustering can effectively distinguish nonedge points from edge points.(a) The distribution of two-dimensional feature vectors corresponds to some typical locations. (b) The distribution of two-dimensional feature vectors corresponds to weak edge points and smooth points.
(a)
(b)Results of feature classification usingK-means clustering. (a) Original image. (b) Low pass subband image of the NSCT. (c)–(f) The results of feature classification of high pass subbands from coarser to finer.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.2. Identifying of the NSCT Modulus Maxima

In this section, we will describe the modulus maxima processing in NSCT domain. The theory of modulus maxima to edge detection lies in the fact that modulus maxima carry a significant degree of information about the position of edges. Mallat and Zhong [7] proved that the wavelet transform modulus maxima detect all the singularities of a function. They calculated the local maxima of discrete wavelet transform at each scale and formed a multiscale edge representation of an image. In what follows, we apply the modulus maxima in the NSCT domain.The NSCT provides a multidirectional structure, which can be used to efficiently obtain the edge direction for each pixel point. Figure8 illustrates an eight directions decomposition of the NSCT. As in the wavelet method, we will select the edge point candidates of an image μ by identifying the local maxima of the function Mj,n at fine scales. The NSCT modulus maxima algorithm consists of the following steps.Figure 8
The gradient map with directional decomposition number is 23.Step 5.
Input the raw edge subbandM(j,n) at each scale. The M(j,n) is obtained in Section 3.1.Step 6.
Estimate the edge directionϕ(j,n) by
(8)ϕj,n=argmax⁡kNSCTμj,k,n.In order to compare modulus values in any direction, we use the interpolation method (in our case, linear interpolation).Step 7.
At each edge point candidate, the magnitude ofM(j,n) is compared with the values of its neighbors along the gradient direction (this is obtained from (8)). If the magnitude is the largest, the point is kept. Otherwise, it is discarded and M(j,n) is set to 0.These produce a set of possible edge points denoted by{B(j,n),j=1,2,…,J}, where J is the decomposition level of the NSCT. The results of the NSCT modulus maxima processing are shown in Figure 9. We can find that the NSCT modulus maxima processing can effectively remove pseudo-edge responses.(a)–(d) Results of the NSCT modulus maxima processing from coarser to finer.
(a)
(b)
(c)
(d)


## 3.3. Edge Tracking

After modulus maxima processing, it is still a challenging problem to extract accurately edge image from the multiscale edge images because very few works propose techniques for combining multiscale information. To this end, we propose to combine the edges appearing at different scales by performing coarse-to-fine edge tracking. This algorithm is used to extract the complete and true edges and improve the robustness against noise. According to the multiscale property of the NSCT, for an imageμ, it provides scale space decomposition in the frequency domain. In general, the fine scales tend to produce results with a higher spatial accuracy, but also to be particularly sensitive to noise. Alternatively, the coarse scales are more robust against noise but tend to suffer from displacements of the edges from their true positions [23]. As shown in Figure 9, we observe that the nonedges disappear and the true edges progressively move away from their true positions and eventually merge or vanish when scale j increases. Besides, we also found that the displacements of the edges are along its gradient direction and the distance d=1 commonly between two consecutive scales. We chose Chebyshev distance as a measure of distance. We can see that the finest scale edge image has the problem of double edges and broken edges as shown in Figure 9(d). Therefore, we replace the finest scale edge image with canny edge image to obtain a complete and continuous edge image. The following sections describe whole edge tracking algorithm process in detail.Let the decomposition level of the NSCT beJ, and B={B1,B2,…,BJ} and Φ={Φ1,Φ2,…,ΦJ} represent the set of edge images and gradient orientation maps, respectively. In our implementation, the tracking starts from the coarsest scale, that is, the edge image BJ. In the first step, we remove isolated pixels by morphological processing of the Bj, where j is currently tracking scale variable. Then for each edge pixel (n1,n2)∈Bj, where Bj(n1,n2)≠0, we obtain the corresponding gradient direction Φj(n1,n2). Next, find the pixel (n1,n2)∈Bj-1, and search nonzero pixels Bj-1(x,y)≠0 in Bj-1, where (x,y) are these pixels satisfying that the maximum displacement d((x,y),(n1,n2))≤Tdist⁡ in direction Φj(n1,n2). If there are nonzero pixels in (x,y), then Bj(n1,n2) is set to 0 and these pixels satisfying Bj-1(x,y)≠0 are set to 255; namely, Bj(n1,n2)=0 and Bj-1(x,y)=255, respectively. We get a strong edge map S satisfying Bj-1(n1,n2)=255 and a weak edge map W satisfying Bj-1(n1,n2)=1. Then, we put S as a high threshold map and put W as a low threshold map, and perform hysteresis algorithm [5]. Therefore, we obtain an edge tracking map Bj-1, and the tracking starts from Bj-1 until the finest scale.We summarize our edge tracking algorithm as shown in Algorithm1.Algorithm 1
Input: the set of edge images B and the set of gradient orientation maps Φ.
Output:the final edge tracking map I.
Begin
Setj=J, where j is currently tracking scale variable and J is decomposition level of the NSCT.
While   j>1
(1) Remove isolated pixels for edge imageBj.
(2) Get gradient directionΦj(n1,n2) for each pixel (n1,n2)∈Bj, where Bj(n1,n2)≠0.
(3) Form the set of edge pixelsP={Bj-1(x,y)} for each pixel (n1,n2)∈Bj-1 corresponds to Bj(n1,n2)≠0,
whereP are these pixels satisfying: the maximum displacement d((x,y),(n1,n2))≤Tdist⁡ in direction
Φj(n1,n2) and Bj-1(x,y)≠0. If P∉⌀, then Bj(n1,n2)=0 and Bj-1(x,y)=255.
(4) Generate a strong edge mapS and a weak edge map W for Bj-1(n1,n2)=255 and Bj-1(n1,n2)=1,
respectively. Then perform edge connector forS and W, as is described in Section 3.3, the
result is assigned toBj-1.
(5) Perform the operationj=j-1 and go back to (1).
End
I
=
B.
EndThe edge tracking results at each scale are shown in Figures10(a)–10(d). We can observe that spurious responses are removed and produce edge results with a higher spatial accuracy. Meanwhile, we can also find that some significant edges are removed at finest scale but it is retained at subfinest scale. Therefore, we integrate next the edges at finest scale and the significant edges by removing small fragments at subfinest. To this end, the final edge map is obtained by the edge subsequent processing (e.g., edge connector and thinning processing). The final edge image using the proposed approach is shown in Figure 10(e).(a)–(d) The results of the edge tracking from coarser to finer. (e) The final edge map on original in Figure7(a).
(a)
(b)
(c)
(d)
(e)


## 4. Experiments and Results

In this section, we intend to compare the performance of the proposed method with some conventional edge detection algorithms, including the Sobel edge detector [2], the LoG edge detector [24], the Canny edge detector [5], and the M-Sobel edge detector [23]. In the experiments we have used the test set of the Berkeley segmentation dataset (BSDS) because it has the advantage of providing multiple human labeled segmentations per image as shown in Figure 11. Thus we use them as ground truth for the quantification of the quality of the edge detection results and these images have a resolution of 481×321 pixels.Berkeley segmentation data set [16]. Top to bottom: image and ground-truth segment boundaries hand drawn by three different human subjects. The BSDS500 consists of 300 training and 200 test images, each with multiple ground-truth segmentations.
(a)
(b)
(c)
(d)We first use different scale decomposition values of the NSCT to generate the final edge images for testing the proposed approach, which are[2,3], [2,2,3], [2,2,2,2], [2,2,3,3], and [2,3,3,3]. The extracted edges using the proposed approach are shown in Figure 12. We can observe that the edge detector can reduce the noised influence with the increase of the number of scale decomposition of the NSCT, especially for noisy image as shown in Figures 12(f), 12(g), and 12(j). As previously described in Section 3.3, with the NSCT scale changing from finer to coarser, the nonedges disappear and the true edges progressively move away from their true positions and eventually merge or vanish. So the slight improvement comes at the cost of a higher computational complexity with the increase of the number of scale decomposition. In addition, we also found that, when the number of scale decomposition is the same, the extracted edges are more complete in the case of 23 direction decomposition. The simulation results are shown in Figures 12(c)–12(e) and Figures 12(h)–12(j). Therefore, we use four scales of decomposition of the NSCT for the following experiments, and the direction numbers of decomposition are 2,3,3,3 from coarser to finer, respectively.(a)–(e) The extracted edges based on the scale decomposition of[2,3], [2,2,3], [2,2,2,2], [2,2,3,3], and [2,3,3,3] for original images. (f)–(j) The extracted edges based on the scale decomposition of [2,3], [2,2,3], [2,2,2,2], [2,2,3,3], and [2,3,3,3] for noisy images.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Then we have compared the performances of the proposed method with the Canny method, the LoG method, the Sobel method, and the M-Sobel method. Figure13(a) is the original images, and the edge detection results using different methods are shown in Figures 13(b)–13(f). We can see that the Canny method and LoG method detect many false edges when they try to obtain more real edges, which can be seen in Figures 13(b) and 13(c). The Sobel method has better edge detection results and less false edges. However, the detected edges are usually incomplete which can be clearly observed in Figure 13(d). Figure 13(f) shows the detection results using the proposed method. It can be seen that our method can remove the false edges and keep the real edges as much as possible. In addition, we can also find that the proposed method compared with the M-Sobel method (Figure 13(e)) has similar or better edge detection results for original images. Meanwhile, the edge detection results of our method are more complete and continuous than the edge images detected using the other three methods. These are mainly because the NSCT has the ability to fully capture directional and other geometrical features for images with edges.Edge detection results using different methods for original images. (a) Original images. (b) Canny edge detector. (c) LoG edge detector. (d) Sobel edge detector. (e) M-Sobel edge detector. (f) Our method.
(a)
(b)
(c)
(d)
(e)
(f)Next, experiments are tested to validate that the proposed approach is efficient to the images with Gaussian white noise, where the noise variance isδn2=0.001 in our experiments. Figure 14(a) is the original images and the extracted results using different methods for noisy images are shown in Figures 14(b)–14(f). We can find that the proposed approach has better robustness against noise than the other four methods. The main reasons are as follows: in general, fine scales of the NSCT provide spatially accurate results but are particularly sensitive to noise; meanwhile, coarse scales of the NSCT have better robustness against noise but tend to suffer from displacements of the edges from their actual position. Thus we achieved better edge detection results by performing coarse-to-fine edge tracking presented in this paper.Edge detection results using different methods for noisy images. (a) Noisy images. (b) Canny edge detector. (c) LoG edge detector. (d) Sobel edge detector. (e) M-Sobel edge detector. (f) Our method.
(a)
(b)
(c)
(d)
(e)
(f)Finally, in order to objectively evaluate the edge detection performances of the proposed approach and other methods, we use the classification inspired methodology by Martin et al. [25]. In essence, edge detection can be seen as a binary classification problem. Therefore, it can be evaluated in terms of success and fallout, comparing the candidate edge image (generated by an edge detection method) with the ground truth (generated by a human). According to the ground truth, the pixels in the candidate edge image can be classified into four different categories: true positive (TP), false positive (FP), true negative (TN), and false negative (FN) [6]. In this way, we also build a confusion matrix as shown in Table 1 [23]. From the confusion matrix we extract the precision and recall evaluations, defined as
(9)PREC=TPTP+FP,REC=TPTP+FN.Table 1
Confusion matrix for the edge detection evaluation.
Reality
Classification
Edge
Nonedge
Edge
TP
FN
Nonedge
FP
TNThe precision and recall measures hold good stability properties and illustrate specific aspects of the problem. Therefore, we use the F-measure by Martin et al. [25] to evaluation of the overall quality of an edge image, defined as
(10)Fα=PREC·RECαPREC+(1-α)REC,
where α is a parameter weighing the relative impact of the precision and recall evaluations. In addition, Pratt’s figure of merit (PFoM) [26] is used for quantitative evaluation of spatial accuracy, defined as
(11)Pl=1max⁡nc,ngt∑k=1nc11+βd2k,
where nc is the number of the candidate edges and ngt is the number of the ground-truth edges. d(k) denotes the Euclidian distance from the kth candidate edge to the corresponding ground-truth edge. β is a constant, generally fixed at 1/9 for keeping the values in [0,1]. A better approach possesses larger Fα and Pl values.In this work, we use a one-to-one pixel matching algorithm to map the edge pixels in the candidate edge image and the ground truth. This matching allows for a certain spatial tolerance (in our case, as much as 0.5% of the diagonal of the image), so that an edge pixel can be slightly moved from its true position, yet being considered as correctly classified. For a given image, the result of each edge detection method is compared with several ground-truth images, each one producing a tuple of evaluations(F0.5,Pl). Then, the comparison producing the best quality assessment, in terms of F0.5, is the one we select. Tables 2 and 3 show the results obtained of each detector in Figures 13(a) and 14(a), respectively. The results show that the proposed approach is quantitatively better than all other methods. More specifically, the proposed method can achieve better results with the typical methods for original images. In particular, it also works well for noisy images and the advantages are obvious.Table 2
TheF0.5 and Pl values, where the 1, 2, 3, and 4 denote top, middle top, middle bottom, and bottom original images in Figure 12(a), respectively.
Original images
Comparison to other methods
Canny
LoG
Sobel
M-Sobel
Our method
1
F
0.5
0.4626
0.5626
0.8452
0.8588
0.9444
P
l
0.4064
0.4691
0.8193
0.8377
0.8351
2
F
0.5
0.2672
0.3240
0.5829
0.6384
0.6271
P
l
0.2046
0.2726
0.5376
0.6243
0.5685
3
F
0.5
0.4474
0.5038
0.6418
0.6560
0.6832
P
l
0.3716
0.5968
0.7227
0.7690
0.7310
4
F
0.5
0.7468
0.7149
0.8484
0.8500
0.8740
P
l
0.6459
0.6839
0.8527
0.8807
0.7792Table 3
TheF0.5 and Pl values, where the 1, 2, 3, and 4 denote top, middle top, middle bottom, and bottom noise images in Figure 13(a), respectively.
Noise images
Comparison to other methods
Canny
LoG
Sobel
M-Sobel
Our method
1
F
0.5
0.1825
0.2004
0.8210
0.8406
0.9404
P
l
0.1561
0.1806
0.7916
0.8343
0.7918
2
F
0.5
0.2002
0.3086
0.5582
0.5848
0.5944
P
l
0.1540
0.2617
0.5611
0.5458
0.5099
3
F
0.5
0.4488
0.5028
0.5368
0.5491
0.6876
P
l
0.3903
0.5920
0.4068
0.5037
0.7349
4
F
0.5
0.3402
0.6992
0.8078
0.8083
0.8686
P
l
0.2535
0.6756
0.8151
0.8405
0.7841


## 5. Conclusions

In this paper, we developed a multiscale edge detection approach based on NSCT and edge tracking. This method is able to distinguish the real edges from the noise and overcome the shortcoming of the limited directions of wavelet transform. In order to produce the final edge image, we have proposed a coarse-to-fine edge tracking algorithm, to improve the edge localization accuracy and yield better edge detection results. Experimental results demonstrate that the proposed approach achieves better edge detection performance compared with the typical methods. In particular, it provides an accurate method for extracting the information about edges and their locations even in presence of noise. Nevertheless, the proposed method has a higher computational complexity due to the NSCT.



---

*Source: 504725-2015-02-22.xml*