MATLAB图像复原系统

MATLAB图像复原系统

一、研究背景

图像复原是数字图像处理中的一个重要课题。它的主要目的是改善给定的图像质量并尽可能恢复原图像。图像在形成、传输和记录过程中，受多种因素的影响，图像的质量都会有不同程度的下降，典型的表现有图像模糊、失真、有噪声等，这一质量下降的过程称为图像的退化。图像复原的目的就是尽可能恢复被退化图像的本来面目。

在成像系统中，引起图像退化的原因很多。例如，成像系统的散焦，成像设备与物体的相对运动，成像器材的固有缺陷以及外部干扰等。成像目标物体的运动，在摄像后所形成的运动模糊。当人们拍摄照片时，由于手持照相机的抖动，结果像片上的景物是一个模糊的图像。由于成像系统的光散射而导致图像的模糊。又如传感器特性的非线性，光学系统的像差，以致在成像后与原来景物发生了不一致的现象，称为畸变。再加上多种环境因素，在成像后造成噪声干扰。人类的视觉系统对于噪声的敏感程度要高于听觉系统，在声音传播中的噪声虽然降低了质量，但时常是感觉不到的。但景物图像的噪声即使很小都很容易被敏锐的视觉系统所感知。图像复原的过程就是为了还原图像的本来面目，即由退化了的图像恢复到能够真实反映景物的图像。

在交通系统、刑事取证中图像的关键信息至关重要，但是在交通、公安、银行、医学、工业监视、军事侦察和日常生活中常常由于摄像设备的光学系统的失真、调焦不准或相对运动等造成图像的模糊，使得信息的提取变得困难。但是相对于散焦模糊，运动模糊图像的复原在日常生活中更为普遍，比如高速运动的违规车辆的车牌辨识，快速运动的人群中识别出嫌疑人、公安刑事影像资料中提取证明或进行技术鉴定等等，这些日常生活中的重要应用都需要通过运动模糊图像复原技术来尽可能地去除失真，恢复图像的原来面目。因此对于运动模糊图像的复原技术研究更具有重要的现实意义。

二、匀速直线运动模糊图像的复原技术

运动模糊现象在数字图像处理实践中经常遇到，在照片曝光期间相机与景物之间的相对运动形成了运动模糊，其中以匀速直线运动最为常见。对于运动模糊图像的复原方法研究非常具有现实意义。因为运动模糊图像在日常生活中普遍存在，给人们的实际生活带来了很多不便甚至危及到安全保障体系。一个典型的例子就是现在很多城市的一些重要交通路口都设置了交通监视系统，它能及时拍摄车辆图像并从图片中分析出该车车牌号。由于车辆在行进中速度较快，所以摄取的画面有时是模糊不清的，这就需要运用运动模糊图像复原技术来进行图像复原，来得到可辨认的车牌图像。运动模糊图像复原技术在银行监视系统中识别经济犯罪、在路况监视系统中监控超速行驶、在刑事侦破中提供线索等方面也发挥着重要的作用。

2.1 匀速直线运动模糊图像的退化模型

在所有的运动模糊中，由匀速直线运动造成的模糊图像的复原问题更具有一般性和普遍意义。因为变速的、非直线运动在某些条件下可以被分解为分段匀速直线运动。

将退化的点扩散函数记作

，在不考虑噪声干扰的情况下，由目标与摄像机相对运动造成的图像模糊模型可以用图3-1来描述：

图2-1匀速直线运动模糊图像的退化模型

得到退化图像如图(b)。

（a）车牌图像（b）退化图像

图2-2 原图像和退化图像

图(c)和(d)分别为原图像和模糊图像的二次傅里叶变化后的图像:

图2-3 二次傅里叶变换后的原图像及退化图像

利用图(d)粗略的计算模糊的方向，可以通过matlab自带的画线工具，选取如下图

所示的三角形，计算a与c之间的夹角。

图2-4模糊方向的获取

