发布时间: 2017-01-25
摘要点击次数:
全文下载次数:
DOI: 10.11834/jig.20170108
2017 | Volumn 22 | Number 1

图像分析和识别

大位移变分光流计算的快速算法

刘博文, 魏伟波, 潘振宽, 王守润

青岛大学计算机科学技术学院, 青岛 266071

收稿日期: 2016-07-14; 修回日期: 2016-10-07

基金项目: 国家科技支撑计划项目（2014BAG03B05）

第一作者简介: 刘博文(1991-), 男, 现为青岛大学计算机科学与技术专业硕士研究生, 研究方向为变分图像处理。E-mail:345469300@qq.com

中图法分类号: TP391

文献标识码: A

文章编号: 1006-8961(2017)01-0066-09

摘要

目的多尺度方法的提出解决了传统HS（Horn Schunck）算法不能计算大位移光流的问题，但同时也增加了迭代运算的步数。为加快迭代收敛速度，研究大位移变分光流计算的快速算法，并分析其性能。方法将用于加快变分图像处理迭代运算的Split Bregman方法、对偶方法和交替方向乘子法应用到大位移光流计算中。结果分别进行了精度、迭代步数、运行时间的对比实验。引入3种快速方法的模型均能够在保证精度的同时，在较少时间内计算出图像序列的光流场，所需时间为传统方法的11%~42%。结论将3种快速方法应用到大位移变分光流计算中，对于不同图像序列均可以较大地提高计算效率。

关键词

光流计算; 大位移光流; 多尺度方法; Split Bregman方法; 对偶方法; 交替方向乘子法

Fast algorithms for large displacement variation optical flow computation

Liu Bowen, Wei Weibo, Pan Zhenkuan, Wang Shourun

College of Computer Science and Technology, Qingdao University, Qingdao 266071, China

Supported by: Supported by:National Key Technology Research and Development Program of the Ministry of Science and Technology of China (2014BAG03B05)

Abstract

Objective The Horn-Schunck (HS) algorithm is one of the most popular optical flow estimation methods. Many scholars have proposed improved HS algorithms to improve accuracy. However, the efficiency of the HS algorithm remains an important problem because the HS algorithm requires much iterative computation. The HS algorithm is based on a differential method, and it only can compute small displacement optical flow. A multi-scale method has been proposed to solve the problem that a differential method cannot compute large displacement optical flow, but the efficiency of this method is slower than before. Fast methods are studied in this research to enhance efficiency. Method In the variation image restoration domain, fast methods for accelerating iteration have yielded good results, and some of the fast methods have been applied to small-displacement optical flow computation domain. In this study, Split Bregman method, dual method, and alternating direction method of multipliers are applied to large displacement optical flow computation for accelerating iteration. Result The accuracy, iteration, and time of different methods are compared quantitatively and qualitatively. The three fast methods all can obtain results with accuracy that is the same as that of the traditional method in lesser time. The time required for fast algorithms is 11%~42% of the time required for the traditional method. Conclusion Computational efficiency can be improved greatly by applying these three fast methods to large displacement variation optical flow computing for different image sequences.

Key words

optical flow computation; large displacement optical flow; multi-scale method; Split Bregman method; dual method; alternating direction method of multipliers

0 引言

光流计算在计算机视觉具有非常重要的作用。所谓光流指图像序列中每个像素点的运动。由于图像序列中目标的运动场难以计算，所以通常用像素点的光流场代替运动场。一般来说，光流计算的方法有以下几种：基于微分的算法(其经典算法为HS (Horn Schunck)算法^[1]和LK (Lucas Kanade)算法^[2-3])、基于区域匹配的方法^[4-5]、基于相位的方法^[6-7]、基于能量的方法^[8]。文献[9]比较了几种通用光流算法的性能，发现HS算法综合性能较好。另外由于变分法和偏微分方程(PDEs)理论的不断完善，HS算法成为研究较为广泛的方法。因为HS算法基于变分法和偏微分方程理论，也可将其称为变分光流计算方法。

大多数学者关注如何得到更准确的光流场。基于变分法的图像处理技术将原问题转化为求解能量泛函的极值问题，能量泛函通常由数据项(data term)和光滑项(smooth term)组成，变分光流计算也不例外。为了提高光流计算的精确度，许多学者做出了努力。在图像序列中，光照变化、目标的运动类型等因素都会影响到光流计算的精度，变分光流模型中直接受到这些因素影响的是数据项，因此部分学者对其数据项进行改进，提出了几种不变假设模型，例如梯度不变假设、拉普拉斯不变假设等等。模型中光滑项决定着光流场的扩散性能，为了让光流场扩散得更准确，部分学者对光滑项进行改进，例如图像(或光流)驱动、图像光流联合驱动^[10]、各向同性(或各向异性)扩散，利用不同的光滑项对光流场进行光滑。文献[11-12]对各种数据项和光滑项进行了整理和介绍。另一部分学者提出了一些其他提高精度的方法，例如文献[13]提出的利用图像不同位置处的特征距离对光流矢量进行匹配调整的方法、文献[14]提出的在模型中添加惩罚项的方法。变分光流模型通常基于HS模型，除此之外，文献[15]将流体与光流结合提出了用于计算流体光流的模型，文献[16]在流体光流模型中添加惩罚项使结果更为精确。文献[17]对近几年光流技术进行了总结。文献[18]将光流技术与背景建模技术结合提出一种新的运动检测方法。由于变分光流计算模型是基于微分的方法，在计算过程中采用了泰勒展开式，而泰勒展开式的前提条件是变化量趋于0，因此微分光流算法在估计大位移光流时会出现较大误差，许多文献[11, 19-20]提到了多尺度方法(也可称为多分辨率分层细化策略、金字塔方法)，多尺度方法较好地解决了这一问题，但该方法需要对图像分层，然后对每一层进行迭代计算，极大地增加了计算复杂度，因此大位移光流计算的加速问题变得尤为重要。

