图像处理实验2代码

2018-07-23 21:07:35 weixin_39569242 阅读数 2550

一、实验目的

(1)了解图像复原的目的及意义,加深对图像复原理论的认识。

(2)掌握维纳滤波复原基本原理。

(3)掌握约束最小二乘方复原方法。

(4)掌握盲解卷积复原方法

二、实验内容

 (1)维纳滤波复原。

 (2)约束最小二乘方复原

 (3)盲解卷积复原

三、实验代码及结果、分析

(1)维纳滤波复原

  • 代码:

I=imread('E:\大三课件\大三下\数字图像处理\实验\实验5\lena.jpg');

I=rgb2gray(I);

noise=5*randn(size(I));      

noise = noise - min(min(noise));

J  = double(I) +noise;

R1=wiener2(J,[10 10]);    %未知噪声

R2=wiener2(J,[10 10],noise);  %已知噪声分布

figure      

subplot(2,2,1),imshow(uint8(I));title('原始图像');      

subplot(2,2,2),imshow(uint8(J));title('退化图像');      

subplot(2,2,3),imshow(uint8(R1));title('盲复原');

subplot(2,2,4),imshow(uint8(R2));title('非盲复原');

  • 结果:

  • 分析:
  1. 维纳滤波又叫最小均方差滤波,它的目标是找到一个原图像f的估计图像f',使得它们之间的均方误差最小。  e^2 = E((f-f')^2)
  2. randa()函数是产生随机噪声,nosize应该是一个矩阵,min(min(noise))表示先取每列中的最小值,再取最小值中的最小值
  3. 关于维纳滤波函数有wiener2和deconvwnr,分别适用于灰度图像和彩色图像
  4. R1为未知噪声的分布,R2为已知噪声分布

 

(2)约束最小二乘方复原方法

  • 代码:

I=imread('E:\大三课件\大三下\数字图像处理\实验\实验5\lena.jpg');

I=rgb2gray(I);I = im2double(I);[hei,wid,~] = size(I);

PSF = fspecial('motion', 21, 11);%摄像物体逆时针方向以11角度运动了21个像素

blurred = imfilter(I, PSF, 'conv', 'circular');%运动模糊图像

% 逆滤波结果

If = fft2(blurred);Pf = psf2otf(PSF,[hei,wid]);deblurred =ifft2(If./Pf);

% 运动模糊+高斯噪声

noise_mean = 0;noise_var = 0.00001;

blurred_noisy = imnoise(blurred, 'gaussian',noise_mean, noise_var);

% 约束最小二乘法

p = [0 -1 0;-1 4 -1;0 -1 0];P = psf2otf(p,[hei,wid]);

gama = 0.002;%γ设置成0.001左右会有比较好的效果If = fft2(blurred_noisy);numerator = conj(Pf);

denominator = Pf.^2 + gama*(P.^2);

deblurred2 = ifft2( numerator.*If./ denominator );

subplot(2,3,1),imshow(I);title('原始图像');

subplot(2,3,2), imshow(blurred); title('运动模糊图像');

subplot(2,3,3), imshow(deblurred); title('运动模糊复原')

subplot(2,3,4), imshow(blurred_noisy),title('运动模糊+噪声图像')

subplot(2,3,5), imshow(deblurred2),title('无约束最小二乘法复原');

subplot(2,3,6); imshow(deconvreg(blurred_noisy, PSF,0)); title('deconvreg函数复原');%自带的去模糊deconvreg函数

  • 结果:

  • 分析:

1、fspecial('type', len,theta)函数,用于预定于滤波算子,type=motion为运动算子,theta参数表示逆时针方向的角度,len表示运动了len个像素

2、约束最小二乘方滤波效果比自带的去模糊deconvreg函数好

3、约束最小二乘方滤波对高噪声和中等噪声产生结果要好于维纳滤波,对于低噪声两种滤波产生结果基本相等

 

(3)盲解卷积复原方

  • 代码:

I=checkerboard(8);

PSF=fspecial('gaussian',7,10);

V=.0001;

BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);

WT=zeros(size(I));

WT(5:end-4,5:end-4)=1;

INITPSF=ones(size(PSF));

FUN=inline('PSF+P1','PSF','P1');

[J P]=deconvblind(BlurredNoisy,INITPSF,20,10*sqrt(V),WT,FUN,0);

subplot(131);imshow(I);title('原图像');

subplot(132);imshow(BlurredNoisy);title('运动模糊+噪声图像');

subplot(133);imshow(J);title('盲去卷积');

  • 结果:

  • 分析:
  1. 盲解卷积复原是指在没有图像退化先验知识、对退化系统了解不足的条件下,通过观察退化图像的多个图像以某种方式抽出退化信息,进行图像复原
  2. 函数imnoise 是表示添加噪声污染一幅图像,叫做噪声污染图像函数,g=imnoise(f,'localvar',V)将均值为0,局部方差为V的高斯噪声添加到图像f上,其中V是与f大小相同的一个数组,它包含了每一个点的理想方差值
  3. Inline()函数:定义一个内置函数,本质上说跟function干的是一样的事,只不过它可以直接内嵌在命令行里,不用另外单独定义function
2017-08-01 12:08:26 luoluonuoyasuolong 阅读数 16078

实验目的

实验一:图像增强

了解图像增强的目的及意义,加深对图像增强的感性认知
1. 掌握直接灰度变换的图像增强方法
2. 掌握灰度直方图的概念及其计算方法
3. 掌握直方图均衡化合直方图规定化得计算过程

实验二:图像平滑和锐化

  1. 利用临域平均法,临域加权平均法,中值滤波对图像进行滤波处理
  2. 掌握彩色图像的滤波处理
  3. 利用梯度方法对图像进行锐化处理
  4. 结合图像的平滑和锐化增强图像

 实验三:二维傅里叶变换

对二维的傅立叶变换有一定的认识,了解图像中傅立叶变换的意义 。

实验原理

图像增强

图像增强的方法是通过一定手段对原图像附加一些信息或变换数据,有选择地突出图像中感兴趣的特征或者抑制(掩盖)图像中某些不需要的特征,使图像与视觉响应特性相匹配。
三. 实验内容及结果分析。图像增强是从像素到像素的操作,是以预定的方式改变图像的灰度直方图。有时又称为对比度增强,灰度变换。点运算不可能改变图像内的空间关系,输出像素的灰度值由输入像素的值决定。其作用:对比度增强:扩展感兴趣特征的对比度。
首先,解释下什么是对比度:通俗的讲,就是明暗的对比程度,通常表现了图像画质的清晰程度,实际上,对比度反映了一幅图像中明暗区域最亮的白和最暗的黑之间不同亮度层级的测量,即指一幅图像灰度反差的大小,差异范围越大代表对比越大,差异范围越小代表对比
越小,好的对比率120:1就可容易地显示生动、丰富的色彩,当对比率高达300:1时,便可支持各阶的颜色。但对比率遭受和亮度相同的困境,现今尚无一套有效又公正的标准来
衡量对比率,所以最好的辨识方式还是依靠使用者眼睛。

