图像处理中的全局优化技术(经典至极).docx

图像处理中的全局优化技术(Global optimization techniques in image processing and computer vision) (一)2013-05-29 14:261659人阅读评论(1)收藏举报算法图像处理计算机视觉imagevisionMulinB按：最近打算好好学习一下几种图像处理和计算机视觉中常用的 global optimization (或 energy minimization) 方法，这里总结一下学习心得。分为以下几篇：1. Discrete Optimization: Graph Cuts and Belief Propagation (本篇)2. Quadratic Optimization : Poisson Equation and Laplacian Matrix3. Variational Methods for Optical Flow Estimation4.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的paper2，以及这篇paper提到的一个max-flow算法3(该max-flow算法最早是发在2001年的一个CVPR Workshop上，后来扩展到TPAMI 3)。需要注意的是，这两篇paper里的Graph Cuts算法，只是针对只有两个label (待求变量是binary variable)的情况。而Boykov的2001 TPAMI paper4提出使用alpha expansion技术将多个label的问题转化成一系列的binary label问题去逼近，这使得Graph Cuts算法开始风靡起来。后来Kolmogorov的2004 TPAMI paper5进一步讨论了什么样的能量方程可以被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算法的步骤很简单，如下：cppview plaincopy1. /alphaexpansionalgorithmpseudo-code2. initializelabeling;3. 4. whilenotconverged5. 6. foreachlabelainL7. 8. constructagraph;9. domax-flowcut;10. ifenergyissmallerthanbefore,acceptit;11. elsedeclineit;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) C+B)4. 5. Valuedelta=A+D-C-B;6. ValuesubtrA=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 paper7。1.2 Belief PropagationBelief Propagation是Graph Cuts的一个强劲对手，其渊源可以追溯到MIT在ICCV 2003上的paper8，比较这两种方法在stereo上的性能，结果貌似不分上下，Graph Cuts略胜一点点。后来大牛Richard Szeliski联合一堆MRF的experts搞了一个benchmark评测这些方法，见TPAMI 2008 paper9，和这个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 paper10，有代码在此).下面简单介绍一下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 paper12。在这种情况下，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算法非常简单，流程如下：cppview plaincopy1. /BeliefPropagationAlgorithmpseudo-code,supposerunTiterations2. for(intt=0;tT;t+)3. 4. updatemessagebetweeneachneighboringpixels;5. 6. 7. selectlabelforeachpixelsuchthatthesumofdatacostandmessagespassingtoitisminimum;最主要的步骤在于每次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 paper10里，针对一些常见的smoothness cost介绍了一些快速算法可以将复杂度减少到O(L)的计算复杂度，基本涵盖了上面表格里列出的cost类型。这样的message update经过T次迭代后，最后通过如下公式计算每个pixel的cost，选择cost最小的label即为最终结果。上面的message passing方式其实比较慢，一般需要很多次iteration才能使得message传递到全图。在Pedro Felzenszwalb的IJCV 2006 paper10里，提出使用multiscale/hierarchical的方式使得message更加快速的进行传递，这使得需要的iteration数量大大减小，使得BP的速度可以很快。1.3 Other Related TechniquesGraph Cuts是从处理Interactive Image Segmentation (或叫Seeded Image Segmentation) 问题开始发家的 2，其他几种可以处理这种问题的方法有：Random Walker, Geodesic Shortest Path, Level Set Method等。这里简单介绍一下。首先值得一提的是，SIGGRAPH 2004井喷了三篇基于Graph Cuts做interactive image editing的文章：GrabCut13, Photomontage14, Lazy Snapping15。其中以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 paper16，注意，该算法在有些文献里叫做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 map17的问题中，这也是后话。Leo Grady后来还将Random Walker和Graph Cuts联系了起来 (ICCV 2007 paper18)，其实从上面的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 paper19。这种算法有个优点，就是所有pixel距离某一种seeds的geodesic distance可以非常快速的计算出来，用N表示pixel的个数，那么用Yatziv的Fast Marching20(在19中使用的算法，是Sethian的Fast Marching Method的快速实现)，或者更适用于并行计算的Parallel Marching21(在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 Reference1 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 CutsA 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, 200810 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, “GrabCutInteractive 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. Natl 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:361556人阅读评论(0)收藏举报算法图像处理计算机视觉imagevisionMulinB按：最近打算好好学习一下几种图像处理和计算机视觉中常用的 global optimization (或 energy minimization) 方法，这里总结一下学习心得。分为以下几篇：1. Discrete Optimization: Graph Cuts and Belief Propagation2. Quadratic Optimization: Poisson Equation and Laplacian Matrix(本篇)3. Variational Methods for Optical Flow Estimation4.TODO: Likelihood Maximization(e.g., Blind Deconvolution)2. Quadratic Optimization: Poisson Equation and Laplacian MatrixQuadratic Optimization (Least Squares Minimization)在图像处理中的魅力要从SIGGRAPH 02和03年的两篇Gradient Domain Image Editing文章说起：Fattal的HDR Compression1和Perez的Poisson Image Editing2。 其后，Levin的两篇文章Colorization3和Closed-Form Matting4更是将其魅力展现的淋漓尽致。而Farbman的基于Weighted Least Squares的WLS filter5也是在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 Pyramid6图像融合(这里有个简洁的中文介绍)，这是题外话。在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形式。以55的图像为例，假设待求图像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的推导过程。如果进一步对上述