MATLAB算子在图像边缘检测中的应用

MATLAB算子在图像边缘检测中的应用

1. 数学方法

边缘检测最通用的方法是检测灰度值的不连续性，这种不连续性用一阶和二阶导数来检测。

1. （1）一阶导数：一阶导数即为梯度，对于平面上的图像来说，我们只需用到二维函数的梯度，即： ,该向量的幅值：,为简化计算，省略上式平方根，得到近似值；或通过取绝对值来近似，得到：

（2）二阶导数：二阶导数通常用拉普拉斯算子来计算，由二阶微分构成：

1. 边缘检测的基本思想：
2. 寻找灰度的一阶导数的幅度大于某个指定阈值的位置；
3. 寻找灰度的二阶导数有零交叉的位置。
4. 几种方法简介
5. Sobel边缘检测器：以差分来代替一阶导数。Sobel边缘检测器使用一个3×3邻域的行和列之间的离散差来计算梯度，其中，每行或每列的中心像素用2来加权，以提供平滑效果。

-1	-2	1
0	0	0
1	2	1

-1	0	1
-2	0	2
-1	0	1

Prewitt边缘检测器:使用下图所示模板来数字化地近似一阶导数。与Sobel检测器相比，计算上简单一些，但产生的结果中噪声可能会稍微大一些。

-1	-1	-1
0	0	0
1	1	1

-1	0	1
-1	0	1
-1	0	1

Roberts边缘检测器：使用下图所示模板来数字化地将一阶导数近似为相邻像素之间的差，它与前述检测器相比功能有限（非对称，且不能检测多种45°倍数的边缘）。

-1	0
0	1

0	-1
1	0

1、Laplace边缘检测器：二维函数的拉普拉斯是一个二阶的微分定义：

模板算子可分为四邻域和八邻域，如下：

0	1	0
1	-4	1
0	1	0

（四邻域）

1	1	1
1	-8	1
1	1	1

（八邻域）

1. LoG边缘检测器

由于噪声点（灰度与周围点相差很大的像素点）对边缘检测有一定的影响，所以效果更好的是LoG算子，即Laplacian-Guass算子。引入高斯函数来平滑噪声：

该函数的Laplace算子：

它把Guass平滑滤波器和Laplace锐化滤波器结合起来，先平滑掉噪声，再进行边缘检测，所以效果比单用Laplace算子要更为平滑，效果更好。

1. Canny边缘检测器

主要分为以下几个步骤：①使用具有指定标准差的一个高斯滤波器来平滑图像，以减少噪声；②在每个点处计算局部梯度和边缘方向；③对步骤②中确定的边缘点产生梯度中的脊线顶部进行追踪，并将实际山不在脊线顶部的像素设置为零，从而在输出中给出一条细线（非最大值抑制），然后使用滞后阈值处理法对这些脊线像素进行阈值处理。最后进行边缘连接。

1. 实验结果

原始图像：

1.Roberts算子

2.Prewitt算子

3.Sobel算子

4.Laplace算子

（4邻域）

（8邻域）

1、LoG

四种结果的比较：（Laplace采用8邻域模板）

三、讨论

