基于MATLAB的图像复原.doc

1、基于 MATLAB 的图像复原摘要随着信息技术的发展，数字图像像已经充斥着人们身边的任意一个角落。由于图像的传送、转换，或者其他原因，可能会造成图像的降质、模糊、变形、质量下降、失真或者其他情况的图像的受损。本设计就针对“图像受损”的问题，在 MATLAB 环境中实现了利用几何失真校正方法来恢复被损坏的图像。几何失真校正要处理的则是在处理的过程，由于成像系统的非线性，成像后的图像与原图像相比，会产生比例失调，甚至扭曲的图像。图像复原从理论到实际的操作的实现，不仅能改善图片的视觉效果和保真程度，还有利于后续的图片处理，这对医疗摄像、文物复原、视频监控等领域都具有很重要的意义。关键字：图像复原；M

2、ATLAB；几何失真校正目录摘要 11 MATLAB 6.x 信号处理 12 图像复原的方法及其应用 132.1 图像复原的方法 132.2 图像复原的应用 143 几何失真校正实现 153.1 空间变换 153.1.1 已知 和 条件下的几何校正 .16yxr,s,3.1.2 和 未知条件下的几何失真 163.2 灰度插值 183.3 结果分析 19参考文献 20附录 2101 MATLAB 6.x 信号处理（1）对 MATLAB 6 进行了简介，包括程序设计环境、基本操作、绘图功能、M 文件以及 MATLAB 6 的稀疏矩阵这五个部分。 MATLAB 的工作环境有命令窗口、启动平台、工作空

3、间、命令历史记录与当前路径窗口这四部分。M文件的编辑调试环境有四个部分的设置，分别是：Editor/Debugger 的参数设置，字体与颜色的设置，显示方式的设置，键盘与缩进的设置。MATLAB 采用路径搜索的方法来查找文件系统的 M 文件，常用的命令文件组在 MATLAB 文件夹中，其他 M 文件组在各种工具箱中。基本操作主要是对一些常用的基本常识、矩阵运算及分解、数据分析与统计这三方面进行阐述。MATLAB 的基本操作对象时矩阵，所以对于矩阵的输入、复数与复数矩阵、固定变量、获取工作空间信息、函数、帮助命令进行了具体的描述。矩阵运算是 MATLAB 的基础，所有参与运算的数都被看做为矩阵。

4、MATLAB 中共有四大矩阵分解函数：三角分解、正交分解、奇异值分解以及特征值分解。数据分析与统计包括面向列的数据分析、数据预处理、协方差矩阵与相关系数矩阵、曲线拟合这四部分。MATLAB 中含有丰富的图形绘制寒素，包括二维图形绘制、三维图像绘制以及通用绘图工具函数等，同时还包括一些专业绘图函数，因此其具有很强大的绘图功能。简单的二维曲线可以用函数 plot 来绘制，而简单的三维曲线图则用 plot3 来绘制。在绘制图形时，MATLAB 自动选择坐标轴表示的数值范围，并用一定的数据间隔标记做标注的数据，当然自己也可以指定坐标轴的范围与数据间隔。专业的绘图函数有绘梯度图制条形图、饼图、三维饼图、

5、箭头图、星点图、阶梯图以及等高线。M 文件时用户自己通过文本编辑器或字处理器生成的，且其之间可以相互调用，用户可以根据自己的需要，自我编写 M 文件。M 文件从功能上可以分为底稿文件与函数文件两类，其中底稿文件是由一系列 MATLAB 语句组成的，而函数文件的第一行必须包含关键字“function”，二者的区别在于函数文件可以接受输入参数，并可返回输出参数，而底稿文件不具备参数传递的功能；在函数文件中定义及使用的变量大都是局部变量，只在本函数的工作区内有效，一旦退出该函数，即为无效变量，而底1稿文件中定义或使用的变量都是全局变量，在退出文件后仍为有效变量。稀疏矩阵是一种特殊类型的矩阵，即矩阵中

6、包括较多的零元素。MATLAB 对稀疏矩阵的存储有两种模式：完全存储和稀疏存储。函数 full 和 sparse 是一对用来对矩阵存储模式进行转换的内部矩阵。函数 sparse 可以用一组非零元素直接创建一个稀疏矩阵，其格式如下：S=Sparse(i，j，s，m， n)其中 i 和 j 都为数组，分别代表矩阵中非零元素的行号和列号；s 是一个全部元素为非零的数组，元素在矩阵中排列的位置为(i,j)；m 为输出矩阵的稀疏矩阵的行数，n 为输出的稀疏矩阵的列数。函数 sparse 还有一种格式为：S=Sparse(i，j，s，m，n，nzmax)其中，参数 i、j、s、m、n 的说明与上面的格式相

