文献阅读 - Poisson Image Editing

泊松图像编辑是一种基于数学偏微分方程（PDE）的方法，用于在目标图像中无缝融合源图像。其核心思想是用源图像的色彩梯度代替目标图像的色彩强度进行融合。该方法通过求解泊松方程来生成梯度匹配的合成图像。具体而言：
核心思想：
泊松图像编辑通过比较源图和目标图的颜色梯度强度，在不改变边界区域颜色强度的前提下进行合成。
关键步骤：

• 确定源图和目标图的定义域及其边界条件。
• 使用拉普拉斯算子提取二阶变化量，并通过最小化二次优化问题求解目标函数。
• 通过离散化方法将连续问题转化为有限差分方程组，并求解得到最终结果。
  应用场景：
  泊松编辑广泛应用于图片修复、无缝克隆等场景中。例如，在平面替换或复杂背景去除时，该方法能够有效保留目标区域的真实感。
  算法特点：
  该方法不仅适用于灰度图处理（单色迁移），还支持多通道（如RGB）处理；同时结合不同克隆方式（常规、混合、单色迁移）以适应不同需求。
  数学基础：
  泊松方程描述了导数匹配约束下的最优问题，其解满足Dirichlet边界条件下的拉普拉斯方程形式。
  实现细节：
  使用高斯-塞德尔迭代法等数值方法求解离散化的线性系统，确保高效稳定地得到结果。
  效果评价：
  通过对合成效果的评估（如对比分析），可以验证该方法在保持边缘细节和颜色一致性的有效性。

1 泊松图像编辑（Poisson Image Editing ）

泊松图像编辑（Poisson Image Editing）的主要思想在于，在图像融合过程中，将图象的色彩亮度替换为其亮度变化率（即使用颜色梯度），从而能够获得更加自然真实的图像融合效果。

经无缝融合（seamless cloning）处理后，在合成图中被遮罩区域（mask region）内目标亮度值与原始输入图景（source image）间的亮度值存在明显差异；然而，在该区域内的亮度变化趋势（intensity gradient）基本吻合。值得注意的是，在该区域边缘处的目标亮度值（target brightness level）与合成图完全匹配。

三位研究者分别为Rez P、Gangnet M以及Blake A在其合作研究中提出了Poisson image editing这一技术，并将其发布于《图形学汇刊》（Acm Transactions on Graphics: a leading platform for graphics-related research）一文中。该技术的相关论文链接为：http://www.irisa.fr/vista/Papers/2003_siggraph_perez.pdf. 论文详细阐述了该方法的核心原理及其在图像编辑领域的应用效果，在《Acm Transactions on Graphics》期刊上于2003年发表于第22卷第3期（第564页）。

1.1 引言（Introduction）

该方法的数学思想是给定狄利克雷边界条件的泊松偏微分方程。

（1）图像中，拉普拉斯算子提取的二阶变化量在感知上最明显；

（2）单一解：在有界区域内的标量函数（见图中标注链接）可以通过其边界面条件以及内部区域内的拉普拉斯运算来精确确定（即该函数\varphi(x,y)能够通过边界条件f(x,y)及其内部区域的拉普拉斯运算来精确确定）。

泊松方程的求解过程等价于一个优化问题：给定边界条件，在基于L_2范数的意义下, 寻找其梯度最逼近指导向量场的那个函数。

1.2 导向插值求解泊松方程（Poisson solution to guided interpolation）

1.2.0 泊松方程（Poisson’s Equation）

• 泊松方程

- \Delta \varphi = f

• 狄利克雷边界条件（Dirichlet boundary condition）

微分方程的“第一类边界条件”，给定微分方程的解在边界处的值。

• 给定狄利克雷边界条件的泊松偏微分方程

\begin{cases} -\Delta\varphi(x,y)=f(x,y) & \text{当 }x,y\in\Omega \\ \varphi(x,y)=g(x,y) & \text{当 }x,y\in\partial{\Omega} \end{cases}