1. 各方法优缺点及比较：
2. Roberts算子：Roberts算子是2×2算子，对具有陡峭的低噪声图像响应最好，并且检测垂直边缘的效果好于斜向边缘，定位精度高。然而，它对噪声敏感，无法抑制噪声的影响。因此，它适用于边缘明显且噪声较少的图像分割。
3. Prewitt算子：Prewitt算子将两个点的各自一定领域内的灰度值求和，并根据两个灰度值和的差来计算x,y的偏导数。它是平均滤波，对噪声有抑制作用，对于灰度渐变的低噪声图像有较好的检测效果，但是像素平均相当于对图像的低通滤波，所以它对边缘的定位不如Roberts算子。对于混合多复杂噪声的图像，效果不太理想。
4. Sobel算子：Sobel算子是滤波算子的形式，用于提取边缘，可以利用快速卷积函数， 简单有效，因此应用广泛。美中不足的是，Sobel算子并没有将图像的主体与背景严格地区分开来，即Sobel算子没有严格地模拟人的视觉生理特征，所以提取的图像轮廓有时并不能令人满意。
5. Laplace算子：Laplace算子是一种各向同性算子，在只关心边缘的位置而不考虑其周围的象素灰度差值时比较合适。Laplace算子对孤立象素的响应要比对边缘或线的响应要更强烈，因此只适用于无噪声图象。存在噪声情况下，使用Laplacian算子检测边缘之前需要先进行低通滤波。
6. LoG算子：由于二阶导数是线性运算，利用LoG卷积一幅图像与首先使用高斯型平滑函数卷积改图像，然后计算所得结果的拉普拉斯是一样的。所以在LoG公式中使用高斯函数的目的就是对图像进行平滑处理，使用Laplacian算子的目的是提供一幅用零交叉确定边缘位置的图像；图像的平滑处理减少了噪声的影响并且它的主要作用还是抵消由Laplacian算子的二阶导数引起的逐渐增加的噪声影响。
7. Canny算子：是一个具有滤波，增强和检测的多阶段的优化算子，在进行处理前，Canny算子先利用高斯平滑滤波器来平滑图像以消除噪声，Canny分割算法采用一阶偏导的有限差分来计算梯度的幅值和方向。在处理过程中，Canny算法还将经过一个非极大值抑制的过程。最后Canny算法将采用两个阈值来连接边缘。高定位精度、低误判率、抑制虚假边缘，适用于高噪声图像。

四、实际应用

1、概述：图像边缘检测广泛应用于车牌识别，人脸识别，地震带检测，生物医疗以及产品外观检测等方面。为了体现其实际应用，在此针对车牌识别来进行实验，实验流程如下：

确定区域

文字分割

模板配对

均值滤波

边缘检测

二值化

灰度处理

1. 实验结果

原图像：

灰度处理：将彩色图像转化为灰度图像，缩小图片占用的空间，减少处理时间

边缘检测：车牌区域的边缘信息一般比其他部分丰富且集中，只要 背景不太复杂，没有过多干扰，边缘特征可将车牌区域与其他部分区别开来，从而消除图像中的无关信息。

腐蚀和膨胀：腐蚀可以分割独立的图像元素，膨胀用于连接相邻的元素，更加利于图像的分割

最后将车牌号分割成七个单个字符，建立模板库，与其进行配对，最终得出结果。

五、程序介绍

1. Roberts边缘检测

 clear;
sourcePic=imread('C:\Users\34899\Desktop\timg.jpg'); %读取原图像
grayPic=mat2gray(sourcePic); %实现图像矩阵的归一化操作
[m,n]=size(grayPic);
newGrayPic=grayPic;%为保留图像的边缘一个像素
robertsNum=0; %经roberts算子计算得到的每个像素的值
robertThreshold=0.2; %设定阈值
for j=1:m-1 %进行边界提取
    for k=1:n-1
        robertsNum = abs(grayPic(j,k)-grayPic(j+1,k+1)) + abs(grayPic(j+1,k)-grayPic(j,k+1));
        if(robertsNum > robertThreshold)
            newGrayPic(j,k)=255;
        else
            newGrayPic(j,k)=0;
        end
    end
end
subplot(1,2,1);imshow(sourcePic);title('原图'); 
subplot(1,2,2);imshow(newGrayPic);title('Robert算子处理后图像'); 

1. Prewitt边缘检测

 clear;
sourcePic=imread('C:\Users\34899\Desktop\timg.jpg'); %读取原图像
grayPic=mat2gray(sourcePic);%实现图像矩阵的归一化操作
[m,n]=size(grayPic);
newGrayPic=grayPic;%为保留图像的边缘一个像素
PrewittNum=0;%经Prewitt算子计算得到的每个像素的值
PrewittThreshold=0.5;%设定阈值
for j=2:m-1 %进行边界提取
		for k=2:n-1    					PrewittNum=abs(grayPic(j-1,k+1)-grayPic(j+1,k+1)+grayPic(j-1,k	)-grayPic(j+1,k)+grayPic(j-1,k-1)-grayPic(j+1,k-1))+abs(grayPic(j-1,k+1)+grayPic(j,k+1)+grayPic(j+1,k+1)-grayPic(j-1,k-1)-grayPic(j,k-1)-grayPic(j+1,k-1));
        if(PrewittNum > PrewittThreshold)
            newGrayPic(j,k)=255;
        else
            newGrayPic(j,k)=0;
        end
    end