7、同，参数 nzmax 用来设置矩阵中非零元素的最大数目。Full 函数可以讲稀疏矩阵变为一般矩阵。将一个矩阵的对角线元素保存在一个稀疏矩阵中，可以使用函数 spdiags 实现，其语法格式为：S=spdiags(B，d，m ，n)创建一个大小为 的稀疏矩阵 S，其非零元素来自矩阵 B 中的元素且按对nm角线排列。参数 d 指定矩阵 B 中用于生成稀疏矩阵 S 的对角线位置。矩阵的主对角线可以认为是第 0 条对角线，每向右移动一条对角线编号加 1，向左下移动一条对角线编号减 1，也就是说 B 中 j 列的元素填充矢量 d 中第 j 个元素所指定的对角线。用外部文件创建的文本文件，如果其中的数据按

8、 3 个列排列，可以将这个文本文件载入工作空间，用于创建一个稀疏矩阵。MATLAB 提供了专门针对稀疏矩阵的函数。处理稀疏矩阵时，计算的复杂程度与稀疏矩阵中的非零元素的个数成正比，计算的复杂程度也与矩阵的行列大小有关，稀疏矩阵的乘法、乘方，包含一定次数的线性方程等，都是比较复杂的运算。稀疏矩阵的行交换与列交换可以用以下两种方法表示：（1）对于交换矩阵 P，对稀疏矩阵 S 的行交换可表示为 ，列交换可以表SP示为 。S（2）对于一个交换矢量 p，p 为一般矢量，包含 1n 个自然数的一个排列。对稀疏矩阵进行行交换，可以表示为 S(p,:)。S(p,:)为列交换形式。对于矩阵 S2的第 i 列进行

9、行交换的形式为 S(p,i)。稀疏矩阵和一般矩阵一样，同样可以进行 LU 分解、Cholesky 分解、QR分解以及一些不完全分解。与一般矩阵的特征值求解函数 eig 不同的是，计算稀疏矩阵的特征值采用函数 eigs。一般矩阵的奇异值分解用函数 svd，对稀疏矩阵额的奇异值分解使用函数 svds。第二章对离散信号进行了详尽的阐述，并就其 MATLAB 的实现作了总结。典型的离散信号有单位抽样序列、单位阶跃系列、正弦序列、复正弦序列、指数序列、随机序列 6 种。单位抽样序列的表达式如下：（1-1）01n又被称为 Kronecker 函数，该信号在离散信号与离散系统的分析与综合中有着重要的作用，在

10、 MATLAB 中可以利用 zeros 函数来实现。如要产生 N 点的单位抽样序列，可通过下列语句实现： ;1;,xNzeros单位阶跃序列的表达式如下：（1-2）0nuMATLAB 中的 ones 函数可以容易实现 N 点单位阶跃序列： 。Nonesx,1正弦序列的表达式如下：（1-3）sfnTAnx2i其 MATLAB 的实现如下所示： faifpixssi;1:0复正弦序列的表达式如下：（1-4）mjen其 MATLAB 的实现如下所示： nwjxNep;1:03指数序列的表达式如下所示：（1-5）nax其 MATLAB 的实现如下所示： ;.:1nxN随机序列在 MATLAB 中是可以

11、很容易实现的，有以下两类：（1）rand(1,N) 产生0 ，1上均匀分布的随机序列；（2）randn(1,N) 产生均值为 0，方差为 1 的高斯随机序列，也就是白噪声序列，其他的分布的随机数可以通过上述随机数的变换而产生的。对离散信号所作的基本运算分别是移位、相加、相乘等等，其 MATLAB 的实现如下所示：（1）信号延迟：给定离散信号 ，若信号 定义为 ，那么，nxnyknxy是信号 在时间轴上右移 k 个抽样周期得到的新序列。nyx（2）信号相加： 。值得注意的是，当序列 和 的长n2112度不等或位置不对应时，首先应使两者位置对齐，然后通过 zeros 函数左右补零使其长度相等后再相

