图像处理-边缘检测 Canny算子的实现 MATLAB&python

一.背景

Canny边缘检测算子是John F. Canny于1986年开发出来的一个多级边缘检测算法。它的目标是找到一个最优的边缘检测算法，最优边缘检测的含义是：①好的检测- 算法能够尽可能多地标识出图像中的实际边缘。②好的定位- 标识出的边缘要与实际图像中的实际边缘尽可能接近。③最小响应- 图像中的边缘只能标识一次，并且可能存在的图像噪声不应标识为边缘。

二.思路

1.Canny算子的基本步骤
①去噪，图像平滑。
任何边缘检测算法都不可能在未经处理的原始数据上很好地處理，一般采用的方法是利用原始数据与高斯核进行卷积，平滑处理。
②用一阶偏导的有限差分来计算梯度的幅值和方向。
计算一阶偏导用于找出灰度值变化较大的像素点，一般用Sobel算子。
③对梯度幅值进行非极大化抑制。
仅仅得到全局的梯度并不足以确定边缘，因此为确定边缘，必须保留局部梯度最大的点，而抑制非极大值。（non-maxima suppression,NMS）
解决方法：利用梯度的方向。
④边缘连接。
采用双阈值的方法进行边缘连接。

2.具体流程
对canny算子的四个步骤分别采用四个函数编写，最后利用主程序按顺序调用四个函数得到一个canny算子，并在每个过程中展示一次处理的结果，用于体现每个步骤的作用。
①图像平滑处理：
利用高斯滤波器对原始的图像，用高斯核函数：

先对x方向进行卷积，再对y方向进行卷积，得到平滑处理后的图像。

②Sobel边缘检测算子计算一阶偏导有限差分：
先将上一步的图像的上下左右各扩展一圈，为了能处理原始图像的边缘像素；然后设置Sobel算子的x和y方向的梯度模板，分别为：

将图像的每个像素分别与xy模板进行相乘叠加，得到该像素点的x,y梯度值，再将扩展的图像复原。最后，每个像素的梯度幅值为：x,y方向的梯度值得平方和开根号得到，梯度方向为反正切值（y梯度/x梯度），将得到得值输出。

③梯度幅值的非极大值抑制
利用下图的线性插值原理计算在8邻域中中间像素点是否为极大值。

线性插值的公式为：

其中，w为根据梯度的角度决定的权重大小。
由梯度角度的大小得到梯度方向的四分类：y梯度是否大于x梯度，以及，y梯度与x梯度符号是否相同。
由y梯度与x梯度的大小比较：
当y大于x时得到权值的计算公式为：绝对值（x梯度/y梯度），g2=左方梯度值，g4=右方梯度值；且当y与x同号，g1=左上角梯度值，g3=右下角梯度值；当y与x不同号，g1=左下方梯度值，g3=右上方梯度值。
当y小于x时得到权值的计算公式为：绝对值（y梯度/x梯度），g2=上方梯度值，g4=下方梯度值；且当y与x同号，g1=左上角梯度值，g3=右下角梯度值；当y与x不同号，g1=左下方梯度值，g3=右上方梯度值。
最后根据计算得到的g1234和权值w计算梯度方向上的两个值，与中间值进行比较，看是否对中间值进行抑制。

④边缘连接
通过双阈值算法检测，先选出大于最大阈值的梯度对应的像素作为边缘，再利用第二阈值，在大于第二阈值小于最大阈值的像素点中，根据邻域的极大值原理选择适当的像素，也作为边缘点，最终将边缘点全部连接起来，称为边缘的图像。

3.流程图
①主程序流程图为：

②高斯滤波平滑处理流程图为：

③Sobel算子流程图为：

④非极大值抑制的流程图为：

⑤边缘连接流程图为：

三.MATLAB代码

①主程序：cannytest.m

clear;clc;
image=imread('Lena.jpg');
image=rgb2gray(image);  %灰度化
%% 步骤一 添加高斯噪声并用高斯滤波平滑化
image_gaussnoise=imnoise(image,'gaussian',0,0.01); %添加高斯噪声
subplot(2,4,1);
imshow(image_gaussnoise);
title('添加高斯噪声');
image_guassfilter=gaussFilter(image_gaussnoise,2); %卷积核sigma值设为2，平滑效果较好，且不模糊
subplot(2,4,2);
imshow(image_guassfilter);
title('高斯滤波平滑处理');
 