end
subplot(1,2,1);imshow(sourcePic);title('原图'); 
subplot(1,2,2);imshow(newGrayPic);title('Prewitt算子处理后图像'); 

1. Sobel边缘检测

      clear all;
	  close all;
	  imag = imread('C:\Users\34899\Desktop\timg.jpg');  %读取关键帧
      imag = rgb2gray(imag);        %转化为灰度图
      subplot(1,2,1);imshow(imag);title('原图'); 
      [high,width] = size(imag);   % 获得图像的高度和宽度
      F2 = double(imag);        
      U = double(imag);       
      uSobel = imag;
      for i = 2:high - 1   %sobel边缘检测
           for j = 2:width - 1
                Gx = (U(i+1,j-1) + 2*U(i+1,j) + F2(i+1,j+1)) - (U(i-1,j-1) + 2*U(i-1,j) + F2(i-1,j+1));
                Gy = (U(i-1,j+1) + 2*U(i,j+1) + F2(i+1,j+1)) - (U(i-1,j-1) + 2*U(i,j-1) + F2(i+1,j-1));
                uSobel(i,j) = sqrt(Gx^2 + Gy^2); 
           end
      end 
      subplot(1,2,2);imshow(im2uint8(uSobel)):title('Sobel边缘检测后'); 

1. Laplace边缘检测

主函数：

f=imread('C:\Users\34899\Desktop\timg.jpg')
fGray=rgb2gray(f);
figure()
imshow(fGray),title('灰度图像');
fGray=double(fGray);
T=60;
LapModType=8;%设置laplace模板方式
fLapEdge=LaplaceEdge(fGray,LapModType,T);
fGrayLapEdge=uint8(fLapEdge);
figure()
imshow(fLapEdge),title('laplace边缘图像');

四邻域算子和八邻域算子

function fLapEdge=LaplaceEdge(fGray,MoldType,Thresh)
%-------------------------参数介绍-----------------
%输入参数：
%         fGray;输入的灰度图像
%         MoldType:模板类型，包括四邻域和八邻域
%         Thresh:边缘检测阈值
%输出参数：
%         fEdge:边缘像素点，存储的是二值化图像
[r,c]=size(fGray);
fLapEdge=zeros(r,c);
 
%四邻域拉普拉斯边缘检测算子
if 4==MoldType
    for i=2:r-1
        for j=2:c-1
            Temp=-4*fGray(i,j)+fGray(i-1,j)+fGray(i+1,j)+fGray(i,j-1)+fGray(i,j+1);
            if Temp>Thresh
                fLapEdge(i,j)=255;
            else
                fLapEdge(i,j)=0;
            end
        end
    end
end
%八邻域拉普拉斯边缘检测算子
if 8==MoldType
    for i=2:r-1
        for j=2:c-1
            Temp=-8*fGray(i,j)+fGray(i-1,j)+fGray(i+1,j)+fGray(i,j-1)+fGray(i,j+1)+fGray(i-1,j-1)+fGray(i+1,j+1)+fGray(i+1,j-1)+fGray(i-1,j+1);
            if Temp>Thresh
                fLapEdge(i,j)=255;
            else
                fLapEdge(i,j)=0;
            end
        end
    end
End

1. 拉普拉斯—高斯检测器

clc; clear all; 
close all;
hsize = [5 5];
sigma = 0.5;
h = fspecial('log', hsize, sigma);
I = imread('C:\Users\34899\Desktop\timg.jpg');
bw = imfilter(I, h, 'replicate');
figure;
subplot(1, 2, 1); imshow(I, []); title('原图像', 'FontWeight', 'Bold');
subplot(1, 2, 2); imshow(bw, []); title('边缘图像', 'FontWeight', 'Bold');

1. Canny算子

