计算机视觉:泊松融合

本周给大家带来的算法是:Poisson Blending!这是图像处理界大名鼎鼎的图像融合算法,自从2003年被发表[1]以后,有很多在此基础上进行改进的研究。

图像融合,就是把不同的图像的不同部分放在一起,形成一张新的图像。融合得越自然,算法就越好。想到融合,最简单的算法就是在融合边界把两张图像的透明度线性相加,形成一张新的图像,比如这幅著名的苹果橘子图:

apple_orange

这种算法叫做alpha blending,它的进阶版本叫做laplacian blending,也就是利用我们上次所说的高斯金字塔来构建拉普拉斯金字塔,在每一层都进行一次alpha blending,从而达到各个频率完全融合的效果。感兴趣可以参考[2]

上面的算法有一个缺点:当两幅图片颜色并非十分相近的时候,融合出的结果并不理想。

想要把左边的图融合成右边的样子,我们就要借助新的工具:泊松融合。首先,依然是目录~

  • Poisson Blending0:预备知识(图像的梯度、泊松方程)
  • Poisson Blending1:泊松方程应用于图像
  • Poisson Blending2:解泊松方程
  • Poisson Blending3:混合柏松融合
  • Poisson Blending4:更多用途

    Poisson Blending0:预备知识(图像的梯度、泊松方程)

    进入正题之前,我们先补充一下基础知识。

      图像的梯度

      什么是图像的梯度?我们可以把图像看成是一个离散的函数,函数在某一点的梯度是它的导数,那么根据导数的定义,我们可以把图像在某一点的梯度定义为:下一个点的像素值-这个点的像素值,即g(i)=f(i+1)-f(i)。至于下一个像素点在那个位置,要根据方向来定。例如行梯度就是右边的相邻点(g(x,y)=f(x+1,y)-f(x,y)),列梯度就是下方的相邻点(g(x,y)=f(x,y+1)-f(x,y))。那么图像的梯度表示什么意思呢?很明显,就是相邻两个像素的差异大小。为了直观表示,我们可以试着表示出一个图像的梯度:

      原图; x方向的梯度; y方向的梯度
      虽然有些费力,但你应该能看出,梯度图像其实是把边缘检测了出来。当然,检测边缘有更加复杂的filter。</p>
      泊松方程

      柏松方程的表示为 equation

      它在物理上有很多应用,例如表示电荷分布与电势的关系,或者重力加速度和密度的关系等等。当然,我们学习泊松融合并不需要了解这些知识,感兴趣的话,你需要学习大学物理(B)第二册:)。图像处理上,则是借鉴了泊松方程以及其他物理分析阐明的函数f、f的梯度(一阶算子)以及f的二阶梯度(拉普拉斯算子)之间的关系。你同样不需要了解这些名词的意义,不过有一点要知道:梯度表示了函数变化的趋势。

      为什么要知道这一点?因为泊松融合的思想并不是让两幅图叠加,而是让目标图像在融合部分“生长”出源图像。也就是说,只提供原图像的斜率,让目标图像根据自己图像的特点,按照对应斜率生成融合部分。由于是按照自己的特点,并没有添加外来的元素,生成的图片就更加自然。

      打个比方来说:如果要融合两座山峰(A山表示为1 2 3 2 2 1, B山表示为12 14 13 14 15 14, 数值表示高度),直接对半拼起来肯定会不自然(取A左边+B右边:1 2 3 14 15 14);如果把其中B山峰生长的趋势记下来(表示为梯度:2 -1 1 1 -1),然后在另一座山峰的融合处上按照此趋势“生长”出另一座山峰,就显得更加自然(结果为C山:1 2 3 4 5 4)。

    Poisson Blending1:泊松方程应用于图像

    在泊松融合中,我们真正要处理的是这样一个函数

    equation

    别犯晕!我来解释一下是什么意思:如果要把图像B融合在图像A上,那么f表示融合的结果图像C(也就是我们要求的图像)(相当于C山),f*表示目标图像A(相当于A山),v表示源图像B的梯度(相当于B山的梯度),▽f表示f的一阶梯度(也就是图像的梯度),Ω表示要融合的区域,由于那个卷曲的字母我打不出来,就用d代替,dΩ代表融合区域的边缘部分。那么这个式子可以解释为:在目标图像A边缘不变的情况下,求融合部分的图像C,使得C在融合部分的梯度与源图像B在融合部分的梯度最为接近(也就是生长趋势最接近)。

    这里有人可能会有疑问:如果梯度最接近的话,直接复制过去不就好了么?我们要注意到,这里有一个条件:C的融合边缘与A的融合边缘相等。如果直接把B复制过去,内部的梯度是与B相等了,但是边缘就会有问题,这个结果往往并不是上式的最小值。

    所以我们从边缘出发,要求得一个生成图,满足边缘不会突变,而且融合部分的变化趋势几乎符合源图像的变化趋势,这样的结果就是我们想要的。

    如果你仍有困惑,可以尝试理解一下一维上的解释,假设我有一个函数f,它在一部分区域上的斜率如下图(1);另一个函数g,它的斜率如下图(2)。我们如何更自然地融合两个函数?如果斜率有突变,就会使得融合的结果不自然(下图(3)),因此我们应该尽量使斜率连续,这里就需要在保持f在区域边界上不变的情况下,生成区域中的斜率使之最为接近g的对应区域中的斜率(下图(4))。

    目标区域f函数的斜率g函数的斜率
    不自然的融合(直接拷贝)自然的融合(泊松融合)

    始终注意,我们这里讨论的是融合二者的斜率(高维上叫 梯度),而不是二者的函数值。斜率连续变化,函数就不会有突兀的改变。


    Poisson Blending2:解泊松方程

    在弄清楚整个算法之后,我们就可以开始求最优函数了。但如何求一个函数的最小值呢?这里,我们继续考虑函数的梯度问题。

    对于一个连续函数f,它的最大/小值点的梯度值一定为0——这应该是众所周知。我们就以此作为切入点,用梯度为0来作为最优函数的条件。

    这里你可能会问,如果是求到最大值怎么办?对于这个函数来讲,由于距离可以无限远,但最小只能为0,因此它拥有最小值,而没有最大值,所以梯度为0的点必定是最小值点。

    这个函数梯度为0的表示可以转化为我们刚刚讲到的泊松方程:

    equation

    再解释一下这个等式:Δ表示拉普拉斯算子,也就是二阶梯度(一阶梯度▽f再求一次梯度啦)。div(Divergence,也就是▽)表示v的梯度,也就是▽B的梯度(B是源图像)。如果你又乱了,就把它理解为:源图像B的二次梯度与新图像C的二次梯度相同。

    好嘞,有了这个等式,我们就有了确定的方法来求解目标图像了。但是由于图像是离散的,我们需要把上面的等式变换一下,符合图像的特点。首先根据梯度的定义,我们先明确下二阶梯度的表示:

    equation

    一个点上图像梯度的简单定义:该像素与周围像素的差值的和(对于不在边缘的像素点,是与四个周围像素的差值)。比较详细的一个推导图像二次梯度的文章可以看这里

    好,那么我们就可以把上方的等式转换为下面的差分方程(总结了每一个Ω内部像素点的所有情况):

    equation

    如果H表示新的图像,一个相邻像素是(1)边界像素,那么它的值是固定的(因为我们要保留边界值);(2)超出选区边界,被排除。N表示H(x,y)邻居像素点数目(< =4,如果此点位于图像边界那么他可能只有2个或3个邻居像素),而左边公式的最后一项则表示如果该像素点的邻居是选区边界的话,那么就直接用源图像A中对应点的值。右边表示目标图像B在此点的二阶梯度。

    于是我们可以开始解方程啦!上式中,与H有关的都是未知数,我们需要做的是把这些像素的值求出来。

    关键问题是,一个融合区域有那么多的未知数,他们互相关联,怎么求呢?这里如果我们把求二阶梯度(也就是式子左边)看成一个n *n(n表示要填充的像素总数)矩阵A,把未知区域的H看做是未知数x,把已知的源图像B的二阶梯度(也就是式子的右边)看做是一个(n*1)列矩阵B,那么问题的本质就是求解Ax=B的问题。

    矩阵形式的表示:

    matrix_equation

    这是优化领域内一个大的学问。很多适用于计算机的方法被提出。因为A巨大(未知数太多了),所以常规求A逆的计算复杂度简直是爆炸级别,所以人们提出了一些迭代方法来求解。最最经典也是最最基础的叫做牛顿迭代法[3],大家可以参考链接学习。在这里我们不讲这些方法的理论基础以及证明(那样我又要写最少十几页),如果感兴趣,可以参考 雅克比迭代法[4] 共轭梯度法[5] 以及十分全面且由浅入深的关于共轭梯度法的 讲解材料[6](虽然很长,但不痛苦,看完绝对收获颇丰)。

    当然,还有tons of better methods,但是由于太complicated,我们在这里使用Jacobi来入门,然后实际运用CG(共轭梯度)方法。

    以下就来详细讲一讲如何把Jacobi应用到这个问题上来。

    通过参照wikipedia上的讲解,我们看到矩阵A可以被分解为一个对角矩阵D以及对角为0矩阵R的和,然后每一步迭代的步骤是X(k+1)=D-1(b-RX(k))。我们可以设置X(0)=0,然后迭代足够多次(例如5000),得到的结果就已经比较满意。或者把停止条件设置成偏差值小于某个数也可以。

    针对特定问题,b代表源图像在融合区域的二次梯度的矩阵,我们可以先把目标图像的融合区域全设置成0(也就是黑色)。怎么表示RX(k)呢?我们注意到,如果按照上面那个最长公式,我们列出当填充区域是3*3矩阵的时候,A应该长着个样子:

    equation

    那么A在分解后,R就应该是这个样子:

    equation

    所以我们只需要对融合区域作用一个这个样子的filter来扫描一遍就可以的到RX(k)

    由于D-1是对角矩阵,而D的对角元素都是-4,因此我们只需要简简单单地让(b-RX(k))乘以-1/4就好,然后把结果赋值给下一个X(k+1),如此循环往复。

    这样下来,融合100*100的图像区域大约会消耗40-50秒的时间,当然,如果运用CG的方法,会达到20-30秒,更快一些。


    Poisson Blending3:混合泊松融合

    在融合的时候,如果目标图像是有纹理的物体,比如棵树干,或者一个橘子,我们希望融合后的图像也具有这样的纹理,应该怎么办呢?如果用我们上面的方法,最好的结果是这样:

    equation

    注意看眼睛周围的皮肤部分,并没有橘子皮一样的纹理。

    再比如如果我们把一幅图融合到一张布匹的图片上,想让融合出的结果具有布匹的纹理,这些都该如何做呢?

    实际上,我们只要对上面公式的参数做小小的改动,就可以得到新的结果。来想一下,纹理在图像的梯度层面表示了什么含义?纹理越多,对应的梯度图像中边缘也就越多,也就是说,图像中梯度不为0的地方也就越多。那么我们只需要对原式中的V做一个小小的改动,原先V表示为源图像的梯度图像,而现在我们把它表示为 源图像梯度图像 与 目标图像梯度图像 对应像素绝对值 的 最大值。也就是Vpq = fp∗−fq∗ if|fp∗−fq∗|>|gp−gq|,或者gp − gq if|fp∗−fq∗|<=|gp−gq|.

    这样,我们就保留了原图的纹理,又合成了新图的图案。

    不过我在做这个步骤的时候,加入了少许创新:我让两个图像在低频部分泊松融合,高频部分直接保留原图(因为高频代表细节、纹理而且不包括颜色信息)。这样的结果似乎效果不错,比如下左是一幅有做旧效果的蒙娜丽莎,注意它的纹理是很多不规则的小点;右边是要融合的比尔盖茨的脸:

    在运用了我的效果之后,融合成新的蒙娜盖茨,依然保留了纹理信息(图挂了,不好意思。。。)

    另外,这种方法也可以融合透明的物体,比如彩虹。

    Poisson Blending4:更多用途

    泊松融合还可以做什么?比如“纹理压平”技术(就是把3D的拍扁成2D),或者对局部颜色的调节。总之,如果你想更深入全面地了解这项技术,最好去参阅下这篇经典的论文[1]


    Reference

  • Poisson Image Editing
  • Image Pyramids and Blending
  • Jacobi method
  • Conjugate gradient method
  • Painless conjugate gradient
  • Poisson Image Editing Tutorial


  • © Mingrui Zhang. All rights reserved.