12、加。（3）信号相乘： ，这是样本与样本之间的点乘运算，在x21MATLAB 中可采用 来实现，但两序列应做如相加运算同样的操作。序列.和 同上，相乘后得到序列 。nx12 n（4）信号标量乘： ，其 MATLAB 很容易实现： 。cxny xcy（5）信号翻转： ，在 MATLAB 中可以直接用 fliplr 函数实现此操作。（6）信号和：对于 N 点信号 ，其和的定义为： ，采用nxNnxy1MATLAB 实现所示： 。sumy（7）信号积：对于 N 点信号 ，其积的定义为： ，MATLABxNnxy1实现如下所示： 。prody4（8）信号能量：有限长信号的能量定义为： ，其NnNnXxx

13、E112MATLAB 实现有两种方法： 或者 。;.cojxsumE;.absum对于0 ，1 上均匀分布的随机噪声可以直接利用 MATLAB 中的 rand 函数实现，均值为 0，方差为 1 的高斯随机噪声即白噪声有函数 randn 产生。对于其他分布（如瑞利分布、对数正态分布等）的随机噪声可以通过上述随机数的变换而产生，这些都是噪声的产生方法。MATLAB 中含有丰富的函数用以生产无线电技术以及通讯等领域广泛采用的信号波形，如方波、三角波和线性调频信号等。其中 MATLAB 内部提供的产生信号波形的函数有五种，分别是：SAWTOOTH 函数、SQUARE 函数、SINC 函数、 DIRIC

14、 函数、CHIRP 函数。第三章对离散系统的基本概念作了描述，然后对离散系统的时域与频域表示方法以及相应的 MATLAB 实现进行了具体的阐述，最后介绍了有关离散系统变换的知识。离散系统的定义是：一个离散系统，可以抽象为一种变换，或是一种映射，即把输入序列 变换为输出序列 ： ，式中，nxnynxTT 代表变换。这样，一个离散系统既可以是一个硬件装置，也可以是一个数字表达式。离散系统有四个重要定义，分别是：线性、移不变性、因果性、稳定性。线性的定义是：设一个离散系统对 的响应是 ，对 的响应是x1y1x2，即ny2（1-6）nxTy2211若该系统对 的响应 ，即nx1y21（1-7）nyTn

15、y 212 那么，该系统是线性的。设一个离散系统对 的响应是 ，即 。若满足xyxTn，则该系统是移不变的，同时具有线性和移不变的离散系kyx统成为线性移不变系统，简称为 LSI 系统。一个 LSI 系统，如果它在任意时刻的输出只决定于现在时刻和过去的输入，而与将来的输入无关，那么，该系统是因果系统。稳定性的定义是一个信号 ，如果存在一个实数 R，使得对所nx5有的 n 都满足 ，那么，称 是有界的。对一个 LSI 系统，若输入Rxnx是有界的，输出 也有界，那么该系统是稳定的。稳定性是一个系统能xny否正常工作的先决条件。对于同一个离散系统，可以从时域和频域两个方面来进行描述，也可以对其进行

16、内部描述。另外，频域描述又可以分为频率响应、转移函数、零极点增益与二次分式几种不同的表示形式。一个 LSI 系统可以用一个常系数线性差分方程来描述：（1-8）NkMrrnxbnyany10式中 ，M 是方程的系数。给定输入信号 以及 ,rbka nx系统的初始条件，可求出该差分方程的解 ，从而得到系统的输出。LSI 系y统的频域表示分为频率响应、转移函数、零极点增益、二次分式四部分。频率响应定义是：任意 LSI 系统都可由单位抽样响应 表示，相应的在频域中可nh以用频率响应 来表示，它是 的离散傅里叶变换。若定义 ，则jweHnhjwezLSI 系统的频率响应变为： ，该式为单位抽样 的 Z

17、变换，nzz nh是系统的转移函数。将转移函数的分子、分母分别作因式分解，便可得z出 LSI 系统的零极点增益表示形式：（1-9）NkkMrrpzgzH1式中， 称为系统的增益因子。使分母多项式等于零的 z 值（即g） ，称为系统的极点，同理，使分子多项式等于零的 z 值（即kp,32,1） ，称为系统的极点。通过系统的零极点增益表示形式，很容Mr易判断一个 LSI 系统的稳定性，也就是说，若所有的极点都位于单位圆内，则系统是稳定的。离散系统的内部描述是用状态方程和输出方程来表示的。离散系统的 MATLAB 实现是用函数来实现的，对这些函数进行简要介绍。能产生单位抽样响应的函数有：zeros

