《数字图像处理：Matlab算法设计与解译》 课件 第5、6章 频率域图像增强、图像恢复

第五章

频率域图像处理

频率域图像是把空间域图像像素的灰度值表示成随位置变化的空间频率，并以频谱（也称为频谱图））的形式表示图像信息分布特征的一种表示方式。

频率域图像处理是指在图像的频率域中对图像进行某种处理的方法，这种方法以傅里叶变换为基础，也即先通过傅里叶变换把图像从空间域变换到频率域，然后用频率域方法对图像进行处理，处理完后再利用傅里叶反变换把图像变换回空间域。5.1二维离散傅里叶变换5.1二维离散傅里叶变换

由于离散傅里叶变换描述了离散信号的时域及空间域表示与频域表示的关系，所以利用基于离散傅里叶变换的时域与频域分析方法可解决大多数图像处理问题，因而离散傅里叶变换在图像处理领域获得了极为广泛的应用。

由于二维离散傅里叶变换对应地可以描述成一个二维函数，所以下面介绍应用于图像处理的。1、二维离散傅里叶变换的定义

5.1.1二维离散傅里叶变换的定义和傅里叶频谱

设f(x,y)是在空间域上等间隔采样得到的M×N的二维离散信号，x和y是离散实变量，u和v为离散频率变量，则二维离散傅里叶变换对一般地定义为：

(u=0,1,…,M-1；v=0,1,…,N-1)(5.1)

(x=0,1,…,M-1；y=0,1,…,N-1)

(5.2)

1、二维离散傅里叶变换的定义

根据:

(5.3)

可知，F(u,v)的运算结果是一个复数结果，用R(u,v)表示其实部，用I(u,v)表示其虚部，可将二维离散傅里叶变换的频谱和相位角分别定义为：(5.4)和欧拉公式:2、图像的傅里叶频谱特性及频谱图傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，从物理效果看，傅立叶变换是将图像从空间域转换到频率域。

5.1.1二维离散傅里叶变换的定义和傅里叶频谱

(a)图像(b)图像的原频谱图

性质包括：线性性、可分离性、平均值性质、周期性、共扼对称性、空间位置和空间频率的平移性、旋转性、尺度变换性、卷积性质等。

本节仅介绍几种比较重要且与书中内容有关的性质。

5.1.2二维离散傅里叶变换的若干重要性质1、变换系数矩阵

5.1.2二维离散傅里叶变换的若干重要性质根据变换式(5.2)，由于u有0,1,…,M-1的M个可能的取值，v有0,1,…,N-1的N个可能的取值，所以f(x,y)由M*N个频率分量组成，所以每个频率分量都与一个特定的(u,v)值相对应；

