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
    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):
        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]:
                j = indices[(*q,)]
                if j > -1:
                    #contribution from inside the domain
                    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)
        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

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: red-green striped image with square region masked out and filled using Poisson image blending But when the mask is a circle shape (run with the --circle and --guide arguments), this is the result: red-green striped image with circular region masked out and filled using Poisson image blending 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:

enter image description here

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?


1 Answer 1


It seems I forgot to accumulate the boundary terms for each relevant entry in the right hand side vector $b$. So

b[i,:] = -image[(*q,)] #boundary term (outside domain)

should have been

b[i,:] -= image[(*q,)]

That fixed it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.