18、函数、filter 函数、impz 函数。零极点增益是用roots 函数来实现的。离散系统变换函数包括：tf2zp 函数、tf2ss 函数、zp2tf 函6数、zp2sos 函数、zp2ss 函数、sos2tfz 函数、sos2zp 函数、sos2ss 函数、ss2tf函数、ss2zp 函数、ss2sos 函数。第四章对离散傅里叶变换 DFT、Chirp z 变换、离散余弦变换DCT、Hilbert 变换进行了阐述，并对其进行了 MATLAB 仿真。有限长序列的离散傅里叶变换公式如下所示：nx（1-10）1010NnWkXNnxkxNkn离散傅里叶 DFT 的性质有线性、正交性、圆周移位、圆周

19、卷积、共轭对称性5 个性质。Z 变换是离散系统与离散信号分析与综合的重要工具。一个离散序列 的 z 变换定义为：nx（1-11）nnzxzXZ 变换有线性、序列移位、与指数序列相乘、 的微分、复序列的共轭、zX序列卷积、序列乘积 7 个特性。Chirp Z 变换即线性调频 Z 变换，可用来计算单位圆上任一段曲线上的 Z 变换。做 DFT 时输入的点数 N 和输出点数 M 可以不相等，从而达到频域”细化“的目的。序列 的 Hilbert 变换是 ，nxnx则 ，求出 ，即可构成 的解析信号： mnxnhx12，也可以用 DFT 求出一个信号 的解析信号及 Hilbert 变jnznx换，步骤是：

20、（1）求 的 DFT ，k=0,1，N-1，其中 对应负数频xkX1,2Nk率。7（2）令 1,2,002NkXkZ（3）对 做逆 DFT，得到 的解析信号 。nxnz（4）由 ，得 。kXjkZxjHilbert 变换具有两个性质，分别是：序列 通过 Hilbert 变换器后，信号频谱的幅度不发生变化；序列 与其 Hilbert 变换 是正交的。nxnx第五章对 IIR 系统、FIR 系统结构、离散系统的 Lattice 结构进行了具体描述。一个 N 阶 IIR 数字滤波器的系统函数可以表示为：（1-12）NkkMzabzH10其差分方程为：（1-13）NkMk nyanxby00无限长单位

21、取样响应 IIR 数字滤波器的主要特点是：（1）单位取样响应是无限长的，即 ， 。h-（2）系数函数 在有限平面 z 上有极点存在。zH（3）结构上存在着输出到输入的反馈网络，即结构式递归的。实现同一个系统函数 ，可以用不同的结构形式，它的主要的结构形式有直接 型、直接 型、级联型与并联型四种。有限长单位取样响应 FIR 滤波器突出特点是单位取样 仅有有限个非零值，即 为一个 N 点序列nhnh，其中系统函数为：10Nn（1-14）10NnnzzH在 Z=0 处有 N-1 阶极点，而没有除 Z 平面原点（Z=0）外的极点。FIRzH8滤波器的结构主要是非递归结构，没有输出反馈。FIR 系统的结

22、构有直接型、级联型两种。 Lattice 结构是一种新的系统结构形式，在功率谱估计、语音处理、自适应滤波等方面已经得到了广泛的应用。这里从零点系统和全极点系统进行讨论。一个 M 阶的 FIR 系统的转移函数 可写为：zH（1-15）MiiiibzbzBH10系数 表示 M 阶 FIR 系统的第 i 个系数。ib第六章是用 MATLAB 来进行 IIR DF（IIR 数字滤波器）的设计。先对数字滤波器进行简介。数字滤波器是数字信号处理的重要基础，在对信号的过滤、检测与参数的估计等信号处理中，数字滤波器是使用最为广泛的一种线性系统。数字滤波器是对数字信号实现滤波的线性是不变系统。设计数字滤波器包括

23、以下几个步骤：（1）按照实际任务的要求，确定滤波器的性能指标。（2）用一个因果、稳定的离散线性时不变的系统函数去逼近这一性能指标。根据不同的要求可以用 IIR 系统函数，也可以用 FIR 系统函数去逼近。（3）利用有限精度算法实现系统函数，这里包括结构的选择、字长选择等。设计一个滤波器，重要的是寻找一个稳定、因果的系统函数去逼近滤波器的技术指标。一个也能过、稳定的模拟滤波器的系统函数 应该满足如下条件：sHa（1）滤波器的单位冲击响应函数 应该是一个实函数，即 是一个具tha a有实系数的 s 的有理函数。（2） 的极点必须分布在 s 平面的左半平面。Ha（3） 的分子多项式的阶数必须小于或者

