Poisson Blending

This project implements poisson blending to blend a masked foreground object into a background image. The method uses this formula $$\boldsymbol{v} = argmin_{\boldsymbol{v}} \sum_{i\in S, j\in N_i \cap S}((v_i - v_j) - (s_i - s_j))^2 + \sum_{i \in S, j \in N_i \cap\neg S}((v_i - t_j) - (s_i - s_j))^2$$ Apart from the required implementations, I have also implemented mixed blending for bells and whistles.

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import lsqr
import cv2
import imageio
import matplotlib.pyplot as plt
from skimage.transform import SimilarityTransform, warp

Below uses a simpler version of blending to reconstruct a small image.

def toy_recon(image):
    h, w = image.shape
    im2var = np.arange(h * w).reshape((h, w)).astype(int)
    i_list = []
    j_list = []
    val_list = []
    b = []
    e = 0
    for i in range(h - 1):
        for j in range(w - 1):
            # Col Constraint
            i_list.append(e)
            j_list.append(im2var[i, j + 1])
            val_list.append(1)
            i_list.append(e)
            j_list.append(im2var[i, j])
            val_list.append(-1)
            b.append(int(image[i, j + 1]) - int(image[i, j]))
            e += 1
            # Row Constraint
            i_list.append(e)
            j_list.append(im2var[i + 1, j])
            val_list.append(1)
            i_list.append(e)
            j_list.append(im2var[i, j])
            val_list.append(-1)
            b.append(int(image[i + 1, j]) - int(image[i, j]))
            e += 1
    i_list.append(e)
    j_list.append(im2var[0, 0])
    val_list.append(1)
    b.append(int(image[0, 0]))
    e += 1
    # This is added because the bottom right corner of image is not restrained.
    i_list.append(e)
    j_list.append(im2var[h - 1, w - 1])
    val_list.append(1)
    b.append(int(image[h - 1, w - 1]))
    e += 1
    i_array = np.array(i_list)
    j_array = np.array(j_list)
    val_array = np.array(val_list)
    b = np.array(b)
    A = csr_matrix((val_array, (i_array, j_array)), shape=(e, im2var.size))
    x, istop, itn, normr = lsqr(A, b)[:4]
    x = x.reshape((h, w))
    return x

image = imageio.imread('demo_data/toy_problem.png')
image_hat = toy_recon(image)

plt.subplot(121)
plt.imshow(image)
plt.title('Input')
plt.subplot(122)
plt.imshow(image_hat)
plt.title('Output')
plt.show()

Below is the implementation of full poisson blending.
The trick that i used to reduce variable number in least square solver is to only list the pixels inside the mask as variables. To achieve this, I used a dictionary to map its pixel location to its variable index.
Because the method and code of mixed blending is very similar to that of poisson blending, it is written in the same code block.

def poisson_blend(fg, mask, bg):
    """
    Poisson Blending.
    :param fg: (H, W, C) source texture / foreground object
    :param mask: (H, W)
    :param bg: (H, W, C) target image / background
    :return: (H, W, C)
    """
    h, w, channels = fg.shape
    bg_copy = np.copy(bg)
#     im2var = np.arange(h * w).reshape((h, w)).astype(int)
    # The next part creates a dictionary of mask pixels
    mask_1c = mask[:, :, 0]
    mask_dict = dict()
    var_ind = 0
    for i in np.argwhere(mask_1c):
        mask_dict[(i[0], i[1])] = var_ind
        var_ind += 1
    for c in range(channels):
        i_list = []
        j_list = []
        val_list = []
        b = []
        e = 0
        for ((row, col), ind) in mask_dict.items():
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if (row - 1, col) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row - 1, col)])
                val_list.append(-1)
                b.append(fg[row, col, c] - fg[row - 1, col, c])
            else:
                b.append(bg[row - 1, col, c] + fg[row, col, c] - fg[row - 1, col, c])
            e += 1
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if (row + 1, col) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row + 1, col)])
                val_list.append(-1)
                b.append(fg[row, col, c] - fg[row + 1, col, c])
            else:
                b.append(bg[row + 1, col, c] + fg[row, col, c] - fg[row + 1, col, c])
            e += 1
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if (row, col - 1) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row, col - 1)])
                val_list.append(-1)
                b.append(fg[row, col, c] - fg[row, col - 1, c])
            else:
                b.append(bg[row, col - 1, c] + fg[row, col, c] - fg[row, col - 1, c])
            e += 1
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if (row, col + 1) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row, col + 1)])
                val_list.append(-1)
                b.append(fg[row, col, c] - fg[row, col + 1, c])
            else:
                b.append(bg[row, col + 1, c] + fg[row, col, c] - fg[row, col + 1, c])
            e += 1
        i_array = np.array(i_list)
        j_array = np.array(j_list)
        val_array = np.array(val_list)
        b = np.array(b)
        A = csr_matrix((val_array, (i_array, j_array)), shape=(e, len(mask_dict)))
        x, istop, itn, normr = lsqr(A, b)[:4]
        for ((row, col), ind) in mask_dict.items():
            bg_copy[row, col, c] = x[ind]
    return bg_copy

