2.5 使用Lucy-Richardson算法的迭代非线性复原
摘要
数字图像处理是将图像信号转换成数字格式，并通过计算机对它们进行处理。图像复原过程往往是对提高图像质量起着重要的作用的数字图像处理方法。图像处理中的一个重要的研究分支是图像重建，其意义在于要检测到获得物体的内部结构图像，而不会其造成任何物体上的损伤。在本文中，先对图像复原与图像重建进行概述，然后介绍几种图像复原技术与图像重建方法。通过MATLAB实验程序获得实际处理效果。
关键词：图像复原；图像重建；MATLAB
Abstract
Digital image processing is to convert the image signal into a digital format and process them through the computer. Image restoration process is often to improve the image quality, it plays an important role in digital image processing methods. Image reconstruction is an important research branch of image processing, in the sense that the object to be detected to obtain images of internal structures without causing objects any damage. In this article, firstly, it will introduce image restoration and reconstruction principle, and then introduce several image restoration techniques and image reconstruction methods. The finally treatment effect obtained by MATLAB experimental procedures.
Key words: image restoration; image reconstruction; MATLAB
基于MATLAB的图像复原与重建设计
前言
随着网络和通信技术的发展，数字图像处理与分析技术已经在科学研究、工业生产、军事技术、医疗卫生、教育等许多领域得到了广泛应用，并产生了巨大的经济效益和社会效益，对推动社会的发展和提高人们生活水平都起到了重要作用[1]。图像复原与重建是数字图像处理的一个重要组成部分，并已被广泛的应用。MATLAB图像处理工具为数字图像处理提供了一个稳定、广泛的软件实现平台。
1 MATLAB的简介
MATLAB是MathWorks公司开发的一款工程数学计算软件。它是集数值符号计算，高质量图形可视化与界面设计为一体。由于其功能强大、操作简单，已成为国际上科学界最具影响力、最有活力的软件。
矩阵是MATLAB的基本数据单位，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解决问题事件比用CFORTRAN等语言简捷方便得多。
MATLAB包括拥有数百个内部函数的主包和三十几种工具包（Toolbox）。工具包又可以分为功能性工具包和科学工具包[12]。功能性工具包用来扩充MATLAB的符号计算，可视化建模仿真，文字处理及实时控制等功能；科学工具包是专业性比较强的工具包，它包括控制工具包、通信工具包、信号处理工具包 [9]。
MATLAB的开放性广受用户欢迎，除去内部函数，MATLAB的所有主包文件和各种工具包都是可读并可以修改的文件，通过对源程序的修改或添加，用户可以构造自己的专用工具包。
MATLAB是一种用高级技术计算语言和交互式环境，它集算法开发、数值计算、数据分析以及数据可视化为一体。有了它，比用传统的编程语言，如C、C++等，更快的解决技术计算问题。
MATLAB高级语言可以用于技术计算；它所形成的开发环境可管理代码、文件和数据；数学函数可用于线性代数、概率统计、傅里叶分析变换、优化、筛选以及积分等；二维和三维图形函数可用于可视化数据；各种工具可用于构建自定义的图形用户界面；各种函数可将基于MATLAB的算法与外部应用程序和语言；它具有非常广泛的应用，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。
一系列支持图像处理操作的函数可以组成图像处理工具。它所支持的图像处理操作有：图像的邻域操作、图像的区域操作、图像的几何操作、图像变换、图像恢复、图像增强，图像复原、图像重建、线性滤波、图像分析和统计等。下面就MATLAB在图像处理中各方面的应用分别进行介绍[6]。
（1）读写和显示数字图像的文件格式。imread()为图像文件读入函数，可以用来读取如：bmp、tif、gif、jpg、png、xwd等格式图像；imwrite()为图像写出函数，而imshow()、image()为图像显示函数。
（2）图像处理的基本运算。加、减等线性运算，卷积、相关、等非线性运算都是MATLAB提供图像处理的基本运算。例如，conv2(X，Y)实现了X，Y两幅图像的卷积。
（3）图像变换。MATLAB提供了离散傅里叶变换（DFT）、快速傅里叶变换（FFT）、离散余弦变换（DCT），连续小波变换（CWT）、离散小波变换（DWT）及其反变换等变换。
（4）图像的分析与图像增强。校正、直方图均衡、中值滤波等都是MATLAB提供的关于图像统计的计算。
（5）图像复原与重建。可用逆滤波复原、维纳滤波复原等方法实现图像复原。傅里叶反投影重建、卷积法重建、代数法重建是常用的图像重建技术。
2图像复原
2.1图像复原的基本概念
图像在形成、传输和记录过程中，由于受到多种原因的影响，图像的质量就会有所下降，典型的表现为图像模糊、失真、有噪声等，这一过程称为图像的退化[5]。
图像复原是试图利用退化过程的先验知识使已退化的图像恢复本来面目，即根据退化的原因，分析引起退化的环境因素，建立相应的数学模型，并沿着使图像降质的逆过程恢复图像[5]。目的在于消除或减轻在图像获取以及传输的过程中造成的图像品质下降，恢复图像的本来面目。因此，复原技术就是把退化模型化，并采用相反的过程进行处理，以便尽可能复原被退化图像的本来面目。
广义上讲，图像复原是一个求逆问题，逆问题经常存在非唯一解，甚至无解。要想恢复全真的景物图像比较困难。为了得到逆问题的有用解，图像复原本身往往需要一个质量标准，即衡量接近全真景物图像的程度，或者说，对图像的估计是否达到最佳的程度。需要有先验知识以及对解的附加约束条件。
典型的图像复原是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用各种逆退化处理的方法进行恢复，使图像质量得到改善。
2.2图像退化的数学模型
一般来说，图像的生成可以简单地被描述为如下数学模型：
f(x,y)是成像景物，H是综合退化因子，g(x,y)是退化图像。
图像f(x,y)可以表示为：
用卷积符号
表示为：
因此还有：
式中，
是像素点的特性函数，
为冲击响应。
假定成像系统是线性移不变系统：
退化模型如图所示
不考虑加性噪声：
考虑加性噪声：
卷积等同于频域内乘积：
2.3逆滤波复原
逆滤波复原法也叫做反向滤波法，其主要过程是首先将要处理的数字图像从空间域转换到傅里叶频域中，进行反向滤波后再由频率域转回到空间域，从而得到复原的图像信号[5]。
1．在不考虑噪声的情况下：
上式两边进行傅里叶变换得
则原始图像
=
然后进行傅里叶逆变换，就可以得到原始图像。由此可看出，如果已知退化图像的傅里叶变换和“滤波”传递函数，则可以求得原始图像的傅里叶变换，经反傅里叶变换就可以求得原始图像f(x,y)，这就是逆滤波法的基本原理。但在实际中用逆滤波法存在病态的情况：当H(u,v)=0时，或非常小的数值点上，F(u,v)将变成无穷大或非常大的数。
2．在有噪声的情况下：
逆滤波原理可以写成：G(u,v)=F(u,v)H(u,v)+N(u,v)
写成逆滤波的方式：
(u,v)=F(u,v)+
但实际用逆滤波存在病态的情况：
噪声存在，当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为特殊常数。 ②
①的维纳滤波要求未退化图像和噪声的功率必须是已知的。虽然用②近似的方法能得到好的结果，但功率谱比常数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
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
end
F2=H1.*F0;
I3=ifft2(fftshift(F2));
subplot(2,2,4);
imshow(uint8(I3));
title(维纳滤波复原图);
运行结果如下：
原始图：
复原后图像：
经过仿真，如上图所示，可以看出逆滤波复原与维纳滤波复原的区别和联系。维纳滤波后虽然仍有一些噪声存在，但已经和原图很接近了。因为原图像和噪声函数都是已知的，可以正确的估算参量。
2.5使用Lucy-Richardson算法的迭代非线性复原
L-R算法是一种迭代非线性复原算法，它是从最大似然公式印出来的，图像用泊松分布加以模型化的。当下面这个迭代收敛时模型的最大似然函数就可以得到一个令人满意的方程：
*代表卷积，
代表未退化图像的估计，g和h和以前定义一样。这个算法的本质是显而易见的。它的非线性本质是在方程右边用
来除产生的[4]。
在IPT中，L-R算法是由名为deconvlucy的函数完成的，此函数的语法为
fr= deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT)
其中，fr代表复原的图像，g代表退化的图像，PSF是点扩散函数，NUMIT为迭代次数（默认为10次），DAMPAR是一个标量，它指定了结果图像与原图像g之间的偏离阈值。WEIGHT是一个与g同样大小的数组，它为每一个像素分配一个权重来反映其重量。
L-R算法程序：
I=imread('C:\ok\苹果.jpg');
PSF=fspecial('gaussian',5,5) ;
Blurred=imfilter(I,PSF,'symmetric','conv');
V=.003;
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方法的一种解释就是将图像数据看成随机量，它们与另外一族可能的随机量之间有着某种似然性。似然函数用
、
和
来加以表达，然后，问题就变成了寻求最大似然函数。在盲去卷积中，最优化问题规定的约束条件并假定收敛时通过迭代来求解，得到的最大
和
就是还原的图像和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');
[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,[]);
模拟实验结果如下：
该算法优点是,同时恢复了图像和点扩张函数,在对失真情况毫无先验知识的情况下, 仍能实现对模糊图像的恢复操作。利用 MATLAB实现的图像恢复, 并对恢复图像的失真情况做了改善。在进行图像恢复时,重建 PSF,对图像进行重建, 得到恢复的图像。
3图像重建
3.1 图像重建的概述
图像处理中的一个重要研究分支是物体图像的重建，它被广泛的应用于检测和观察中，而这种重建方法一般是根据物体的一些横截面部分的投影而进行的。在一些应用中，某个物体的内部结构图像的检测只能通过这种重建才不会有任何物理上的损伤。由于这种无损检测技术的显著优点，因此，它的适用面非常广泛，它在各个不同的领域都显示出独特的重要性。例如：医疗放射学、核医学、电子显微、无线和雷达天文学、光显微和全息成像学及理论视觉等等领域都有应用。在医学影像处理中是医学图像获取的重要方法。如医疗放射学、核医学，电子显微等领域是必不可少的技术，在工业生产中的无损检测技术图像重建也扮演重要的角色。
图像重建经过多年研究已取得巨大进展，产生了许多有效的算法，如：傅里叶反投影法、卷积反投影法、代数法、迭代法等，其中以卷积反投影法运用最为广泛。近年来，由于与计算机图形学相结合，把多个二维图像合成三维图像，并加以光照模型和各种渲染技术，已能生成各种具有强烈真实感的高质量三维人工合成图像。
3.2傅里叶反投影重建
傅里叶反投影重建方法可以说是最简单的一种变换重建方法，一个三维（或二维）物体，它的二维（或一维）投影的傅里叶变换恰好与此物体的傅里叶变换的主题部分相等，傅里叶变换的重建方法正是以此为基础的。该重建方法最早于1974年由Shepp和Logan提出，该方法是建立在“投影切片定理”这一理论基础之上的。根据此结论，可以先对投影进行旋转和傅里叶变换以构造整个傅里叶变换域中的各个方向的切片数据，然后再对其进行傅里叶反变换即可得到重建后的目标（原空间域中的图像）。
傅里叶变换重建的原理如下：
令f(x,y)代表一图像函数，则此二维函数的傅里叶变换为：
而图像在x轴上的投影为：
投影的一维傅氏变换为：
，它恰与二维傅氏变换的表达式一致。
即：
现在假设将函数投影到一条经过旋转的直线上，该直线的旋转角度为
。
定义旋转坐标为：
而将函数投影的直线选为x轴。投影点通过对距离t轴为
处的一平行线进行函数积分，因此，该投影可如下表示：
这里，积分路径是沿着
直线进行的。此投影的一维傅氏变换为：
，
展开后为：
为使展开式与投影的二维傅里叶变换相等，把指数项做某种代换得到：
因而，若点(u,v)在一条
角一定而距原点距离为r的直线上，投影变换将于二维变换中的一直线有相同的傅氏变换，即：
若投影变换
中的所有
及
值都是已知的，则图像的二维变换也是可以确定的。为得到图像函数，我们必须进行反变换运算，即：
这些结论很容易推广到三维情形中。令：
表示一物体，这里f可以为实数或负数。它的三维傅氏变换由下式给出
而变换的核心部分是：
通过定义，纵剖面或在
面上的投影是：
注意到
的二维傅里叶变换正好等于上述三维变换的核心部分。这也说明了如果投影在
平面上旋转了
角度，相应的傅里叶变换部分也将在变换域内的
平面内转过
角。这样，投影可以采用不同的方向角
插入到三维变换域中。建立一个傅里叶变换空间需要很多的投影。最后通过傅里叶反变换重建图像
。既然在三维空间中的任意平面都可以被重建，那么，一个二维图像
的重建也不失一般性。
我们可以重写二维投影方程，定出
即投影平面
：
，这里ds是光线几何路径中的微分长度。
傅里叶变换的结论由下面给出：
若已知无数的投影，从极坐标
中计算得到的投影变换推出在矩形平面
中的傅里叶变换并不困难。但是，若只有有限个投影是有效的，则可能需要在变换中插入一些数据。另外需要注意的是，虽然只需一维傅里叶变换的投影数据就可构成变换空间，但图像重建则需要二维反变换。由此，我们得出一个结论，即：三维图像不能在得到部分投影数据的过程中局部的重建，而必须延迟到所有投影数据都获得之后才能重建。
3.3 卷积法重建
卷积法重建的基本思路和方法：
逆投影原理：从各个方向得到的投影逆向返回到该方向的各个位置，如果对多个投影方向中的每个方向都进行这样的逆投影，就可能建立平面上的一个部分。典型的方法是卷积逆投影重建。
卷积重建法是一种变换重建法，可以根据傅里叶变换投影定理推出。
按照二维傅里叶反变换标准定义，有
作代换：
，
写成极坐标（R,
）的形式：
①
利用傅里叶变换共轭对称性，有：
②
令：
③
则③可以表示为：
， ④
在②中，当用FFT计算投影数据的傅里叶变换F(R,
)时，投影数据g(
)总被有限截断。当
的采样间隔为d时，在变换域R的变化范围为
到
，于是投影反变换重建公式可以近似写成：
采用标记：
⑤
根据前式③，结合傅里叶投影定理可知：
由上式可以得出，要实现对已经得到的投影数据实现图像重建，则可以采取两步：首先将投影数据
和相应脉冲滤波器⑤进行卷积，然后由式④对不同旋转角
求和，就能实现图像重建。
卷积可以看成一种滤波手段，卷积投影相当于先对数据滤波再将结果逆投影回来，这样可以使模糊得到校正。
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
1004
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