24、等于分母多项式的阶数。实际上设计模拟原型滤波器是要寻求一个逼近理想低通滤波器的函数，所谓原型低通滤波器是指低通模拟或数字滤波器。模拟低通滤波器的设计方法有：巴特沃斯、切比雪夫和椭圆滤波器。模拟低通巴特沃斯滤波器是以巴特沃斯函数作为滤波器的系统函数，它的幅度平方函数表示为：9（1-16）NcajjH2221式中的 N 为正整数，表示滤波器的阶数。 为通带的截止频率，或 的c sHa3 分贝带宽。模拟的低通巴特沃斯滤波器的设计过程包括以下两个过程：（1）按照给定的通带和阻带指标确定阶数 N。（2）从幅度平方函数确定系统函数 。sHa切比雪夫 型滤波器在通带内幅度特性是等波纹的，在阻带是单调的，而切

25、比雪夫 则相反，它在通带内是单调的，在阻带内是等波纹的。切比雪夫滤波器由 3 个参数需要确定： 、 和 N，其中 是给定的截止频率， 由容许cc 的通带波纹或通带的幅度误差 确定。椭圆滤波器是采样有限零点设计的滤波器，它能更好的逼近理想的低通滤波器。它的幅度平方函数为：（1-17）221NaJjH式中的 是雅可比椭圆函数， 是与通带衰减有关的函数，阶数 N 等NJ于通带和阻带内最大点和最小点的总和。脉冲响应不变法的设计原理是使得数字滤波器的单位取样响应序列 模仿模拟滤波器的冲激响应 。双线性变nhtha换化是使得数字滤波器的频率响应模仿模拟滤波器的频率响应模仿模拟滤波器的频率响应的一种方法。这

26、种方法的基本思路是：首先将整个 s 平面压缩到平面的一条带宽为 （从 到 ）的横带里，然后通过标准的变换1sT2T关系 将横带变换成整个 z 平面上去，这便得到 s 平面与 z 平面之间的一Tsez1一对应的单值关系。第 7 章对有限长单位脉冲响应数字滤波器进行介绍。窗函数在设计 FIR数字滤波器中有很重要的作用，正确的选择窗函数可以提高所设计的数字滤波器的性能，或者在满足设计要求的情况下，减小 FIR 数字滤波器的阶数。常见的窗函数有矩形窗、三角窗、布拉克曼窗、汉宁窗、海明窗、凯塞窗、巴特里特窗、切比雪夫窗 8 种。设计 FIR 数字滤波器最简单的方法是窗函数法，通常也称为傅里叶级数法。FI

27、R 滤波器的设计同 IIR 数字滤波器的设计一样，10首先给出要求的理想滤波器的频率响应 ，设计一个 FIR 数字滤波器频jdeH率响应 ，去逼近理想的频率响应 。然而窗函数法设计 FIR 数字jeHj滤波器是在时域进行的，因而必须由理想的频率响应 推导出对应的单jde位取样响应 ，在设计一个 FIR 数字滤波器的单位取样响应 去逼近nhd nh。频率取样法设计 FIR 数字滤波器是从频域出发，根据频域的采样定理，d对给定的理想滤波器的频率响应 进行等间隔的采样：jdeH（1-18） 1,0/2 NkeHNkjd 等波纹切比雪夫法是采用最大误差最小准则得到最佳数字滤波器，并且最佳解是唯一的。切

28、比雪夫逼近法设计 FIR 数字滤波器的过程是：（1）规定所需的频率响应 ，加权函数 和滤波器的单位脉冲响应jdejeW的长度 N。（2）形成 、 和 。jePjWjP（3）利用雷米兹多重交换算法求解逼近问题。（4）计算滤波器的单位脉冲响应 。nh采用切比雪夫逼近设计方法能够得到既有严格线性相位，又有很好衰减特性的滤波器，因此，切比雪夫逼近法在滤波器设计中占有很重要的位置。第 8 章讲述了功率谱估计的问题。功率谱估计方法可以分为经典谱估计和现代谱估法两种，经典谱估计法又可分为直接法和间接法，Bartlett 法和Welch 法是直接法的改进。随机序列 自相关函数估计的两种形式为：nx（1-19）

