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源碼，可以直接拿來運行。
附錄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');