图像处理中的全局优化技术（经典至极）

图像处理中的全局优化技术(Global optimization techniques in image processing and computer vision) (一)

2013-05-29 14:26 1659人阅读 评论(1) 收藏 举报

算法图像处理计算机视觉imagevision

MulinB按：最近打算好好学习一下几种图像处理和计算机视觉中常用的 global optimization (或 energy minimization) 方法，这里总结一下学习心得。分为以下几篇： 1. Discrete Optimization: Graph Cuts and Belief Propagation (本篇) 2. Quadratic Optimization : Poisson Equation and Laplacian Matrix 3. Variational Methods for Optical Flow Estimation

4. TODO: Likelihood Maximization (e.g., Blind Deconvolution)

1. Discrete Optimization: Graph Cuts and Belief Propagation

很多图像处理和视觉的问题可以看成是pixel-labeling问题，用energy minimization framework可以formulate成以下能量方程：

其中第一项叫data term (或叫unary term)，是将label l_p赋给像素p时的cost，第二项叫smoothness term (或叫pairwise term)，是每两个相邻pixel的labeling不同时产生的cost (w_pq是spatial varying weight，比如跟p和q的difference相关)。传统的smoothness term一般只考虑两两(pairwise)相邻的pixel，最近几年的CVPR和ICCV上开始出现很多higher-order MRF的研究，比如这位大牛的paper，这是题外话。这种energy minimization framework其实从概率的角度看，等价于求Markov Random Field (MRF) 的maximum a posteriori (MAP) 概率。求解这类energy minimization的方法，最流行的有两个，Graph Cuts和Belief Propagation。

1.1 Graph Cuts

刚开始学习Graph Cuts时，不知道到底这方法是从哪篇paper最早提出来的，因为在后来的paper里引用的参考文献一般指向好几个来源，这里先大致梳理一下参考文献。一般认为，将Graph Cuts引入图像处理领域的先驱paper是Greig等人1989年发表的[1]，不过貌似没有引起太大的注意。真正使Graph Cuts在图像领域风靡起来的是两个俄罗斯人(貌似是在Cornell做出的成果)，Yuri Boykov和Vladimir Kolmogorov，以及他们的导师Ramin Zabih。首先引起大家注意的是Boykov在ICCV 2001上的使用Graph Cuts解决Interactive Image Segmentation的paper[2]，以及这篇paper提到的一个max-flow算法[3] (该max-flow算法最早是发在2001年的一个CVPR Workshop上，后来扩展到TPAMI [3])。需要注意的是，这两篇paper里的Graph Cuts算法，只是针对只有两个label (待求变量是binary variable)的情况。而Boykov的2001 TPAMI paper[4]提出使用alpha expansion技术将多个label的问题转化成一系列的binary label问题去逼近，这使得Graph Cuts算法开始风靡起来。后来Kolmogorov的2004 TPAMI paper[5] 进一步讨论了什么样的能量方程可以被Graph Cuts优化，并给出了一个简单清晰的根据能量方程构造相应graph的算法，该算法基本成为被大家广泛使用的Graph Cuts算法。Boykov和Kolmogorov的代码可以从这里找到。

下面简单介绍一下Graph Cuts算法，先从binary label开始(见参考文献[2] [3])。顾名思义，Graph Cuts是将图像中的所有pixel以及两个label作为node，根据data cost和smoothness cost建立node之间的edge，这样构造一个(无向)graph，然后通过cut算法将整个graph切成两个分离的部分。如下图所示：

注意，图中的cut会切断它经过的所有edge(包括蓝色、红色、和土黄色的edge)。如果将两个label的node看成两个特殊的terminal，这样的一个cut会阻断所有s连往t的路径(edge)。在Boykov的ICCV 2001 paper [2]中，他证明了通过简单的方式构造这样的一个graph，如果能找到一个min-cut (即该cut经过的edge cost加起来在所有possible的cut中最小)，其实就是上面的能量方程的最小解 (见paper中的Theorem 1)。那么如何找到min-cut呢？

在图论里，有证明找到min-cut其实等价于找到max-flow，即从s流往t的最大流量。其实，min-cut等价于max-flow的intuition很简单，从一个terminal流往另一个terminal的最大流量，其瓶颈肯定是min-cut的位置。这里有个有意思的介绍，关于网络里s-t flow的计算。计算max-flow的经典算法主要有两种，一种是基于augmenting-path，一种是基于push-relabeling。在Boykov和Kolmogorov的TPAMI 2004 paper [3]里，介绍了一种基于augmenting-path，为了图像这种扁平graph量身定制的max-flow算法，通过实验证明了其效率，这里有他们的代码。

在解决multi-labeling问题时 (其实是更为普遍和常见的问题)，在能量方程满足某些特定的条件下(注意：该限定条件其实挺难满足，后面讨论！)，可以使用alpha expansion算法将其转化为一系列binary-labeling问题来逼近，参见Boykov TPAMI 2001 paper [4]。见下图：

这种alpha expansion思路很简单，当处理每一个label时(假设其为a)，将其他所有的label看成一个label package(假如称之为b)，这时问题就变成了binary-labeling。此时在进行cut时，如果一个原来是a的pixel被cut给b，将无法确定到底给该pixel具体哪一个label(由

于b是个大杂烩)。所以在进行cut时，只允许原来是b的pixel被cut给a，也就是标记为a的pixel在expanding，这就是算法名字的来源。需要注意的是，为了使得这样的一次alpha expansion可以被max-flow算法计算出来，graph的构造比之前的binary-labeling要稍微复杂一些 (比如仅仅允许alpha expansion的话，有些跟b相连的edge weight要设成无穷大)。使用alpha expansion算法的步骤很简单，如下：

[cpp] view plaincopy

1. // alpha expansion algorithm pseudo-code 2. initialize labeling; 3.

4. while not converged 5. {

6. for each label a in L 7. {

8. construct a graph; 9. do max-flow cut;

10. if energy is smaller than before, accept it; 11. else decline it; 12. } 13. }

值得注意的是，这种alpha expansion只是multi-labeling问题的近似求解，而之前的max-flow算法是binary-labeling问题的exact求解方法。而且，为了使得这种alpha expansion时的graph可以被构造出来，能量方程需要满足一定的限定条件，具体来说，是能量方程中的pairwise term函数V_pq需要满足某些限定条件。在Boykov TPAMI 2001 paper [4]中称之为V_pq必须是一个metric(类似于满足“距离”的定义，比如：可交换、满足三角不等式)。在Kolmogorov TPAMI 2004 paper [5] 中将其推广为V_pq必须是一个submodular函数(文中称之为regular，其实后来都称之为submodular)，即函数V_pq必须满足V(0,0) + V(1,1) <= V(0,1) + V(1,0)条件。该条件乍一看貌似很容易满足，特别是对于binary-labeling来说。然而，注意，在alpha expansion中，该条件变成了如下，

注意，其中l_p和l_q可能不一样。这样一来其实该条件没那么容易满足！比如常用的quadratic cost就不满足！常用的pairwise term的V_pq函数如下表所示 (Potts Model中delta函数是unit impulse函数)：

例如，quadratic cost时，l_p=3, l_q=10, a=5, 这时，上面条件里左侧第一项的值是49，第二项是0，不等式右侧第一项是4，第二项是25，显然不等式不成立！