Matlab命令窗口输入：ginput，选中三个顶角后回车，可得：

ans =

136.6009 146.3977

137.4205 166.8895

145.2074 166.8895

粗略取值后，经Matlab计算得：

atan(20/9)*180/pi

ans = 65.7723

则运动方向为90-ans≈25（误差为5）。

四、匀速直线运动模糊图像的几种复原方法

图像恢复也称图像复原是指根据相应的退化模型和先验知识，把品质下降了的图像加以重建．从退化图像中重构出原始图像。匀速直线运动模糊图像恢复技术的研究具有一般性和普遍意义，因为非匀速直线运动在某些条件下可以近似地视为匀速直线运动或者可以看作是多个匀速直线运动的合成。而在匀速直线运动模糊的所有模型中，水平方向的模糊更具有代表性和简单性，因为其它方向上的模糊可以由水平方向上的恢复方法推广得到，或者直接将图像旋转就可以把其它方向的运动模糊转化为水平方向上的匀速直线运动模糊。

对于退化的复原，一般采用两种方法，一种方法适用于对图像缺乏先验知识的情况，此时可对退化过程(模糊和噪声)建立模型，进行描述，进而寻找一种去除或削弱其影响的过程。另一方面，若对于原始图像有足够的先验知识，则对原始图像建立一个数学模型，并根据它对退化图像进行恢复会更有效。本文采用的是已知退化模型和噪声的情况下，对图像复原的常用方法进行研究。

下面介绍逆滤波、有约束的最小二乘方滤波、维纳滤波和Richardson-Lucy迭代算法进行图像复原，并进行讨论分析。

4.1 逆滤波复原

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

如果退化图像为

，原始图像为

，在不考虑噪声的情况下，其退化模型用式（4-1）表示，现在将其重写如下：

(4-1)

上式两边进行傅里叶变换得式（4-2）

(4-2)

式中

，分别是退化图像

，点扩散函数

，原始图像

的傅里叶变换。

由式子（4-2）以及傅里叶逆变换公式可得式（4-3）和式（4-4）

(4-3)

(4-4)

式中，

可以理解为成像系统的“滤波”传递函数。在频域中系统的传递函数与原图像信号相乘实现“正向滤波”，这里

除以

起到了“反向滤波”的作用。这意味着，如果已知退化图像的傅里叶变换和“滤波”传递函数，则可以求得原始图像的傅里叶变换，经反傅里叶变换就可求原始图像

。这就是逆滤波复原的基本原理。前面为了分析问题的简化，没有考虑噪声的影响，在有噪声的情况下，逆滤波复原的基本原理可以写成如（4-5）和（4-6）形式：

(4-5)

(4-6)

式中，

是噪声

。

由于在逆滤波复原公式（4-4）中，

处于分母的位置上，利用式（4-4）和（4-6）进行图像复原处理可能会发生下列情况：即在

平面上有些点或域会产

=0或

非常小的情况，在这种情况，即使没有噪声，也无法精确的恢复

。另外，在有噪声存在时，在

的领域内,

的值可能比

的值小得多，因此有式（4-6）得到的噪声项可能会非常大，这样也会使

不能正确恢复。一般来说，逆滤波复原法不能正确地估计

的零点，因此必须采用一个折中的方法进行解决。实际上，逆滤波不是用

，而是采用另外一个关于

的函数

。它的处理框图如图4-1所示：

图4-1 实际的逆滤波处理框图

在没有零点并且也不存在噪声的情况，有式（4-7）：

(4-7)

图4-1的模型包括了退化和恢复运算。退化和恢复总的传递函数可以用

，

来表示。此时有式（4-8）

(4-8)

一般情况，可以将图像的退化过程视为一个具有一定带宽的带通滤波器，随着频率的升高，该滤波器的带通特性很快下降，即

的幅度随着

平面原点的距离的增加而迅速下降，而噪声项

的幅度变化是比较平缓的。在远离

平面的原点

的值就会变得很大，而对于大多数图像来说，