tic  % 记录CPU使用时间，开始计时
I = imread('C:\Users\34899\Desktop\timg.jpg');    % 读入图像
I=rgb2gray(I);               % 转化为灰色图像
% figure,imshow(I);title('原图')        % 显示原图像
% a = im2single(I); % 将图像矩阵转换为单精度类型
% a = im2uint8(a);  % 将图像矩阵转换为8位无符号整数类型
% I = uint16(I);  % 8位无符号整数转化为16位无符号整数
 
%**************************** matlab的高斯平滑 ****************************
[row,line] = size(I);        % 行列
temp = zeros(row+8,line+8);  % 申请空间
imagedata = zeros(row,line); % 申请空间
% 高斯核计算,核大小为9*9
GaussianDieOff = .0001;
sigma = 1;
pw = 1:30; % possible widths
ssq = sigma^2;
width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');
if isempty(width)
  width = 1;  % the user entered a really small sigma
end
t = (-width:width);
gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq); 
kernel = gau' * gau;  % 高斯核,大小为9*9
kernel = im2single(kernel);
 
% 图像边缘复制
for q = 5:row+4  % 首先载入源图像信息
    for p = 5:line+4
        temp(q,p) = I(q-4,p-4);
    end
end
for q = 1:row+8  % 复制列5:line+4
    for p = 5:line+4
        if(q<5) 
            temp(q,p) = temp(5,p);
        end
        if(q>row+4) 
            temp(q,p) = temp(row+4,p); 
        end
    end
end
for q = 1:row+8  % 复制行1:row+8
    for p = 1:line+8
        if(p<5) 
            temp(q,p) = temp(q,5); 
        end
        if(p>line+4) 
            temp(q,p) = temp(q,line+4); 
        end
    end
end
for q = 1:row  % 高斯卷积
    for p = 1:line
        for n = 1:9
            for m = 1:9
                imagedata(q,p) = imagedata(q,p) + kernel(n,m) * temp(q+n-1,p+m-1);
            end
        end
    end
end
%************************ matlab的高斯平滑说明结束 ************************
 
%****************************** matlab求梯度 ******************************
dx = zeros(size(I));         % 申请空间,存放x方向的幅值
dy = zeros(size(I));         % 申请空间,存放y方向的幅值
[x,y]=meshgrid(-width:width,-width:width);
dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq); % 核
% 图像边缘复制
for q = 5:row+4  % 首先载入源图像信息
    for p = 5:line+4
        temp(q,p) = imagedata(q-4,p-4);
    end
end
for q = 1:row+8  % 复制列5:line+4
    for p = 5:line+4
        if(q<5) 
            temp(q,p) = temp(5,p); 
        end
        if(q>row+4) 
            temp(q,p) = temp(row+4,p); 
        end
    end
end
for q = 1:row+8  % 复制行1:row+8
    for p = 1:line+8
        if(p<5) 
            temp(q,p) = temp(q,5); 
        end
        if(p>line+4) 
            temp(q,p) = temp(q,line+4); 
        end
    end
end
for q = 1:row  % x方向卷积
    for p = 1:line
        for n = 1:9
            for m = 1:9
                dx(q,p) = dx(q,p) + dgau2D(n,m) * temp(q+n-1,p+m-1);
            end
        end
    end
end
dgau2D = dgau2D';
for q = 1:row  % y方向卷积
    for p = 1:line
        for n = 1:9
            for m = 1:9
                dy(q,p) = dy(q,p) + dgau2D(n,m) * temp(q+n-1,p+m-1);
            end
        end
    end
end
%**************************** matlab求梯度结束 ****************************
 
%***************************** 自适应阈值计算 *****************************
value = sqrt(dx.^2 + dy.^2);  % 幅值
value = double(value);
value = value/max(max(value));% 归一化 
temp_value = value;
counts=imhist(value, 64);
high_threshold = find(cumsum(counts) > 0.7*row*line,1,'first') / 64;
low_threshold = 0.4 * high_threshold;  % 低阈值
clear counts
%*************************** 自适应阈值计算结束 ***************************
 
