在使用相机进行拍摄的过程中, 由于相机相对于被拍摄对象间存在相对运动, 经常会造成图像模糊.在实际应用中, 相机的运动形式多种多样, 不同的运动形式产生的图像运动模糊机理也不尽相同, 尽管在视觉上对模糊的感觉是相同的.例如, 传统的(基于图像平移运动的) 卷积模型与以前向运动为主的相机运动模型之间存在很大的差别:前者相当于图像中各个像素点以相同的方向、速度进行运动, 表现出一致性;而后者相当于图像中各个像素点以不同的方向、速度进行运动, 表现出“非一致”性, 致使各个像素点的模糊是“非一致”的.

去除图像模糊是一个典型的逆问题, 其相应的正问题为图像的模糊过程.对正问题的正确理解, 是求解逆问题的关键之一.

为得到有效去除运动模糊的算法, 需要在分析其基本原理的基础上建立正确的图像运动模糊模型.

本文焦点放在一个具有典型代表性的运动模糊情况——相机以前向平移方式运动的情况, 即, 前向运动.该运动导致图像中所有像素点有各自不同的移动方向和速度.对于这种非一致性的图像模糊, 称为前向运动模糊.其他运动模糊, 如平移、抖动等常见的运动模糊均属一致性运动模糊, 是非一致性运动模糊的简化情况.

对于图像的前向运动模糊, 目前的方法都不能解决这一问题.我们的研究也表明:本文给出的模型并非对以前研究的简单扩展, 而是针对前向运动模糊成像的一种恢复方法.实际上, 在计算机视觉、近景摄影测量等领域, 也常遇到这种成像环境.

1 相关工作

实际得到的图像可以看作是“理想”图像经过某种退化过程所得到的结果, 称为退化图像.图像复原就是要通过退化图像估计出理想图像^[1].长久以来, 不管是静态模糊复原模型(例如TV(total variation)模型^[2]), 还是一些经典的去除运动模糊的模型(例如Levin^[3]和Cho等人^[4]的工作), 图像的模糊过程都被描述为一个卷积过程.这样做的原因是模型比较简单, 易于分析.事实上, 模型本身就是实际情况经过理想化假设而得到的近似结果.对于静态模糊(例如图像失焦、Gauss模糊等), 使用卷积模型所表现出的缺陷并不十分明显, 但是对于运动模糊, 卷积模型表现出较大的局限性.这是由于相机的运动(例如抖动、前向运动等)所引起的图像模糊并不是一个卷积过程, 图像上不同位置的像素点的模糊过程和方式不同.因此在前向运动成像条件下, 要取得好的图像复原效果, 需要了解图像运动模糊的正过程.

对非卷积模型, 也称为非均匀模糊, 已经有很多研究成果.如天文图像, 由于大气湍流等因素的影响, 即使是对于静态模糊, 光斑的大小和形状也会随着像素点位置的不同而发生很大的变化.Fusco等人^[5]对天文图像中的不一致模糊进行了深入研究.Whyte等人^[6]指出, 相机抖动所引起的图像模糊是非均匀的(他们用相机的转动模型来描述相机的抖动).围绕着图像的非均匀模糊这个热点研究问题, 提出了许多不同的方法.例如, Portz等人^[7]提出了一种有效的光流估计方法, Gupta等人^[8]提出了“运动密度函数”, 并尝试用它来实现图像去模糊问题. Hirsch等人^[9]提出了一种快速算法, 用于去除相机抖动所引起的非均匀模糊.特别值得提出的是, Zheng等人^[10]针对相机前向平动所引起的运动模糊进行了研究, 给出了解决方法.然而, Zheng处理的目标物体是在一个平面上, 因此可以忽略景深的影响, 相对于需要考虑景深的情况, 在难度上大为减小.

上面提到的这些方法尽管关注重点不同, 但其基本思路依然使用卷积模型来估计模糊核, 只不过是先将图像划分为一系列局部小区域, 对不同的区域使用不同的模糊核.因此, 这些方法对模糊核的估计速度较慢.有的方法, 如Zheng等人^[10]的工作, 还需要人来参与.

