【数字图像处理】Canny边缘检测详解及编程实现

原博客网址：http://blog.csdn.net/xiajun07061225/article/details/6926108

Canny边缘检测算法一直是边缘检测的经典算法。下面详细介绍Canny边缘检测算法的原理以及编程实现。

Canny边缘检测基本原理：
(1)图象边缘检测必须满足两个条件：一能有效地抑制噪声；二必须尽量精确确定边缘的位置。
 (2)根据对信噪比与定位乘积进行测度，得到最优化逼近算子。这就是Canny边缘检测算子。
 (3)类似与Marr（LoG）边缘检测方法，也属于先平滑后求导数的方法。

Canny 的目标是找到一个最优的边缘检测算法，最优边缘检测的含义是：
(1)好的检测 - 算法能够尽可能多地标识出图像中的实际边缘。
(2)好的定位 - 标识出的边缘要尽可能与实际图像中的实际边缘尽可能接近。
(3)最小响应 - 图像中的边缘只能标识一次，并且可能存在的图像雜訊不应标识为边缘。

Canny边缘检测算法的步骤：

(1)去噪

任何边缘检测算法都不可能在未经处理的原始数据上很好地處理，所以第一步是对原始数据与高斯 mask 作卷积，得到的图像与原始图像相比有些轻微的模糊（blurred）。这样，单独的一个像素雜訊在经过高斯平滑的图像上变得几乎没有影响。

(2)用一阶偏导的有限差分来计算梯度的幅值和方向。

(3)对梯度幅值进行非极大值抑制。

仅仅得到全局的梯度并不足以确定边缘，因此为确定边缘，必须保留局部梯度最大的点，而抑制非极大值。（non-maxima suppression,NMS）
解决方法：利用梯度的方向。

四个扇区的标号为0到3，对应3*3邻域的四种可能组合。在每一点上，邻域的中心象素M与沿着梯度线的两个象素相比。如果M的梯度值不比沿梯度线的两个相邻象素梯度值大，则令M=0。

(4)用双阈值算法检测和连接边缘。

减少假边缘段数量的典型方法是对N[i，j]使用一个阈值。将低于阈值的所有值赋零值。但问题是如何选取阈值？
 解决方法：双阈值算法。双阈值算法对非极大值抑制图象作用两个阈值τ1和τ2，且2τ1≈τ2，从而可以得到两个阈值边缘图象N1［i,j］和N2［i，j］。由于N2［i，j］使用高阈值得到，因而含有很少的假边缘，但有间断(不闭合)。双阈值法要在N2［i，j］中把边缘连接成轮廓，当到达轮廓的端点时，该算法就在N1［i,j］的8邻点位置寻找可以连接到轮廓上的边缘，这样，算法不断地在N1［i,j］中收集边缘，直到将N2［i,j］连接起来为止。

在连接边缘的时候，用数组模拟队列的实现。以进行8-连通域搜索。

更详细的资料请参考维基百科：http://zh.wikipedia.org/wiki/Canny%E7%AE%97%E5%AD%90

下面是我编程实现的Canny边缘检测代码，如有错误，请大家包涵、指正：

[cpp]  view plain copy print ?