设Ω为有界开区域，则称Δ=∂²/∂x²+∂²/∂y²为拉普拉斯微分算子（Laplacian operator）

1.2.1 引导插值（Guided Interpolation）

• 引导插值：利用引导向量场（guidance vector field）进行图像插值。

S：R^2的闭子集，表示图像的定义域；

\Omega：S的闭子集，其边界为\partial \Omega；

f^{\ast}：其定义域是在\Omega域边界及其外部区域上的已知数量场（quantity field）；

f：待求解的无量纲标量场；该场限定在\Omega内部区域（即\Omega - \partial \Omega）。

\mathbf{v}：定义在\Omega域上的矢量场（vector field）

• 最小化问题：

We aim to find a function f that minimizes the integral \iint_{\Omega}\|\nabla f -\mathbf{v}\|^2 over the domain \Omega, subject to the boundary condition \left.f\right|_{\partial{\Omega}}=\left.f^{\ast}\right|_{\partial{\Omega}}.

（3）

在本问题中，在公式\nabla = \left[ \frac{\partial }{\partial x}, \frac{\partial }{\partial x} \right]中定义了梯度算子（gradient operator），该运算符用于描述矢量场\mathbf{v}的变化方向，在此语境下被定义为引导场（guidance field）。这个问题的解是满足给定狄利克雷边界条件的Poisson partial differential equation的唯一解：

\Delta f equals \nabla\cdot\mathbf{v} within \Omega, subject to the boundary condition that f on \partial\Omega is equal to f^* on its surface.

（4）

此外该方程组\nabla \cdot \mathbf{v} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}代表向量场\mathbf{v} = (u, v)的散度（divergence）。

对于三通道彩色图像来说，泊松方程需要在三个通道上分别独立求解。

当矢量场\mathbf{v}被视为保守场时（conservative），若满足\mathbf{v} = \nabla g的关系，则可以在区域\Omega内构造一个修正函数\tilde{f}。此时，在该区域上通过引入修正函数\tilde{f}到基函数f中进行调整后，则原问题中的泊松偏微分方程将被转化为拉普拉斯方程。

The Laplacian of \tilde{f} is zero within the domain \Omega, and it satisfies the boundary condition \tilde{f} = f^* - g on the boundary \partial\Omega.

由此可见，在区域Ω内部定义的附加校正\tilde{f}等于源图像与目标图像在边界\partial \Omega处的失配(f^{\ast} - g)所形成的膜插值表面。

1.2.2 离散泊松解算器（Discrete Poisson solver）

S：表示图像的一部分区域；

p：表示一个像素点；

N_p：表示与像素p相邻的四个像素的集合；

\left< p, q \right>：表示满足q \in N_p的像素对；

域\Omega的边界定义为：\partial \Omega等于所有属于S\setminus\Omega且满足N_p与\Omega有非空交集的所有点p

f_p表示f在p处的值；

v_{pq}表示\mathbf{v} \left( \frac{p + q}{2} \right)沿着有向边\left[ p, q \right]的投影为：

v_{pq} = \mathbf{v} \left( \frac{p + q}{2} \right) \cdot \vec{pq}

该任务为计算灰度强度的集合\left. f \right|_{\Omega} = \left\{ f_p, p \in \Omega \right\}

针对定义于任意边界的狄利克雷边界值问题, 通常建议采用有限差分法进行离散化处理, 将方程(3)转换为一个离散化的二次优化模型, 进而实现对其的有效求解

\text{minimum}_{f|\Omega} \sum_{\langle p,q\rangle:\Omega\cap\langle p,q\rangle\neq\emptyset} (f_p - f_q - v_{pq})^2, 其中 f 在 \Omega 上满足 f_p = f_p^{\ast} 对于所有位于 \Omega 边界的点 p

（6）

【公式推导】

方程（6）的推导过程中需要注意的是对N_p这一变量的定义域设定。在分析向量\vec{pq}的方向特性时发现其可能取值仅限于0、π/2、π、3π/2四个方向范围内，在此约束下该向量长度不会超过4个单位长度。

