基于MATLA的图像复原系统

基于MATLA的图像复原系统

二、图像复原的基本概念

图像在形成、传输和记录过程中，由于受到多种原因的影响，图像的质量就会有所下降，典型的表现为图像模糊、失真、有噪声等，这一过程称为图像的退化[5]。

图像复原是试图利用退化过程的先验知识使已退化的图像恢复本来面目，即根据退化的原因，分析引起退化的环境因素，建立相应的数学模型，并沿着使图像降质的逆过程恢复图像[5]。目的在于消除或减轻在图像获取以及传输的过程中造成的图像品质下降，恢复图像的本来面目。因此，复原技术就是把退化模型化，并采用相反的过程进行处理，以便尽可能复原被退化图像的本来面目。

广义上讲，图像复原是一个求逆问题，逆问题经常存在非唯一解，甚至无解。要想恢复全真的景物图像比较困难。为了得到逆问题的有用解，图像复原本身往往需要一个质量标准，即衡量接近全真景物图像的程度，或者说，对图像的估计是否达到最佳的程度。需要有先验知识以及对解的附加约束条件。

典型的图像复原是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用各种逆退化处理的方法进行恢复，使图像质量得到改善。

2.2图像退化的数学模型

2.3逆滤波复原

逆滤波复原法也叫做反向滤波法，其主要过程是首先将要处理的数字图像从空间域转换到傅里叶频域中，进行反向滤波后再由频率域转回到空间域，从而得到复原的图像信号[5]。

但实际用逆滤波存在病态的情况：

噪声存在，当H(u,v)很小或为零时，则噪声被放大。这意味着退化图像中小噪声的干扰在H(u,v)较小时，会对逆滤波恢复的图像产生很大的影响，有可能使恢复的图像和f(x,y)相差很大，甚至面目全非。

实验证明，当退化图像的噪声较小，即轻度降质时，采用逆滤波复原的方法可以获得较好的结果。通常，在离频率平面原点较远的地方数值较小或为零，因此图像复原在原点周围的有限区域内进行，即将退化图像的傅里叶频谱限制在没出零点而且数值又不是太小的有限范围内。

2.4 维纳滤波复原

逆滤波比较简单，但没有清楚地说明如何处理噪声，而维纳滤波综合了退化函数和噪声统计特性两个方面进行复原处理。维纳滤波是维纳在1949年提出的，并应用于一维平稳时间序列，获得了满意的结果。这是最早也是最著名的线性滤波技术。

采用维纳滤波是假设图像信号可以近似看成平稳随机过程的前提下，按照使f(x,y)和之间的均方误差达到最小的准则函数来实现图像复原的，即

式中，E（）代表求期望值。因此维纳滤波又称为最小均方误差滤波器。

维纳滤波需要假定下述条件成立：

系统为线性空间移不变系统。
退化图像、原始图像、噪声都是均匀随机场，噪声的均值为零，且与图像不相关。

维纳滤波的复原滤波函数，即滤波器的传递函数为：

没有噪声时，维纳滤波退化为逆滤波。有噪声时，维纳滤波利用信噪功率比对恢复过程进行修正，在信噪功率比很小的区域内，P(u,v)的值也很小，这使恢复图像较小地依赖于退化图像。在H(u,v)很小或等于零时，P(u,v)的分母不为零，维纳滤波没有病态问题[7]。

在实际系统中，维纳滤波经常用下式近似：

①的维纳滤波要求未退化图像和噪声的功率必须是已知的。虽然用②近似的方法能得到好的结果，但功率谱比常数K的估计一般没有合适的解[8]。

MATLAB语言程序

clear;

I=imread('C:\ok\原始图.jpg');

imshow(I);

I=rgb2gray(I); %将原图像转化为黑白图

figure;

subplot(2,2,1);

imshow(I);

title('转成黑白图像');

[m,n]=size(I);

F=fftshift(fft2(I));

k=0.0025;

for u=1:m

for v=1:n

H(u,v)=exp((-k)*(((u-m/2)^2+(v-n/2)^2)^(5/6)));

end

G=F.*H;

I0=real(ifft2(fftshift(G)));

I1=imnoise(uint8(I0),'gaussian',0,0.001)

subplot(2,2,2);

imshow(uint8(I1));

title('模糊退化且添加高斯噪声的图像');

F0=fftshift(fft2(I1));

F1=F0./H;

I2=ifft2(fftshift(F1));

subplot(2,2,3);

imshow(uint8(I2));

title('全逆滤波复原图');

K=0.1;

for u=1:m

for v=1:n