直方图均衡化:

直方图均衡化处理的“中心思想”是把原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。直方图均衡化就是对图像进行非线性拉伸,重新分配图像像素值,使一定灰度范围内的像素数量大致相同。直方图均衡化就是把给定图像的直方图分布改变成“均匀”分布直方图分布。缺点:
  1)变换后图像的灰度级减少,某些细节消失;
  2)某些图像,如直方图有高峰,经处理后对比度不自然的过分增强。
  直方图均衡化是图像处理领域中利用图像直方图对对比度进行调整的方法。
  这种方法通常用来增加许多图像的局部对比度,尤其是当图像的有用数据的对比度相当接近的时候。通过这种方法,亮度可以更好地在直方图上分布。这样就可以用于增强局部的对比度而不影响整体的对比度,直方图均衡化通过有效地扩展常用的亮度来实现这种功能。
  这种方法对于背景和前景都太亮或者太暗的图像非常有用,这种方法尤其是可以带来X光图像中更好的骨骼结构显示以及曝光过度或者曝光不足照片中更好的细节。这种方法的一个主要优势是它是一个相当直观的技术并且是可逆操作,如果已知均衡化函数,那么就可以恢复原始的直方图,并且计算量也不大。这种方法的一个缺点是它对处理的数据不加选择,它可能会增加背景杂讯的对比度并且降低有用信号的对比度。
  直方图均衡化的基本思想是把原始图的直方图变换为均匀分布的形式,这样就增加了象素灰度值的动态范围从而可达到增强图像整体对比度的效果。
设r和s人别表示原图像灰度级和经过直方图均衡化后的图像灰度级。为便于讨论,对r和s进行归一化,使0r,s1 对于一幅给定的图像归一化后灰度级分布在0r1 范围内。对[0, 1]区间内的任何一个值进行如下变换:

s=T(r)

该变换应该满足条件
1)对于0r1 ,有 0s1
2)在0r1 区间内设r的概率密度为 Prr,s的概率密度可由Prr 求出
Ps(s)=(Pr(r)drds|r=T1(s))

假设变换函数为
s=T(r)=r0Pr(w)dw

式中:w是积分变量,而 r0Pr(w)dw就是r的累积分布函数。
离散形式的直方图均衡化
设一幅图像的像元数为n,公有m个灰度级, nk代表灰度级为r的像元的数目,则第k个灰度级出现的概率可表示为:
Pr(rk)=nkn,0rk1,k=0,1,...,m1
变换函数 改写为:
sk=T(rk)=j=0kPr(rj)=j=0knjn0rk1,k=0,1,...,m1

均衡化后的图像灰度值可直接由原图像的直方图算出。

直方图规定化

在实际应用中,希望能够有目的地增强某个灰度区间的图像, 即能够人为地修正直方图的形状, 使之与期望的形状相匹配,这就是直方图规定化的基本思想。换句话说,希望可以人为地改变直方图形状,使之成为某个特定的形状,直方图规定化就是针对上述要求提出来的一种增强技术,它可以按照预先设定的某个形状来调整图像的直方图。直方图规定化是在运用均衡化原理的基础上,通过建立原始图像和期望图像之间的关系,选择地控制直方图,使原始图像的直方图变成规定的形状,从而弥补了直方图均衡不具备交互作用的特性。
直方图规定化步骤:
1. 其增强原理是先对原始的直方图均衡化:S = T(r)
2. 找到一种映射有原来的直方图映射到新的目的直方图。
3. 将规定化的直方图进行均衡化。v = G(z)
3. 由于都是均衡化,故令 S = v,则:z = G-1(v) = G-1[T(r)]

图像平滑和锐化

图像平滑是指用于突出图像的宽大区域、低频成分、主干部分或抑制图像噪声和干扰高频成分,使图像亮度平缓渐变,减小突变梯度,改善图像质量的图像处理方法。图像平滑的方法包括:插值方法,线性平滑方法,卷积法等等。这样的处理方法根据图像噪声的不同进行平滑,比如椒盐噪声,就采用线性平滑方法!

均值滤波器

均值滤波器(averaging filter)是消除噪声的最简单的方法。原理:使用某像素周围mxn像素范围内的平均值来置换该像素值。通过使图像模糊,达到看不到细小噪声的目的。不良反应:使用这种方法,在噪声被消除的同时,目标图像也变模糊了。

中值滤波

消除噪声最好的结果是,在消除噪声的同时,图像边缘完好的保留。中值滤波能够比较好的实现这一点。原理:查看mxn邻域内的像素灰度,按照从小到大的顺序进行排列,结果取中间值。
中值滤波器与均值滤波器比较的优势:在均值滤波器中,由于噪声成分被放入平均计算中,所以输出受到了噪声的影响,但是在中值滤波器中,由于噪声成分很难选上,所以几乎不会影响到输出。因此同样用3x3区域进行处理,中值滤波消除的噪声能力更胜一筹。中值滤波无论是在消除噪声还是保存边缘方面都是一个不错的方法。

锐化原理

在图像增强过程中,通常利用各类图像平滑算法消除噪声,图像的常见噪声主要有加性噪声、乘性噪声和量化噪声等。一般来说,图像的能量主要集中在其低频部分,噪声所在的频段主要在高频段,同时图像边缘信息也主要集中在其高频部分。这将导致原始图像在平滑处理之后,图像边缘和图像轮廓模糊的情况出现。为了减少这类不利效果的影响,就需要利用图像锐化技术,使图像的边缘变得清晰。图像锐化处理的目的是为了使图像的边缘、轮廓线以及图像的细节变得清晰,经过平滑的图像变得模糊的根本原因是因为图像受到了平均或积分运算,因此可以对其进行逆运算(如微分运算)就可以使图像变得清晰。微分运算是求信号的变化率,由傅立叶变换的微分性质可知,微分运算具有较强高频分量作用。从频率域来考虑,图像模糊的实质是因为其高频分量被衰减,因此可以用高通滤波器来使图像清晰。但要注意能够进行锐化处理的图像必须有较高的性噪比,否则锐化后图像性噪比反而更低,从而使得噪声增加的比信号还要多,因此一般是先去除或减轻噪声后再进行锐化处理。
2.4 图像锐化的几种方式
1) 水平竖直直梯度
2) Robert 梯度
3) Sobel 算子
4) 二阶梯度
在图像处理中,一阶微分是通过梯度算法来实现的,对于一幅图像用函数f(x,y)表示,定义在f(x,y)在点(x,y)处的梯度是一个矢量,定义为:

G[f(x,y)]=[dfdx,dfdy]

梯度的方向在函数f(x,y)最大变化率的方向上,梯度的幅度G[f(x,y)]可以由以下公式算出:
G[f(x,y)]=([(dfdx)2+(dfdy)2])
由上式可知:梯度的数值就是f(x,y)在其最大变化率方向上的单位距离所增加的量。对于数字图像而言,微分可用差分来近似。因此上式可以写成:
G(f(i,j))=|f(i,j)f(i+1,j)|+|f(i,j)f(i,j+1)|
这种方法就是上面的水平竖直梯度法。

