# Image Edge Detection Based on Gaussian Mixture Model in Nonsubsampled Contourlet Domain

Authors: Li Yang; Chang Xia; Chang Juan
Journal: Journal of Electrical and Computer Engineering
 (2016)

---

## Abstract

In order to get accurate location and continuous edges, Gaussian mixture model and local direction modulus nonmaxima suppression are used in high frequency subbands of nonsubsampled Contourlet transform. The distribution of NSCT high frequency subbands coefficients has the “high spikes, long tail” non-Gaussian statistical characteristic. Gaussian mixture model (GMM) is used to distinguish the linear singular signal and the nonlinear singular signal on the high frequency subbands. Local direction modulus nonmaxima suppression is used to refine the linear singular signal. An appropriate threshold is used to distinguish edge pixels and nonedge pixels to get binary image. The experimental results demonstrate that the proposed method can capture more continuous edges in multiple directions and has accurate edge location. And the edges are with great convenience for the image recognition.


---

## Body

## 1. Introduction

Image edge detection plays a very important role in image processing [1]. The results of image edge detection will directly affect the performance of image understanding, analysis, and recognition. The classical edge detectors based on gradient are concise and easy to understand, such as Sobel and Prewitt operators. Because of doing the average operation in these algorithms, they will lose some detailed information and make the edges blur. So the accuracy of location is not high. Canny edge detection usually gives good performance. Canny algorithm uses two different thresholds to detect strong edges and weak edges, respectively. Weak edges will be detected when weak edges are connected to strong edges. But weak edges will be ignored when weak edges are not connected to strong edges. Due to using Gaussian filter in the Canny edge detection method, some weak edges will be smoothed. And the detected edges obtained are more miscellaneous, so it is difficult to identify the main contours.As the multiscale analysis is introduced, the image processing methods based on wavelet have obtained good results. Wavelet transform is famous for its time frequency localization, multiscale, and multiresolution. It is a very effective tool for image edge detection in image processing [2], such as the wavelet modulus maxima edge detection method. But the basis functions of wavelet transform are isotropic, so it cannot capture more directional information in images [3]. In the recent years, multiscale geometric analysis has grown constantly, and new methods based on the new multiscale transforms are proposed, such as the curvelet and Contourlet transform [4]. The Contourlet transform was proposed to address the lack of geometrical structure in the separable two-dimensional wavelet transform. The research results show that the performance of the image edge detection based on the new multiscale transforms is better than those based on the wavelet transform. The Contourlet transform is a real 2D image representation using a directional filter and a Laplacian pyramid, which can effectively capture contours in an image and can achieve better expression of image than the wavelet transform. It has the anisotropic characteristics. Because of the process of downsampling and upsampling, the Contourlet transform is shift-variance and always Gibbs phenomena around singularities. In 2006, the nonsubsampled Contourlet transform (NSCT) [5] was proposed by Cunha et al. The NSCT has better performance of the NSCT in image denoising and enhancement applications. NSCT is composed of the nonsubsampled pyramid filter and the nonsubsampled direction filter, so it is shift-invariance [6].In this paper, a method of edge detection based on nonsubsampled Contourlet transform and Gaussian mixture model is proposed. The probability and statistics analysis is applied in the NSCT coefficients. The Gaussian mixture model is used to distinguish the linear singular signal and nonlinear singular signal on the high frequency subbands automatically. Local direction modulus nonmaxima suppression is used to refine the linear singular signal. An appropriate threshold is used to distinguish edge pixels and nonedge pixels to get binary image. We compare the proposed method with the classical Canny operator, Sobel operator, the wavelet method [2], and the method based on cellular neural network (CNN) [7]. The proposed method can capture more continuous edges in multiple directions and has accurate edge location.


## 2. Modeling NSCT High Frequency Coefficients

### 2.1. Nonsubsampled Contourlet Transform