用于解决图像处理加速迭代问题的快速方法在图像恢复、图像分解等领域已经发挥了重要作用，例如Split Bregman算法^[21]、对偶(dual)算法^[22-23]、交替方向乘子法(ADMM)算法^[24-25]等。其中一部分算法已经应用到小位移光流计算领域^{[23, 26-27]}，并在提升算法效率方面取得了突破。其中，Bregman算法最早用于解决求解变分模型极值迭代步数较多的问题，文献[21]将其与Split思想结合提出了Split Bregman算法。对偶方法定义了一个对偶空间，用求解对偶项极值代替求解原模型极值，避免了原模型求解过程中分母可能为0的问题，同时也加快了迭代速度。增广拉格朗日(ALM)方法^[24]将多种快速方法纳入了一个体系，由于快速方法通常需要求解多个未知量，因此需要采用ADMM方法对未知量进行交替求解。

本文重点研究了几种快速算法对减少大位移光流计算迭代步骤的影响。首先介绍了变分光流计算的经典模型--HS模型和多尺度方法，然后将Split Bregman算法、Dual算法、ADMM算法应用到大位移光流计算模型中，最后利用数值实验对几种方法的性能进行比较。

1 变分光流算法介绍

1.1 HS模型

在图像序列中，可以用$f (x, y, t)$表示t时刻的图像的灰度值, 用$\boldsymbol{u}={[{u_1}, {u_2}]^{\rm{T}}}$表示光流，光流分量${u_1}, {u_2}$分别表示光流的水平方向分量和垂直方向分量。灰度一致假设可表示为

$ f(x + {u_1}, y + {u_2}, t + 1)-f(x, y, t) = 0 $

(1)

将式(1)用泰勒展开式展开得到光流约束(OFC)方程，HS模型将OFC方程作为变分模型的数据项，OFC方程为

$ {f_x}{u_1} + {f_y}{u_2} + {f_t} = 0 $

(2)

式中

$ \begin{array}{l} {f_x} = \frac{{\partial f}}{{\partial x}}, {f_y} = \frac{{\partial f}}{{\partial y}}\\ {f_t} = f(t + 1)-f(t) \end{array} $

由于泰勒展开式的前提条件是变化量趋于0，因此HS模型只能计算小位移光流。

式(2)中含有两个未知量，因此有无穷多个解。HS算法假设图像序列中光流场是连续光滑的，将式(2)与光滑项结合，得到HS模型为

$ E(\boldsymbol{u})=\frac{1}{2}\left[\begin{align} & \underbrace{\iint_{\Omega }{{{({{f}_{x}}{{u}_{1}}+{{f}_{y}}{{u}_{2}}+{{f}_{t}})}^{2}}\text{d}x\text{d}y}}_{\text{data}\ \text{term}}+ \\ & \underbrace{\lambda \iint_{\Omega }{({{\left| \nabla {{u}_{1}} \right|}^{2}}+{{\left| \nabla {{u}_{2}} \right|}^{2}})\text{d}x\text{d}y}}_{\text{smooth}\ \text{term}} \\ \end{align} \right] $

(3)

式中，$\lambda $为参数，调节光滑项所占比重。

1.2 多尺度方法

由于基于微分的光流算法利用了泰勒展开式，前提条件是变化量趋于0，因此微分光流算法只能计算小位移光流。而在图像中，位移的单位为像素，位移量的大小表现为相隔的像素的个数，像素又与图像分辨率直接相关，因此，可对原图像进行缩小，将原来的大位移运动变为小位移运动。将原图像缩小不同尺寸，得到许多个小分辨率的图像，原来整个的大位移变为一个个的小位移，然后再将小位移合并得到大位移光流。这种方法称为多尺度方法。多尺度方法具有通用性，为了表示方便，本文用HS模型为例。

将原图像缩小不同分辨率，从低分辨率到高分辨率依次分为第1层到第$k$层，将原来的大位移光流矢量$\boldsymbol{u}$分为多个小位移光流之和，在每一层分辨率下分别计算小位移光流然后相加。第$k$层的光流可以分为前$k$-1层光流之和与第$k$层光流两部分，前$k$-1层光流是大位移光流，第$k$层光流是小位移光流，用公式表示为

$ {{\boldsymbol{u}}^{k}}={{\boldsymbol{u}}^{k-1}}+\text{d}{{\boldsymbol{u}}^{k-1}} $

(4)

为了表述方便，定义

