CN103201765A - 用于从所观察到的数字图像序列恢复数字图像的方法与设备 - Google Patents

Abstract

本公开涉及用于从所观察到的数字图像序列恢复数字图像的方法与设备。根据本发明，用于从所观察到的数字图像序列（y₁,…,y_T）恢复数字图像（x）的计算机实现方法包括步骤：获得所观察的数字图像（y_t）；基于所观察到的图像（y_t）估计点扩散函数（f_t）；基于所估计的点扩散函数（f_t）和所观察到的图像（y_t）来估计恢复的数字图像（x）并且重复以上步骤。为了校正透镜的光学像差，可以使用透镜的点扩散函数。

Description

用于从所观察到的数字图像序列恢复数字图像的方法与设备

技术领域

本发明涉及用于从一个或多个观察到的数字图像重构真实图像或图像序列的方法与设备。所观察到的图像可能由于透镜的光学像差、不期望的运动（相机抖动）或类似的现象而模糊。

背景技术

在如理论上由近轴光学器件所描述的理想光学***中，点源发射出的所有光线都汇聚到焦平面中的单个点，形成清楚而清晰（sharp）的图像。光学***与这种行为的偏离被称为像差，会造成图像不期望的模糊。

照相透镜的制造商试图通过组合若干个透镜来最小化光学像差。复合透镜的设计与复杂性依赖于各种因素，例如，光圈大小、焦距和对失真的限定。光学像差是不可避免的而且透镜的设计总是在包括成本的各个参数之间进行折中。

类似地，当从地球表面观察天文对象时，其发射或反射的光必须穿过大气，这会产生模糊的观察图像。因此，在天文学中，给所观察到的天体的图像去模糊是一个基本的问题。由于模糊不仅未知而且由于大气湍流引起的折射指数的波动还随时间和空间持续变化的事实，这个问题很复杂。众所周知，按湍流平稳的时间量程的曝光时间（即，比十分之一秒还要短）产生包含所研究天体的高频信息的图像。这个事实在散斑成像（Speckle Imaging）中被利用，其中散斑成像是用于恢复封装在短时间曝光图像中的高频信息的各种技术与算法的集合。由于大气湍流的随机本质，散斑成像技术必须考虑大量的图像来实际恢复衍射受限的信息。

在散斑干涉成像中，短曝光图片在傅立叶域中被处理。通过CCD技术中提供优越信噪比的新发展（电子倍增CCD照相机），近年来所谓的幸运成像（Lucky Imaging）方法变得很流行。通常要收集大量的图像（>10000），其中只有非常小的一个百分比用于产生一个高质量的图像。对于所选图像的组合，已经建议了不同的注册与超分辨率方法。

另一方面，相机抖动是手持式、更长曝光照片的一个常见问题，尤其发生在低照明情况下，例如建筑物内部。除了例如平移摄影的少数例外，相机抖动都是不期望的，因为它常常破坏细节并模糊图像。

透镜上光学像差、大气模糊或者相机抖动及许多其它因素的影响可以表达为（或者有时候至少近似为）实时图像的真实清晰的变换，这种变换可以写成矩阵矢量乘法（MVM）

y=Ax                 (1)

其中x是代表兼容或真实图像的矢量，y是观察到的图像，而A是表示线性图像变换的矩阵。

但是，在一般的情况下，在给定x和y的情况下估计A是计算代价高昂的。因此，期望能适当地限定解空间，即，不为了有效地执行MVM而过分牺牲可表达性，因为这些对于盲解卷积算法是基本的操作。

因此，本发明的一个目标是提供用于处理光学像差、相机抖动或大气模糊的影响的方法，所述方法可以有效地实现并且有可能在线执行。

发明内容

这个目标是根据本发明由根据独立权利要求的方法与设备实现的。有利的实施方式在从属权利要求中定义。

本发明的关键思想是采用在维持计算可行的同时对许多现实世界应用能够充分表达的一类线性过滤器。

根据本发明，满足这种需求的线性过滤器是所谓的点扩散函数（PSF）。本发明范围内的点扩散函数可以是空间不变的或者空间变化的，即，依赖于要应用其的图像片的位置。

在空间不变PSF的情况下，这些过滤器可以由（具有某个长度k的）矢量a表示，而且它们的操作可以定义为x与a的卷积，

$y_i=\underset{j=0}{\overset{k-1}{\mathrm{\&Sigma;}}}{a}_j{x}_{i-j},$ 对于0≤i<m         (2)

其中，x的长度是n而y的长度是m。有效部分可以从全卷积中裁剪，使得m=n-k+1。由于变换（2）在x中是线性的而且可以写成对于0≤i<m和0≤j<k，对于A_i,i+j=aj，有y=Ax。换句话说，A在每一行中都包含a的偏移副本。对于这种结构的A，MVM可以利用具有适当零填充的FET在O(n log n)次乘法中执行。如果信号比PSF长得多，即，n>>k，则通过把信号切片并且使用Stockham的重叠相加（OLA）方法（1966年ACM春季联合计算机会议，1966年4月26-28日会议录第229-233页，T.STOCKHAM J_R.High-speedconvolution and correlation），MVM可以甚至更快地被处理。如果q是用于片的FFT的尺寸，则OLA花费O(n log q)。

在空间变化过滤的情况下，图像可以分成不同的片，其中每一片都可以用其自己的过滤器来处理。

如果片重叠，那么这种空间变化的过滤器可以有效地通过根据本发明的重叠相加（OLA）方法的修改来实现。通常的OLA方法把写图像裁剪成重叠的片、利用某种开窗函数抑制（dampen）每一片的边界、利用相同的过滤器卷积每一片，并且随后把变换后的片相加，以获得输出图像。

通过用其自己的过滤器处理每一片，可以获得空间变化的线性过滤器。更特别地，对于一个输入图像x的多个p个重叠的片中的每一片r（0≤r<p），可以使用只对x中对应片的间隔非零的开窗函数w(r)（表示为与x相同长度的矢量）。w(r)与x关于输入的乘法把x中在第r个片之外的所有输入都设置成零，给剩余的输入开窗。用a(r)（0=r<p）指示长度为k的p个过滤器，根据本发明的空间变化过滤器可以定义为

$y_i=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{k-1}{\mathrm{\&Sigma;}}}{a}_j^{\left(r\right)}{w}_{i-j}^{\left(r\right)}{x}_{i-j},0\&le;i<m.---\left(3\right)$

对于固定的片尺寸和固定个数的片，可以计算重叠的量和片的位置。空间变化过滤器(0)(k-1)的参数是k个依赖位置的点扩散函数a⁽⁰⁾,…,a^(k-1)。

为了避免输出图像y中片的重叠区域的假象，所选的窗口函数必须加起来为一，即

$\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{w}_i^r=1\mathrm{for}0\&le;i<m.---\left(4\right)$

在实践当中，规格化窗口函数可能总是加强这个属性。

由于x在等式（3）中看起来只是线性的，因此它可以写成y=Ax。过滤器矩阵A是由下式给出的

$A=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{A}^{\left(r\right)}\mathrm{Diag}\left(w^{\left(r\right)}\right)---\left(5\right)$

其中A^(r)是对应于卷积a^(r)（0=r<p）的矩阵，而Diag(v)是具有沿其对角线的矩阵v的对角矩阵。但是，这种表示没有描述对于A可以如何有效地计算MVM，也没有描述对于A^T可以如何有效地计算MVM。例如，A可以等效地表示为以下矩阵乘积之和：