NSCT consisted of the nonsubsampled pyramid filter banks (NSPFB) and the nonsubsampled directional filter banks (NSDFB) [6]. The two independent parts are shift-invariant. NSPFB is used for multiscale decomposition. NSDFB uses singular points of the same direction to synthesize the NSCT coefficients. NSDFB is shown in Figure 1(a). Frequency resolution is shown in Figure 1(b).Figure 1
The structure of nonsubsampled Contourlet transform.
(a)
Nonsubsampled directional filter bank (b)
Frequency resolutionThe distribution of NSCT high frequency subbands coefficients has the “high spikes, long tail” non-Gaussian statistical characteristic, as shown in Figure2(a). High frequency subband coefficients consist of a few “big” coefficients and most “small” coefficients. The “small” coefficients contain less information and the “big” coefficients contain main information. Other high frequency directional subbands also have the similar distribution. Based on these characteristics to distinguish the linear singular signal and the nonlinear singular signal, Gaussian mixture model is used to model subbands, as shown in Figure 2(b); we use one Gaussian density function with big variance and one Gaussian density function with small variance to approximate the distribution. The “big” coefficients represent the linear singular signal. The linear singular signal is corresponding to the image edges.Figure 2
Statistics of the image NSCT coefficients.
(a)
“High spikes, long tail” non-Gaussian statistical characteristic (b)
The two classes of the coefficients


### 2.2. Parameter Estimation of Gaussian Mixture Model

We model each high frequency directional subband through the Gaussian mixture model.xv,w is the coefficient which is in the v scale and the w high frequency directional subband; N is the number of all coefficients in xv,w. Let xm,nv,w be in xv,w matrix, where m is the column and n is the row. Assume that each sample in xm,nv,w is independent of each other and obeys the K Gaussian distribution [8], K=2. xm,nv,w will be divided into two classes: one class is “big” state; the other is “small” state.Density function ofxm,nv,w is fv,w(x)=∑k=0Kπkv,wNv,wxm,nv,w∣μkv,w,Σkv,w. There, ∑k=0Kπkv,w=1, where k=0 or 1: 0 means “big” state and 1 means “small” state; μkv,w is the k class means; Σkv,w is the k class covariance, where v is the decomposition scale and w is the filter direction.Expectation maximization algorithm can be used to calculate the optimal parameters of GMM by E-Step and M-Step.E-Step: calculating the probability ofxm,nv,w generated by each component(1)γv,wi,k=πkv,wNxm,nv,w∣μkv,w,Σkv,w∑j=0Kπjv,wNxm,nv,w∣μjv,w,Σjv,w.In formula (1), γv,w(i,k) is the probability of xm,nv,w and πkv,w is the weighting factor.M-Step: calculatingNkv,w, ukv,w, and πkv,w, iteratively(2)Nkv,w=∑i=1Nγv,wi,k,(3)ukv,w=1Nkv,w∑i=1Nγv,wi,kxm,nv,w,(4)πkv,w=Nkv,wN.


### 2.3. Algorithm1: Data Screening Algorithm by Gaussian Mixture Model

See Algorithm1.Algorithm 1:Data screening algorithm by Gaussian mixture model.
Input: high-frequency directional sub-bandsxm,nv,w (1) Employ (1)~(3) to calculate the probability of xm,nv,w generated by each component (2) Employ (4)~(7) to calculate the ukv,w, πkv,w, Nkv,w (3) Repeated iteration, until the values of the likelihood function achieves convergence. (4) Calculate the probability of xi generated by each component according to xm,nv,w, ukv,w, π k v , w, Nkv,w. Let pro0 be the probability of xm,nv,w generated by the “big” variance of Gaussian model; let pro1 be the probability ofxm,nv,w generated by the “small” variance of Gaussian model. (5) If pro0 < pro1, set 0 to xm,nv,w Output:xm,nv,w (new high-frequency directional sub-bands xm,nv,w)


### 2.4. Local Direction Modulus Nonmaxima Suppression

After executing Algorithm1, we use local direction modulus nonmaxima suppression to refine the linear signal. High frequency directional subbands contain the directional information, so there is no need to calculate the direction of the modulus; just compare the modulus of pixel and the modulus of the two gradient directional pixels of this pixel [9]. Determine whether it is a local maximum modulus and retain if it is.For example, there are eight directions in some new high frequency directional subbands. So there are eight sectors, which areπ/16, 3π/16, 5π/16, 7π/16, 9π/16, 11π/16, 13π/16, and 15π/16. Mod (xm,nv,w) is the modulus of xm,nv,w. argv,w is the equivalent gradient of xv,w. Figure 3 is the directional subband equivalent gradient.Figure 3
Directional subband equivalent gradient.