却变小，在这种情况下，噪声反而占优势，自然无法满意的恢复出原图像。这一规律说明，应用逆滤波时仅在原点领域采用

方能有效。换句话说，应使

在下述范围内选择式（4-9）

(4-9)

式中，

的选择应该将

的零点排除在此领域之外。

4.2 有约束最小二乘方复原

为了在数学上更容易处理，通常在无约束复原的基础上附加一定的约束条件，从而在多个可能结果选择一个最佳结果，这便是有约束复原方法。

最小二乘方约束复原是指除了要求了解关于退化系统的传输函数H之外，还需要知道某种噪声的统计特性和噪声与图像的某些相关情况。根据所了解的噪声先验知识的不同，应该采用不同的约束条件，可得到不同的图像复原技术。

在最小二乘约束复原中，复原问题表现为在满足约束条件下，对于这类问题的有约束最小化问题，通常采用拉格朗日乘法进行处理。即寻找一个

，使得如下准则函数（4-10）最小。

(4-10)

式中，

为的线性算子，

为一个常数（称为拉格朗日乘子）。对式（4-10）求导，可得式（4-11）和（4-12）

(4-11)

(4-12)

令

，得式（4-13）

(4-13)

常数必须反复迭代齐整，直到满足约束条件。求解（4-12）的关键就是如何选用一个合适的变换矩阵

。

相对于无约束问题，有约束条件的图像复原更符合图像退化的实际情况，因此其适应面更加广泛。对式（4-13），若选择不同形式的矩阵

，则可得到不同类型的有约束最小二乘方类图像的复原方法。

4.3 维纳滤波复原

一般地，噪声源往往具有平坦的功率谱，即使不是如此，其随着频率的升高而下降的趋势也要比典型图像的功率谱慢得多。因此，可以认为功率谱的低频部分以信号为主，而高频部分则主要被噪声所占据。由于逆滤波滤波器的幅值通常随着频率的升高而升高，会增强高频部分的噪声。因此，为了克服以上缺点，提出了采用最小均方误差的方法(维纳滤波)进行模糊图像恢复。

维纳(wiener)滤波可以归于反卷积(或反转滤波)算法一类，它是由Wiener首先提出的，并应用于一维信号，取得了很好的效果。此后算法又被引入二维信号处理中，也取得了相当满意的效果，尤其是在图像复原领域，由于维纳滤波器的复原效果良好，计算量较低，并且抗噪性能优良，因而在图像复原领域得到了广泛的应用，并不断得到改进发展。

如果取式（4-14）为：

(4-14)

和

分别是图像和噪声的自相关矩阵，即得式（4-15）和式（4-16）：

(4-15)

(4-16)

并且它们都是正定对称矩阵，则根据约束条件有式（4-17）：

(4-17)

实际上就意味着使噪声与信号的比对复原图像的影响最小。图像和噪声的相关矩阵都是把图像当作随机过程来研究，从而描述其统计特性的量，在这里最小二乘方的最佳已经演变成均方误差最小准则下的最佳。

根据式(4-17)可求得频域维纳滤波公式（4-18）如下：

(4-18)

其中，

表示退化函数，

表示噪声功率谱，

表示为退化函数的噪声功率谱。

=1时，为标准维纳滤波器；

时，为含参数维纳滤波器。当没有噪声时，维纳滤波器则退化成理想反滤波器。

因此可以用一个比值K代替两者之比，从而得到简化的维纳滤波公式（4-19）如下：

(4-19)

维纳滤波是一种综合考虑了退化函数和噪声统计特征两个方面进行恢复处理的方法，它建立在认为图像和噪声是随机过程的基础上，而目标是找一个未污染图像的估计值，使它们的均方误差最小。维纳滤波复原法不存在极点，即当很小或变为零时，分母至少为K，而且

的零点也转换成了维纳滤波器的零点，抑制了噪声，所以它在一定程度上克服了逆滤波复原方法的缺点。