本文从泛函运动的角度给出广义运动模糊模型.在此基础上, 从光流场的角度出发, 建立了一种前向运动模糊模型.基本出发点是, 对于前向运动运动模糊, 模糊核所描述的是像素点的运动轨迹, 而光流所描述的是像素点的运动速度^{[11, 12]}.如果曝光时间很短, 可以近似认为像素点的运动轨迹等于运动速度乘以曝光时间.因此, 可以通过估计光流来计算前向运动模糊核.

光流估计是机器视觉中的一个经典问题, 已经有很多成熟算法^{[5, 11, 12, 13]}.为提高光流估计的速度, 可以先估计出相机的运动速度和景深(即无源导航^{[14, 15]}), 然后计算出光流, 并最终生成前向运动模糊核.这一过程就是本文提出的“三步法”:无源导航计算光流生成前向运动模糊核.

2 广义运动模糊和广义运动模糊核

相机的运动会造成像平面上图像亮度模式的变化, 而图像中像素点的灰度值取决于像平面上相应位置的图像感知器对曝光时间Dt内的光子数的积累量的测量结果^[11].因此, 这个问题可以分解为两个子问题:

1) 图像亮度模式的变化如何造成图像的模糊?

2) 相机的前向运动将使得图像亮度模式如何发生变化?

本节主要解决第1个问题, 下一节解决第2个问题, 从而建立起“广义运动”所对应的“前向运动”模糊模型.

2.1 广义运动模糊

图像是场景中物体在像平面上的二维投影, 图像中的某一个像素点对应于场景中的某一个物体“小块”.由于场景和相机之间的相对运动, 在不同时刻, 这个“小块”将被投影到像平面上的不同位置, 在相机快门曝光的瞬间Dt=T₀-0内, 这个“小块”将被投影到不同细分像平面上的不同位置(所谓细分像平面, 是假设的一种没有模糊的理想像平面), 从而引起图像亮度模式的变化.也就是说, 细分像平面上某一固定位置的亮度将随着时间而发生变化, 其中, T₀为曝光时间.

由于相机的运动, 场景中某一个“小块”的像在像平面内形成一条曲线γ(x₁, y₁, t)={x(t), y(t)^T}, 如图 1所示, 原来的x-y平面扩展为x-y-t空间.其中, (x(0) , y(0) )^T=(x₁, y₁)^T表示曲线的起始点位置.为确定亮度的变化, 需要确定场景中某一个“小块”的像在细分像平面上如何运动, 即, 确定“小块”的像的运动轨迹函数γ(x₁, y₁, t);同时, 也要知道场景中某一个“小块”的像在细分像平面上运动过程中其亮度值如何变化?也就是要知道 函数I(γ(x₁, y₁, t))的具体形式.

2.2 广义运动模糊核

运动模糊图像I_b(x, y)是图像序列I(x, y, t)的线性叠加结果(如图 1所示), 即

${{I}_{b}}(x,y)=\frac{1}{{{T}_{0}}}\int_{0}^{{{T}_{0}}}{I(x,y,t)\text{d}t}$

(1)

图像序列I(x, y, t)是由图像序列中的第一张图像I₀(x₁, y₁)(通过像素点的运动而)生成的, 因此, 模糊图像也可以看作是I₀(x₁, y₁)中各个像素点的灰度值的线性叠加结果.于是, 公式(1) 可以被进一步写为

${{I}_{b}}(x,y)=\iint{w(x,y,{{x}_{1}},{{y}_{1}}){{I}_{0}}({{x}_{1}},{{y}_{1}})\text{d}{{x}_{1}}\text{d}{{y}_{1}}}$

(2)

本文中, 公式(2) 被称为广义运动模糊模型, 核函数w(x, y, x₁, y₁)的具体形式为(证明参见附录)

$w(x,y,{{x}_{1}},{{y}_{1}})=\frac{1}{{{T}_{0}}}\frac{T(l({{x}_{1}},{{y}_{1}},t)\in A(x,y))}{\delta x\delta y}$

(3)

其中, A(x, y)表示以(x, y)^T为中心的像素点区域;而dx和dy分别表示该像素区域的长和宽;T(γ(x₁, y₁, t)A(x, y))表示曲线γ(x₁, y₁, t)在像素区域A(x, y)中的停留时间.