1. I = imread('rice.png');  
2. I = double(I);  
3. [height,width] = size(I);  
4. J = I;  
5.   
6. conv = zeros(5,5);%高斯卷积核  
7. sigma = 1;%方差  
8. sigma_2 = sigma * sigma;%临时变量  
9. sum = 0;  
10. for i = 1:5  
11.     for j = 1:5  
12.         conv(i,j) = exp((-(i - 3) * (i - 3) - (j - 3) * (j - 3)) / (2 * sigma_2)) / (2 * 3.14 * sigma_2);%高斯公式  
13.         sum = sum + conv(i,j);  
14.     end  
15. end  
16. conv = conv./sum;%标准化  
17.   
18. %对图像实施高斯滤波  
19. for i = 1:height  
20.     for j = 1:width  
21.         sum = 0;%临时变量  
22.         for k = 1:5  
23.             for m = 1:5  
24.                 if (i - 3 + k) > 0 && (i - 3 + k) <= height && (j - 3 + m) > 0 && (j - 3 + m) < width  
25.                     sum = sum + conv(k,m) * I(i - 3 + k,j - 3 + m);  
26.                 end  
27.             end  
28.         end  
29.         J(i,j) = sum;  
30.     end  
31. end  
32. figure,imshow(J,[])  
33. title('高斯滤波后的结果')  
34. %求梯度  
35. dx = zeros(height,width);%x方向梯度  
36. dy = zeros(height,width);%y方向梯度  
37. d = zeros(height,width);  
38. for i = 1:height - 1  
39.     for j = 1:width - 1  
40.         dx(i,j) = J(i,j + 1) - J(i,j);  
41.         dy(i,j) = J(i + 1,j) - J(i,j);  
42.         d(i,j) = sqrt(dx(i,j) * dx(i,j) + dy(i,j) * dy(i,j));  
43.     end  
44. end  
45. figure,imshow(d,[])  
46. title('求梯度后的结果')  
47.   
48. %局部非极大值抑制  
49. K = d;%记录进行非极大值抑制后的梯度  
50. %设置图像边缘为不可能的边缘点  
51. for j = 1:width  
52.     K(1,j) = 0;  
53. end  
54. for j = 1:width  
55.     K(height,j) = 0;  
56. end  
57. for i = 2:width - 1  
58.     K(i,1) = 0;  
59. end  
60. for i = 2:width - 1  
61.     K(i,width) = 0;  
62. end  
63.   
64. for i = 2:height - 1  
65.     for j = 2:width - 1  
66.         %当前像素点的梯度值为0，则一定不是边缘点  
67.         if d(i,j) == 0  
68.             K(i,j) = 0;  
69.         else  
70.             gradX = dx(i,j);%当前点x方向导数  
71.             gradY = dy(i,j);%当前点y方向导数  
72.             gradTemp = d(i,j);%当前点梯度  
73.             %如果Y方向幅度值较大  
74.             if abs(gradY) > abs(gradX)  
75.                 weight = abs(gradX) / abs(gradY);%权重  
76.                 grad2 = d(i - 1,j);  
77.                 grad4 = d(i + 1,j);  
78.                 %如果x、y方向导数符号相同  
79.                 %像素点位置关系  
80.                 %g1 g2  
81.                 %   C  
82.                 %   g4 g3  
83.                 if gradX * gradY > 0  
84.                     grad1 = d(i - 1,j - 1);  
85.                     grad3 = d(i + 1,j + 1);  
86.                 else  
87.                     %如果x、y方向导数符号反  
88.                     %像素点位置关系  
89.                     %   g2 g1  
90.                     %   C  
91.                     %g3 g4  
92.                     grad1 = d(i - 1,j + 1);  
93.                     grad3 = d(i + 1,j - 1);  
94.                 end  
95.             %如果X方向幅度值较大  
96.             else  
97.                 weight = abs(gradY) / abs(gradX);%权重  
98.                 grad2 = d(i,j - 1);  
99.                 grad4 = d(i,j + 1);  
100.                 %如果x、y方向导数符号相同  
101.                 %像素点位置关系  
102.                 %g3  
103.                 %g4 C g2  
104.                 %     g1  
105.                 if gradX * gradY > 0  
106.                     grad1 = d(i + 1,j + 1);  
107.                     grad3 = d(i - 1,j - 1);  
108.                 else  
109.                     %如果x、y方向导数符号反  
110.                     %像素点位置关系  
111.                     %     g1  
112.                     %g4 C g2  
113.                     %g3  
114.                     grad1 = d(i - 1,j + 1);  
115.                     grad3 = d(i + 1,j - 1);  
116.                 end  
117.             end  
118.             %利用grad1-grad4对梯度进行插值  
119.             gradTemp1 = weight * grad1 + (1 - weight) * grad2;  
120.             gradTemp2 = weight * grad3 + (1 - weight) * grad4;  
121.             %当前像素的梯度是局部的最大值，可能是边缘点  
122.             if gradTemp >= gradTemp1 && gradTemp >= gradTemp2  
123.                 K(i,j) = gradTemp;  
124.             else  
125.                 %不可能是边缘点  
126.                 K(i,j) = 0;  
127.             end  
128.         end  
129.     end  
130. end  
131. figure,imshow(K,[])  
132. title('非极大值抑制后的结果')  
133.   
134. %定义双阈值：EP_MIN、EP_MAX，且EP_MAX = 2 * EP_MIN  
135. EP_MIN = 12;  
136. EP_MAX = EP_MIN * 2;  
137. EdgeLarge = zeros(height,width);%记录真边缘  
138. EdgeBetween = zeros(height,width);%记录可能的边缘点  
139. for i = 1:height  
140.     for j = 1:width  
141.         if K(i,j) >= EP_MAX%小于小阈值，不可能为边缘点  
142.             EdgeLarge(i,j) = K(i,j);  
143.         else if K(i,j) >= EP_MIN  
144.                 EdgeBetween(i,j) = K(i,j);  
145.             end  
146.         end  
147.     end  
148. end  
149. %把EdgeLarge的边缘连成连续的轮廓  
150. MAXSIZE = 999999;  
151. Queue = zeros(MAXSIZE,2);%用数组模拟队列  
152. front = 1;%队头  
153. rear = 1;%队尾  
154. edge = zeros(height,width);  
155. for i = 1:height  
156.     for j = 1:width  
157.         if EdgeLarge(i,j) > 0  
158.             %强点入队  
159.             Queue(rear,1) = i;  
160.             Queue(rear,2) = j;  
161.             rear = rear + 1;  
162.             edge(i,j) = EdgeLarge(i,j);  
163.             EdgeLarge(i,j) = 0;%避免重复计算  
164.         end  
165.         while front ~= rear%队不空  
166.             %队头出队  
167.             temp_i = Queue(front,1);  
168.             temp_j = Queue(front,2);  
169.             front = front + 1;  
170.             %8-连通域寻找可能的边缘点  
171.             %左上方  
172.             if EdgeBetween(temp_i - 1,temp_j - 1) > 0%把在强点周围的弱点变为强点  
173.                 EdgeLarge(temp_i - 1,temp_j - 1) = K(temp_i - 1,temp_j - 1);  
174.                 EdgeBetween(temp_i - 1,temp_j - 1) = 0;%避免重复计算  
175.                 %入队  
176.                 Queue(rear,1) = temp_i - 1;  
177.                 Queue(rear,2) = temp_j - 1;  
178.                 rear = rear + 1;  
179.             end  
180.             %正上方  
181.             if EdgeBetween(temp_i - 1,temp_j) > 0%把在强点周围的弱点变为强点  
182.                 EdgeLarge(temp_i - 1,temp_j) = K(temp_i - 1,temp_j);  
183.                 EdgeBetween(temp_i - 1,temp_j) = 0;  
184.                 %入队  
185.                 Queue(rear,1) = temp_i - 1;  
186.                 Queue(rear,2) = temp_j;  
187.                 rear = rear + 1;  
188.             end  
189.             %右上方  
190.             if EdgeBetween(temp_i - 1,temp_j + 1) > 0%把在强点周围的弱点变为强点  
191.                 EdgeLarge(temp_i - 1,temp_j + 1) = K(temp_i - 1,temp_j + 1);  
192.                 EdgeBetween(temp_i - 1,temp_j + 1) = 0;  
193.                 %入队  
194.                 Queue(rear,1) = temp_i - 1;  
195.                 Queue(rear,2) = temp_j + 1;  
196.                 rear = rear + 1;  
197.             end  
198.             %正左方  
199.             if EdgeBetween(temp_i,temp_j - 1) > 0%把在强点周围的弱点变为强点  
200.                 EdgeLarge(temp_i,temp_j - 1) = K(temp_i,temp_j - 1);  
201.                 EdgeBetween(temp_i,temp_j - 1) = 0;  
202.                 %入队  
203.                 Queue(rear,1) = temp_i;  
204.                 Queue(rear,2) = temp_j - 1;  
205.                 rear = rear + 1;  
206.             end  
207.             %正右方  
208.             if EdgeBetween(temp_i,temp_j + 1) > 0%把在强点周围的弱点变为强点  
209.                 EdgeLarge(temp_i,temp_j + 1) = K(temp_i,temp_j + 1);  
210.                 EdgeBetween(temp_i,temp_j + 1) = 0;  
211.                 %入队  
212.                 Queue(rear,1) = temp_i;  
213.                 Queue(rear,2) = temp_j + 1;  
214.                 rear = rear + 1;  
215.             end  
216.             %左下方  
217.             if EdgeBetween(temp_i + 1,temp_j - 1) > 0%把在强点周围的弱点变为强点  
218.                 EdgeLarge(temp_i + 1,temp_j - 1) = K(temp_i + 1,temp_j - 1);  
219.                 EdgeBetween(temp_i + 1,temp_j - 1) = 0;  
220.                 %入队  
221.                 Queue(rear,1) = temp_i + 1;  
222.                 Queue(rear,2) = temp_j - 1;  
223.                 rear = rear + 1;  
224.             end  
225.             %正下方  
226.             if EdgeBetween(temp_i + 1,temp_j) > 0%把在强点周围的弱点变为强点  
227.                 EdgeLarge(temp_i + 1,temp_j) = K(temp_i + 1,temp_j);  
228.                 EdgeBetween(temp_i + 1,temp_j) = 0;  
229.                 %入队  
230.                 Queue(rear,1) = temp_i + 1;  
231.                 Queue(rear,2) = temp_j;  
232.                 rear = rear + 1;  
233.             end  
234.             %右下方  
235.             if EdgeBetween(temp_i + 1,temp_j + 1) > 0%把在强点周围的弱点变为强点  
236.                 EdgeLarge(temp_i + 1,temp_j + 1) = K(temp_i + 1,temp_j + 1);  
237.                 EdgeBetween(temp_i + 1,temp_j + 1) = 0;  
238.                 %入队  
239.                 Queue(rear,1) = temp_i + 1;  
240.                 Queue(rear,2) = temp_j + 1;  
241.                 rear = rear + 1;  
242.             end  
243.         end  
244.         %下面2行用于观察程序运行的状况  
245.         i  
246.         j  
247.     end  
248. end  
249.   
250. figure,imshow(edge,[])  
251. title('双阈值后的结果')  

对图片rice.png进行处理后的 结果 如下：