另外需要一提的是，其实即使上述的条件不成立，仍然可以使用alpha expansion，使用的时候可以对一些不满足条件的项进行修改，这种技术在CVPR 05的一篇paper里提出[6]，叫做truncating，在这里的代码(文件energy.h，函数add_term2里)可以找到例子，代码非常简单：

[cpp] view plaincopy

1. //Truncating for non-submodualr term, code by Kolmogorov 2. //A, D, C, B分别是上述不等式的四项 3. if ( A+D > C+B) 4. {

5. Value delta = A+D-C-B; 6. Value subtrA = delta/3; 7.

8. A = A-subtrA; 9. C = C+subtrA;

10. B = B+(delta-subtrA*2); 11. }

这种truncating之后的算法其实并不能保证最终结果是strong local minimum (注意，没有truncating的alpha expansion只是能保证找到strong local minimum，不能保证是global mininum)，但是实际使用中效果不错。另外专门针对non-submodular进行优化的算法有QBPO，见Kolgoromov TPAMI 2007 paper[7]。

1.2 Belief Propagation

Belief Propagation是Graph Cuts的一个强劲对手，其渊源可以追溯到MIT在ICCV 2003上的paper [8]，比较这两种方法在stereo上的性能，结果貌似不分上下，Graph Cuts略胜一点点。后来大牛Richard Szeliski联合一堆MRF的experts搞了一个benchmark评测这些方法，见TPAMI 2008 paper[9]，和这个benchmark，上面有code可以下载，集成了很多算法，非常值得研究。结论是，Graph Cuts (alpha expansion)算法表现比较出色，另外Belief Propagation的一个改进算法叫Tree ReWeighted Message Passing (TRW-S)也表现出色，而Belief Propagation似乎表现不那么理想。不过其实其中除了Photomontage之外的其他比较中，基本都难分高下(在Photomontage中，alpha expansion完爆其他算法)。在实际使用中，至少在stereo的Middlebury benchmark上，Belief Propagation还是占了大多数席位。跟Graph Cuts相比，Belief Propagation至少有以下几个优点： (1). 对能量方程不存在submodular限制;

(2). implement非常简单 (虽然理论比较难懂，但是算法implement非常之简单);

(3). 快速算法很快 (比如Hierarchical Belief Propagation，见Pedro Felzenszwalb的IJCV 2006 paper[10]，有代码在此).

下面简单介绍一下Belief Propagation算法，其背景知识可以参考Coursera上Daphne Koller的Probabilistic Graphical Model。Belief Propagation (BP)算法最早可以追溯到Pearl的书 [11]，该算法最初是设计在tree或者graph without cycles上通过message passing的方式计算MAP的(也就是上面提到的energy minimization)，在这种情况下，BP可以计算出exact的结果(其实等价于dynamic programming算法)。而在更为普遍的情况下，即graph中有cycle时，BP并不能保证可以converge，但是在很多情况下，可以converge到strong local minimum (有点类似alpha expansion)，见MIT大牛的TIT 2001 paper [12]。在这种情况下，BP一般被称为Loopy Belief Propagation，即graph中存在cycle使得message passing存在loop。另外，一般提到BP有两种，一种是max-product，一种是sum-product，前者是用来计算MAP的(max of probability product)，后者是用来计算marginal propability的，我们这里只讨论max-product (其实sum-product算法类似，只是计算message时稍有不同)。BP算法非常简单，流程如下：

[cpp] view plaincopy

1. //Belief Propagation Algorithm pseudo-code, suppose run T iterations 2. for (int t=0; t

4. update message between each neighboring pixels; 5. } 6.

7. select label for each pixel such that the sum of data cost and messages pass

ing to it is minimum;

最主要的步骤在于每次iteration中得message update。在计算某个pixel传递给另一个pixel的message时，要考虑上一个iteration中其他邻居pixel传给该pixel的message，如下图所示：

其中，蓝色箭头表示上一个iteration的message，红色箭头表示当前要计算的message。那么，message到底如何计算呢？假设label的candidate有L个(回忆一下，我们要解决的是pixel-labeling问题)，那么每条message包含有L个数值entry，表示的信息是“嘿，邻居，从我当前的状态看，你选择每个label的代价分别是...”。如果用f_q表示每个label，那么p传给q的message中的每个entry如下计算：

其中，t表示迭代数，f_p表示p的possible label(也有L个)，s是p的其余三个neighbors。非常简单，当要计算从p传给q的消息中关于f_q的entry时，考察每一个f_p，把data cost，smoothness cost (从f_p到f_q的cost), 其余neighbors传来的message统统加起来，在其中选择最小值即可(类似于dynamic programming算法)。图示如下：

可以看到，这个过程计算复杂度是O(L*L)。在Pedro Felzenszwalb的IJCV 2006 paper[10]里，针对一些常见的smoothness cost介绍了一些快速算法可以将复杂度减少到O(L)的计算复杂度，基本涵盖了上面表格里列出的cost类型。这样的message update经过T次迭代后，最后通过如下公式计算每个pixel的cost，

选择cost最小的label即为最终结果。

上面的message passing方式其实比较慢，一般需要很多次iteration才能使得message传递到全图。在Pedro Felzenszwalb的IJCV 2006 paper[10]里，提出使用

multiscale/hierarchical的方式使得message更加快速的进行传递，这使得需要的iteration数量大大减小，使得BP的速度可以很快。

1.3 Other Related Techniques

Graph Cuts是从处理Interactive Image Segmentation (或叫Seeded Image Segmentation) 问题开始发家的 [2]，其他几种可以处理这种问题的方法有：Random Walker, Geodesic Shortest Path, Level Set Method等。这里简单介绍一下。

首先值得一提的是，SIGGRAPH 2004井喷了三篇基于Graph Cuts做interactive image editing的文章：GrabCut[13], Photomontage [14], Lazy Snapping [15]。其中以GrabCut影响力最大。GrabCut其实只是对Graph Cuts (Boykov ICCV 2001 [2]) 的一个扩展，使其变得更加robust。第一项扩展是将计算data cost从histogram扩展成用Gaussian Mixture Model，第二项扩展是将一次Graph Cuts计算扩展成迭代多次Graph Cuts。由于这样算法变得更加robust，其用户输入也可以变得较少，直接一个框就够了(当然使用fg/bg brush也可以，而且他们的system也可以用matting brush进行边界matting)，比如下图。

Random Walker图像分割(segmentation或matting)算法是另一种基于用户输入seeds (或者scribble) 的分割算法，见Leo Grady TPAMI 2006 paper[16]，注意，该算法在有些文献里叫做Random Walks图像分割算法，不过后来其作者倾向于叫其Random Walker算法，估计为了和随机游走算法(Random Walks)区分。其实该算法只是以Random Walks/Walker为出发点导出的，最终算法是一个quadratic optimization，有closed-form solution，可以通过求解一个线性方程组解决，所以算法可以很快(虽然数值求解大型线性方程组可能需要用到迭代算法，但是相当于其他迭代算法如Graph Cuts和Belief Propagation而言，求解线性方程组还是非常快的)。具体来讲，Random Walker的Motivation是：将图像看成是一个graph，给定用户输入的seeds，假设在每个pixel放置一个random walker，那么该random walker最先到达哪一个seed(概率上)，就将该pixel标记为那个seed的label。有相关理论证明(貌似一些关于potential theory或者circuit theory，whatever...)，求解这样的random walker模拟等价于求解Dirichlet problem，通过minimize相应的Dirichlet Integral，可以求得其解(即称为harmonic function)。而在graph上formulate该Dirichlet Integral，最后形式非常简单(Grady称之为combinatorial formulation，貌似其他的文献里有称之为discrete finite differential formulation)，这样，问题最终归结求解下面的非常简单的quadratic minimization问题：