广义运动模糊模型代表了运动模糊的普遍模式, 可用以描述所有运动模糊.例如, 当所有的曲线γ(x₁, y₁, t)具有完全相同的形状时, 如果将曲线γ(x₁, y₁, t)的起始点位置(x₁, y₁)^T和像素区域A(x, y)的中心位置(x, y)^T同时移动(x₀, y₀)^T, 曲线在像素区域中的停留时间并不发生变化.不失一般性, 令(x₁, y₁)^T=(x₀, y₀)^T, 于是, 广义模糊核函数为

$w\left( x,y,{{x}_{1}},{{y}_{1}} \right)=w(x-{{x}_{1}},y-{{y}_{1}},0,0)=k(x-{{x}_{1}},y-{{y}_{1}})$

(4)

而广义运动模糊模型公式(2) 就变为卷积模型为

${{I}_{b}}(x,y)=\iint{k(x-{{x}_{1}},y-{{y}_{1}}){{I}_{0}}({{x}_{1}},{{y}_{1}})\text{d}{{x}_{1}}\text{d}{{y}_{1}}}$

(5)

对于绝大多数的运动模糊情况, 曲线γ(x₁, y₁, t)的形状并不相同, 因此, 卷积模型公式(5) 不成立.本文将依据广义运动模糊模型得到前向运动模糊模型.

3 无源导航算法

公式(3) 中的曲线γ(x₁, y₁, t)的具体表达式为

$l({{x}_{1}},{{y}_{1}},t):\left\{ \begin{array}{*{35}{l}} x(t)={{x}_{1}}+u(x,y)\cdot t \\ y(t)={{y}_{1}}+v(x,y)\cdot t \\ \end{array} \right.(0\le t\le {{T}_{0}})$

(6)

其中, 像素点的速度(u, v)^T又被称为光流^{[11, 12]}.假设曝光时间T₀非常短, 于是, 公式(6) 可以近似为

$l({{x}_{1}},{{y}_{1}},t):\left\{ \begin{array}{*{35}{l}} x(t)={{x}_{1}}+u(x,y)\cdot t \\ y(t)={{y}_{1}}+v(x,y)\cdot t \\ \end{array} \right.(0\le t\le {{T}_{0}})$

(7)

因此, 计算曲线γ(x₁, y₁, t)等价于估计光流(u, v)^T.用(U, V, W)^T表示相机沿着三维场景空间X轴、Y轴和Z轴的运动速度, 根据透视投影模型^[11], 场景中的一点(X, Y, Z)^T被投影成像平面上的一点(x, y)^T, 即

$x=\frac{f}{Z}X;y=\frac{f}{Z}Y$

(8)

而光流(u, v)^T和相机速度(U, V, W)^T之间的关系满足^[14]:

$u=\frac{-Uf+xW}{Z};v=\frac{-Vf+yW}{Z}$

(9)

若W≠0, 那么(u, v)^T会随着(x, y)^T的不同而发生变化(如图 2所示), 因此, 对于图 2的模糊并不是一个卷积过程.

直接估计(u, v)^T很复杂, 并且精度不高^[12].对于相机前向运动的情况, 可以先估计(U, V, W)^T和Z(x, y);然后, 通过公式(9) 计算光流场(u, v)^T.在机器视觉领域, 估计参数(U, V, W)^T和Z(x, y)也被称为无源导航^{[14, 15]}.

3.1 估计相机运动速度

在亮度不变的假设条件下, 光流约束方程为^[12]

${{E}_{x}}u+{{E}_{y}}v+{{E}_{t}}=0$

(10)

其中, (E_x, E_y)^T表示图像的亮度梯度, 而E_t表示像素点亮度的时间变化率.将公式(9) 带入公式(10) 得到:

$(-Uf+xW){{E}_{x}}+(-Vf+yW){{E}_{y}}+Z{{E}_{t}}=0$

(11)

公式(11) 中同时包含相机的运动速度(U, V, W)^T和场景的深度Z(x, y), 直接求解较为困难, 而无源导航^[15]是求解该问题的一种快速算法, 其基本过程为:

1) 选择点集Ω={(x, y)^T|(E_x, E_y)^T≠(0, 0) ^T并且E_t=0}.

2) 根据点集Ω, 使用最小二乘法^[16]求出U和V(令W=1) .

