一、Sobel算子
sobel算子的核心是像素矩阵的卷积,而卷积本质就是对指定的图像区域的像素值进行加权求和的过程,其计算过程为图像区域中的每个像素值分别与卷积模板的每个元素对应相乘,将卷积的结果作求和运算,运算到的和就是卷积运算的结果。
sobel算子包含垂直(左)和水平(右)两个方向的卷积模板
若A为原始图像,则
改变后的灰度值有两种计算方式:
最后设置一个阈值,运算后的像素值大于该阈值输出为0,小于该阈值输出为255。
sobel算子matlab代码:
%RGB_YCbCr
clc;
clear all;
close all;
RGB_data = imread('微信图片_20191227124400.jpg');%
R_data = RGB_data(:,:,1);
G_data = RGB_data(:,:,2);
B_data = RGB_data(:,:,3);
imshow(RGB_data);
title('原图像');
[ROW,COL, DIM] = size(RGB_data);
Y_data = zeros(ROW,COL);
Cb_data = zeros(ROW,COL);
Cr_data = zeros(ROW,COL);
Gray_data = RGB_data;
for r = 1:ROW
for c = 1:COL
Y_data(r, c) = 0.299*R_data(r, c) + 0.587*G_data(r, c) + 0.114*B_data(r, c);%彩色图像变为灰度图像公式Gray = R*0.299 + G*0.587 + B*0.114
Cb_data(r, c) = -0.172*R_data(r, c) - 0.339*G_data(r, c) + 0.511*B_data(r, c) + 128;
Cr_data(r, c) = 0.511*R_data(r, c) - 0.428*G_data(r, c) - 0.083*B_data(r, c) + 128;
end
end
Gray_data(:,:,1)=Y_data;
Gray_data(:,:,2)=Y_data;
Gray_data(:,:,3)=Y_data;
figure;
imshow(Gray_data);
title('灰度图像');
%Median Filter
imgn = imnoise(Gray_data,'salt & pepper',0.02); %当噪声类型是’salt & pepper’的时候,第三个参数的意思是噪声密度,比如0.1,那么总像素个数的10%为黑白点,当然是黑点还是白点都是随机的。
figure;
imshow(imgn);
title('椒盐噪声图像');
Median_Img = Gray_data;
for r = 2:ROW-1
for c = 2:COL-1
median3x3 =[imgn(r-1,c-1) imgn(r-1,c) imgn(r-1,c+1)
imgn(r,c-1) imgn(r,c) imgn(r,c+1)
imgn(r+1,c-1) imgn(r+1,c) imgn(r+1,c+1)];
sort1 = sort(median3x3, 2, 'descend');%对矩阵的每行进行降序排列
sort2 = sort([sort1(1), sort1(4), sort1(7)], 'descend');%对矩阵第一列进行降序排列
sort3 = sort([sort1(2), sort1(5), sort1(8)], 'descend');%对矩阵第二列进行降序排列
sort4 = sort([sort1(3), sort1(6), sort1(9)], 'descend');%对矩阵第三列进行降序排列
mid_num = sort([sort2(3), sort3(2), sort4(1)], 'descend');%对第一列最小值,第二列中间值,第三列最大值进行降序排列
Median_Img(r,c) = mid_num(2);%取中间值
end
end
figure;
imshow(Median_Img);
title('中值滤波后图像');
%Sobel_Edge_Detect
Median_Img = double(Median_Img);
Sobel_Threshold = 160;
Sobel_Img = zeros(ROW,COL);
for r = 2:ROW-1
for c = 2:COL-1
Sobel_x = Median_Img(r-1,c+1) + 2*Median_Img(r,c+1) + Median_Img(r+1,c+1) - Median_Img(r-1,c-1) - 2*Median_Img(r,c-1) - Median_Img(r+1,c-1);
Sobel_y = Median_Img(r-1,c-1) + 2*Median_Img(r-1,c) + Median_Img(r-1,c+1) - Median_Img(r+1,c-1) - 2*Median_Img(r+1,c) - Median_Img(r+1,c+1);
Sobel_Num = abs(Sobel_x) + abs(Sobel_y);
%Sobel_Num = sqrt(Sobel_x^2 + Sobel_y^2);
if(Sobel_Num > Sobel_Threshold)
Sobel_Img(r,c)=0;
else
Sobel_Img(r,c)=255;
end
end
end
figure;
imshow(Sobel_Img);
title('Sobel边缘检测图像');
Sobel算法只采用水平和垂直2个方向模板,对其他方向的边缘不敏感,对于纹理较复杂、斜向边缘较多的图像轮廓提取不是很理想,抗噪声能力也比较低。
改进方法:
在原有水平和垂直方向上增加了两个对角线方向上的模板
二、Canny 算子
Canny 算子进行边缘检测步骤:
1)对图像进行二维高斯滤波;
通过使用高斯滤波器与图像进行卷积,来减少噪声对边缘检测结果的影响和平滑图像,以减少边缘检测器上明显的噪声影响。
2)通过一阶微分计算图像的灰度梯度幅值和方向;
通过梯度方向确定边缘的方向,就可以把边缘的梯度方向大略分成几种角度(如0°、45°、90°和135 ), 并可以找到这个像素梯度方向的邻接像素。 该梯度值的计算和sobel计算过程一样,改变后的灰度值有两种计算方式:
3)对计算出的梯度幅值进行非极大值抑制(Non-MaximaSuppression, NMS);
由于梯度值最大不一定为边缘点,需要对一阶微分计算后的图像数据进行非极大值抑制,只保留幅值局部变化最大的点。将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。 如果当前像素的梯度强度与另外两个像素相比最大,则该像素点保留为边缘点,否则该像素点将被抑制。
4)通过人为设定的高低阈值确定图像的边缘。
为了解决噪声和颜色变化引起的一些边缘像素,必须用弱梯度值过滤边缘像素,并保留具有高梯度值的边缘像素,可以通过选择高低阈值来实现。如果边缘像素的梯度值高于高阈值,则将其标记为强边缘像素;如果边缘像素的梯度值小于高阈值并且大于低阈值,则将其标记为弱边缘像素;如果边缘像素的梯度值小于低阈值,则会被抑制。阈值的选择取决于给定输入图像的内容。
canny算子matlab代码:
clc;
clear all;
image = imread('微信图片_20191227124400.jpg');
image = rgb2gray(image);
subplot(221);
imshow(image);
title('原始图像');
image = double(image)/256;
[m,n] = size(image);
w = fspecial('gaussian');
image_1 = imfilter(image,w,'replicate');%滤波
subplot(222);
imshow(int8(256*image_1));
title('高斯滤波后的图像'); % 梯度计算
op = fspecial('sobel')/4; % 用Sobel算子来求导数
x = op'; y =op;
bx = imfilter(image_1,x,'replicate');
by = imfilter(image_1,y,'replicate');
b_abs = sqrt(bx.*bx+by.*by); % 求梯度的幅值
b_angle = angle(by-i*bx);
b_ang = b_angle/3.1416*180; % 求梯度的方向
% 梯度方向确定
for r = 1:m
for c = 1:n
if((b_ang(r,c)>=22.5 && b_ang(r,c)<67.5) || (b_ang(r,c)>=-157.5 && b_ang(r,c)<-112.5))
dir(r,c) = 1;
else if ((b_ang(r,c)>=67.5 && b_ang(r,c)<112.5)|| (b_ang(r,c)>=-112.5 && b_ang(r,c)<-67.5))
dir(r,c) = 2;
else if ((b_ang(r,c)>=112.5 && b_ang(r,c)<157.5)|| (b_ang(r,c)>=-67.5 && b_ang(r,c)<-22.5))
dir(r,c) = 3;
else
dir(r,c) = 0;
end
end
end
end
end
% 遍历图像
b_ab = [zeros(m,1),b_abs,zeros(m,1)]; % 串联矩阵
b_ab = [zeros(1,n+2);b_ab;zeros(1,n+2)]; %相当于在求梯度后的像素值加了一周0
for r = 2:m+1
for c = 2:n+1
switch dir(r-1,c-1)
case 0
if((b_ab(r,c)<b_ab(r+1,c))| (b_ab(r,c)<b_ab(r-1,c)))
b1(r-1,c-1) = 0;
else
b1(r-1,c-1) = b_ab(r,c);
end
case 1
if((b_ab(r,c)<b_ab(r+1,c-1))| (b_ab(r,c)<b_ab(r-1,c+1)))
b1(r-1,c-1) = 0;
else
b1(r-1,c-1) = b_ab(r,c);
end
case 2
if((b_ab(r,c)<b_ab(r,c-1))| (b_ab(r,c)<b_ab(r,c+1)))
b1(r-1,c-1) = 0;
else
b1(r-1,c-1) = b_ab(r,c);
end
case 3
if((b_ab(r,c)<b_ab(r-1,c-1))| (b_ab(r,c)<b_ab(r+1,c+1)))
b1(r-1,c-1) = 0;
else
b1(r-1,c-1) = b_ab(r,c);
end
end
end
end
for r = 1:m
for c = 1:n
if (b1(r,c)>0.24)
b2(r,c) = 1;
else
b2(r,c) = 0;
end
end
end
for r = 1:m
for c = 1:n
if(b1(r,c)>0.36)
b3(r,c)=1;
else
b3(r,c) = 0;
end
end
end
image_2 = b3;
for k = 1:10
for r = 2:m-1
for c = 2: n-1
if (b2(r,c)==1 && (image_2(r,c)==1||image_2(r+1,c)==1 ...
||image_2(r+1,c-1)==1||image_2(r+1,c+1)==1||image_2(r,c-1)==1 ...
||image_2(r,c+1)==1||image_2(r-1,c-1)==1||image_2(r-1,c)==1 ...
||image_2(r-1,c+1)==1))
image_2(r,c) = 1;
else
image_2(r,c) = 0;
end
end
end
end
subplot(223);
imshow(image_2);
title('Canny算子检测后的图像');
由于梯度计算对噪声非常敏感,并且采用高斯滤波需要人为设定方差,所以容易出现漏检现象。
改进方法:
利用双线性滤波代替高斯滤波,同时在非极大值抑制时将梯度幅值极大 值点进行保留并作为候选边缘点,并且进行自适应阈值选取。
clear all;
close all;
%%%%%%%%%%%%%%载入图像%%%%%%%%%%%%%%
X = zeros(200); % 定义一个灰度块矩阵A
X(1:100,1:100) = 85;
X(101:200,101:200) = 170;
X(101:200,1:100) = 255;
[dy, dx]=gradient(X,.1,.1); % 求梯度
figure('NumberTitle', 'off', 'Name', '原灰度块1');
imshow(uint8(X)) % 显示原灰度块1
[Ix,Iy]=size(X); %图像长宽
%%%%%%%%%%%%%%%%%%双边滤波器%%%%%%%%%%%%%%%%%%%%%%%
r=15; %模板半径
imgn=zeros(Ix+2*r+1,Iy+2*r+1);
imgn(r+1:Ix+r,r+1:Iy+r)=X;
imgn(1:r,r+1:Iy+r)=X(1:r,1:Iy); %扩展上边界
imgn(1:Ix+r,Iy+r+1:Iy+2*r+1)=imgn(1:Ix+r,Iy:Iy+r); %扩展右边界
imgn(Ix+r+1:Ix+2*r+1,r+1:Iy+2*r+1)=imgn(Ix:Ix+r,r+1:Iy+2*r+1); %扩展下边界
imgn(1:Ix+2*r+1,1:r)=imgn(1:Ix+2*r+1,r+1:2*r); %扩展左边界
sigma_d=7; sigma_r=0.2;
[x,y] = meshgrid(-r:r,-r:r);
w1=exp(-(x.^2+y.^2)/(2*sigma_d^2)); %以距离作为自变量高斯滤波器
for i=r+1:Ix+r
for j=r+1:Iy+r
w2=exp(-(imgn(i-r:i+r,j-r:j+r)-imgn(i,j)).^2/(2*sigma_r^2)); %以周围和当前像素灰度差值作为自变量的高斯滤波器
w=w1.*w2;
s=imgn(i-r:i+r,j-r:j+r).*w;
imgn(i,j)=sum(sum(s))/sum(sum(w));
end
end
figure;
imshow(mat2gray(imgn(r+1:Ix+r,r+1:Iy+r)));
title('双边滤波器图像');%显示滤波后的图像
%%%%%%%%%%%%%%二级尺度%%%%%%%%%%%%%%
m=1.0;
delta=2^m;
%%%%%%%%%高斯函数X,Y方向的偏导%%%%%%%
N=20; %长度和幅度
A=-1/sqrt(2*pi);
phi_x=zeros(N,N);
phi_y=zeros(N,N);
for i=1:N %偏导
for j=1:N
x=i-(N+1)/2;
y=j-(N+1)/2;
phi_x(i,j)=A*(x/delta^2).*exp(-(x*x+y*y)/(2*delta^2));
phi_y(i,j)=A*(y/delta^2).*exp(-(x*x+y*y)/(2*delta^2));
end
end
%%%%%%%%%图象行、列卷积%%%%%%%%%%%%%%
W_x=conv2(mat2gray(imgn(r+1:Ix+r,r+1:Iy+r)),phi_x,'same');
W_y=conv2(mat2gray(imgn(r+1:Ix+r,r+1:Iy+r)),phi_y,'same');
%%%%%%%%%%%%求梯度%%%%%%%%%%%%%%%%%%%
Mod_G=zeros(Ix,Iy);
for i2=1:Ix
for j2=1:Iy
Mod_G(i2,j2)=sqrt(W_x(i2,j2)*W_x(i2,j2)+W_y(i2,j2)*W_y(i2,j2));
end
end
M=fix(Mod_G);
%%%%%%%%%%%%求幅角%%%%%%%%%%%%%%%%%%%
Arg=zeros(Ix,Iy);
for i3=1:Ix
for j3=1:Iy
angle=atan(W_y(i3,j3)/W_x(i3,j3))*180/pi; %求反正切
if ((W_x(i3,j3)>0)&&(W_y(i3,j3)>0)) %第一象限
angle=angle+0;
elseif ((W_x(i3,j3)<0)&&(W_y(i3,j3)>0)) %第二象限
angle=angle+180;
elseif ((W_x(i3,j3)<0)&&(W_y(i3,j3)<0)) %第三象限
angle=angle+180;
elseif ((W_x(i3,j3)>0)&&(W_y(i3,j3)<0)) %第四象限
angle=angle+360;
elseif (abs(W_y(i3,j3)/W_x(i3,j3))>100000000) %大值90度
angle=90;
elseif (abs(W_y(i3,j3)/W_x(i3,j3))==0) %零值0度
angle=0;
end
if(angle>180) %幅角限制在0-180
angle=angle-180;
end
Arg(i3,j3)=angle;
end
end
%%%%%%%%%%%%图像边界%%%%%%%%%%%%%%%%%%%%
%%%%%%沿着梯度方向检测小波变换系数模的局部极大值点即可得到图像的边缘点%%%%%%%%%
Edge=zeros(Ix,Iy);
for i4=2:(Ix-1)
for j4=2:(Iy-1)
if ((Arg(i4,j4)>=0&&Arg(i4,j4)<=22.5))
if (Mod_G(i4,j4)>Mod_G(i4+1,j4)&&Mod_G(i4,j4)>Mod_G(i4-1,j4))%水平
Edge(i4,j4)=Mod_G(i4,j4);
end
elseif ((Arg(i4,j4)>=(90-22.5)&&Arg(i4,j4)<=(90+22.5)))
if (Mod_G(i4,j4)>Mod_G(i4,j4+1)&&Mod_G(i4,j4)>Mod_G(i4,j4-1))%45度
Edge(i4,j4)=Mod_G(i4,j4);
end
elseif ((Arg(i4,j4)>=(45-22.5)&&Arg(i4,j4)<=(45+22.5)))
if (Mod_G(i4,j4)>Mod_G(i4+1,j4+1)&&Mod_G(i4,j4)>Mod_G(i4-1,j4-1))%垂直
Edge(i4,j4)=Mod_G(i4,j4);
end
else
if (Mod_G(i4,j4)>Mod_G(i4+1,j4-1)&&Mod_G(i4,j4)>Mod_G(i4-1,j4+1))%135度
Edge(i4,j4)=Mod_G(i4,j4);
end
end
end
end
MAX_E=max(Edge(:));
if (MAX_E>0)
Edge=Edge/MAX_E; %归一化
end
%%%%%%%%%%%根据直方图获取高低阈值%%%%%%%%%%
count=0; %累计直方图统计可能是边缘点的个数
hist=zeros(1,1024);
for k1=1:Ix
for k2=1:Iy
if (Edge(k1,k2)~=0)
hist(1,M(k1,k2)+1)=hist(1,M(k1,k2)+1)+1;
count=count+1;
end
end
end
for k3=1:1024
if hist(1,k3)~=0
nmaxmag=k3;
end
end
p=0.7;
dot=ceil(p*count);
k4=1;
nedgenb=hist(1,1);
while (k4<(nmaxmag-1)&&(nedgenb<dot))
k4=k4+1;
nedgenb=nedgenb+hist(1,k4);
end
k4=k4/1000;
high=k4; %高阈值
low=0.4*high; %低阈值,一般高为低的2.5倍
%%%%%%%%%%%%边缘点判定%%%%%%%%%%%%%%%%%%
for i5=1:Ix
for j5=1:Iy
if (Edge(i5,j5)>=high) %高于高阈值一定是边缘
Edge(i5,j5)=255;
elseif (Edge(i5,j5)<=low) %低于高阈值一定是边缘
Edge(i5,j5)=0;
else %高低阈值之间判断8邻域是否高于阈值否则不是边缘
if((Edge(i5-1,j5)>high)&&(Edge(i5,j5-1)>high)&&(Edge(i5+1,j5)>=high)&&(Edge(i5,j5+1)>=high)...
&&(Edge(i5-1,j5-1)>high)&&(Edge(i5+1,j5-1)>high)&&(Edge(i5-1,j5+1)>high)&&(Edge(i5+1,j5+1)>high))
Edge(i5,j5)=255;
else
Edge(i5,j5)=0;
end
end
end
end
figure;
imshow(Edge);
title('边缘检测图像')
三、Roberts算子
Roberts算子和Sobel算子原理基本一样。只是Roberts算子使用的偶数模板,如2*2,通过对角线上的值对图像像素进行偏移,其垂直(左)和水平(右)卷积模板为:
若A为原始图像,则
梯度两种计算方式为:
通常为了简化计算,而采用第二种。
最后设置一个阈值,运算后的像素值大于该阈值输出为0,小于该阈值输出为255。
roberts算子matlab代码:
clc;
clear all;
img = imread('微信图片_20191227124400.jpg'); % 读取图像
figure;
imshow(img);
title('原始图像');
img_gray = rgb2gray(img); % 转换成灰度图
figure;
imshow(img_gray);
title('灰度图像');
[m,n] = size(img_gray); % 得到图像的大小
new_img_gray = img_gray; % 新建一个一样大的图像
pxValue = 0; % roberts计算所得到的像素值
% 对边界象素操作
threshold_value=30;
for i=1:m-1
for j=1:n-1
pxValue = abs(img_gray(i,j)-img_gray(i+1,j+1))+...
abs(img_gray(i+1,j)-img_gray(i,j+1));
if(pxValue > threshold_value)
new_img_gray(i,j) = 0;
else
new_img_gray(i,j) =255;
end
end
end
figure;
imshow(new_img_gray);
title('roberts边缘检测图像');
roberts算子对图像边缘检测而言定位精度高, 在水平和垂直方向效果较好。但该算子对噪声较 敏感,无法消除局部干扰,对图像中目标和背景灰度 差异表现并不显著的弱边缘却很难检测识别,这将 导致提取的目标边缘出现间断。并且该算子的阈值 需要人为设定,使得传统的 roberts算子在提取不 同目标轮廓时,具有很大的局限性。
改进方法:
采用3×3邻域代替roberts算法中2×2邻域来计算梯度幅值;并利用图像块之间相似性的 三维块匹配的去噪模型,提高roberts算子的检测精度和抗噪性能;通过最佳阈值迭代方法代替人为指定阈值 来获取最佳分割阈值,有效地提取图中目标轮廓。
四、Prewitt算子
Prewitt算子采用 3*3卷积 模板对区域内的像素值进行计算。
其垂直(左)和水平(右)卷积模板为:
若A为原始图像,则
梯度两种计算方式为:
最后设置一个阈值,运算后的像素值大于该阈值输出为0,小于该阈值输出为255。
prewitt算子matlab代码:
I=imread('微信图片_20191227124400.jpg');
I=im2double(I);
figure;
imshow(I);
title('原始图像');
img_gray = rgb2gray(I); % 转换成灰度图
figure;
imshow(img_gray);
title('灰度图像');
[height width R]=size(img_gray);
threshold_value=0.5;
for i=2:height-1
for j=2:width-1
Dx=[I(i+1,j-1)-I(i-1,j-1)]+[I(i+1,j)-I(i-1,j)]+[I(i+1,j+1)-I(i-1,j+1)];
Dy=[I(i-1,j+1)-I(i-1,j-1)]+[I(i,j+1)-I(i,j-1)]+[I(i+1,j+1)-I(i+1,j-1)];
P(i,j)=sqrt(Dx^2+Dy^2);
if (P(i,j)>threshold_value)
P(i,j)=0;
else
P(i,j)=255;
end
end
end
figure;
imshow(P);
title('prewitt边缘检测图像');
由于Prewitt 算子模板只能对水平和垂直方向的边缘进行检测,无法检测斜侧方向,使边缘检测不完整。
改进方法:
可以通过添加左右斜侧方向的模板,使边缘检测更加完整。
参考文献:[1]https://blog.csdn.net/ahilll/article/details/82050358
[2]https://blog.csdn.net/Jacky_Ponder/article/details/64126411?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.nonecase&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.nonecase
[3]https://blog.csdn.net/Prince_IT/article/details/90205259
[4]https://blog.csdn.net/weixin_43788499/article/details/84502478?locationNum=7&fps=1
[5]https://blog.csdn.net/baidu_21578557/article/details/51786249?locationNum=2&fps=1
[6]https://blog.51cto.com/3754839/1401776