Canny边缘检测-白红宇

Canny边缘检测

阅读量：5169 次

发布时间：2019-06-13

本文共 8564 字，大约阅读时间需要 28 分钟。

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-连通域搜索。

更详细的资料请参考维基百科：

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

[cpp]

I = imread('rice.png');

I = double(I);

[height,width] = size(I);

J = I;

conv = zeros(5,5);%高斯卷积核

sigma = 1;%方差

sigma_2 = sigma * sigma;%临时变量

sum = 0;

for i = 1:5

for j = 1:5

conv(i,j) = exp((-(i - 3) * (i - 3) - (j - 3) * (j - 3)) / (2 * sigma_2)) / (2 * 3.14 * sigma_2);%高斯公式

sum = sum + conv(i,j);

conv = conv./sum;%标准化

%对图像实施高斯滤波

for i = 1:height

for j = 1:width

sum = 0;%临时变量

for k = 1:5

for m = 1:5

if (i - 3 + k) > 0 && (i - 3 + k) <= height && (j - 3 + m) > 0 && (j - 3 + m) < width

sum = sum + conv(k,m) * I(i - 3 + k,j - 3 + m);

J(i,j) = sum;

figure,imshow(J,[])

title('高斯滤波后的结果')

%求梯度

dx = zeros(height,width);%x方向梯度

dy = zeros(height,width);%y方向梯度

d = zeros(height,width);

for i = 1:height - 1

for j = 1:width - 1

dx(i,j) = J(i,j + 1) - J(i,j);

dy(i,j) = J(i + 1,j) - J(i,j);

d(i,j) = sqrt(dx(i,j) * dx(i,j) + dy(i,j) * dy(i,j));

figure,imshow(d,[])

title('求梯度后的结果')

%局部非极大值抑制

K = d;%记录进行非极大值抑制后的梯度

%设置图像边缘为不可能的边缘点

for j = 1:width

K(1,j) = 0;

for j = 1:width

K(height,j) = 0;

for i = 2:width - 1

K(i,1) = 0;

for i = 2:width - 1

K(i,width) = 0;

for i = 2:height - 1

for j = 2:width - 1

%当前像素点的梯度值为0，则一定不是边缘点

if d(i,j) == 0

K(i,j) = 0;

else

gradX = dx(i,j);%当前点x方向导数

gradY = dy(i,j);%当前点y方向导数

gradTemp = d(i,j);%当前点梯度

%如果Y方向幅度值较大

if abs(gradY) > abs(gradX)

weight = abs(gradX) / abs(gradY);%权重

grad2 = d(i - 1,j);

grad4 = d(i + 1,j);

%如果x、y方向导数符号相同

%像素点位置关系

%g1 g2

% g4 g3

if gradX * gradY > 0

grad1 = d(i - 1,j - 1);

grad3 = d(i + 1,j + 1);

else

%如果x、y方向导数符号反

%像素点位置关系

% g2 g1

%g3 g4

grad1 = d(i - 1,j + 1);

grad3 = d(i + 1,j - 1);

%如果X方向幅度值较大

else

weight = abs(gradY) / abs(gradX);%权重

grad2 = d(i,j - 1);

grad4 = d(i,j + 1);

%如果x、y方向导数符号相同

%像素点位置关系

%g4 C g2

% g1

if gradX * gradY > 0

grad1 = d(i + 1,j + 1);

grad3 = d(i - 1,j - 1);

else

%如果x、y方向导数符号反

%像素点位置关系

% g1

%g4 C g2

grad1 = d(i - 1,j + 1);

grad3 = d(i + 1,j - 1);

%利用grad1-grad4对梯度进行插值

gradTemp1 = weight * grad1 + (1 - weight) * grad2;

gradTemp2 = weight * grad3 + (1 - weight) * grad4;

%当前像素的梯度是局部的最大值，可能是边缘点

if gradTemp >= gradTemp1 && gradTemp >= gradTemp2

K(i,j) = gradTemp;

else

%不可能是边缘点

K(i,j) = 0;

figure,imshow(K,[])

title('非极大值抑制后的结果')

%定义双阈值：EP_MIN、EP_MAX，且EP_MAX = 2 * EP_MIN

EP_MIN = 12;

EP_MAX = EP_MIN * 2;