其解满足联立线性方程组（simultaneous linear equations）：

左边部分：\left|\,N_p\,\right|乘以f_p减去\sum\limits_{q\in N_p\cap\Omega}的f_q之和；右边部分：等于\sum\limits_{q\in N_p}\!\!\bigcap{\partial{\Omega}}}的f^*_q加上\sum\nolimits_{{}}v_{{pq}}}之总和，并且对于任意属于Ω中的元素p都成立

（7）

对于属于区域边界的点p，在其邻域内元素的数量满足\left| N_p \right| < 4；若p位于区域Ω的内部区域，则方程组(7)仅涉及内部节点。

\left| N_p \right| f_p - \sum_{q \in N_p \bigcap \Omega} f_q = \sum_{q \in N_p} v_{pq}

（8）

方程（7）具有带状稀疏性（banded sparse）、对称性（symmetric）、正定性（positive-definite）特征，并满足高斯-塞德尔迭代条件

1.3 无缝融合（seamless cloning）

1.3.1 导入梯度（Importing gradients）

源图像中的梯度场可以作为最基本的引导场\mathbf{v}：

\mathbf{v} = \nabla g

其中，g表示源图像。代入方程（4）可得：

The Laplacian of function f, denoted as \nabla^2 f, is equal to the Laplacian of function g within the domain \Omega, subject to the boundary condition where f equals f^* on the boundary \partial\Omega.

（10）

v_{pq} = g_p - g_q, \ \forall \left< p, q \right>

导入梯度在 \Omega域内完全不考虑目标图像。

1.3.2 混合梯度（Mixing gradients）

混合梯度会将 f^{\ast}与g的特性组合起来。

该研究采用Poisson方法通过非保守引导场实现更好的效果。在\Omega域内的任一点上 选取具有较大梯度的函数f^{\ast}和g中的一个作为引导向量\mathbf{v}的取值依据

\bm{v =}\begin{cases}\nabla f^*(\bm{x}) &\text{当 }|\nabla f^*(\bm{x})| > |\nabla g(\bm{x})|\\[4pt]\nabla g(\bm{x}) &\text{其他情况}\end{cases},\quad\forall\bm{x}\in\Omega

（12）

引导向量\mathbf{v}的离散化表示为：

该算法中的相似度计算采用以下规则：当两个特征值f_p^{\ast}与f_q^{\ast}之间的差值\left| f_p^{\ast} - f_q^{\ast} \right|显著超过目标函数g_p - g_q的变化范围时，则取其差值作为v_{pq}的估计值；否则，则采用目标函数的变化量作为v_{pq}的估计值。对于所有有序对\langle p, q\rangle均适用此计算规则

1.3.3 单色迁移（Monochrome transfer）

若无需求携带色彩信息，则可将源图像转为单色形式进行处理

1.4 选区编辑（Selection editing）

无间隙的图像融合问题涉及将目标区域从源图像中提取并嵌入到目标图像中。其中引导场的方向或全部由源图像g的梯度决定。

选区编辑处理的核心问题是基于单一图像实现特定区域的操作，并通过引导场建立该区域与原始图像的空间关系以实现变换效果。其主要操作形式包括：基于原始图像构建引导场以指导选区图像的去扁平化操作；通过调整该区域的局部亮度以实现亮度分布的一致性；利用颜色空间模型对目标区域进行色彩调校；以及通过无缝拼接技术保证整体视觉效果的一致性和连贯性。

纹理被扁平化以及局部呈现明暗变化是以对原始梯度场\nabla f^{\ast}进行非线性调整来完成的；背景与前景色调的调整以及边缘融合则是采用就地无缝整合的方式完成的。原始图像经过域内处理（提供新的基础图像）或域外处理（补充新的边界信息）后便被用作新的来源。

1.4.1 纹理扁平化（Texture flattening）