虽然维纳滤波避免了频域处理的病态问题，对噪声放大有自动抑制作用，且噪声越强，作用越明显，但是对具体问题，有时得到的结果不能令人满意。这是因为

（1）维纳滤波假设是线性系统，但实际上图像的记录和评价图像的人的视觉系统往往都是非线性的。

（2）维纳滤波是根据最小均方误差准则设计的最佳滤波器。这个准则不见得与人类视觉判断准则相符合。均方误差(MSE)准则对所有的误差，不管其处在图像中的位置如何，都赋以同样的权值。由于使均方误差最小化，图像进行了平滑。而人眼对暗处和高梯度区域的误差具有较大维纳滤波以一种并非最适合人眼的方式对图进行评价。

（3）维纳滤波基于平稳随机过程模型。但实际存在的图像并不一定都符合这个模型，大多数图像都是高度非平稳的，而且它只利用了图像的协方差信息，可能还有大量的有用信息没有充分利用。

4.4 Lucy-Richardson复原

Lucy-Richardson算法能够按照泊松噪声统计标准求出给定的PSF卷积后，最有可能成为输入模糊图像的图像。当PSF已知，但图像噪声信息未知时，也可以使用这个函数进行有效的工作。

从成像方程和poissian统计可以有（4-20）推导：

(4-20)

式中，是原始图像；

是PSF()函数；

是无噪声模糊图像。在已知

时，在每个像素点估计

的联合似然函数为式（4-21）：

(4-21)

当式(4-21)存在时，最大联合似然函数的解存在。解为式（4-22）：

(4-22)

则可以得到Lucy-Richardson迭代式，得式（4-23）：

(4-23)

可以看出每次迭代时，都可以提高解的似然性，随着迭代次数的增加，最终会收敛在具有最大似然性的解处。

MATLAB提供的deconvlucy( )函数，就是利用加速收敛的Lucy-Richardson算法对图像进行复原。deconvlucy( )函数还能够用于实现复杂图像重建的多钟算法中。这些重建算法都是基于原始Lucy-Richardson最大化可能性算法。

deconvlucy( )函数的调用方式如下：

J=deconvlucy( I，PSF，NUMIT, DAMPAR, WEIGHT, READOUT, SUBSMPL)

其中，I表示输入图像。PSF表示点扩散函数。其他参数都是可选参数：NUMIT表示算法的重复次数，默认值为10；DAMPAR表示偏差阈值，默认值为0（无偏差）；WEIGHT表示像素加权值，默认值为原始图像的数值；READOUT表示噪声矩阵，默认值为0；SUBSMPL表示子采样时间，默认值为1。

噪声痕迹是最大化可能性数据逼近算法的常见问题。经过多次重复处理，尤其是在低信噪比条件下，重建图像可能会出现一些斑点，这些斑点并不代表图像的真实结构，是输出图像过于逼近噪声所产生结果。要有效控制这些痕迹，可以使用deconvlucy( )函数的收敛参数DAMPAR，该参数指定了收敛过程中结果图像与原始图像背离程度的阈值。对于那些超过阈值的数据，将不再允许进行反复计算。

如果对采样不足的数据进行重建，而重建过程建立在一个较好的网格操作基础上，则重建效果就可以大大提高。如已知PSF具有较高的分辨率，则deconvlucy( )函数使用SUBSMPL参数指定采样不足的比例。另外，PSF还可以通过观察像素偏移或者光学模型技术获得。

五、几种方法下的图像复原效果图

本文利用MATLAB软件对运动模糊和复原进行仿真，研究使用MATLAB把退化现象模糊化及将复原的解决方案公式化。实验证明该软件功能强大，语言简洁易学，人机界面友好，具有丰富的技术支持，应用简单而效果良好。

5.1 逆滤波复原

图5-1逆滤波复原过程图，图(a)是选取的原始图像，图(b)是利用MATLAB对原始图像进行运动模糊和加噪声仿真而生成的仿真图像，模糊长度为10个像素。经过逆滤波复原图像为图（c）。