### 2.5. Algorithm2: Local Mold Maxima Suppression Algorithm

See Algorithm2.Algorithm 2:Local mold maxima suppression algorithm.
Input: new high-frequency directional sub-bandsxm,nv,w If3π/16<argv,w<5π/16∥11π/16<argv,w<13π/16 Ifmod⁡(xm,nv,w)<mod⁡(xm,n-1v,w)∥mod⁡(xm,nv,w)<mod(xm,n+1v,w) xm,nv,w=0; End End ElseIfπ/16<argv,w<3π/16∥9π/16<argv,w<11π/16 Ifmod⁡(xm,nv,w)<mod⁡(xm-1,n+1v,w)∥mod⁡(xm,nv,w)<mod(xm+1,n-1v,w) xm,nv,w=0; End End ElseIf5π/16<argv,w<7π/16∥13π/16<argv,w<15π/16 Ifmod⁡(xm,nv,w)<mod⁡(xm-1,nv,w)∥mod⁡(xm,nv,w)<mod⁡(xm+1,nv,w) xm,nv,w=0; End End Else Ifmod⁡(xm,nv,w)<mod⁡(xm-1,nv,w)∥mod⁡(xm,nv,w)<mod⁡(xm+1,nv,w) xm,nv,w=0; End End Output:xm,nv,w (new high-frequency directional sub-bands xm,nv,w)


### 2.6. Our Proposed Edge Detection Algorithm

The proposed edge detection algorithm is summarized as follows:Input: original image.Step 1.
Decompose the original image into high frequency directional subbandsxm,nv,w and low frequency subband by NSCT.Step 2.
Let low frequency subband coefficients be 0.Step 3.
Employ Algorithm1 to process xm,nv,w, we get new xm,nv,w, and new xm,nv,w are the linear singular signal here.Step 4.
Employ Algorithm2 to process xm,nv,w, which is got by Step 3, we get new xm,nv,w, and new xm,nv,w are the linear singular signal refined here.Step 5.
Reconstruct image using high frequency coefficientsxm,nv,w processed by Step 3 and low frequency coefficients processed by Step 1.Step 6.
After NSCT inverse transform, there will be a data matrix. There will be two types of data in that matrix, the modulus of one type is big, and the modulus of the other type is small. An appropriate threshold can be easily set to separate them. Set threshold to distinguish edge pixels and nonedge pixels to get binary image.Output: edge detection result.


## 2.1. Nonsubsampled Contourlet Transform

NSCT consisted of the nonsubsampled pyramid filter banks (NSPFB) and the nonsubsampled directional filter banks (NSDFB) [6]. The two independent parts are shift-invariant. NSPFB is used for multiscale decomposition. NSDFB uses singular points of the same direction to synthesize the NSCT coefficients. NSDFB is shown in Figure 1(a). Frequency resolution is shown in Figure 1(b).Figure 1
The structure of nonsubsampled Contourlet transform.
(a)
Nonsubsampled directional filter bank (b)
Frequency resolutionThe distribution of NSCT high frequency subbands coefficients has the “high spikes, long tail” non-Gaussian statistical characteristic, as shown in Figure2(a). High frequency subband coefficients consist of a few “big” coefficients and most “small” coefficients. The “small” coefficients contain less information and the “big” coefficients contain main information. Other high frequency directional subbands also have the similar distribution. Based on these characteristics to distinguish the linear singular signal and the nonlinear singular signal, Gaussian mixture model is used to model subbands, as shown in Figure 2(b); we use one Gaussian density function with big variance and one Gaussian density function with small variance to approximate the distribution. The “big” coefficients represent the linear singular signal. The linear singular signal is corresponding to the image edges.Figure 2
Statistics of the image NSCT coefficients.
(a)
“High spikes, long tail” non-Gaussian statistical characteristic (b)
The two classes of the coefficients


## 2.2. Parameter Estimation of Gaussian Mixture Model