其中，l是待求的labeling，l_F和l_B是用户标记的seeds处的label。其实该问题与后面我要谈的quadratic optimization里很多formulation很相似，尤其是Levin的Colorization和Closed-form Matting，这是后话。另外，该formulation还出现在diffusion map[17] 的问题中，这也是后话。Leo Grady后来还将Random Walker和Graph Cuts联系了起来 (ICCV 2007 paper[18])，其实从上面的formulation也可以大概看出来，上面的objective function只有smoothness term，而将其constraint移到objective function里就成为data term，跟

Graph Cuts的唯一区别就是，Random Walker要解决的能量方程smoothness term是quadratic，而Graph Cuts可以用在不同的smoothness term的能量方程中(见上文的表格)。Random Walker有一个引人注意的特点，就是在weak boundary的表现非常出色(估计是因为其motivation跟potential有关系)，见下图例子。

从Random Walker的motivation可以想到，如果直接计算从每个pixel出发到seeds的最短路径(在路径上加入跟图像相关的weight时，这种路径为geodesic)，然后看pixel距离哪个seed近，就将pixel标记为那个seed的label。这就是基于Geodesic Distance的分割，见ICCV 2007 paper[19]。这种算法有个优点，就是所有pixel距离某一种seeds的geodesic distance可以非常快速的计算出来，用N表示pixel的个数，那么用Yatziv的Fast

Marching[20] (在[19]中使用的算法，是Sethian的Fast Marching Method的快速实现)，或者更适用于并行计算的Parallel Marching [21](在[22]中使用的算法，这里可以找到代码)，计算复杂度都可以达到O(N)。跟Random Walker相比，基于Geodesic Distance的方法速度更快，不过对seeds的位置依赖比较严重，而且在weak boundary表现也不会好，比如下面的例子。

提到Fast Marching Method，就不得不提另一个在图像领域影响很大的方法：Level Set Method。Level Set Method最早提出是为了解决boundary evolve的问题，比如用户在一个图像上画一个圈，然后这个圈根据图像内容进行演化，最终将一个物体圈出来，比如下图：

这种boundary evolve的问题最原始的方法是explicitly去跟踪boundary上的一些control points，而如果boudanry变化的过程中发生topological结构变化时(比如一个圈圈分裂成两个圈圈)，这种跟踪方式会很变得很难。Level Set Method就是为了解决这种难题。简单来讲，Level Set Method大概思路就是把原来二维的boundary evolve问题重新参数化为三维的surface evolve问题，新参数为笛卡尔坐标x-y以及时间t，这样新的问题成为一个在笛卡尔grid(即pixel grid)中的PDE，一般称为Hamilton-Jacob Equation，可以使用数值解法求出。其原理介绍可以参考Coursera上Guillermo Sapiro的Image and Video Processing，个人认为他讲的非常清晰易懂(可以直接看Section 6: Geometric PDEs)。

最后，再提一下Fast Marching Method。当上面的boundary evolve中曲线法向速度符号不变时(比如说boundary一直朝外扩张)，我们可以用更快速的方法来求解这个PDE：从boundary出发，在其所在的笛卡尔grid里计算并记录从boundary到每个pixel的距离(即到达时间)，如此得到一个类似于等高线图的map，而从这个map上，其实可以得到任意t时的boundary (取记录了某一个相同时间/距离的pixels)，如下图所示。这种map可以使用类似于Dijkstra最短路径算法的方法计算(一般称之为weve-front propagation方法)，或者用raster-scan的方法计算(更易于并行化)。Sethian的Fast Marching Method以及上文提到的Yatziv的快速实现[20]属于前者，而并行算法[21]属于后者。

1.4 Reference

[1] D. Greig, B. Porteous, and A. Seheult, “Exact Maximum A Posteriori Estimation for Binary Images,” J. Royal Statistical Soc., 1989.

[2] Y. Boykov and M.-P. Jolly, “Interactive Graph Cuts for Optimal Boundary and Region Segmentation of Objects in N-D Images,” ICCV, 2001.

[3] Y. Boykov and V. Kolmogorov, “An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision,” TPAMI, 2004.

[4] Y. Boykov, O. Veksler, and R. Zabih, “Fast Approximate Energy Minimization via Graph Cuts,” TPAMI, 2001.

[5] V. Kolmogorov and R. Zabih, “What Energy Functions can be Minimized via Graph Cuts,” TPAMI, 2004.

[6] C. Rother, S. Kumar, V. Kolmogorov, and A. Blake, “Digital Tapestry,” CVPR, 2005. [7] V. Kolmogorov and C. Rother, “Minimizing Nonsubmodular Functions with Graph Cuts—A Review,” TPAMI, 2007.

[8] M. Tappen and W. Freeman, “Comparison of Graph Cuts with Belief Propagation for Stereo, Using Identical MRF Parameters,” ICCV, 2003.

[9] R. Szeliski, R. Zabih, D. Scharstein, O. Veksler, V. Kolmogorov, A. Agarwala, M. Tappen, and C. Rother. “A comparative study of energy minimization methods for markov random fields with smoothness-based priors.” TPAMI, 2008

[10] P. Felzenszwalb and D. Huttenlocher, “Efficient Belief Propagation for Early Vision,” IJCV, 2006.

[11] J. Pearl, Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann, 1988.

[12] Weiss, Y. and Freeman,W.T. 2001. “On the optimality of solutions of the max-product belief propagation algorithm in arbitrary graphs.” IEEE Transactions on Information Theory, 2001.

[13] C. Rother, V. Kolmogorov, and A. Blake, “?GrabCut?—Interactive Foreground Extraction Using Iterated Graph Cuts,” SIGGRAPH, 2004.

[14] A. Agarwala,M. Dontcheva, M. Agrawala, S. Drucker, A. Colburn, B. Curless, D. Salesin, and M. Cohen, “Interactive Digital Photomontage,” SIGGRAPH, 2004. [15] Y. Li, J. Sun, C.-K. Tang, and H.-Y. Shum, “Lazy snapping,” SIGGRAPH, 2004. [16] L. Grady, “Random Walks for Image Segmentation,” TPAMI, 2006.

[17] R.R. Coifman, S. Lafon, A.B. Lee, M. Maggioni, B. Nadler, F. Warner, and S.W. Zucker, “Geometric Diffusions as a Tool for Harmonic Analysis and Structure Definition of Data: Diffusion Maps,” Proc. Nat?l Academy of Sciences USA, 2005.

[18] A.K. Sinop and L. Grady, “A Seeded Image Segmentation Framework Unifying Graph Cuts and Random Walker Which Yields a New Algorithm,” ICCV, 2007.

[19] X. Bai and G. Sapiro, “A Geodesic Framework for Fast Interactive Image and Video Segmentation and Matting,” ICCV, 2007.

[20] L. Yatziv, A. Bartesaghi, and G. Sapiro, “O(n) implementation of the fast marching algorithm,” Journal of Computational Physics, 2006.