%***************************** 方向与阈值判别 *****************************
num = 0;      % 当前堆栈个数
flag = zeros(80000,2);       % 堆栈
temp_flag = zeros(80000,2);  % 临时堆栈
imagedata = zeros(row,line); % 初始化
% 方向：
% 0为0°~45°
% 1为45°~90°
% 2为90°~135°
% 3为135°~180°
direction = zeros(size(I));  % 申请空间,存放小于高阈值而大于低阈值
for q = 2:row-1
    for p = 2:line-1
        if((dy(q,p)<=0 && dx(q,p)>-dy(q,p))  || (dy(q,p)>=0 && dx(q,p)<-dy(q,p)))
            d = abs(dy(q,p)/dx(q,p));  % 0°~45°方向
            gradmag = temp_value(q,p);
            gradmag1 = temp_value(q,p+1)*(1-d) + temp_value(q-1,p+1)*d; 
            gradmag2 = temp_value(q,p-1)*(1-d) + temp_value(q+1,p-1)*d;
        elseif((dx(q,p)>0 && -dy(q,p)>=dx(q,p))  || (dx(q,p)<0 && -dy(q,p)<=dx(q,p)))
            d = abs(dx(q,p)/dy(q,p));  % 45°~90°方向
            gradmag = temp_value(q,p);
            gradmag1 = temp_value(q-1,p)*(1-d) + temp_value(q-1,p+1)*d; 
            gradmag2 = temp_value(q+1,p)*(1-d) + temp_value(q+1,p-1)*d;
        elseif((dx(q,p)<=0 && dx(q,p)>dy(q,p)) || (dx(q,p)>=0 && dx(q,p)<dy(q,p)))
            d = abs(dx(q,p)/dy(q,p));  % 90°~135°方向
            gradmag = temp_value(q,p);
            gradmag1 = temp_value(q-1,p)*(1-d) + temp_value(q-1,p-1)*d; 
            gradmag2 = temp_value(q+1,p)*(1-d) + temp_value(q+1,p+1)*d;
        elseif((dy(q,p)<0 && dx(q,p)<=dy(q,p)) || (dy(q,p)>0 && dx(q,p)>=dy(q,p)))
            d = abs(dy(q,p)/dx(q,p));  % 135°~180°方向
            gradmag = temp_value(q,p);
            gradmag1 = temp_value(q,p-1)*(1-d) + temp_value(q-1,p-1)*d; 
            gradmag2 = temp_value(q,p+1)*(1-d) + temp_value(q+1,p+1)*d;
        end
        if(gradmag>=gradmag1 && gradmag>=gradmag2)
            if(gradmag >= high_threshold)   % 高阈值判别
                value(q,p) = 255;
            elseif(gradmag >= low_threshold)% 低阈值判别
                value(q,p) = 125;
            else
                value(q,p) = 0;  % 小于低阈值
            end
        else
            value(q,p) = 0;  % 非极大值抑制
        end
    end 
end
for q = 2:row-1      % 高阈值判别，检查高阈值邻域8个方向范围内是否存在低阈值
    for p = 2:line-1
        if(value(q,p) == 255)
            imagedata(q,p) = 255;
            if(value(q-1,p-1) == 125)
                value(q-1,p-1) = 255;
                imagedata(q-1,p-1) = 255;
                if((q-1 > 1) && (p-1 > 1))
                    num = num + 1;
                    flag(num,1) = q-1;
                    flag(num,2) = p-1;
                end
            end
            if(value(q-1,p) == 125)
                value(q-1,p) = 255;
                imagedata(q-1,p) = 255;
                if(q-1 > 1)
                    num = num + 1;
                    flag(num,1) = q-1;
                    flag(num,2) = p;
                end
            end
            if(value(q-1,p+1) == 125)
                value(q-1,p+1) = 255;
                imagedata(q-1,p+1) = 255;
                if((q-1 > 1) && (p+1 < line))
                    num = num + 1;
                    flag(num,1) = q-1;
                    flag(num,2) = p+1;
                end
            end
            if(value(q,p-1) == 125)
                value(q,p-1) = 255;
                imagedata(q,p-1) = 255;
                if(p-1 > 1)
                    num = num + 1;
                    flag(num,1) = q;
                    flag(num,2) = p-1;
                end
            end
            if(value(q,p+1) == 125)
                value(q,p+1) = 255;
                imagedata(q,p+1) = 255;
                if(p+1 < line)
                    num = num + 1;
                    flag(num,1) = q;
                    flag(num,2) = p+1;
                end
            end
            if(value(q+1,p-1) == 125)
                value(q+1,p-1) = 255;
                imagedata(q+1,p-1) = 255;
                if((q+1 < row) && (p-1 > 1))
                    num = num + 1;
                    flag(num,1) = q+1;
                    flag(num,2) = p-1;
                end
            end
            if(value(q+1,p) == 125)
                value(q+1,p) = 255;
                imagedata(q+1,p) = 255;
                if(q+1 < row)
                    num = num + 1;
                    flag(num,1) = q+1;
                    flag(num,2) = p;
                end
            end
            if(value(q+1,p+1) == 125)
                value(q+1,p+1) = 255;
                imagedata(q+1,p+1) = 255;
                if((q+1 < row) && (p+1 < line))
                    num = num + 1;
                    flag(num,1) = q+1;
                    flag(num,2) = p+1;
                end
            end
        end
    end