29、10mNnxxxR用 FFT 计算自相关函数的一般步骤：（1）对 补 N 个零，得 ，对 做 DFT 得 ，k=0，1，2N-nxxkX1；11（2）求 的幅平方，然后除以 N，得 ；kX 21kX（3）对 做逆变换，的 。21N mRx在 MATLAB 中，函数 xcorr 用来进行自相关函数估计，且为局域上述 FFT 的快速算法，其格式为：C=XCORR(A,flag)，该函数返回长度为 2N-1 的自相关序列。AR 模型功率谱估计是现代谱估计的主要内容。 AR 模型又称为自回归模型，它是一个全极点的模型，该模型现在的输出时现在的输入和过去输出的加权和。AR 模型谱估计的性质有：AR 谱的

30、平滑特性、AR 谱的分辨率。基于矩阵特征分解的功率谱估计包括特征向量估计与 MUSIC 估计，这两种估计方法均为非参数估计方法，特征向量估计主要使用混有白噪声的正弦信号的功率皮估计，而 MUSIC 估计适合更为普遍情况下正弦信号参数估计的方法。函数 Pmusic 为 MUSIC 估计，而函数 Peig 为特征向量估计。122 图像复原的方法及其应用2.1 图像复原的方法在图像的获取、传输和保存的过程中，由于各种原因，如大气的湍流效应、传感器特性的非线性、成像设备与物体之间的相对运动、摄像设备中光学系统的衍射、胶片颗粒噪声和感光胶卷的非线性、光学系统的像差以及电视摄像扫描的非线性等所引起的几何失

31、真，都可能造成图像的畸变和失真。通常，由于这些因素引起的质量下降称为图像退化 1。图像退化的表现是图像出现模糊、失真以及附加噪声等。由于图像的退化，在图像接收端显示的图像已经不再是发送的原始图像，图像的效果明显变差，因此我们可以采用一些技术手段来尽可能的减轻甚至消除图像质量的下降，还原图像的本来面目，这就是图像复原 2。图像复原是图像处理领域中一种非常重要的处理技术，与图像增强等其他基本图像处理技术类似，也是以获取视觉质量程度的改善为目的，所不一样的地方是图像复原过程实际上是一个估计的过程，需要根据一些特定的退化模型，对退化图像进行校正。简而言之，图像复原的处理过程就是提升退化图像的质量，并通

32、过图像质量的提高来达到图像在视觉上的改善。因为引起图像退化的原因很多，且性质各不相同，目前还没有统一的复原方法，许多研究人员根据不同的应用环境，采取了不同的退化模型 3-4、估计准则和处理技巧，所以得到了各种复原方法。图像复原的算法是整个技术的核心内容。目前，国内在这方面的研究才刚起步，而国外已经取得了较好的成果。早期的图像恢复是用光学的方法对失真的测试图像进行复原，但是关于数字图像复原技术的研究则是从对天文观测图像的后期处理中才逐渐发展起来的。其中一个成功例子是在 1964 年，NASA的喷气推进实验室用计算机处理有关月球的照片。照片是用电视摄像机在空间飞行器上拍摄的，图像的复原包括消除噪声

33、和干扰等因素，校正几何失真和对比度损失以及反卷积。另一个典型的例子则是对肯尼迪遇刺事件现场照片的处理。由于事情发生的太突然，照片是在相机移动地过程中拍摄的。图像复原的主要目的就是消除移动造成的失真。还有其在在医学领域，图像复原的技术能13广泛的应用于 X 光，CT，B 超等成像系统，用来抑制各种医学成像系统或图像获取系统的噪声，改善医学图像的分辨率 5。现今的复原方法有：维纳滤波、逆滤波、最小二乘滤波、几何失真校正等多种复原方法。随着数信号处理和图像处理的发展，新的复原算法不断出现，不同的复原方法所需的条件是不相同的，所以在应用中可以根据具体情况加以选择。2.2 图像复原的应用目前国内外对于图