def mixed_blending(fg, mask, bg):
    """
    Poisson Blending.
    :param fg: (H, W, C) source texture / foreground object
    :param mask: (H, W)
    :param bg: (H, W, C) target image / background
    :return: (H, W, C)
    """
    h, w, channels = fg.shape
    bg_copy = np.copy(bg)
#     im2var = np.arange(h * w).reshape((h, w)).astype(int)
    # The next part creates a dictionary of mask pixels
    mask_1c = mask[:, :, 0]
    mask_dict = dict()
    var_ind = 0
    diff = 0
    for i in np.argwhere(mask_1c):
        mask_dict[(i[0], i[1])] = var_ind
        var_ind += 1
    for c in range(channels):
        i_list = []
        j_list = []
        val_list = []
        b = []
        e = 0
        for ((row, col), ind) in mask_dict.items():
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if abs(fg[row, col, c] - fg[row - 1, col, c]) >= abs(bg[row, col, c] - bg[row - 1, col, c]):
                diff = fg[row, col, c] - fg[row - 1, col, c]
            else:
                diff = bg[row, col, c] - bg[row - 1, col, c]
            if (row - 1, col) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row - 1, col)])
                val_list.append(-1)
                b.append(diff)
            else:
                b.append(bg[row - 1, col, c] + diff)
            e += 1
            
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if abs(fg[row, col, c] - fg[row + 1, col, c]) >= abs(bg[row, col, c] - bg[row + 1, col, c]):
                diff = fg[row, col, c] - fg[row + 1, col, c]
            else:
                diff = bg[row, col, c] - bg[row + 1, col, c]
            if (row + 1, col) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row + 1, col)])
                val_list.append(-1)
                b.append(diff)
            else:
                b.append(bg[row + 1, col, c] + diff)
            e += 1
            
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if abs(fg[row, col, c] - fg[row, col - 1, c]) >= abs(bg[row, col, c] - bg[row, col - 1, c]):
                diff = fg[row, col, c] - fg[row, col - 1, c]
            else:
                diff = bg[row, col, c] - bg[row, col - 1, c]
            if (row, col - 1) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row, col - 1)])
                val_list.append(-1)
                b.append(diff)
            else:
                b.append(bg[row, col - 1, c] + diff)
            e += 1
            
            i_list.append(e)
            j_list.append(ind)
            val_list.append(1)
            if abs(fg[row, col, c] - fg[row, col + 1, c]) >= abs(bg[row, col, c] - bg[row, col + 1, c]):
                diff = fg[row, col, c] - fg[row, col + 1, c]
            else:
                diff = bg[row, col, c] - bg[row, col + 1, c]
            if (row, col + 1) in mask_dict:
                i_list.append(e)
                j_list.append(mask_dict[(row, col + 1)])
                val_list.append(-1)
                b.append(diff)
            else:
                b.append(bg[row, col + 1, c] + diff)
            e += 1
            
        i_array = np.array(i_list)
        j_array = np.array(j_list)
        val_array = np.array(val_list)
        b = np.array(b)
        A = csr_matrix((val_array, (i_array, j_array)), shape=(e, len(mask_dict)))
        x, istop, itn, normr = lsqr(A, b)[:4]
        for ((row, col), ind) in mask_dict.items():
            bg_copy[row, col, c] = x[ind]
    return bg_copy

Below is the default demo of poisson blending.

# after alignment (masking_code.py)
ratio = 1
fg = cv2.resize(imageio.imread('demo_data/source_01.jpg'), (0, 0), fx=ratio, fy=ratio)
bg = cv2.resize(imageio.imread('demo_data/target_01.jpg'), (0, 0), fx=ratio, fy=ratio)
mask = cv2.resize(imageio.imread('demo_data/mask_01.jpg'), (0, 0), fx=ratio, fy=ratio)

offset = -np.array([5, 200])
fg, mask = align_source(fg, bg, mask, offset)

fg = fg / fg.max()
bg = bg / bg.max()
mask = (mask > 0)

blend_img = poisson_blend(fg, mask, bg)

plt.rcParams["figure.figsize"] = (18,10)
plt.subplot(221)
plt.imshow(fg)
plt.title('Foreground')
plt.subplot(222)
plt.imshow(bg)
plt.title('Background')
plt.subplot(223)
plt.imshow(fg * mask + bg * (1 - mask))
plt.title('Naive Blend')
plt.subplot(224)
plt.imshow(blend_img)
plt.title('Poisson Blend')
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Below is my favorite blending result, as the blending problem itself is not trivial and the blending result makes the lighting look natural.

