matlab数字图像处理汇编.docx
- 文档编号:4890266
- 上传时间:2022-12-11
- 格式:DOCX
- 页数:22
- 大小:1.04MB
matlab数字图像处理汇编.docx
《matlab数字图像处理汇编.docx》由会员分享,可在线阅读,更多相关《matlab数字图像处理汇编.docx(22页珍藏版)》请在冰豆网上搜索。
matlab数字图像处理汇编
中国地质大学
课程设计
课程名称数字图像处理
教师姓名傅华明
学生姓名
学生学号
学生班级
数字图像处理报告
1.实验要求:
已知一图象p2-04-01,经过低通滤波得到其退化图象p2-04-02。
采用逆滤波方式将其重新复原。
低通滤波器采用两种:
(1)巴特沃茨滤波器
(2)高斯滤波器,其截止频率自行设定。
原理分析:
此题先进行傅立叶变换后高通滤波,再傅立叶反变换得到退化图象,再傅立叶变换,最后逆滤波。
这里要用到傅立叶变换,低通滤波及图象复原知识。
傅立叶变换
傅立叶变换是数字图像处理中应用最广的一种变换,其中图像增强、图像复原和图像分析与描述等,每一类处理方法都要用到图像变换,尤其是图像的傅立
叶变换。
离散傅立叶(Fourier)变换的定义:
二维离散傅立叶变换(DFT)为:
逆变换为:
式中,
在DFT变换对中,称为离散信号的频谱,而称为幅度谱,为相位角,功率谱为频谱的平方,它们之间的关系为:
图像的傅立叶变换有快速算法。
逆滤波:
则
1/H(u,v)称为逆滤波器。
对上式再进行傅立叶反变换可得到f(x,y)。
但实际上碰到的问题都是有噪声的,因而只能求F(u,v)的估计值:
然后再作傅立叶逆变换,得
程序:
%**********************显示图像********************%
f=imread('p2-04-01.bmp');
subplot(3,4,1);
imshow(f);title('原始图像');
c=imread('p2-04-02.bmp');
subplot(3,4,2);
imshow(c);title('低通滤波后退化图像');
%*******************二维傅里叶变化******************%
F=fft2(f);%求原始图像的傅里叶变化
F=fftshift(F);%将傅里叶变化零频率搬移到频谱中间
subplot(3,4,3);
imshow(log(abs(F)),[-110]);title('原始图像的傅里叶变换');
Q=fft2(c);%求退化图像的傅里叶变化
Q=fftshift(Q);%将傅里叶变化零频率搬移到频谱中间
subplot(3,4,4);
imshow(log(abs(Q)),[-110]);title('退化图像的傅里叶变换');
%******************理想频域低通滤波*****************%
clearall;
f=imread('p2-04-01.bmp');
F=fft2(f);%求原始图像的傅里叶变化
F=fftshift(F);%将傅里叶变化零频率搬移到频谱中间
[M,N]=size(F);%确定图像大小,M为行数,N为列数
m=fix(M/2);n=fix(N/2);%确定傅里叶变化原点(即直流部分),并且数据向0取整
d0=20;%设定低通滤波上限频率
fori=1:
M
forj=1:
N
d=sqrt((i-m)^2+(j-n)^2);
if(d<=d0)
h=1;
elseh=0;
end
H(i,j)=h*F(i,j);%滤波矩阵点乘
end
end
subplot(3,4,5);
imshow(H);title('理想频域低通滤波器处理');
H=ifftshift(H);
J1=ifft2(H);
J2=uint8(real(J1));%提取J1的实部,并将该数据定义8位无符号整数
subplot(3,4,6);
imshow(J2);title('理想低通滤波退化后图像');
%******************理想低通滤波复原*****************%
FR=fft2(J1);
FR=fftshift(FR);%退化图像经过二维傅里叶变换
[M,N]=size(FR);%确定图像大小,M为行数,N为列数
m=fix(M/2);n=fix(N/2);%确定傅里叶变化原点(即直流部分),并且数据向0取整
d1=20;%设定低通滤波上限频率
fori=1:
M
forj=1:
N
d=sqrt((i-m)^2+(j-n)^2);
if(d>=d1)
M(i,j)=F(i,j);
else
M(i,j)=FR(i,j);
end
end
end
M=ifftshift(M);
J3=ifft2(M);%逆滤波复原图像
J4=uint8(real(J3));%提取J3的实部,并将该数据定义8位无符号整数
subplot(3,4,9);
imshow(J4);title('理想低通滤波退化后图像的复原图');
%******************巴特沃斯低通滤波*****************%
clearall;
f=imread('p2-04-01.bmp');
F=fft2(f);%求原始图像的傅里叶变化
F=fftshift(F);%将傅里叶变化零频率搬移到频谱中间
n=2;%设定巴特沃斯滤波器阶数
[M,N]=size(F);%确定图像大小,M为行数,N为列数
m=fix(M/2);n=fix(N/2);%确定傅里叶变化原点(即直流部分),并且数据向0取整
d0=30;%设定巴特沃斯低通滤波上限频率
fori=1:
M
forj=1:
N
d=sqrt((i-m)^2+(j-n)^2);
h1(i,j)=1/(1+(d/d0)^(2*n));
H1(i,j)=h1(i,j)*F(i,j);%滤波矩阵点乘
end
end
subplot(3,4,7);
imshow(H1);title('巴特沃斯低通滤波器处理');
H1=ifftshift(H1);
J1=ifft2(H1);
J2=uint8(real(J1));%提取J1的实部,并将该数据定义8位无符号整数
subplot(3,4,8);
imshow(J2);title('巴特沃斯低通滤波退化后图像');
%****************巴特沃斯低通滤波复原***************%
FR=fft2(J1);
FR=fftshift(FR);%退化图像经过二维傅里叶变换
[M,N]=size(FR);%确定图像大小,M为行数,N为列数
m=fix(M/2);n=fix(N/2);%确定傅里叶变化原点(即直流部分),并且数据向0取整
d1=20;%设定低通滤波上限频率
fori=1:
M
forj=1:
N
d=sqrt((i-m)^2+(j-n)^2);
if(d>=d1)
M(i,j)=F(i,j);
else
M(i,j)=FR(i,j);
end
end
end
M=ifftshift(M);
J3=ifft2(M);%逆滤波复原图像
J4=uint8(real(J3));%提取J3的实部,并将该数据定义8位无符号整数
subplot(3,4,11);
imshow(J4);title('巴特沃斯低通滤波退化后图像的复原图');
%************退化图像巴特沃斯低通滤波复原************%
clearall;
f=imread('p2-04-01.bmp');
F=fft2(f);%求原始图像的傅里叶变化
F=fftshift(F);%将傅里叶变化零频率搬移到频谱中间
c=imread('p2-04-02.bmp');
Q=fft2(c);%求退化图像的傅里叶变化
Q=fftshift(Q);%将傅里叶变化零频率搬移到频谱中间
[M,N]=size(Q);%确定图像大小,M为行数,N为列数
m=fix(M/2);n=fix(N/2);%确定傅里叶变化原点(即直流部分),并且数据向0取整
d0=2;%设定低通滤波上限频率
fori=1:
M
forj=1:
N
d=sqrt((i-m)^2+(j-n)^2);
if(d>=d0)
M(i,j)=F(i,j);
else
M(i,j)=Q(i,j);
end
end
end
M=ifftshift(M);
J3=ifft2(M);%逆滤波复原图像
J4=uint8(real(J3));%提取J3的实部,并将该数据定义8位无符号整数
subplot(3,4,12);
imshow(J4);title('退化图像巴特沃斯低通滤波的复原图');
运行结果如下
2.实验要求:
已知p2-08为原图象,将该图象与平滑函数
卷积产生模糊,再叠加零均值,方差分别为8,16和32的高斯随机噪声得到的1组待复原的图象。
采用
(1)维纳滤波方法进行复原;
(2)约束最小二乘方滤波方法进行复原;并进行比较。
代码:
closeall;
clear;
clc;
f=imread('p2-08.tif');
[row,line]=size(f);
F=fft2(f);
%figure;
%imshow(F);
%title('原图像频谱');
r=1:
row;
c=1:
line;
[C,R]=meshgrid(c,r);
%%%%%%%%%%%%%%用向量优化法,加快计算速度;实质是以空间换时间
%%%%%%%%%%%%%%%%%%下面的RC得到的均是二维数组,计算需满足二维数组计算法则
h=exp(sqrt(R.*R+C.*C)/240);
%%%%%%%%%h=fftshift(h);
H=fft2(h);
%imshow(H);
%title('传递函数频谱');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%方法1imfilter卷积方法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%gimfilter=abs(imfilter(f,h));
%%%gimfilter=mat2gray((gimfilter));
%%%%%%%%%%%%%方法2g1=(ifft2(F.*fft2(h)));%%%%%%%%%%%%%%%%%%%%%
g=(ifft2((F.*H)));
%%%g=mat2gray(abs(g1));
%%%g=ifft2(fft2(g1)./fft2(h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%方法3c=conv2(a,b)%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%方法4有扩展的卷积计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%h(2*row-1,2*line-1)=0;
%%H=fft2(h);
%%g=dftfilt(f,H);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%倘若这句加进去后,下面得到的信噪比将相当大
%g=mat2gray(g)*255;
M=0;%%%%%%平均值
%%%%%%%%%%%%V=8;%%%%%%方差s
%%%%%%%%%%%%n
noise8=imnoise(zeros(size(f)),'gaussian',M,8);
noise16=imnoise(zeros(size(f)),'gaussian',M,16);
noise32=imnoise(zeros(size(f)),'gaussian',M,32);
%%%%%%%%%%%%%%%%%%%%加入高斯噪声后的图像%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fnoise8=g+(noise8);
fnoise16=g+noise16;
fnoise32=g+noise32;
%%%%%%%%%%%%%%滤波%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NPnoise8=abs(fftn(noise8)).^2;%%%%%%%%%%%噪声功率
NPnoise16=abs(fftn(noise16)).^2;
NPnoise32=abs(fftn(noise32)).^2;
%%%%%%%%%%%%%%%%原始图像功率谱%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
IP=abs(fftn(f)).^2;
IPnoise8=fft2(fnoise8);
IPnoise16=fft2(fnoise16);
%%%%%%%%%%%%%%%%对图像的复原%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fnoise8=((((abs(H)).^2)./H)./(abs(H)).^2+NPnoise8./IP).*IPnoise8;
Fnoise16=((((abs(H)).^2)./H)./(abs(H)).^2+NPnoise16./IP).*IPnoise16;
Fnoise32=((((abs(H)).^2)./H)./(abs(H)).^2+NPnoise32./IP).*fft2(fnoise32);
%%max(g2(:
))
%%min(g2(:
))
subplot(231);
imshow(f)
title('原图像');
subplot(232);
imshow(mat2gray(g))%%%%灰度归一化处理
title('模糊化的图像');
subplot(233);
imshow(mat2gray(fnoise8))
title('加方差为8的高斯噪声的图像');
subplot(234);
imshow(mat2gray(fnoise16))
title('加方差为16的高斯噪声的图像');
subplot(235);
imshow(mat2gray(fnoise32))
title('加方差为32的高斯噪声的图像');
figure;
subplot(221);
imshow(f)
title('原图像');
subplot(222);
imshow(mat2gray(ifft2(Fnoise8)))
title('复原加方差为8的高斯噪声的图像');
subplot(223);
imshow(mat2gray(ifft2(Fnoise16)))
title('复原加方差为16的高斯噪声的图像');
subplot(224);
imshow(mat2gray(ifft2(Fnoise32)))
title('复原加方差为32的高斯噪声的图像');
可得如下图形:
在做完卷积后,产生的模糊图像。
图像已经变得相当模糊,难以分辨。
下图是按题目所给要求,进行的对方差为8,16,32等等几种图像进行的复原。
若把g=mat2gray(g)加上产生的恢复效果,富原出来的效果并不好,难以分辨。
下面是去掉g=mat2gray(g)时,进行的图像复原,得到的复原图效果较好。
实验时发现即使噪声的方差达到10的7次方仍能较好的复原原图像。
下面是方差为16时的高斯噪声的平均功率谱大小,和原图像的功率谱大小。
可知其信噪比相当大。
3.图象几何教正
1,将p3-01校正为黑底白三角形。
2,将p3-02校正为黑底白菱形。
1.选取校正点
orthophoto=J;
unregistered=J;
cpselect(unregistered,orthophoto);
其中,orthophoto是基准图像,即正投影图像,unregistered为倾斜图像,cpelect函数为选取校正点,如下图,选取4对校正点
2.几何校正
input_points=evalin('base','input_points');
base_points=evalin('base','base_points');
mytform=cp2tform(input_points,base_points,'projective');
set(handles.axes1,'HandleVisibility','ON');
axes(handles.axes1);
imshow(imtransform(J,mytform));
Q=imtransform(J,mytform);
其中,cp2tform函数为通过四对点来构建投影变换结构mytform,imtransform函数将变形的图像J通过mytform结构校正,得到校正后的图像Q;如下图
程序
clc
clear
f=imread('p3-01.tif');
basepoints=[103;10259;2670;267259];;%校正图像匹配点坐标
inputpoints=[103;61230;21527;267259]%失真图像控制点坐标
tform=cp2tform(inputpoints,basepoints,'projective');%估计线性变换方程
g=imtransform(f,tform);
figure
(1)
subplot(1,2,1),
imshow(f);title('原始图像');
subplot(1,2,2),
imshow(g);title('校正图像');
结果
方法二
选取的校正点不合适,导致校正后的图像严重变形。
因为不规则的三角形可利用MPP算法去把细微的曲边拉直,所以可以考虑提取不规则三角形的轮廓运用MPP算法进行拉直而后再将得到的轮廓填充白色。
程序
clc;
B=imread('p3-01.tif');
[x,y]=minperpoly(B,2);
b2=connectpoly(x,y);
b=boundaries(B);
b=b{1};
[M,N]=size(B);
xmin=min(b(:
1));
ymin=min(b(:
2));
B2=bound2im(b2,M,N,xmin,ymin);
imshow(B);
title('原图像');
figure
imshow(B2);
title('修正图像');
结果
4.根据傅立叶反变换实施图象重建。
原理:
图像f(x,y)的傅立
其对x轴的投影可表示为
对其进行傅立叶变换,得
可见f(x,y)向x轴投影的傅立叶变换是和f(x,y)的傅立叶变换沿v=0的断面一致的。
因为傅立叶变换是正交变换,所以这一点对任意的方向都是成立的。
若对多个方向上的投影数据分别进行傅立叶变换,就可求出沿着与这个方向相同的直线上的F(u,v)。
如果把由它们计算出的F(u,v)进行傅立叶逆变换,就得到重建的图像f(x,y)。
因为投影数据的傅立叶变换得到的是极坐标形式的F(u,v),因此为了求得在直角坐标系中的F(u,v),就必须在(u,v)空间进行内插,或者按照极坐标进行傅立叶逆变换,在图像空间进行内插,得到重建的图像。
代码:
clc;
closseall;
f=zeros(15,15);
f(3:
3,3:
12)=1;
f(4:
4,3:
13)=1;
f(5:
13,3:
5)=1;
f(5:
5,12:
13)=1;
f(6:
6,6:
13)=1;
f(7:
7,6:
12)=1;
f(8:
8,6:
8)=1;
f(9:
9,7:
9)=1;
f(10:
10,8:
10)=1;
f(11:
11,9:
11)=1;
f(12:
12,10:
12)=1;
f(13:
13,11:
13)=1;
subplot(1,2,1);
imshow(f);
title('构建图像');
theta1=0:
179;
[R1,xp]=radon(f,theta1);
I1=iradon(R1,1);
subplot(1,2,2);
imshow((I1));
title('重构图像');
运行结果
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 数字图像 处理 汇编