[21] O. Weber, Y. Devir, A. Bronstein, M. Bronstein, and R. Kimmel,“Parallel algorithms for approximation of distance maps on parametric surfaces,” SIGGRAPH, 2008. [22] A. Criminisi, T. Sharp, and C. Rother, “Geodesic Image and Video Editing,” TOG, 2010.

图像处理中的全局优化技术(Global optimization techniques in image processing and computer vision) (二)

2013-06-12 16:36 1556人阅读 评论(0) 收藏 举报

算法图像处理计算机视觉imagevision

MulinB按：最近打算好好学习一下几种图像处理和计算机视觉中常用的 global optimization (或 energy minimization) 方法，这里总结一下学习心得。分为以下几篇： 1. Discrete Optimization: Graph Cuts and Belief Propagation

2. Quadratic Optimization: Poisson Equation and Laplacian Matrix (本篇) 3. Variational Methods for Optical Flow Estimation

4. TODO: Likelihood Maximization (e.g., Blind Deconvolution)

2. Quadratic Optimization: Poisson Equation and Laplacian Matrix

Quadratic Optimization (Least Squares Minimization)在图像处理中的魅力要从

SIGGRAPH 02和03年的两篇Gradient Domain Image Editing文章说起：Fattal的HDR Compression[1] 和Perez的Poisson Image Editing [2]。 其后，Levin的两篇文章Colorization [3]和Closed-Form Matting [4]更是将其魅力展现的淋漓尽致。而Farbman的基于Weighted Least Squares的WLS filter [5]也是在Edge-preserving Filter领域名声大噪。由于目标函数是quadratic, 这类问题的求解一般比较容易，大多都可以最终归结为求解一个大型稀疏线性方程组。而数值求解大型线性方程组是一个由来已久的问题，有着各种现成的solver，更是有着为以上这类问题量身定做的solver，见下文solver小节。

题外话：以色列的耶路撒冷希伯来大学(The Hebrew University of Jerusalem)的Lischinski教授貌似很偏爱这类方法，上面提到的这些文章大多有他的署名。

2.1 Problem I: Gradient Domain Image Editing

有心理学为证(见Poisson Image Editing [2]文章的introduction部分)，对图像的gradient进行修改可以产生比较不容易感知到的artifacts，这使得很多图像编辑的工作可以放到gradient domain使得效果很逼真，比如下图的图像拼合例子(图例来自[2])：

其实对gradient domain进行修改而获得逼真的编辑效果由来已久，最早见于1983年Burt-Adelson的Laplacian Pyramid [6]图像融合(这里有个简洁的中文介绍)，这是题外话。在gradient domain进行图像编辑的pipeline一般如下(图例修改自ICCV 2007 Course -- Gradient Domain Manipulation Techniques，顺便赞一个，nice ppt!)：

其中第一步的gradient processing根据不同的需求有具体的操作，比如HDR Compression里是将较大的gradient value进行削弱，而上面的图像拼合例子(Seamless Clone)则是将源图像的gradient拷贝到目标区域。而其中第二步中由gradient重建出新图像并非那么容易，因为经过编辑后的gradient一般是不可积分的，这时Quadratic Optimization粉墨登场。

假设待求图像为I，修改后的已知gradient是G，则通过Least Squares Minimization可以将问题formulate成如下(使得待求图像I的gradient在L2 norm下尽量接近G)：

注意其中的约束条件，比如，在图像拼合例子中，非编辑区域的像素是已知的，在求解编辑区域的像素时，边界上的像素值是约束条件。上面的formulation是假设图像I是定义在x-y连续空间的函数，所以其实上述目标函数是关于I的functional(泛函，也就是“函数的函数”)。使用calculus of variation(变分法)中的Euler-Lagrange Equation (one unknown function, two variables)可以将其转化为一个非常经典的偏微分方程形式，这就是Poisson Equation：

注意其中G是已知的，I是未知的，Δ是Laplacian operator，div是divergence operator。当已知边界像素值时，该偏微分方程具有第一类边界条件(Dirichlet boundary condition)，比如图像拼合；当处理整个图像时，该偏微分方程具有第二类边界条件(Neumann boundary condition)，即已知边界导数值(设为0)，比如HDR Compression。

上面的formulation是在x-y连续空间(像素坐标是连续的)，而用于图像处理时，一般需要将其离散化(因为事实上像素坐标(x,y)是整数)，上面相应的偏微分形式可以使用有限差分(finite difference)形式近似代替。具体来讲，离散化的discrete Laplacian operator如下，

而divergence operator中的一阶偏导可以用前向或者后向差分近似(由于G本身是由gradient得来，一般如果之前计算gradient使用前向差分，那么这里计算div就使用后向差分，这样使得两次差分的结果等价于一次二阶中心差分，具体参考[1])，比如这里的divG可以由以下后向差分近似，

于是，整个Poisson Equation离散化之后，每个pixel都有一个线性方程，假设图像有N个pixel，那么整个Poisson Equation就成了一个包含N个方程的大型方程组。如果将这个大型方程组写成矩阵形式(假设将待求图像I拉成一个长的vector，用x表示，将已知的divG也拉成一个长的vector，用b表示)，离散化的Poisson Equation变成了经典的Ax=b形式。以5×5的图像为例，假设待求图像I为如下形式(每个pixel的值都是未知)：

将其拉成列向量x(按列展开)，则整个Discrete Poisson Equation (with Neumann boundary condition)写成Ax=b形式即，

该矩阵A可以直接从Laplacian operator得来，一般称为Laplacian Matrix(其实如果将图像看成graph，该矩阵即为graph theory中的Laplacian Matrix)。注意，对角线上值为2和3的元素是对应在图像边界上的pixel(因为其discrete Laplacian operator无法完整展开，包含了一些不存在的neighboring pixel)，如果将边界条件改为Dirichlet boundary condition并且未知区域周边的pixle都是已知的话，对角线上的元素就全为4，比如下面的例子。假设待求图像I为如下形式(未知pixel的周边pixel是已知)：

则未知向量x包含9个元素，整个Discrete Poisson Equation (with Dirichlet boundary condition)变成以下形式，

注意等式右侧包含了边界已知pixel的值。这里的Laplacian Matrix较为规整，主要是因为所有的未知pixel处的discrete Laplacian operator可以完整的展开。

总之，上面的discrete Poisson Equation都可以归结为求解一个大型线性方程组，其中的Laplacian Matrix具有以下特点： 1. 对角线元素为non-negative； 2. 非对角线元素为non-positive； 3. 对称(symmetric)；

4. 正定或半正定(positive semidefinite)； 5. 属于分块对角阵。

下面的solver小节再讨论如何求解这样的线性方程组。另外，其实如果直接对上面Least Squares Minimization的目标函数进行离散化，也可以得到相同的方程组(省去了应用Euler-Lagrange Equation的步骤)，参见文章[2]的推导过程。

其中C_desc是descriptor matching的cost。辅助变量u_d和v_d的作用是可以将该目标函数decouple成两个子问题：descriptor matching和variational求解。迭代计算这两个子问题即可求得最终结果，比如，先descriptor matching得到u_d和v_d，然后再用上述variational framework计算，然后再descriptor matching，然后再variational framework…

该integrated framework其实在benchmark上成绩并不突出，只是后来基于此framework做的DeepFlow(ICCV 2013 oral paper)，效果亮瞎了(见Sintel和KITTI benchmark)，值得期待(目前paper还没放出)。