通过稀疏采样算法对多通道数字图像的梯度矩阵\nabla f^{\ast}进行处理后，筛选出最具代表性的特征。

该向量由矩阵\mathbf{M}\left(\mathbf{x}\right)与函数f^{*}在点\mathbf{x}处的梯度相乘得到。

（14）

其中所述的\mathbf{M(x)}被定义为一种二进制蒙版矩阵，在其应用过程中可实现对稀疏筛选区域的精确识别。在特定条件下，在当前分析框架下当所述的\mathbf{M(x)}被配置作为边缘检测器时。

v_{pq} is defined as the difference between f_p^{\ast} and f_q^{\ast} when an edge connects nodes p and q, otherwise it is zero. This holds for every pair of nodes \left< p, q \right>.

在该模型中，在边缘检测器所选的边缘数量越少（即选择能力越强）的情况下，在目标区域进行的映射点就越稀疏（即分布越分散），从而呈现出的扁平化程度就越高。

1.4.2 局部明亮变化（Local illumination changes）

Fattal et al. 2002，给出一种平滑修改原始图像动态范围的方法：

（1）对图像对数的梯度场进行变换，用于减小大的梯度值、增加小的梯度

在全局图像域内，在基于诺伊曼边界条件（Neumann boundary conditions）的约束下，通过求解泊松方程\Delta f = \nabla \cdot \mathbf{v}来实现利用变换后的向量场\mathbf{v}来重建图像的对数值f。

基于适当的狄利克雷边界条件，
该方法被用于校正选定区域的图像。
采用基于简化的Fattal et al变换方案，
引导场被定义于对数域（log-domain），其特性具有特定的数学属性

\mathbf{V} = \alpha^{\beta} \left| \nabla f^{\ast} \right|^{- \beta} \nabla f^{\ast}

其中，
\alpha 表示为 0.2 \left| \nabla f^{\ast} \right|_{\mathrm{avg}}（在 \Omega 域内 f^{\ast} 的平均梯度范数的 20%），而 \beta 则被设定为 0.2。

1.4.3 局部色彩变化（Local color changes）

基于输入的原始图像I及其选定区域Ω域，在此基础上, 通过无缝融合技术对原始图像进行前景色与背景色的修改。

（1）优化颜色设置：对于定义在\Omega域外的目标函数f^{\ast}而言，在图中中间位置示例：通过对选区外区域的去色处理操作后所获得的图像即为所求的目标函数f^{\ast}。

（2）调整前景色效果：在\Omega域内的源函数g中，请参考下图右侧部分：将选择区域内的图像各通道灰度值分别乘以1.5、0.5及0.5后所获得的图像数据即可定义为源函数g。

1.4.4 无缝拼接（seamless tiling）

选区（\Omega域）为矩形：

选取区域内的原始图像作为基准图g后，在其基础上构建边值条件

2 OpenCV实现

2.1 函数原型

seamlessClone(src, dst, mask, center, flags)

* 常规融合（Normal Cloning）

源图像的纹理（梯度）保留在融合区域

* 混合融合（Mixed Cloning）

融合区域的纹理（梯度）由源图像与目标图像共同构成。混合融合过程会倾向于选择源图像与目标图像中的主要纹理（梯度），从而避免生成平滑区域。

* 单色迁移（Monochrome Transfer）

提取源图像的纹理特征（梯度），去除其色彩信息；通过这一处理过程可使融合区域的颜色与目标图像保持一致，并且该方法特别适用于皮肤质地填充。

colorChange(src, mask, result, red_mul, green_mul, blue_mul)

* 局部颜色改变（Local Color Change）

在相同的一幅图像中进行色彩融合时，其后三个参数对应RGB三个通道的权重系数，在数值上通常设置于0.5至2.5之间。其中数值越大则会对图像进行锐化处理。

illuminationChange(src, mask, result, alpha=0.2f, beta=0.4f)

* 局部亮度改变（Local Illumination Change）