end
done = num;  % 完成标志,等于0表示当前连线已完成
while done ~= 0
    num = 0;
    for temp_num = 1:done
        q = flag(temp_num,1);
        p = flag(temp_num,2);
        if(value(q-1,p-1) == 125)
            value(q-1,p-1) = 255;
            imagedata(q-1,p-1) = 255;
            if((q-1 > 1) && (p-1 > 1))
                num = num + 1;
                temp_flag(num,1) = q-1;
                temp_flag(num,2) = p-1;
            end
        end
        if(value(q-1,p) == 125)
            value(q-1,p) = 255;
            imagedata(q-1,p) = 255;
            if(q-1 > 1)
                num = num + 1;
                temp_flag(num,1) = q-1;
                temp_flag(num,2) = p;
            end
        end
        if(value(q-1,p+1) == 125)
            value(q-1,p+1) = 255;
            imagedata(q-1,p+1) = 255;
            if((q-1 > 1) && (p+1 < line))
                num = num + 1;
                temp_flag(num,1) = q-1;
                temp_flag(num,2) = p+1;
            end
        end
        if(value(q,p-1) == 125)
            value(q,p-1) = 255;
            imagedata(q,p-1) = 255;
            if(p-1 > 1)
                num = num + 1;
                temp_flag(num,1) = q;
                temp_flag(num,2) = p-1;
            end
        end
        if(value(q,p+1) == 125)
            value(q,p+1) = 255;
            imagedata(q,p+1) = 255;
            if(p+1 < line)
                num = num + 1;
                temp_flag(num,1) = q;
                temp_flag(num,2) = p+1;
            end
        end
        if(value(q+1,p-1) == 125)
            value(q+1,p-1) = 255;
            imagedata(q+1,p-1) = 255;
            if((q+1 < row) && (p-1 > 1))
                num = num + 1;
                temp_flag(num,1) = q+1;
                temp_flag(num,2) = p-1;
            end
        end
        if(value(q+1,p) == 125)
            value(q+1,p) = 255;
            imagedata(q+1,p) = 255;
            if(q+1 < row)
                num = num + 1;
                temp_flag(num,1) = q+1;
                temp_flag(num,2) = p;
            end
        end
        if(value(q+1,p+1) == 125)
            value(q+1,p+1) = 255;
            imagedata(q+1,p+1) = 255;
            if((q+1 < row) && (p+1 < line))
                num = num + 1;
                temp_flag(num,1) = q+1;
                temp_flag(num,2) = p+1;
            end
        end
    end
    done = num;
    flag = temp_flag;
end
%*************************** 方向与阈值判别结束 ***************************
 
%****************************** 目标图像显示 ******************************
imagedata = uint8(imagedata);  % 16位无符号整数转化为8位无符号整数
figure,imshow(imagedata);title('canny算子检测后')     % 显示非极大值抑制后的图像
clear all
%**************************** 目标图像显示结束 ****************************
 
toc  % 记录CPU使用时间，结束计时