3) 选择点集Φ={(x, y)^T|(E_x, E_y)^T≠(0, 0) ^T并且E_t≠0}.

4) 根据点集Φ和步骤2) 中求出的U和V, 使用最小二乘法^[16]求出Z(x, y)(令W=1) .

我们将详细讨论步骤2) 和步骤4) 的具体实现过程.首先, 对于点集Ω, 公式(11) 进一步写为

$(-Uf+x){{E}_{x}}+(-Vf+y){{E}_{y}}=0$

(12)

最小二乘拟合通过最小化如下目标函数:

${{\Gamma }_{1}}(U,V)=\iint_{\Omega }{{{((-Uf+x){{E}_{x}}+(-Vf+y){{E}_{y}})}^{2}}\text{d}x\text{d}y}$

(13)

来求得相应的U和V.令公式(13) 关于U和V的偏导数等于0, 进一步得到:

$\left\{ \begin{array}{*{35}{l}} \left( \iint_{\Omega }{E_{x}^{2}\text{d}x\text{d}y} \right)U+\left( \iint_{\Omega }{{{E}_{x}}{{E}_{y}}\text{d}x\text{d}y} \right)V=\iint_{\Omega }{\frac{{{E}_{x}}}{f}(x{{E}_{x}}+y{{E}_{y}})\text{d}x\text{d}y} \\ \left( \iint_{\Omega }{{{E}_{x}}{{E}_{y}}\text{d}x\text{d}y} \right)U+\left( \iint_{\Omega }{E_{y}^{2}\text{d}x\text{d}y} \right)V=\iint_{\Omega }{\frac{{{E}_{y}}}{f}(x{{E}_{x}}+y{{E}_{y}})\text{d}x\text{d}y} \\ \end{array} \right.$

(14)

求解上式, 即可得到关于U和V的解析解.

3.2 估计场景的深度

在前向运动视频的实际应用中, 一个较为合理的假设是将场景视为平面, 如路面、护栏、整齐排列的电杆都位于平面范围^[17], 即Z=-n₁X-n₂Y+Z₀.根据公式(8) 整理得到:Z=Z₀/(1+n₁X/f+n₂Y/f), 于是, 公式(11) 整理为

$g\left( x,y \right)f+g\left( x,y \right)x{{n}_{1}}+g\left( x,y \right)y{{n}_{2}}+f{{E}_{t}}{{Z}_{0}}=0$

(15)

其中, g(x, y)=(-Uf+x)E_x+(-Vf+y)E_y.最小二乘拟合通过最小化如下目标函数:

$\iint_{\Phi }{{{\left( g(x,y)+g(x,y)\frac{x}{f}{{n}_{1}}+g(x,y)\frac{y}{f}{{n}_{2}}+{{E}_{t}}{{Z}_{0}} \right)}^{2}}\text{d}x\text{d}y}$

(16)

来求得相应的n₁, n₂和Z₀.令公式(16) 关于n₁, n₂和Z₀的偏导数等于0, 得到:

$\left\{ \begin{array}{*{35}{l}} a{{n}_{1}}+b{{n}_{2}}+c{{Z}_{0}}={{h}_{1}} \\ b{{n}_{1}}+d{{n}_{2}}+e{{Z}_{0}}={{h}_{2}} \\ c{{n}_{1}}+e{{n}_{2}}+g{{Z}_{0}}={{h}_{3}} \\ \end{array} \right.$

(17)

其中, 公式(17) 中的参数为

$\begin{align} & a=\iint_{\Phi }{{{g}^{2}}(x,y){{x}^{2}}\text{d}x\text{d}y,}\text{ }b=\iint_{\Phi }{{{g}^{2}}(x,y)xy\text{d}x\text{d}y,} \\ & c=f\iint_{\Phi }{g(x,y)x{{E}_{t}}\text{d}x\text{d}y}\text{, }d=\iint_{\Phi }{{{g}^{2}}(x,y){{y}^{2}}\text{d}x\text{d}y}, \\ & e=f\iint_{\Phi }{g(x,y)y{{E}_{t}}\text{d}x\text{d}y}\text{, }g={{f}^{2}}\iint_{\Phi }{E_{t}^{2}\text{d}x\text{d}y}, \\ & {{h}_{1}}=f\iint_{\Phi }{{{g}^{2}}(x,y)x\text{d}x\text{d}y}\text{, }{{h}_{2}}=f\iint_{\Phi }{{{g}^{2}}(x,y)y\text{d}x\text{d}y}, \\ & {{h}_{3}}={{f}^{2}}\iint_{\Phi }{g(x,y){{E}_{t}}\text{d}x\text{d}y}. \\ \end{align}$

通过求解公式(16) , 可以得到关于n₁, n₂和Z₀的解析解, 进而得到平面场景的表达Z(x, y).将U, V以及Z(x, y)带入公式(9) (令W=1) , 就得到了光流场(u, v)^T的具体形式.

4 生成前向运动模糊核

光流(u, v)^T是基于“前向运动”这个假设估计出来的, 因此, 我们将根据公式(7) 计算出的模糊核称为前向运动模糊核, 记为W={w_{i, j, k, γ}}.具体计算包括两个子过程.

1) 确定出从像素点A_{k, γ}出发的线段(即公式(7) )所“扫过”的像素点集合R_{k, γ}(称为离散线段);

2) 计算R_{k, γ}所包含的各个像素点(i, j)的权重系数w_{i, j, k, γ}.

