Canny边缘检测

Posted emg0818

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Canny边缘检测相关的知识,希望对你有一定的参考价值。

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进行处理后的结果如下:

 

canny边缘检测一共四个部分:

  1.对原图像高斯平滑

  2.对高斯平滑后的图像进行sobel边缘检测。这里需要求横的和竖的还有联合的,所以一共三个需要sobel边缘检测图像。

  3.对联合的sobel检测图像进行非极大抑制

  4.连接边缘点并进行滞后阈值处理。

下面是代码:

main.m

复制代码
clear all;
close all;
clc;

img=imread(\'lena.jpg\');
imshow(img);
[m n]=size(img);
img=double(img);

%%canny边缘检测的前两步相对不复杂,所以我就直接调用系统函数了
%%高斯滤波
w=fspecial(\'gaussian\',[5 5]);
img=imfilter(img,w,\'replicate\');
figure;
imshow(uint8(img))

%%sobel边缘检测
w=fspecial(\'sobel\');
img_w=imfilter(img,w,\'replicate\');      %求横边缘
w=w\';
img_h=imfilter(img,w,\'replicate\');      %求竖边缘
img=sqrt(img_w.^2+img_h.^2);        %注意这里不是简单的求平均,而是平方和在开方。我曾经好长一段时间都搞错了
figure;
imshow(uint8(img))

%%下面是非极大抑制
new_edge=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        Mx=img_w(i,j);
        My=img_h(i,j);
        
        if My~=0
            o=atan(Mx/My);      %边缘的法线弧度
        elseif My==0 && Mx>0
            o=pi/2;
        else
            o=-pi/2;            
        end
        
        %Mx处用My和img进行插值
        adds=get_coords(o);      %边缘像素法线一侧求得的两点坐标,插值需要       
        M1=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %插值后得到的像素,用此像素和当前像素比较 
        adds=get_coords(o+pi);  %边缘法线另一侧求得的两点坐标,插值需要
        M2=My*img(i+adds(2),j+adds(1))+(Mx-My)*img(i+adds(4),j+adds(3));   %另一侧插值得到的像素,同样和当前像素比较
        
        isbigger=(Mx*img(i,j)>M1)*(Mx*img(i,j)>=M2)+(Mx*img(i,j)<M1)*(Mx*img(i,j)<=M2); %如果当前点比两边点都大置1
        
        if isbigger
           new_edge(i,j)=img(i,j); 
        end        
    end
end
figure;
imshow(uint8(new_edge))

%%下面是滞后阈值处理
up=120;     %上阈值
low=100;    %下阈值
set(0,\'RecursionLimit\',10000);  %设置最大递归深度
for i=1:m
    for j=1:n
      if new_edge(i,j)>up &&new_edge(i,j)~=255  %判断上阈值
            new_edge(i,j)=255;
            new_edge=connect(new_edge,i,j,low);
      end
    end
end
figure;
imshow(new_edge==255)
复制代码

get_coords.m

复制代码
function re=get_coords(angle)       %angle是边缘法线角度,返回法线前后两点
    sigma=0.000000001;
    x1=ceil(cos(angle+pi/8)*sqrt(2)-0.5-sigma);
    y1=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);
    x2=ceil(cos(angle-pi/8)*sqrt(2)-0.5-sigma);
    y2=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);
    re=[x1 y1 x2 y2];

end
复制代码

connect.m

复制代码
function nedge=connect(nedge,y,x,low)       %种子定位后的连通分析
    neighbour=[-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];  %八连通搜寻
    [m n]=size(nedge);
    for k=1:8
        yy=y+neighbour(k,1);
        xx=x+neighbour(k,2);
        if yy>=1 &&yy<=m &&xx>=1 && xx<=n  
            if nedge(yy,xx)>=low && nedge(yy,xx)~=255   %判断下阈值
                nedge(yy,xx)=255;
                nedge=connect(nedge,yy,xx,low);
            end
        end        
    end 

end
复制代码

每步运行效果:

原图

高斯模糊后

sobel边缘检测后

非极大抑制后

上阈值120,下阈值100检测结果。

以上是关于Canny边缘检测的主要内容,如果未能解决你的问题,请参考以下文章

sobel边缘检测优缺点与canny算子的优缺点?

Canny边缘检测

sobel边缘检测优缺点与canny算子的优缺点?

OpenCV入门指南第三篇Canny边缘检测

Python,Opencv cv2.Canny()边缘检测

一些关于Canny边缘检测算法的改进