<jupyter_start><jupyter_text>Poisson matting Algorithm1. Initialize alpha, foreground and background values. Alpha = 1 foreground, 0.5 unknown, 0 background. Approximate foreground and background  in unknown region by correspponding nearest known pixel. Smooth the F_B image with Gaussian filter to supress significant changes due to noise and inaccurate estimation.2. Reconstruct alpha by solving Poisson equation using many iterations of Gauss-Seidel, with overrelaxation. 3. F,B refinement. Set new definite foreground and background based on alpha > 0.95 and alpha < 0.05  respectively. Update unknown F, B using new nearest known F, B pixels and apply another Gaussian filter. Iterate steps 2,3 untill change in matting results is small or no new Foreground or background pixels are found in step 3.<jupyter_code>from __future__ import division
import numpy as np
import scipy.misc
import scipy.ndimage
import matplotlib.pyplot as plt
import cv2
import time
from numba import jit
img = scipy.misc.imread('troll.png')
gray_img = scipy.misc.imread('troll.png', flatten='True')
trimap = scipy.misc.imread('trollTrimap.bmp', flatten='True')
# scale = 0.1
# img = scipy.misc.imresize(img, scale)
# gray_img = scipy.misc.imresize(gray_img, scale)
# trimap = scipy.misc.imresize(trimap, scale)
h, w ,c  = img.shape
fg = trimap == 255
bg = trimap == 0
unknown = True ^ np.logical_or(fg,bg)
img = img
fg_img = gray_img*fg
bg_img = gray_img*bg
alphaEstimate = fg + 0.5 * unknown
approx_bg = cv2.inpaint(bg_img.astype(np.uint8),(unknown +fg ).astype(np.uint8)*255,3,cv2.INPAINT_TELEA)*(np.logical_not(fg)).astype(np.float32)
approx_fg = cv2.inpaint(fg_img.astype(np.uint8),(unknown +bg ).astype(np.uint8)*255,3,cv2.INPAINT_TELEA)*(np.logical_not(bg)).astype(np.float32)
plt.imshow(approx_fg, cmap='gray')
plt.show()
plt.imshow(approx_bg, cmap='gray')
plt.show()
# Smooth F - B image
approx_diff = approx_fg - approx_bg
approx_diff = scipy.ndimage.filters.gaussian_filter(approx_diff, 0.9)
plt.imshow(np.minimum(np.maximum(approx_diff,0),255), cmap='gray')
plt.show()<jupyter_output><empty_output><jupyter_text>Poisson Equation$$\Delta \alpha = div \frac{\nabla I}{F-B}$$![](https://wikimedia.org/api/rest_v1/media/math/render/svg/9ae5cade73b31dee0e04eb21b1e6e13b5944ef41)Gauss seidel method![ig](https://wikimedia.org/api/rest_v1/media/math/render/svg/02f26444f4da61ae84a9824a6b5292949a483fcc)<jupyter_code>dy, dx = np.gradient(gray_img)
d2y, _ = np.gradient(dy/approx_diff)
_, d2x = np.gradient(dx/approx_diff)
b = d2y + d2x
def computeAlphaSlow(alpha, b):
    alphaNew = alpha.copy()
    alphaOld = np.zeros(alphaNew.shape)
    threshold = 0.1
    n = 1
    while (n < 50 and np.sum(np.abs(alphaNew - alphaOld)) > threshold):
        alphaOld = alphaNew.copy()
        for i in range(1, h-1):
            for j in range(1, w-1):
                if(unknown[i,j]):
                    alphaNew[i,j] = 1/4  * (alphaNew[i-1 ,j] + alphaNew[i,j-1] + alphaOld[i, j+1] + alphaOld[i+1,j] - b[i,j])
        n +=1
    return alphaNew
@jit
def computeAlphaJit(alpha, b):
    alphaNew = alpha.copy()
    alphaOld = np.zeros(alphaNew.shape)
    threshold = 0.1
    n = 1
    while (n < 50 and np.sum(np.abs(alphaNew - alphaOld)) > threshold):
        alphaOld = alphaNew.copy()
        for i in range(1, h-1):
            for j in range(1, w-1):
                if(unknown[i,j]):
                    alphaNew[i,j] = 1/4  * (alphaNew[i-1 ,j] + alphaNew[i,j-1] + alphaOld[i, j+1] + alphaOld[i+1,j] - b[i,j])
        n +=1
    return alphaNew
def computeAlpha(alpha, b):
    alphaNew = alpha.copy()
    alphaOld = np.zeros(alphaNew.shape)
    threshold = 0.1
    n = 1
    ca = np.transpose(np.nonzero(unknown))
    c = ca[ (ca[:,0]>=1) & (ca[:,0]<h-1) &  (ca[:,1]>=1) & (ca[:,1]<w-1)]
    c0 = c[:,0]
    c1 = c[:,1]
    while (n < 50 and np.sum(np.abs(alphaNew - alphaOld)) > threshold):
        alphaOld = alphaNew.copy()
        alphaNew[c0, c1] = 1/4  * (alphaNew[c0 -1, c1] + alphaNew[c0, c1 -1] + alphaOld[c0,c1+1] + alphaOld[c0+1,c1] - b[c0,c1])
        n +=1
    return alphaNew
%timeit computeAlphaSlow(alphaEstimate, b)
%timeit computeAlphaJit(alphaEstimate, b)
%timeit computeAlpha(alphaEstimate, b)
alpha = computeAlphaJit(alphaEstimate, b)
alpha = np.minimum(np.maximum(alphaNew,0),1).reshape(h,w)
plt.imshow(alpha, cmap='gray')
plt.show()
scipy.misc.imsave('trollAlpha.png', alpha)
plt.imshow((alpha.reshape(h,w,1).repeat(3,2)*img).astype(np.uint8))
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1462979482604817, 'fraction_numerical': 0.03991971454058876, 'mean_word_length': 4.022396416573348, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '109949', 'license': 'Apache-2.0', 'n_tokens_mistral': 1822, 'n_tokens_neox': 1709, 'n_words': 379, 'path': 'misc/Poisson Matting.ipynb', 'repo_name': 'juxiangwu/image-processing', 'size': '4484'}