彩色图像的平滑

彩色图像平滑处理和灰度图像的平滑处理类似,差别在于彩色图像需要对各个通道依次进行平滑操作,使得每一个通道对应的图像都经过平滑处理,时间复杂度要远比灰度图像的高。还有就是彩色图像的噪声和灰度图的噪声有所差别,有时候需要区别的对待。

二维傅立叶变换

傅里叶提出任何周期函数都可以表示为不同频率的正弦和/或余弦和的形式,每个正弦和/或余弦乘以不同的系数(傅里叶级数)。图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯度.在噪声点和图像边缘处的频率为高频。
傅立叶变换在实际中有非常明显的物理意义,设f是一个能量有限的模拟信号,则其傅立叶变换就表示f的谱。从纯粹的数学意义上看,傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看,傅立叶变换是将图像从空间域转换到频率域,其逆变换是将图像从频率域转换到空间域。换句话说,傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数,傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数.
傅立叶变换以前,图像(未压缩的位图)是由对在连续空间(现实空间)上的采样得到一系列点的集合,我们习惯用一个二维矩阵表示空间上各点,则图像可由z=f(x,y)来表示。由于空间是三维的,图像是二维的,因此空间中物体在另一个维度上的关系就由梯度来表示,这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度?因为实际上对图像进行二维傅立叶变换得到频谱图,就是图像梯度的分布图,当然频谱图上的各点与图像上各点并不存在一一对应的关系,即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点,实际上图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的频率的大小(可以这么理解,图像中的低频部分指低梯度的点,高频部分相反)。一般来讲,梯度大则该点的亮度强,否则该点亮度弱。这样通过观察傅立叶变换后的频谱图,也叫功率图,我们首先就可以看出,图像的能量分布,如果频谱图中暗的点数更多,那么实际图像是比较柔和的(因为各点与邻域差异都不大,梯度相对较小),反之,如果频谱图中亮的点数多,那么实际图像一定是尖锐的,边界分明且边界两边像素差异较大的。对频谱移频到原点以后,可以看出图像的频率分布是以原点为圆心,对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外,还有一个好处,它可以分离出有周期性规律的干扰信号,比如正弦干扰,一副带有正弦干扰,移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心,对称分布的亮点集合,这个集合就是干扰噪音产生的,这时可以很直观的通过在该位置放置带阻滤波器消除干扰。
图像是两个参数的函数,通过一组正交函数的线性组合可以将其分解,而傅里叶就是通过谐波函数来分解的。而对于离散傅里叶变换,傅里叶变换的条件是存在的。
傅里叶变换进行图像处理有几个特点
1. 直流成分F(0,0)等于图像的平均值;
2. 能量频谱|F(u,v)|^2完全对称于原点;其中F=PfQ, f表示原图,P和Q都是对称的实正交矩阵,这个公式表示傅里叶变换就是个正交矩阵的正交变换
3.图像f平移(a,b)后,F只有exp[-2pij(au/M+bv/M)]的相位变化,能量频谱不发生变化。
4. 图像f自乘平均等于能量频谱的总和,f的分散等于能量频谱中除直流成分后的总和。
5.图像f(x,y)和g(x,y)的卷积h(x,y)=f(x,y)*g(x,y)的傅里叶变换H(u,v)等于f(x,y)和g(x,y)各自的傅里叶变换的乘积。

图像中的每个点通过傅里叶变换都成了谐波函数的组合,也就有了频率,这个频率则是在这一点上所有产生这个灰度的频率之和,也就是说傅里叶变换可以将这些频率分开来。当想除去图像背景时,只要去掉背景的频率就可以了。

在进行傅里叶变换时,实际上在某一特定的频率下,计算每个图像位置上的乘积。什么乘积呢,就是f(x,y)exp[-j2pi(ux+vy)],然后计算下一个频率。这样就得到了频率函数。
也就是说,我们看到傅里叶变换的每一项(对每对频率u,v,F(u,v)的值)是由f(x)函数所有值的和组成。f(x)的值与各种频率的正弦值和余弦值相乘。因此,频率u, v决定了变换的频率成分(x, y也作用于频率,但是它们相加,对频率有相同的贡献)。

通常在进行傅里叶变换之前用(-1)^(x+y)乘以输入的图像函数,这样就可以将傅里叶变换的原点F(0,0)移到(M/2,N/2)上。

每个F(u,v)项包含了被指数修正的f(x,y)的所有值,因而一般不可能建立图像特定分量和其变换之间的联系。然而,一般文献通常会有关于傅里叶变换的频率分量和图像空间特征之间联系的阐述。变换最慢的频率成分(u=v=0)对应一幅图像的平均灰度级。当从变换的原点移开时,低频对应着图像的慢变换分量,较高的频率开始对应图像中变化越来越快的灰度级。这些事物体的边缘和由灰度级的突发改变(如噪声)标志的图像成分。

实验内容及结果分析

直方图均衡化

从文件中读取一幅图像,转换成灰度图,将其进行均衡化处理,最后对比原始灰度图和均衡化后的结果。
下面的函数是我自己写的均衡化操作函数
核心函数:

1.  function re = HistgramEqualization(gray)  
2.  %直方图均衡化的函数  
3.  [m, n] = size(gray);  
4.  count = zeros(1, 256);  
5.  for i = 1 : m  
6.      for j = 1 : n  
7.        count(gray(i, j) + 1) = count(gray(i, j) + 1) + 1;    
8.      end  
9.  end  
10. ratio = count ./ (m * n);  
11. new_ratio = zeros(1, 256);  
12. % 统计各个灰度级的数量  
13. for i = 1 : 256  
14.     for j = 1 : i  
15.         new_ratio(i) = new_ratio(i) + ratio(j);  
16.     end  
17. end  
18. %新的比例  
19. new_count = floor(new_ratio .* 256);  
20. for i = 1 : m  
21.     for j = 1 : n  
22.       re(i, j) = new_count(gray(i, j) + 1);    
23.     end  
24. end  
25. end  

这个函数是自己编写的一个将原始图像进行直方图均衡化的函数。
原始灰度图
原始灰度直方图

均衡化后的图像
均衡化后的直方图

从均衡化后的结果可以看出来,图像的对比度略有提升,对应的直方图的分布也相对均匀。但缺点是某些像素值是没有的,也就是灰度值少了一些。

直方图规定化

直方图规定化是给出一个映射使直方图按所需映射要求进行变换, 测试程序的代码如下所示:

