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
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]的推导过程。
如果进一步对上述的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。
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)处的元素为：
可以看出，只有当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都不为空)：
注意其中形如
另外值得一提的是，使用和这里同样的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-Seidel | N^2 | Diagonal dominant |
| SOR | N^(3/2) | Diagonal dominant |
| Conjugate Gradient | N^(3/2) | Symmetric positive-definite |
| ABF[17] | N^(3/2) | Symmetric positive-definite |
| HSC[13] | N | M-matrix (off-diagonal non-positive) |
| Multigrid method | N | General |
| FFT-based method | NlogN | Poisson Equation |
2.4 Reference
[2] P. Perez, M. Gangnet, and A. Blake, "Poisson Image Editing," SIGGRAPH, 2003.
[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," TOG, 2009.
[9] A. Orzan, A. Bousseau, H. Winnemöller, P. Barla, J. Thollot, and D. Salesin, "Diffusion Curves: A Vector Representation for Smooth-Shaded Images," SIGGRAPH, 2008.
[10] J. McCann and N. Pollard, "Real-Time Gradient-Domain Painting," SIGGRAPH, 2008.
[12] K. He, J. Sun, and X. Tang, "Guided Image Filtering," ECCV, 2010. ->TPAMI 2012 extension.
[13] D. Krishnan, R. Fattal, and R. Szeliski, "Efficient Preconditioning of Laplacian Matrices for Computer Graphics," SIGGRAPH, 2013.
[14] R. Szeliski, "Locally Adapted Hierarchical Basis Preconditioning," SIGGRAPH, 2006.
[15] T. Davis, "Direct Methods for Sparse Linear Systems," 2006.
[16] Y. Saad, "Iterative Methods for Sparse Linear Systems," Second Edition, 2003.
[17] D. Krishnan and R. Szeliski, "Multigrid and Multilevel Preconditioners for Computational Photography," SIGGRAPH Asia, 2011.