34、像复原技术的研究和应用主要集中在空间探索、天文观测、物质研究、遥感遥测、军事科学、生物科学、医学影象、刑事侦察、交通监控等领域。如在生物方面，主要是用于生物活体细胞内部组织的三维再现和重构，通过复原荧光显微镜所收集的细胞内部的逐层切片图，来重现细胞内部构成；医学方面，如对肿瘤周围组织进行显微观察，以获取肿瘤安全切缘与癌肿原先部位之间关系的定量数据；天文方面则采用迭代盲反卷积进行气动光学效应图像复原研究等。143 几何失真校正实现在图像的获取或者显示的过程中通常会产生几何失真，比如成像系统有一定的几何非线性，其主要原因是视像管摄像机和阴极射线管显示器的扫描偏转系统有一定的非线性，所以会造成图像的

35、枕形失真或桶性失真。另外，由卫星拍摄的地球表面的图像往将往覆盖较大的面积，由于地球表面呈球形，这样拍摄的图像同样会有比较大的几何失真。几何失真一般分为系统失真和非线性失真。系统失真是有规律的、能预测的；而非线性系统的失真则是随机的。当对图像进行定量分析的时候，就要先对畸变的图像进行准确的几何失真校正 12（即把几何失真的图像校正成与原图像相差无几的图像） ，以免影响分析精度。基本的方法是：首先，要建立几何校正的数学模型；其次，用已知的条件来确定模型参数；最后，根据畸变模型对图像进行几何校正。通常可以分为以下两步：（1）图像空间坐标的变换；（2）确定校正空间各像素的灰度值（灰度内插） 。3.1

36、空间变换假设一幅图像为 ，经过几何失真变成了 ，其中， 表示yxf, vug,vu,的是失真图像的坐标，而不是原图像的坐标 。上述变化可表示为：yx,（3-1）ru,（3-2）yxsv这里， 和 是空间变换，产生了几何失真图像 。yxr,s, vug,若函数 和 已知，则可以依据一个坐标系统的像素坐标算出另一坐标系统的对应像素的坐标。在已知情况下，通常 和 可用多项yxr,s,式来近似：（3-3）10Nijjiyxau（3-4）10ijjibv15式中， 为多项式的次数， 和 为各项系数。Nijaijb3.1.1 已知 和 条件下的几何校正yxr,s,若我们具备先验知识 、 ，则希望将几何畸变

37、图像 恢复yxr,s, vug,为基准几何坐标的图像 。几何校正通常分为直接和间接两种方法。f（1）直接法。先通过 来推出 ，然后计算每个像素yxsvu,vusyrx,的校正坐标值，保持各像素灰度值不变，然后生成一幅校正图像，但其像素的分布是不太规则的，可能会出现疏密不均、像素挤压等现象。（2）间接法。假设复原图像的像素在基准坐标系统内是等距网格的交叉点，从网格交叉点的坐标 出发，可以算出在已知几何畸变图像上的yx,坐标 ，即：vu,（3-5）srvu,虽然点 坐标为整数，但 一般不为整数，不会刚好就处于畸变图yx, v像像素的中心，所以不能直接明确该点的灰度值，而只能由其在几何失真图像的周围

38、像素灰度内插求出，作为对应像素 的灰度值，据此可以获得校正yx,之后的图像。由于间接法内插灰度比较容易，所以在通常情况下采用的是间接法来进行几何失真校正。3.1.2 和 未知条件下的几何失真yxr,s,在这种情况下，通常用原始图像和几何畸变图像上多对“连接点”的坐标来确定 和 。r,s,假设基准图像像素的空间坐标 和待校正图像的对应像素的空间坐标yx,之间的关系用二元多项式来表示，表达式如下：vu,（3-6）10,Nijjiaru（3-7）10,ijjiyxbyxsv16式中， 为多项式的次数， 和 为各项待定系数。Nijaijb对于线性失真：（3-8）yaxyxru010,（3-9）sv1对

39、于一般的（非线性）二次失真：（3-10）xyxayxr1010,（3-11）bbsv利用“连接点”建立失真图像与校正图像之间其他像素空间位置的相应关系，而“连接点”在失真图像和校正图像中的位置是精确已知的。失真图像和校正图像中有四边形区域，这两个四边形的顶点就是相应的“连接点”。假定四边形区域中的几何畸变的过程可以用二次失真方程来表示，即：（3-12）xyaxayxr1010, （3-13）bbs将以上两个表达式代入式（3-1）和式（3-2）中，得：（3-14）xyxu1010（3-15）v因为一共有四对“连接点”，代入式（3-14）和式（3-15）可得 8 个联立方程，由这些方程可以解出 8

