算法解剖系列-Canny边缘检测原理及实现「建议收藏」

2022-09-02 11:01:04 浏览数 (1)

大家好,又见面了,我是你们的朋友全栈君。

基本原理

  • 须满足条件:抑制噪声;精确定位边缘。
  • 从数学上表达了三个准则[信噪比准则(低错误率)、定位精度准则、单边缘响应准则],并寻找表达式的最佳解。
  • 属于先平滑后求导的方法。

算法基本步骤

1 、 1、 1、使用高斯滤波平滑图像

令 f ( x , y ) 表 示 数 据 ( 输 入 源 数 据 ) , G ( x , y ) 表 示 二 维 高 斯 函 数 ( 卷 积 操 作 数 ) , f s ( x , y ) 令f(x,y)表示数据(输入源数据),G(x,y)表示二维高斯函数(卷积操作数),f_s(x,y) 令f(x,y)表示数据(输入源数据),G(x,y)表示二维高斯函数(卷积操作数),fs​(x,y)为卷积平滑后的图像。 G ( x , y ) = 1 2 π σ 2 e − ( x 2 y 2 ) 2 σ 2 G(x,y)=frac{1}{2pisigma^2}e^frac{-(x^2 y^2)}{2sigma^2} G(x,y)=2πσ21​e2σ2−(x2 y2)​ f s ( x , y ) = f ( x , y ) ∗ G ( x , y ) f_s(x,y)=f(x,y)ast G(x,y) fs​(x,y)=f(x,y)∗G(x,y)

G u e s s 过 程 Guess过程 Guess过程
  • 用坐标点 ( x , y ) 表 示 一 个 3 × 3 的 邻 域 , 设 中 心 点 的 坐 标 为 ( 0 , 0 ) (x,y)表示一个3times3的邻域,设中心点的坐标为(0,0) (x,y)表示一个3×3的邻域,设中心点的坐标为(0,0),相邻的点以此类推。
  • 计算权重矩阵。设定方差 σ 2 = 0.64 sigma^2=0.64 σ2=0.64的值,将对应各个坐标点 ( x , y ) (x,y) (x,y)带入二维高斯公式 G ( x , y ) G(x,y) G(x,y)中,得到一个权重矩阵,归一化权重矩阵(矩阵中各个点除以权重之和),得到标准的权重矩阵,即高斯模板。
  • 计算高斯模糊。设在一幅图像中的 3 × 3 3times3 3×3区域内,用各像素点的灰度值乘以对应点的权重。

– 将得到的9个值求和,就是中心点的高斯模糊值。

**具体过程**:

  • 简单来说就是使用 G u e s s Guess Guess模板在原始图像中进行移位、相乘、相加的过程。

2 、 2、 2、计算幅值图像,角度图像 补充:求变化率时,对于一元函数,即求导;对于二元函数,求偏导。

数字图像处理中,用一阶有限差分近似求取灰度值的梯度值(变化率)。 (即:使差商( Δ f / Δ x Delta {f} / Delta {x} Δf/Δx) 近似取代微商( ∂ f / ∂ x partial{f} / partial{x} ∂f/∂x)。求灰度的变化率,分别取x和y方向上相邻像素做差,代替求取x和y方向一阶偏导)

例:计算一点 x 方 向 和 y x方向和y x方向和y方向的梯度幅值和方向

上图中显示一段直的边缘线段放大后一部分,每个方块代表一个像素点,用一个方框强调点处边缘的幅值和方向。令灰色像素值为0,白色像素值为1.

如图关于一点为中心的 3 × 3 3 times 3 3×3 邻域,使用 P r e w i t t Prewitt Prewitt卷积模板进行计算:

p x = [ − 1 − 1 − 1 0 0 0 1 1 1 ] , p y = [ − 1 0 1 − 1 0 1 − 1 0 1 ] p_x= begin{bmatrix} -1&-1&-1\ 0&0&0\ 1&1&1 end{bmatrix} , p_y= begin{bmatrix} -1&0&1\ -1&0&1\ -1&0&1 end{bmatrix} px​=⎣⎡​−101​−101​−101​⎦⎤​,py​=⎣⎡​−1−1−1​000​111​⎦⎤​ 根据 x 方 向 和 y x方向和y x方向和y方向的卷积模板,可知,在 3 × 3 3 times 3 3×3邻域中从底部一行像素值减去顶部一行的像素,得到 x x x方向的偏导数(梯度);同样,从右边一列像素值减去左边一列的像素,得到 y y y方向的偏导数。

x x x方向的梯度: g x = ∂ f / ∂ x = ( 0 − 0 ) ( 0 − 1 ) ( 0 − 1 ) = − 2 g_x=partial{f}/partial{x}=(0-0) (0-1) (0-1)=-2 gx​=∂f/∂x=(0−0) (0−1) (0−1)=−2 y y y方向的梯度: g y = ∂ f / ∂ x = ( 1 − 0 ) ( 1 − 0 ) ( 0 − 0 ) = 2 g_y=partial{f}/partial{x}=(1-0) (1-0) (0-0)=2 gy​=∂f/∂x=(1−0) (1−0) (0−0)=2