且对于某个特定的(u,v)值来说，当(x,y)取遍所有可能的值(x=0，1，…，N-1；y=0，1，…，N-1）时，就可得到对应于该特定的(u,v)值的一个变换系数矩阵：(5.2)

5.1.2二维离散傅里叶变换的若干重要性质1、变换系数矩阵

可见，该矩阵的值仅与N有关，与f(x,y)无关。

2、可分离性

式(5.1)和式(5.2)的二维离散傅里叶变换对可写成如下的分离形式：

上述的可分离表示形式说明，可以连续运用两次一维DFT来实现一个二维DFT。

5.1.2二维离散傅里叶变换的若干重要性质5.1.2二维离散傅里叶变换的若干重要性质2、可分离性

2、可分离性

行变换列变换5.1.2二维离散傅里叶变换的若干重要性质3、平均值

一幅图像的灰度平均值可表示为：

(5.10)5.1.2二维离散傅里叶变换的若干重要性质4、周期性

对于M×N的图像和二维离散傅里叶变换对的一般定义式(5.1)和(5.2)，F(u,v)的周期性定义为：

(m,n=0,±1,±2,…）(5.13)5.1.2二维离散傅里叶变换的若干重要性质5、共轭对称性

设f(x,y)为实函数，则其傅里叶变换F(u,v)具有共轭对称性：(5.14)(5.15)5.1.2二维离散傅里叶变换的若干重要性质6、平移性

对于M×N的图像f(x,y)和二维离散傅里叶变换对的一般定义式(5.1)和(5.2)，若设用符号表示函数与其傅里叶变换的对应性，则傅里叶变换的平移性可表示为：

(5.16)

(5.17)其中，式(5.16)说明，给函数乘以一个指数项，就相当于把其变换后的傅里叶频谱在频率域进行平移。

式(5.17)说明，给傅里叶频谱乘以一个指数项，就相当于把其反变换后得到的函数在空间域进行平移。

5.1.2二维离散傅里叶变换的若干重要性质6、平移性

可见，当空域中f(x,y)产生移动时，在频域中只发生相移，而傅立叶变换的幅值不变，即：

同理，当频域中F(u,v)产生移动时，相应的f(x,y)在空域中也只发生相移，而幅值不变。5.1.2二维离散傅里叶变换的若干重要性质1、图像傅里叶频谱关于(M/2，N/2)的对称性

5.1.3图像的傅里叶频谱特性分析设f(x,y)是一幅大小为M×N的图像，根据离散傅立叶变换的周期性公式(5.13)，有：根据(5.22)，对于u=0，M-U=M当u=0、v=0时：当u=0、v=1时：当u=0、v=2时：┆┆当u=0、v=N/2时：0N/2NMM/2(M,N)(M/2,N/2)ABCDvu(M/2,N)(M,N/2)1、图像傅里叶频谱关于(M/2，N/2)的对称性

(5.19)5.1.3图像的傅里叶频谱特性分析1、图像傅里叶频谱关于(M/2，N/2)的对称性

(5.19)同理，对于v=0，N–v=N：当u=0、v=0时：当u=1、v=0时：当u=2、v=0时：┆┆当u=M/2、v=0时：由此可得：

频谱图A区与D区和B区与C区关于坐标(M/2,N/2)对称。5.1.3图像的傅里叶频谱特性分析1、图像傅里叶频谱关于(M/2，N/2)的对称性

下图是原点坐标位于(0,0)的图像的傅里叶变换频谱关于(M/2,N/2)对称的两个例子。

关于(M/2,N/2)对称示例1/示例2

(a)图像(b)图像的原频谱图(a)图像(b)图像的原频谱图5.1.3图像的傅里叶频谱特性分析2、原点坐标平移到(M/2，N/2)后的傅里叶频谱

(0,0）(M/2,N/2）vuvu0NM(M,N）yx0NM(M,N）vu5.1.3图像的傅里叶频谱特性分析2、原点坐标平移到(M/2，N/2)后的傅里叶频谱

即：基于原频谱图，只要把其左上角的A区向下和向右平移到结果频谱图的D区,对应地把其右下角的D区向上和向左平移到结果频谱图的A区，同理把其右上角的B区向下和向左平移到结果频谱图的C区,对应地把其左下角的C区向上和向右平移到结果频谱图的B区，就可以得到中心对称的频谱图。

5.1.3图像的傅里叶频谱特性分析2、原点坐标平移到(M/2，N/2)后的傅里叶频谱

坐标位于(0,0)的频谱图坐标位于(0,0)的频谱图图5.7原点平移到(M/2，N/2)后的频谱图

原点在(0,0)时的频谱图:5.1.3图像的傅里叶频谱特性分析(a)(b)2、原点坐标平移到(M/2，N/2)后的傅里叶频谱

（a）原图像（b）移动前的幅度谱（c）移动后幅度谱5.1.3图像的傅里叶频谱特性分析坐标位于(0,0)的频谱图坐标位于(0,0)的频谱图图5.7原点平移到(M/2，N/2)后的频谱图

原点在(0,0)时的频谱图:的傅里叶频谱图。5.1.3图像的傅里叶频谱特性分析(a)(b)离散傅里叶变换的-频谱示例

离散傅里叶变换的-频谱示例

离散傅里叶变换的-频谱示例

离散傅里叶变换的-频谱示例

3、对图像进行傅里叶变换的意义

（1）简化计算，在空间域中处理图像时所进行的复杂的卷积运算，等同于在频率域中简单的乘积运算。5.1.3图像的傅里叶频谱特性分析4、对图像进行傅里叶变换的意义

（2）由于在用频谱图表示的频率域图像中，中心部位是能量集中的低频特征，反映的是图像的平滑部分；随着不断远离频谱图的中心位置，对应于空间图像中变化越来越快的细节、边缘、结构复杂区域、突变部位和噪声等高频成分逐渐加强。所以，在频率域中滤波的概念更为直观，更容易理解；也即，某些在空间域中难于处理或处理起来比较复杂的问题，在频率域却比较容易处理。5.1.3图像的傅里叶频谱特性分析4、对图像进行傅里叶变换的意义

（3）某些只能在频率域处理的特定应用需求，比如频率域图像特征提取、数据压缩、纹理分析、水印嵌入等。

5.1.3图像的傅里叶频谱特性分析5.1.4快速离散傅里叶变换及其实现5、二维DFT的Matlab编程

相关函数：（1）F=fft2(f1);%二维傅里叶变换

（2）FS=fftshift(F);%将变换的频率图像四角移动到中心

（3）S=log(1+abs(FS));%对频谱进行对数运算，提升中心低频部分，为了更好地显示

（4）fr=real(ifft2(ifftshift(FS)));%二维傅里叶逆变换5.1.4快速离散傅里叶变换及其实现5、二维DFT的Matlab编程

（5）ret=im2uint8(mat2gray(fr))其中：①mat2gray的功能是实现图像矩阵的归一化操作。②im2uint8的功能是把图像数据类型转换为无符号八位整型。

【例】傅里叶变换（Fouriertransform）matlab程序(Fourier_T51.m)clc;clearall;closeall;gray_img=imread('d:\0_matlab图像课编程\girl.jpg')%读灰度图像FT_img=fft2(gray_img);%二维傅里叶正变换FTS_img=fftshift(FT_img);%分4块平移频谱图，实现频谱图的中心对称FTS_log_img=log(1+abs(FTS_img));%频谱图对数运算Inverse_I=ifftshift(FTS_img);%反移频谱中心Inverse_img=real(ifft2(Inverse_I));%傅里叶逆变换，并取变换结果的实部Inverse_FT_img=uint8(Inverse_img);%转换成8位图像subplot(1,3,1);imshow(gray_img);title('原灰度图像');subplot(1,3,2);imshow(FTS_log_img,[]);title('傅立叶频谱图像');subplot(1,3,3);imshow(Inverse_FT_img);title('反变换结果图像');

5.2频率域图像处理的实现思路

二维离散傅里叶变换很好地描述了二维离散信号的空间域与频域之间的关系，所以对于那些在空间域中表述起来比较困难，甚至是不太可能实现的图像处理问题，可以先通过对图像进行离散傅里叶变换把图像变换到频率域，然后利用适当的频率域图像处理方式对图像进行处理，处理完后再把它转换回空间域中，就可解决那些在空间域不便于解决的图像处理问题。由傅立叶频谱的特性可知，u和v同时为0时的频率成分对应于图像的平均灰度级。当从(傅立叶)变换的原点离开时，低频对应着图像的慢变化分量，比如一幅图像中较平坦的区域；当进一步离开原点时，较高的频率开始对应图像中变化越来越快的灰度级，它们反映了一幅图像中物体的边缘和灰度级突发改变(如噪声)部分的图像成分。

频率域图像增强正是基于这种机理，通过对图像的傅立叶频谱进行低通滤波(使低频通过，使高频衰减)来虑除噪声，通过对图像的傅立叶频谱进行高通滤波(使高频通过，使低频衰减)突出图像中的边缘和轮廓。

5.2.1基本实现思想5.2.1基本实现思想5.2.1基本实现思想以上过程可简要地描述为图5.10。5.2.1基本实现思想5.2.1基本实现思想

频率域（高低痛）滤波的原理：转移函数H(u,v)的设计:

一是先凭直观感觉选择一个理想的滤波器模型，然后通过反复的滤波实验和参数修正来逼近并设计出实际的滤波器。

二是利用频率成分和图像外表之间的对应关系选择频率域滤波器。

三是基于数学和统计准则设计频率域滤波器。5.2.2转移函数的设计

对于大小为M×N的函数f(x,y)和h(x,y)，其卷积形式表示为：

用F(u,v)和H(u,v)分别表示f(x,y)和h(x,y)的傅立叶变换，则有傅立叶变换变换对：（5.23）

（5.24）

也即，空间域的卷积在频率域简化为相乘，频率域的卷积在空间域简化为相乘；有时也可以将频率域的滤波器函数变换到空间域，然后再在空间域对图像进行滤波运算。

5.2.2转移函数的设计

5.3基于频率域的图像噪声消除——频率域低通滤波

在频率域中，图像中的噪声和边缘对应于傅立叶频谱的高频部分，选择能使低频通过、使高频衰减的转移函数，就可以实现低通滤波，达到虑除噪声的目的。

5.3.1理想低通滤波器

1.

理想低通滤波器的转移函数定义

其中，D0是1个非负整数，D(u,v)为频率平面从原点到点(u,v)的距离。（5.25）

D(u,v)的值：

设已经将傅里叶频谱的原点平移到(M/2,N/2），则点（u,v）到频率平面原点(M/2,N/2)的距离为：（5.26）

(0,0）(M/2,N/2）vu5.3.1理想低通滤波器

1.

理想低通滤波器的转移函数定义

2.

理想低通滤波器的含义在半径为D0的圆内，所有的频率没有衰减地通过该滤波器；而在此半径的圆之外的所有频率完全被衰减掉。所以称D0为截至频率。5.3.1理想低通滤波器

3.

理想低通滤波器的转移函数横截面图和透视图

1（a）转移函数（b）透视图该透视图的含义是：

只有那些位于该圆柱体内的频率范围的信号才能通过，而位于圆柱体外的频率成分都将被虑除掉。频谱幅度谱5.3.1理想低通滤波器

例5.4频率域理想低通滤波器的滤波效果及低频特性分析若一般地设R为截止频率的圆周半径，EB为圆周内能量（图像功率）与原图像总能量（总功率）的百分比，根据图像信号能量在频率域上的分布有：(5.27)

5.3.1理想低通滤波器

例5.4（续1）(a)原图像(b)频谱图

(c)截止频率半径10

(d)截止频率半径20(e)截止频率半径40(f)截止频率半径80

5.3.1理想低通滤波器

5.3.1理想低通滤波器

1.巴特沃斯低通滤波器的转移函数定义（5.28）

其中，D0为截至频率，D(u,v)为频率平面从原点到点(u,v)的距离，且D(u,v)由下式给出。即：5.3.2巴特沃斯低通滤波器（5.26）

2.转移函数横截面图和透视图（阶数n为1～4）

（a）转移函数（b）透视图透视图的含义是：只有那些位于该草帽型体内的频率范围的信号才能通过，而位于草帽型体外的频率成分都将被虑除掉。频谱幅度谱5.3.2巴特沃斯低通滤波器图5.14是利用巴特沃斯低通滤波器进行图像去噪的实验结果。5.3.2巴特沃斯低通滤波器【例】利用巴特沃斯低通滤波器进行图像去噪matlab编程。clc;clearall;closeall;img0=imread('d:\0_matlab图像课编程\lena.jpg');Noisy_img=imnoise(img0,‘salt&pepper’,0.02);%噪声密度为0.02的椒盐噪声f=double(Noisy_img);FT_img=fft2(f);%傅里叶变换FTS_img=fftshift(FT_img);%平移频谱图为中心对称D0=input('\n请输入非负的截止频率值(10/20/30)D0=')%正整数n=input('\n请输入巴特奥斯滤波器的阶数(1/2/3)n=');%滤波器阶数1至3[h,w]=size(img0);M=fix(h/2);N=fix(w/2);%将h/2和w/2分别向零方向(向下)取整foru=1:hforv=1:wduv=sqrt((u-M)^2+(v-N)^2);huv=1/(1+(duv/D0)^(2*n));%计算巴特奥斯低通滤波器转移函数值

huv_img(u,v)=huv*FTS_img(u,v);%对傅里叶频谱进行滤波

endend【例】利用巴特沃斯低通滤波器进行图像去噪matlab编程。Inverse_I=ifftshift(huv_img);%反移频谱中心Inverse_img=real(ifft2(Inverse_I));%傅里叶逆变换，并取变换结果的实部Inverse_H_img=uint8(Inverse_img);%转换成8位图像figure;subplot(1,3,1);imshow(img0);title('原图像');subplot(1,3,2);imshow(Noisy_img);title('加入椒盐噪声图像');subplot(1,3,3);imshow(Inverse_H_img);title('低通滤波结果图像');截至频率D0=10，滤波器阶数n=2

1.二维高斯低通滤波器的转移函数定义（5.29）

其中，D(u,v)为频率平面从原点到点(u,v)的距离，σ

表示高斯曲线扩展的程度，且D(u,v)由下式给出：当σ

=D0时，可得到高斯低通滤波器的一种更为标准的表示形式：（5.30）

5.3.3高斯低通滤波器（5.26）

2.转移函数横截面图和透视图（D=10,20,40,100）

透视图的含义是：只有那些位于该草帽型体内的频率范围的信号才能通过，而位于草帽型体外的频率成分都将被虑除掉。（a）转移函数（b）透视图频谱幅度谱5.3.3高斯低通滤波器5.4基于频率域的图像增强——频率域高通滤波

在频率域中，图像中的边缘和灰度的陡峭变化对应于傅立叶频谱的高频部分，选择能使高频通过、使低频衰减的转移函数，就可以实现高通滤波，达到突出图像的高频边缘成分，实现图像增强的效果。

1.理想高通滤波器的转移函数定义（5.31）

其中，D0为截至频率；D(u,v)为频率平面从原点到点(u,v)的距离，且D(u,v)由下式给出：5.4.1理想高通滤波器

（5.26）

2.理想高通滤波器的含义将以半径为D0的圆周内的所有频率置零，而让圆周外的所有频率毫不衰减地通过。

5.4.1理想高通滤波器

3.转移函数的横截面图和透视图1透视图的含义是：只有那些位于该圆柱体外的频率范围的信号才能通过，而位于圆柱体内的频率成分都将被虑除掉。频谱幅度谱5.4.1理想高通滤波器

1.巴特沃斯高通滤波器的转移函数定义（5.32）

其中，D0为截至频率；D(u,v)为频率平面从原点到点(u,v)的距离，且D(u,v)由下式给出：5.4.2巴特沃斯高通滤波器（5.26）

2.

转移函数的横截面图和透视图透视图的含义是：只有那些位于该倒立型草帽体外的频率范围的信号才能通过,而位于倒立型草帽体内的频率成分都将被虑除掉。频谱幅度谱5.4.2巴特沃斯高通滤波器

1.高斯高通滤波器的转移函数定义（5.33）

其中，D(u,v)为频率平面从原点到点(u,v)的距离，且D(u,v)由下式给出：5.4.3高斯高通滤波器（5.26）

2.转移函数的横截面图和透视图透视图的含义是：只有那些位于该倒立型草帽体外的频率范围的信号才能通过,而位于倒立型草帽体内的频率成分都将被虑除掉。频谱幅度谱5.4.3高斯高通滤波器5.4.3高斯高通滤波器

2.转移函数的横截面图和透视图（a）原图（b）D0=30（c）D0=60由上图可以看出：随着D0值的增大，增强效果更加明显，即使对于微小的物体和细线条，用高斯滤波器滤波后也比较清晰。5.4.3高斯高通滤波器【例】利用高斯高通滤波器进行图像增强matlab编程。clc;clearall;closeall;img0=imread('d:\0_matlab图像课编程\girl.jpg');f=double(img0);FT_img=fft2(f);%傅里叶变换FTS_img=fftshift(FT_img);%平移频谱图为中心对称D0=input('\n请输入非负的截止频率值(10/20/30/40)D0=')%正整数[h,w]=size(img0);M=fix(h/2);N=fix(w/2);%将h/2和w/2分别向零方向(向下)取整d=2*D0^2;foru=1:hforv=1:wduv=sqrt((u-M)^2+(v-N)^2);huv=1-exp(-duv^2/d);%计算高斯高通滤波器转移函数值

huv_img(u,v)=huv*FTS_img(u,v);%对傅里叶频谱进行滤波

endend【例】利用巴特沃斯低通滤波器进行图像去噪matlab编程。Inverse_I=ifftshift(huv_img);%反移频谱中心Inverse_img=real(ifft2(Inverse_I));%傅里叶逆变换，并取变换结果的实部img1=f+Inverse_img;%形成结果图像result_img=uint8(img1);%转换成8位图像figure;subplot(1,2,1);imshow(img0);title('原图像');subplot(1,2,2);imshow(result_img);title('高斯高通滤波结果图像');截至频率D0=305.5带阻滤波和带通滤波

在某些应用中，图像的质量可能受到带有一定规律的结构性噪声的影响。比如，图像上叠加有正弦干扰图案就是这类噪声的一个典型情况。当正弦干扰图案比较明显时，会在图像的频谱平面上出现2个比较明显的对称点（由于傅立叶变换的共轭对称性所致）。象这种用于消除以某点为对称中心的给定区域内的频率，或用于阻止以原点为对称中心的一定频率范围内信号通过的问题，就可以用带阻滤波器实现。

5.5.1带阻滤波器

1.带阻滤波的传递函数定义1一个用于消除以某点为中心，以D为半径的圆域上的带阻滤波器，可以通过把以原点为中心的高通滤波器平移到该点得到，设该带阻滤波器的中心为点(u0,v0)，半径为D0，则其传递函数定义为：

（5.34）

其中：

（5.35）

5.5.1带阻滤波器由于傅立叶变换的共轭对称性，要求带阻滤波器必须成对出现，所以一个用于消除以（u0,v0）为中心，以D0为半径的对称区域内的所有频率的理想带阻滤波器的转移函数定义为：（5.36）

其中：

（5.37）

（5.38）

1.带阻滤波的传递函数定义25.5.1带阻滤波器

2.一种n阶径向对称的巴特沃斯带阻滤波器的传递函数定义（5.39）

其中，W为阻带带宽，

D0为阻带中心半径。

5.5.1带阻滤波器

3.带阻滤波转移函数的透视图透视图的含义是：

只有那些位于两个立方体外的频率范围的信号才能通过，而位于两个立方体内的频率成分都将被虑除掉。5.5.1带阻滤波器

1.带通滤波器转移函数的定义（5.40）

带通滤波器也可以通过对相应的带阻滤波器进行“翻转”获得。若设H’(u,v)为带阻滤波器的传递函数，则对应的带通滤波器的传递函数H(u,v)可定义为：

（5.41）

5.5.2带通滤波器

2.带通滤波器转移函数的透视图透视图的含义是：

只有那些位于两个立方体内的频率范围的信号会被通过，而位于两个立方体外的频率成分都将被虑除掉。5.5.2带通滤波器作业：

谢谢第五章

图像恢复

所谓图像恢复，就是使退化了的图像去除退化因素，并以最大保真度恢复成原来图像的一种技术。

图像恢复与图像增强的研究内容有一定的交叉性。一般认为，图像增强是一种改进图像视觉效果的技术；而图像恢复是一种对退化（或品质下降）了的图像去除退化因素，并进而复原或重建被退化了的图像的技术。

根据以上定义，通过去模糊函数去除图像模糊应属于一种图像恢复技术。6.1图像的退化模型

6.1.1常见退化现象的物理模型1、造成图像退化的因素造成图像退化的原因是多方面的。

下面是一些具体的退化因素的例子：透镜象差/色差聚焦不准（失焦，限制了图象锐度）模糊（限制频谱宽度）噪声（是一个统计过程）抖动（机械、电子）

非线性退化空间模糊退化平移退化叠加随机噪声退化图6.1常见的4种退化现象的物理模型示意图

6.1.1常见退化现象的物理模型2、常见的4种退化现象的物理模型成像模糊恢复后结果运动成像模糊恢复后结果3、图像退化举例16.1.1常见退化现象的物理模型光脉冲退化的光脉冲6.1.1常见退化现象的物理模型3、图像退化举例26.1.1常见退化现象的物理模型3、图像退化举例3原始图像模糊图像抖动模糊图像6.1.1常见退化现象的物理模型3、图像退化举例4

设f(x,y)是一幅原图像，图像的退化过程可以理解为一个作用于原图像f(x,y)的系统H，或理解为施加于原图像f(x,y)上的一个运算H；同时数字图像也常会因受一些随机误差，也即噪声n(x,y)而退化。由此可得退化模型:并可以表示为：

（6.1）

6.1.2图像退化模型的表示6.2逆滤波图像恢复

图像恢复按是否对图像恢复施加约束条件分为：

无约束恢复方法；

有约束恢复方法。

逆滤波图像恢复方法是一种典型的无约束最小二乘方恢复方法；维纳滤波是一种典型的有约束图像恢复方法。

接下来先介绍逆滤波图像恢复方法，然后介绍维纳滤波图像恢复方法。为了分析噪声的统计特征，以及噪声与图像的相关情况，需要引入最小二乘方恢复。由式（6.1）有：6.2.1无约束最小二乘方恢复（6.2）

当叠加噪声n无法知道时，显然可从g–Hf

获得n。

由于g是已知的退化图像，所以如果取为f的估计，就可使在最小均方误差的意义下代替Hf，并通过求退化后的实际图像g与退化图像的估值

的模(或范数)平方，也即从而可把图像的恢复问题看作是对求下式的最小值：

如果在求最小值的过程中，不施加任何约束，称这种复原为无约束复原，或称为非约束复原。（6.2）（6.4）6.2.1无约束最小二乘方恢复（6.3）由于有：∴

（6.8）6.2.1无约束的最小二乘方恢复根据极值条件

（6.8）给式（6.8）两端同乘以得则有（6.10）

当图像矩阵的尺寸满足N=M，且和满秩非奇异（即可逆）时，则有

式（6.11）说明：当已知H时，便可由g

求出

的估值。（6.11）

如果对式（6.11）两边取傅里叶变换，可以证明有：

式对式（6.7）的结果求傅立叶反变换，就可得到恢复后的图像：（6.11）6.2.2

逆滤波图像恢复方法u=0,1,…,M-1；v=0,1,…,N-1）（6.12）u=0,1,…,M-1；v=0,1,…,N-1）（6.13）

【例】设有退化函数对原图像进行退化和逆滤波图像恢复的结果如图6.3所示。6.2.2

逆滤波图像恢复方法

【例】逆滤波（Inversefilter）图像复原matlab编程。clc;clearall;closeall;img0=imread('d:\0_matlab图像课编程\lena.jpg');subplot(1,3,1);imshow(img0);title('原图像');%显示原图像%1.频率域方法生成退化图像f=im2double(img0);FTS_img=fftshift(fft2(f));%傅里叶变换-平移频谱图为中心对称%执行退化：退化函数H(u,v)=exp(-0.0025*((u-M/2).^2+(v-N/2).^2).^(5/6))[M,N]=size(img0);foru=1:Mforv=1:NH(u,v)=exp(-0.025*((u-M/2).^2+(v-N/2).^2).^(5/6));endendG_FTS_img=FTS_img.*H;%G=F*Hdegraded_I=ifftshift(G_FTS_img);%反移频谱中心degraded_img0=real(ifft2(degraded_I));

%傅里叶逆变换，并取变换结果的实部degraded_img=uint8(abs(degraded_img0)*256);%生成的退化图像subplot(1,3,2);imshow(degraded_img);title(‘生成的退化图像’);

%2.对退化图像进行逆滤波Filter_R=80;%设置逆滤波半径ifFilter_R>M/2

F_Id_img=G_FTS_img./(H+eps);%全滤波F=G/Helse%对滤波半径范围内进行滤波

fori=1:M

forj=1:N

ifsqrt((i-M/2).^2+(j-N/2).^2)<Filter_R

F_Id_img(i,j)=G_FTS_img(i,j)./(H(i,j)+eps);

end

end

endendInverse_f_I=ifftshift(F_Id_img);%反移频谱中心Inverse_f_img0=ifft2(Inverse_f_I);%傅里叶反变换Inverse_f_img=uint8(abs(Inverse_f_img0)*255);subplot(1,3,3);imshow(Inverse_f_img);title('滤波半径为80的逆滤波复原图像');6.2.2

逆滤波图像恢复方法【例】逆滤波（Inversefilter）图像复原matlab编程。6.2.3无约束图像恢复的病态性（6.12）6.2.3无约束图像恢复的病态性6.3维纳滤波图像恢复6.3.1

有约束最小二乘方恢复

为了克服图像恢复过程中的病态性，常常会在图像的恢复过程施加某种约束，于是就引入了有约束的最小二乘方恢复方法。6.3.1

有约束最小二乘方恢复1、有约束的最小二乘方恢复方法

有约束的最小二乘方恢复方法需要知道噪声的模平方，有学者已经证明，

能用噪声的均值

和方差表示为：（6.16）

也就是说，只要知道噪声的均值和方差，就可实现有约束的最小二乘方图像恢复方法。1、有约束的最小二乘方恢复方法

下面先讨论有约束恢复的一般表示形式。

设对原图像施加某一线性运算，求在约束条件下，使

为最小的原图像f

的最佳估计

。

这一问题实际上是求极值问题，常采用拉格朗日乘数法来实现。也就说，要寻找一个，使得构造的辅助函数（准则函数）（6.17）（6.18）为最小。6.3.1

有约束最小二乘方恢复

寻找一个，使得构造的辅助函数（准则函数）（6.18）为最小，也即令（6.19）其中：是拉格朗日乘子，是约束项，如果找到

为最小的原图像f的最佳估值

时，就为0。设，并带入式（6.19）可得（6.21）（6.19）所以有6.3.1

有约束最小二乘方恢复（6.21）6.3.1

有约束最小二乘方恢复1、有约束的最小二乘方恢复方法6.3.2

维纳滤波图像恢复方法

维纳滤波的总体思路是寻找图像

的一种估计，使得

和之间的均方误差最小。6.3.2

维纳滤波图像恢复方法设和分别表示原图像和噪声的自相关矩阵，对图像取线性运算

如果用E{}一般地表示自相关矩阵，则有（6.22）（6.23）（6.24）6.3.2

维纳滤波图像恢复方法

将式（6.22）代入式（6.21）可得

如果式（6.25）可使

的模最小，也即使噪声和信号的比对复原图像的影响最小。

式（6.25）即是最小均方误差滤波恢复方法的表示式。（6.22）（6.25）6.3.2

维纳滤波图像恢复方法

有学者已经证明，当式（6.25）中r=1时，即可得下式的（标准）维纳滤波器公式:（6.26）其中：为噪声的功率谱，为图像的功率谱。

由式(6.26)可知，当没有噪声时，

，维纳滤波器就可简化成逆滤波器；当有噪声时，维纳滤波器也可用信噪功率比作为修正函数对逆滤波器进行修正，可在均方误差最小的意义上提供最佳恢复。6.3.2

维纳滤波图像恢复方法

通常将噪声假设为白噪声，则噪声的功率谱

为常数，也即认为（6.27）由于通常难以估计，一种近似的解决方法是用一个系数K

来代替，，这样（6.26）就可用下式来近似（6.28）其中，K

是根据信噪比的某些先验知识来预先设定的一个常数。

【例】设有退化函数对原图像进行退化和逆滤波图像恢复的结果如图6.4所示。6.3.2

维纳滤波图像恢复方法

【例】维纳滤波（WienerFiltering）图像复原matlab程序。clc;clearall;closeall;img0=imread('d:\0_matlab图像课编程\lena.jpg');subplot(1,4,1);imshow(img0);title('原图像');%1.频率域方法生成退化图像f=double(img0);FTS_img=fftshift(fft2(f));%傅里叶变换-平移频谱图为中心对称%执行退化：退化函数H(u,v)=exp(-0.005*((u-M/2).^2+(v-N/2).^2).^(5/6))[M,N]=size(img0);foru=1:Mforv=1:NH(u,v)=exp(-0.005*((u-M/2).^2+(v-N/2).^2).^(5/6));endendG=FTS_img.*H;%G=F?HG_I=ifftshift(G);

%反移频谱中心G_img=real(ifft2(G_I));

%傅里叶反变换，并取变换结果的实部degraded_img=uint8(G_img);

%生成的退化图像subplot(1,4,2);imshow(degraded_img);title(‘退化图像’);

%2.对退化图像叠加高斯噪声f1=imnoise(uint8(G_img),'gaussian',0,0.001)subplot(1,4,3);imshow(uint8(f1));title('模糊退化且添加高斯噪声的图像');%3.对退化图像进行维纳滤波F1=fftshift(fft2(f1));%对退化图像进行傅里叶变换，移频谱中心对称K=0.1;foru=1:Mforv=1:NH(u,v)=exp(-0.005*(((u-M/2)^2+(v-N/2)^2)^(5/6)));H0(u,v)=(abs(H(u,v)))^2;H1(u,v)=H0(u,v)/(H(u,v)*(H0(u,v)+K));endendF2=F1.*H1;%matlab中的点乘，也即卷积乘法

f2=ifft2(ifftshift(F2));%反移频谱中心-傅里叶反变换f2_img=uint8(f2);%转换成0-255图像数据subplot(1,4,4);imshow(f2_img);title('维纳滤波复原图');6.2.2

逆滤波图像恢复方法

【例】维纳滤波（WienerFiltering）图像复原matlab程序6.4匀速直线运动模糊的恢复

匀速直线运动模糊

6.4匀速直线运动模糊的恢复在图像的运动分析中，比较简单的情况就是对由于相机镜头和对象之间在曝光瞬间的相对运动而造成的图像模糊的恢复。这种情况或者发生在相机处于静止状态而目标在场景中运动，或者发生在相机移动而目标处于静止状态这两种情况。

其中，最简单的相机和目标的相对运动，即可以看成是匀速直线运动而造成的模糊图像的恢复问题。

6.4匀速直线运动模糊的恢复基于式(6.31)和(6.32)，就可推出匀速直线运动模糊的恢复模型{自学}6.5图像噪声与被噪声污染图像的恢复

图像噪声去除是图像处理技术中，图像增强与图像恢复的交叉研究问题，一般认为是一种图像预处理技术。

为了在有噪声的情况下恢复图像，就需要了解噪声的统计性质，以及噪声与图像之间的相关性质。

数字图像常会因受一些随机误差而退化，这种退化通常称为噪声（noise）。图像噪声通常是一种空间上不相联系的离散和孤立的像素的变化现象。图像噪声也是一种图像退化因素。对图像来说，噪声是一种外部干扰。但噪声本身也是一种信号（携带了噪声源的信息）。6.5.1图像噪声6.5.1图像噪声图像噪声是一个随机量，所以噪声一般用其概率特征来描述。1.常用图像噪声的概率密度函数

对相关概念的理解：1、概率分布函数的几何意义概率分布函数是对一个概率事件的描述，反映了事件的全貌；用于统计有多少情况可以发生。

2、概率密度函数的几何意义

概率密度函数可以理解为每个事件出现的机会大小，反应的是每一种情况有多大的机会发生。6.5.1图像噪声对相关概念的理解：3、概率分布函数与概率密度函数区别与联系?

解答：一元函数下：概率分布函数是概率密度函数的变上限积分,就是原函数；概率密度函数是概率分布函数的一阶导函数。

多元函数下：联合分布函数是联合密度函数的重积分，联合密度函数是联合分布函数关于每个变量的偏导。6.5.1图像噪声

（1）高斯噪声高斯噪声是一种源于电子电路噪声和由低照明度或高温带来的传感器噪声。高斯噪声也称为正态噪声，其概率密度函数为：（6.59）其中，高斯随机变量z表示灰度值；μ表示z的平均值或期望值；σ表示z的标准差，而标准差的平方σ2称为z的方差。6.5.1图像噪声—常用图像噪声的概率密度函数

注意：高斯随机变量Z的值不一定要和原图像的范围一致，例如可以是0-10。

高斯噪声是白噪声的一个特例（幅度符合高斯分布）。

白噪声：当图像面上不同点的噪声互不相关时，称为白噪声。

白噪声的特点：频率均匀覆盖整个频谱，功率谱为常量，也即其强度不随频率的增加而衰减。

高斯噪声的形成：电子噪声、弱光照/温度条件下的传感器噪声

（1）高斯噪声（续1）6.5.1图像噪声—常用图像噪声的概率密度函数

（2）瑞利噪声瑞利噪声的概率密度函数为：

概率密度的均值和方差分别为：

（6.60）

（6.61）

（6.62）

瑞利噪声的形成：深度成像、超声波图像。

6.5.1图像噪声—常用图像噪声的概率密度函数

（3）均匀分布噪声均匀分布噪声的概率密度函数为：概率密度的期望值和方差分别为：

（6.63）

（6.64）

（6.65）

6.5.1图像噪声—常用图像噪声的概率密度函数

（4）脉冲噪声（椒盐噪声）噪声脉冲可以是正的或负的，一般假设为a和b。都是“饱和”值的双极性脉冲噪声也称椒盐噪声。(双极)脉冲噪声的概率密度为：（6.66）

含义:

脉冲噪声在Pa或Pb均不可能为零，且在脉冲可能是正的，也可能是负值的情况下，称为双极脉冲噪声。

椒盐噪声的形成：快速瞬变、误切换。6.5.1图像噪声—常用图像噪声的概率密度函数

（4）脉冲噪声（椒盐噪声）(续1)如果b>a，灰度b的值在图像中将显示一个亮点，而灰度a的值在图像中将显示一个暗点。如果Pa或Pb均不可能为零时，脉冲噪声值就类似于随机分布在图像上的胡椒和盐粉微粒，所以双极脉冲噪声也称为椒盐噪声.式（6.66）表示的脉冲噪声如果Pa或Pb为零，则脉冲噪声称为单极脉冲噪声。（6.42）6.5.1图像噪声—常用图像噪声的概率密度函数

（4）脉冲噪声（椒盐噪声）(续1)通常情况下脉冲噪声总是数字化为允许的最大值或最小值，所以负脉冲以黑点(胡椒点)出现在图像中，正脉冲以白点（盐点）出现在图像中。

6.5.1图像噪声—常用图像噪声的概率密度函数

（4）脉冲噪声（椒盐噪声）(续2)

（5）其他噪声周期噪声：图像采集过程中的电子或电磁干扰指数和Gamma分布：激光成像6.5.1图像噪声—常用图像噪声的概率密度函数

图像中噪声的概率密度函数举例1：

加高斯噪声图像加瑞利噪声图像

加gamma噪声图像

叠加噪声后图像的直方图原图对应的直方图6.5.1图像噪声按噪声信号与图像信号的相关性可以把噪声分为两类：加性噪声和乘性噪声。2.图像噪声的分类

6.5.1图像噪声—图像噪声的分类

（1）加性噪声

加性噪声是指叠加在图像上的噪声，也即它们与信号的关系是相加的，它与图像信号的有无及灰度值大小无关，即使信号为零，它也会存在。这种在图像通过信道传输时，独立于图像信号的噪声称为加性噪声（additivenoise）；含有这类噪声的图像一般表示为（6.67）

其中，噪声

和输入图像

是相互独立的变量。6.5.1图像噪声—图像噪声的分类

（2）乘性噪声

乘性噪声是指对有用信号有调幅作用的噪声，也即它们与信号的关系是相乘的，该类噪声的幅值与图像本身的灰度值有关；但当有用信号为零时，该噪声的干扰影响就不存在了，也即信号在它在，信号不在它也就不在了。这种噪声称为乘性噪声（multiplicativenoise），含有这类噪声的图像一般表示为（6.68）

比如，电视光栅退化和胶片材料的退化都是乘性噪声。6.5.1图像噪声—图像噪声的分类

（2）乘性噪声

由于乘性噪声的处理是比较复杂的，所以通常总是假定信号或图像和噪声是互相独立的，也即一般都假设噪声是加性噪声。

在红外图像的成像过程中，由于红外波的相互干涉作用，往往存在有散斑噪声，也即这种噪声在图像上呈斑点分布状。由于散斑噪声既包含乘性噪声的成分，也包含加性噪声的成分，所以含有这类噪声的图像一般表示为（6.69）

6.5.1图像噪声

下面以叠加加性噪声为例进行说明。

若假设输入图像f(x,y)的灰度级取值范围为[0，L-1]，则产生加性零均值高斯噪声的具体步骤为：3.给图像叠加噪声的方法①取图像灰度值的标准差＞0。②对每一水平相邻的像素（x，y）、（x，y+1）产生一对位于[0，1]范围的独立的随机数，。③计算（6.70）④计算

和

⑤置

⑥跳转到③，直到扫描完所有像素为止。

（6.71）（6.72）6.5.1图像噪声【例】给原图像叠加噪声密度为0.05的椒盐噪声的实例。3.给图像叠加噪声的方法6.5.1图像噪声

【例】

给图像叠加生成椒盐噪声的matlab程序。3.给图像叠加噪声的方法6.5.1图像噪声

◆典型的Matlab叠加噪声函数

（1）Matlab中用椒盐噪声污染图像f的函数及编程

函数：

imnoise(f,'salt&pepper',d)

其中，d是噪声密度（即包括噪声值的图像区域的百分比），含义是大约有d×numel（f）个像素受到影响，也即有d×100%的像素受到污染。一般的默认噪声密度为0.01—0.05。

Matlab用椒盐噪声污染图像的编程实现

clc;clearall;closeall;img0=imread('d:\0_matlab图像课编程\lena.jpg');f=im2double(img0);imgnoise=(rand(size(f)));%生成随机数组d=0.05;%需要的椒盐噪声密度result_f=f;[h,w]=size(f);fori=1:h

forj=1:w

if(imgnoise(i,j)>(1-d/2))

result_f(i,j)=1;

else

if(imgnoise(i,j)<d/2)

result_f(i,j)=0;

else

continue;

end

end

endendsubplot(1,2,1);imshow(img0);title('原图像');%显示原图像subplot(1,2,2);imshow(result_f);title(‘加0.05椒盐噪声图像’);

6.5.1图像噪声

◆典型的Matlab叠加噪声函数

（2）Matlab中给图像叠加加性零均值高斯噪声的函数imnoise(img0,'gaussian',0.02)；

给图像加噪声密度为0.02的高斯噪声imnoise(img0,'speckle',0.02)；

给图像加噪声密度为0.02的斑点噪声6.5.1图像噪声4.图像的信噪比

设

是含有噪声的图像,是没有被噪声污染的图像，则图像的信噪比（signal-to-noiseratio，SNR）的定义为：

SNR=

设信噪比的对数表示形式如式（6.50）所示，单位为分贝。

（6.73）

（6.74）

设f(x,y)是一幅原图像，经过退化过程H(x,y)后，形成的退化图像为g(x,y)。当一幅图像中存在的唯一退化因素是噪声n(x,y)，并且噪声与图像不相关时，则在空间域中的退化图像就可以表示为：

在图像中仅存