$ \left\{ \begin{array}{l} {f_x} = {\partial _x}f(\boldsymbol{x} + \boldsymbol{u}, t + 1)\\ {f_y} = {\partial _y}f(\boldsymbol{x} + \boldsymbol{u}, t + 1)\\ {f_z} = f(\boldsymbol{x} + \boldsymbol{u}, t + 1)-f(\boldsymbol{x}, t) \end{array} \right. $

(5)

式中，$f (\boldsymbol{x} + \boldsymbol{u}, t + 1)$利用$f (\boldsymbol{x}, t + 1)$和$\boldsymbol{u}$插值得到，插值算法可采用双线性插值。大位移光流模型的灰度一致方程为

$ \begin{array}{l} f(x + {u_1} + {\rm{d}}{u_1}, y + {u_2} + {\rm{d}}{u_2}, t + 1)-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;f(x, y, t) = 0 \end{array} $

(6)

大位移光流的OFC方程为

$ {f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_z} = 0 $

(7)

根据大位移光流的OFC方程可得到引入多尺度方法的HS模型，即

$ \begin{array}{l} E(\boldsymbol{u}) = \\ \frac{1}{2}\left[ \begin{array}{l} \int {\int_\Omega {{{({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_z})}^2}{\rm{d}}x{\rm{d}}y} } \\ \lambda \int {\int_\Omega {({{\left| {\nabla ({u_1} + {\rm{d}}{u_1})} \right|}^2} + {{\left| {\nabla ({u_2} + {\rm{d}}{u_2})} \right|}^2}){\rm{d}}x{\rm{d}}y} } \end{array} \right] \end{array} $

(8)

若要将多尺度方法引入其他变分光流模型，只需将式(4)代入原模型。

对式(8)求欧拉-拉格朗日方程，得到

$ \left\{ \begin{array}{l} ({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_z}){f_x}-\\ \;\;\;\lambda \Delta ({u_1} + {\rm{d}}{u_1}) = 0\\ ({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_z}){f_y}-\\ \;\;\;\lambda \Delta ({u_2} + {\rm{d}}{u_2}) = 0 \end{array} \right. $

(9)

在每一层分辨率下首先采用式(9)计算得到${\rm{d}}\boldsymbol{u}$，然后利用式(4)得到$\boldsymbol{u}$，第1层的初始光流设为0，即${\boldsymbol{u}^0}=\boldsymbol{0}$。采用改进的变分光流模型计算时，由于计算成本较大，为了节约时间可在低分辨率下使用各向同性扩散项计算，高分辨率下采用改进的扩散项计算。

2 大位移变分光流计算的快速算法

选用光滑项为TV正则项的引入多尺度方法的HS模型(M-HS-TV model)作为举例，TV项为全变分去噪模型^[28]的光滑项，具有边缘保持的特点，M-HS-TV模型为

$ \begin{align} & E(\text{d}\boldsymbol{u})=\frac{1}{2}{{\iint_{\Omega }{({{f}_{x}}\text{d}{{u}_{1}}+{{f}_{y}}\text{d}{{u}_{2}}+{{f}_{z}})}}^{2}}\text{d}x\text{d}y+ \\ & \lambda \iint_{\Omega }{(\left| \nabla ({{u}_{1}}+\text{d}{{u}_{1}}) \right|+\left| \nabla ({{u}_{2}}+\text{d}{{u}_{2}}) \right|)\text{d}x\text{d}y} \\ \end{align} $

(10)

2.1 Split Bregman方法

Split Bregman算法是将分裂思想与Bregman算法结合的方法，首先需要介绍Bregman算法。

设求解问题为

$ \mathop {\min }\limits_u \{ J(u) + \theta H(u)\} $

(11)

式中，$H (u)$为约束条件，$\theta $为参数。

定义$u$到$v$点的Bregman距离为

$ D_J^P(u, v) = J(u)-J(v)-< p, u-v > $

(12)

式中，$p$表示函数$J$中$v$点处的子梯度(subgradient)。函数$J$的子梯度的集合用$\partial J$表示，则$p \in \partial J$。利用Bregman距离，可将原问题式(11)转化为问题

$ \mathop {\arg \;\min }\limits_u \{ D_J^P(u, {u^k}) + \theta H(u)\} $

(13)

得到

$ \left\{ \begin{array}{l} {u^{k + 1}} = \mathop {\arg \;\min }\limits_u \{ J(u)-< p, u-{u^k} > + \theta H(u)\} \\ {p^{k + 1}} = {p^k}-\theta \nabla H \end{array} \right. $

(14)

变分图像处理模型中难以求解的是含有$\left| {\nabla u} \right|$的项，分裂(Split)方法是引入辅助变量$\boldsymbol{w}$，用$\boldsymbol{w}$代替${\nabla u}$，令

$ \boldsymbol{w} = \nabla u $

(15)

此时约束条件变为

$ H(u) = {\left| {\boldsymbol{w}-\nabla u} \right|^2} $

(16)

根据式(14)可得到

$ {p^{k + 1}} = {p^k}-\theta ({\boldsymbol{w}^k}-\nabla {u^k}) $

(17)

根据文献[21]方法，约束条件为线性时，可用公式

$ \left\{ \begin{array}{l} H(u) = {\left| {\boldsymbol{w}-\nabla u-{b^k}} \right|^2}\\ {b^{k + 1}} = {b^k}-({\boldsymbol{w}^k} - \nabla {u^k}) \end{array} \right. $

(18)

代替式(16)(17)。

对于模型式(10)，引入Split Bregman方法后得到

$ \begin{align} & E(\text{d}\boldsymbol{u})=\frac{1}{2}{{\iint_{\Omega }{({{f}_{x}}\text{d}{{u}_{1}}+{{f}_{y}}\text{d}{{u}_{2}}+{{f}_{t}})}}^{2}}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \sum\limits_{i=1}^{2}{[\lambda \iint_{\Omega }{\left| {{\boldsymbol{w}}_{i}} \right|\text{d}x\text{d}y}} \\ & \ \ \ \ \ \frac{\theta }{2}\iint_{\Omega }{{{({{\boldsymbol{w}}_{i}}-\nabla ({{u}_{i}}+\text{d}{{u}_{i}})-\boldsymbol{b}_{i}^{k})}^{2}}\text{d}x\text{d}y}] \\ \end{align} $

(19)

对式(19)求欧拉-拉格朗日方程得

$ \left\{ \begin{array}{l} ({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_z}){f_x} + \theta \nabla \cdot \\ \;\;\;({\boldsymbol{w}_1}-\nabla ({u_1}{\rm{ + d}}{u_1})-{\boldsymbol{b}_1}) = 0\\ ({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_z}){f_y} + \theta \nabla \cdot \\ \;\;\;({\boldsymbol{w}_2}-\nabla ({u_2}{\rm{ + d}}{u_2}) - {\boldsymbol{b}_2}) = 0\\ {\boldsymbol{w}_1} = \max (\left| {\nabla ({u_1} + {\rm{d}}{u_1}) + {\boldsymbol{b}_1}} \right| - \frac{\lambda }{\theta }, 0)\\ \;\;\;\;\frac{{\nabla ({u_1} + {\rm{d}}{u_1}) + {\boldsymbol{b}_1}}}{{\left| {\nabla ({u_1} + {\rm{d}}{u_1}) + {\boldsymbol{b}_1}} \right|}}, 0\frac{\boldsymbol{0}}{{\left| \boldsymbol{0} \right|}} = \boldsymbol{0}\\ {\boldsymbol{w}_2} = \max (\left| {\nabla ({u_2} + {\rm{d}}{u_2}) + {\boldsymbol{b}_2}} \right| - \frac{\lambda }{\theta }, 0)\\ \;\;\;\;\frac{{\nabla ({u_2} + {\rm{d}}{u_2}) + {\boldsymbol{b}_2}}}{{\left| {\nabla ({u_2} + {\rm{d}}{u_2}) + {\boldsymbol{b}_2}} \right|}}, 0\frac{\boldsymbol{0}}{{\left| \boldsymbol{0} \right|}} = \boldsymbol{0} \end{array} \right. $

(20)

式中，$\boldsymbol{b}$的计算公式为

$ \left\{ \begin{array}{l} \boldsymbol{b}_1^{k + 1} = \boldsymbol{b}_1^k + \nabla (u_1^k + {\rm{d}}u_1^k)-\boldsymbol{w}_1^k\\ \boldsymbol{b}_2^{k + 1} = \boldsymbol{b}_2^k + \nabla (u_2^k + {\rm{d}}u_2^k)-\boldsymbol{w}_2^k \end{array} \right. $

(21)

式(20)中计算$\boldsymbol{w}$的公式称为广义软阈值公式，避免了分母为0的情况，可表示为

$ \left\{ \begin{array}{l} {\boldsymbol{w}_1} = {\rm{shrink}}(\nabla ({u_1} + {\rm{d}}{u_1}) + {\boldsymbol{b}_1}, \frac{\lambda }{\theta })\\ {\boldsymbol{w}_2} = {\rm{shrink}}(\nabla ({u_2} + {\rm{d}}{u_2}) + {\boldsymbol{b}_2}, \frac{\lambda }{\theta }) \end{array} \right. $

(22)

式中

$ {\rm{shrink}}(x, \boldsymbol{\lambda }) = \left\{ \begin{array}{l} x-\boldsymbol{\lambda} \frac{\boldsymbol{x}}{{\left| \boldsymbol{x} \right|}}\;\;\;\;\;\left| \boldsymbol{x} \right| > \lambda \\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left| \boldsymbol{x} \right| \le \lambda \end{array} \right. $

2.2 Dual方法

对于变分模型中计算复杂的光滑项，Split Bregman方法采用分裂思想，用辅助变量$\boldsymbol{w}$代替计算$\nabla u$。而对偶方法定义了对偶空间

$ \sup \text{ }\!\!\{\!\!\text{ }\iint_{\Omega }{u(\nabla \cdot \boldsymbol{p})}\text{d}x\text{d}y:\left| \boldsymbol{p} \right|\le \text{1 }\!\!\}\!\!\text{ } $

(23)

使用对偶项

$ E(\boldsymbol{p})=\iint_{\Omega }{u(\nabla \cdot \boldsymbol{p})}\text{d}x\text{d}y $

(24)

代替光滑项

$ {{E}_{\text{smooth}}}(u)=\iint_{\Omega }{\left| \nabla u \right|}\text{d}x\text{d}y $

(25)

对于模型式(10)，需要引入辅助变量$\text{d}\boldsymbol{v}$，利用$\text{d}v$计算光流矢量$\text{d}u$，得到

$ \begin{align} & E(\text{d}\boldsymbol{u}, \text{d}\boldsymbol{v}, \boldsymbol{p})=\frac{1}{2}\iint_{\Omega }{{{({{f}_{x}}\text{d}{{v}_{1}}+{{f}_{y}}\text{d}{{v}_{2}}+{{f}_{t}})}^{2}}}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \sum\limits_{i=1}^{2}{[}\lambda \iint_{\Omega }{({{u}_{i}}+\text{d}{{u}_{i}})\nabla \cdot {{\boldsymbol{p}}_{i}}}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{\theta }{2}\iint_{\Omega }{{{(\text{d}{{u}_{i}}+\text{d}{{v}_{i}})}^{2}}\text{d}x\text{d}y}] \\ \end{align} $

(26)

固定$\text{d}\boldsymbol{v}$，对$\boldsymbol{p}$和$\text{d}\boldsymbol{u}$进行求解，得到问题

$ \begin{align} & \underset{u}{\mathop{\min }}\, \underset{\boldsymbol{p}}{\mathop{\operatorname{m}\text{ax}}}\, \{E(\text{d}\boldsymbol{v}, \boldsymbol{p})=\sum\limits_{i=1}^{2}{[}\lambda \iint_{\Omega }{({{u}_{i}}+\text{d}{{u}_{i}})\nabla \cdot {{\boldsymbol{p}}_{i}}}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \ \frac{\theta }{2}\iint_{\Omega }{{{(\text{d}{{u}_{i}}+\text{d}{{v}_{i}})}^{2}}\text{d}x\text{d}y}]:\left| {{\boldsymbol{p}}_{i}} \right|\le 1\}, i=1, 2 \\ \end{align} $

(27)

式(27)对$\text{d}\boldsymbol{u}$的变分导(函)数为

$ \text{d}{{u}_{i}}=\text{d}{{v}_{i}}-\frac{\lambda }{\theta }\nabla \cdot {{\boldsymbol{p}}_{i}}, i=1, 2 $

(28)

计算得到

(29)

2.2.1 原始解法

问题式(29)是一个约束问题，可将其转为无约束问题

$ \begin{align} & \underset{\boldsymbol{p}}{\mathop{\min }}\, \{E(\boldsymbol{p})=\sum\limits_{i=1}^{2}{[}\iint_{\Omega }{{{(\nabla \cdot {{\boldsymbol{p}}_{i}}-\frac{\theta ({{u}_{i}}+\text{d}{{v}_{i}})}{\lambda })}^{2}}}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\alpha }_{i}}\iint_{\Omega }{({{\left| {{\boldsymbol{p}}_{i}} \right|}^{2}}-1)\text{d}x\text{d}y}]:{{\alpha }_{i}}>0, \\ & \ \ \ \ \ \ \ \ \ \ \ \left| {{\boldsymbol{p}}_{i}} \right|=1H\alpha =0, \left| {{\boldsymbol{p}}_{i}} \right|<1\}, i=1, 2 \\ \end{align} $

(30)

对式(30)求欧拉-拉格朗日方程得到

$ \begin{array}{l} \nabla (\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda })-{\alpha _i}{\boldsymbol{p}_i} = \boldsymbol{0}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, 2 \end{array} $

(31)

当${\alpha _i} > 0, \left| {{\boldsymbol{p}_i}} \right|=1, i=1, 2$时，根据式(31)得到

$ \begin{array}{l} {\alpha _i} = \pm \nabla (\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda }) = \\ \left| {\nabla \left( {\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda }} \right)} \right|, i = 1, 2 \end{array} $

(32)

当${\alpha _i}=0, \left| {{\boldsymbol{p}_i}} \right| < 1, i=1, 2$时，根据式(31)得到

$ \nabla \left( {\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda }} \right) = \boldsymbol{0}, i = 1, 2 $

(33)

可令

$ \begin{array}{l} {\alpha _i} = \left| {\nabla (\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda })} \right|\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, 2 \end{array} $

(34)

将式(34)代入式(31)得到

$ \begin{array}{l} \frac{{\partial {\boldsymbol{p}_i}}}{{\partial t}} = \nabla \left( {\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda }} \right)-\\ \;\;\;\left| {\nabla \left( {\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({u_i} + {\rm{d}}{v_i})}}{\lambda }} \right)} \right|{\boldsymbol{p}_i}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, 2 \end{array} $

