基于领域梯度的泊松图像融合 (Gradient domain fusion)

实验目的

1. 熟悉泊松融合的意义和用途，掌握泊松融合的基本方法；
2. 理解泊松融合的原理，熟悉基于梯度的图像融合的实现思路；
3. 掌握运用编程实现基于梯度的泊松融合将一个图像粘贴到目标图像中。

实验原理

泊松融合概述

泊松融合可以解决原图像中的部分区域复制到目标图像的区域$\Omega$中，边界过渡不自然的问题。主要思想是：融合需要保证使前景部分区域与背景而言，尽量平滑并且保持边界一致。首先将该部分区域的梯度场覆盖到目标图像的梯度场上，得到融合图像的梯度场，对其求偏导得到散度b，然后通过非边界点与拉普拉斯卷积核进行卷积，并且通过边界点的约束条件，得到稀疏矩阵A，最后通过求解$Af=b$的方程得到$f$,$f$就是融合图像的每个像素点RGB值。

(1)平滑可表示为：$F={\min }_f{\iint }_{\Omega }|\nabla f{|}^2$
其中$\Delta f$表示f二阶微分$|\nabla f{|}^2$，即为直角坐标系下的散度$div$，可通过拉普拉斯算子计算

。因为在一阶微分$\nabla f$取最小值时，二阶微分$\Delta f$的值为0，所以平滑条件可以转化为$div=\Delta f=0$\
(2)保持边界一致表示为：$f{|}_{\partial \Omega }=f^{\ast }{|}_{\partial \Omega }$
$\partial \Omega$表示区域$\Omega$的边界,$f^{\ast }$表示在目标图像中的像素。
最后得到的满足上面条件的泊松方程如下：$$\Delta f=0over\Omega ,f{|}_{\partial \Omega }=f^{\ast }{|}_{\partial \Omega }$$
将该泊松方程转化为线性向量形式如下图：

计算稀疏矩阵A

假设除去边界后要融合的区域像素点总数是N，则矩阵A的大小为N*N。对于每一个像素点，他的拉普拉斯算子的是邻居像素值的和减去4倍的该点像素值。对于邻居像素点不是边界点时，因此每一行有五个值不是0，分别对应拉普拉斯算子的值。如果邻居像素点时边界点是，则该邻居像素点对应的值为0而不是为1。\
例如内容如下所示：

计算向量b

对于每个区域内的像素点,该邻居像素点不是边界点时，则该点的$b_i=-4g_{i,j}+g_{i-1,j}+g_{i+1,j}+g_{i,j-1}+g_{i,j+1}$；如果该像素点的邻居p是边界点时，则该点的$b_i=-4g_{i,j}+g_{i-1,j}+g_{i+1,j}+g_{i,j-1}+g_{i,j+1}-f_p^{\ast }$，其中g表示源图像像素，$f_p^{\ast }$表示目标图像的像素值。
具体实现如下：

from cv2 import cv2
import numpy as np

def neighbors(h, w):
    return (h + 1, w), (h - 1, w), (h, w + 1), (h, w - 1)

def map_Omega(src,mask):#获取源图像需要粘贴的像素点位置以及编号
    h,w=src.shape[0],src.shape[1]
    coordinate_map=[]
    index_map=np.zeros([h,w],np.uint16)
    idx=0
    for i in range(h):
        for j in range(w):
            if mask[i,j].all()==np.array([255, 255, 255]).all():
                coordinate_map.append((i,j)) #记录原图像需要粘贴到目标图像的像素点位置
                index_map[i,j]=idx #将每个需要粘贴的像素点进行编号
                idx+=1
    return index_map,coordinate_map

if __name__ == '__main__':
    center = (220, 240)  #将源图像部分区域移到目标图片的偏移位置

    g= cv2.imread('shake.jpeg').astype(np.float64)  # 源图像
    s = cv2.imread('img_background1.jpg').astype(np.float64)  # 目标图像
    #获取源图像要复制的区域
    poly = np.array([[9,0],[3,11],[17,39],[1,78],[18,71],[33,59],[56,91],[52,122],
                     [87,101],[112,110],[145,107],[170,135],[160,94],[179,45],[168,27],
                     [97,9],[54,44],[36,44],[12,0]], np.int32)

    g_height, g_width, _ = g.shape
    g_mask = np.zeros([g_height, g_width, 3])
    cv2.fillPoly(g_mask, [poly], (255, 255, 255))
    #g_mask[:, :] = [255, 255, 255]
    cv2.imwrite('ig_mask.jpg', g_mask.astype(np.uint8))

    #获取要复制区域的边界，由于边界等于目标图像像素，所以将要复制区域的边界区域mask设为0
    iindex_map,icoordinate_map=map_Omega(g, g_mask)#获取要复制区域的像素点位置并且编号
    back = np.zeros([g_height, g_width, 3])#获取复制部分在目标图像的位置的目标区域的像素
    tmp=[]#记录边界的像素点位置
    for i in range(len(icoordinate_map)):
        y,x=icoordinate_map[i]
        back[y,x]=s[y + center[0], x + center[1]]
        if (y == 0 or y == g_height - 1) or (x == 0 or x == g_width - 1):#判断是否图像边界
            tmp.append((y, x))
            continue
        for neighbor in neighbors(y, x):
            if g_mask[neighbor].all() == np.array([0, 0, 0]).all():  # 判断是否是复制部分的边界
                tmp.append((y, x))
                break
    for i in range(len(tmp)):#将边界的mask设置为0
        y,x=tmp[i]
        g_mask[y, x] = np.array([0, 0, 0])
    cv2.imwrite('back.jpg', back.astype(np.uint8))
    cv2.imwrite('g_mask.jpg', g_mask.astype(np.uint8))

    index_map,coordinate_map = map_Omega(g, g_mask)#获取去除边界后的要复制区域的像素点位置并且编号
    N= len(coordinate_map)
    matrix_A = np.zeros([N, N, 3], np.int8)
    vector_b = np.zeros([N, 3])

    #计算矩阵A和向量b
    for i in range(N):
        y,x=coordinate_map[i]
        temp_b=-4*g[y,x]+g[y+1,x]+g[y-1,x]+g[y,x+1]+g[y,x-1]#计算像素点的梯度，用拉普拉斯计算
        for neighbor in neighbors(y, x):
            if g_mask[neighbor].all() ==np.array([0, 0, 0]).all():#判断邻居像素点是否是边界
                temp_b -= back[neighbor]#如果是邻居像素点是边界点，这该像素点梯度要减去邻居像素点在目标图像的像素值
            else:  # if not
                matrix_A[index_map[y, x], index_map[neighbor]] = np.array([1, 1, 1])
        matrix_A[index_map[y, x], index_map[y, x]] = np.array([-4, -4, -4])#如果邻居不是边界，直接用拉普拉斯计算梯度
        vector_b[index_map[y, x]] = temp_b
    #求解f=b/A
    channels = []
    for _ in range(3):
        a_solution = np.linalg.solve(matrix_A[:, :, _], vector_b[:, _])#求解线性方程矩阵
        a_solution[a_solution < 0] = 0#将小于0的解变为0
        a_solution[a_solution > 255] = 255#将大于255的解变为255
        channels.append(a_solution)#分别获取r，g，三个通道的解
    final_solution = np.dstack(channels).astype(np.uint8)#将channels数组在第三维上进行堆叠


    for i in range(N):
        y,x=coordinate_map[i]
        s[y+ center[0], x + center[1]]=final_solution[0, i]
    cv2.imwrite('test2.jpg', s.astype(np.uint8))

实验结果