边缘提取与锐化增强
基于单方向梯度算子,Robert算子,Sobel算子,Pwitter算子,Laplacian算子,多方向模板,LoG算子编写的matlab增强程序。
原理
1.单方向梯度算子
用差分定义一元函数f(x)一阶微分:
水平方向梯度计算方法为:
垂直方向梯度计算方法为:
若想得到常规梯度模板,用绝对值近似平方和平方根操作以方便计算:
2.Robert算子
利用以下图中的符号表示一个3*3区域内图像点的灰度:
Robert交叉差分算子中:
则梯度图像计算为:
为了方便运算,可利用以下公式:
故Robert交叉梯度算子计算方法为
3.Sobel算子
利用图(a)中的像素点,当以z5为对称中心的3*3领域内时,我们可得:
即:
故Sobel算子的计算方法为:
同样,为了在实际问题中方便运算,我们用绝对值近似平方和平方根:
4.Priwitt算子
Priwitt算子的计算方法为:
同样,为了在实际问题中方便运算,我们用绝对值近似平方和平方根:
5.Laplacian算子
我们此处用的是Laplacian算子四邻域模板:
也可用Laplacian算子的八邻域模板:
故四邻域Laplacian算子计算方法为:
八邻域Laplacian算子计算方法为:
6.多方向模板
采用Kirsch算子
Krisch 算子8个方向的模板为:
取这8个方向最大值,即产生最大输出值的模板方向,表示该点边缘的方向。如果所有方向上的边缘模板接近于0,该点没有边缘;如果所有方向的边缘模板输出值都近似相等,则没有可靠边缘。
7.LoG算子
G(x,y)是对图像处理的高斯函数:
对平滑后的图像做拉普拉斯变换,得:
即先对图像平滑。后对图像做拉普拉斯变换求二阶微分,再与原来的图像做卷积,此处我所用模板为5*5LoG模板:
边缘提取结果
LOG算子放大图如下:
锐化增强结果
对比分析
单方向模板具有单一性,对单方向边缘的提取效果较好,即水平、垂直方向;Robert算子利用局部差分寻找边缘,边缘定位较为准确,尤其是45度方向的边缘,但容易丢失一部分边缘,不具备抑制噪声的能力;Sobel算子和Prewitt算子对图像的处理更好,但容易出现多像素宽度,两者的区别在于权重大小;Laplacian 算子对噪声敏感,使噪声加强,也容易丢失部分边缘信息;LoG能够检测出较细边缘,且能有较好的去噪效果,边缘提取效果较为理想。LoG算子放大图如下所示:
代码
锐化增强
%Write by 长安 Rjex
function enhance(name)Iold = imread(name);
I = double(Iold);
[x,y,z] = size(I);
A = [-1 1];
B = [1
-1];
I3 = zeros(x,y,z);
for k = 1:z
for i = 1:x
for j = 1:y-1
p = I(i,j:j+1,k);
I3(i,j,k) = I(i,j,k) + abs(sum(sum(p.*A)));
end
end
end
I3 = uint8(I3);
subplot(2,4,1);imshow(I3),title('水平梯度算子');
for k = 1:z
for i = 1:x-1
for j = 1:y
p = I(i:i+1,j,k);
I3(i,j,k) =I(i,j,k)+abs(sum(sum(p.*B)));
end
end
end
I3 = uint8(I3);
subplot(2,4,2);imshow(I3),title('垂直梯度算子');
A = [-1 0
0 1];
B = [0 -1
1 0];I3 = zeros(x,y,z);
for k = 1:z
for i = 1:x-1
for j = 1:y-1
p = I(i:i+1,j:j+1,k);
I3(i,j,k) = I(i,j,k)+abs(sum(sum(p.*A)))+abs(sum(sum(p.*B)));
end
end
end
I3 = uint8(I3);
subplot(2,4,3);
imshow(I3),title('Robert算子');
I3 = double(I3);
C = [-1 -2 -1
0 0 0
1 2 1];
D = [-1 0 1
-2 0 2
-1 0 1];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
p = I(i-1:i+1,j-1:j+1,k);
I3(i,j,k) = I(i,j,k)+max(abs(sum(sum(p.*C))),abs(sum(sum(p.*D))));
end
end
end
I3 = uint8(I3);
subplot(2,4,4);imshow(I3),title('Sobel算子');
I3 = double(I3);
C = [-1 -1 -1
0 0 0
1 1 1];
D = [-1 0 1
-1 0 1
-1 0 1];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
p = I(i-1:i+1,j-1:j+1,k);
I3(i,j,k) = I(i,j,k)+max(abs(sum(sum(p.*C))),abs(sum(sum(p.*D))));
end
end
end
I3 = uint8(I3);
subplot(2,4,5);imshow(I3),title('Priwitt算子');
C = [0 -1 0
-1 5 -1
0 -1 0];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
p = I(i-1:i+1,j-1:j+1,k);
I3(i,j,k) = sum(sum(p.*C));
end
end
end
I3 = uint8(I3);
subplot(2,4,6);imshow(I3),title('Lap算子');
E(:,:,1) = [5 5 5
-3 0 -3
-3 -3 -3];
E(:,:,2) = [-3 5 5
-3 0 5
-3 -3 -3];
E(:,:,3) = [-3 -3 5
-3 0 5
-3 -3 5];
E(:,:,4) = [-3 -3 -3
-3 0 5
-3 5 5];
(:,:,5) = [-3 -3 -3
-3 0 -3
5 5 5];
E(:,:,6) = [-3 -3 -3
5 0 -3
5 5 -3];
E(:,:,7) = [5 -3 -3
5 0 -3
5 -3 -3];
E(:,:,8) = [5 5 -3
5 0 -3
-3 -3 -3];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
for l =1:8
p = I(i-1:i+1,j-1:j+1,k);
flag(l) = sum(sum(p.*E(:,:,l)));
end
I3(i,j,k) = I(i,j,k)+max(flag);
end
end
end
I3 = uint8(I3);
subplot(2,4,7);imshow(I3),title('多方向模板');
C = [-2 -4 -4 -4 -2
-4 0 8 0 -4
-4 8 24 8 -4
-4 0 8 0 -4
-2 -4 -4 -4 -2];
for k = 1:z
for i = 3:x-2
for j = 3:y-2
p = I(i-2:i+2,j-2:j+2,k);
I3(i,j,k) = I(i,j,k)+sum(sum(p.*C));
end
end
end
I3 = uint8(I3);
subplot(2,4,8);imshow(I3),title('LoG算子');
锐化增强
%Write by 长安 Rjex
3.edgeImage.m
function edgeImage(name)
I = double(imread(name));
[x,y,z] = size(I);
A = [-1 1];B = [1 -1];
I3 = zeros(x,y,z);
for k = 1:z
for i = 1:x
for j = 1:y-1
p = I(i,j:j+1,k);
I3(i,j,k) = abs(sum(p.*A));
end
end
end
I3 = uint8(I3);
subplot(2,4,1);imshow(I3),title('水平梯度算子');
for k = 1:z
for i = 1:x-1
for j = 1:y
p = I(i:i+1,j,k);
I3(i,j,k) = abs(sum(p.*B));
end
end
end
I3 = uint8(I3);
subplot(2,4,2);imshow(I3),title('垂直梯度算子');
A = [-1 0
0 1];
B = [0 -1
1 0];
I3 = zeros(x,y,z);
for k = 1:z
for i = 1:x-1
for j = 1:y-1
p = I(i:i+1,j:j+1,k);
I3(i,j,k) = abs(sum(sum(p.*A)))+abs(sum(sum(p.*B)));
end
end
end
I3 = uint8(I3);
subplot(2,4,3);imshow(I3),title('Robert算子');
I3 = double(I3);
C = [-1 -2 -1
0 0 0
1 2 1];
D = [-1 0 1
-2 0 2
-1 0 1];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
p = I(i-1:i+1,j-1:j+1,k);
I3(i,j,k) = abs(sum(sum(p.*C)))+abs(sum(sum(p.*D)));
end
end
end
I3 = uint8(I3);
subplot(2,4,4);imshow(I3),title('Sobel算子');
I3 = double(I3);
C = [-1 -1 -1
0 0 0
1 1 1];
D = [-1 0 1
-1 0 1
-1 0 1];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
p = I(i-1:i+1,j-1:j+1,k);
I3(i,j,k) = abs(sum(sum(p.*C)))+abs(sum(sum(p.*D)));
end
end
end
I3 = uint8(I3);
subplot(2,4,5);
imshow(I3),title('Priwitt算子');
C = [0 -1 0
-1 4 -1
0 -1 0];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
p = I(i-1:i+1,j-1:j+1,k);
I3(i,j,k) = sum(sum(p.*C));
end
end
end
I3 = uint8(I3);
subplot(2,4,6);imshow(I3),title('Lap算子');
E(:,:,1) = [5 5 5
-3 0 -3
-3 -3 -3];
E(:,:,2) = [-3 5 5
-3 0 5
-3 -3 -3];
E(:,:,3) = [-3 -3 5
-3 0 5
-3 -3 5];
E(:,:,4) = [-3 -3 -3
-3 0 5
-3 5 5];
E(:,:,5) = [-3 -3 -3
-3 0 -3
5 5 5];
E(:,:,6) = [-3 -3 -3
5 0 -3
5 5 -3];
E(:,:,7) = [5 -3 -3
5 0 -3
5 -3 -3];
E(:,:,8) = [5 5 -3
5 0 -3
-3 -3 -3];
for k = 1:z
for i = 2:x-1
for j = 2:y-1
for l =1:8
p = I(i-1:i+1,j-1:j+1,k);
flag(l) = sum(sum(p.*E(:,:,l)));
end
I3(i,j,k) = max(flag);
end
end
end
I3 = uint8(I3);
subplot(2,4,7);imshow(I3),title('多方向模板');
C = [-2 -4 -4 -4 -2
-4 0 8 0 -4
-4 8 24 8 -4
-4 0 8 0 -4
-2 -4 -4 -4 -2];
for k = 1:z
for i = 3:x-2
for j = 3:y-2
p = I(i-2:i+2,j-2:j+2,k);
I3(i,j,k) = sum(sum(p.*C));
end
end
end
I3 = uint8(I3);
subplot(2,4,8);imshow(I3),title('LoG算子');