4.1 生成离散线段

首先, 根据公式(9) 计算出像素区域A_{k, l}的中心点(x_{k, l}, y_{k, l})^T的运动速度(u_{k, l}, v_{k, l})^T.首先考虑(u_{k, l}, v_{k, l})^T的方向在0°和45°之间的情况(如图 3所示).根据公式(7) , 可以得到相应的直线段方程:

${{l}_{k}}_{,l}:(x-{{x}_{k}}_{,l}){{v}_{k}}_{,l}=(y-{{y}_{k}}_{,l}){{u}_{k}}_{,l}({{x}_{k}}_{,l}\le x\le {{x}_{k}}_{,l}+{{u}_{k}}_{,l}{{T}_{0}})$

(18)

将该直线段分别向上和向下平移1/2个像素点尺寸, 可以得到两条新的直线段:

$\left\{ \begin{align} & l_{k,l}^{+}:(x-{{x}_{k,l}}){{v}_{k,l}}=(y-{{y}_{k,l}}-1/2){{u}_{k,l}} \\ & l_{k,l}^{-}:(x-{{x}_{k,l}}){{v}_{k,l}}=(y-{{y}_{k,l}}+1/2){{u}_{k,l}} \\ \end{align} \right.,\text{ }{{x}_{k,l}}\le x\le {{x}_{k,l}}+{{u}_{k,l}}{{T}_{0}}$

(19)

本文中, 我们将这两条直线段之间的区域定义为像素点A_{k, l}所“扫过”的区域(如图 3所示).进一步, 需要确定出哪些像素点包含在该区域之中, 并且, 将所有这些被检测出的像素点保存成一个数组R_{k, l}={(i, j)}.

对于(u_{k, l}, v_{k, l})^T的方向在0°和45°之间的情况, 我们考虑像素区域A_{k, l}的右上角点(如图 3所示).找出R_{k, l}中所包含的所有像素点的右上角点是一个迭代过程.如果R_{k, l}中所包含的第m个右上角点的坐标为(x_{i, j}+1/2, y_{i, j}+1/2) ^T, 那么第m+1个右上角点只有两种可能, 即(x_{i, j}+3/2, y_{i, j}+1/2) ^T或者(x_{i, j}+3/2, y_{i, j}+3/2) ^T.如果第m+1个右上角点的坐标为(x_{i, j}+3/2, y_{i, j}+1/2) ^T, 我们将其称为前行点(图 4中的B₅到B₆).如果第m+1个右上角点的坐标为(x_{i, j}+3/2, y_{i, j}+3/2) ^T, 我们将其称为上跳点(图 4中的A₃到A₄).

确定第m+1个右上角点到底是前行点还是上跳点的关键在于第m个右上角点到直线区域下边界的竖直距离$d_{k,l}^{(m)}$(即图 4中的线段A₂A₃或B₂B₅).令r_{k, l}=u_{k, l}/v_{k, l}, 可以得到如下的判断准则(如图 4所示):

1) 第m+1个右上角点为