$A=Z_y^T\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{C}_r^T{F}^H\mathrm{Diag}\left(F{Z}_{\mathrm{\&alpha;}}{a}^{\left(r\right)}\right)F{C}_r\mathrm{Diag}\left(w^{\left(r\right)}\right)---\left(6\right)$

其中（i）C_r是从长度为n的矢量裁剪掉第r个片的矩阵；（ii）Z_a是把零附加到a(r)的零填充矩阵，使得其大小与片的大小匹配；（iii）F是离散傅立叶变换（DFT）矩阵（由FFT实现）；（iv）F^H是DFT矩阵的厄密共轭（由逆FFT实现）；及（v）Z_v是向矢量预加零的零填充矩阵，使得其大小与从求和得到的矢量的大小匹配。等式（6）忠实地代表A的证据是直接从卷积的FFT实现得到的。

从右向左读等式（6），这个表达式简洁地描述了有效地计算Ax所需的步骤。也可以把A^T读作

$A^T=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}\mathrm{Diag}\left(w^{\left(r\right)}\right)C_r^T{F}^T\mathrm{Diag}\left(F{Z}_a{a}^{\left(r\right)}\right)\stackrel{\&OverBar;}{F}{C}_r{Z}_y$              (7)

$=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}\mathrm{Diag}\left(w^{\left(r\right)}\right)C_r^T{F}^H\stackrel{\&OverBar;}{\mathrm{Diag}\left(F{Z}_a{a}^{\left(r\right)}\right)}{\mathrm{FC}}_r{Z}_y$

其中是A的关于组成部分的共轭复数。等式（7）是正确的，因为A是实数值矩阵，因此

从右向左读等式（7）描述了为A^T有效计算MVM所需的步骤。换句话说，步骤与等式（6）中的相似，但是是在最后而不是在开始开窗，而且对于PSF的FFT的共轭复数，导致计算关于片的相关，而不是卷积。

对于利用空间变化过滤器A的非盲解卷积，对于A和A^T的有效MVM就足够了。

对于盲解卷积，由于等式（3）在k个PSF中也是线性的，因此可以定义一个矩阵X，使得y=Ax=Xa，其中a指示PSF的堆栈序列a⁽⁰⁾,…,a^(p-1)。现在等式（6）可以利用D(v)w=D(w)v和从矢量a裁剪掉第r个PSF的某个矩阵B_r重写，

$X=Z_y^T\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{C}_v^T{F}^H\mathrm{Diag}\left({\mathrm{FC}}_r\mathrm{Diag}\left(w^{\left(r\right)}\right)x\right){\mathrm{FZ}}_a{B}_r---\left(8\right)$

这个表达式允许通过简单地取用于X的表达式的转置来导出对于X^T的有效MVM的算法，

$X^T=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{B}_r^T{Z}_a^T{F}^T\mathrm{Diag}\left({\mathrm{FC}}_r\mathrm{Diag}\left(w^{\left(r\right)}\right)x\right)\stackrel{\&OverBar;}{F}{C}_r{Z}_y$                   (9)

$=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}{B}_r^T{Z}_a^T{F}^H\stackrel{\&OverBar;}{\mathrm{Diag}\left(F{C}_r\mathrm{Diag}\left(w^{\left(r\right)}\right)x\right)}{\mathrm{FC}}_r{Z}_y$

其中

可以再次使用，因为X是实数。换句话说，由等式（9）暗示的用于X^Tv的算法包括把v分成片、把它们与来自x的片关联，并最终把结果求和。

以上所述对A、A^T、X和X^T的MVM的计算复杂性与用于空间不变过滤的重叠相加方法的相同，大约是O(n log q)，其中q是用于片的FFT的尺寸。换句话说，根据本发明的空间变化线性过滤器可以与空间不变过滤器一样有效地实现。由于PSF a(r)比输入图像x小得多，因此存储空间变化过滤器的存储器需求比通用线性变换所需的O(mn)小得多。但是，为上述变换（例如，在盲解卷积中）估计参数通常更加昂贵，因为参数的个数随着片的个数而增加，而对于空间不变过滤器，只需要估计单个过滤器。

至于表达能力，由一个PSF a表示的空间不变过滤器是根据本发明的空间变化过滤器的一个特例。当所有PSF a(r)=a时，利用（4），可以看到（3）简化成（2）。

在另一个极端，根据本发明的空间变化过滤器可以利用m个图像片实现任何线性变换A，对于A的每一行有一个图像片。于是，所有的窗口函数都可以设置成常量1/m，而且特征化空间变化过滤的PSF设置成A的行。这种情况完全退化成片重叠而且与信号一样长的PSF只评估一次。但是，显示根据本发明的空间变化过滤实际上覆盖了从空间不变过滤到任意线性变换的整个范围，为了更具表达性而折中了一些计算效率。

解决该问题的本发明性方法具有以下主要优点：

-它具有低资源需求，因为图像可以流方式处理；这避免了存储图像的需求，并且导致更有效的算法（关于运行时间和存储器使用）。

-典型的在线方法通常需要学习速率的设置。本发明性方法是基于对图像估计的倍增更新，有效地绕过了设置学习速率的需求，这简化了算法及其实现。

-该方法实现起来非常简单，而且，如果期望的话，它可以在成像硬件本身当中实现，以便提供动态去模糊。

附图说明

当联系附图研究以下本发明不同实施方式的具体描述时，本发明的这些及其它方面、优点和特征将更容易理解，附图中：

图1示出了具有电机驱动平台的实验装备和用于测量透镜的点扩散函数的照相机。

图2示出了用于模拟的点扩散函数。

图3示出了一个小的PSF集合如何可以参数化平滑变化的模糊：（左边）具有实际相机抖动拍摄的网格，（中间）被由九个PSF（右边）参数化的EFF框架模糊的网格。

图4示出了如何同时捕捉被实际相机抖动及其空间变化的PSF模糊的图像；（a）真实的图像与点的网格组合成（b）RGB图像，该图像是带相机抖动拍摄的，而且（d）分成蓝和红通道，以便分离描述模糊的PSF和被模糊的图像。

图5示出了如何对单像素点只应用一次全部可能的单应（homography）和可能的相机抖动如何通过线性组合点网格的不同单应来生成。

图6示出了模糊估计阶段的概述。

图7示出了对于天文成像的在线盲解卷积的算法。

具体实施方式

为了简化，本发明将对矢量值的图像进行描述。对矩阵值图像的通用性是直接的。

I.在本发明的第一种实施方式中，将显示如何减轻由于不完美光学器件造成的图像降级。校准步骤可以用于在空间变化点扩散函数（PSF）中编码光学***的光学像差。然后，校正后的图像可以利用那个函数通过不平稳的解卷积获得。通过在图像形成模型中包括拜耳（Bayer）阵列，去马赛克可以作为解卷积的一部分来执行。