“`
1. I = imread(‘123.jpg’);
2. I = rgb2gray(I);
3. hgram=0:255;
4. J = histeq(I,hgram);
5. figure(),imshow(I),title(‘原始图像’);
6. figure(), imshow(J),title(‘直方图规定化后图像’);
7. figure(), imhist(I),title(‘原始图像直方图’);
8. figure(), imhist(J),title(‘规定化后图像直方图’);
““
效果如下图所示
规定化后的效果
规定化后的直方图

邻域平均法对图像进行平滑

首先在原始图像上加上一些噪声,然后利用邻域平均法对图像进行平滑处理。
下面是我自己写的对灰度图进行邻域平均法进行平滑的函数,参数gray是待输入的灰度图像矩阵,参数N表示模板的维数,N只能为3或者5,当N为3时,模板为3×3的,程序运行时可根据不同情况选取不同的模板。

1.  function re = AverageBoxFilter(gray, N)  
2.  %均值滤波函数  N 是 模板的维数  
3.  [m, n] = size(gray);  
4.  Answer = zeros(m, n);  
5.  out = zeros(m + (N - 1), n + (N - 1));  
6.  out((N - 1) / 2 + 1 : (N - 1) / 2 + m, (N - 1) / 2 + 1 : (N - 1) / 2 + n) = gray(:, :);  
7.  if N == 3  
8.      neighborIndex =  [     
9.      -1, -1; -1, 0; -1, 1;  
10.      0, -1;   0, 0;   0, 1;  
11.      1, -1;   1, 0;   1, 1 ];  
12. out(1, 2 : n + 1) = gray(1,1 : n);  
13. out(m + 2, 2 : n + 1) = gray(1, 1 : n);  
14. out(2 : m + 1, 1) = gray(1 : m, 1);  
15. out(2 : m + 1, n + 2) = gray(1 : m, n);  
16. out(1, 1) = gray(1, 1);  
17. out(1, n + 2) = gray(1, n);  
18. out(m + 2, 1) = gray(m, 1);  
19. out(m + 2, n + 2) = gray(m, n);  
20. else  
21.     if N == 5  
22.         neighborIndex =  [     
23.          -2, -2; -2, -1; -2, 0; -2, 1; -2, 2;  
24.          -1, -2; -1, -1; -1, 0; -1, 1; -1, 2;  
25.          0, -2; 0, -1; 0, 0; 0, 1; 0, 2;  
26.          1, -2; 1, -1; 1, 0; 1, 1; 1, 2;  
27.          2, -2; 2, -1; 2, 0; 2, 1; 2, 2; ];  
28.     end  
29. end  
30. neighbor = zeros(N * N, 1);  
31. %neighborIndex = zeros(9, 2);  
32. for i = (N - 1) / 2 + 1 : (N - 1) / 2 + m  
33.     for j =  (N - 1) / 2 + 1 : (N - 1) / 2 + n  
34.     for k = 1 : N * N  
35.         neighbor(k)  = out(i + neighborIndex(k, 1),  j + neighborIndex(k, 2));  
36.        % sort(neighbor);  
37.     end     
38.      Answer(i - 1, j - 1) = uint8(mean(neighbor));  
39.     end  
40. end  
41. re = Answer;  end 

邻域均值滤波效果如下:
加入高斯噪声的图像
3×3模板均值滤波
5×5模板均值滤波
从上述三幅图像中我们可以看出利用均值滤波函数进行去除高斯噪声是有一定效果的,图2.2和图2.3 都明显有去噪效果。此外,3×3模板和5×5模板效果也会略有差异,从图中看出5×5模板明显失去了一些图像细节。

中值滤波

图像中值滤波一般用于图像去除椒盐噪声,类似于电视机出现的雪花,手机屏幕出现明显的噪点等等。
下面的函数是我自己写的去除椒盐噪声的中值滤波函数,参数gray为待输入的原始图像,N为可选参数表示模板的维数。

1.  function re = MiddleBoxFilter(gray, N)  
2.  %中值滤波函数, N表示模板的尺寸 N*N  
3.  [m, n] = size(gray);  
4.  Answer = zeros(m, n);  
5.  out = zeros(m + (N - 1), n + (N - 1));  
6.  out((N - 1) / 2 + 1 : (N - 1) / 2 + m, (N - 1) / 2 + 1 : (N - 1) / 2 + n) = gray(:, :);  
7.  if N == 3  
8.      %边界处理  
9.      neighborIndex =  [     
10.     -1, -1; -1, 0; -1, 1;  
11.      0, -1;   0, 0;   0, 1;  
12.      1, -1;   1, 0;   1, 1  ];  
13. out(1, 2 : n + 1) = gray(1,1 : n);  
14. out(m + 2, 2 : n + 1) = gray(1, 1 : n);  
15. out(2 : m + 1, 1) = gray(1 : m, 1);  
16. out(2 : m + 1, n + 2) = gray(1 : m, n);  
17. out(1, 1) = gray(1, 1);  
18. out(1, n + 2) = gray(1, n);  
19. out(m + 2, 1) = gray(m, 1);  
20. out(m + 2, n + 2) = gray(m, n);  
21. else  
22.     if N == 5  
23.         neighborIndex =  [     
24.          -2, -2; -2, -1; -2, 0; -2, 1; -2, 2;  
25.          -1, -2; -1, -1; -1, 0; -1, 1; -1, 2;  
26.          0, -2; 0, -1; 0, 0; 0, 1; 0, 2;  
27.          1, -2; 1, -1; 1, 0; 1, 1; 1, 2;  
28.          2, -2; 2, -1; 2, 0; 2, 1; 2, 2; ];  
29.     end  
30. end  
31. neighbor = zeros(N * N, 1);  
32. %neighborIndex = zeros(9, 2);  
33. for i = (N - 1) / 2 + 1 : (N - 1) / 2 + m  
34.     for j =  (N - 1) / 2 + 1 : (N - 1) / 2 + n  
35.     for k = 1 : N * N  
36.         neighbor(k)  = out(i + neighborIndex(k, 1),  j + neighborIndex(k, 2));    
37.     end     
38.      Answer(i - 1, j - 1) = median(neighbor);  
39.     end  
40. end  
41. re = Answer;  
42. re = uint8(re); 

中值滤波效果如下
加入椒盐噪声的图像
中值滤波结果
通过两幅图片的对比我们可以看出来,中值滤波对于椒盐噪声的去除效果非常明显,经过中值滤波操作,像基本不会有椒盐噪声的存在。

图像锐化利用一阶水平竖直梯度

水平数值一阶梯度上面原理已经说过,下面列出我写的代码,输入参数pic为待处理的图像矩阵,varargin表示可选的阈值参数,得到的梯度值大于这个值则梯度值赋为255,否则赋为0,这样可以更清楚地对图像的边缘进行表现。

1.  function re = Gradient(pic, varargin)  
2.  %计算水平数值梯度  
3.  [m, n] = size(pic);  
4.  if numel(varargin) == 0  
5.      threshold = 40;  
6.  else  
7.      threshold = varargin{1};  
8.  end  
9.  grad = zeros(m - 1, n - 1);  
10. for i = 1 : m - 1  
11.     for j = 1 : n -1  
12.         G = abs(pic(i + 1, j) - pic(i,  j)) + abs(pic(i, j) - pic(i, j + 1));  
13.         if G > threshold;  
14.             grad(i, j) = 255;  
15.         else  
16.             grad(i, j) = 0;  
17.         end  
18.     end  
19. end  
20. re = grad;  
21. end  

用一阶的水平竖直梯度进行图像锐化的效果如下图所示
原始灰度图
水平竖直梯度
从图中的对比中我们发现图像中具有明显边界的区域很明显地被识别出来,即在图像中轮廓的边界往往蕴含着梯度信息。

图像锐化之Robert梯度

Robert梯度旨在寻找不仅仅水平和竖直两个方向的梯度,所以Robert表示为斜向的梯度。下面是我写的关于Robert梯度计算的函数RobertGradient(pic, varargin),同样的,pic表示原始图像矩阵,varargin表示可选的阈值参数,默认为40,当得到的梯度模大于40,对应的像素点的值为255,否则为0。

1.  function re = RobertGradient(pic, varargin)  
2.  %计算Robert梯度  
3.  [m, n] = size(pic);  
4.  if size(varargin) == 0  
5.      threshold = 40;  
6.  else  
7.      threshold = varargin{1};  
8.  end  
9.  grad = pic(1 : m -1, 1 : n - 1);  
10. for i = 1 : m - 1  
11.     for j = 1 : n - 1   
12.         G = abs(pic(i + 1, j + 1) - pic(i, j)) + abs(pic(i + 1, j) - pic(i, j + 1));  
13.         if G > threshold  
14.             grad(i, j) = 255;  
15.         else  
16.             grad(i, j) = 0;  
17.         end  
18.     end  
19. end  
20. re = grad;  
21. end  

Robert梯度的算法效果如下:
Robert梯度

图像锐化之Sobel梯度

Sobel算子根据像素点上下、左右邻点灰度加权差,在边缘处达到极值这一现象检测边缘。对噪声具有平滑作用,提供较为精确的边缘方向信息,边缘定位精度不够高。当对精度要求不是很高时,是一种较为常用的边缘检测方法。下面是我写的Sobel算子进行锐化的函数。

1.  function re = Sobel_da(gray, varargin)  
2.  [m, n] = size(gray);  
3.  if numel(varargin) == 0  
4.      threshold = 20;  
5.  else  
6.      threshold = varargin{1};  
7.  end  
8.  gray = double(gray);  
9.  S = gray;  
10. NIndex = [  
11.     -1, -1; -1, 0; -1, 1;  
12.     0, -1; 0, 0; 0, 1;  
13.     1, -1; 1, 0; 1, 1];  
14. NxMask = [1, 0, -1, 2, 0, -2, 1, 0, -1];  
15. NyMask = [-1, -2, -1, 0, 0, 0, 1, 2, 1];  
16. for i = 2 : m - 1  
17.     for j = 2 : n - 1  
18.         Sx = 0;  
19.         Sy = 0;  
20.         for k = 1 : 9  
21.             Sx = NxMask(k) * gray(i + NIndex(k, 1), j + NIndex(k, 2)) + Sx; 
22.             Sy = NyMask(k) * gray(i + NIndex(k, 1), j + NIndex(k, 2)) + Sy; 
23.         end  
24.         G = sqrt(double(Sx^2 + Sy^2));  
25.         S(i, j) = uint8(G);  
26.     end  
27. end  
28. re = S;  
29. end  
30. function re = Sobel_da(gray, varargin)  
31. [m, n] = size(gray);  
32. if numel(varargin) == 0  
33.     threshold = 20;  
34. else  
35.     threshold = varargin{1};  
36. end  
37. gray = double(gray);  
38. S = gray;  
39. NIndex = [  
40.     -1, -1; -1, 0; -1, 1;  
41.     0, -1; 0, 0; 0, 1;  
42.     1, -1; 1, 0; 1, 1];  
43. NxMask = [1, 0, -1, 2, 0, -2, 1, 0, -1];  
44. NyMask = [-1, -2, -1, 0, 0, 0, 1, 2, 1];  
45. for i = 2 : m - 1  
46.     for j = 2 : n - 1  
47.         Sx = 0;  
48.         Sy = 0;  
49.         for k = 1 : 9  
50.             Sx = NxMask(k) * gray(i + NIndex(k, 1), j + NIndex(k, 2)) + Sx;  
51.             Sy = NyMask(k) * gray(i + NIndex(k, 1), j + NIndex(k, 2)) + Sy;
52.         end  
53.         G = sqrt(double(Sx^2 + Sy^2));  
54.         S(i, j) = uint8(G);  
55.     end  
56. end  
57. re = S;  
58. end 

Sobel的算法效果如下图所示:
Sobel算子
Sobel边缘检测的结果能用线条较好地描述图像,源图像中对比度较高的区域在结果图中体现为高灰度的像素,简单地说就是对源图像进行了“描边”操作。在进行图像细化时,最直接的思路就是对线条宽度进行缩减,然而图像中的线条往往是不规则的,这使得直接对线条进行处理的方法很难实现。
从上面三种图像锐化操作的结果对比中我们可以看出水平竖直梯度和Robert梯度的锐化效果近似,但对比于Sobel算子,他们在细节上有很明显的缺失。Sobel感觉上图像的边缘信息不仅仅被提取出来而且进一步被细化。

图像傅立叶变换

图像的傅立叶变换
图像的傅立叶变换,原始图像由N行N列构成,N必须是基2的,把这个N*N个包含图像的点称为实部,另外还需要N*N个点称为虚部,因为FFT是基于复数的。计算图像傅立叶变换的过程很简单:首先对每一行做一维FFT,然后对每一列做一维FFT。具体来说,先对第0行的N个点做FFT(实部有值,虚部为0),将FFT输出的实部放回原来第0行的实部,FFT输出的虚部放回第0行的虚部,这样计算完全部行之后,图像的实部和虚部包含的是中间数据,然后用相同的办法进行列方向上的FFT变换,这样N*N的图像经过FFT得到一个N*N的频谱。
实验内容分析
下面我们进行在频域进行图像压缩的实验,实验方法是将图像进行二维的傅立叶变换,然后在频域矩阵中去掉某些频率的值这样使得原本的图像失去一些对应频率的基向量,从而达到数据压缩的效果。这里我们在频域中利用图像像素乘法的思想对频域图像进行压缩,让ROI区域以外的幅度值变为0,ROI部分的区域值不变,最后进行反傅立叶变化,实现图像的压缩变换。这种压缩变换往往是不可逆的,所以并不属于无失真压缩方法。本实验中ROI的选取是以图像中心做为圆心,以指定数值为半径做的圆作为实验中的ROI。
实验的代码

1.  I_in=imread('./cat.jpg');  
2.  I_t=im2double(I_in);  
3.  I_fft=fft2(I_t);  
4.  figure,imshow(abs(I_fft));  
5.    
6.    
7.  I_fft=fftshift(I_fft);  
8.  figure,imshow(abs(I_fft));  
9.  figure,imshow(real(I_fft));  
10. figure,imshow(imag(I_fft));  
11. n=250;%圆形半径  
12. [height,width,t] = size(I_fft);  
13. a=zeros(height,width,t);   
14.     for i=1:height  
15.       for j=1:width  
16.           for k=1:3  
17.           if sqrt((i-height/2)^2+(j-width/2)^2)<=n % 如果在圆内  
18.             a(i,j,k)=1;   
19.           end  
20.           end  
21.     end  
22. end  
23. I_filter1=a;  
24. figure;imshow(abs(I_filter1));  
25.   
26. I_filter=I_filter1.*I_fft;  %获取某圆内的频率分量  
27. figure;imshow(abs(I_filter));  %  
28.   
29. I_o=ifftshift(I_filter);  
30. I_o=ifft2(I_o);  
31. figure('name','最终'),imshow(abs(I_o)); 

实验结果
首先原始图像是一只可爱的猫咪如下图:
cat
将上述图像进行傅立叶变换得到
Fourier变换实部
Fourier变换虚部

上图分别表示傅立叶变换后图像的实部和傅立叶变换后的虚部。
下面将通过选取不同的ROI区域对原图像进行压缩变换,最后对比结果。
压缩变换的区域
压缩变换后的图像
从结果图中可以看出图像经过这样的压缩变换后仅仅能大致看清事物的轮廓,其他细节部分已经完全分辨不出。通过上述实验我们实践了通过傅立叶变换进行压缩变换的原理,并且认识到傅立叶变换在信号及图像处理中具有重要意义。

总结

已经做过四次实验报告了,看着自己写的一些完成一些图像处理的基础算法还是有一点小成就感呢,似乎有点小激动,不过最大的感慨就是,知识和实践相结合是一件非常美妙的事,实践能帮助增强对知识本身的理解,知识推动者实践的热情。无论学习也好,科研也罢,尽自己最大的努力去争取,去探索,我相信我们成功的概率一定会越来越大的。

2018-07-23 21:08:29 weixin_39569242 阅读数 4035

一、实验目的

(1)了解图像增强的目的及意义,加深对图像增强的

     感性认识,巩固所学的图像增强的理论知识和相

     关算法。

(2)熟练掌握直方图均衡化和直方图规定化的计算过

     程。

(3)熟练掌握空域滤波中常用的平滑和锐化滤波器。

(4)熟练掌握低通和高通滤波器的使用方法,明确不

     同性质的滤波器对图像的影响和作用。

(5)掌握最简单的伪彩色变换方法。

二、实验内容

(1)任意选择几幅图像,对其进行平滑处理,用

     不同的平滑模板,对结果进行分析。

(2)任意选择几幅图像,对其进行中值滤波,用

     不同的滤波模板对结果进行分析。

(3)任意选择几幅图像,对其进行梯度锐化,选

     择不同的阈值参数,观察图像有何变化。

(4)对图像进行伪彩色变换,比较彩色增强后的

     图像与原图像有何不同。

三、实验代码及结果、分析

(1)平滑滤波----邻域平均法

  • 代码:

I=imread('E:\JZ数字图像处理\实验3\shiyansan.jpg');

I=rgb2gray(I);

J=imnoise(I,'salt & pepper',0.03);       %加均值为0,方差为0.03的椒盐噪声G=imnoise(I,'gaussian',0.03);             %加均值为0,方差为0.03的高斯噪声。

A=fspecial('average',[4,4]); %4x4均值滤波

A1=fspecial('average',[6,6]); %6x6均值滤波

A2=fspecial('average',[8,8]); %8x8均值滤波

A3=fspecial('average',[10,10]); %10x10均值滤波

J1= imfilter(J,A);  %imfilter均值滤波函数

J2= imfilter(J,A1);  

J3= imfilter(J,A2);  

J4= imfilter(J,A3);   

G1= imfilter(G,A);  

G2= imfilter(G,A1);  

G3= imfilter(G,A2);  

G4= imfilter(G,A3);  

figure(1)

subplot(2,3,1),imshow(I),title('原图像');

subplot(2,3,4),imshow(I),title('原图像');

subplot(2,3,2),imshow(J),title('加入椒盐噪声');         %显示有椒盐噪声图像

subplot(2,3,5),imshow(G),title('加入高斯噪声');         %显示有高斯噪声图像

subplot(2,3,3),imshow(J1),title('4×4均值滤波');%显示有椒盐噪声图像的滤波

subplot(2,3,6),imshow(G1),title('4×4均值滤波'); %显示有高斯噪声图像的滤波

figure(2);

subplot(2,3,1),imshow(I),title('原图像');

subplot(2,3,2),imshow(J),title('加入椒盐噪声');  %显示有椒盐噪声图像

subplot(2,3,4),imshow(J1),title('4×4');

subplot(2,3,3),imshow(J2),title('6×6');

subplot(2,3,5),imshow(J3),title('8×8');

subplot(2,3,6),imshow(J4),title('10×10');

figure(3)

subplot(2,3,1),imshow(I),title('原图像');

subplot(2,3,2),imshow(G),title('加入高斯噪声');       %显示有高斯噪声图像

subplot(2,3,4),imshow(G1),title('4×4');

subplot(2,3,3),imshow(G2),title('6×6');

subplot(2,3,5),imshow(G3),title('8×8');

subplot(2,3,6),imshow(G4),title('10×10');

  • 结果:

  • 分析:
  1. 平滑滤波的作用是对高频分量进行削弱或消除,增强图像的低频分量
  2. 领域平均法是线性运算、中值滤波是非线性运算
  3. 在该题中加的椒盐噪声、高斯噪声方差均为0.03,可以由图1-2、图1-3看出原图像在加入两种噪声之后,画质都有所损失。但是椒盐噪声只是在原图像的某些地方存在,而高斯噪声几乎是分布在原图像的每个地方,造成了图像比之前模糊
  4. 椒盐噪声:又称脉冲噪声,它随机改变一些像素值,是由图像传感器,传输信道,解码处理等产生的黑白相间的亮暗点噪声。椒盐噪声往往由图像切割引起;
  5. 高斯分布:也称正态分布,有均值和方差两个参数,均值反应了对称轴的方位,方差表示了正态分布曲线的胖瘦。高斯分布是最普通的噪声分布
  6. 平滑滤波的模板需要选择合适,模板选择太小,噪声无法滤除,模板选择太大,造成了过度处理,图像会变得逐渐模糊,不管是椒盐噪声、还是高斯噪声如图1-4、1-5模板逐渐增大,图像逐渐模糊。

(2)平滑滤波----中值法

  • 代码:

I=imread('E:\大三课件\大三下\数字图像处理\实验\实验3\shiyansan.jpg');

I=rgb2gray(I);

temp = I;

I = double(I);

%Robert梯度

w1 = [-1 0; 0 1]; 

w2 = [0 -1;1 0];

G1 = imfilter(I,w1,'corr','replicate');

G2 = imfilter(I,w2,'corr','replicate');

G = abs(G1)+abs(G2);

figure(1);

subplot(2,2,1);imshow(temp),title('原图像');

subplot(2,2,2);imshow(abs(G1),[]),title('w1滤波');

subplot(2,2,3);imshow(abs(G2),[]),title('w2滤波');

subplot(2,2,4);imshow(G,[]),title('Robert梯度');

%sobel梯度

w11 = fspecial('sobel');

w21 = w11';

G11 = imfilter(I,w11);

G21 = imfilter(I,w21);

G1 = abs(G11)+abs(G21);

figure(2);

subplot(2,2,1);imshow(temp),title('原图像');

subplot(2,2,2);imshow(G11,[]),title('水平sobel');  

subplot(2,2,3);imshow(G21,[]),title('竖直sobel');

subplot(2,2,4);imshow(G1,[]),title('sobel');

%拉普拉斯滤波

w12 = [0 -1 0;-1 4 -1;0 -1 0];

L12 = imfilter(I,w12,'corr','replicate');

figure(3);

subplot(1,2,1);imshow(temp),title('原图像');

subplot(1,2,2);imshow(abs(L12),[]);

  • 结果:

 

  • 分析:
  1. 锐化滤波器可以消除或减弱图像的低频分量从而增强图像中物体的边缘轮廓信息,使得除边缘以外的像素点的灰度值趋向于零
  2. 我们要将图像类型从uint8转换为double.因为锐化模板计算时常常使输出产生负值,如果采用无符号的 uint8 型, 则负位会被截断
  3. 图3-2Robert梯度分别以w1和w2为模板,w1对接近正45 度边缘有较强响应: w2对接近负45 度边缘有较强响应,对原图像进行滤波就可得到GI和G2,最终的Robert交叉梯度图像(b)为:G = |G1| + |G2|.  
  4. 图3-3,第一个计算了一幅图像的竖直梯度,第2个计算了一幅图像的水平梯度 它们的和可以作为完整的Sobel梯度
  5. 拉普拉斯滤波输出图像中的双边缘,拉普拉斯锐化对一些离散点有较强的响应
  6.  
  7. (3)伪彩色变换
  • 代码:

im=imread('E:\JZ数字图像处理\实验3\shiyansan.jpg');

gray=rgb2gray(im);

I=double(gray);

[m,n]=size(I);

L=256;

for i=1:m

    for j=1:n

if I(i,j)<=L/2    %绿色通道

    R(i,j)=0;

    G(i,j)=2*I(i,j);

    B(i,j)=L;

else if I(i,j)<=L/4   %蓝色通道

        R(i,j)=0;

        G(i,j)=L;

        B(i,j)=-4*I(i,j)+2*L;

    else if I(i,j)<=3*L/2   %红色通道

            R(i,j)=2*I(i,j)-2*L;

            G(i,j)=L;

            B(i,j)=0;

        else

            R(i,j)=L;

            G(i,j)=-4*I(i,j)+4*L;

            B(i,j)=0;

        end

    end

end

    end

end

for i=1:m

    for j=1:n

        rgbim(i,j,1)=R(i,j);

        rgbim(i,j,2)=G(i,j);

        rgbim(i,j,3)=B(i,j);

    end

end

rgbim=rgbim/256;

figure(1);

subplot(1,2,1),imshow(gray),title('原图');

subplot(1,2,2),imshow(rgbim),title('伪彩色变换');

  • 结果:

  • 分析:
  1. 伪彩色处理是指通过将每一个灰度级匹配到彩色空间上的一点,将单色图像映射为一副彩色图像的彩色图像
  2. 最小的灰度值0映射为红色,中间的灰度值L/2映射为绿色,最高的灰度值L映射为蓝色
2017-05-31 10:30:10 hongbin_xu 阅读数 7927

以下这些实验中的代码全部是我自己编写调试通过的,到此,最后进行一下汇总。

数字图像处理实验(1):PROJECT 02-01, Image Printing Program Based on Halftoning (基于半色调技术的图像打印技术
链接:http://blog.csdn.net/hongbin_xu/article/details/70340458

数字图像处理实验(2):PROJECT 02-02, Reducing the Number of Gray Levels in an Image(减小图像的灰度级
链接:http://blog.csdn.net/hongbin_xu/article/details/70344579

数字图像处理实验(3):PROJECT 02-03, Zooming and Shrinking Images by Pixel Replication(通过像素复制缩放图片
链接:http://blog.csdn.net/hongbin_xu/article/details/70344600

数字图像处理实验(4):PROJECT 02-04 [Multiple Uses],Zooming and Shrinking Images by Bilinear Interpolation(通过双线性内插法缩放图片
链接:http://blog.csdn.net/hongbin_xu/article/details/70460542

数字图像处理实验(5):Proj03-01 ~ Proj03-06(这里有6个实验
链接:http://blog.csdn.net/hongbin_xu/article/details/70994340

数字图像处理实验(5):PROJECT 04-01 [Multiple Uses],Two-Dimensional Fast Fourier Transform(二维快速傅里叶变换
链接:http://blog.csdn.net/hongbin_xu/article/details/71366094

数字图像处理实验(6):PROJECT 04-02,Fourier Spectrum and Average Value(傅立叶频谱和平均值
链接:http://blog.csdn.net/hongbin_xu/article/details/71375543

数字图像处理实验(7):PROJECT 04-03 , Lowpass Filtering(低通滤波
链接:http://blog.csdn.net/hongbin_xu/article/details/72722264

数字图像处理实验(8):PROJECT 04-04,Highpass Filtering Using a Lowpass Image(使用低通图像进行高通滤波
链接:http://blog.csdn.net/hongbin_xu/article/details/72722475

数字图像处理实验(9):PROJECT 04-05,Correlation in the Frequency Domain(频域的相关性
链接:http://blog.csdn.net/hongbin_xu/article/details/72722842

数字图像处理实验(10):PROJECT 05-01 [Multiple Uses],Noise Generators(噪音发生器
链接:http://blog.csdn.net/hongbin_xu/article/details/72775845

数字图像处理实验(11):PROJECT 05-02,Noise Reduction Using a Median Filter(使用中值滤波器减少噪声
链接:http://blog.csdn.net/hongbin_xu/article/details/72775908

数字图像处理实验(12):PROJECT 05-03,Periodic Noise Reduction Using a Notch Filter(使用陷波滤波器减少周期噪声
链接:http://blog.csdn.net/hongbin_xu/article/details/72779188

数字图像处理实验(13):PROJECT 05-04,Parametric Wiener Filter(参数维纳滤波器
链接:http://blog.csdn.net/hongbin_xu/article/details/72779509

数字图像处理实验(14):PROJECT 06-01,Web-Safe Colors(网页安全色
链接:http://blog.csdn.net/hongbin_xu/article/details/72784278

数字图像处理实验(15):PROJECT 06-02,Pseudo-Color Image Processing(伪彩色图像处理
链接:http://blog.csdn.net/hongbin_xu/article/details/72784321

数字图像处理实验(16):PROJECT 06-03,Color Image Enhancement by Histogram Processing(通过直方图处理的彩色图像增强
链接:http://blog.csdn.net/hongbin_xu/article/details/72784377

数字图像处理实验(17):PROJECT 06-04,Color Image Segmentation(彩色图像分割
链接:http://blog.csdn.net/hongbin_xu/article/details/72784421

2018-04-24 20:36:38 weixin_39569242 阅读数 2861

 

一、实验目的

1)了解图像变换的意义和手段。

2)熟悉傅立叶变换的基本性质。

3)通过实验了解二维频谱的分布特点。

4)了解余弦变换或WalshHadamard变换

二、实验内容

    任意选择几幅图像,对其进行傅立叶变换,观察图像的频谱特征。

   1)对图像进行傅里叶变换(包括移位和未移位),观察频谱信号

2)观察频谱的三维图形

① 移位前:

v 代码:

f=zeros(64,64);

f(15:50,15:50)=1;%输入64*64的黑色图像矩阵

F=fft2(f);%平移前傅立叶变换

F2=fftshift(abs(F));%频谱中心化

x=1:64;y=1:64;

%绘制平移前

figure(1)

subplot(1,3,1),imshow(f);title('平移前原始图像');

subplot(1,3,2),imshow(abs(F));title('图像傅里叶变换');

subplot(1,3,3),imshow(abs(F2));title('图像中心化傅里叶变换');

figure(2)

subplot(2,1,1),title('三维图像');

mesh(abs(real(F)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,F2(x,y),'FaceColor','red');

结果:


                        图1

                 

2:彩色三维图像                          图3FFT变换

② 移位后:

代码:

ff=circshift(f,[-10 -5]);%将原始图像向左向上移,及X轴,Y轴负方向

FF=fft2(ff);%平移后傅立叶变换

FF2=fftshift(abs(FF));%频谱中心化

%绘制平移后

figure(3)

subplot(2,3,1),imshow(ff);title('平移后原始图像');

subplot(2,3,2),imshow(abs(FF));title('图像傅里叶变换');

subplot(2,3,3),imshow(abs(FF2));title('图像中心化傅里叶变换');

figure(4)

subplot(2,1,1),title('三维图像');

mesh(abs(real(FF)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,FF2(x,y),'FaceColor','red');

结果:

                      图4

 

                        

  图5:彩色三维图像                            图6FFT变换

   3)进行傅里叶逆变换,重建图像

   4)根据幅值谱重建图像信息

   5)根据相位谱重建图像信息

代码:

%平移之前的频谱、相位谱、频谱逆变换、相谱逆变换、傅里叶逆变换

xf1=log(abs(F));%频谱

xf2=angle(F);%相位谱

xf12=ifft2(xf1);

xf13=ifftshift(abs(xf12));%频谱逆变换

xf22=ifft2(xf2);

xf23=ifftshift(abs(xf22));%相位逆变换

xr1=ixff1=log(abs(FF));%平移后频谱

%平移之后频谱、相位谱

xff2=angle(FF);%平移后相位谱

fft2(F);%进行傅里叶逆变换

%平移前平移后相位谱、频谱对比

figure(5)

subplot(2,2,1),imshow(abs(xf2));title('傅里叶相位谱');

subplot(2,2,2),imshow(abs(xff2));title('平移后傅里叶相位谱');

subplot(2,2,3),imshow(abs(xf1));title('傅里叶幅度谱');

subplot(2,2,4),imshow(abs(xff1));title('平移后傅里叶幅度谱');

%平移前逆变换、相位谱逆变换、幅度谱逆变换

figure(6)

subplot(2,2,1),imshow(abs(xf23));title('傅里叶相位谱逆变换');

subplot(2,2,2),imshow(abs(xf13));title('傅里叶幅度谱逆变换');

subplot(2,2,3);imshow(abs(xr1));title('逆变换');

结果:

 

7:平移前后相位谱、幅度谱

 

8:平移前逆变换、相位谱逆变换、幅度谱逆变换

 

三、实验分析:

1、图像平移之后的傅里叶幅度谱不会发生变化,而仅仅是相位谱产生了一定的相移特性,如图7

2、图像进行傅里叶变换后,进行傅里叶逆变换会还原为原来的图像,如图8

3、相位谱逆变换之后会大致还原原图像的轮廓、幅度谱逆变换之后还原不出原图像,只会确定图像的精度。如图8

附代码:

f=zeros(64,64);

f(15:50,15:50)=1;%输入64*64的黑色图像矩阵

F=fft2(f);%平移前傅立叶变换

F2=fftshift(abs(F));%频谱中心化

x=1:64;

y=1:64;

 

%平移之前的频谱、相位谱、频谱逆变换、相谱逆变换、傅里叶逆变换

xf1=log(abs(F));%频谱

xf2=angle(F);%相位谱

xf12=ifft2(xf1);

xf13=ifftshift(abs(xf12));%频谱逆变换

xf22=ifft2(xf2);

xf23=ifftshift(abs(xf22));%相位逆变换

xr1=ifft2(F);%进行傅里叶逆变换

 

%平移之后的傅里叶变换频谱、相位谱

ff=circshift(f,[-10 -5]);%将原始图像向左向上移,及X轴,Y轴负方向

FF=fft2(ff);%平移后傅立叶变换

FF2=fftshift(abs(FF));%频谱中心化

xff1=log(abs(FF));%平移后频谱

xff2=angle(FF);%平移后相位谱

 

%绘制平移前

figure(1)

subplot(1,3,1),imshow(f);title('平移前原始图像');

subplot(1,3,2),imshow(abs(F));title('图像傅里叶变换');

subplot(1,3,3),imshow(abs(F2));title('图像中心化傅里叶变换');

 

figure(2)

subplot(2,1,1),title('三维图像');

mesh(abs(real(F)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,F2(x,y),'FaceColor','red');

 

%绘制平移后

figure(3)

subplot(2,3,1),imshow(ff);title('平移后原始图像');

subplot(2,3,2),imshow(abs(FF));title('图像傅里叶变换');

subplot(2,3,3),imshow(abs(FF2));title('图像中心化傅里叶变换');

 

figure(4)

subplot(2,1,1),title('三维图像');

mesh(abs(real(FF)),'FaceColor','white');

subplot(2,1,2),title('FFT变换');

mesh(x,y,FF2(x,y),'FaceColor','red');

 

%平移前平移后相位谱、频谱对比

figure(5)

subplot(2,2,1),imshow(abs(xf2));title('傅里叶相位谱');

subplot(2,2,2),imshow(abs(xff2));title('平移后傅里叶相位谱');

subplot(2,2,3),imshow(abs(xf1));title('傅里叶幅度谱');

subplot(2,2,4),imshow(abs(xff1));title('平移后傅里叶幅度谱');

 

%平移前逆变换、相位谱逆变换、幅度谱逆变换

figure(6)

subplot(2,2,1),imshow(abs(xf23));title('傅里叶相位谱逆变换');

subplot(2,2,2),imshow(abs(xf13));title('傅里叶幅度谱逆变换');

subplot(2,2,3);imshow(abs(xr1));title('逆变换');