$\left\{ \begin{array}{*{35}{l}} {{({{x}_{k,l}}+3/2,{{y}_{k,l}}+1/2)}^{T}}, & 如果d_{k,l}^{(m)}-{{r}_{k.l}}>0 \\ {{({{x}_{k,l}}+3/2,{{y}_{k,l}}+3/2)}^{T}}, & 如果d_{k,l}^{(m)}-{{r}_{k.l}}\le 0 \\ \end{array} \right.$

(20)

2) 第m+1个右上角点到直线区域下边界的竖直距离为

$\left\{ \begin{array}{*{35}{l}} d_{k,l}^{(m+1)}=d_{k,l}^{(m)}-{{r}_{k.l}},\text{ 如果}d_{k,l}^{(m)}-{{r}_{k.l}}>0 \\ d_{k,l}^{(m+1)}=d_{k,l}^{(m)}-{{r}_{k.l}}+1,\text{ 如果}d_{k,l}^{(m)}-{{r}_{k.l}}\le 0 \\ \end{array} \right.$

(21)

事实上, 公式(20) 和公式(21) 给出了一个交叉迭代过程, 迭代的初始值为:1) 第1个右上角点(x_{i, j}+1/2, y_{i, j}+1/2) ^T以及2) $d_{k,l}^{(0)}=1-1/2{{r}_{k,l}}$.于是就求出了R_{k, l}中所包含的所有的右上角点(即图 4中的“小圆圈”).

对于(u_{k, l}, v_{k, l})^T的方向在其他角度范围内的情况, 只需进行微小的调整即可.例如, 如果(u_{k, l}, v_{k, l})^T的方向在45°~90°之间, 仍然考虑像素区域的右上角点, 只是此时应该不断增加y值的大小, 然后判断x的值是否需要加1.该方法最早由Bresenham提出^[18], 是计算机图形学中经典的离散直线(或曲线)生成算法^[19].

4.2 计算相交面积

在找到所有的右上角点以后(仍然以(u_{k, l}, v_{k, l})^T的方向在0°和45°之间的情况为例), 需要进一步计算权重系数, 即, 相应的像素区域和A_{k, l}所“扫过”的区域直线段区域之间的相交面积.仔细观察图 3不难发现, 对于右上角点, 存在如下两种不同的情况(如图 4所示)。

1) 如果该右上角点是上跳点(例如图 4中的A₄), 那么在该右上角点左侧, 有3个像素点属于R_{k, l};

2) 如果该右上角点是前行点(例如图 4中的B₆), 那么在该右上角点左侧, 只有两个像素点属于R_{k, l}.

在由公式(20) 和公式(21) 所构成的交叉迭代过程中, 记录下每一个新得到的右上角点是上跳点还是前行点, 于是就得到了R_{k, l}中的所有像素点, 然后计算R_{k, l}所包含的每一个像素点(i, j)所对应的权重系数w_{i, j, k, l}, 即

${{w}_{i,j,k,l}}={{{S}_{i,j}}}/{\left( \sum\limits_{(i,j)\in {{R}_{k,l}}}{{{S}_{i,j}}} \right),}\;\text{ 对于}(i,j)\in {{R}_{k,l}}$

(22)

其中, S_{i, j}表示像素点区域A_{k, l}与直线段区域(即, 图 3中两条虚线之间的区域)之间的相交面积.一个重要的观察结果是, 图 4中的四边形A₁A₂A₈A₆和B₁B₂B₃B₄的面积都是1个像素区域.

对于上跳点的情况, 三角形A₂A₃A₇的面积为

$S({{A}_{2}}{{A}_{3}}{{A}_{7}})=\frac{1}{2}{{A}_{2}}{{A}_{3}}\times {{A}_{3}}{{A}_{7}}=\frac{1}{2}{{{(d_{k,l}^{(m)})}^{2}}}/{{{r}_{k,l}}}\;$

(23)

而三角形A₅A₄A₆的面积为

$S({{A}_{5}}{{A}_{4}}{{A}_{6}})=\frac{1}{2}{{A}_{5}}{{A}_{4}}\times {{A}_{4}}{{A}_{6}}=\frac{1}{2}{{{(1-d_{k,l}^{(m+1)})}^{2}}}/{{{r}_{k,l}}}\;$