在第一方面，本发明适于校正透镜或光学***中的单色像差，包括球面像差（在球面透镜中，焦距是从轴开始的距离的函数）及多种离轴像差：当射线的交点关于其轴偏移时，在倾斜的光束中发生彗差（coma）；当焦面非平面时，发生场曲（field curvature）；当矢状面和切向的焦面不一致时（即，***对于离轴光束不是旋转对称的），散光（astigmatism）指示这种情况。所有这些单色像差都导致跨图像变化的模糊。任何这种模糊都可以在EFF框架中通过适当地选择局部模糊过滤器f⁽⁰⁾、…、f^(p-1)来表示。

在第二方面，本发明还可以用于校正透镜或光学***中的色像差。包括玻璃在内的大多数材料的折射率依赖于所透射的光的波长。在轴的方向，这导致透镜的焦点是该波长的函数（纵向色像差）；离轴，我们观察到由于不同波长的不同焦距直接暗示图像比例随波长稍有变化的事实所造成的横向色。例如，蓝光的模糊与红光的模糊在形状与位置上都稍有不同。根据本发明，这种色像差可以通过利用独立的空间变化PSF给三个颜色通道建模来描述。更具体而言，假定x_R、x_G和x_B是真实图像的三个颜色通道，为了简化，它们都是列矢量，而且对于每个颜色通道，空间变化模糊可以用来自等式（3）的EFF矩阵F_R、F_G、F_B来建模。于是，进入图像传感器的颜色过滤器阵列（CFA）的，即，在拜耳过滤/马赛克之前的，图像z可以写成

$z=\left[\begin{array}{c}{z}_R\\ z_G\\ z_B\end{array}\right]=\left[\begin{array}{c}{F}_R{x}_R\\ F_G{x}_G\\ F_B{x}_B\end{array}\right]=\left[\begin{array}{ccc}{F}_R& 0& 0\\ 0& F_G& 0\\ 0& 0& F_B\end{array}\right]x---\left(10\right)$

在第三方面，本发明还包括校正所谓渐晕的有效途径。因为倾斜的光束不完全到达焦平面，所以图像的强度朝图像的角落减小。这可以通过拍摄平场帧，即，均匀背景的图像，并且用它分割图像来校正。虽然这是直接了当的，但是本发明还可以通过忽略节能限定来处理渐晕。在那种情况下，等式（1）中的过滤器f^(r)不必加起来等于一，即，只需要对于所有的j和r都有

而且

通过允许调光过滤器，本发明的过程对渐晕进行校正。应当指出，等式（2）不受放松节能限定的影响。

到目前为止，已经描述了从包括三个颜色通道x_R、x_G和x_B的真实全色图像x到由于透镜像差而模糊的全色图像z的正演模型。就是该图像刚好在被马赛克之前将进入CFA，即，z包含所有三个颜色通道的完整色度信息。CFA的操作可以描述为由某个矩阵D表示的线性图，其结果y=Dz将是命中该CFA后面的光敏传感器的图像，它可以利用存储在照相机中的原始图像来识别。矩阵D是列比行多三倍的矩形矩阵。

正演模型把透镜像差和拜耳过滤组合到单个矩阵A中并且添加了噪声，即

$y=D\left[\begin{array}{ccc}{F}_R& 0& 0\\ 0& F_G& 0\\ 0& 0& F_B\end{array}\right]x+n=\mathrm{Ax}+n---\left(11\right)$

假定拜耳矩阵D中的权重是固定的而且是已知的（简单的拜耳矩阵可以无视颜色通道之间的串扰而使用），线性变换A，即，PSF，是由确定用于三个颜色通道的EFF矩阵F_R、F_G、F_B的一组过滤器来参数化的。这些过滤器依赖于所使用的透镜和照相机。

假定等式（5）中的噪声是高斯噪声，那么从原理上讲，通过求解最小二乘问题，即，通过关于x最小化

未知的全色图像x可以从测量到的原始图像y恢复。但是，由EFF框架参数化的PSF只是真实PSF的近似，因此，当最小化L₂范数时，这导致次优的结果。因此，最小化L₁范数||y-Ax||₁是优选的，这导致更好的结果。

处理实际照片的一个挑战是像素可能是饱和的，即，它们的真实的值由于有限的动态范围而被截断了。因而，测量到的被截断的像素的值不与模糊的物理模型一致。通过只对不饱和的像素求和，饱和的像素可以被排除在数据保真项||y-Ax||₁之外。这个项对应于隐式底层概率模型的可能性项（或者数据拟合）。但是，因为本发明性方法试图从单个原始图像估计三个颜色通道，即，存在三倍于观察的未知数，所以我们去模糊问题的提法是不妥当的。为了使其正则化，可以包括关于自然图像的现有知识：已经显示图像梯度近似地遵循超拉普拉斯分布。通过把形式为

的正则化项加到目标函数，这可以结合到优化问题中。这种正则化的效果是罚（penalize）强梯度并且因此来平滑图像。遵循Farsiu等（S.Farsiu、M.Elad和P.Milanfar等人在2006年IEEE Trans.Image Process第15卷第1期141-159页上所写的Multiframe demosaicing and super-resolution of color image.），在应用正则化之前，RGB图像可以变换到亮度/色度彩色空间中（在本发明的本实施方式中，使用YUV）。这允许在色度空间中更强地正则化，而在亮度中较弱。人眼对亮度比对色度的差别更加敏感，即，视觉上令人愉快的结果必须在亮度通道中是清晰的。从RGB到YUV的变换是利用适当选择的矩阵C的矩阵矢量乘法，

$\left[\begin{array}{c}{x}_Y\\ x_U\\ x_V\end{array}\right]=C\left[\begin{array}{c}{x}_R\\ x_G\\ x_B\end{array}\right]$