H(u,v)=exp(-k*(((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));

end

F2=H1.*F0;

I3=ifft2(fftshift(F2));

subplot(2,2,4);

imshow(uint8(I3));

title(维纳滤波复原图);

运行结果如下：

原始图：

复原后图像：

经过仿真，如上图所示，可以看出逆滤波复原与维纳滤波复原的区别和联系。维纳滤波后虽然仍有一些噪声存在，但已经和原图很接近了。因为原图像和噪声函数都是已知的，可以正确的估算参量。

2.5使用Lucy-Richardson算法的迭代非线性复原

L-R算法程序：

I=imread('C:\ok\苹果.jpg');

PSF=fspecial('gaussian',5,5) ;

BN=imnoise(Blurred,'gaussian',0,V);

luc=deconvlucy(BN,PSF,5);

figure

subplot(2,2,1);

imshow(I);

title('原始图像');

subplot(2,2,2);

imshow (Blurred);

title('模糊后的图像');

subplot(2,2,3);

imshow (BN);

title('加噪后的图像');

subplot(2,2,4);

imshow (luc);

title('恢复后的图像');

模拟实验结果如下：

用Lucy-Richardson算法可以较好的恢复图像[1]。但由于迭代产生的噪声痕迹是最大化可能性数据逼近法的常见问题，在低信噪比条件下，恢复图像可能会出现一些斑点，这些斑点并不代表图像的真实结构，只不过是恢复图像过于逼近噪声所产生的结果。另外此法存在一些较严重的缺陷，一是噪声放大问题，二是对于恢复图像中的不同部分，分别执行多少迭代才合适的问题。因为图像中信噪比高的部分可能需要数百次迭代才能获得满意的结果；而另一些光滑的对象可能只需很少次数即可达到满意的结果，所以适当选择迭代次数对图像恢复也很重要。这两个问题若得不到解决，将会对最终结果产生不利影响。

2.6 盲去卷积

通常图像恢复方法均在成像系统的点扩展函数PSF已知下进行, 实际上它通常是未知的. 在 PSF未知的情况下, 盲去卷积是实现图像恢复的有效方法。因此，盲去卷积算法就是那些不以PSF知识为基础的图像复原的方法。

在过去的20年里，一种盲去卷积的方法已经受到了人们的极大重视，它是以最大似然估计(MLE)为基础，即一种用被随机噪声所干扰的量进行估计的最优化策略。简要的说，关于MLE方法的一种解释就是将图像数据看成随机量，它们与另外一族可能的随机量之间有着某种似然性。似然函数用g(x,y),f(x,y)和h(x,y)加以表达，然后，问题就变成了寻求最大似然函数。在盲去卷积中，最优化问题规定的约束条件并假定收敛时通过迭代来求解，得到的最大f(x,y)和h(x,y)就是还原的图像和PSF[3]。

工具箱通过函数deconvblind来执行盲去卷积，它有如下语法：

[f,PSFe]=deconvblind(g,INITPSF)

其中，g代表退化函数，INITPSF是点扩散函数的出事估计。PSFe是这个函数最终计算到的估计值，fr是利用估计的PSF复原的图像。用来去的复原图像的算法是L-R迭代复原算法[13]。PSF估计受其初始推测尺寸的巨大影响，而很少受其值的影响。

盲去卷积程序：

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');

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,[]);

模拟实验结果如下：

该算法优点是,同时恢复了图像和点扩张函数,在对失真情况毫无先验知识的情况下, 仍能实现对模糊图像的恢复操作。利用 MATLAB实现的图像恢复, 并对恢复图像的失真情况做了改善。在进行图像恢复时,重建 PSF,对图像进行重建, 得到恢复的图像。

三、图像重建

3.1 图像重建的概述

图像处理中的一个重要研究分支是物体图像的重建，它被广泛的应用于检测和观察中，而这种重建方法一般是根据物体的一些横截面部分的投影而进行的。在一些应用中，某个物体的内部结构图像的检测只能通过这种重建才不会有任何物理上的损伤。由于这种无损检测技术的显著优点，因此，它的适用面非常广泛，它在各个不同的领域都显示出独特的重要性。例如：医疗放射学、核医学、电子显微、无线和雷达天文学、光显微和全息成像学及理论视觉等等领域都有应用。在医学影像处理中是医学图像获取的重要方法。如医疗放射学、核医学，电子显微等领域是必不可少的技术，在工业生产中的无损检测技术图像重建也扮演重要的角色。

图像重建经过多年研究已取得巨大进展，产生了许多有效的算法，如：傅里叶反投影法、卷积反投影法、代数法、迭代法等，其中以卷积反投影法运用最为广泛。近年来，由于与计算机图形学相结合，把多个二维图像合成三维图像，并加以光照模型和各种渲染技术，已能生成各种具有强烈真实感的高质量三维人工合成图像。

3.2傅里叶反投影重建

傅里叶反投影重建方法可以说是最简单的一种变换重建方法，一个三维（或二维）物体，它的二维（或一维）投影的傅里叶变换恰好与此物体的傅里叶变换的主题部分相等，傅里叶变换的重建方法正是以此为基础的。该重建方法最早于1974年由Shepp和Logan提出，该方法是建立在“投影切片定理”这一理论基础之上的。根据此结论，可以先对投影进行旋转和傅里叶变换以构造整个傅里叶变换域中的各个方向的切片数据，然后再对其进行傅里叶反变换即可得到重建后的目标（原空间域中的图像）。

傅里叶变换重建的原理如下：

若已知无数的投影，从极坐标中计算得到的投影变换推出在矩形平面F(u,v)中的傅里叶变换并不困难。但是，若只有有限个投影是有效的，则可能需要在变换中插入一些数据。另外需要注意的是，虽然只需一维傅里叶变换的投影数据就可构成变换空间，但图像重建则需要二维反变换。由此，我们得出一个结论，即：三维图像不能在得到部分投影数据的过程中局部的重建，而必须延迟到所有投影数据都获得之后才能重建。

3.3 卷积法重建

卷积法重建的基本思路和方法：

逆投影原理：从各个方向得到的投影逆向返回到该方向的各个位置，如果对多个投影方向中的每个方向都进行这样的逆投影，就可能建立平面上的一个部分。典型的方法是卷积逆投影重建。

卷积重建法是一种变换重建法，可以根据傅里叶变换投影定理推出。

按照二维傅里叶反变换标准定义，有

MATLAB程序：

P=phantom(256);

subplot(2,2,1);

imshow(P);

title('大脑的幻影图');

theta1=0:10:170;

[R1,xp]=radon(P,theta1);

theta2=0:5:175;

[R2,xp]=radon(P,theta2);

theta3=0:2:178;

[R3,xp]=radon(P,theta3);

figure;

imagesc(theta3,xp,R3);

title('大脑幻影图的90条投影光速的Radon变换');

colormap(hot);

colorbar;

xlabel('\theta');

ylabel('x\prime');

I1=iradon(R1,10);%R1有18条投影光速

I2=iradon(R2,5);% R2有36条投影光速

I3=iradon(R3,2);%R3有90条投影光速

figure(1);

subplot(2,2,2);

imshow(I1);

title('用R1重建图像');

subplot(2,2,3);

imshow(I2);

title('用R2重建图像');

subplot(2,2,4);

imshow(I3);

title('用R3重建图像');

通过MATLAB演示我们可以看出：用卷积逆投影来实现，数据质量高的情况下可重建出准确清晰的图像，而用傅里叶反变换重建法实现，不是很容易，且重建的图像质量很差。

3.4代数重建方法

代数重建技术就是事先对未知图像的各像素给予一个初始估值，然后利用这些假设数据去计算各射线穿过对象时可能得到的投影值（射影和），再用它们和实测投影值进行比较，根据差异获得一个修正值，利用这些修正值，修正各对应射线穿过的诸像素值。如此反复迭代，直到计算值和实测值接近到要求的精确度为止[5]。

具体实施步骤：

（1）、对于未知图像各像素均给予一个假定的初始值，从而得到一组初始计算图像；

（2）、根据假设图像，计算对应各射线穿过时，应得到的各个相应投影值，，……，；

（3）、将计算值和对应的实测值进行比较，然后取对应差值作为修正值；

（4）、用每条射线的修正值修正和该射线相交的诸像素值；

（5）、用修正后的像素值重复1~4各步，直到计算值和实测值之差，即修正值小到所期望的值为止。

只要所测的射线投影值组成一个独立的集合，那么代数重建便将收敛于唯一解。

参考文献

[1] 龚声蓉，刘纯平，王强.数字图像处理[M].北京:清华大学出版社,2006.46-84,123-144[2] 张德丰，雷晓平.MATLAB基础与工程应用.北京：清华大学出版社，2012.

[3] 冈萨雷斯等.数字图像处理（MATLAB版）[M].北京:电子工业出版社,2005.5:123-133

[4] 张铮,倪红霞,苑春苗，杨立红.精通MATLAB数字图像处理与识别[M].北京:人民邮电出版社,2013.8:194-196

[5] 陈传波,金先级.数字图像处理[M].北京：机械工业出版社,2004.119-124,161-176

[6] 林旭梅，葛广英.MATLAB实用教程[M].山东：中国石油大学出版社，2010.3:49-75

[7] 张德丰,张葡青.维纳滤波图像恢复的理论分析与实现[J].中山大学学报,2006.45(6)：44-47

[8] 何东健.数字图像处理[M].西安:西安电子科技大学出版社,2003.261-279

[9] 徐昕.MATLAB工具箱应用指南：控制工程篇[M].北京:电子工业出版社,2000.1-12

[10] 罗军辉，冯平.MATLAB7.0在图像处理中的应用[M].北京:机械工业出版社,2005.257-277

[11] 赵怡红,张常年编著.数字信号处理及其MATLAB实现.北京:化学工业出版社,2002.1:38

-125

[12]徐晰.MATLAB工具箱应用指南：控制工程篇[M].北京：电子工业出版社，2000.1-12

[13]刘维一，于德月等.用迭代法消除数字图像放大后的模糊[J].光电子.激光，2002,13（4）：389-400

[14]陶洪.数字共焦显微技术及其图像复原算法研究四川大学硕士学位论文CNKI::CDMD:10610.2.2003.6632

[15]苏开亮.运动模糊图像的恢复及恢复质量评价[D].西安：西安电子科技大学出版社，2010,21-25