(24)

进一步求得:

$S({{A}_{1}}{{A}_{3}}{{A}_{7}}{{A}_{8}}{{A}_{4}}{{A}_{5}})=1-S({{A}_{2}}{{A}_{3}}{{A}_{7}})-S({{A}_{5}}{{A}_{4}}{{A}_{6}})$

(25)

对于前行点的情况, 计算过程更加简单.首先计算:

$S({{B}_{5}}{{B}_{2}}{{B}_{3}}{{B}_{6}})=\frac{1}{2}({{B}_{2}}{{B}_{5}}+{{B}_{3}}{{B}_{6}})\times 1=\frac{1}{2}(d_{k,l}^{(m)}+d_{k,l}^{(m+1)})$

(26)

然后, 计算梯形B₁B₅B₆B₄的面积:

$S({{B}_{1}}{{B}_{5}}{{B}_{6}}{{B}_{4}})=1-S({{B}_{5}}{{B}_{2}}{{B}_{3}}{{B}_{6}})$

(27)

5 数值求解

为了进一步验证上述理论分析结果与实际的吻合程度, 本文以真实场景为对象设计了实验.高速铁路检测列车前端安装的摄像机获取的图像是典型的前向运动视频.列车的高速前向运动, 使得图像中含有很强的运动模糊.如果在去模糊过程中采用前向运动模糊核所获得的去模糊效果优于卷积模糊核的效果, 无疑可以进一步表明本文的研究思路以及所提出前向运动模糊核是正确的.为此, 本文建立依据卷积模糊核与前向运动模糊核的图像复原模型, 并对两者进行比较.

5.1 非线性热扩散方程

我们通过求解一个泛函优化模型^[20]:

$\underset{{{I}_{0}}}{\mathop{\min }}\,\left\{ \frac{1}{2}||{{I}_{b}}-\iint{w(x,y,{{x}_{1}},{{y}_{1}}){{I}_{0}}({{x}_{1}},{{y}_{1}})\text{d}{{x}_{1}}\text{d}{{y}_{1}}}||_{2}^{2}+\lambda |||\nabla {{I}_{0}}||{{|}_{1}} \right\}$

(28)

来实现图像复原, 其中, I₀(x, y)的全变分的表达式为

$|||\nabla {{I}_{0}}||{{|}_{1}}=\iint{|\nabla {{I}_{0}}|\text{d}x\text{d}y}=\iint{\sqrt{{{(\partial {{I}_{0}}/\partial x)}^{2}}+{{(\partial {{I}_{0}}/\partial y)}^{2}}}\text{d}{{x}_{1}}\text{d}{{y}_{1}}}$

(29)

通过变分法^[20], 可以得到求解公式(28) 的演化方程:

$\frac{\partial }{\partial t}\nabla {{I}_{0}}=div(\nabla {{I}_{0}}/|\nabla {{I}_{0}}|)+\lambda \cdot s(x,y)$

(30)

即, 公式(28) 所对应的变分梯度下降流^[20].其中, $\nabla $表示梯度算子, div表示散度算子^[16], 而s(x, y)的表达式为

$s(x,y)=\iint{w(x,y,{{x}_{1}},{{y}_{1}})({{I}_{b}}({{x}_{1}},{{y}_{1}})-\iint{w(x,y,{{x}_{1}},{{y}_{1}}){{I}_{0}}(x,y)\text{d}x\text{d}y})\text{d}{{x}_{1}}\text{d}{{y}_{1}}}$

(31)

公式(30) 是一个各向异性的非线性热扩散方程^[20].

5.2 迭代算法

求解公式(30) 的离散数值方法(例如有限元和有限差分法^[16])就是相应的广义运动模糊复原算法.本文中, 我们采用有限差分的显式迭代格式来进行计算.令${{I}^{(n)}}=\{I_{i,j}^{(n)}\}$表示第n步迭代中复原出的图像, ${{I}^{(b)}}=\{I_{i,j}^{(b)}\}$表示模糊图像, D_x={D_{i, j;x}}和D_y={D_{i, j;y}}表示梯度矩阵, 即

$D_{i,j;x}^{(n)}=I_{i,j}^{(n)}-I_{i,j-1}^{(n)};D_{i,j;y}^{(n)}=I_{i,j}^{(n)}-I_{i-1,j}^{(n)}$