EdgeLarge = zeros(height,width);%记录真边缘

EdgeBetween = zeros(height,width);%记录可能的边缘点

for i = 1:height

for j = 1:width

if K(i,j) >= EP_MAX%小于小阈值，不可能为边缘点

EdgeLarge(i,j) = K(i,j);

else if K(i,j) >= EP_MIN

EdgeBetween(i,j) = K(i,j);

%把EdgeLarge的边缘连成连续的轮廓

MAXSIZE = 999999;

Queue = zeros(MAXSIZE,2);%用数组模拟队列

front = 1;%队头

rear = 1;%队尾

edge = zeros(height,width);

for i = 1:height

for j = 1:width

if EdgeLarge(i,j) > 0

%强点入队

Queue(rear,1) = i;

Queue(rear,2) = j;

rear = rear + 1;

edge(i,j) = EdgeLarge(i,j);

EdgeLarge(i,j) = 0;%避免重复计算

while front ~= rear%队不空

%队头出队

temp_i = Queue(front,1);

temp_j = Queue(front,2);

front = front + 1;

%8-连通域寻找可能的边缘点

%左上方

if EdgeBetween(temp_i - 1,temp_j - 1) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i - 1,temp_j - 1) = K(temp_i - 1,temp_j - 1);

EdgeBetween(temp_i - 1,temp_j - 1) = 0;%避免重复计算

%入队

Queue(rear,1) = temp_i - 1;

Queue(rear,2) = temp_j - 1;

rear = rear + 1;

%正上方

if EdgeBetween(temp_i - 1,temp_j) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i - 1,temp_j) = K(temp_i - 1,temp_j);

EdgeBetween(temp_i - 1,temp_j) = 0;

%入队

Queue(rear,1) = temp_i - 1;

Queue(rear,2) = temp_j;

rear = rear + 1;

%右上方

if EdgeBetween(temp_i - 1,temp_j + 1) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i - 1,temp_j + 1) = K(temp_i - 1,temp_j + 1);

EdgeBetween(temp_i - 1,temp_j + 1) = 0;

%入队

Queue(rear,1) = temp_i - 1;

Queue(rear,2) = temp_j + 1;

rear = rear + 1;

%正左方

if EdgeBetween(temp_i,temp_j - 1) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i,temp_j - 1) = K(temp_i,temp_j - 1);

EdgeBetween(temp_i,temp_j - 1) = 0;

%入队

Queue(rear,1) = temp_i;

Queue(rear,2) = temp_j - 1;

rear = rear + 1;

%正右方

if EdgeBetween(temp_i,temp_j + 1) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i,temp_j + 1) = K(temp_i,temp_j + 1);

EdgeBetween(temp_i,temp_j + 1) = 0;

%入队

Queue(rear,1) = temp_i;

Queue(rear,2) = temp_j + 1;

rear = rear + 1;

%左下方

if EdgeBetween(temp_i + 1,temp_j - 1) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i + 1,temp_j - 1) = K(temp_i + 1,temp_j - 1);

EdgeBetween(temp_i + 1,temp_j - 1) = 0;

%入队

Queue(rear,1) = temp_i + 1;

Queue(rear,2) = temp_j - 1;

rear = rear + 1;

%正下方

if EdgeBetween(temp_i + 1,temp_j) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i + 1,temp_j) = K(temp_i + 1,temp_j);

EdgeBetween(temp_i + 1,temp_j) = 0;

%入队

Queue(rear,1) = temp_i + 1;

Queue(rear,2) = temp_j;

rear = rear + 1;

%右下方

if EdgeBetween(temp_i + 1,temp_j + 1) > 0%把在强点周围的弱点变为强点

EdgeLarge(temp_i + 1,temp_j + 1) = K(temp_i + 1,temp_j + 1);

EdgeBetween(temp_i + 1,temp_j + 1) = 0;

%入队

Queue(rear,1) = temp_i + 1;

Queue(rear,2) = temp_j + 1;

rear = rear + 1;

%下面2行用于观察程序运行的状况

figure,imshow(edge,[])

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)
     
      <=M2); %如果当前点比两边点都大置1                if isbigger           new_edge(i,j)=img(i,j);         end            endendfigure;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    endendfigure;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

每步运行效果：

原图