由此,可以得到该点梯度的幅值和方向: M ( x , y ) = g x 2 g y 2 = 2 2 M(x,y)=sqrt{g_x^2 g_y^2}=2sqrt{2} M(x,y)=gx2​ gy2​ ​=22 ​ α ( x , y ) = a r c t a n [ g y / g x ] = 135 ° alpha(x,y)=arctan[g_y/g_x]=135° α(x,y)=arctan[gy​/gx​]=135°

如下图表示了中心点的梯度向量、方位角以及边缘方向。(任一点的边缘与梯度向量正交)

注意

  • g x , g y , M ( x , y ) 和 α ( x , y ) g_x,g_y,M(x,y)和alpha(x,y) gx​,gy​,M(x,y)和α(x,y)都是与原图像大小相同的图像。
  • α ( x , y ) alpha(x,y) α(x,y)是为了确定梯度的方向,以及确定非极大值抑制时的对比方向;
  • M ( x , y ) M(x,y) M(x,y)用来确定是否为边缘点,幅值越大,表示像素区域灰度值变化越明显,更能代表边缘点.

Canny 算子卷积模板: p x = [ − 1 − 1 1 1 ] , p y = [ 1 − 1 1 − 1 ] p_x= begin{bmatrix} -1&-1\ 1&1 end{bmatrix} , p_y= begin{bmatrix} 1&-1\ 1&-1 end{bmatrix} px​=[−11​−11​],py​=[11​−1−1​]

p [ i , j ] = ( f [ i , j 1 ] − f [ i , j ] f [ i 1 , j 1 ] − f [ i 1 , j ] ) / 2 p[i,j]=(f[i,j 1]-f[i,j] f[i 1,j 1]-f[i 1,j])/2 p[i,j]=(f[i,j 1]−f[i,j] f[i 1,j 1]−f[i 1,j])/2 Q [ i , j ] = ( f [ i 1 , j ] − f [ i , j ] f [ i 1 , j 1 ] − f [ i , j 1 ] ) / 2 Q[i,j]=(f[i 1,j]-f[i,j] f[i 1,j 1]-f[i,j 1])/2 Q[i,j]=(f[i 1,j]−f[i,j] f[i 1,j 1]−f[i,j 1])/2 M [ i , j ] = p [ i , j ] 2 Q [ i , j ] 2 M[i,j]=sqrt{p[i,j]^2 Q[i,j]^2} M[i,j]=p[i,j]2 Q[i,j]2 ​ θ = a r c t a n ( Q [ i , j ] / p [ i , j ] ) theta=arctan(Q[i,j]/p[i,j]) θ=arctan(Q[i,j]/p[i,j])

3 、 3、 3、对幅值图像进行非极大值抑制

代码语言:javascript复制
步骤:
1. 对角度在四个方向进行划分,形成新的角度图;
2. 根据角度图(表征着梯度的方向),对幅值进行非极大值抑制。
  • 首先将角度划分成四个方向范围:水平 ( 0 ° ) 、 − 45 ° 、 垂 直 ( 90 ° ) 、 45 ° (0°)、-45°、垂直(90°)、 45° (0°)、−45°、垂直(90°)、 45°。如下图:

– 接着讨论对3×3区域的四个基本边缘方向进行非极大值抑制。

做法:若中心点(即:访问点)在沿其方向上邻域的梯度幅值最大,则保留;否则,抑制。