# after alignment (masking_code.py)
ratio = 1
fg = cv2.resize(imageio.imread('demo_data/porsche_newsource.png'), (0, 0), fx=ratio, fy=ratio)
bg = cv2.resize(imageio.imread('demo_data/track.jpeg'), (0, 0), fx=ratio, fy=ratio)
mask = cv2.resize(imageio.imread('demo_data/track._mask.png'), (0, 0), fx=ratio, fy=ratio)

# offset = -np.array([5, 200])
# fg, mask = align_source(fg, bg, mask, offset)

fg = fg / fg.max()
bg = bg / bg.max()
mask = (mask > 0)

blend_img = poisson_blend(fg, mask, bg)

plt.rcParams["figure.figsize"] = (18,10)
plt.subplot(221)
plt.imshow(fg)
plt.title('Foreground')
plt.subplot(222)
plt.imshow(bg)
plt.title('Background')
plt.subplot(223)
plt.imshow(fg * mask + bg * (1 - mask))
plt.title('Naive Blend')
plt.subplot(224)
plt.imshow(blend_img)
plt.title('Poisson Blend')
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Below are some more results of poisson blending, the first one produce very good result, however the problem itself is not very difficult.
The second one does not produce aesthetic results, as the grass in the foreground is too bright, when it is blended to the background, it changed the color tone of the robomower entirely, thus the result is not saticfactory.

# after alignment (masking_code.py)
ratio = 1
fg = cv2.resize(imageio.imread('demo_data/letters_newsource.png'), (0, 0), fx=ratio, fy=ratio)
bg = cv2.resize(imageio.imread('demo_data/postit.jpg'), (0, 0), fx=ratio, fy=ratio)
mask = cv2.resize(imageio.imread('demo_data/postit_mask.png'), (0, 0), fx=ratio, fy=ratio)

# offset = -np.array([5, 200])
# fg, mask = align_source(fg, bg, mask, offset)

fg = fg / fg.max()
bg = bg / bg.max()
mask = (mask > 0)

blend_img = poisson_blend(fg, mask, bg)

plt.rcParams["figure.figsize"] = (18,10)
plt.subplot(221)
plt.imshow(fg)
plt.title('Foreground')
plt.subplot(222)
plt.imshow(bg)
plt.title('Background')
plt.subplot(223)
plt.imshow(fg * mask + bg * (1 - mask))
plt.title('Naive Blend')
plt.subplot(224)
plt.imshow(blend_img)
plt.title('Poisson Blend')
plt.show()


# after alignment (masking_code.py)
ratio = 1
fg = cv2.resize(imageio.imread('demo_data/robomower_newsource.png'), (0, 0), fx=ratio, fy=ratio)
bg = cv2.resize(imageio.imread('demo_data/cmu.jpg'), (0, 0), fx=ratio, fy=ratio)
mask = cv2.resize(imageio.imread('demo_data/cmu_mask.png'), (0, 0), fx=ratio, fy=ratio)

# offset = -np.array([5, 200])
# fg, mask = align_source(fg, bg, mask, offset)

fg = fg / fg.max()
bg = bg / bg.max()
mask = (mask > 0)

blend_img = poisson_blend(fg, mask, bg)

plt.rcParams["figure.figsize"] = (18,10)
plt.subplot(221)
plt.imshow(fg)
plt.title('Foreground')
plt.subplot(222)
plt.imshow(bg)
plt.title('Background')
plt.subplot(223)
plt.imshow(fg * mask + bg * (1 - mask))
plt.title('Naive Blend')
plt.subplot(224)
plt.imshow(blend_img)
plt.title('Poisson Blend')
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Below is a demonstration of mixed blending.
The benefit of mixed blending is that it preserves some of the background texture while blending the image. The final result looks like some sort of mural painted on the wall.

# after alignment (masking_code.py)
ratio = 1
fg = cv2.resize(imageio.imread('demo_data/monalisa_newsource.png'), (0, 0), fx=ratio, fy=ratio)
bg = cv2.resize(imageio.imread('demo_data/wall.jpg'), (0, 0), fx=ratio, fy=ratio)
mask = cv2.resize(imageio.imread('demo_data/wall_mask.png'), (0, 0), fx=ratio, fy=ratio)

# offset = -np.array([5, 200])
# fg, mask = align_source(fg, bg, mask, offset)

fg = fg / fg.max()
bg = bg / bg.max()
mask = (mask > 0)

blend_img = mixed_blending(fg, mask, bg)

plt.rcParams["figure.figsize"] = (18,10)
plt.subplot(221)
plt.imshow(fg)
plt.title('Foreground')
plt.subplot(222)
plt.imshow(bg)
plt.title('Background')
plt.subplot(223)
plt.imshow(fg * mask + bg * (1 - mask))
plt.title('Naive Blend')
plt.subplot(224)
plt.imshow(blend_img)
plt.title('Poisson Blend')
plt.show()

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).