3.4 Reference

[1] B.D. Lucas and T. Kanade, “An Iterative Image Registration Technique with an Application to Stereo Vision,” IJCAI, 1981. [2] B.K.P. Horn and B.G. Schunck, “Determining Optical Flow,” Artificial Intelligence, 1981. [3] M.J. Black and P. Anandan, “The Robust Estimation of Multiple Motions: Parametric and Piecewise-Smooth Flow Fields,” CVIU, 1996. [4] T. Brox, A. Bruhn, N. Papenberg, and J. Weickert, “High Accuracy Optical Flow Estimation Based on a Theory for Warping,” ECCV, 2004.

[5] C. Zach, T. Pock, and H. Bischof, “A Duality Based Approach for Realtime TV-L1 Optical Flow,” DAGM, 2007. [6] A. Wedel, T. Pock, C. Zach, H. Bischof, and D. Cremers, “An Improved Algorithm for TV-L1 Optical Flow Computation,”Dagstuhl Motion Workshop, 2008. [7] D. Sun, S. Roth, and M.J. Black, “Secrets of Optical Flow Estimation and Their Principles,” CVPR, 2010. [8] D. Sun, S. Roth, and M.J. Black, “A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them,” IJCV, 2013. [9] C. Liu, “Beyond Pixels: Exploring New Representations and Applications for Motion Analysis,” PhD Thesis, MIT, 2009. [10] T. Brox and J. Malik, “Large Displacement Optical Flow: Descriptor Matching in Variational Motion Estimation,” TPAMI, 2011.

如果进一步对上述的Gradient Domain Image Editing Pipeline进行推广，将原始图像考虑进来(注意，前面的方法是完全由新的gradient重建图像，不考虑原始图像的影响)，这样可以得到一个更为强大的framework，这就是GradientShop[7] (其实之前Levin的

Colorization[3]，以及Lischinski的Interactive Local Adjustment[8]，和Farbman的WLS filter[5]都可以算作这个framework的特例)，新的框架如下：

假设constraint image为图像S(可以是原始图像本身或者经过处理后的图像、用户输入的像素值等)，新的gradient为G，待求图像为I，则Least Squares Minimization的目标函数推广为如下：

其中lambda为weighting factor，用来权衡constraint image和new gradient对待求图像的影响，可以整个image用统一的lambda值，也可以用一个weighting map为每个pixel赋上不同的lambda值。比如，在Levin的Colorization[3]中，待求的图像I是chrominance通道(比如YUV中的U和V通道)，G=0，S是用户输入的scribbles，lambda_2是scribble

的mask(有scribble为1，否则0)，lambda_1是根据已知灰度图中相邻pixel的affinity进行赋值(pixel灰度值越相似，值越大)，这样获得的Laplacian Matrix可以使得用户输入的scribble根据pixel affinity进行有效的propagation。再以5×5图像为例，假设图像M(元素为m_xy)为S的reverse mask(有输入为0，没有输入为1)，S的元素用s_xy表示(无用户输入处的pixel为0)，用w_{x1y1-x2y2}表示pixel(x1,y1)和pixel(x2,y2)的affinity weight(已经归一化)，则Colorization的整个线性方程组可以表示为如下(右键查看大图)：

可见，其矩阵跟上述Poisson Equation中的矩阵形式一模一样(满足Laplacian Matrix的几个特点，其实也可以称之为Laplacian Matrix)，只是矩阵元素的值不再是固定值，而是跟输入有关，而且是spatially varying。注意：这里使用的neighborhood是四邻域，Levin的代码中使用的是8邻域。

GradientShop framework的另一个特例是WLS filter [5]。令constraint image S为输入图像本身，G=0，而lambda_1是根据输入图像的gradient进行spatially varying weighting：某个pixel的gradient越大，lambda_1越小(期待与输入图像较为接近)；gradient越小，lambda_1越大(进行smoothing)。这样的效果就是进行edge-preserving smoothing，在decomposition-manipulation-recombination时可以避免artifacts的产生。

最后值得一提的是，在Gradient Domain进行操作的其他work有：Diffusion Curves [9],Gradient Domain Painting [10] 等，都是很有趣的东东。另外，没有通过解线性方程组而是通过操作Laplacian Pyramid的local Laplacian filter [11]也是很有趣的work，个人非常喜欢，能生成类似于梦幻般的图像效果并且没有什么artifacts，例如下图。

2.2 Problem II: Closed-Form Matting

Image Matting就是精确抠图问题，主要是用来抠出毛发等具有半透明边界的物体，可以用以下方程描述：

给一个输入图像I，以及部分foreground和background的标记(用Trimap或者scribbles标记，如下图)，目标是求出alpha matte, F和B。

Trimap中白色部分和黑色部分分别是已经标记为F和B的区域，所以matting的主要目的是求出灰色未知区域的alpha matte值(介于0到1之间的值)。基于scribbles的输入相当于稀疏的Trimap，除了白色和黑色的scribbles外，都是未知区域。较为经典的几种算法包括Bayesian Matting，Poisson Matting，Closed-Form Matting，Robust Matting等，这里还有个网站专门评测每年新提出的matting算法，这里有篇专门介绍matting算法的2007年的survey (by Jue Wang and Michael Cohen)。近年来的很多新算法其实大多是基于Closed-Form Matting和Robust Matting做improvement的，比如基于Closed-Form Matting的有Non-local Matting和KNN Matting，基于Robust Matting的有Shared-Matting，Global Sampling Matting，Weighted Color and Texture Matting，以及各种带有sampling字样的算法。这里主要介绍一下Closed-Form Matting。

Levin的Closed-Form Matting [4] 主要基于一个非常有见地的假设 -- Local Linear Model，大概意思就是，在一个很小的local window里，alpha matte应当可以用图像的三个颜色通道线性表出(linear combination)，即下图所示：

这个模型的精妙之处在于：第一，将原有的composition equation中的三个未知变量alpha，F，B解耦合，变成了未知变量alpha，a，b，而三个未知变量之间不再有相乘的关系，也就是nonlinear model转化成了linear model；第二，由于a和b是在一个local window里不变的，而每个pixel都被好几个重叠的local window覆盖，这些互相重叠的local window使得用户的输入可以有效进行propagation。有了这个模型，可以直接对其进行Regularized Least Squares Minimization来求解alpha，即：

其中w_j是每个pixel j附近的local window，S是用户输入的scribbles，epsilon是

regularization parameter(该regularization term相当于prior：希望alpha matte较为smooth)。该optimization中包含三个未知变量(每个pixel：alpha, a, b)，但幸运的是，目标函数相对于每个未知变量都是quadratic的，因此可以推导出其closed-form的解(分别对alpha，a，b求偏导并令式子等于0)。仔细观察目标函数，发现其关于a和b比较简单：每个pixelj可以单独计算。如果先将alpha看做已知，把目标函数对a和b求偏导并令其等于0，可以很容易计算出每个pixelj处a和b的表达式(包含alpha)，然后再将其代入目标函数中，可以得到一个只包含一个未知变量alpha的目标函数(仍然是quadratic)，写成矩阵形式如下(注意，这时求alpha不像刚才求a和b那么简单，无法对每个pixel分别求出alpha，因为每个pixel的未知量alpha与其他pixel的未知量有关联，只能通过求解线性方程组才能求出alpha)：