对于x_Y、x_U和x_V，根据本发明的一种可能的组合目标函数可以写成 ${\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+{\mathrm{\&alpha;}\left|\right|\&dtri;x_Y\left|\right|}_{\mathrm{\&gamma;}}^{\mathrm{\&gamma;}}+\mathrm{\&beta;}{\left|\right|\&dtri;x_U\left|\right|}_{\mathrm{\&gamma;}}^{\mathrm{\&gamma;}}+\mathrm{\&beta;}{\left|\right|{\&dtri;x}_V\left|\right|}_{\mathrm{\&gamma;}}^{\mathrm{\&gamma;}}+n{\mathrm{\&sigma;}}^2({\left|\right|x_R\left|\right|}_2^2/4+{\left|\right|x_G\left|\right|}_2^2/2+{\left|\right|x_B\left|\right|}_2^2/4)---\left(2\right)$

在本发明人的模拟实验中，好的结果是通过设置α=10^-4、β=10^-3及γ=0.65和σ=10^-3来获得的。对于实际的图像，用于α和β的最优值小一个为十的系数。

目标函数可以对本发明的装备采用Krishnan和Fergus的（D.Krishnan和R.Fergus.于2009年在Advances in Neural Inform.Process.Syst.上发表的Fast image deconvolution using hyper-Laplacian priors.）方法、在凸起和非凸起相之间交替来关于x最小化，其中非凸起相是通过查找表加速的。

为了提高执行速度，利用本发明性目标函数稍做修改的版本，本发明还可以用于傅立叶空间中的直接解卷积。由于去马赛克算子在傅立叶空间中不是对角线的，因此可以单独地对每个已经去马赛克的颜色通道工作并且求解问题

${\left|\right|y-\mathrm{Bx}\left|\right|}_2^2+\mathrm{\&alpha;}{\left|\right|{\&dtri;x}_Y\left|\right|}_{\mathrm{\&gamma;}}^{\mathrm{\&gamma;}}+{\mathrm{n\&sigma;}}^2{\left|\right|x\left|\right|}_2^2---\left(13\right)$

但是，B的倒置是必要的。利用上述关系，B在本发明性框架中的应用可以写成

$\underset{r}{\mathrm{\&Sigma;}}{L}_r^T{F}^H\mathrm{Diag}\left({\mathrm{FPf}}^{\left(r\right)}\right){\mathrm{FK}}_r\mathrm{Diag}\left(w^{\left(r\right)}\right)---\left(14\right)$

在这个对所有片的求和当中，矩阵P零填充每一片，矩阵K_r和L_r是被裁剪的矩阵。F应用离散傅立叶变换。这个表达式可以近似地倒置成

其中|z|和分别指示关于输入的绝对值与共轭复数。矩阵R正则化结果，例如，离散拉普拉斯算子。权重N是通过对固定图像应用倒置获得的并且对于除去由于倒置窗口而产生的假象是必要的。尽管损失了少量的图像质量，但是对于21.1兆像素图像的运行时间只有2分钟。

撇开衍射效应（例如，通过确保像素尺寸大于艾里斑（Airydisk）），点光源应当只影响数字照相机的成像传感器上的单个像素。但是，这只有在数字照相机是理想的光学***时才会发生。实践当中，以上所讨论的各种透镜像差将会在成像传感器的更大区域上展开点光源。这种局部模式特征化PSF，因此，通过记录跨图像平面的这些模式，不平稳卷积的过滤器可以如上所述地进行设置。

为了自动化设置点扩散函数（PSF）所需的测量，照相机可以安装到具有两个旋转自由度的电机驱动的平台上，如图1中所示。透镜测量过程可以在黑暗的房间中通过远程改变照相机朝点光源的角度来进行（在12米的距离内发射光通过100μm的光圈的汽灯），使得在后续的曝光中，光点可以在传感器上等距的位置被捕捉。在实验中，本发明人对EFF框架使用18乘27的支撑点网格。模糊内核是通过平均点光源的三个图像并阈值化噪声来记录的。

II.根据本发明的第二种实施方式，上述思想也可以应用到相机抖动的校正。相机抖动是手持式、更长曝光照片的一个常见问题，尤其发生在低照明情况下，例如建筑物内部。除了例如平移摄影的少数例外，相机抖动都是不期望的，因为它常常破坏细节并模糊图像。特定相机抖动的影响可以通过对清晰图像，即，利用三脚架记录的图像，的线性变换来描述。为了简化而把图像指示为列矢量，所记录的模糊图像y可以写成清晰图像x的线性变换，即，写成y=Ax，其中A是描述相机抖动的未知矩阵。盲图像去模糊的任务是在只给定模糊图像y而没有A的情况下恢复x。

为了获得通用的图像去模糊方法，线性变换y=Ax可以由上述可以处理平滑变化的PSF的有效过滤器流来表示。

本发明性方法不是简单地把不同的PSF应用到不同的图像区域，而是为每个像素产生不同的PSF。其原因是通常片被选择成重叠至少50%，使得位于一个像素的PSF是几个过滤器的某种线性组合，其中权重选择成平滑地融入和融出过滤器，并且因此导致在每个像素的PSF都趋于不同。图3显示了小到3×3的PSF阵列，对应于p=9和九个重叠的片（底部一行的右面板），可以参数化精密模仿实际相机抖动（左列）的平滑变化的模糊（中间列）。

根据本发明用于校正相机抖动的单图像盲解卷积算法估计EFF变换的参数矢量a并且执行空间变化的非盲解卷积，以便获得校正后的图像。

更特别地，如果x是利用模糊图像y初始化的，那么参数化为EFF的线性变换A的估计可以如下执行：

在步骤X10，通过边缘保留双边过滤可以在图像x的平区域中除去噪声，并且通过冲击过滤过分强调边缘。为了遏制由于冲击过滤而增强的噪声，可以应用空间自适应梯度幅值阈值。

本发明的预测步骤类似于由Cho和Lee所给出的预测步骤。它避免了非线性的优化，如果被非线性过滤操作（即，冲击和双边过滤及梯度幅值阈值化）强调的图像特征必须由图像事先对x实现的话，那么非线性优化将是必需的。根据本发明本实施方式的方法也受益于这个窍门而且超参数可以像Cho和Lee那样精确地设置（对于关于非线性过滤操作的细节，见[2]）。但是，梯度阈值化是空间自适应地应用的，即，单独地对每一片应用。否则，有些区域中的大梯度可能完全擦掉该区域中较少纹理的梯度，从而导致那些区域中差的PSF估计。

在步骤X20，给定模糊图像y与对被预测的x的当前估计，只利用（在预处理效果中导致的）x的梯度图像并且加强邻接PSF之间的平滑度，可以更新对PSF的估计。

给定非线性过滤的图像x的阈值梯度图像作为预测步骤的输出，PSF估计最小化正则化之后的最小二乘成本函数，

$\underset{z}{\mathrm{\&Sigma;}}{\left|\right|{\&PartialD;}_{z}y-{A\&PartialD;}_{z}x\left|\right|}^2+\mathrm{\&lambda;}{\left|\right|a\left|\right|}^2+\mathrm{vg}\left(a\right)---\left(16\right)$

其中z的范围在集合{h,v,hh,υυ,hυ}之上，即，考虑y和x的一阶和二阶、水平和垂直导数。忽略零阶导数（即，图像x和y本身）具有像在Cho和Lee著作中所讨论的预处理效果。矩阵A也依赖于PSF a的矢量。

可以添加正则化项g(a)，这促进了邻接PSF之间的相似性，

$g\left(a\right)=\underset{r=0}{\overset{p-1}{\mathrm{\&Sigma;}}}\underset{s\&Element;N\left(\text{r}\right)}{\mathrm{\&Sigma;}}{\left|\right|a^{\left(r\right)}-a^{\left(s\right)}\left|\right|}^2---\left(17\right)$

其中，如果片r和s是邻居的话，有s∈N(r)。

在步骤X30中，可以可选地识别估计得不好的PSF的区域并且用邻接的PSF代替。

由于高频信息，即，图像细节，对于PSF估计是必需的，因此，对于具有较少结构化区域的图像（例如天空），可能无法在每个地方都估计出合理的PSF。这个问题源自于发现，尽管有些区域可能关于局部PSF信息不足，但是它可以看起来是模糊的，并且因此将需要解卷积。这些区域可以通过阈值化对应PSF的熵来识别。被拒绝的PSF可以用它们的邻接PSF的平均来代替。由于可能有其邻接PSF也被拒绝的区域，因此可以执行递归过程，该递归过程把被接受的PSF传播到被拒绝的PSF。

在步骤X40中，当前去模糊的图像x可以通过事先对梯度图像利用平滑以便最小化最小二乘成本函数来更新。

在Cho和Lee及Krishnan和Fergus的著作中，图像估计步骤都涉及直接解卷积，这对应于通过傅立叶域中的零填充PSF对模糊图像的简单的关于像素的分割。不幸的是，直接解卷积对于表示为EFF的线性变换通常是不存在的，因为它涉及对片的求和。但是，直接解卷积可以用某个正则化后的最小二乘成本函数的优化

来代替。

尽管在（i）中估计了线性变换，但是正则化器是关于梯度图像的Tikhonov，即，p=2。由于估计出的x随后在预测步骤中被处理，因此可以认为正则化在（i）的图像估计步骤中是冗余的。但是，由于对a的不充分估计，因此调整对于抑制振铃（ringing）是关键的。

图像估计步骤与[2]和[3]的主要区别在于线性变换A不再是卷积，而是由EFF框架实现的空间变化过滤器。

在本实施方式中，上述步骤对多量程图像金字塔中的每一量程重复七次。它们以大小为3×3像素的PSF和对应缩减采样的观察图像开始。对于增加和缩减采样，可以采用简单的线性插值方案。a中的结果PSF和在每一量程的结果图像x被增加采样并且初始化下一量程。这种迭代过程的最终结果是参数化空间变化线性变换的PSF。

在获得了形式为PSF阵列的用于线性变换的估计之后，交替（ii）的步骤事先利用自然图像统计对所记录的图像y执行空间变化非盲解卷积（就像[13]中一样）。为此，最近由Krishnan和Fergus建议的方法可以用来处理表示为EFF的线性变换。

在步骤X50中，估计事先利用稀疏性正则化的本征变量，该本征变量近似x的梯度。这可以有效地利用查找表来求解。

在步骤X60中，当前去模糊的图像x是通过在对前一步骤的本征变量罚梯度图像的欧几里德范数的同时最小化最小二乘成本函数来更新的。

III.在本发明的一种备选实施方式中，给出了用于校正数字图像的相机抖动的方法，如下：

令g是模糊的照片，我们要为其恢复清晰的版本f。EFF把不平稳的模糊近似为R个不同模糊的片之和，

$g=\underset{r=1}{\overset{R}{\mathrm{\&Sigma;}}}{a}^{\left(r\right)}*(w^{\left(r\right)}\&times;f)---\left(18\right)$

其中权重图像w^(r)具有与清晰图像f相同的大小，因此第r个片可以写成w^(r)×f，其中×表示关于像素的（Hadamard）乘积。具有第r个模糊内核a^(r)的平稳模糊（aka卷积）表示为a^(r)*。由于w^(r)≥0选择成只在f的一个小区域内非零，因此每一片的卷积可以有效地利用短期快速傅立叶变换（FFT）实现。如果片选择成具有充分重叠而且a^(r)是截然不同的，那么由EFF表示的模糊将在每个图像位置不同，从一个像素到另一个像素梯度变化。

为了把EFF的可能模糊限定到相机抖动，用于模糊内核a^(r)的基础可以利用单应来创建。但是，代替把单应应用到清晰图像f（就像在[25，6]中所进行的），本发明对单像素点的网络只应用一次所有可能的单应。于是，可能的相机抖动会是由于线性组合点网格的不同单应而产生的，见图1。应当指出，这些单应可以预先计算，而不需要知道模糊图像g。

更特别地，令p是德耳塔峰值的图像，其中该峰值确切地位于片的中心，而且p具有与清晰图像f相同的尺寸。片的中心是由对应权重图像w^(r)的支撑的中心确定的。令θ索引一组单应并且由H_θ来指示特定的单应。于是，点网格p的不同视图p_θ=H_θ(p)可以通过应用单应H_θ来生成，其中p_θ同样是与f有相同尺寸的图像。通过近似地把这些视图p_θ裁剪成局部模糊内核

每个片一个，可以获得用于局部模糊内核a^(r)的基础，现在它可以被限定到基础模糊内核的线性组合，

$a^{\left(r\right)}=\underset{\mathrm{\&theta;}}{\mathrm{\&Sigma;}}{\mathrm{\&mu;}}_{\mathrm{\&theta;}}{b}_{\mathrm{\&theta;}}^{\left(r\right)}---\left(19\right)$

其中，μ_θ为整体模糊确定对应单应的相关性，即，μ_θ不依赖片索引r。模糊内核基础的构造确保由a^(r)参数化的结果产生的整体模糊对应于可能的相机抖动。这种线性组合不对应于单个单应，而是线性组合了几个单应，这些单应构成可能的相机抖动。因此，根据本发明的用于相机抖动的快速正演模型可以写成：

在这里，引入了符号表示μ◇f来指示根据本发明的用于相机抖动的快速正演模型。这种模型利用参数μ_θ来参数化相机抖动，其中μ_θ在矢量μ中汇总。这些参数看起来只是线性的，这对于有效优化是关键的。同样，清晰图像f只对固定的μ看起来是线性的。因而，存在矩阵M和A，使得用于矢量值的正演模型可以表达为矩阵矢量乘法（MVM），

g=μ◇f=Mf=Aμ                  (21)

为了方便，下文中将使用这种符号表示；用于矢量值图像的所有公式都可以直接了当地一般化到矩阵值图像。

为了有效地评估正演模型，或者换句话说，为了对M或X计算MVM，可以首先利用模糊参数μ和模糊内核基础（利用等式（2））计算所有的a^(r)，然后运行以上具体描述的EFF的快速实现。类似地，我们可以对M^T和A^T获得MVM的快速实现。关于点网格p的单应计算可以预先计算，而且在更新模糊参数μ_θ之后和更新清晰图像估计的估计之后都不再需要。这个事实对于本发明性方法的快速运行时间是基本的。

现在以被相机抖动模糊的照片g开始，未知的清晰图像f可以在两个阶段中恢复：（i）对不平稳PSF的模糊估计估计，及（ii）利用适合（tailored to）不平稳模糊的非盲解卷积过程的清晰图像恢复。以下，将具体描述这两个阶段的单个步骤而且在适当的地方明确地包括超参数的值，确定所涉及的项之间的权重。这些值在本发明人的所有实验过程中都是固定的。

在第一阶段，本发明性方法试图在只给出模糊照片的情况下恢复在曝光过程中由照相机所采取的运动：

在步骤Y10中，为了减少模糊并增强图像质量，可以使用冲击与双边过滤的组合。模糊参数是通过找出把当前图像估计f变换成所记录的模糊图像g的最佳不平稳模糊来估计的。但是，一开始，当前的估计f可能甚至无法接近清晰图像，而且在几次迭代后，它可能仍然是稍微模糊的。为了加速模糊参数估计的收敛，预测步骤可以通过冲击过滤来强调f中的边缘并且通过双边过滤来降低平噪声区域的重要性。

在步骤Y20中，估计相机运动参数，这最佳地解释来自步骤Y10的预测图像的模糊图片。为了符号表示的简化，假定g、μ和f是矢量值的。模糊或者相机运动参数μ可以通过最小化下式来更新

其中g的离散导数可以用符号写成

即，

对于矩阵值的图像，考虑水平和垂直导数。此外，

指示预测步骤（i）的结果而M_s是只选择提供信息的边缘并且方便内核估计的加权掩码。特别地，它忽略了尺寸比局部内核小的结构，这可能会误导内核估计。

等式（22）中的项可以如下推动：如果假定加性高斯噪声n的话，第一个项与对数相似性

成比例。考虑g和μ◇f的导数带来了几个好处：首先，Shan等人的[20]已经显示这种具有图像导数的项有助于通过对边缘设置权重来减少振铃假象。其次，它减少了优化问题等式5的条件个数并由此导致更快的收敛。通过抑制μ中的高强度值，第二个被加数罚μ的L₂范数并有助于避免平凡解。第三个项

加强了μ的平滑度，并因此有利于相机运动空间中的连通性。

在步骤Y30中，经非盲解卷积来估计本征图像。在模糊估计阶段中重复更新的清晰图像估计f不需要完美地恢复真实的清晰图像。但是，它应当在交替的更新过程中，即，步骤Y10、Y20和Y30中，引导PSF估计。由于在这第一个阶段中花费了大部分计算时间，因此清晰图像更新步骤应当很快。这促进对清晰图像采用基于L₂的正则化项，虽然结果估计有可能显示出一些振铃而且有可能还有其它假象（这些假象在预测步骤中处理）。因而愿意关于f最小化

通过用傅立叶空间中关于像素的分割来代替f中的迭代优化，Cho和Lee获得了大的加速。尽管有已知的修复假象，但是它们显示这种非迭代的更新步骤足以引导PSF估计。我们把傅立叶空间中这种关于像素的分割称为直接解卷积（DD）并且为我们对相机抖动的快速正演模型提供类似的更新。

首先，可以修改在[8]中给出的矩阵表达式，以便获得在部分3中介绍的用于M的明确的表达式，

其中，对于变化的θ，B^(r)是具有列矢量

的矩阵，即，

矩阵C_r和E_r近似地选择成裁剪矩阵，F是离散傅立叶变换矩阵，而Z_a是零填充的矩阵。此外，Diag(υ)指示具有沿其对角线的矢量υ的对角矩阵。

用于图像估计的直接更新步骤的基本思想是组合傅立叶空间中关于片的、关于像素的分割与重新加权和边缘衰落，来最小化振铃假象。以下表达式可以用于近似“倒置”本发明性正演模型g=Mf：

$f\&ap;\mathrm{Diag}\left(v\right)\underset{r}{\mathrm{\&Sigma;}}\mathrm{Diag}{\left(w^{\left(r\right)}\right)}^{1/2}{C}_r^T{F}^H\frac{\stackrel{\&OverBar;}{F{Z}_a{B}^{\left(r\right)}\mathrm{\&mu;}}\&times;\left({\mathrm{FE}}_r\mathrm{Diag}{\left(w^{\left(r\right)}\right)}^{1/2}g\right)}{{\left|{\mathrm{FZ}}_a{B}^{\left(r\right)}\mathrm{\&mu;}\right|}^2+\frac{1}{2}{\left|{\mathrm{FZ}}_{l}l\right|}^2}---\left(25\right)$

其中，用于具有复数输入的z的|z|计算关于输入的绝对值，而

计算关于输入的共轭复数。平方根是关于像素来取的。项Diag(υ)是某个附加的权重，我们在下一段中实验证明。这一部分必须实现成关于像素的。这一部分的分母中的项|FZ_ll|²源自等式（6）中利用对应于离散拉普拉斯算子的l=[-1,2,-1]^T的正则化。

给定模糊照片g和模糊参数μ，等式（24）中的更新等式近似真正的清晰图像f，而且可以通过从右向左读它来有效地实现。等式（24），即，直接解卷积，而没有附加的加权项（即，υ=1）近似真实图像，但是也显露了由于开窗造成的假象。通过应用附加的加权项υ，这些假象可以被有效抑制。加权υ可以通过应用等式（24）来计算，而对于与模糊图像g有相同尺寸的固定图像没有附加的加权项。解卷积的固定图像显露了相同的假象。通过关于像素取其倒置，它充当校正的加权项，这能够除去由于开窗造成的大部分假象并且同时计算也是很快的。

为了避免局部最小值并且使第一阶段中的模糊估计健壮，可以对几个量程采用由粗到细的迭代。开始，所记录的图像g的分辨率可以减小而且执行模糊估计阶段。然后，较低分辨率的模糊估计初始化下一量程，等等，直到所记录图像的全分辨率。在最粗的量程，未知的清晰图像f可以通过模糊图像g的缩减采样版本初始化（在零和一之间缩放）。对于固定的量程，在图6中可视化步骤（i）-（iii），并且将在下面具体描述。

在估计出并固定模糊参数μ之后，最终的清晰图像f可以通过在清晰图像更新步骤之前用之前基于梯度图像的稀疏性的自然图像代替L₂图像来恢复（例如，Fergus等人的[4]）即，最小化该项。

其中L₁项代表之前用于某个α≤1的自然图像。

为了最小化等式（9），Krishnan与Fergus的[10]方法可以修改成适用于不平稳情况下的平稳非盲解卷积：在引入辅助变量υ之后，交替地最小化f和υ中的

在f和υ中九次交替的更新过程中，权重2^t从1增加到256，其中t=0，1，...，8。选择a=2/3允许对υ中更新的解析公式，见Krishnan与Fergus的[10]。

GPU实现

以上具体描述的算法本身导致关于图形处理单元（GPU）的并行化。我们在PyCUDA中重新实现算法的所有步骤，PyCUDA是对NVIDIA的CUDA API的Python包装。为了评估加速，对Intel核心i5上的单核MATLAB的运行时间与Tesla C2050进行比较。表1显示去模糊逼真尺寸的实际图像在GPU上可以比在通常的（单核）处理器上快大约一百倍地执行。

IV.在本发明的另一种实施方式中，将显示过滤器如何可以应用到重构天文图像的问题。

假定在每个时间点t=1,2,...，通过用某个未知的模糊内核f_t模糊并添加噪声n_t，望远镜记录从共同底层图像x获得的模糊图像y_t。模糊图像{y₁,y₂,...}以流化方式到达，而且为了建模真实性，假定只存在一个大小为W（窗口尺寸）的用于保存它们的有限缓冲区。如果一种方法一次只访问一个图像的话，即，W=1，就说该方法是在线的，同而如果图像流终止并且所有观察到的（即，T个）图像都可用的话，就称该方法为批处理方法（W=T）。自然，当1<W<T时，就是混合过程。

为了符号表示的简化，该方法将关于表示为列矢量的一维图像来描述。但是，所述算法和结果可以很容易地关于两维图像重新用公式表示。事实上，本发明人所做的所有实验实际上都是在两维（矩阵）装备中执行的。

令每个所记录的图像y_t都是长度为l_y的矢量，而且未知的真实图像x具有长度l_x。x的大气模糊建模为具有非负模糊内核的非环形卷积（点扩散函数（PSF））f_t；这种卷积由f_t*x指示，而且f_t可以表示为长度l_f=l_x-l_y+1的列矢量。

为了流化来自望远镜的数据，每个观察到的图像y_t将由不同的和未知的PSF f_t来卷积。

非环形的卷积可以等价地（利用适当的零填充）写成矩阵矢量乘法，即，f_t*x=F_tx，其中ly×lx矩阵F_t只依赖于f_t。变量f_t和F_t在下文中将可以互换使用。对于有效的实现，F_t的具体形式允许经快速傅立叶变换（FFT）对矩阵矢量乘积F_tx的快速计算。矢量F_tx的单独成分是由(F_tx)_j来指示的。矢量的部分应当理解为关于组成部分的部分。

对于高强度的天文数据，y_t、f_t和x之间的关系可以修改为：

y_t=f_t*x+n_t                   (28)

其中噪声n_t假定为是高斯噪声，其平均数为0而且对角线协方差矩阵为　²I。假定固定但未知的σ，可以通过最小化以下损失来去模糊

$L_t(x,f_t;y_t)={\left|\right|y_t-F_{t}x\left|\right|}_2^2---\left(29\right)$

这个损失必须接受关于x和f_t的非负限定来最小化。尽管最小二乘模型看起来更加健壮，但是对于低强度数据，泊松模型被认为更加合适，

y_t～Poisson(f_t*x)                  (30)

由此，本发明将把其注意力单独地集中到（28）和（29）。但是，如果期望的话，本发明性在线重构方案可以扩展成处理（30）。

如果已经访问了所有的图像{y_t}_t≥1，那么理想地愿意最小化整体损失

$\begin{array}{c}\min \\ x,f_1,f_2,...,f_T\end{array}\underset{t=1}{\overset{T}{\mathrm{\&Sigma;}}}{L}_t(x,f_t;y_t)---\left(31\right)$

其中L_t(x,f_t;y_t)是在如由（29）给出的tth步骤引起的损失。但是，图像y_t是以尺寸为W的有限缓冲区的流化方式，同时最小化整体损失（31）需要保存所有的y_t。

由于x和f_t之间经卷积f_t*x的联合依赖性，整体损失（31）是非凸的。接下来，对于W=1的纯在线情况描述近似地解决这些难题的在线算法，应当指出，它可以很容易地一般化，来处理W>1的混合场景。

从每个观察到的图像y_t，需要估计PSF f_t及底层图像x。令x_t-1是真实图像的当前估计。现在，可以获得新的观察y_t并且可以更新真实图像的估计，同时估计PSF f_t。为此，我们可以最小化在时间步长t引入的损失。正式声明，以下非负限定的问题可以求解

$\begin{array}{c}\min \\ f_t\&GreaterEqual;0,x\&GreaterEqual;0\end{array}{L}_t(x,f_t;y_t)={\left|\right|y_t-F_{t}x\left|\right|}^2---\left(32\right)$

以便获得更新后的估计x_t和PSF f_t。

因而，通过在每个步长t只求解问题（32），我们获得了解决以上提到的第一个问题的在线过程。但是，第二个问题仍然存在，因为甚至连问题（32）都是非凸的，由此通常对其全局最优的解决方案不能有效地找到。幸运的是，目标函数表现充分良好，因为，如果其它保持固定的话，在每个变量中单独看它是凸起的。这个关键属性是用于解决（32）的简单交替最小化或者下行（descent）方案的要点。这种方案构成确保下行的迭代序列

即

$L_t(x_t^{k+1},f_t^{k+1};y_t)\&le;L_t(x_t^k,f_t^k\text{;}{y}_t)---\left(33\right)$

并且在一些弱假设下，它可以示为收敛到（32）的平稳点。

但是，太好地执行交替下行是要付出代价的，即，y_t的过度拟合。为了避免过度拟合，可以只执行交替下行的几次迭代，而且，就象我们的实验一样，将显示一或两次迭代就足以得到好的结果。现在要回答的问题是如何执行这种迭代，使得可以满足下行条件（6）。我们建议以下两个步骤：

给定模糊图像的输入流{y_t}_t≥1，以下步骤是为恢复真实的底层图像x而建议的。附带地还获得了对单个模糊内核f_t的估计。

现在，该方法执行以下两步：

1.f_t=argmin_f≥0||y_t-Fx_t-1||²

2.x_t=argmin_x≥0g_t(x,x_t-1)

步骤1是非负最小二乘（NNLS）问题而且它可以利用已知的方法，优选地是[3]的LBFGS-B算法，有效地解决。

在步骤2中，g_t(x,x_t-1)是必须满足 $\&ForAll;x,\tilde{x}:g_t(x,\tilde{x})\&GreaterEqual;g_t(x,x)=L_t(x,f_t;y_t).$ 的辅助函数。

由于通过定义g_t(x_t-1;x_t-1)≥g_t(x_t,x_t-1)≥g_t(x_t;x_t)，因此这个条件确保步骤2中x_t的选择只会减小我们的损失函数。

很容易看到，对于二次损失，这种函数是由 $g_t(x,\tilde{x})=y_t^T{y}_t-{2y}_t^T{F}_{t}x+{\tilde{x}}^T{F}_{\&prime;t}^T{F}_t\left(\frac{x\&times;x}{\tilde{x}}\right)$ 给出的，

其中×指示关于元素的乘积而且两个矢量的分割应当理解为是关于组成部分的。对步骤2的解答是通过以封闭的形式对x_t求解

获得的，这产生以下的倍增更新，

$x_t=x_{t-1}\&times;\frac{F_t^T{y}_t}{F_t^T{F}_t{x}_{t-1}}---\left(34\right)$

这种更新考虑对x的非负限定，因为所有涉及的量都是非负的。

利用NNLS求解f_t的原因是对于x我们仅仅是利用辅助函数的技术执行下行。原因如下：

未知的图像x对于所有观察到的图像是公共的，尽管PSF对于每个图像都是不同的。因此，x和f_t需要区别对待。为此，给定图像x_t-1的当前估计和新的输入图像y_t，损失L_t可以关于f_t最小化，而对于x仅仅是关于L_t下行。这么做基本上确保f_t很好地拟合，而且x_t与x_t-1不会有太大区别。

在执行步骤1和2之后，很清楚f_t和x的估计已经得到了改进，因为步骤1对f最小化，而在步骤2中辅助函数g_t(x;x_t-1)确保我们进行了下行。因而，下行保证（6）保留，并且对于任意x_t-1和f采取形式

L_t(x_t,f_t;y_t)≤L_t(x_t-1,f;y_t)             (35)

根据以上讨论用于天文成像的在线盲解卷积方法在图7中用伪码示出。

对于泊松情况，本发明性方法采取形式：

1.f_t=argmin_f≥0KL(y_t;f*x_t-1)

2.x_t=argmin_x≥0g_t(x,x_t-1)

其中KL指示一般化的Kullback-Leibler发散，而g_t(x,x_t-1)是用于关联损失的适当辅助函数。步骤1同样可以利用LBFGS-B方法来求解，而[17]的EM导数产生

$g_t(x,\tilde{x})=1^T{F}_{t}x-y_t^T\mathrm{\&Lambda;}\log \left(x\right)+c$

其中c是独立于x的常量，而Λ是矩阵

以上所述的方法可以适当地修改成适合这些改变。

给定一系列模糊图像的话，以上所述的基本算法可以恢复底层的清楚图像。现在，该基本方法扩展成结合超级分辨率和饱和度校正。本发明性方法很简单，但是非常有效，因为它利用了丰富的数据但是没有牺牲效率。

由于我们只对单个图像x感兴趣但是却有几个观察y_t，虽然有模糊，但是我们可能有可能推导出一个超级分辨率的图像-假定分辨率的变化结合到正演模型中的话。为此，调整大小矩阵可以定义为

$D_m^n=(I_m\&CircleTimes;1_n^T)(I_n\&CircleTimes;1_m)/m,$

其中I_m是m×m的单位矩阵，1_m是单位矩阵的m维列矢量，而

指示Kronecker乘积。矩阵

把长度为n的的矢量υ变换成长度为m的矢量。υ的输入之和

υ被保留。这对于图像来说是一个有利的属性，因为所观察到的光子的数量不应当依赖于分辨率。即使m和n彼此不是倍数，但是

将适当地插值。

令n是y的长度。对于k倍的超级分辨率，x和f可选择成足够大，使得f*x具有长度kn。现在，可以使用图像模型

这将导致修改后的更新步骤：

正缩放因子k不限于是整数。只有kn需要是整数。

饱和的（过度曝光的）像素会不利地影响图像的恢复，在天文成像中尤其是这样，其中人可能想捕捉暗淡的星体连同亮几倍的星体。现实的解卷积方法应当能够处理饱和的像素，即，那些命中（或者靠近）最大可能像素值的像素。

处理饱和像素的一种合理途径是给它们加权。由于每个帧y_t会具有饱和的不同像素（不同的帧是以不同方式对准的），因此可以在每次迭代时检查哪些像素是饱和的。为了忽略这些像素，可以定义加权矩阵，

C=Diag([y_t＜ρ])

其中ρ指示最大像素强度，而Iverson括号[·]应用是关于组成部分的。于是，通过用加权的范数

代替欧几里德范数，更新可以忽略饱和像素来写。这种方法等效于从概率模型除去饱和的像素。

根据本发明，x中饱和的像素在大部分帧y_t中都可以恢复。由于PSF f_t，对应于x中这种像素的光子已经在每个帧y_t中跨整个像素集合散开了。因而，如果不是y_t中所有这些像素都是饱和的，那么x中对应像素的真实值将被恢复。

本发明所考虑的另一种应用是星体亮度的估计（光度测定），为此需要线性传感器响应（为了本发明的目的，所使用的CCD传感器可以假定为是线性的）。于是，强度计数可以翻译成星体的幅值。很清楚，这么做对于饱和CCD的星体不是直接可能的。但是，如果像素不是饱和的话，就可以使用用于具有饱和校正的解卷积的建议方法来重构我们可能已经记录的光子计数（图像强度）。这些计数可以转换成星体的幅值。对于四边形星团，获得了令人鼓舞的结果-见表1。

第二种应用解决磁共振成像（MRI）中对象运动的常见问题。MRI是在临床前研究中所使用的用于可视化人体和动物内部结构与功能的医学成像模型。与计算断层扫描（CT）相比，MRI提供了不同软组织之间的更大对比度，这使得它在神经***（脑），肌肉骨骼，心血管，肿瘤（癌症）成像中尤其有用。

图6示出了在临床前研究中用于对比度MRI的老鼠胸腔的两个图像序列的典型帧。这些序列对应于不同高度的两个轴向切片。二者都是利用Bruker的7次Tesla ClinScan取得的并且包含200个帧，每个帧的大小是128×128像素。如可以从这些帧看到的，对象运动产生大的模糊并且导致图像质量的显著损失。除全局对象运动之外，心跳也造成局部失真。全局和局部失真都可以通过用于空间变化过滤的本发明来描述。本发明人应用具有大小为20×20像素的4×4PSF的空间变化盲解卷积算法（选择大小为64×64像素的Bartlett-Hanning窗口，具有50%的重叠）。对于内核估计，可以采用附加的Tikhonov正则化。图6示出了我们方法的估计图像。更有趣的是，我们的方法可以用于有效的运动校正，因为在每个时间步长所估计的对象图像都保留在相同的位置。通过对内核施加能量限定，所估计出的PSF不仅给出关于对象运动的信息，而且还给出关于强度变化的信息，这是在对比度或功能性MRI中主要兴趣所在。可以看到，本发明性方法恢复了更多的图像细节。

Claims (22)

1.一种用于从由图像获取设备获取的图像（y）重构真实图像（x）的计算机实现方法，其中图像获取设备具有透镜，该方法包括步骤：

-获得数字形式的所获取图像（y）；

-基于数字形式的所获取图像和透镜的点扩散函数（A）来估计真实图像（x）；及

-输出所估计的真实图像（x）。

2.如权利要求1所述的方法，其中透镜的点扩散函数（A）是空间变化点扩散函数（A）。

3.如权利要求2所述的方法，其中真实图像（x）的每个颜色通道是通过单独的空间变化点扩散函数（FR、FG、FB）来估计的。

4.如权利要求2所述的方法，其中估计真实图像（x）还对渐晕进行校正。

5.如权利要求3所述的方法，其中估计步骤还基于图像获取设备的颜色过滤器阵列（CFA）的模型。

6.如权利要求2所述的方法，其中透镜的点扩散函数（A）是通过利用点光源自动测量透镜来确定的。

7.如权利要求2所述的方法，其中透镜是复合透镜。

8.如权利要求1所述的方法，其中图像获取设备是电子显微镜。

9.一种用于从由图像获取设备获取的图像（_y）重构真实图像（x）的设备，其中图像获取设备具有透镜，该设备包括：

-用于获得数字形式的所获取图像（y）的装置；

-用于基于数字形式的所获取图像和透镜的点扩散函数（A）来估计真实图像（x）的装置；及

-用于输出所估计的真实图像（x）的装置。

10.如权利要求9所述的设备，其中点扩散函数（A）是空间变化点扩散函数（A）。

11.一种用于从所观察到的数字图像序列（y₁,…,y_T）恢复数字图像（x）的计算机实现方法，该方法包括步骤：

-获得所观察到的数字图像（y_t）；

-基于所观察到的图像（y_t）估计点扩散函数（f_t）；

-基于所估计的点扩散函数（f_t）和所观察到的图像（y_t）来估计被恢复的数字图像（x）；

-重复以上步骤。

12.如权利要求11所述的方法，其中被恢复的数字图像（x）具有比所观察到的数字图像（y_t）更高的分辨率。

13.如权利要求11或12所述的方法，其中所观察到的数字图像中的饱和像素至少部分地被忽略。

14.如权利要求11所述的方法，应用到所观察到的数字图像序列（y₁,…,y_T）的单独颜色通道。

15.如权利要求11所述的方法，其中被恢复的图像（x）的全部颜色通道被同时估计。

16.如权利要求11所述的方法，其中：

-所观察到的图像被分成片；

-多个点扩散函数（f_t）被联合估计，每个片一个点扩散函数（f_t）；及

-被恢复的数字图像（x）是基于联合估计的点扩散函数（f_t）与所观察到的图像来估计的。

17.如权利要求16所述的方法，还包括步骤：

-通过所观察到的数字图像（y_t）的双边过滤除去噪声；及

-通过所观察到的数字图像（y_t）的冲击过滤强调边缘。

18.如权利要求11所述的方法，其中该方法被在线执行。

19.如权利要求11所述的方法，其中所观察到的数字图像序列（y₁,…,y_T）是天文图像序列。

20.如权利要求11所述的方法，其中所观察到的数字图像序列（y₁,…,y_T）是磁共振图像（MRI）序列。

21.一种用于从所观察到的数字图像序列（y₁,…,y_T）恢复数字图像（x）的设备，该设备包括：

-用于获得所观察到的数字图像（y_t）的模块；

-用于基于所观察到的图像（y_t）估计点扩散函数（f_t）的模块；

-用于基于所估计的点扩散函数（f_t）和所观察到的图像（y_t）来估计被恢复的数字图像（x）的模块；

其中

该设备适于执行如权利要求11所述的方法。

22.一种计算机可读介质，存储如权利要求11所述的方法恢复的数字图像。

Images