(35)

计算得到

$ \begin{array}{l} \boldsymbol{p}_i^{k + 1} = \frac{{\boldsymbol{p}_i^k + t\nabla (\nabla \cdot \boldsymbol{p}_i^k-\frac{{\theta (u_i^k + {\rm{d}}v_i^k)}}{\lambda })}}{{1 + t\left| {\nabla (\nabla \cdot \boldsymbol{p}_i^k-\frac{{\theta (u_i^k + {\rm{d}}v_i^k)}}{\lambda })} \right|}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, 2 \end{array} $

(36)

式中，$t$为步长，通常$t$≤0.125。

式(26)对${{\rm{d}}\boldsymbol{v}}$的变分导数为

$ \left\{ \begin{array}{l} ({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_t}){f_x}-\\ \;\;\;\theta ({\rm{d}}{u_1}-{\rm{d}}{v_1}) = 0\\ ({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_t}){f_y}-\\ \;\;\;\theta ({\rm{d}}{u_2} - {\rm{d}}{v_2}) = 0 \end{array} \right. $

(37)

式(36)和式(37)是求解模型式(26)的公式。

2.2.2 梯度投影法

梯度投影(gradient projection)法^[29]是求解约束问题式(29)的另一种解法。方法是直接对式(29)求欧拉-拉格朗日方程得到

$ \frac{{\partial {\boldsymbol{p}_i}}}{{\partial t}} = \nabla (\nabla \cdot {\boldsymbol{p}_i}-\frac{{\theta ({\boldsymbol{u}_i} + {\rm{d}}{\boldsymbol{v}_i})}}{\lambda }), i = 1, 2 $

(38)

再添加约束条件

$ {\boldsymbol{p}_i} = \frac{{{\boldsymbol{p}_i}}}{{\max (\left| {{\boldsymbol{p}_i}} \right|, 1)}}, i = 1, 2 $

(39)

使辅助变量$\boldsymbol{p}$满足$\left| {{\boldsymbol{p}_i}} \right| \le 1$。

2.3 ADMM方法

ADMM (交替方向乘子法)与分裂方法类似，区别在于添加了新的约束项

$ \mathit{\boldsymbol{\beta }} \cdot (\mathit{\boldsymbol{w}}-\nabla \mathit{\boldsymbol{u}}-\nabla {\rm{d}}\mathit{\boldsymbol{u}}) = 0 $

(40)

引入ADMM方法的大位移光流模型为

$ \begin{align} & E(\text{d}\boldsymbol{u})=\frac{1}{2}\iint_{\Omega }{{{({{f}_{x}}\text{d}{{u}_{1}}+{{f}_{y}}\text{d}{{u}_{2}}+{{f}_{t}})}^{2}}\text{d}x\text{d}y} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \sum\limits_{i=1}^{2}{[}\lambda \iint_{\Omega }{\left| {{\boldsymbol{w}}_{i}} \right|}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \iint_{\Omega }{{{\mathit{\boldsymbol{\beta }}}_{i}}\cdot ({{\boldsymbol{w}}_{i}}-\nabla ({{u}_{i}}+\text{d}{{u}_{i}}))}\text{d}x\text{d}y+ \\ & \ \ \ \ \ \ \frac{\theta }{2}\iint_{\Omega }{{{({{\boldsymbol{w}}_{i}}-\nabla ({{u}_{i}}+\text{d}{{u}_{i}}))}^{2}}\text{d}x\text{d}y}] \\ \end{align} $

(41)

式中，$\mathit{\boldsymbol{\beta }}$的迭代公式为

$ \begin{array}{l} \mathit{\boldsymbol{\beta }}_i^{k + 1} = \mathit{\boldsymbol{\beta }}_i^k + \theta (\boldsymbol{w}_i^k-\nabla u_i^k-\nabla {\rm{d}}u_i^k)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, 2 \end{array} $

(42)

其他未知量的求解公式为

$ \left\{ {\begin{array}{*{20}{l}} {({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_t}){f_x} + \nabla \cdot {\mathit{\boldsymbol{\beta }}_1} + }\\ {\;\;\;\;\;\theta \nabla \cdot ({\mathit{\boldsymbol{w}}_1}-\nabla ({u_1} + {\rm{d}}{u_1})) = 0}\\ {({f_x}{\rm{d}}{u_1} + {f_y}{\rm{d}}{u_2} + {f_t}){f_y} + \nabla \cdot {\mathit{\boldsymbol{\beta }}_2} + }\\ {\;\;\;\;\;\theta \nabla \cdot ({\mathit{\boldsymbol{w}}_2}-\nabla ({u_2} + {\rm{d}}{u_2})) = 0}\\ \begin{array}{l} {\mathit{\boldsymbol{w}}_1} = {\rm{shrink}}(\nabla ({u_1} + {\rm{d}}{u_1})-\frac{{{\mathit{\boldsymbol{\beta }}_1}}}{\theta }, \frac{\lambda }{\theta })\\ {\mathit{\boldsymbol{w}}_2} = {\rm{shrink}}(\nabla ({u_2} + {\rm{d}}{u_2}) - \frac{{{\mathit{\boldsymbol{\beta }}_2}}}{\theta }, \frac{\lambda }{\theta }) \end{array} \end{array}} \right. $

(43)

3 数值实验

为了说明快速方法对加速迭代收敛的效果，本文分别采用模型(10)和引入3种快速方法的模型对图像序列进行计算，采用的实验图像是来自Middlebury光流数据库^[30]的Dimetrodon图像序列和Grove2图像序列。实验环境为Intel (R) Pentium (R) 4 CPU 2.66 GHz, 1.99 GB内存，采用的编程环境为Matlab R2010b。

Dimetrodon图像为固定场景图像，Grove2图像为转动场景图像，分别如图 1所示，图 1(a)(b)均为图像序列的相邻两帧图像，图 1(c)表示对应的用彩色空间表示的标准光流场，不同颜色表示不同方向的光流，颜色亮度大小表示光流模的大小。

图 1 Dimetrodon图像序列与标准光流场

Fig. 1 Dimetrodon, Grove2 sequence and standard optical flow (color plot) ((a) frame 10; (b) frame 11; (c) standard optical flow)

按照多尺度方法共分30层，采样因子0.95。为了使对比实验更为简洁和保证条件一致，低分辨率下均采用线性各向同性模型，最高分辨率下采用模型(10)和引入快速方法的模型，当能量函数变化量小于一定阈值时停止迭代。如果快速方法在最高分层可以发挥加速作用，说明在所有分层都可以发挥加速作用。平均角误差(AAE)计算公式为

$ \begin{array}{l} {\rm{AAE = }}\frac{1}{M}\sum\limits_{I = 2}^M {\arccos ((u_i^c{u_i} + v_i^c{v_i} + {k^2})/} \\ (\sqrt {{{(u_i^c)}^2} + {{(v_i^c)}^2} + {k^2}} \times \sqrt {{{({u_i})}^2} + {{({v_i})}^2} + {k^2}} )) \end{array} $

(44)

式中，$M$表示图像总像素，${u^c}$、${v^c}$表示标准光流，$k$表示图像序列相隔帧数。

图 2表示对图 1 Dimetrodon图像序列使用不同方法得到的光流场，各向同性扩散项参数为10。其中，图 2(a)表示使用未引入快速方法的模型的结果，参数$\lambda $=10，图 2(b)表示引入Split Bregman方法的结果，参数$\lambda $=5, $\theta $=400，图 2(c)表示引入Dual方法的结果，参数$\lambda $=1, $\theta $=400，图 2(d)表示引入ADMM方法的结果，参数$\lambda $=4, $\theta $=400。由于引入不同快速方法后模型改变，为了使能量函数收敛需要重新选取最优的光滑项参数$\lambda $，因此不同方法光滑项参数不同。从视觉效果上看，Split Bregman方法与ADMM方法结果基本一致，并且与原模型结果相差不大，而Dual方法得到的光流较为“粗糙”。引入快速方法所得光流场与原模型所得光流场边缘基本一致。与标准光流相比可知，4种方法均能得到较为准确的光流场。

图 2 Dimetrodon序列使用不同方法得到的光流场

Fig. 2 Dimetrodon sequence and optical flow got by different methods ((a) M-HS-TV; (b) Split Bregman; (c) Dual; (d) ADMM)

图 3表示对Dimetrodon图像序列使用不同方法得到的能量函数值随迭代步数变化曲线。由于所用模型不同，比较能量函数最小值并无较大意义，需要比较的是不同模型到达收敛条件时的最少步数。观察图 3可知，Split Bregman方法和ADMM方法函数变化曲线基本重合，并且3种快速方法在100步以内能量函数值变化就已经趋于平稳。而M-HS-TV模型在迭代100步之后仍在明显下降。说明快速方法的加速迭代效果非常明显。

图 3 Dimetrodon图像不同方法的能量函数变化曲线

Fig. 3 Decrease of different methods' function value on increasing iteration

表 1为最高分辨率分层使用不同方法对图 1 Dimetrodon序列运算占用的CPU时间和AAE比较。从AAE结果上看，几种方法的AAE相差不大，说明快速方法不会对结果的精度产生较大影响。从时间上看，引入快速方法后，所需时间均少于原方法，Dual方法所用时间为原方法的30%左右，另外两方法所用时间为原方法的20%以内。这是因为引入快速方法后，到达收敛状态所需迭代步数大幅度下降，虽然快速方法会增加单步迭代内的计算成本，但是减少迭代步数所减少的时间大于增加的时间，因此所需总时间也减少。变分光流计算问题为线性问题，几种算法的空间复杂度均为O (${N^2}$)。

表 1 不同方法各项数据的比较
Table 1 Comparison between different methods

下载CSV

方法	AAE	迭代步数	CPU时间/s
M-HS-TV	11.828 2	400	47
Split Bregman	11.647 4	30	7
Dual	11.810 7	70	15
ADMM	11.528 2	30	5

图 4表示对Grove2图像序列使用不同方法得到的光流场，各向同性扩散项参数为10。其中，图 4(a)表示使用未引入快速方法的模型的结果，参数$\lambda $=10，图 4(b)表示引入Split Bregman方法的结果，参数$\lambda $=5, $\theta $=400，图 4(c)表示引入Dual方法的结果，参数$\lambda $=1, $\theta $=500，图 4(d)表示引入ADMM方法的结果，参数$\lambda $=4, $\theta $=400。从视觉效果上看，几种方法结果基本一致，目标与背景光流场边缘保持较好。背景处岩石、草地的光流场边缘可以明显的区分，图 4中近处的树木光流场边缘非常明显，较远处的树木光流场边缘也可以与背景区分。说明引入快速方法后，结果不会受到较大影响，仍可以计算出准确光流。

图 4 Grove2序列使用不同方法得到的光流场

Fig. 4 Grove2 sequence and optical flow got by different methods ((a) M-HS-TV; (b) Split Bregman; (c) Dual; (d) ADMM)

图 5表示对Grove2图像序列使用不同方法得到的能量函数值随迭代步数变化曲线。观察图 5可知，Split Bregman方法和ADMM方法函数变化曲线几乎完全重合，3种快速方法在50步以内能量函数值变化就趋于平稳状态。而M-HS-TV模型在50步之后仍未收敛。图 3和图 5说明引用快速方法后能量函数值会迅速下降。

图 5 Grove2图像不同方法的能量函数变化曲线

Fig. 5 Decrease of different methods' function value on increasing iteration

表 2表示最高分辨率分层使用不同方法对图 1 Grove2序列运算占用的CPU时间和AAE。从结果上可以看出，3种快速方法均可以在较短时间内得到AAE较低的结果，快速方法中ADMM方法迭代步数较少、AAE较低，综合性能最好。表 1和表 2说明快速方法可以通过减少迭代步骤从而减少算法总时间。

表 2 不同方法各项数据的比较
Table 2 Comparison between different methods

下载CSV

方法	AAE	迭代步数	CPU时间/s
M-HS-TV	5.246 9	370	49
Split Bregman	5.112 8	20	14
Dual	5.231 0	50	21
ADMM	5.101 4	20	13

通过表 1和表 2数据可知，使用传统方法通常需要上百步迭代才能得到最终结果，而引入快速方法仅需几十步迭代。不同图像序列计算时所需迭代步数不同，迭代步数越多快速方法性能越好。

4 结论

本文将几种快速算法应用于大位移变分光流计算中，由于变分光流模型较多，且快速方法具有通用性，因此选取较为简单的M-HS-TV模型为例。实验结果证明，对于固定场景和转动场景图像序列，几种快速方法均能在较短时间内得到较好的光流场，适用范围较广，3种方法中ADMM方法综合性能最好。在实际应用中，为了得到误差较低的光流场，通常需要使用复杂模型，而复杂模型求解困难，且计算效率较低，使用快速方法可以解决这一问题。然而，快速方法的性能会受到一些因素的影响。使用传统方法计算光流时通常需要几百步的迭代步数，对于不同图像序列所需总迭代步数不同。使用快速方法可以迭代步数缩短至几十步，但是减少的时间与使用传统方法时所需的总迭代步数有关。总迭代步数越多快速方法性能越好。总迭代步数减少时，虽然快速方法也能提升一倍以上效率，但是相对来说性能有所下降。

本文主要从算法角度提升模型效率，下一步可以从数值实现方式等角度进行分析和研究。

参考文献

[1] Horn B K P, Schunck B G. Determining optical flow[J]. Artificial Intelligence , 1980, 17 (1-3) : 185–203. DOI:10.1016/0004-3702(81)90024-2

[2] Lucas B D, Kanade T. An iterative image registration technique with an application to stereo vision[C]//Proceedings of 7th International Joint Conference on Artificial Intelligence. San Francisco, CA:Morgan Kaufmann Publishers Inc, 1981:674-679.

[3] Yuan G W, Chen Z Q, Gong J, et al. A moving object detection algorithm based on a combination of optical flow and three-frame difference[J]. Journal of Chinese Computer Systems , 2013, 34 (3) : 668–671. [ 袁国武, 陈志强, 龚健, 等. 一种结合光流法与三帧差分法的运动目标检测算法[J]. 小型微型计算机系统 , 2013, 34 (3) : 668–671. DOI:10.3969/j.issn.1000-1220.2013.03.047 ]

[4] Anandan P. A computational framework and an algorithm for the measurement of visual motion[J]. International Journal of Computer Vision , 1989, 2 (3) : 283–310. DOI:10.1007/BF00158167

[5] Yu C. A research of optical flow algorithm based on block matching[D]. Nanchang:Nanchang Hangkong University, 2014. [余聪.基于块匹配的光流计算方法研究[D].南昌:南昌航空大学, 2014.] http://cdmd.cnki.com.cn/Article/CDMD-10406-1014059146.htm

[6] Wu Y T, Kanade T, Cohn J, et al. Optical flow estimation using wavelet motion model[C]//Proceedings of the 6th International Conference on Computer Vision. Bombay, India:IEEE, 1998:992-992.[DOI:10.1109/ICCV.1998.710837]

[7] Duan X H, Wang Y Q, Wang P A, et al. Estimation of optical flow for phase image using Gabor filter banks[J]. Journal of Computer-Aided Design & Computer Graphics , 2005, 17 (10) : 2157–2161. [ 段先华, 王元全, 王平安, 等. 利用Gabor滤波的相位图像进行光流估计[J]. 计算机辅助设计与图形学学报 , 2005, 17 (10) : 2157–2161. DOI:10.3321/j.issn:1003-9775.2005.10.003 ]

[8] Kim J, Sikora T. Hybrid recursive energy-based method for robust optical flow on large motion fields[C]//Proceedings of the 2005 IEEE International Conference on Image Processing. Genova, Italy:IEEE, 2005:I-129-32.[DOI:10.1109/ICIP.2005.1529704]

[9] Barron J L, Fleet D J, Beauchemin S S. Performance of optical flow techniques[J]. International Journal of Computer Vision , 1994, 12 (1) : 43–77. DOI:10.1109/CVPR.1992.223269

[10] Mei G H, Chen Z, Wei S G, et al. New algorithm for estimation of variational optical flow with image-and flow-driven[J]. Journal of Image and Graphics , 2011, 16 (12) : 2159–2168. [ 梅广辉, 陈震, 危水根, 等. 图像光流联合驱动的变分光流计算新方法[J]. 中国图象图形学报 , 2011, 16 (12) : 2159–2168. DOI:10.11834/jig.20111210 ]

[11] Papenberg N, Bruhn A, Brox T, et al. Highly accurate optic flow computation with theoretically justified warping[J]. International Journal of Computer Vision , 2006, 67 (2) : 141–158. DOI:10.1007/s11263-005-3960-y

[12] Bruhn A. Variational optic flow computation:accurate modelling and efficient numerics[D]. Germany:Saarland University, 2006. http://www.oalib.com/references/17190094

[13] Brox T, Bregler C, Malik J. Large displacement optical flow[C]//Proceedings of 2009 IEEE Conference on Computer Vision and Pattern Recognition. Piscataway:IEEE, 2009:41-48.[DOI:10.1109/CVPR.2009.5206697]

[14] Aubert G, Deriche R, Kornprobst P. Computing optical flow via variational techniques[J]. SIAM Journal on Applied Mathematics , 1999, 60 (1) : 156–182. DOI:10.1137/S0036139998340170

[15] Liu T S, Shen L X. Fluid flow and optical flow[J]. Journal of Fluid Mechanics , 2008, 614 : 253–291. DOI:10.1017/S0022112008003273

[16] Wang B, Cai Z M, Shen L X, et al. An analysis of physics-based optical flow[J]. Journal of Computational and Applied Mathematics , 2015, 276 : 62–80. DOI:10.1016/j.cam.2014.08.020

[17] Fortun D, Bouthemy P, Kervrann C. Optical flow modeling and computation:a survey[J]. Computer Vision and Image Understanding , 2015, 134 : 1–21. DOI:10.1016/j.cviu.2015.02.008

[18] Zhang S F, Zhang W S, Ding H, et al. Background modeling and object detecting based on optical flow velocity field[J]. Journal of Image and Graphics , 2011, 16 (2) : 236–243. [ 张水发, 张文生, 丁欢, 等. 融合光流速度与背景建模的目标检测方法[J]. 中国图象图形学报 , 2011, 16 (2) : 236–243. DOI:10.11834/jig.20110220 ]

[19] Lauze F, Kornprobst P, Mémin E. A course to fine multiscale approach for linear least squares optical flow estimation[C]//Proceedings of the 2004 British Machine Vision Conference. London, UK:BMVA Press, 2004:1-10.[DOI:10.5244/C.18.79]

[20] Yan W J. Research on estimation and performance evaluation of optical flow from image sequence[D]. Nanchang:Nanchang Hangkong University, 2012. [晏文敬.图像序列光流计算及其性能评估的研究[D].南昌:南昌航空大学, 2012.] http://cdmd.cnki.com.cn/Article/CDMD-10406-1015968491.htm

[21] Goldstein T, Osher S. The split Bregman method for L1-regularized problems[J]. SIAM Journal on Imaging Sciences , 2009, 2 (2) : 323–343. DOI:10.1137/080725891

[22] Chambolle A. An algorithm for total variation minimization and applications[J]. Journal of Mathematical Imaging and Vision , 2004, 20 (1-2) : 89–97. DOI:10.1023/B:JMIV.0000011325.36760.1e

[23] Zach C, Pock T, Bischof H. A duality based approach for realtime TV-L¹ optical flow[C]//29th DAGM Symposium. Berlin Heidelberg:Springer, 2007:214-223.[DOI:10.1007/978-3-540-74936-3_22]

[24] Wu C L, Tai X C. Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models[J]. SIAM Journal on Imaging Sciences , 2010, 3 (3) : 300–339. DOI:10.1137/090767558

[25] Goldstein T, O'Donoghue B, Setzer S, et al. Fast alternating direction optimization methods[J]. SIAM Journal on Imaging Sciences , 2014, 7 (3) : 1588–1623. DOI:10.1137/120896219

[26] Wei W B, Pan Z K, Li Y Y, et al. Median formula for optic flow computation of small displacement[J]. Chinese Journal of Scientific Instrument , 2011, 32 (10) : 2256–2260. [ 魏伟波, 潘振宽, 李媛媛, 等. 小位移光流计算的中值公式[J]. 仪器仪表学报 , 2011, 32 (10) : 2256–2260. ]

[27] Yu J J, Pan Z K, Wei W B. Dual method for small displacement optic flow computation[J]. Computer Engineering , 2010, 36 (7) : 260–261. [ 于晶晶, 潘振宽, 魏伟波. 小位移光流计算的对偶方法[J]. 计算机工程 , 2010, 36 (7) : 260–261. DOI:10.3969/j.issn.1000-3428.2010.07.089 ]

[28] Rudin L I, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms[J]. Physica D:Nonlinear Phenomena , 1992, 60 (1-4) : 259–268. DOI:10.1016/0167-2789(92)90242-F

[29] Duval V, Aujol J F, Vese L. A projected gradient algorithm for color image decomposition[EB/OL].[2016-09-30]. https://www.math.u-bordeaux.fr/~jaujol/PAPERS/CMLA2008-21.pdf

[30] Middlebury.Flowdatabase[DB/OL].[2016-09-30]. http://vision.middlebury.edu/flow/data/. [Middlebury光流数据库[DB/OL].[2016-09-30]. http://vision.middlebury.edu/flow/data/.]