%% 步骤二 求梯度
[grad,gradx,grady,angle]=Sobel_filter(image_guassfilter,3);
subplot(2,4,3);
imshow(gradx);
title('x方向梯度');
subplot(2,4,4);
imshow(grady);
title('y方向梯度');
subplot(2,4,5);
imshow(grad);
title('梯度');
 
%% 步骤三 梯度幅值非极大化抑制
nmsgrad=nms(grad,angle,gradx,grady);
subplot(2,4,6);
imshow(nmsgrad);
title('非极大值抑制后梯度');
 
%% 步骤四 双阈值算法检测和连接边缘
edgee=edge_correct(nmsgrad);
subplot(2,4,7);
imshow(edgee);
title('边缘连接后的canny算子最终结果');
 
%% 用matlab自带的canny算子函数实现边缘检测作为结果的对比
ed = edge(image_gaussnoise,'canny', 0.5);
subplot(2,4,8);
imshow(ed);
title('matlab自带的边缘检测函数');

②步骤一：高斯平滑:gaussFilter.m

%% 高斯平滑滤波函数（利用matlab自带的函数设计）
function output=gaussFilter(image,sigma)
output=image;
ksize=double(uint8(3*sigma)*2+1);%窗口大小一半为3*sigma 
window = fspecial('gaussian', [1,ksize], sigma); %使用1行ksize列的高斯核对图像先进行x方向卷积，再进行y方向卷积
for i = 1:size(image,3)
    ret = imfilter(image(:,:,i),window,'replicate');
    ret = imfilter(ret,window','replicate');
    output(:,:,i) = ret;
end 
end

③步骤二：sobel算子求梯度：Sobel_filter.m

%% Sobel算子计算梯度得到边缘图像
function [output,outputx,outputy,outputangle]=Sobel_filter(image,n)
image=double(image);
[h,w]=size(image);
%% 原灰度图像一圈扩展1个像素
big_image = zeros(h+2,w+2);
for i = 2:h+1
    for j = 2:w+1
        big_image(i,j) = image(i-1,j-1);
    end
end
%% 设置Sobel算子x,y方向模板
Hx = [-1,-2,-1;0,0,0;1,2,1];
Hy = Hx';
%% x,y方向梯度计算
gradx_image = zeros(h+2,w+2);
grady_image = zeros(h+2,w+2);
W = zeros(n,n);%n阶的卷积核
for i = 1:h
    for j = 1:w
        %模板移动窗口
        W = big_image(i:i+2,j:j+2);
        Sx = Hx .* W;
        Sy = Hy .* W;
        gradx_image(i+1,j+1) = sum(sum(Sx));
        grady_image(i+1,j+1) = sum(sum(Sy));
    end
end
%% 将一圈扩展1个像素的图像复原
gradx = zeros(h,w);
grady = zeros(h,w);
for i = 1:h
    for j = 1:w
        gradx(i,j) = gradx_image(i+1,j+1);
        grady(i,j) = grady_image(i+1,j+1);
    end
end
%sobel梯度
angle_grad=atand(grady./gradx); %求梯度方向(角度制）
grad = sqrt(gradx.^2 + grady.^2);  %得到图像的sobel梯度
outputx=uint8(gradx);
outputy=uint8(grady);
output=uint8(grad);
outputangle=angle_grad;
end

④步骤三：非极大值抑制处理：nms.m

%% 梯度幅值非极大值抑制
function output=nms(grad,gradangle,gradx,grady)
grad=double(grad);
[h,w]=size(grad);
%% 处理角度
sector=zeros(h,w);
canny1=zeros(h,w);
sector((gradangle>=45)&(gradangle<=-45))=0;%y方向梯度大于x梯度
sector((gradangle<45)&(gradangle>-45))=1;%y方向梯度小于x梯度
sector(gradangle>=0)=2;%x，y方向梯度符号相同
sector(gradangle<0)=3;%x，y方向梯度符号不同