(32)

首先计算:

$p_{i,j;x}^{(n)}=\frac{D_{i,j;x}^{(n)}}{\sqrt{{{(D_{i,j;x}^{(n)})}^{2}}+{{(D_{i,j;y}^{(n)})}^{2}}}};p_{i,j;y}^{(n)}=\frac{D_{i,j;y}^{(n)}}{\sqrt{{{(D_{i,j;x}^{(n)})}^{2}}+{{(D_{i,j;y}^{(n)})}^{2}}}}$

(33)

然后, 进一步计算求得:

$q_{i,j}^{(n)}=p_{i,j+1;x}^{(n)}-p_{i,j;x}^{(n)}+p_{i+1,j;y}^{(n)}-p_{i,j;x}^{(n)}$

(34)

$s_{i,j}^{(n)}=\sum\limits_{k}{\sum\limits_{l}{\left( {{w}_{i,j,k,l}}\cdot \left( I_{k,l}^{\left( b \right)}-\left( \sum\limits_{{{i}'}}{\sum\limits_{{{j}'}}{{{w}_{{i}',{j}',k,l}}\cdot I_{{i}',{j}'}^{(n)}}} \right) \right) \right)}}$

(35)

最终, 得到了公式(30) 所对应的显式迭代格式:

$I_{i,j}^{(n+1)}=I_{i,j}^{(n)}+{{\lambda }_{1}}\cdot q_{i,j}^{(n)}+{{\lambda }_{2}}\cdot s_{i,j}^{(n)}$

(36)

迭代的初始值选为模糊图像(即I⁽⁰⁾ =I^(b)).

6 实验结果

现在我们来考虑实际情况——高速检测列车所拍摄的图像进行去模糊处理.车头前置摄像机获取图像为前向运动图像, 因此, 使用传统的反卷积模型无法取得好的图像复原效果.图 5给出了模糊图像的复原结果.

图 6是根据图 5(a)~图 5(d)计算出的模糊核结果.图 6(a)给出了一部分前向运动模糊核w_{i, j, k, l}, 每一个“小方块”的大小为39x39个像素点, 对应于一个固定的(k, l).图6(b)中给出了反卷积模型^[4]所估计出的模糊核.

从图6(b)可以看出, 反卷积模型“认为”模糊是由图像的“上下晃动”产生的.基于这个错误的图像模糊假设, 复原方法将无法有效地去除图像中的运动模糊, 甚至引入了严重的伪影.相比之下, 我们的方法取得了较为理想的效果, 有效地去除了图像中的模糊, 并且在图像的中心区域, 即模糊程度很小的区域, 并没有产生“伪影”(反卷积模型的复原结果中存在明显的伪影(如图 7所示)).

为了便于进一步对比观察, 图 7给出了图 5中的模糊图像、反卷积模型的复原结果, 即, 图 5(a)~图 5(d)中方框圈出来的小图像块.

7 总结

图像复原模型是一个逆问题, 对于求解逆问题来说, 对正问题的正确理解是极其重要的.经典的图像复原方法(即反卷积模型)将图像的模糊过程模拟为一个卷积过程.对于前向运动模糊, 这个近似假设是不成立的.本文提出了一种广义运动模糊模型, 在此基础上, 给出了前向运动模糊模型和前向运动模糊核, 其目的就是为恢复前向运动模糊找到理论依据.在这里, 首先给出了基于光流场的广义运动模糊模型的理论分析;接下来给出了前向运动模糊模型的快速的运动估计方法、平面场景朝向的快速估计方法、前向运动模糊核的数值生成方法;最后, 通过对实际的前向运动模糊图像的去模糊实验, 验证了本文结论的正确性.

本文研究的前向运动模糊只是广义运动模糊的一种特殊情况, 即便是对于相机相对于固定场景运动的特殊情况, 相机的运动也可能同时包含平动和转动^{[14, 15]}.此时, 估计光流和生成广义运动模糊核将变得更加困难.如何结合更加复杂的相机运动假设以及更加复杂的正则化方法^[21]设计出针对其他类型的广义运动模糊的图像复原算法, 将是我们进一步研究和探索的问题.