（a）原图像（b）模糊加噪图像

（c）复原图像

图5-1 逆滤波复原过程

5.2 有约束最小二乘方复原的实现

通过MATLAB仿真来实现有约束的最小二乘方复原，图5-3是有约束的最小二乘方复原图。分别取参数为0、1、0.001、0.01、0.05、0.1对应图5-3里面的(a)-(f)。

(a)

=0 (b)

(c)

= 0.3 (d)

=0.05

(e)

= 0.01 (f)

=0.001

图5-3 有约束最小二乘方在不同参数下的恢复情况图

5.3 维纳滤波复原的实现

图(a)是选取的原始图像，图(b)是利用MATLAB对原始图像进行运动模糊和加噪声仿真而生成的仿真图像，模糊长度为5个像素。

采用维纳滤波恢复算法对模糊图像进行恢复，在加噪声的情况下，参数k分别选取0．0001、0．001、0．005、0．01、0.1和1。各种图中(c)-(h)为对应参数下的恢复图像。图5-4有噪声下维纳滤波在不同参数下的恢复情况。

（a）原始图像（b）含噪声运动模糊图像d=5，v=0.001

(e) K=0.005 (f) K=0.01

(g) K=0.1 (h) K=1

图5-4 有噪声下维纳滤波在不同参数下的恢复情况

5.4 Richardson-Lucy复原的实现

图(a)是原始图像，图(b)是对原图进行运动模糊仿真而生成的仿真图像，模糊长度为10个像素，模糊方向为水平方向。

采用Richardson-Lucy恢复算法对模糊图像进行恢复，迭代次数参数分别选取20次、50次、100次、150次、200次和300次。所有图的(c-h)为对应迭代次数下的复原图像。

(a) 原图像 (b) 水平运动10像素加噪声图像

(e) 迭代100次 (f) 迭代150次

(g) 迭代200次 (h) 迭代300次

图5-5 R-L算法在不同参数下的复原图像

六、源码

(1)模糊图像的傅里叶频谱获取

j=imread('车牌.jpg');

figure(1);

imshow(j);

title('原图像');

len=20;

theta=30;

psf=fspecial('motion',len,theta);

j1=imfilter(j,psf,'circular','conv');

figure,imshow(j1);

title('PSF模糊图像');

J=rgb2gray(j);

K=fft2(J); %傅里叶变换

M=fftshift(K); %直流分量移到频谱中心

N=abs(M); %计算频谱幅值

P=(N-min(min(N)))/(max(max(N))-min(min(N)))*225;%归一化

Figure；

imshow(P);

title('原图像的傅里叶变换频谱');

J1=rgb2gray(j1);

K1=fft2(J1); %傅里叶变换

M1=fftshift(K1); %直流分量移到频谱中心

N1=abs(M1);%计算频谱幅值

P1=(N1-min(min(N1)))/(max(max(N1))-min(min(N1)))*225;%归一化

Figure；

imshow(P1);

title('模糊图像的傅里叶变换频谱');

(2)模糊长度获取程序

f1=rgb2gray(j1);

f1=im2double(f1);

h = fspecial('Sobel'); %Sobel算子

J = conv2(f1,h,'same'); %Sobel算子微分

IP=abs(fft2(J)); %图像能量谱密度

S=fftshift(real(ifft2(IP)));

figure,plot(S);

title('自相关图')

(3)逆滤波复原实现程序

%对运动模糊图像进行逆滤波复原

clear all

I=imread('lena.bmp');

LEN=10;

THETA=5;

PSF=fspecial('motion',LEN,THETA);

g=imnoise(MF,'gaussian',0,0.0001);

MFN=imadd(MF,im2uint8(noise));

NSR=sum(noise(:).^2)/sum(MFN(:).^2);

wnr1=deconvwnr(g,PSF);

figure(1),subplot;imshow(I);

figure(2),subplot;imshow(MF);

