I am trying to implement Poisson image blending as in the paper Poisson Image Editing. This is the task of filling in a masked region of an image by minimizing $$\min_f\int_\Omega \left | \nabla f - \mathbf v\right |^2$$ with $$f|_{\partial\Omega} = f^*_{\partial\Omega}$$ which is equivalent to the Poisson problem $$\Delta f = \mathbf \nabla\cdot \mathbf v$$ I use the same discretization of the problem as presented there, which is equivalent to using the 4-neighbour discrete Laplacian kernel $$\begin{bmatrix} 0 & 1 & 0\\ 1 & -4 & 1\\ 0 & 1 & 0 \end{bmatrix}$$ and solving the linear system $$Lx = b$$ where $x$ is the vector of unknown image values, $L$ is the sparse matrix that applies the above kernel to $x$, approximating the Laplacian, and $b$ is the sum of boundary terms and gradient guides. Specifically, this is the system of equations $$\sum_{q\in N_p\cap\Omega}f_q-\left |N_p\right | f_p=-\sum_{q\in N_p\cap\partial\Omega}f^*_q-\sum_{q\in N_p}v_{pq}$$ where $N_p$ is the set of 4 neighbors of point $p$ (this is the same as the equation 7 in the paper, albeit with opposite signs).
My implementation in Python (scipy) is below. I use a sparse factorization of L to solve the system for each color channel in the image.
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
import argparse
def poisson_problem(image, mask, guide=None, threshold = 0.5):
indices = np.full(mask.shape, -1) #map of coordinates inside the domain to unique indices
invdices = [] #list of coordinates inside the domain
ind = 0
for p,m in np.ndenumerate(mask):
if m > threshold:
indices[p] = ind
ind += 1
invdices.append(p)
N = len(invdices)
b = np.zeros([N, image.shape[2]]) #for RHS of equation
#build sparse matrix
data = []
I = []
J = []
for i,p in enumerate(invdices):
data.append(-4)
I.append(i)
J.append(i)
for dim in (0, 1):
for dir in (-1,1):
q = [*p]
q[dim] += dir
if q[dim] < 0 or q[dim] >= mask.shape[dim]:
continue
j = indices[(*q,)]
if j > -1:
#contribution from inside the domain
data.append(1)
I.append(i)
J.append(j)
else:
b[i,:] = -image[(*q,)] #boundary term (outside domain)
if guide is not None:
b[i,:] -= image[p] - image[(*q,)] #vector guide term
L = sp.csc_matrix((data, (I,J)), shape=(N,N))
return L, b, invdices
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('--circle', action='store_true')
parser.add_argument('--guide', action='store_true')
args = parser.parse_args()
img = np.zeros([256, 256, 3])
img[:128,:,0] = 1
img[129:,:,1] = 1
if args.circle:
inds = np.indices(img.shape, dtype=np.float)[0:2,:,:,0] - np.array([[[img.shape[0]//2]], [[img.shape[1]//2]]])
mask = ((inds[0]*inds[0]+inds[1]*inds[1]) < min(*img.shape[0:2])**2/16).astype(np.float)
else:
mask = np.zeros(img.shape[0:2])
mask[img.shape[0]//4:-img.shape[0]//4, img.shape[1]//4:-img.shape[1]//4] = 1
#solve poisson problem for each color channel (only BCs change)
L, b, I = poisson_problem(img, mask, img if args.guide else None)
factor = linalg.factorized(L)
xs = np.stack([factor(b[:,i]) for i in range(b.shape[1])], 1)
#composite and display results
img2 = np.zeros(img.shape)
for i, p in enumerate(I):
img2[p] = xs[i]
mask = np.expand_dims(mask,2)
img3 = img * (1-mask) + img2 * mask
plt.imshow(img3)
plt.show()
To make sure it works, I'm pasting the original image into the masked region, hoping to get the original image back. When the mask is square and the gradient guides are provided (--guide), it seems to work without noticeable errors:
But when the mask is a circle shape (run with the --circle and --guide arguments), this is the result:
This does not even appear to be preserving the boundary condition except where the mask boundary is tangent to the axes. This tendency towards 0 seems to occur even when solving the simple Laplace equation, with no guide:
I suspect that this is some kind of error resulting from the discretization, but as far as I can tell, this is the same method as in the paper which claims to work on arbitrary domain shapes. What could be causing this?