其中alpha是所有pixel的alpha value组成的列向量，其中L是可以由输入图像计算出的已知矩阵，Levin称之为Matting Laplacian Matrix，与前面介绍的Poisson Equation中的Laplacian Matrix有着很一致的形式和功能(propagation user input)。假设图像有N个pixel，Matting Laplacian是一个N×N的矩阵，其中在(i,j)处的元素为：

其中w_k是local window，|w_k|是local window内pixel的个数，注意I_i是3×1的vector(即rgb三个channel的pixel值)，I_3表示3×3单位矩阵，其他几个符号如下表示：

可以看出，只有当i和j代表的pixel处于同一个local window时，元素(i,j)才不为0。将δ_ij单独提出来(其实只有对角线上的元素有这一项)，令

Levin称之为Matting Affinity，跟前面colorization中用到的affinity相似，只是这里的affinity值可能为负。如果将local window size为2×2(只是为了方便直观显示，事实上一般用奇数边长的window size)，以5×5的图像为例，Matting Laplacian可以表示成如下形式(注意：

只要两个相邻的pixel能放入一个2×2的window，其affinity就不为空，因此事实上每个pixel周围8邻域的pixel的affinity都不为空)：

注意其中形如的元素表示的是输入图像中pixel(3,3)和pixel(3,4)的affinity。由上可

见，Matting Laplacian的形式与colorization中的Laplacian Matrix形式基本一致(如果colorization也用8邻域，那么形式一模一样)。考虑进去目标函数的constraint(即input scribbles)，可以使用跟上面colorization类似的mask列出方程组，或者使用Levin的TPAMI文章里介绍的方法。这里有Levin提供的Matlab代码。需要注意的是：这里的对角线元素不一定非负(non-negative)，非对角线元素也不一定非正(non-positive)，但是整个矩阵仍然是半正定(positive semidefinite)(Levin的paper [4]里有证明)，当然也是对称阵。(注意有些为前面的Laplacian Matrix特别定制的解线性方程组的算法不一定适用于这里的Matting Laplacian Matrix，比如[13]，针对的是非对角元素必须non-positive)。

另外值得一提的是，使用和这里同样的Local Linear Model，假设alpha matte是一个给定的guidance图像，将目标改成要计算出从原始图像到guidance图像的local linear transformation(即求出a和b的值)，将global optimization变为local window内的optimization，可以导出一个计算起来非常快的edge-preserving filter，这就是guided filter [12]。该filter由于其计算速度快，可以作为经典的bilateral filter的近似，很受欢迎。

2.3 Solvers: Solving Large Sparse Linear System

解大型稀疏线性方程组是一个很well-studied的问题，一般来讲，解法分为两大类：直接法(direct solver)和迭代法(iterative solver)。直接解法一般指高斯消元法(Gaussian

Elimination)，比如对于普通矩阵用LU Decomposition，对于对称正定阵用Cholesky Decomposition。这类算法的问题是，计算的过程中稀疏矩阵会变成稠密矩阵，如果矩阵规模较大(对于图像处理来讲很正常，比如对于1000×1000的图像，矩阵规模达到

1000000×1000000)，存储计算过程中的稠密矩阵都是问题，这就是文献中经常提到的fill-in问题[14]。这类方法的另一个问题是不太容易并行化。当然，也有很多新的direct算法克服了这些问题，参见Tim Davis的大作Direct Methods for Sparse Linear Systems [15](貌似Matlab里的backslash用的就是direct算法)。不过在computer vision中较受欢迎的还是iterative算法，推荐一本很棒的教科书，就是Yousef Saad的Iterative Methods for Sparse Linear Systems [16]。在ICCV 2007 gradient course中，有一个很棒的总结(section3 ppt)。这里也有个不错的总结。这里有个很棒的介绍Jacobi/Gauss-Seidel/SOR的ppt。

先将上述涉及到的方程组的系数矩阵及可以使用的算法进行一个分类，注意有些算法是对所有矩阵通用的，有些算法只能用于一些特定矩阵。以下分类基本是从上往下呈包含关系： 1. 普通矩阵(general)：LU Decomposition，Multigrid method；

2. 对角优势阵(diagonal dominant)：Jacobi/Gauss-Seidel/SOR(在diagonal dominant的情况下才收敛)；

2. 对称正定阵(symmetric positive-definite, SPD)：Cholesky Decomposition，Conjugate Gradient，Hierarchical Preconditioning (ABF) [14][17]；

3. M-matrix(非对角元素是non-positive的Laplacian Matrix，包括上述的Poisson Equation，colorization以及WLS filter中的Laplacian Matrix，但是不包括matting中的矩阵)：Hierarchical Preconditioning (HSC)[13]；

4. Spatially homogeneous Laplacian Matrix in Poisson Equation(在[13]中被如此称呼，注意跟homogeneous Poisson Equation不同)：FFT-based Poisson Solver (or see Chapter 2.2.6, Saad [16])；

这些涉及到的算法的复杂度如下(参见这里)(N是未知元素的个数，即图像像素个数)：

算法 LU Decomposition时间复杂度 N^3 适用矩阵类型 General Jacobi/Gauss-SeidelSORN^2 N^(3/2) Diagonal dominant Diagonal dominant Symmetric positive-definite Symmetric positive-definite M-matrix (off-diagonal non-positive) General Poisson Equation Conjugate GradientABF[17]N^(3/2) N^(3/2) HSC[13]N Multigrid methodN FFT-based methodNlogN 由上面可以看出，一般想要速度快的solver都会用Multigrid method，只是对于irregular的矩阵，quality不是很好。另外由于computer vision中的problem一般都是对称正定阵(甚至大多都是M-matrix)，Preconditioned Conjugate Gradient和SOR也很受欢迎。对于M-matrix，今年SIGGRAPH出现的HSC[13]算法也值得一试。

2.4 Reference

[1] R. Fattal, D. Lischinski, and M. Werman, \Gradient Domain High Dynamic Range Compression,\

[2] P. Perez, M. Gangnet, and A. Blake, \Poisson Image Editing,\[3] A. Levin, D. Lischinski, and Y. Weiss, “Colorization Using Optimization,” SIGGRAPH, 2004.

[4] A. Levin, D. Lischinski, and Y. Weiss, “A Closed-Form Solution to Natural Image Matting,” CVPR, 2006. ->TPAMI2008 extension.

[5] Z. Farbman, R. Fattal, D. Lischinski, and R. Szeliski, “Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation,” SIGGRAPH, 2008. [6] P. Burt and E. Adelson, \A Multiresolution Spline with Application to Image Mosaics,\TOG, 1983.

[7] P. Bhat, C.L. Zitnick, M. Cohen, and B. Curless, \GradientShop: A Gradient-Domain Optimization Framework for Image and Video Filtering,\

[8] D. Lischinski, Z. Farbman, M. Uyttendaelem, and R. Szeliski, “Interactive Local Adjustment of Tonal Values,” SIGGRAPH, 2006.

[9] A. Orzan, A. Bousseau, H. Winnem?ller, P. Barla, J. Thollot, and D. Salesin, \Diffusion Curves: A Vector Representation for Smooth-Shaded Images,\[10] J. McCann and N. Pollard, \Real-Time Gradient-Domain Painting,\2008.