figure(3),subplot;imshow(wnr1);

(4)有约束最小二乘方复原的实现程序

clear all;

clc;

%通过模拟水平运动模糊建立退化函数

d=2;

h=zeros(2*d+1,2*d+1);

h(d+1,1:2*d)=1/(2*d);

%模糊原图像并加入噪声

fig1=imread('lena.bmp');

[m n]=size(fig1);

fe=zeros(m+2*d,n+2*d);

fe(1:m,1:n)=fig1;

he=zeros(m+2*d,n+2*d);

he(1:2*d+1,1:2*d+1)=h;

F=fft2(fe);

H=fft2(he);

g=imnoise(uint8(ifft2(F.*H)),'gaussian',0,0.0001);

G=fft2(double(g));

%最小二乘平方滤波器

p=[0 1 0;1 -4 1; 0 1 0];

pp=zeros(m+2*d,n+2*d);

pp(1:3,1:3)=p;

p=fft2(pp);

r=0.001;

F_est=(conj(H)./(abs(H).^2+r.*abs(p).^2)).*G;

fig_est=real(ifft2(F_est));

%显示结果

figure(1),imshow(fig1);

figure(2),imshow(uint8(g(d+1:m+d,d+1:n+d)),[min(g(:)) max(g(:))]);

figure(3),imshow(uint8(fig_est(1:m,1:n)),[min(fig_est(:)) max(fig_est(:))]);

(5)维纳滤波复原实现程序

clear all;

clc;

%通过模拟水平运动模糊建立退化函数

d=2;

h=zeros(2*d,2*d+1);

h(d+1,1:2*d+1)=1/(2*d);

%模糊原图像并加入噪声

fig1=imread('lena.bmp');

[m n]=size(fig1);

fe=zeros(m+2*d,n+2*d);

fe(1:m,1:n)=fig1;

he=zeros(m+2*d,n+2*d);

he(1:2*d+1,1:2*d+1)=h;

F=fft2(fe);

H=fft2(he);

g=imnoise(uint8(ifft2(F.*H)),'gaussian',0,0.0001);

G=fft2(double(g));

%维纳滤波器设计

k=0.1;

F_est=((abs(H).^2)./(abs(H).^2+k)).*G./H;

fig_est=real(ifft2(F_est));

figure(1),imshow(fig1);

figure(2),imshow(g);

figure(3),imshow(fig_est);

(6)Lucy-Richardson复原的实现程序

%Lucy_Richardson滤波复原实现

clear all;

clc;

g=imread('lena.bmp');

figure(1),imshow(g);title('原图像');

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

Blurred=imfilter(PSF,'symmetric','conv');

V=.0001;

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

figure(2),imshow(BlurredNoisy);title('模糊和加噪图');

luc1=deconvlucy(BlurredNoisy,PSF,50);

figure(3),imshow(luc1);title('复原图像，NUMIT=50');

（6）客观评价三种方法实现程序（平均平方误差、峰值信噪比、归一化互相关）

f0=fig1;

g1=x;

fima=y;

[x_size,y_size]=size(g1) % 模糊噪声图像参数

fima=double(fima);

fprintf(1,'MSE__blurnoisy=%4.5e,PSNR__blurnoisy=%4.5e,

NCR__blurnoisy=%4.5e\n', MSE_blur_noisy,PSNR_blur_noisy,NCR_blur_noisy);

%---------------------去模糊后的图像参数----------------------------------%

MSE_deblurring=sum(sum((fima-f0).^2))/(x_size*y_size); % MSE

PSNR_deblurring=10*log10(255*255/MSE_deblurring); % PSNR

NCR_deblurring=sum(sum(fima.*f0))/sum(sum((f0).^2)); %NCR

fprintf(1,'MSE_deblurring=%4.5e,PSNR_deblurring=%4.5e,

NCR_deblurring=%4.5e\n', MSE_deblurring, PSNR_deblurring, NCR_deblurring);