We model each high frequency directional subband through the Gaussian mixture model.xv,w is the coefficient which is in the v scale and the w high frequency directional subband; N is the number of all coefficients in xv,w. Let xm,nv,w be in xv,w matrix, where m is the column and n is the row. Assume that each sample in xm,nv,w is independent of each other and obeys the K Gaussian distribution [8], K=2. xm,nv,w will be divided into two classes: one class is “big” state; the other is “small” state.Density function ofxm,nv,w is fv,w(x)=∑k=0Kπkv,wNv,wxm,nv,w∣μkv,w,Σkv,w. There, ∑k=0Kπkv,w=1, where k=0 or 1: 0 means “big” state and 1 means “small” state; μkv,w is the k class means; Σkv,w is the k class covariance, where v is the decomposition scale and w is the filter direction.Expectation maximization algorithm can be used to calculate the optimal parameters of GMM by E-Step and M-Step.E-Step: calculating the probability ofxm,nv,w generated by each component(1)γv,wi,k=πkv,wNxm,nv,w∣μkv,w,Σkv,w∑j=0Kπjv,wNxm,nv,w∣μjv,w,Σjv,w.In formula (1), γv,w(i,k) is the probability of xm,nv,w and πkv,w is the weighting factor.M-Step: calculatingNkv,w, ukv,w, and πkv,w, iteratively(2)Nkv,w=∑i=1Nγv,wi,k,(3)ukv,w=1Nkv,w∑i=1Nγv,wi,kxm,nv,w,(4)πkv,w=Nkv,wN.


## 2.3. Algorithm1: Data Screening Algorithm by Gaussian Mixture Model

See Algorithm1.Algorithm 1:Data screening algorithm by Gaussian mixture model.
Input: high-frequency directional sub-bandsxm,nv,w (1) Employ (1)~(3) to calculate the probability of xm,nv,w generated by each component (2) Employ (4)~(7) to calculate the ukv,w, πkv,w, Nkv,w (3) Repeated iteration, until the values of the likelihood function achieves convergence. (4) Calculate the probability of xi generated by each component according to xm,nv,w, ukv,w, π k v , w, Nkv,w. Let pro0 be the probability of xm,nv,w generated by the “big” variance of Gaussian model; let pro1 be the probability ofxm,nv,w generated by the “small” variance of Gaussian model. (5) If pro0 < pro1, set 0 to xm,nv,w Output:xm,nv,w (new high-frequency directional sub-bands xm,nv,w)


## 2.4. Local Direction Modulus Nonmaxima Suppression

After executing Algorithm1, we use local direction modulus nonmaxima suppression to refine the linear signal. High frequency directional subbands contain the directional information, so there is no need to calculate the direction of the modulus; just compare the modulus of pixel and the modulus of the two gradient directional pixels of this pixel [9]. Determine whether it is a local maximum modulus and retain if it is.For example, there are eight directions in some new high frequency directional subbands. So there are eight sectors, which areπ/16, 3π/16, 5π/16, 7π/16, 9π/16, 11π/16, 13π/16, and 15π/16. Mod (xm,nv,w) is the modulus of xm,nv,w. argv,w is the equivalent gradient of xv,w. Figure 3 is the directional subband equivalent gradient.Figure 3
Directional subband equivalent gradient.


## 2.5. Algorithm2: Local Mold Maxima Suppression Algorithm

See Algorithm2.Algorithm 2:Local mold maxima suppression algorithm.
Input: new high-frequency directional sub-bandsxm,nv,w If3π/16<argv,w<5π/16∥11π/16<argv,w<13π/16 Ifmod⁡(xm,nv,w)<mod⁡(xm,n-1v,w)∥mod⁡(xm,nv,w)<mod(xm,n+1v,w) xm,nv,w=0; End End ElseIfπ/16<argv,w<3π/16∥9π/16<argv,w<11π/16 Ifmod⁡(xm,nv,w)<mod⁡(xm-1,n+1v,w)∥mod⁡(xm,nv,w)<mod(xm+1,n-1v,w) xm,nv,w=0; End End ElseIf5π/16<argv,w<7π/16∥13π/16<argv,w<15π/16 Ifmod⁡(xm,nv,w)<mod⁡(xm-1,nv,w)∥mod⁡(xm,nv,w)<mod⁡(xm+1,nv,w) xm,nv,w=0; End End Else Ifmod⁡(xm,nv,w)<mod⁡(xm-1,nv,w)∥mod⁡(xm,nv,w)<mod⁡(xm+1,nv,w) xm,nv,w=0; End End Output:xm,nv,w (new high-frequency directional sub-bands xm,nv,w)


