MATLAB在图像复原中的应用研究
摘 要：图像复原是图象处理的一个重要课题。图像复原也称图象恢复，是图象处理中的一大类技术。它的主要目的是改善给定的图像质量。当给定了一幅退化了的或者受到噪声污染了的图像后，利用退化现象的某种先验知识来重建或恢复原有图像是复原处理的基本过程。可能的退化有光学系统中的衍射，传感器非线性畸变，光学系统的像差，摄影胶片的非线性，大气湍流的扰动效应，图像运动造成的模糊及几何畸变等等。噪声干扰可以由电子成像系统传感器、信号传输 过程或者胶片颗粒性造成。各种退化图像的复原都可归结为一种过程，具体地说就是把退化模型化，并且采用相反的过程进行处理，以便恢复出原图像。文章介绍了图象退化的原因，几种常用的图像滤波复原技术，以及用MATLAB实现图像复原的方法。
关键词：退化模型 ；噪声干扰；图像滤波；图像复原
1．图像复原的概念
1.1图像复原的定义
图像复原也称图象恢复，是图象处理中的一大类技术。所谓图像复原，是指去除或减轻在获取数字图像过程中发生的图像质量下降（退化）这些退化包括由光学系统、运动等等造成图像的模糊，以及源自电路和光度学因素的噪声。 图像复原的目标是对退化的图像进行处理，使它趋向于复原成没有退化的理想图像。成像过程的每一个环节（透镜，感光片，数字化等等）都会引起退化。在进行图像复原时，既可以用连续数学，也可以用离散数学进行处理。其次，处理既可在空间域，也可在频域进行。
1.2图象恢复与图象增强的异同
相同点：改进输入图像的视觉质量 。
不同点：图象增强目的是取得较好的视觉结果(不考虑退化原因)； 图象恢复根据相应的退化模型和知识重建或恢复原始的图像(考虑退化原因)。
1.3 图象退化的原因
图象退化指由场景得到的图像没能完全地反映场景的真实内容，产生了失真等问题。其原因是多方面的。如：
透镜象差/色差
聚焦不准（失焦，限制了图像锐度）
模糊（限制频谱宽度）
噪声（是一个统计过程）
抖动（机械、电子）
1.4图象退化举例
如图1所示是两个图象退化的例子。
图1 退化图像与原始图像
2．退化模型
2.1图象退化模型概述
图像复原处理的关键问题在于建立退化模型。在用数学方法描述图像时，它的最普遍的数学表达式为
这样一个表达式可以代表一幅活动的、彩色的立体图像。当研究的是静止的、单色的、平面的图像时，则其数学表达式就简化为
基于这样的数学表达式，可建立如图2所示的退化模型。由图2的模型可见，一幅纯净的图像
是由于通过了一个系统H及加性噪声
而使其退化为一幅图像
的。
图2 图像退化模型
图像复原可以看成是一个估计过程。如果已经给出了退化图像
并估计出系统参数H，从而可近似地恢复
。这里，
是一种统计性质的噪声信息。当然，为了对处理结果做出某种最佳的估计，一般应首先明确一个质量标准。根据图像的退化模型及复原的基本过程可见，复原处理的关键在于对系统H的基本了解。就一般而言，系统是某些元件或部件以某种方式构造而成的整体。退化模型可分为连续函数退化模型和离散函数退化模型。
2.2连续函数退化模型
假定系统H对坐标为(,)处的冲激函数(x-,y-)的冲激响应为h(x,,y,)，则
此式说明，如果系统H对冲激函数的响应为已知，则对任意输入的响应可用上式求得，即，线性系统H完全可以由冲激响应来表征。图像中冲激响应也称为点扩散函数。
在有噪音的情况下：
2.3离散函数退化模型
对和进行均匀取样后，就可引伸出离散函数的退化模型。用一维的来说明。如果f (x)和h(x)周期分别A和B的序列，为避免卷积周期重叠需要对它们进行周期扩展为周期为M ≥ A + B – 1。
f(x)0 ≤ x ≤ A-1h(x)0≤ x≤ B-1
fe(x)=he(x)=
0A-1≤ x ≤ M-10 B-1< x≤ M-1
那么它们的时域离散卷积可定义为下式：
显然，上式也是具有周期M的序列。
如果用矩阵来表示上述离散退化模型，可写成下式之形式:
退化过程为：
图像f(x,y)被线性操作h(x,y)所模糊，并叠加上噪声n(x,y)，构成了退化后的图像g(x,y)。退化后的图像与复原滤波器卷积得到复原的f(x,y)图像。
图3 图像的退化/复原过程模型
3.图象复原技术
3.1无约束恢复
由退化模型得：
最小均方误差准：
在最小二乘方意义上说，希望找到一个f使下式的值最小:
3.2逆滤波
设M = N，则：
退化函数H (u, v)与F (u, v)相乘为退化过程，用H (u, v)去除G (u, v) 是复原过程，称其为逆滤波。可描述为：
记M (u, v)为复原转移函数，则其等于1 / H (u, v).
3.3 维纳（Wiener）滤波器
它一种最小均方误差滤波器。
设 Rf 是 f 的相关矩阵：
Rf的第 ij 元素是E{fi fj}，代表 f 的第 i 和第 j 元素的相关。
设 Rn是n 的相关矩阵：
根据两个象素间的相关只是它们相互距离而不是位置的函数的假设，可将Rf和Rn都用块循环矩阵表达，并借助矩阵W来对角化：
fe(x, y)的功率谱，记为Sf(u, v) ；ne(x, y)的功率谱，记为Sn(u, v)。D是1个对角矩阵，D(k, k) = (k)，则有：
定义：
代入：
两边同乘以W –1，有：
最后整理得：
3.4图像复原例图
以下的几幅图是用MATLAB软件根据不同的复原方法进行的图像复原。根据图4例图可看出不同复原方法的区别。
原图退化图像
全逆滤波 半径受限逆滤波 维纳滤波结果
图4复原例图
4．图像复原的MATLAB实现实例
维纳滤波复原规则化滤波复原
Lucy-Richardson复原盲目去卷积复原
图5 图像复原实例
5.结束语
本文简要介绍了图像退化的原因，图像退化的模型，图像复原的概念，几种常用的图像复原的方法，以及利用MATLAB实现图像复原的几个例子。简单的讲述了MATLAB在图像复原中的应用。
参考文献：
[1]阮秋琦编著.—2版。：电子工业，2007.2.
[2](美)卡斯尔曼(castleman，k．R)著；X志刚等译.数字图像处理.：电子工业，2002.2.
[3] 孙家广等主编.计算机图形学.第3版.：清华大学，1998.9.
[4] 罗军辉等主编.MATLAB7.0在图像处理中的应用.第1版.：机械工业，2007.7.
附录：
（1）.维纳滤波复原源代码：
I=checkerboard(8); noise=0.1*randn(size(I));
PSF=fspecial('motion',21,11);
Blurred=imfilter(I,PSF,'circular');
BlurredNoisy=im2uint8(Blurred+noise);
NP=abs(fftn(noise)).^2;
NPOW=sum(NP(:)/numel(noise));
NCORR=fftshift(real(ifftn(NP)));
IP=abs(fftn(I)).^2;
IPOW=sum(IP(:)/numel(noise));
ICORR=fftshift(real(ifftn(IP)));
ICORR1=ICORR(:,ceil(size(I,1)/2));
NSR=NPOW/IPOW;
subplot(221);imshow(BlurredNoisy,[]);
title('模糊和噪声图像');
subplot(222);imshow(deconvwnr(BlurredNoisy,PSF,NSR),[]);
title('deconbwnr(A,PSF,NSR)');
subplot(223);imshow(deconvwnr(BlurredNoisy,PSF,NCORR,ICORR),[]);
title('deconbwnr(A,PSF,NCORR,ICORR)');
subplot(224);imshow(deconvwnr(BlurredNoisy,PSF,NPOW,ICORR1),[]);
title('deconbwnr(A,PSF,NPOW,ICORR_1_D)');
（2）.规则化滤波复原程序源代码：
I=checkerboard(8);
PSF=fspecial('gaussian',7,10);
V=.01;
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);
NOISEPOWER=V*numel(I);
subplot(221);imshow(BlurredNoisy);
title('A=Blurred and Noisy');
subplot(222);imshow(J);
title('[J LAGRA]=deconvreg(A,PSF,NP)');
subplot(223);imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA/10));
title('deconvreg(A,PSF,[],0.1*LAGRA)');
subplot(224);imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA/10))
title('deconvreg(A,PSF,[],10*LAGRA');
（3）.Lucy-Richardson复原滤波源代码：
I=checkerboard(8);
PSF=fspecial('gaussian',7,10);
V=.0001;
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);
WT(5:end-4,5:end-4)=1;
J1=deconvlucy(BlurredNoisy,PSF);
J2=deconvlucy(BlurredNoisy,PSF,20,sqrt(V));
J3=deconvlucy(BlurredNoisy,PSF,20,sqrt(V),[],WT);
subplot(221);imshow(BlurredNoisy);
title('A=Blurred and Noisy');
subplot(222);imshow(J1);
title('deconvlucy(A,PSF)');
subplot(223);imshow(J2);
title('deconvlucy(A,PSF,NI,DP)');
subplot(224);imshow(J3);
title('deconvlucy(A,PSF,NI,DP,[],WT)');
（4）.盲目去卷积复原源代码：
I=checkerboard(8);
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(221);imshow(BlurredNoisy);title('A=Blurred and Noisy');
subplot(222);imshow(PSF,[]);title('True PSF');
subplot(223);imshow(J);title('Deblured Image');
subplot(224);imshow(P,[]);
