[人工智能]改进Canny算法分析暨2021年11月13号笔记

????????最近迷上了图像边缘检测的研究，在这个机器学习横行的时代，还在专研图像处理的传统方法，可以说有那么点任性，哈哈。今天我们赏析一篇高引用的中文论文——《应用Otsu改进Canny算子的图像边缘检测方法》。

????????是的，这是一篇2014年发表于<计算机与数字工程>的论文，作者是顺德职业技术学院的张志强教授和中国科学院研究生院的宋海生教授，知网上引用次数高达48次。

? ? ? ? 也可以说这是一篇考古文(下次来个更老的)，不过我坚信，没有最好的方法，只有最合适用的地方。废话不多说，ACTION!

原文摘要：

????????针对传统Canny算子进行边缘检测时易丢失边缘细节的缺陷，提出了一种改进的Canny边缘 检测算法。首先从数学形态学角度分析传统Canny思想和缺陷；接着提出应用尺度自适应调整的高斯滤波器改进传统高斯平滑滤波；然后使用最大类间分差法（Otsu）替代传统双阈值检测算法求出最佳阈值，有效平衡去噪能力和保留边缘细节信息二者之间的矛盾；为解决Otsu遍历时间长，实时性差的缺点，提出应用Kalman滤波器最小误差估计算法预先估计后续图像序列的阈值范围，以加快遍历过程；最后进行实验分析比较，证明该算法可以更快速有效地检测出图像的边缘。

关键词：Canny；自适应滤波器；Otsu；Kalman

摘要赏析：

????????其实这个摘要确实就是摘要(不是废话么)，介绍的比较简洁。Canny算子又叫Canny算法，因为比较经典，所以带个“子”(比如，孔子，孟子，劳资。。)。在边缘检测过程中，Canny算子主要分为四个部分，分别是：

1，图像去噪。因为不清楚噪声的来源，所以一般按照随机噪声来处理，也就是常说的高斯滤波。2，方向和梯度的提取。在去噪以后我们认为图像已经干净了，所以根据像素值确定像素点的方向和相邻像素点的梯度。

3，非线性抑制。原始的Canny算子将图像分为若干个3*3区域(不够的边缘填充)，如果沿中心点梯度方向的插值比周围8个方向的像素大，那么这个中心点就有可能处于边缘上，也就是局部最大值。

4，双阈值边缘追踪。经过非线性抑制，可以说和边缘沾点边的像素点都被保留了下来，但是呢，有些不是。所以Canny算子玩了点花的，给图像设置了上下两个界限，只有同时满足a<x<b的条件，像素点才能被保留。这样一来，一些过明或者过暗的像素点就被筛掉了，完成了边缘点的提取。最后经过像素点的连通，就形成了边缘检测图像。

? ? ? ? 这篇文章呢，有三个部分进行了改进。

1，采用尺度自适应高斯滤波器替代了传统滤波器。

2，使用Otsu算法替代了传统凭经验获取双阈值的方法。

3，Otsu算法计算最佳阈值的过程是全部遍历的，使用Kalman滤波中的最小误差估计算法提前估计出图像的阈值范围，可以加快遍历的速度。

? ? ? ? 现在看起来这些改进不算什么，但是在2014年，图像处理领域刚刚蓬勃发展的时候，还是难能可贵的，也怪不得引用率这么高。

????????等等，你以为这期已经结束了么，嘿嘿，这期有代码实现福利！！！当然，是单独的

MTALAB实现：

Canny算法

clc
clear all
close all

img=imread('lena.jpg');%提取图像
Radius=5;
[ img ] = multiScaleSharpen( img, Radius);
img=rgb2gray(img);
sigma=sqrt(2);
img = im2single(img);
[m,n] = size(img);
% 用于输出的边界位图
e = false(m,n);

PercentOfPixelsNotEdges = .7; % Used for selecting thresholds
ThresholdRatio = .4;          % Low thresh is this fraction of the high.
    
% Calculate gradients using a derivative of Gaussian filter
%使用高斯滤波器的导数计算梯度
[dx, dy] = smoothGradient(img, sigma);

%计算梯度的大小
% Calculate Magnitude of Gradient
magGrad = hypot(dx, dy);
    
% Normalize for threshold selection
%用于阈值选择的规格化
magmax = max(magGrad(:));
if magmax > 0
   magGrad = magGrad / magmax;
end
   
% Determine Hysteresis Thresholds
counts=imhist(magGrad, 64);
highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,1,'first') / 64;
lowThresh = ThresholdRatio*highThresh;
%lowThresh =highThresh/3;
    
% Perform Non-Maximum Suppression Thining and Hysteresis Thresholding of Edge

[rstrong,cstrong] = find(magGrad>highThresh & E);
if ~isempty(rstrong) % result is all zeros if idxStrong is empty
rstrong = rem(idxStrong-1, m)+1;
cstrong = floor((idxStrong-1)/m)+1;
H = bwselect(E, cstrong, rstrong, 8);
else
       H = false(size(E));
end
thresh = [lowThresh highThresh];
end

figure(1)
imshow(H)

Otsu算法：

clc
clear all
close all

I1 = imread('lena.jpg');  %变为双精度，即0-1
I2 = rgb2gray(I1);
I3 = im2double(I2);

%%
% k2=graythresh(I3);              %得到最优阈值
tic
I4 = im2uint8(I3(:));
num_bins = 256;
counts = imhist(I3,num_bins);

% level = otsuthresh(counts);
%%  
num_bins = numel(counts);
counts = double( counts(:) );
p = counts / sum(counts);
omega = cumsum(p);
mu = cumsum(p .* (1:num_bins)');
mu_t = mu(end);
sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega));
%找到sigma_b_平方的最大值的位置。
%最大值可能会延伸到多个箱子，因此将这些位置平均起来。如果maxval是NaN，意味着sigma_b_平方都是NaN，那么返回0。
maxval = max(sigma_b_squared);

idx = mean(find(sigma_b_squared == maxval));
% Normalize the threshold to the range [0, 1].
t = (idx - 1) / (num_bins - 1);
%%
k=t*(255)+1;
toc
J1=im2bw(I3,t);                  %转换成二值图，k为分割阈值
figure(2),imshow(J1);

就这，白白