邊緣檢測之Canny
1. 寫在前面算法
最近在作邊緣檢測方面的一些工做，在網絡上也找了不少有用的資料，感謝那些積極分享知識的先輩們，本身在理解Canny邊緣檢測算法的過程當中也走了一些彎路，在編程實現的過程當中，也遇到了一個讓我懷疑人生的BUG（日了狗狗）。就此寫下此文，做爲後記，也但願此篇文章能夠幫助那些在理解Canny算法的道路上暫入迷途的童鞋。廢話少說，上乾貨。編程
2. Canny邊緣檢測算法的發展歷史網絡
Canny邊緣檢測於1986年由JOHN CANNY首次在論文《A Computational Approach to Edge Detection》中提出，就此拉開了Canny邊緣檢測算法的序幕。函數
Canny邊緣檢測是從不一樣視覺對象中提取有用的結構信息並大大減小要處理的數據量的一種技術，目前已普遍應用於各類計算機視覺系統。Canny發現，在不一樣視覺系統上對邊緣檢測的要求較爲相似，所以，能夠實現一種具備普遍應用意義的邊緣檢測技術。邊緣檢測的通常標準包括：性能
1) 以低的錯誤率檢測邊緣，也即意味着須要儘量準確的捕獲圖像中儘量多的邊緣。this
2) 檢測到的邊緣應精肯定位在真實邊緣的中心。spa
3) 圖像中給定的邊緣應只被標記一次，而且在可能的狀況下，圖像的噪聲不該產生假的邊緣。orm
爲了知足這些要求，Canny使用了變分法。Canny檢測器中的最優函數使用四個指數項的和來描述，它能夠由高斯函數的一階導數來近似。對象
在目前經常使用的邊緣檢測方法中，Canny邊緣檢測算法是具備嚴格定義的，能夠提供良好可靠檢測的方法之一。因爲它具備知足邊緣檢測的三個標準和實現過程簡單的優點，成爲邊緣檢測最流行的算法之一。blog
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源碼，能夠直接拿來運行。
附錄1
% -------------------------------------------------------------------------
% 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) == 0
error('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 > 1
imgGray = rgb2gray(I);
else
imgGray = I;
end
imgCopy = imgReplicate(imgGray,(m-1)/2,(n-1)/2);
% Gaussian smoothing filter.
imgData = zeros(row,col);
for ii = 1:row
for jj = 1:col
window = 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:row
for jj = 1:col
window = 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;
else
dir(ii,jj) = 1;
end
end
end
% Normalize for threshold selection.
magMax = max(gradm(:));
if magMax ~= 0
gradm = 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-1
for jj = 2:col-1
E = 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) == 0
d = abs(grady(ii,jj)/gradx(ii,jj));
gradm1 = E*(1-d) + NE*d;
gradm2 = W*(1-d) + SW*d;
elseif dir(ii,jj) == 1
d = abs(gradx(ii,jj)/grady(ii,jj));
gradm1 = N*(1-d) + NE*d;
gradm2 = S*(1-d) + SW*d;
elseif dir(ii,jj) == 2
d = abs(gradx(ii,jj)/grady(ii,jj));
gradm1 = N*(1-d) + NW*d;
gradm2 = S*(1-d) + SE*d;
elseif dir(ii,jj) == 3
d = abs(grady(ii,jj)/gradx(ii,jj));
gradm1 = W*(1-d) + NW*d;
gradm2 = E*(1-d) + SE*d;
else
gradm1 = highThresh;
gradm2 = highThresh;
end
% Non-Maxima Suppression.
if gradValue >= gradm1 && gradValue >= gradm2
if gradValue >= highThresh
imgBW(ii,jj) = 1;
gradmCopy(ii,jj) = highThresh;
elseif gradValue >= lowThresh
gradmCopy(ii,jj) = lowThresh;
else
gradmCopy(ii,jj) = 0;
end
else
gradmCopy(ii,jj) = 0;
end
end
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-1
for jj = 2:col-1
if gradmCopy(ii,jj) == lowThresh
neighbors = [...
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;
end
end
end
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:rExt
for jj = 1:cExt
imgCopy(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 > 1
imgCanny1 = edge(rgb2gray(img),'canny'); % system function.
else
imgCanny1 = 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');
- 1. 邊緣檢測之Canny
- 2. Canny邊緣檢測
- 3. Canny邊緣檢測原理
- 4. Canny算法 邊緣檢測
- 5. Canny邊緣檢測算法
- 6. Canny 邊緣檢測算法
- 7. 邊緣檢測---Canny算子
- 8. OpenCV—python 邊緣檢測（Canny）
- 9. 邊緣檢測算法--------canny
- 10. OpenCV——邊緣檢測入門、Canny邊緣檢測
- 更多相關文章...
- • 使用TCP協議檢測防火牆 - TCP/IP教程
- • ping命令檢測遠程主機是否可用 - TCP/IP教程
- • 互聯網組織的未來：剖析GitHub員工的任性之源
- • Github 簡明教程
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. 邊緣檢測之Canny
- 2. Canny邊緣檢測
- 3. Canny邊緣檢測原理
- 4. Canny算法 邊緣檢測
- 5. Canny邊緣檢測算法
- 6. Canny 邊緣檢測算法
- 7. 邊緣檢測---Canny算子
- 8. OpenCV—python 邊緣檢測（Canny）
- 9. 邊緣檢測算法--------canny
- 10. OpenCV——邊緣檢測入門、Canny邊緣檢測