Canny算子–邊緣檢測[通俗易懂]

大家好，又見面了，我是你們的朋友全棧君。

Canny邊緣檢測算法的發展曆史

Canny邊緣檢測于1986年由JOHN CANNY首次在論文《A Computational Approach to Edge Detection》中提出，就此拉開了Canny邊緣檢測算法的序幕。

Canny邊緣檢測是從不同視覺對象中提取有用的結構資訊并大大減少要處理的資料量的一種技術，目前已廣泛應用于各種計算機視覺系統。Canny發現，在不同視覺系統上對邊緣檢測的要求較為類似，是以，可以實作一種具有廣泛應用意義的邊緣檢測技術。邊緣檢測的一般标準包括：

1) 以低的錯誤率檢測邊緣，也即意味着需要盡可能準确的捕獲圖像中盡可能多的邊緣。

2) 檢測到的邊緣應精确定位在真實邊緣的中心。

3) 圖像中給定的邊緣應隻被标記一次，并且在可能的情況下，圖像的噪聲不應産生假的邊緣。

為了滿足這些要求，Canny使用了變分法。Canny檢測器中的最優函數使用四個指數項的和來描述，它可以由高斯函數的一階導數來近似。

在目前常用的邊緣檢測方法中，Canny邊緣檢測算法是具有嚴格定義的，可以提供良好可靠檢測的方法之一。由于它具有滿足邊緣檢測的三個标準和實作過程簡單的優勢，成為邊緣檢測最流行的算法之一。

Canny邊緣檢測算法的處理流程

Canny邊緣檢測算法可以分為以下5個步驟：

1) 使用高斯濾波器，以平滑圖像，濾除噪聲。

2) 計算圖像中每個像素點的梯度強度和方向。

3) 應用非極大值（Non-Maximum Suppression）抑制，以消除邊緣檢測帶來的雜散響應。

4) 應用雙門檻值（Double-Threshold）檢測來确定真實的和潛在的邊緣。

5) 通過抑制孤立的弱邊緣最終完成邊緣檢測。

下面詳細介紹每一步的實作思路。

高斯平滑濾波

為了盡可能減少噪聲對邊緣檢測結果的影響，是以必須濾除噪聲以防止由噪聲引起的錯誤檢測。為了平滑圖像，使用高斯濾波器與圖像進行卷積，該步驟将平滑圖像，以減少邊緣檢測器上明顯的噪聲影響。大小為(2k+1)x(2k+1)的高斯濾波器核的生成方程式由下式給出：

下面是一個sigma = 1.4，尺寸為3×3的高斯卷積核的例子（需要注意歸一化）：

若圖像中一個3×3的視窗為A，要濾波的像素點為e，則經過高斯濾波之後，像素點e的亮度值為：

其中*為卷積符号，sum表示矩陣中所有元素相加求和。

重要的是需要了解，高斯卷積核大小的選擇将影響Canny檢測器的性能。尺寸越大，檢測器對噪聲的敏感度越低，但是邊緣檢測的定位誤差也将略有增加。一般5×5是一個比較不錯的trade off。

計算梯度強度和方向

圖像中的邊緣可以指向各個方向，是以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算子的方向

若圖像中一個3×3的視窗為A，要計算梯度的像素點為e，則和Sobel算子進行卷積之後，像素點e在x和y方向的梯度值分别為：

其中*為卷積符号，sum表示矩陣中所有元素相加求和。根據公式（3-2）便可以計算出像素點e的梯度和方向。

非極大值抑制

非極大值抑制是一種邊緣稀疏技術，非極大值抑制的作用在于“瘦”邊。對圖像進行梯度計算後，僅僅基于梯度值提取的邊緣仍然很模糊。對于标準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.5 抑制孤立低門檻值點

到目前為止，被劃分為強邊緣的像素點已經被确定為邊緣，因為它們是從圖像中的真實邊緣中提取出來的。然而，對于弱邊緣像素，将會有一些争論，因為這些像素可以從真實邊緣提取也可以是因噪聲或顔色變化引起的。為了獲得準确的結果，應該抑制由後者引起的弱邊緣。通常，由真實邊緣引起的弱邊緣像素将連接配接到強邊緣像素，而噪聲響應未連接配接。為了跟蹤邊緣連接配接，通過檢視弱邊緣像素及其8個鄰域像素，隻要其中一個為強邊緣像素，則該弱邊緣點就可以保留為真實的邊緣。

抑制孤立邊緣點的僞代碼描述如下：

總結

通過以上5個步驟即可完成基于Canny算法的邊緣提取，圖5-1是該算法的檢測效果圖，希望對大家有所幫助。

　　　　　　　　　　　　　　　　　　　　　　　　圖5-1 Canny邊緣檢測效果

附錄是完整的MATLAB源碼，可以直接拿來運作。

附錄1

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

% -------------------------------------------------------------------------

% Edge Detection Using Canny Algorithm.

% Auther: Yongli Yan.

% Mail: [email protected]

% 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

%% 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

% 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

%% 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

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

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');

原部落格轉載自https://www.cnblogs.com/techyan1990/p/7291771.html

釋出者：全棧程式員棧長，轉載請注明出處：https://javaforall.cn/143786.html原文連結：https://javaforall.cn