## 2.6. Our Proposed Edge Detection Algorithm

The proposed edge detection algorithm is summarized as follows:Input: original image.Step 1.
Decompose the original image into high frequency directional subbandsxm,nv,w and low frequency subband by NSCT.Step 2.
Let low frequency subband coefficients be 0.Step 3.
Employ Algorithm1 to process xm,nv,w, we get new xm,nv,w, and new xm,nv,w are the linear singular signal here.Step 4.
Employ Algorithm2 to process xm,nv,w, which is got by Step 3, we get new xm,nv,w, and new xm,nv,w are the linear singular signal refined here.Step 5.
Reconstruct image using high frequency coefficientsxm,nv,w processed by Step 3 and low frequency coefficients processed by Step 1.Step 6.
After NSCT inverse transform, there will be a data matrix. There will be two types of data in that matrix, the modulus of one type is big, and the modulus of the other type is small. An appropriate threshold can be easily set to separate them. Set threshold to distinguish edge pixels and nonedge pixels to get binary image.Output: edge detection result.


## 3. Experimental Results

In order to quantify the performance of the proposed method, FOM (Pratt’s Figure of Merit) [10] is used to compare the performance of edge detection methods, because the edges of the standard image must be obtained while using FOM. Because of this, we can just use the synthetic images. The FOM is computed as follows:(5)FOM=1max⁡NI,NA∑i=1NA11+αdi,where NI and NA are the number of ideal and detected edge pixels, respectively, the parameter α is a constant, and di is the vertical distance between an actual edge pixel and the nearest ideal edge pixel. The FOM measures the accuracy of localization. The more the FOM is, the more the quality of the edge detection is. In order to objectively discuss the accuracy of the localization, the false edge pixels are not covered while computing the FOM. A synthetic image is used to calculate FOM.The proposed method, CNN, Canny, Sobel operators, NSCT, and the wavelet method are used to detect the edge of the synthetic image. The NSCT method is used to detect edges without using GMM. Figure4 shows the results detected by the different methods. Table 1 shows the FOM calculations by different methods. The results demonstrate that the accuracy of the localization of the proposed method is more accurate than Canny operator, Sobel operator, wavelet method, NSCT method, and the method based on CNN. The accuracy of the localization of the proposed method is the best. Because of not using GMM to distinguish the linear singular signal and the nonlinear singular signal on the high frequency subbands, the edges detected by the NSCT method are interrupted by the nonlinear singular signal. So the accuracy of the localization of the NSCT method is very low. The Canny, Sobel operators, and CNN method are sensitive to noise, so their accuracy of the localization is lower than the proposed method.Table 1
The FOM calculations with different methods.
Figure4 The proposed method The method based on CNN Canny operator Sobel operator NSCT Wavelet method FOM 0.9631 0.8934 0.8256 0.7937 0.6432 0.7231Figure 4
The results of the edge detection: (a) the synthetic image, (b) the proposed method, (c) the method based on CNN, (d) Canny operator, (e) Sobel operator, (f) wavelet method, and (g) NSCT method.
(a) (b) (c) (d) (e) (f) (g)Three standard images, which are Barbara (size256∗256), Lenna (size 256∗256), and boat (size 512∗512), are used to test in the experiments. We select “9-7” pyramid decomposition and “pkva” directional filter bank for NSCT decomposition. Decomposition level is 3 and each level has eight high frequency directions. Experimental results are shown in Figures 5, 6, and 7.Figure 5
Edge detection results: (a) original image; (b) the result by Canny; (c) the result by our method; (d) the result by Sobel; (e) the result by wavelet method; (f) the method based on CNN; (g) NSCT method.
(a) (b) (c) (d) (e) (f) (g)Figure 6
Edge detection results: (a) original image; (b) the result by Canny; (c) the result by our method; (d) the result by Sobel; (e) the result by wavelet method; (f) the method based on CNN; (g) NSCT method.
(a) (b) (c) (d) (e) (f) (g)Figure 7
Edge detection results: (a) original image; (b) the result by Canny; (c) the result by our method; (d) the result by Sobel; (e) the result by wavelet method; (f) the method based on CNN; (g) NSCT method.
(a) (b) (c) (d) (e) (f) (g)The linear singular edges of the table cloth detected in Figure5(b) are clearer than those in Figures 5(c), 5(d), 5(e), and 5(f). The texture edges of the tassels in Figure 6(b) are more continuous than those in Figures 6(c), 6(d), 6(e), and 6(f). Contour edges detected of the hull in Figure 7(b) are more complete than those in Figures 7(c), 7(d), and 7(f). The edges are too miscellaneous to distinguish the main contours in Figures 5(c), 6(c), 7(c), 5(d), 6(d), and 7(d). Because the wavelet just has two directions, it cannot capture more directional details. Because of not using GMM to distinguish the linear singular signal and the nonlinear singular signal on the high frequency subbands, the edges detected by the NSCT method are interrupted by the nonlinear singular signal. From Figures 5(g), 6(g), and 7(g), we also cannot distinguish the main contours. Because the NSCT has good directionality, anisotropy, and decorrelation, the proposed method can capture more directional details. And the proposed method can capture the main contours because of using the GMM. From Figures 5(b), 6(b), and 7(b), we can clearly see the main contours in the images.In order to evaluate the algorithm, the quantitative analysis is adopted [11]. In Table 2, CEN denotes the number of continuous edge pixels in edge image. TEN represents the number of all pixels in edge image. R denotes the ratio of TEN to CEN. The more R is, the more image edge is continuous. R is computed as follows:(6)R=TENCEN.Table 2
The comparison of the different detection methods.
Image Detection method CEN TEN R Barbara Canny 9896 8713 0.8805 Our method 16113 14984 0.9299 Sobel 8646 7163 0.82848 Wavelet 7556 6388 0.8454 CNN 12703 10793 0.8496 NSCT 21822 19790 0.9069  Lenna Canny 7935 7184 0.9041 Our method 8462 7924 0.9364 Sobel 6731 5831 0.86629 Wavelet 5470 4822 0.8815 CNN 4544 3627 0.7982 NSCT 23823 22043 0.9253  Boat Canny 32011 29030 0.9069 Our method 56781 54791 0.965 Sobel 34536 30172 0.87364 Wavelet 25920 23479 0.9058 CNN 18805 13647 0.7257 NSCT 80251 73809 0.9197Table2 shows that the edges obtained by the proposed method are more continuous than other methods. Because the wavelet just has two directions, it cannot capture more directional details. Though the NSCT method can capture more edges than the Canny and Sobel operators, the pixels detected by the NSCT method contain so much the nonlinear singular signal. The proposed method can capture more details than other methods.MSE (Mean Squared Error) is a standard to measure the accuracy of localization. Generally, when MSE is greater, the accuracy of localization is less. When MSE is less, the accuracy of localization is higher. The MSE is computed as follows:(7)MSE=1mn∑i=0m-1∑j=0n-1Ii,j-Ki,j2.Table3 shows the MSE results of different methods. We can see that the proposed method has a lower MSE than other methods. It means that the proposed method has high accuracy of localization. Because of not using GMM to distinguish the linear singular signal and the nonlinear singular signal on the high frequency subbands, the edges detected by the NSCT method are interrupted by the nonlinear singular signal. So the accuracy of localization is low in the NSCT method.Table 3
MSE of the different detection methods.
Image Sobel MSE Canny MSE CNN MSE NSCT Wavelet MSE Proposed method MSE Barbara 0.64891 0.4351 0.2574 0.7859 0.5874 0.1989 Boat 0.5945 0.3946 0.2415 0.6845 0.6842 0.1826 Lenna 0.84657 0.5874 0.4059 0.7548 0.3897 0.2834The results show that, compared with traditional edge detection methods, the proposed method possessed better edge locating ability and keeps more edge details.


## 4. Conclusions

We have presented an edge detection method based on nonsubsampled Contourlet transform and Gaussian mixture model. The proposed method is based on nonsubsampled Contourlet transform and Gaussian mixture model. It can capture more directional information and continuous edges. It overcomes the shortage of the limited directions of wavelet transform. The main contours in the image can be clearly expressed. This method can effectively locate the edge. The thick edges detected by the proposed method express the main edges in the images, which are with great convenience for the image recognition. The proposed method runs slowly; it needs to be improved in the future.



---

*Source: 4125909-2016-07-28.xml*