对区域内光照增强（illuminaton）效果显著的地方进行优化调整，在具体参数设置上发现当beta值增大时（即β增大时），图像整体平滑程度也随之提高；而当alpha值增大时（即α增大时），图像恢复的程度与原图细节更为接近；其中alpha和beta的取值范围均限定在0至2之间；该方法在加亮或曝光不足的区域表现出较好的效果

textureFlattening(src, mask, result, low_thresh, high_thresh, kernel_size)

* 纹理平整（Texture Flattening）

仅仅保留边缘处的梯度，选中区域的质地没有保留。

2.2 Seamless Cloning Examples （Learn OpenCV）

    %matplotlib inline
    import matplotlib.pyplot as plt
    import math
    
    def plotImgs(idxFig, lstImgs, figsize=(12, 24)):
    
    fig = plt.figure(idxFig, figsize=figsize)
    fig.clf()
    
    intNumImgs = len(lstImgs)
    intNumCols = int(math.sqrt(intNumImgs))
    intNumRows = int(math.ceil(intNumImgs / intNumCols))
    
    idxImg = 0
    for idxRow in range(intNumRows):
        for idxCol in range(intNumCols):
    
            ax = fig.add_subplot(intNumRows, intNumCols, idxImg + 1, frameon=False)
            if lstImgs[idxImg].dtype == "uint8":
                ax.imshow(cv2.cvtColor(lstImgs[idxImg], cv2.COLOR_BGR2RGB))
            else:
                ax.imshow(lstImgs[idxImg])
            ax.get_yaxis().set_visible(False)
            ax.get_xaxis().set_visible(False)
            idxImg += 1
            if idxImg >= intNumImgs:
                break
    
    plt.show()

    # Standard imports
    import cv2
    import numpy as np 
     
    # Read images
    src = cv2.imread("./img/airplane.jpg")
    dst = cv2.imread("./img/sky.jpg")
     
    
    # Create a rough mask around the airplane.
    src_mask = np.zeros(src.shape, src.dtype)
    poly = np.array([[4,80], [30,54], [151,63], [254,37], [298,90], [272,134], [43,122]], np.int32)
    cv2.fillPoly(src_mask, [poly], (255, 255, 255))
     
    # This is where the CENTER of the airplane will be placed
    center = (800,100)
     
    # Clone seamlessly.
    # Normal Cloning
    output1 = cv2.seamlessClone(src, dst, src_mask, center, cv2.NORMAL_CLONE)
    # Mixed Cloning
    output2 = cv2.seamlessClone(src, dst, src_mask, center, cv2.MIXED_CLONE)
    output3 = cv2.seamlessClone(src, dst, src_mask, center, cv2.MONOCHROME_TRANSFER)
    
    # Save result
    #cv2.imwrite("./output/opencv-seamless-cloning-example.jpg", output)

    idxFig = 1
    plotImgs(1, [src, src_mask, dst, output1, output2, output3], figsize=(12, 8))

    # Read images : src image will be cloned into dst
    im = cv2.imread("./img/wood-texture.jpg")
    obj = cv2.imread("./img/iloveyouticket.jpg")
     
    # Create an all white mask
    mask = 255 * np.ones(obj.shape, obj.dtype)
     
    # The location of the center of the src in the dst
    width, height, channels = im.shape
    center = (height // 2, width // 2)
     
    # Seamlessly clone src into dst and put the results in output
    normal_clone = cv2.seamlessClone(obj, im, mask, center, cv2.NORMAL_CLONE)
    mixed_clone = cv2.seamlessClone(obj, im, mask, center, cv2.MIXED_CLONE)
    mono_clone = cv2.seamlessClone(obj, im, mask, center, cv2.MONOCHROME_TRANSFER)
     
    # Write results
    #cv2.imwrite("images/opencv-normal-clone-example.jpg", normal_clone)
    #cv2.imwrite("images/opencv-mixed-clone-example.jpg", mixed_clone)

    idxFig += 1
    plotImgs(1, [obj, im, normal_clone, mixed_clone, mono_clone], figsize=(12, 8))