[11] S. Paris, S.W. Hasinoff, J. Kautz, \Local Laplacian Filters: Edge-aware Image Processing with a Laplacian Pyramid,\

[12] K. He, J. Sun, and X. Tang, \Guided Image Filtering,\TPAMI 2012 extension.

[13] D. Krishnan, R. Fattal, and R. Szeliski, \Efficient Preconditioning of Laplacian Matrices for Computer Graphics,\

[14] R. Szeliski, \Locally Adapted Hierarchical Basis Preconditioning,\[15] T. Davis, \Direct Methods for Sparse Linear Systems,\

[16] Y. Saad, \Iterative Methods for Sparse Linear Systems,\[17] D. Krishnan and R. Szeliski, \Multigrid and Multilevel Preconditioners for Computational Photography,\

图像处理中的全局优化技术(Global optimization techniques in image processing and computer vision) (三)

2013-09-25 11:22 1083人阅读 评论(1) 收藏 举报

算法图像处理计算机视觉imagevision

MulinB按：最近打算好好学习一下几种图像处理和计算机视觉中常用的 global optimization (或 energy minimization) 方法，这里总结一下学习心得。分为以下几篇： 1. Discrete Optimization: Graph Cuts and Belief Propagation 2. Quadratic Optimization : Poisson Equation and Laplacian Matrix 3. Variational Methods for Optical Flow Estimation (本篇) 4. TODO: Likelihood Maximization (e.g., Blind Deconvolution)

3. Variational Methods for Optical Flow Estimation

Optical Flow(光流法)这个词乍一看很能唬住人，啥东东这么高级，是追踪光的流动轨迹吗？没这么玄乎。其实optical flow是一个很朴实的low-level vision的东西，就是每个pixel从一帧图像到另一帧图像的位置偏移(displacement)。如下图所示，

Two Input Frames Optical Flow (Vector Plot) Optical Flow (Color Plot)

上面的color plot，其实是通过一个二维的color wheel将2D的motion vector用颜色show出来，常用的color wheel如下所示(中心点表示横向和纵向的位移都为0，用白色表示)：

在上面的例子中，可以看出背景中大多数pixel是往左上方偏移一点(相对于镜头)，因此在optical flow中，背景呈现浅蓝色(在color wheel的第二象限)。至于计算optical flow这个东东到底有啥用途，可以说绝对是视频编辑的基石，参见这里(Art of Optical Flow，被墙的可以参看这里的英文原版)有一篇有趣的介绍optical flow在电影编辑中的作用(比如合成《黑客帝国》中的超级慢镜头)。

正是由于其重要的作用，如何计算optical flow从1980s就开始成为计算机视觉的研究热门。早期的计算主要focus在计算subpixel级的displacement，随着视频分辨率的增加，最近的很多算法开始考虑较大的displacement，通常需要先计算帧与帧之间pixel级别的correspondence，然后进行warping后再使用传统算法计算subpixel精度的optical flow结果。至于评测optical flow算法的accuracy，最经典的一个benchmark是Middlebury，但其用于排名的dataset只有12张图，比较容易overfitting，最近两年又出现两个比较popular的benchmark，Sintel和KITTI。

3.1 Classic Framework

我们先从optical flow算法的目的说起。令I(x,y,t)表示在t时刻frame上像素(x,y)处的亮度(或颜色)，那么optical flow的目的就是求出在t+1时刻的frame上，该像素相对于原来(x,y)的位移量(u,v)，用方程表示即：

其中u和v是未知量。这就是optical flow中著名的Brightness Constancy Model。

不过，数字图像毕竟是离散化pixel by pixel的，如果只给出两帧图像，如何计算出每个pixel的subpixel级的displacement(位移)呢？如果从correspondence的角度去考虑，在frame1中的某个pixel，只能找到其在frame2中对应的pixel整数位置，这样只能得到pixel级的displacement，而无法精确到subpixel精度。在经典的optical flow算法中，一般利用一阶泰勒展开作为工具来建立图像gradient和displacement之间的关系，这一步骤通常被称为Linearization。具体原理如下(这里有个不错的tutorial介绍的也很详细，by David Fleet and Yair Weiss)：

假设图像的intensity是连续的，如下图所示(1D case)，黑色的曲线表示frame1中的图像，蓝色的曲线表示frame2中的图像，待求的displacement是深蓝色的箭头dx/dt。对曲线进行一阶泰勒展开，其实就是假设曲线的局部是线性的，这样可以考察如图的红绿蓝组成的三角形。注意，dI/dx并非是红色线段的长度，而是其斜率。这样可以得到图中所示的关系，注意负号是因为斜率其实表示的是钝角的tan。这样一来就建立了图像的derivative和displacement之间的关系，注意dI/dx是图像在spatial上的derivative，dI/dt是图像在temporal上的derivative。

将上面的这个方程用在2D的图像上，对于每个pixel，可以写出以下方程：

其中，I_x和I_y是图像沿spatial的x和y两个方向的derivative(即图像gradient的两个分量)，I_t是图像沿temporal的derivative(可以用两帧图像的difference来近似)，u和v是该pixel沿x和y方向的位移，也就是待求的optical flow未知量。这其实就是对上面的Brightness Constancy Model的Lineariazation的结果，也就是optical flow中著名的Gradient Constraint Equation。

需要注意的是，上面这个模型是建立在两个假设基础上：第一，图像变化是连续的；第二，位移不是很大。如果这两点假设不成立，那么泰勒展开的近似是很差的。直觉想象一下，如果图像不是连续的，而且displacement很大，无论如何是无法将displacement和图像gradient联系起来的。不过在实际中，上述两个假设可以很容易使其成立。关于第一点，一般可以预先对图像进行Gaussian Smoothing，使其变化较为平缓；关于第二点，一般是对图像降分辨率建立金字塔，通过Coarse-To-Fine的方式一步一步去求解。

基于上面的方程，最经典的两个optical flow算法，Lucas-Kanade方法[1]和Horn-Schunck方法[2]，分别从不同的角度增加了求解该方程的稳定性。Lucas-Kanade方法是将每个pixel周围的一些pixels考虑进来，但每个pixel的未知量是单独求解的，所以是一种local方法；而Horn-Schunck方法是将上面的方程纳入到一个regularization的framework中，加入optical flow是locally smooth的prior，所有的pixel的未知量之间相互依赖，需要用global optimization的方法求解。目前state-of-the-art的方法，一般都是基于Horn-Schunck的framework做的，这里要介绍的经典optical flow算法也是基于这个framework经过几次改良得到的。

石器时代：令小写的u和v分别表示每个pixel处沿x和y方向的displacement，大写的U和V分别表示由所有pixel的u和v组成的map，那么Horn-Schunck的目标函数可以写作如下形式：

可以看出，第一项data term其实就是Brightness Constancy Model的Least Square形式，第二项regularization term目的是令结果的U和V较为smooth。对该目标函数进行minimization即可求出U和V。由于该目标函数是quadratic，根据calculus of variation(变分法)中的 Euler-Lagrange Equation (two unknown function -- U and V, two variables -- x and y) 可以将其转化为偏微分方程形式，再进行离散化，最终可以归结为求解一个大型线性方程组，利用Conjugate Gradient或者SOR可以很容易求解，参见上一篇文章。

