摘要
随着信息技术的发展，数字图像像已经充斥着人们身边的任意一个角落。由于图像的传送、转换，或者其他原因，可能会造成图像的降质、模糊、变形、质量下降、失真或者其他情况的图像的受损。本设计就针对“图像受损”的问题，在MATLAB环境中实现了利用几何失真校正方法来恢复被损坏的图像。几何失真校正要处理的则是在处理的过程，由于成像系统的非线性，成像后的图像与原图像相比，会产生比例失调，甚至扭曲的图像。
图像复原从理论到实际的操作的实现，不仅能改善图片的视觉效果和保真程度，还有利于后续的图片处理，这对医疗摄像、文物复原、视频监控等领域都具有很重要的意义。
关键字：图像复原；MATLAB；几何失真校正
目录
- 对MATLAB 6 进行了简介，包括程序设计环境、基本操作、绘图功能、M文件以及MATLAB 6 的稀疏矩阵这五个部分。MATLAB的工作环境有命令窗口、启动平台、工作空间、命令历史记录与当前路径窗口这四部分。M文件的编辑调试环境有四个部分的设置，分别是：Editor/Debugger的参数设置，字体与颜色的设置，显示方式的设置，键盘与缩进的设置。MATLAB采用路径搜索的方法来查找文件系统的M文件，常用的命令文件组在MATLAB文件夹中，其他M文件组在各种工具箱中。基本操作主要是对一些常用的基本常识、矩阵运算及分解、数据分析与统计这三方面进行阐述。MATLAB的基本操作对象时矩阵，所以对于矩阵的输入、复数与复数矩阵、固定变量、获取工作空间信息、函数、帮助命令进行了具体的描述。矩阵运算是MATLAB的基础，所有参与运算的数都被看做为矩阵。MATLAB中共有四大矩阵分解函数：三角分解、正交分解、奇异值分解以及特征值分解。数据分析与统计包括面向列的数据分析、数据预处理、协方差矩阵与相关系数矩阵、曲线拟合这四部分。MATLAB 中含有丰富的图形绘制寒素，包括二维图形绘制、三维图像绘制以及通用绘图工具函数等，同时还包括一些专业绘图函数，因此其具有很强大的绘图功能。简单的二维曲线可以用函数plot来绘制，而简单的三维曲线图则用plot3来绘制。在绘制图形时，MATLAB自动选择坐标轴表示的数值范围，并用一定的数据间隔标记做标注的数据，当然自己也可以指定坐标轴的范围与数据间隔。专业的绘图函数有绘梯度图制条形图、饼图、三维饼图、箭头图、星点图、阶梯图以及等高线。M文件时用户自己通过文本编辑器或字处理器生成的，且其之间可以相互调用，用户可以根据自己的需要，自我编写M文件。M文件从功能上可以分为底稿文件与函数文件两类，其中底稿文件是由一系列MATLAB语句组成的，而函数文件的第一行必须包含关键字“function”，二者的区别在于函数文件可以接受输入参数，并可返回输出参数，而底稿文件不具备参数传递的功能；在函数文件中定义及使用的变量大都是局部变量，只在本函数的工作区内有效，一旦退出该函数，即为无效变量，而底稿文件中定义或使用的变量都是全局变量，在退出文件后仍为有效变量。稀疏矩阵是一种特殊类型的矩阵，即矩阵中包括较多的零元素。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的说明与上面的格式相同，参数nzmax用来设置矩阵中非零元素的最大数目。Full函数可以讲稀疏矩阵变为一般矩阵。将一个矩阵的对角线元素保存在一个稀疏矩阵中，可以使用函数spdiags实现，其语法格式为：
S=spdiags(B，d，m，n)
创建一个大小为
的稀疏矩阵S，其非零元素来自矩阵B中的元素且按对角线排列。参数d指定矩阵B中用于生成稀疏矩阵S的对角线位置。矩阵的主对角线可以认为是第0条对角线，每向右移动一条对角线编号加1，向左下移动一条对角线编号减1，也就是说B中j列的元素填充矢量d中第j个元素所指定的对角线。用外部文件创建的文本文件，如果其中的数据按3个列排列，可以将这个文本文件载入工作空间，用于创建一个稀疏矩阵。MATLAB提供了专门针对稀疏矩阵的函数。处理稀疏矩阵时，计算的复杂程度与稀疏矩阵中的非零元素的个数成正比，计算的复杂程度也与矩阵的行列大小有关，稀疏矩阵的乘法、乘方，包含一定次数的线性方程等，都是比较复杂的运算。稀疏矩阵的行交换与列交换可以用以下两种方法表示：
- 对于交换矩阵P，对稀疏矩阵S的行交换可表示为
，列交换可以表示为
。
- 对于一个交换矢量p，p为一般矢量，包含1~n个自然数的一个排列。对稀疏矩阵进行行交换，可以表示为S(p,:)。S(p,:)为列交换形式。对于矩阵S的第i列进行行交换的形式为S(p,i)。
稀疏矩阵和一般矩阵一样，同样可以进行LU分解、Cholesky分解、QR分解以及一些不完全分解。与一般矩阵的特征值求解函数eig不同的是，计算稀疏矩阵的特征值采用函数eigs。一般矩阵的奇异值分解用函数svd，对稀疏矩阵额的奇异值分解使用函数svds。
第二章对离散信号进行了详尽的阐述，并就其MATLAB的实现作了总结。典型的离散信号有单位抽样序列、单位阶跃系列、正弦序列、复正弦序列、指数序列、随机序列6种。单位抽样序列的表达式如下：
（1-1）
又被称为Kronecker函数，该信号在离散信号与离散系统的分析与综合中有着重要的作用，在MATLAB中可以利用zeros函数来实现。如要产生N点的单位抽样序列，可通过下列语句实现：
单位阶跃序列的表达式如下：
（1-2）
MATLAB中的ones函数可以容易实现N点单位阶跃序列：
。
正弦序列的表达式如下：
（1-3）
其MATLAB的实现如下所示：
复正弦序列的表达式如下：
（1-4）
其MATLAB的实现如下所示：
指数序列的表达式如下所示：
（1-5）
其MATLAB的实现如下所示：
随机序列在MATLAB中是可以很容易实现的，有以下两类：
- rand(1,N)产生[0，1]上均匀分布的随机序列；
- randn(1,N)产生均值为0，方差为1的高斯随机序列，也就是白噪声序列，其他的分布的随机数可以通过上述随机数的变换而产生的。
对离散信号所作的基本运算分别是移位、相加、相乘等等，其MATLAB的实现如下所示：
- 信号延迟：给定离散信号
，若信号
定义为
，那么，
是信号
在时间轴上右移k个抽样周期得到的新序列。
- 信号相加：
。值得注意的是，当序列
和
的长度不等或位置不对应时，首先应使两者位置对齐，然后通过zeros函数左右补零使其长度相等后再相加。
- 信号相乘：
，这是样本与样本之间的点乘运算，在MATLAB中可采用
来实现，但两序列应做如相加运算同样的操作。序列
和
同上，相乘后得到序列
。
- 信号标量乘：
，其MATLAB很容易实现：
。
- 信号翻转：
，在MATLAB中可以直接用fliplr函数实现此操作。
- 信号和：对于N点信号
，其和的定义为：
，采用MATLAB实现所示：
。
- 信号积：对于N点信号
，其积的定义为：
，MATLAB实现如下所示：
。
- 信号能量：有限长信号的能量定义为：
，其MATLAB实现有两种方法：
或者
。
对于[0，1]上均匀分布的随机噪声可以直接利用MATLAB中的rand函数实现，均值为0，方差为1的高斯随机噪声即白噪声有函数randn产生。对于其他分布（如瑞利分布、对数正态分布等）的随机噪声可以通过上述随机数的变换而产生，这些都是噪声的产生方法。MATLAB中含有丰富的函数用以生产无线电技术以及通讯等领域广泛采用的信号波形，如方波、三角波和线性调频信号等。其中MATLAB内部提供的产生信号波形的函数有五种，分别是：SAWTOOTH函数、SQUARE函数、SINC函数、DIRIC函数、CHIRP函数。
第三章对离散系统的基本概念作了描述，然后对离散系统的时域与频域表示方法以及相应的MATLAB实现进行了具体的阐述，最后介绍了有关离散系统变换的知识。离散系统的定义是：一个离散系统，可以抽象为一种变换，或是一种映射，即把输入序列
变换为输出序列
：
，式中，T代表变换。这样，一个离散系统既可以是一个硬件装置，也可以是一个数字表达式。离散系统有四个重要定义，分别是：线性、移不变性、因果性、稳定性。线性的定义是：设一个离散系统对
的响应是
，对
的响应是
，即
（1-6）
若该系统对
的响应
，即
（1-7）
那么，该系统是线性的。
设一个离散系统对
的响应是
，即
。若满足
，则该系统是移不变的，同时具有线性和移不变的离散系统成为线性移不变系统，简称为LSI系统。一个LSI系统，如果它在任意时刻的输出只决定于现在时刻和过去的输入，而与将来的输入无关，那么，该系统是因果系统。稳定性的定义是一个信号
，如果存在一个实数R，使得对所有的n都满足
，那么，称
是有界的。对一个LSI系统，若输入
是有界的，输出
也有界，那么该系统是稳定的。稳定性是一个系统能否正常工作的先决条件。对于同一个离散系统，可以从时域和频域两个方面来进行描述，也可以对其进行内部描述。另外，频域描述又可以分为频率响应、转移函数、零极点增益与二次分式几种不同的表示形式。一个LSI系统可以用一个常系数线性差分方程来描述：
（1-8）
式中
，M是方程的系数。给定输入信号
以及系统的初始条件，可求出该差分方程的解
，从而得到系统的输出。LSI系统的频域表示分为频率响应、转移函数、零极点增益、二次分式四部分。频率响应定义是：任意LSI系统都可由单位抽样响应
表示，相应的在频域中可以用频率响应
来表示，它是
的离散傅里叶变换。若定义
，则LSI系统的频率响应变为：
，该式为单位抽样
的Z变换，
是系统的转移函数。将转移函数的分子、分母分别作因式分解，便可得出LSI系统的零极点增益表示形式：
（1-9）
式中，
称为系统的增益因子。使分母多项式等于零的z值（即
），称为系统的极点，同理，使分子多项式等于零的z值（即
），称为系统的极点。通过系统的零极点增益表示形式，很容易判断一个LSI系统的稳定性，也就是说，若所有的极点都位于单位圆内，则系统是稳定的。离散系统的内部描述是用状态方程和输出方程来表示的。离散系统的MATLAB实现是用函数来实现的，对这些函数进行简要介绍。能产生单位抽样响应的函数有：zeros函数、filter函数、impz函数。零极点增益是用roots函数来实现的。离散系统变换函数包括：tf2zp函数、tf2ss函数、zp2tf函数、zp2sos函数、zp2ss函数、sos2tfz函数、sos2zp函数、sos2ss函数、ss2tf函数、ss2zp函数、ss2sos函数。
第四章对离散傅里叶变换DFT、Chirp z 变换、离散余弦变换DCT、Hilbert变换进行了阐述，并对其进行了MATLAB仿真。有限长序列
的离散傅里叶变换公式如下所示：
（1-10）
离散傅里叶DFT的性质有线性、正交性、圆周移位、圆周卷积、共轭对称性5个性质。Z变换是离散系统与离散信号分析与综合的重要工具。一个离散序列
的z变换定义为：
（1-11）
Z变换有线性、序列移位、与指数序列相乘、
的微分、复序列的共轭、序列卷积、序列乘积7个特性。Chirp Z变换即线性调频Z变换，可用来计算单位圆上任一段曲线上的Z变换。做DFT时输入的点数N和输出点数M可以不相等，从而达到频域”细化“的目的。序列
的Hilbert变换是
，则
，求出
，即可构成
的解析信号：
，也可以用DFT求出一个信号
的解析信号及Hilbert变换，步骤是：
- 求
的DFT
，k=0,1，…，N-1，其中
对应负数频率。
- 令
- 对
做逆DFT，得到
的解析信号
。
- 由
，得
。
Hilbert变换具有两个性质，分别是：序列
通过Hilbert变换器后，信号频谱的幅度不发生变化；序列
与其Hilbert变换
是正交的。
第五章对IIR系统、FIR系统结构、离散系统的Lattice结构进行了具体描述。一个N阶IIR数字滤波器的系统函数可以表示为：
（1-12）
其差分方程为：
（1-13）
无限长单位取样响应IIR数字滤波器的主要特点是：
- 单位取样响应是无限长的，即
，
。
- 系数函数
在有限平面z上有极点存在。
- 结构上存在着输出到输入的反馈网络，即结构式递归的。
实现同一个系统函数
，可以用不同的结构形式，它的主要的结构形式有直接
型、直接
型、级联型与并联型四种。有限长单位取样响应FIR滤波器突出特点是单位取样
仅有有限个非零值，即
为一个N点序列
，其中系统函数为：
（1-14）
在Z=0处有N-1阶极点，而没有除Z平面原点（Z=0）外的极点。FIR滤波器的结构主要是非递归结构，没有输出反馈。FIR系统的结构有直接型、级联型两种。 Lattice结构是一种新的系统结构形式，在功率谱估计、语音处理、自适应滤波等方面已经得到了广泛的应用。这里从零点系统和全极点系统进行讨论。一个M阶的FIR系统的转移函数
可写为：
（1-15）
系数
表示M阶FIR系统的第i个系数。
第六章是用MATLAB来进行IIR DF（IIR数字滤波器）的设计。先对数字滤波器进行简介。数字滤波器是数字信号处理的重要基础，在对信号的过滤、检测与参数的估计等信号处理中，数字滤波器是使用最为广泛的一种线性系统。数字滤波器是对数字信号实现滤波的线性是不变系统。设计数字滤波器包括以下几个步骤：
- 按照实际任务的要求，确定滤波器的性能指标。
- 用一个因果、稳定的离散线性时不变的系统函数去逼近这一性能指标。根据不同的要求可以用IIR系统函数，也可以用FIR系统函数去逼近。
- 利用有限精度算法实现系统函数，这里包括结构的选择、字长选择等。
设计一个滤波器，重要的是寻找一个稳定、因果的系统函数去逼近滤波器的技术指标。一个也能过、稳定的模拟滤波器的系统函数
应该满足如下条件：
- 滤波器的单位冲击响应函数
应该是一个实函数，即
是一个具有实系数的s的有理函数。
-
的极点必须分布在s平面的左半平面。
-
的分子多项式的阶数必须小于或者等于分母多项式的阶数。
实际上设计模拟原型滤波器是要寻求一个逼近理想低通滤波器的函数，所谓原型低通滤波器是指低通模拟或数字滤波器。模拟低通滤波器的设计方法有：巴特沃斯、切比雪夫和椭圆滤波器。模拟低通巴特沃斯滤波器是以巴特沃斯函数作为滤波器的系统函数，它的幅度平方函数表示为：
（1-16）
式中的N为正整数，表示滤波器的阶数。
为通带的截止频率，或
的3分贝带宽。模拟的低通巴特沃斯滤波器的设计过程包括以下两个过程：
- 按照给定的通带和阻带指标确定阶数N。
- 从幅度平方函数确定系统函数
。
切比雪夫
型滤波器在通带内幅度特性是等波纹的，在阻带是单调的，而切比雪夫
则相反，它在通带内是单调的，在阻带内是等波纹的。切比雪夫滤波器由3个参数需要确定：
、
和N，其中
是给定的截止频率，
由容许的通带波纹或通带的幅度误差
确定。椭圆滤波器是采样有限零点设计的滤波器，它能更好的逼近理想的低通滤波器。它的幅度平方函数为：
（1-17）
式中的
是雅可比椭圆函数，
是与通带衰减有关的函数，阶数N等于通带和阻带内最大点和最小点的总和。脉冲响应不变法的设计原理是使得数字滤波器的单位取样响应序列
模仿模拟滤波器的冲激响应
。双线性变换化是使得数字滤波器的频率响应模仿模拟滤波器的频率响应模仿模拟滤波器的频率响应的一种方法。这种方法的基本思路是：首先将整个s平面压缩到
平面的一条带宽为
（从
到
）的横带里，然后通过标准的变换关系
将横带变换成整个z平面上去，这便得到s平面与z平面之间的一一对应的单值关系。
- 对有限长单位脉冲响应数字滤波器进行介绍。窗函数在设计FIR数字滤波器中有很重要的作用，正确的选择窗函数可以提高所设计的数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶数。常见的窗函数有矩形窗、三角窗、布拉克曼窗、汉宁窗、海明窗、凯塞窗、巴特里特窗、切比雪夫窗8种。设计FIR数字滤波器最简单的方法是窗函数法，通常也称为傅里叶级数法。FIR滤波器的设计同IIR数字滤波器的设计一样，首先给出要求的理想滤波器的频率响应
，设计一个FIR数字滤波器频率响应
，去逼近理想的频率响应
。然而窗函数法设计FIR数字滤波器是在时域进行的，因而必须由理想的频率响应
推导出对应的单位取样响应
，在设计一个FIR数字滤波器的单位取样响应
去逼近
。频率取样法设计FIR数字滤波器是从频域出发，根据频域的采样定理，对给定的理想滤波器的频率响应
进行等间隔的采样：
（1-18）
等波纹切比雪夫法是采用最大误差最小准则得到最佳数字滤波器，并且最佳解是唯一的。切比雪夫逼近法设计FIR数字滤波器的过程是：
- 规定所需的频率响应
，加权函数
和滤波器的单位脉冲响应的长度N。
- 形成
、
和
。
- 利用雷米兹多重交换算法求解逼近问题。
- 计算滤波器的单位脉冲响应
。
采用切比雪夫逼近设计方法能够得到既有严格线性相位，又有很好衰减特性的滤波器，因此，切比雪夫逼近法在滤波器设计中占有很重要的位置。
- 讲述了功率谱估计的问题。功率谱估计方法可以分为经典谱估计和现代谱估法两种，经典谱估计法又可分为直接法和间接法，Bartlett法和Welch法是直接法的改进。随机序列
自相关函数估计的两种形式为：
（1-19）
用FFT计算自相关函数的一般步骤：
- 对
补N个零，得
，对
做DFT得
，k=0，1，…，2N-1；
- 求
的幅平方，然后除以N，得
；
- 对
做逆变换，的
。
在MATLAB中，函数xcorr用来进行自相关函数估计，且为局域上述FFT的快速算法，其格式为：C=XCORR(A,’flag’)，该函数返回长度为2N-1的自相关序列。AR模型功率谱估计是现代谱估计的主要内容。AR模型又称为自回归模型，它是一个全极点的模型，该模型现在的输出时现在的输入和过去输出的加权和。AR模型谱估计的性质有：AR谱的平滑特性、AR谱的分辨率。基于矩阵特征分解的功率谱估计包括特征向量估计与MUSIC估计，这两种估计方法均为非参数估计方法，特征向量估计主要使用混有白噪声的正弦信号的功率皮估计，而MUSIC估计适合更为普遍情况下正弦信号参数估计的方法。函数Pmusic为MUSIC估计，而函数Peig为特征向量估计。
2.1 图像复原的方法
在图像的获取、传输和保存的过程中，由于各种原因，如大气的湍流效应、传感器特性的非线性、成像设备与物体之间的相对运动、摄像设备中光学系统的衍射、胶片颗粒噪声和感光胶卷的非线性、光学系统的像差以及电视摄像扫描的非线性等所引起的几何失真，都可能造成图像的畸变和失真。通常，由于这些因素引起的质量下降称为图像退化[1]。
图像退化的表现是图像出现模糊、失真以及附加噪声等。由于图像的退化，在图像接收端显示的图像已经不再是发送的原始图像，图像的效果明显变差，因此我们可以采用一些技术手段来尽可能的减轻甚至消除图像质量的下降，还原图像的本来面目，这就是图像复原[2]。
图像复原是图像处理领域中一种非常重要的处理技术，与图像增强等其他基本图像处理技术类似，也是以获取视觉质量程度的改善为目的，所不一样的地方是图像复原过程实际上是一个估计的过程，需要根据一些特定的退化模型，对退化图像进行校正。简而言之，图像复原的处理过程就是提升退化图像的质量，并通过图像质量的提高来达到图像在视觉上的改善。
因为引起图像退化的原因很多，且性质各不相同，目前还没有统一的复原方法，许多研究人员根据不同的应用环境，采取了不同的退化模型[3-4]、估计准则和处理技巧，所以得到了各种复原方法。
图像复原的算法是整个技术的核心内容。目前，国内在这方面的研究才刚起步，而国外已经取得了较好的成果。早期的图像恢复是用光学的方法对失真的测试图像进行复原，但是关于数字图像复原技术的研究则是从对天文观测图像的后期处理中才逐渐发展起来的。其中一个成功例子是在1964年，NASA的喷气推进实验室用计算机处理有关月球的照片。照片是用电视摄像机在空间飞行器上拍摄的，图像的复原包括消除噪声和干扰等因素，校正几何失真和对比度损失以及反卷积。另一个典型的例子则是对肯尼迪遇刺事件现场照片的处理。由于事情发生的太突然，照片是在相机移动地过程中拍摄的。图像复原的主要目的就是消除移动造成的失真。还有其在在医学领域，图像复原的技术能广泛的应用于X光，CT，B超等成像系统，用来抑制各种医学成像系统或图像获取系统的噪声，改善医学图像的分辨率[5]。
现今的复原方法有：维纳滤波、逆滤波、最小二乘滤波、几何失真校正等多种复原方法。随着数信号处理和图像处理的发展，新的复原算法不断出现，不同的复原方法所需的条件是不相同的，所以在应用中可以根据具体情况加以选择。
目前国内外对于图像复原技术的研究和应用主要集中在空间探索、天文观测、物质研究、遥感遥测、军事科学、生物科学、医学影象、刑事侦察、交通监控等领域。如在生物方面，主要是用于生物活体细胞内部组织的三维再现和重构，通过复原荧光显微镜所收集的细胞内部的逐层切片图，来重现细胞内部构成；医学方面，如对肿瘤周围组织进行显微观察，以获取肿瘤安全切缘与癌肿原先部位之间关系的定量数据；天文方面则采用迭代盲反卷积进行气动光学效应图像复原研究等。
在图像的获取或者显示的过程中通常会产生几何失真，比如成像系统有一定的几何非线性，其主要原因是视像管摄像机和阴极射线管显示器的扫描偏转系统有一定的非线性，所以会造成图像的枕形失真或桶性失真。另外，由卫星拍摄的地球表面的图像往将往覆盖较大的面积，由于地球表面呈球形，这样拍摄的图像同样会有比较大的几何失真。
几何失真一般分为系统失真和非线性失真。系统失真是有规律的、能预测的；而非线性系统的失真则是随机的。当对图像进行定量分析的时候，就要先对畸变的图像进行准确的几何失真校正[12]（即把几何失真的图像校正成与原图像相差无几的图像），以免影响分析精度。基本的方法是：首先，要建立几何校正的数学模型；其次，用已知的条件来确定模型参数；最后，根据畸变模型对图像进行几何校正。通常可以分为以下两步：
（1）图像空间坐标的变换；
（2）确定校正空间各像素的灰度值（灰度内插）。
3.1 空间变换
，经过几何失真变成了
，其中，
表示的是失真图像的坐标，而不是原图像的坐标
。上述变化可表示为：
（3-1）
（3-2）
和
是空间变换，产生了几何失真图像
。
和
已知，则可以依据一个坐标系统的像素坐标算出另一坐标系统的对应像素的坐标。在已知情况下，通常
和
可用多项式来近似：
（3-3）
（3-4）
式中，
为多项式的次数，
和
为各项系数。
3.1.1 已知
和
条件下的几何校正
若我们具备先验知识
、
，则希望将几何畸变图像
恢复为基准几何坐标的图像
。几何校正通常分为直接和间接两种方法。
- 直接法。先通过
来推出
，然后计算每个像素的校正坐标值，保持各像素灰度值不变，然后生成一幅校正图像，但其像素的分布是不太规则的，可能会出现疏密不均、像素挤压等现象。
- 间接法。假设复原图像的像素在基准坐标系统内是等距网格的交叉点，从网格交叉点的坐标
出发，可以算出在已知几何畸变图像上的坐标
，即：
（3-5）
虽然点
坐标为整数，但
一般不为整数，不会刚好就处于畸变图像像素的中心，所以不能直接明确该点的灰度值，而只能由其在几何失真图像的周围像素灰度内插求出，作为对应像素
的灰度值，据此可以获得校正之后的图像。由于间接法内插灰度比较容易，所以在通常情况下采用的是间接法来进行几何失真校正。
3.1.2
和
未知条件下的几何失真
在这种情况下，通常用原始图像和几何畸变图像上多对“连接点”的坐标来确定
和
。
假设基准图像像素的空间坐标
和待校正图像的对应像素的空间坐标
之间的关系用二元多项式来表示，表达式如下：
（3-6）
（3-7）
式中，
为多项式的次数，
和
为各项待定系数。
对于线性失真：
（3-8）
（3-9）
对于一般的（非线性）二次失真：
（3-10）
（3-11）
利用“连接点”建立失真图像与校正图像之间其他像素空间位置的相应关系，而“连接点”在失真图像和校正图像中的位置是精确已知的。失真图像和校正图像中有四边形区域，这两个四边形的顶点就是相应的“连接点”。假定四边形区域中的几何畸变的过程可以用二次失真方程来表示，即：
（3-12）
（3-13）
将以上两个表达式代入式（3-1）和式（3-2）中，得：
（3-14）
（3-15）
因为一共有四对“连接点”，代入式（3-14）和式（3-15）可得8个联立方程，由这些方程可以解出8个系数
、
、
、
、
、
、
、
。这些系数就可以形成用于变换四边形区域内所有像素的几何失真模型，即空间映射公式。通常来说，可以把一幅图像分割成一系列覆盖全图的四边形区域的集合，在每个区域寻找足够的“连接点”用来计算进行映射所需的系数。
只要有了系数，校正（即复原）图像就不那么困难了。如果想找到非线性失真图像在任一点的
的值，需要简单地知道
在失真图像中的什么地方被映射。为此，可以把
代入式（3-14）和式（3-15）得到几何失真坐标
。在无失真图像中被映射到
点的值是
。这样简单地令
，就得到了复原图像的值。
3.2 灰度插值
最简单的灰度灰度插值是最近邻插值（也称为零阶插值），该方法实现起来相对简单，但是有时不够精确，甚至经常产生不希望的认为疵点，如高分辨率图像直边的扭曲；对于通常的图像处理，双线性插值很实用；更完善的技术如样条插值、立方卷积内插等可以得到较平滑的结果，但是更平滑的仿真的代价就是增加计算开销。
在MATLAB[13]软件中生成如图3.1所示的figure图像，然后在figure图像上的失真图像和原始图像上选取九对点，最后依据选取的九对点来进行校正。
图3.1 选取连接点图像
图3.2 几何失真及校正后图像
图3.1是调用cpselect函数，使系统启动交互连接点工具，然后在figure文件上的图像交替选择对应的坐标点，将其输出到Workspace中，进行坐标点的统计；图3.2是校正的结果，第二张是几何失真之后的图像，第三张是几何校正后的图像，校正后的图像和原图像相比，几乎没有什么差别，所以几何失真校正复原出来的结果很好，但是它的适用条件是图像的几何失真，其他的失真方式，比如说图像的质量退化，就不可以采用这种复原方法。
[1] 朱冠南.基于MATLAB的图像复原设计[J].技术交流,2009:87.
[2] 王婷.退化图像的复原改进算法研究与实现[D].哈尔滨:哈尔滨工程大学,2007.
[3] 何东健.数字图像处理[M].西安:西安电子科技大学出版社,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.
几何失真校正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);
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.7500;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 424.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('几何失真校正');
901
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