40、 个系数 、0a、 、 、 、 、 、 。这些系数就可以形成用于变换四边形区10a10b101b域内所有像素的几何失真模型，即空间映射公式。通常来说，可以把一幅图像分割成一系列覆盖全图的四边形区域的集合，在每个区域寻找足够的“连接点”用来计算进行映射所需的系数。只要有了系数，校正（即复原）图像就不那么困难了。如果想找到非线性失真图像在任一点的 的值，需要简单地知道 在失真图像中的什0,yx0,yxf么地方被映射。为此，可以把 代入式（3-14）和式（3-15）得到几何失0,真坐标 。在无失真图像中被映射到 点的值是 。这样简单0,vu0,vu0,vug地令 ，就得到了复原图像的值。0,gyxf

41、173.2 灰度插值最简单的灰度灰度插值是最近邻插值（也称为零阶插值），该方法实现起来相对简单，但是有时不够精确，甚至经常产生不希望的认为疵点，如高分辨率图像直边的扭曲；对于通常的图像处理，双线性插值很实用；更完善的技术如样条插值、立方卷积内插等可以得到较平滑的结果，但是更平滑的仿真的代价就是增加计算开销。在 MATLAB13软件中生成如图 3.1 所示的 figure 图像，然后在 figure 图像上的失真图像和原始图像上选取九对点，最后依据选取的九对点来进行校正。图 3.1 选取连接点图像18图 3.2 几何失真及校正后图像3.3 结果分析图 3.1 是调用 cpselect 函数，使系

42、统启动交互连接点工具，然后在 figure文件上的图像交替选择对应的坐标点，将其输出到 Workspace 中，进行坐标点的统计；图 3.2 是校正的结果，第二张是几何失真之后的图像，第三张是几何校正后的图像，校正后的图像和原图像相比，几乎没有什么差别，所以几何失真校正复原出来的结果很好，但是它的适用条件是图像的几何失真，其他的失真方式，比如说图像的质量退化，就不可以采用这种复原方法。19参考文献1 朱冠南.基于 MATLAB的图像复原设计J.技术交流,2009:87.2 王婷 .退化图像的复原改进算法研究与实现D.哈尔滨:哈尔滨工程大学,2007.3 何东健.数字图像处理M.西安:西安电子科

43、技大学出版社,2003:263-264.4 胡学龙,许开宇 .数字图像处理M.北京:电子工业出版社,2006:111-132.5 吴亚东.图像复原算法研究D. 成都:电子科技大学,2006.6 杨杰. 数字图像处理及MATLAB实现M.北京:电子工业出版社,2013:127-131. 7 Rafael C Gonzalez,Richard E Woods,Steven L Eddins.Digital image processing(MATLAB)M.Addison-Weslwy, Publishing Company Inc. New York, 1993:189-207.20附录几何失真

44、校正 9 对连接点的选取代码：I=imread(F:素材2.jpg);theta=pi/6; %变换旋转角度T=cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1; %变换矩阵tform1=maketform(affine,T); %几何变换结构x=imtransform(I,tform1,nearest); %以最近邻插值进 行仿射变换F=cpselect(x,I); %调用 cpselect 函数来启动交互 选择连接点工具几何失真校正代码如下所示：I=imread(F:素材2.jpg);subplot(131);imshow(I);

45、title(原图像);k=0.7;theta=pi/6;T=k*cos(theta) k*sin(theta) 0;-k*sin(theta) k*cos(theta) 0;0 0 1;tform1=maketform(affine,T);x=imtransform(I,tform1,nearest);subplot(132);imshow(x);title(最近邻插值的放射变换 );base_points=1.7500 1.2500;188.2500 2.7500;332.7500 1.2500;327.7500 194.7500;177.7500 189.2500;2.7500 200.7

46、500;2.7500 396.2500;166.2500 396.2500;331.2500 394.2500; %从原始图像 I 中选择的 9 个连接点的坐标矩阵input_points=200.9787 3.1263;360.0565 96.1718;486.1182 168.2070;389.3209 336.2898;257.2562 255.2496;103.4310 174.2100;4.3825 346.0441;144.7012 21424.0822;288.0213 508.1233; %从几何失真图像 x 中选择的 9 个连接点的坐标矩阵tform=cp2tform(input_points,base_points,affine); %由 9 对连接点坐标矩阵建立几何变换结构F=imtransform(x,tform,XData,1 333,YData,1 398);subplot(133);imshow(F);title(几何失真校正);