4 、 4、 4、双阈值检测和连接边缘

  • 选取高阈值 T H 和 低 阈 值 T L , 比 率 为 2 : 1 或 3 : 1 。 ( 一 般 取 T H = 0.3 或 0.2 , T L = 0.1 T_H和低阈值T_L,比率为2:1或3:1。(一般取T_H=0.3或0.2,T_L=0.1 TH​和低阈值TL​,比率为2:1或3:1。(一般取TH​=0.3或0.2,TL​=0.1)
  • 取出非极大值抑制后的图像中的最大梯度幅值,重新定义高低阈值。即: T H × M a x , T L × M a x T_Htimes{Max},T_Ltimes{Max} TH​×Max,TL​×Max。(当然可以自己给定)
  • 将 小 于 T L 的 点 抛 弃 , 赋 0 ; 将 大 于 T H 的 点 立 即 标 记 ( 这 些 点 就 是 边 缘 点 ) , 赋 1 。 将小于T_L的点抛弃,赋0;将大于T_H的点立即标记(这些点就是边缘点),赋1。 将小于TL​的点抛弃,赋0;将大于TH​的点立即标记(这些点就是边缘点),赋1。
  • 将 大 于 T L , 小 于 T H 的 点 使 用 8 连 通 区 域 确 定 ( 即 : 只 有 与 T H 像 素 连 接 时 才 会 被 接 受 , 成 为 边 缘 点 , 赋 1 ) 将大于T_L,小于T_H的点使用8连通区域确定(即:只有与T_H像素连接时才会被接受,成为边缘点,赋1) 将大于TL​,小于TH​的点使用8连通区域确定(即:只有与TH​像素连接时才会被接受,成为边缘点,赋1)。

注意:双阈值做法是将候选像素点拼接成轮廓,轮廓的形成时对这些像素运用滞后性阈值。

算法实现

Matlab代码
代码语言:javascript复制
clear all;
clc;

I = imread('rice.png');%读图
% I = rgb2gray(I);%灰度转换
I = double(I);%转化为双精度
[H,W] = size(I);%获取图像大小

%%  Step1:使用高斯滤波平滑图像

B = [1 2 1;2 4 2;1 2 1];%高斯滤波系数
B = 1/16.*B;%高斯滤波模板 方差=0.8
A = conv2(I,B,'same');%使用高斯模板进行卷积.计算二维卷积,结果与原图像大小相同 

%%  Step2:计算梯度的幅值图像,角度图像.

%Prewitt梯度模板
dx = [-1 0 1;-1 0 1;-1 0 1];%x方向的梯度模板
dy = [1 1 1; 0 0 0;-1 -1 -1];%y方向的梯度模板
gx = conv2(A,dx,'same');%获取x方向的梯度图像.使用梯度模板进行二维卷积,结果与原图像大小相同
gy = conv2(A,dy,'same');%获取y方向的梯度图像.使用梯度模板进行二维卷积,结果与原图像大小相同
M = sqrt((gx.^2)   (gy.^2));%获取幅值图像.大小与原图像相等.(.^)表示数组乘方
a = atan2(gy,gx);%获取弧度,范围:-pi~pi
a = a*180/pi;%将弧度转换为角度,得到角度图像,与原图像大小相等.

%%  Step3:对幅值图像进行应用非极大值抑制

%首先将角度划分成四个方向范围:水平(0°)、-45°、垂直(90°)、 45°
for i = 1:H
    for j = 1:W
        if((a(i,j) >= -22.5) && (a(i,j) < 0)||(a(i,j) >= 0) && (a(i,j) < 22.5) || (a(i,j) <= -157.5) && (a(i,j) >= -180)||(a(i,j) >= 157.5)&&(a(i,j) <= 180))
            a(i,j) = 0;
        elseif((a(i,j) >= 22.5) && (a(i,j) < 67.5) || (a(i,j) <= -112.5) && (a(i,j) > -157.5))
            a(i,j) = -45;
        elseif((a(i,j) >= 67.5) && (a(i,j) < 112.5) || (a(i,j) <= -67.5) && (a(i,j) >- 112.5))
            a(i,j) = 90;
        elseif((a(i,j) >= 112.5) && (a(i,j) < 157.5) || (a(i,j) <= -22.5) && (a(i,j) > -67.5))
            a(i,j) = 45;  
        end
    end
end
%讨论在 3x3 区域内在该点梯度方向上进行非极大值抑制.获取非极大值抑制图像
Nms = zeros(H,W);%定义一个非极大值图像
for i = 2:H-1
    for j= 2:W-1
        if (a(i,j) == 0 && M(i,j) == max([M(i,j), M(i,j 1), M(i,j-1)]))
            Nms(i,j) = M(i,j);
        elseif (a(i,j) == -45 && M(i,j) == max([M(i,j), M(i 1,j-1), M(i-1,j 1)]))
            Nms(i,j) = M(i,j);
        elseif (a(i,j) == 90 && M(i,j) == max([M(i,j), M(i 1,j), M(i-1,j)]))
            Nms(i,j) = M(i,j);
        elseif (a(i,j) == 45 && M(i,j) == max([M(i,j), M(i 1,j 1), M(i-1,j-1)]))
            Nms(i,j) = M(i,j);
        end;
    end;
end;

%%  Step4:双阈值检测和连接边缘

DT = zeros(H,W);%定义一个双阈值图像
TL = 0.1 * max(max(Nms));%低阈值
TH = 0.3 * max(max(Nms));%高阈值
for i = 1  : H
    for j = 1 : W
        if (Nms(i, j) < TL)
            DT(i,j) = 0;
        elseif (Nms(i, j) > TH)
            DT(i,j) = 1 ;
        %对TL < Nms(i, j) < TH 使用8连通区域确定
        elseif ( Nms(i 1,j) > TH || Nms(i-1,j) >TH || Nms(i,j 1) > TH || Nms(i,j-1) > TH || Nms(i-1, j-1) > TH || Nms(i-1, j 1) > TH || Nms(i 1, j 1) > TH || Nms(i 1, j-1) > TH)
            DT(i,j) = 1;
        end;
    end;
end;
figure, imshow(DT); %最终的边缘检测为二值图像

效果图

发布者:全栈程序员栈长,转载请注明出处:https://javaforall.cn/139429.html原文链接:https://javaforall.cn

0 人点赞