Canny边缘检测算法
1. 写在前面
最近在做边缘检测方面的一些工作,在网络上也找了很多有用的资料,感谢那些积极分享知识的先辈们,自己在理解Canny边缘检测算法的过程中也走了一些弯路,在编程实现的过程中,也遇到了一个让我怀疑人生的BUG(日了狗狗)。就此写下此文,作为后记,也希望此篇文章可以帮助那些在理解Canny算法的道路上暂入迷途的童鞋。废话少说,上干货。
2. Canny边缘检测算法的发展历史
Canny边缘检测于1986年由JOHN CANNY首次在论文《A Computational Approach to Edge Detection》中提出,就此拉开了Canny边缘检测算法的序幕。
Canny边缘检测是从不同视觉对象中提取有用的结构信息并大大减少要处理的数据量的一种技术,目前已广泛应用于各种计算机视觉系统。Canny发现,在不同视觉系统上对边缘检测的要求较为类似,因此,可以实现一种具有广泛应用意义的边缘检测技术。边缘检测的一般标准包括:
1) 以低的错误率检测边缘,也即意味着需要尽可能准确的捕获图像中尽可能多的边缘。
2) 检测到的边缘应精确定位在真实边缘的中心。
3) 图像中给定的边缘应只被标记一次,并且在可能的情况下,图像的噪声不应产生假的边缘。
为了满足这些要求,Canny使用了变分法。Canny检测器中的最优函数使用四个指数项的和来描述,它可以由高斯函数的一阶导数来近似。
在目前常用的边缘检测方法中,Canny边缘检测算法是具有严格定义的,可以提供良好可靠检测的方法之一。由于它具有满足边缘检测的三个标准和实现过程简单的优势,成为边缘检测最流行的算法之一。
3. Canny边缘检测算法的处理流程
Canny边缘检测算法可以分为以下5个步骤:
1) 使用高斯滤波器,以平滑图像,滤除噪声。
2) 计算图像中每个像素点的梯度强度和方向。
3) 应用非极大值(Non-Maximum Suppression)抑制,以消除边缘检测带来的杂散响应。
4) 应用双阈值(Double-Threshold)检测来确定真实的和潜在的边缘。
5) 通过抑制孤立的弱边缘最终完成边缘检测。
下面详细介绍每一步的实现思路。
3.1 高斯平滑滤波
为了尽可能减少噪声对边缘检测结果的影响,所以必须滤除噪声以防止由噪声引起的错误检测。为了平滑图像,使用高斯滤波器与图像进行卷积,该步骤将平滑图像,以减少边缘检测器上明显的噪声影响。大小为(2k+1)x(2k+1)的高斯滤波器核的生成方程式由下式给出:
下面是一个sigma = 1.4,尺寸为3x3的高斯卷积核的例子(需要注意归一化):
若图像中一个3x3的窗口为A,要滤波的像素点为e,则经过高斯滤波之后,像素点e的亮度值为:
其中*为卷积符号,sum表示矩阵中所有元素相加求和。
重要的是需要理解,高斯卷积核大小的选择将影响Canny检测器的性能。尺寸越大,检测器对噪声的敏感度越低,但是边缘检测的定位误差也将略有增加。一般5x5是一个比较不错的trade off。
3.2 计算梯度强度和方向
图像中的边缘可以指向各个方向,因此Canny算法使用四个算子来检测图像中的水平、垂直和对角边缘。边缘检测的算子(如Roberts,Prewitt,Sobel等)返回水平Gx和垂直Gy方向的一阶导数值,由此便可以确定像素点的梯度G和方向theta 。
其中G为梯度强度, theta表示梯度方向,arctan为反正切函数。下面以Sobel算子为例讲述如何计算梯度强度和方向。
x和y方向的Sobel算子分别为:
其中Sx表示x方向的Sobel算子,用于检测y方向的边缘; Sy表示y方向的Sobel算子,用于检测x方向的边缘(边缘方向和梯度方向垂直)。在直角坐标系中,Sobel算子的方向如下图所示。
图3-1 Sobel算子的方向
若图像中一个3x3的窗口为A,要计算梯度的像素点为e,则和Sobel算子进行卷积之后,像素点e在x和y方向的梯度值分别为:
其中*为卷积符号,sum表示矩阵中所有元素相加求和。根据公式(3-2)便可以计算出像素点e的梯度和方向。
3.3 非极大值抑制
非极大值抑制是一种边缘稀疏技术,非极大值抑制的作用在于“瘦”边。对图像进行梯度计算后,仅仅基于梯度值提取的边缘仍然很模糊。对于标准3,对边缘有且应当只有一个准确的响应。而非极大值抑制则可以帮助将局部最大值之外的所有梯度值抑制为0,对梯度图像中每个像素进行非极大值抑制的算法是:
1) 将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。
2) 如果当前像素的梯度强度与另外两个像素相比最大,则该像素点保留为边缘点,否则该像素点将被抑制。
通常为了更加精确的计算,在跨越梯度方向的两个相邻像素之间使用线性插值来得到要比较的像素梯度,现举例如下:
图3-2 梯度方向分割
如图3-2所示,将梯度分为8个方向,分别为E、NE、N、NW、W、SW、S、SE,其中0代表00~45o,1代表450~90o,2代表-900~-45o,3代表-450~0o。像素点P的梯度方向为theta,则像素点P1和P2的梯度线性插值为:
因此非极大值抑制的伪代码描写如下:
需要注意的是,如何标志方向并不重要,重要的是梯度方向的计算要和梯度算子的选取保持一致。
3.4 双阈值检测
在施加非极大值抑制之后,剩余的像素可以更准确地表示图像中的实际边缘。然而,仍然存在由于噪声和颜色变化引起的一些边缘像素。为了解决这些杂散响应,必须用弱梯度值过滤边缘像素,并保留具有高梯度值的边缘像素,可以通过选择高低阈值来实现。如果边缘像素的梯度值高于高阈值,则将其标记为强边缘像素;如果边缘像素的梯度值小于高阈值并且大于低阈值,则将其标记为弱边缘像素;如果边缘像素的梯度值小于低阈值,则会被抑制。阈值的选择取决于给定输入图像的内容。
双阈值检测的伪代码描写如下:
3.5 抑制孤立低阈值点
到目前为止,被划分为强边缘的像素点已经被确定为边缘,因为它们是从图像中的真实边缘中提取出来的。然而,对于弱边缘像素,将会有一些争论,因为这些像素可以从真实边缘提取也可以是因噪声或颜色变化引起的。为了获得准确的结果,应该抑制由后者引起的弱边缘。通常,由真实边缘引起的弱边缘像素将连接到强边缘像素,而噪声响应未连接。为了跟踪边缘连接,通过查看弱边缘像素及其8个邻域像素,只要其中一个为强边缘像素,则该弱边缘点就可以保留为真实的边缘。
抑制孤立边缘点的伪代码描述如下:
4 总结
通过以上5个步骤即可完成基于Canny算法的边缘提取,图5-1是该算法的检测效果图,希望对大家有所帮助。
图5-1 Canny边缘检测效果
附录是完整的MATLAB源码,可以直接拿来运行。
% -------------------------------------------------------------------------
% Edge Detection Using Canny Algorithm.
% Auther: Yongli Yan.
% Mail: yanyongli@ime.ac.cn
% Date: 2017.08.01.
% The direction of Sobel operator.
% ^(y)
% |
% |
% |
% 0--------->(x)
% Direction of Gradient:
% 3 2 1
% 0 P 0
% 1 2 3
% P = Current Point.
% NW N NE
% W P E
% SW S SE
% Point Index:
% f(x-1,y-1) f(x-1,y) f(x-1,y+1)
% f(x, y-1) f(x, y) f(x, y+1)
% f(x+1,y-1) f(x+1,y) f(x+1,y+1)
% Parameters:
% percentOfPixelsNotEdges: Used for selecting thresholds.
% thresholdRatio: Low thresh is this fraction of the high.
% -------------------------------------------------------------------------
function imgCanny = edge_canny(I,gaussDim,sigma,percentOfPixelsNotEdges,thresholdRatio)
%% Gaussian smoothing filter.
m = gaussDim(1);
n = gaussDim(2);
if mod(m,2) == 0 || mod(n,2) == 0error('The dimensionality of Gaussian must be odd!');
end
% Generate gaussian convolution kernel.
gaussKernel = fspecial('gaussian', [m,n], sigma);
% Image edge copy.
[m,n] = size(gaussKernel);
[row,col,dim] = size(I);
if dim > 1imgGray = rgb2gray(I);
elseimgGray = I;
end
imgCopy = imgReplicate(imgGray,(m-1)/2,(n-1)/2);
% Gaussian smoothing filter.
imgData = zeros(row,col);
for ii = 1:rowfor jj = 1:colwindow = imgCopy(ii:ii+m-1,jj:jj+n-1);GSF = window.*gaussKernel;imgData(ii,jj) = sum(GSF(:));end
end
%% Calculate the gradient values for each pixel.
% Sobel operator.
dgau2Dx = [-1 0 1;-2 0 2;-1 0 1];
dgau2Dy = [1 2 1;0 0 0;-1 -2 -1];
[m,n] = size(dgau2Dx);
% Image edge copy.
imgCopy = imgReplicate(imgData,(m-1)/2,(n-1)/2);
% To store the gradient and direction information.
gradx = zeros(row,col);
grady = zeros(row,col);
gradm = zeros(row,col);
dir = zeros(row,col); % Direction of gradient.
% Calculate the gradient values for each pixel.
for ii = 1:rowfor jj = 1:colwindow = imgCopy(ii:ii+m-1,jj:jj+n-1);dx = window.*dgau2Dx;dy = window.*dgau2Dy;dx = dx'; % Make the sum more accurate.dx = sum(dx(:));dy = sum(dy(:));gradx(ii,jj) = dx;grady(ii,jj) = dy;gradm(ii,jj) = sqrt(dx^2 + dy^2);% Calculate the angle of the gradient.theta = atand(dy/dx) + 90; % 0~180.% Determine the direction of the gradient.if (theta >= 0 && theta < 45)dir(ii,jj) = 2;elseif (theta >= 45 && theta < 90)dir(ii,jj) = 3;elseif (theta >= 90 && theta < 135)dir(ii,jj) = 0;elsedir(ii,jj) = 1;endend
end
% Normalize for threshold selection.
magMax = max(gradm(:));
if magMax ~= 0gradm = gradm / magMax;
end
%% Plot 3D gradient graph.
% [xx, yy] = meshgrid(1:col, 1:row);
% figure;
% surf(xx,yy,gradm);
%% Threshold selection.
counts = imhist(gradm, 64);
highThresh = find(cumsum(counts) > percentOfPixelsNotEdges*row*col,1,'first') / 64;
lowThresh = thresholdRatio*highThresh;
%% Non-Maxima Suppression(NMS) Using Linear Interpolation.
gradmCopy = zeros(row,col);
imgBW = zeros(row,col);
for ii = 2:row-1for jj = 2:col-1E = gradm(ii,jj+1);S = gradm(ii+1,jj);W = gradm(ii,jj-1);N = gradm(ii-1,jj);NE = gradm(ii-1,jj+1);NW = gradm(ii-1,jj-1);SW = gradm(ii+1,jj-1);SE = gradm(ii+1,jj+1);% Linear interpolation.% dy/dx = tan(theta).% dx/dy = tan(90-theta).gradValue = gradm(ii,jj);if dir(ii,jj) == 0d = abs(grady(ii,jj)/gradx(ii,jj));gradm1 = E*(1-d) + NE*d;gradm2 = W*(1-d) + SW*d;elseif dir(ii,jj) == 1d = abs(gradx(ii,jj)/grady(ii,jj));gradm1 = N*(1-d) + NE*d;gradm2 = S*(1-d) + SW*d;elseif dir(ii,jj) == 2d = abs(gradx(ii,jj)/grady(ii,jj));gradm1 = N*(1-d) + NW*d;gradm2 = S*(1-d) + SE*d;elseif dir(ii,jj) == 3d = abs(grady(ii,jj)/gradx(ii,jj));gradm1 = W*(1-d) + NW*d;gradm2 = E*(1-d) + SE*d;elsegradm1 = highThresh;gradm2 = highThresh;end% Non-Maxima Suppression.if gradValue >= gradm1 && gradValue >= gradm2if gradValue >= highThreshimgBW(ii,jj) = 1;gradmCopy(ii,jj) = highThresh;elseif gradValue >= lowThreshgradmCopy(ii,jj) = lowThresh;elsegradmCopy(ii,jj) = 0;endelsegradmCopy(ii,jj) = 0;endend
end
%% High-Low threshold detection.Double-Threshold.
% If the 8 pixels around the low threshold point have high threshold, then
% the low threshold pixel should be retained.
for ii = 2:row-1for jj = 2:col-1if gradmCopy(ii,jj) == lowThreshneighbors = [...gradmCopy(ii-1,jj-1), gradmCopy(ii-1,jj), gradmCopy(ii-1,jj+1),...gradmCopy(ii, jj-1), gradmCopy(ii, jj+1),...gradmCopy(ii+1,jj-1), gradmCopy(ii+1,jj), gradmCopy(ii+1,jj+1)...];if ~isempty(find(neighbors) == highThresh)imgBW(ii,jj) = 1;endendend
end
imgCanny = logical(imgBW);
end
%% Local functions. Image Replicate.
function imgRep = imgReplicate(I,rExt,cExt)
[row,col] = size(I);
imgCopy = zeros(row+2*rExt,col+2*cExt);
% 4 edges and 4 corners pixels.
top = I(1,:);
bottom = I(row,:);
left = I(:,1);
right = I(:,col);
topLeftCorner = I(1,1);
topRightCorner = I(1,col);
bottomLeftCorner = I(row,1);
bottomRightCorner = I(row,col);
% The coordinates of the oroginal image after the expansion in the new graph.
topLeftR = rExt+1;
topLeftC = cExt+1;
bottomLeftR = topLeftR+row-1;
bottomLeftC = topLeftC;
topRightR = topLeftR;
topRightC = topLeftC+col-1;
bottomRightR = topLeftR+row-1;
bottomRightC = topLeftC+col-1;
% Copy original image and 4 edges.
imgCopy(topLeftR:bottomLeftR,topLeftC:topRightC) = I;
imgCopy(1:rExt,topLeftC:topRightC) = repmat(top,[rExt,1]);
imgCopy(bottomLeftR+1:end,bottomLeftC:bottomRightC) = repmat(bottom,[rExt,1]);
imgCopy(topLeftR:bottomLeftR,1:cExt) = repmat(left,[1,cExt]);
imgCopy(topRightR:bottomRightR,topRightC+1:end) = repmat(right,[1,cExt]);
% Copy 4 corners.
for ii = 1:rExtfor jj = 1:cExtimgCopy(ii,jj) = topLeftCorner;imgCopy(ii,jj+topRightC) = topRightCorner;imgCopy(ii+bottomLeftR,jj) = bottomLeftCorner;imgCopy(ii+bottomRightR,jj+bottomRightC) = bottomRightCorner;end
end
imgRep = imgCopy;
end
%% End of file.
附錄2
%% test file of canny algorithm.
close all;
clear all;
clc;
%%
img = imread('./picture/lena.jpg');
% img = imnoise(img,'salt & pepper',0.01);
%%
[~,~,dim] = size(img);
if dim > 1imgCanny1 = edge(rgb2gray(img),'canny'); % system function.
elseimgCanny1 = edge(img,'canny'); % system function.
end
imgCanny2 = edge_canny(img,[5,5],1.4,0.9,0.4); % my function.
%%
figure;
subplot(2,2,1);
imshow(img);
title('original');
%%
subplot(2,2,3);
imshow(imgCanny1);
title('system function');
%%
subplot(2,2,4);
imshow(imgCanny2);
title('my function');
Canny边缘检测算法相关推荐
- OpenCV+python:Canny边缘检测算法
1,边缘处理 图像边缘信息主要集中在高频段,通常说图像锐化或检测边缘,实质就是高频滤波.我们知道微分运算是求信号的变化率,具有加强高频分量的作用. 在空域运算中来说,对图像的锐化就是计算微分.由于数字 ...
- 【算法随记一】Canny边缘检测算法实现和优化分析。
以前的博文大部分都写的非常详细,有很多分析过程,不过写起来确实很累人,一般一篇好的文章要整理个三四天,但是,时间越来越紧张,后续的一些算法可能就以随记的方式,把实现过程的一些比较容易出错和有价值的细节 ...
- python canny检测_【数字图像分析】基于Python实现 Canny Edge Detection(Canny 边缘检测算法)...
Canny 边缘检测算法 Steps: 高斯滤波平滑 计算梯度大小和方向 非极大值抑制 双阈值检测和连接 代码结构: Canny Edge Detection |Gaussian_Smoothing ...
- Canny边缘检测算法原理及其VC实现详解(一)
原文地址:http://blog.csdn.net/likezhaobin/article/details/6892176 图象的边缘是指图象局部区域亮度变化显著的部分,该区域的灰度剖面一般可以看作是 ...
- python canny算法_Python 实现 Canny 边缘检测算法
Canny 边缘检测算法由计算机科学家 John F. Canny 于 1986 年提出的.其不仅提供了算法,还带来了一套边缘检测的理论,分阶段的解释如何实现边缘检测.Canny 检测算法包含下面几个 ...
- 图像处理:Canny边缘检测算法原理(一)
图象的边缘是指图象局部区域亮度变化显著的部分,该区域的灰度剖面一般可以看作是一个阶跃,既从一个灰度值在很小的缓冲区域内急剧变化到另一个灰度相差较大的灰度值.图象的边缘部分集中了图象的大部分信息,图象边 ...
- opencv 图像边缘检测 Canny边缘检测算法使用
图解边缘检测 opencv 应用Canny算法进行边缘检测 import cv2 as cv import numpy as npimg = cv.imread('baby_g.jpg', 0) # ...
- 图像处理:推导Canny边缘检测算法
目录 概述 最优边缘检测 算法实现的步骤 1.灰度化与高斯滤波 2.计算图像的梯度和梯度方向 3.非极大值抑制 4.双阈值筛选边缘 5.利用滞后的边界跟踪 6.在图像中跟踪边缘 数学推导 Opencv ...
- Canny边缘检测算法(python 实现)
文章目录 最优边缘准则 算法实现步骤 1. 应用高斯滤波来平滑(模糊)图像,目的是去除噪声 2. 计算梯度强度和方向 3. 应用非最大抑制技术NMS来消除边误检 4. 应用双阈值的方法来决定可能的(潜 ...
最新文章
- hive olap 数据仓库_数据仓库那些事儿
- 解析 Java 类和对象的初始化过程 由一个单态模式引出的问题谈起
- DELPHI的DBGRID有两个难点
- python --> Python初阶 --> 基础语法 --> 条件和分支
- Zookeeper_安全认证讲解
- Failed to issue method call: Unit mysql.service failed to load: No such file or directory解决的方式...
- 微信小程序运行报错---invoke event
- android源码包下载
- 安徽省智慧政务新模式及典型应用
- FPGA内部硬件结构简介
- leetcode题库5-- 最长回文子串
- MD5校验工具的使用
- RT-Thread 入门学习笔记:使用虚拟GPS设备验证GPS nmea解析
- java实训文献_java实训论文参考文献写作指导
- Html5 Canvas 百行代码实现扫雷
- GAN网络学习笔记系列2-Cluster GAN
- 好听的摇滚_好听的摇滚歌曲有哪些 十大最好听中国摇滚歌曲
- 字符数组动态开辟空间和静态开辟空间
- 【显著性检测】基于HC算法实现图像显著性检测附MATLAB代码
- 微信小程序之数独挑战九宫格
热门文章
- 贪便宜买了减价香蕉之后
- volunteer is great
- 我妈妈的优点:做事情特别细致
- 【l转】VS2015下解决:无法解析的外部符号 __imp___vsnprintf 及__iob_func
- 数据蒋堂 | 怎样看待存储过程的移植困难
- Elasticsearch 简介
- redis的数据类型及设置方法
- Chrome插件会干坏事儿的
- 4)线性表[顺序表和链表]
- glance was not installed properly