青铜时代：根据robust statistics理论，quadratic的目标函数对outlier太敏感。而其实在locally smooth assumption下，最理想的情况是一个场景中只有一种motion。如果有多种motion，那么就相互构成outlier，会严重影响结果，这显然是quadratic目标函数的软肋。而这种情况恰恰是现实中最常发生的。因此，Michael J. Black和P. Anandan[3] 提出了将上述目标函数中的L2 norm换成更为robust的函数(比如L1 norm或者truncated L1 norm)，形成了目前较为常用的形式：

目前最常用的robust函数是Charbonnier函数psi(x^2)=sqrt(x^2+1e-6)，即L1 norm (增加了一个很小的数字1e-6为了使其convex，容易求解optimization)。这样一来，data term是L1 norm，regularization term使用了L1 norm后其实就成了Total Variation，因此这个目标函数又被称为TV-L1 formulation。另外，目前流行的coarse-to-fine、warping的framework也是在这篇paper里基本定形的。

白银时代：为了避免光照对Brightness Constancy Model的影响，Thomas Brox等人[4] 提出在data term里引入Gradient Constancy Model(不妨看做在原来的3个color channel之外多加两个gradient channel)，大大提高了算法的精确度。值得一提的是，这篇paper里提到的fixed point iteration解法(这里有Ce Liu的C++ code)比之前的graduated non-convexity算法(这里有Deqing Sun的Matlab code)貌似速度要快不少(或许是implementation的区别)。

另外值得一提的是，在Middlebury上有个叫improved TV-L1[6] 的算法，速度貌似很快，效果也不错。其实之前的很多算法都是基于TV-L1 formulation的，不过这篇paper [6] 主要是介绍一种快速算法的。其实，该文的主要contribution在于：第一，该文提出用texture decomposition的方式先提取image的texture，然后用texture作为data term，旨在避免光照等的影响，跟前面的加入gradient data term思路类似。第二，该文提出对每一个iteration求解出的flow进行median filter去掉outlier，后来这个小trick被证实是很有效的一步 [7]。另外，该文 [6] 中对implementation的描述也比较细致，如果要实现optical flow算法，这篇paper很值得一看(可以直接忽略其之前的版本[5])。

黄金时代：对经典算法中使用的各种trick进行总结的集大成者是Deqing Sun等人的CVPR 2010 paper [7]及其扩展IJCV 2013 paper [8]。该猛人将之前paper中五花八门的trick统统试了一遍，总结出哪些是精华，哪些是糟粕。最难能可贵的是，他将所有的matlab代码公布于众，使得后来者省去了不少揣测如何做implementation的时间，强赞一个。美中不足的是，这些实验是在Middlebury的training data上完成的，但是Middlebury的dataset有

点小，只有几张图，难免让人觉得有可能overfitting，调查的trick是否真的如paper中得出的结论那般，值得商榷(在其IJCV 2013 paper [8]中也发现，算法在KITTI的benchmark上和在Middlebury上performance很不一样)。

3.2 Algorithm and Implementation Details

这里简单地对求解上述framework的算法理一理，主要参考Thomas Brox的paper [4] 以及Ce Liu的Implementation及他对Brox算法的重新推导[9](Appendix A)。为了简化推导过程，下面假设图像I是由3个color channel加上2个gradient channel组成的5个channel的图像(paper [4]中推导时将gradient channel显式的表示了出来，看起来过于复杂)。由上面的TV-L1 formulation，经过Linearization后，目标函数如下所示(下标表示偏导)：

根据calculus of variation(变分法)中的 Euler-Lagrange Equation (two unknown function -- U and V, two variables -- x and y) ，可以得到以下两个偏微分方程：

这个方程有点麻烦，主要是因为robust函数psi引入的non-linearity。我们不妨先考虑没有psi时的情况，即Horn-Schunck的quadratic目标函数，这时上述偏微分方程很简单：

用divergence运算符及gradient运算符表示，即为：

注意这两个方程是对于每个pixel建立的，其中u和v是未知量，如果将divergence和gradient离散化(见上一篇中的基础介绍)，可以表示成当前点未知量uv与周围几个pixel未知量uv之间的线性关系。这样以来，相当于每个pixel可以列出2个线性方程，假设图像一共有N个pixel，那么一共可以建立2N个线性方程，未知量也是2N个，由于每个pixel与周围一些pixel存在依赖关系，这2N个方程组成的方程组需要联立求解。于是归结为求解一个大型稀疏线性方程组的问题(参见上一篇文章)，较为简单，可以用Conjugate Gradient或SOR求解。

回到上面那个复杂点的nonlinear偏微分方程，一个简单而有效的方法来消除non-linearity是，将nonlinear factor看做已知，用上一次迭代的未知量直接代入计算！这样多迭代几次，如果算法可以收敛，最终求出的未知量是可以满足整个nonlinear方程的！这就是fixed point iteration算法，一个可以化复杂为简单的有效算法。注意，每一次迭代其实都需要求解上述的线性方程组，在实际应用中，一般固定迭代次数(5~10次)即可。如果用大写Psi'_D表示上面方程中data term中的nonlinear factor，用大写Psi'_S表示smooth term中的nonlinear factor，那么复杂的方程可以简化为如下形式：

由于Psi'_D和Psi'_S可以直接使用上一次迭代的结果进行计算，该方程和上面的linear system解法一样。

更进一步，由于Linearization在u和v比较小时才比较精确，为了提高精度，在上面的迭代中，可以用上一个迭代的uv结果对第二帧图像进行warp，用warp之后的图像代替原来的第二帧图像，这样每次求的未知量uv其实只是增量，一般用du和dv表示，每次迭代将其累加到uv的结果上即可得到最终结果。具体来讲，上面的方程中涉及到第二帧图像的主要是I_t，将其用warp之后的第二帧图像计算，即作以下替换：

注意由于u和v是float point类型，上面的warp过程中一般需要进行插值计算，可以用bilinear或者bicubic插值算法。另外要注意的是，Psi'_D中也包含I_t项，如果也将其替换，那么需要作以下调整(不需要再考虑u和v)：

这样一来，方程就变为：

其中u和v用上一次iteration的结果，只有du和dv是未知量。注意，divergence和gradient operator离散化后其实是线性运算符，所以sooth term可以进行分离，整理后可以得到：

注意这两个方程是对于每个pixel建立的，等式左侧的du和dv是未知量，右侧是已知量(可以用上一个迭代的结果计算)。等式左侧含有divergence和gradient operator的项其实耦合了每个pixel与周围pixel的关系，这是regularization的功能所在。像上面介绍的一样，将该方程离散化后，可以建立2N个线性方程，联立可以求得每个pixel的du和dv。

另外，这个求解过程一般从图像金字塔上coarse-to-fine地进行，从最coarse的level开始，每个level计算完成再迁移到下一个level，使用上一个level的uv结果进行初始化。总之，整个算法可以简单的描述为以下流程(具体参见Ce Liu的C++ code)：

3.3 Large Displacement Optical Flow

在Thomas Brox和牛魔王Jitendra Malik近来的一篇TPAMI paper [10]里，提出将descriptor matching融合到optical flow的framework中，增加其对一些小的object的motion捕捉(这种小的object在之前的coarse-to-fine的方法中一般会丢失)。主要思想是增加一对辅助未知量u_d,v_d，将descriptor matching的cost和之前的cost function关联起来，这时的variational framework目标函数变为：

本文来源：https://www.bwwdw.com/article/873f.html