实验二 图像直方图及灰度变换
一、实验目的与要求
- 掌握图像灰度直方图的概念及其计算方法，编写灰度直方图统计程序。
- 通过对图像直方图的分析，学习应用直方图法解决诸如图像二值化等具体问题。
- 熟悉直方图均衡化的计算过程及其应用。
- 掌握图像灰度变换技术，通过调整图像的对比度和亮度等参数，改善视觉效果。
二、实验相关知识
直方图是图像最基本的统计特征，是图像亮度分布的概率密度函数，反映了图像灰度值的分布情况。直方图是多种空间域处理技术的基础。直方图操作能有效地用于图像增强，如通过直方图均衡化处理，可使图像在整个灰度级范围内的分布均匀化，即在每个灰度级上都具有相同的像素点数，从而获得较好的视觉效果。另外，直方图固有的信息也可用在图像分割等其它图像处理的应用中。
灰度级变换技术可用g(x,y)=T[f(x,y)]的形式表示，其中f(x,y)为输入图像，g(x,y)为输出图像，T是对图像f进行某种处理的操作。由于(x,y)处的g值仅由f在该点处的亮度决定，T也称为一个亮度或灰度级变换函数，它与位置坐标(x,y)无关，所以通常写成如下的简化形式s=T®，其中r和s分别表示图像f和g在相应点(x,y)的亮度。利用灰度变换可调整一幅图像的明暗、对比度等。
三、实验内容
-
编写一个图像灰度直方图统计函数my_imhist，选择一幅图像利用my_imhist显示其直方图，将结果与MATLAB图像处理工具箱中提供的灰度直方图函数imhist的处理结果进行比较，并在同一窗口中显示出来。
实验代码：
a = imread('img1.jpg');
L = 256;
my = zeros(L,1);
[row,col] = size(a);
n = row*col;
for i = 1:row
for j = 1:col
num = a(i,j)+1;
my(num) = my(num)+1;
end
end
subplot(1,2,1);bar(my);title("my_imhist");
subplot(1,2,2);imhist(a);title("imhist");
实验结果：
- 利用以上编写的函数my_imhist或imhist，估算图像iris.tif中瞳孔的半径（以像素为单位）。
答：因为瞳孔是黑色的，可以大致认为图像中所有黑色的像素块都为瞳孔。根据my_imhist可以算出0-50之间的像素个数大约为13665个，然后圆形面积公式的R=（S/π）^（1/2），得R≈66。但此结果并不精确，因为眉毛也是黑色的。
- 、按照教材68页上的公式（4.1.6），编程实现图像的分段线性灰度变换.
实验代码：
Img = imread( 'img1.jpg' );
[m,n,k] = size(Img);
J = zeros(m,n,k);
a = 255/3;
b = 2*a;
M = 255;
c = a + 20;
d=M-c;
for i = 1 : m*n*k
if (Img(i) >= 0 && Img(i) < a)
J(i)=(c/a) * Img(i);
else
if (Img(i) >= a && Img(i) < b)
J(i) = (d-c)/(b-a) * (Img(i)-a) + c;
else
J(i)= (d-b)/(d-a) * (Img(i)-b) + d;
end
end
end
figure(1);
subplot(2,2,1);
imshow(Img);
title( '原图像' );
subplot(2,2,2);
imshow(uint8(J));
title( '变换后' );
subplot(2,2,3);
imhist(Img);
subplot(2,2,4);
imhist(uint8(J));
- 编写一个灰度图像的直方图均衡化函数（不可使用库函数）。（可使用的灰度等级数量不变即可），并对下图例题进行测试，给出测试结果。
实验代码：
A = [5,5,4,4,2,0,7;
5,5,4,2,0,0,7;
4,4,4,2,0,7,7;
1,3,3,0,7,7,6;
1,0,0,0,7,7,6;
7,7,7,7,7,7,6;
7,7,7,7,7,7,6];
[m,n] = size(A);
A1 = zeros(m,n);
for i = 1:m*n
A1(i) = A(i)/7;
end
subplot(2,2,1);
imshow(A1);
title('原图像');
subplot(2,2,3)
imhist(A1);
nk = zeros(m+1);
for i = 1:m
for j = 1:n
a = A(i,j)+1;
nk(a) = nk(a) + 1;
end
end
pr = zeros(m+1);
for i = 1:m+1
pr(i) = nk(i)/(m*n);
end
Sk = zeros(m+1);
Sk(1) = pr(1);
for i = 2:m+1
Sk(i) = Sk(i-1) + pr(i);
end
Sk_write_as = zeros(m+1);
for i = 1:m+1
a = 1;
for j = 0:m
if a >= abs(Sk(i) - (j/7))
Sk_write_as(i) = j/7;
a = Sk(i) - j/7;
end
end
end
a = Sk_write_as(3)*7;
A2 = zeros(m,n);
for i = 1:m
for j = 1:n
A2(i,j) = Sk_write_as(A(i,j)+1) * 7;
end
end
A3 = zeros(m,n);
for i = 1:m*n
A3(i) = A2(i) / 7;
end
subplot(2,2,2);
imshow(A3);
title('新图像');
subplot(2,2,4);
imhist(A3);
实验三 图像空间域平滑与锐化
一、实验目的与要求
- 加深对图像增强及边缘检测技术的感性认识，应用MATLAB工具箱自带的处理函数或自己编程完成相关的工作，分析处理结果，巩固所学理论知识。
- 熟练掌握空域滤波中常用的平滑和锐化滤波器，针对不同类型和强度的噪声，进行滤波处理，体会并正确评价滤波效果，了解不同滤波方式的使用场合，能够从理论上作出合理的解释。
二、实验相关知识
图像增强是指按特定的需要突出一幅图像中的某些有用信息，同时消弱或去除某些不需要的信息的处理方法，其主要目的是使处理后的图像对某些特定的应用比原来的图像更加有效。图像平滑与锐化处理是图像增强的主要研究内容。
图像分割是由图像处理过渡到图像分析的关键步骤，一般是按照一定的规则，把图像分成互不重叠的若干区域或子集，并提取出感兴趣目标的技术和过程。本实验主要研究利用边缘提取算子的分割方法。
和本实验有关的几个常用Matlab函数：
(1) imnoise：用于对图像生成模拟噪声，如：
j=imnoise(i,‘gaussian’,0,0.02) %在图像i上叠加均值为0、方差为0.02的高斯噪声，得到含噪图像j
j=imnoise(i,‘salt & pepper’,0.04) %在图像i上叠加密度为0.04的椒盐噪声，得到含噪图像j
(2) fspecial：用于产生预定义滤波器，如：
h=fspecial(‘average’,3); %产生3×3模板的均值滤波器
h=fspecial(‘sobel’); %产生sobel水平边缘增强的滤波器
可选项还有：'gaussian’高斯低通滤波器、'laplacian’拉普拉斯滤波器、'log’高斯拉普拉斯滤波器等
(3) imfilter、filter2、conv2：均是基于卷积的图像滤波函数，都可用于图像滤波，用法类似，如：
i=imread(‘p1.tif’);
h=fspecial(‘prewitt’); %产生prewitt算子的水平方向模板
j1=imfilter(i,h); %或者j2=filter2(h,i); 或者j3=conv2(i,h);
(4) medfilt2：用于图像的中值滤波，如：
j=medfilt2(i, [M N]); %对矩阵i进行二维中值滤波，邻域为M×N，缺省值为3×3
三、实验内容
- 图像平滑（去噪）：编写超限像素平滑法，灰度最相近的K个邻点平均法（函数名称可以自定义），并对上面实验1选择的噪声图片进行处理，显示处理前后的各个图像，分析不同方法对不同噪声的处理效果及其优缺点。
实验代码：
clear all;
I = imread('1.png');
[m,n] = size(I);
I1 = imnoise(I,'salt & pepper',0.04);
I2 = imnoise(I,'gaussian',0,0.02);
subplot(3,3,1);
imshow(I);
title('原图像');
subplot(3,3,4);
imshow(I1);
title('原图像添加椒盐噪声');
subplot(3,3,7);
imshow(I2);
title('原图像添加高斯噪声');
g1 = zeros(m,n);
g2 = zeros(m,n);
% 求出g（x）
for i = 1:m
for j = 1:n
if(i == 1 || j == 1 || i == m || j == n)
g1(i,j) = I1(i,j);
g2(i,j) = I2(i,j);
else
for x = -1:1
for y = -1:1
g1(i,j) = g1(i,j) + double(I1(i+x,j+y));
g2(i,j) = g2(i,j) + double(I2(i+x,j+y));
end
end
g1(i,j) = g1(i,j) / 9;
g2(i,j) = g2(i,j) / 9;
end
end
end
% 求g'(x)
for i = 1:m
for j = 1:n
if(abs(g1(i,j)-double(I1(i,j)))>64)
g1(i,j) = g1(i,j);
else
g1(i,j) = I1(i,j);
end
if(abs(g2(i,j)-double(I2(i,j)))>64)
g2(i,j) = g2(i,j);
else
g2(i,j) = I2(i,j);
end
end
end
subplot(3,3,5)
imshow(uint8(g1));
title('超限像素平滑法')
subplot(3,3,8)
imshow(uint8(g2));
title('超限像素平滑法')
%最相近的K个邻点平均法
A1 = zeros(1,9);
A2 = zeros(1,9);
K = 5;
for i = 1:m
for j = 1:n
if(i == 1 || j == 1 || i == m || j == n)
g1(i,j) = I1(i,j);
g2(i,j) = I2(i,j);
else
for x = -1:1 %找到当前像素点的3×3邻域（A)
for y = -1:1
A1((x+1)*3 + y+2) = I1(i+x,j+y);
A2((x+1)*3 + y+2) = I2(i+x,j+y);
end
end
[~,index] = sort(abs(A1-double(I1(i,j))));
AA = A1(index(1:K)); %%AA矩阵为最相近的K个点
g1(i,j) = sum(AA) / K;
[~,index] = sort(abs(A2-double(I2(i,j))));
AA = A2(index(1:K)); %%AA矩阵为最相近的K个点
g2(i,j) = sum(AA) / K;
end
end
end
subplot(3,3,6)
imshow(uint8(g1));
title('最相近K个邻点平均法')
subplot(3,3,9)
imshow(uint8(g2));
title('最相近K个邻点平均法');
实验结果：
思考题：
超限像素平滑算法对抑制椒盐噪声比较有效，对保护仅有微小灰度差的细节及纹理也有效。随着邻域增大，该算法去噪能力增强，但模糊程度也加大。
- 图像锐化:编写梯度锐化算法函数my_grad(method，T)；参数method可以是梯度算子、Roberts算子、Prewitt和Sobel算子，T是梯度阈值，根据参数method和T选用不同的算子和阈值进行锐化处理，用二值图像表示处理结果，小于T的用黑色表示，大于等于T的用白色表示。显示处理前后的各个图像，分析不同方法对锐化效果及其优缺点。
实验代码：
main.m
f = my_grad("tidu",0);
subplot(2,2,1);imshow(f);title("梯度算子");
f = my_grad("Roberts",0);
subplot(2,2,2);imshow(f);title("Roberts算子");
f = my_grad("Prewitt",0);
subplot(2,2,3);imshow(f);title("Prewitt算子");
f = my_grad("Sobel",0);
subplot(2,2,4);imshow(f);title("Sobel算子");
my_grad.m
function [f] = my_grad(method,T)
switch(method)
case "tidu"
f = tidu(T);
case "Roberts"
f = Roberts(T);
case "Prewitt"
f = Prewitt(T);
case "Sobel"
f = Sobel(T);
end
tidu.m
function [f] = tidu(T)
I = imread('wire.png');
%梯度算子
w1 = [-1,1;0,0];
w2 = [-1,0;1,0];
[m,n] = size(I);
fx = zeros(m,n);
fy = zeros(m,n);
for i = 1:m-1
for j = 1:n-1
for x = 1:2
for y = 1:2
fx(i,j) = fx(i,j) + double(I(i+x-1,j+y-1)) * w1(x,y);
fy(i,j) = fy(i,j) + double(I(i+x-1,j+y-1)) * w2(x,y);
end
end
end
end
f = abs(fx) + abs(fy);
for i = 1:m-1
for j = 1:n-1
if f(i,j) >= T
else
f(i,j) = I(i,j);
end
end
end
Prewitt.m
function [f] = Prewitt(T)
I = imread('wire.png');
%Prewitt算子
w1 = [-1,0,1;-1,0,1;-1,0,1];
w2 = [-1,-1,-1;0,0,0;1,1,1];
[m,n] = size(I);
fx = zeros(m,n);
fy = zeros(m,n);
for i = 2:m-1
for j = 2:n-1
for x = 1:3
for y = 1:3
fx(i,j) = fx(i,j) + double(I(i+x-2,j+y-2)) * w1(x,y);
fy(i,j) = fy(i,j) + double(I(i+x-2,j+y-2)) * w2(x,y);
end
end
end
end
f = abs(fx) + abs(fy);
for i = 1:m-1
for j = 1:n-1
if f(i,j) >= T
else
f(i,j) = I(i,j);
end
end
end
Roberts.m
function [f] = Roberts(T)
I = imread('wire.png');
%Roberts算子
w1 = [-1,0;0,1];
w2 = [0,-1;1,0];
[m,n] = size(I);
fx = zeros(m,n);
fy = zeros(m,n);
for i = 1:m-1
for j = 1:n-1
for x = 1:2
for y = 1:2
fx(i,j) = fx(i,j) + double(I(i+x-1,j+y-1)) * w1(x,y);
fy(i,j) = fy(i,j) + double(I(i+x-1,j+y-1)) * w2(x,y);
end
end
end
end
f = abs(fx) + abs(fy);
for i = 1:m-1
for j = 1:n-1
if f(i,j) >= T
else
f(i,j) = I(i,j);
end
end
end
Sobel.m
function [f] = Sobel(T)
I = imread('wire.png');
%Sobel算子
w1 = [-1,0,1;-2,0,2;-1,0,1];
w2 = [-1,-2,-1;0,0,0;1,2,1];
[m,n] = size(I);
fx = zeros(m,n);
fy = zeros(m,n);
for i = 2:m-1
for j = 2:n-1
for x = 1:3
for y = 1:3
fx(i,j) = fx(i,j) + double(I(i+x-2,j+y-2)) * w1(x,y);
fy(i,j) = fy(i,j) + double(I(i+x-2,j+y-2)) * w2(x,y);
end
end
end
end
f = abs(fx) + abs(fy);
for i = 1:m-1
for j = 1:n-1
if f(i,j) >= T
else
f(i,j) = I(i,j);
end
end
end
subplot(2,2,4);imshow(uint8(f));title("Sobel算子");
实验结果：
思考题：
Roberts算子优点是方法简单，缺点是对噪声敏感，常用于不含噪声的图像边缘增强。
Prewitt算子与Sobel都在检测边缘点时的同时具有抑制噪声的能力，检测出的边缘宽度至少为二像素缺点是会丢失一些细节信息，是边缘有一定的模糊。但由于Sobel算子的加权作用，其使边缘的模糊程度要稍低于Prewitt算子。
实验四 图像频域平滑与锐化
一、实验目的与要求
- 了解频域变换过程，掌握频域变换特点
- 熟练掌握频域滤波中常用的平滑和锐化滤波器，能够对不同要求的图像进行滤波处理，体会并正确评价滤波效果，了解不同滤波方式的使用场合，能够从理论上作出合理的解释。
二、实验相关知识
- 图像增强是指按特定的需要突出一幅图像中的某些有用信息，同时消弱或去除某些不需要的信息的处理方法，其主要目的是使处理后的图像对某些特定的应用比原来的图像更加有效。图像平滑与锐化处理是图像增强的主要研究内容。
- 和本实验有关的几个常用Matlab函数：
(1) imnoise：用于对图像生成模拟噪声，如：
j=imnoise(i,‘gaussian’,0,0.02) %在图像i上叠加均值为0、方差为0.02的高斯噪声，得到含噪图像j
j=imnoise(i,‘salt & pepper’,0.04) %在图像i上叠加密度为0.04的椒盐噪声，得到含噪图像j
(2) fspecial：用于产生预定义滤波器，如：
h=fspecial(‘average’,3); %产生3×3模板的均值滤波器
h=fspecial(‘sobel’); %产生sobel水平边缘增强的滤波器
可选项还有：'gaussian’高斯低通滤波器、'laplacian’拉普拉斯滤波器、'log’高斯拉普拉斯滤波器等
(3) imfilter、filter2、conv2：均是基于卷积的图像滤波函数，都可用于图像滤波，用法类似，如：
i=imread(‘p1.tif’);
h=fspecial(‘prewitt’); %产生prewitt算子的水平方向模板
j1=imfilter(i,h); %或者j2=filter2(h,i); 或者j3=conv2(i,h);
(4) fft2:二维快速傅里叶变换函数
(5) fftshift:中心变换函数
(6) abs:取绝对值或复数取幅值
三、实验内容
- 图像频域平滑（去噪）：使用自生成图像（包含白色区域，黑色区域，并且部分区域添加椒盐噪声），然后进行傅里叶变换，并且分别使用理想低通滤波器、巴特沃斯低通滤波器、指数低通滤波器和梯形低通滤波器（至少使用两种低通滤波器），显示滤波前后的频域能量分布图，空间图像。分析不同滤波器对噪声、边缘的处理效果及其优缺点。
实验代码：
img = imread('ren.png');
img = rgb2gray(img);
img_noise = imnoise(img,'salt & pepper',0.04);
subplot(4,2,1);imshow(img_noise);title('原图像');
F = fft2(double(img_noise));
F = fftshift(F);
[m,n] = size(F);
F1 = abs(F);
T = log(F1+1);
subplot(4,2,2);imshow(T,[]);title('傅里叶变换-log(abs(f)+1)');
m_mid = fix(m/2);
n_mid = fix(n/2);
IF = zeros(m,n);
h = zeros(m,n);
d0 = fix(m/5);
d = zeros(m,n);
for i = 1:m
for j = 1:j
d(i,j) = sqrt((i-m_mid)^2+(j-n_mid)^2);
end
end
%理想低通滤波器
for i = 1:m
for j = 1:n
if d(i,j) <= d0
h(i,j) = 1;
else
h(i,j) = 0;
end
IF(i,j) = h(i,j)*F(i,j);
end
end
F1 = abs(IF);
T = log(F1+1);
subplot(4,2,3);imshow(T,[]);title('理想低通滤波器-log(abs(f)+1)');
IF = ifftshift(IF);
IF = uint8(real(ifft2(IF)));
subplot(4,2,4);imshow(IF);title('理想低通滤波器,D0=行数/5');
IF = zeros(m,n);
h = zeros(m,n);
%巴特沃斯低通滤波器
n_btws = 2;
for i = 1:m
for j = 1:n
h(i,j) = 1/(1+0.414*(d(i,j)/d0)^(2*n_btws));
IF(i,j) = h(i,j)*F(i,j);
end
end
F1 = abs(IF);
T = log(F1+1);
subplot(4,2,5);imshow(T,[]);title('巴特沃斯低通滤波器-log(abs(f)+1)');
IF = ifftshift(IF);
IF = uint8(real(ifft2(IF)));
subplot(4,2,6);imshow(IF);title('巴特沃斯低通滤波器,n=2');
%梯形低通滤波器
d1 = 2*d0;
IF = zeros(m,n);
h = zeros(m,n);
for i = 1:m
for j = 1:n
if d(i,j) <= d0
h(i,j) = 1;
else
if d(i,j) > d1
h(i,j) = 0;
else
h(i,j) = (d(i,j)-d1)/(d0-d1);
end
end
IF(i,j) = h(i,j)*F(i,j);
end
end
F1 = abs(IF);
T = log(F1+1);
subplot(4,2,7);imshow(T,[]);title('梯形低通滤波器-log(abs(f)+1)');
IF = ifftshift(IF);
IF = uint8(real(ifft2(IF)));
subplot(4,2,8);imshow(IF);title('梯形低通滤波器,d1=2*d0');
实验结果：
思考题：
由于高频成分包含有大量的边缘信息，因此采用理想低通滤波器在去噪声的同时将会导致边缘信息损失而使图像边模糊。巴特沃斯低通滤波器的特性是连续性衰减，而不像理想滤波器那样陡峭变化，即明显的不连续性。因此采用该滤波器滤波在抑制噪声的同时，图像边缘的模糊程度大大减小，没有振铃效应产生。采用梯形低通滤波器滤波在抑制噪声的同时，图像边缘的模糊程度较用Butterworth滤波产生的大些，无明显的振铃效应。
- 图像频域平滑（锐化）：选择一幅图像，例如rice.png，分别使用理想高通滤波器、巴特沃斯高通滤波器、指数高通滤波器和梯形高通滤波器（至少使用两种高通滤波器），显示滤波前后的频域能量分布图，空间图像。分析不同滤波器处理效果及其优缺点。
实验代码：
clear all;
img = imread('ren.png');
img = rgb2gray(img);
subplot(4,2,1);imshow(img);title('原图像');
F = fft2(double(img));
F = fftshift(F);
[m,n] = size(F);
F1 = abs(F);
T = log(F1+1);
subplot(4,2,2);imshow(T,[]);title('傅里叶变换-log(abs(f)+1)');
m_mid = fix(m/2);
n_mid = fix(n/2);
IF = zeros(m,n);
h = zeros(m,n);
d0 = fix(m/10);
d = zeros(m,n);
for i = 1:m
for j = 1:n
d(i,j) = sqrt((i-m_mid)^2+(j-n_mid)^2);
end
end
%理想低通滤波器
for i = 1:m
for j = 1:n
if d(i,j) <= d0
h(i,j) = 0;
else
h(i,j) = 1;
end
IF(i,j) = h(i,j)*F(i,j);
end
end
F1 = abs(IF);
T = log(F1+1);
subplot(4,2,3);imshow(T,[]);title('理想高通滤波器-log(abs(f)+1)');
IF = ifftshift(IF);
IF = uint8(real(ifft2(IF)));
subplot(4,2,4);imshow(IF);title('理想高通滤波器,D0=行数/20');
实验结果：
1290
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
