Poisson Blending

By Eric | Published: October 15, 2013

Poisson blending is one of the topics that spent me days trying to understand recently (not fully understand yet), it is a wonderful method, and using wonderful maths, as well. Now I’ll try to explain this method, with as less Math formulae as I can.

ABOUT IMAGE BLENDING

As an image blending method, our goal is to seamlessly blend an object or texture from a source image into a target image, here is an example:


we are going to put the dog (or bear??) in 1st image, the two children in 2nd image into the regions in the 3rd image. The simplest way to do that is copy and paste the pixels from 1st and 2nd images into the 3rd image, the result will like this:

Weird, we can see very noticeable difference in these regions, because the color of water in the three sources are different, even if the backgrounds are matched, we can also see seams between these regions. Why?
Our people sight is much more sensitive to gradient than the overall intensity of an image,  that means, we don’t know the exact value of each pixel, but if value changes somewhere, we can easily know that it changes, and this in fact leads to our goal: clone pixels to another image, but let color not change abruptly, i.e. maximally preserve the gradient of the source region. This result is much better than the above one:

POISSON BLENDING OVERVIEW

Poisson blending is one of gradient domain image processing methods, maybe this figure can explain it well:

• v – Gradient of a region in an image (as vector)
• g – Selected region of source (dog, children…)
• f* – Known functions that exist in domain S (water in 3rd source above)
• f – Unknown functions that exist in domain Ω
• Ω – Region g that is now placed on domain S (target background)
• ∂Ω – Boundaries between the source and target regions

So Actually what we are going to do is nothing but given vector field v, find the value of f in unknown region that optimize.


whose solution is the solution of Poisson equation:

where div v = ∂v/∂x + ∂v/∂y,  is the Divergence Property of gradient field, and ∆ is the Laplacian operator:


Applying the Laplacian operator to the Poisson equation, we can get:

div G = -4 f(x,y) + f(x-1,y) + f(x,y-1) + f(x+1,y) + f(x,y+1)

1-D EXAMPLE

Say we are going to copy the red part of left image (4, 3, 5, 4) into the right image.
Fist, we got the gradient of left image, as showed above the red bars, now we are going to move these pixels to the right hand side image, without abruptly change of value, so we want:
    With   f1 = 6, f6 = 1.
that is,
  Let’s denote it Q, so we can get:
 we can solve it by using matrix:

f2 = 6, f3 = 4, f4 = 5, f5 = 3

We can find that, the matrix above is a convolution of kernel (-2 4 -2), because every pixel value is influenced by the two adjacent pixels, this is very similar to the Laplacian operator in 2D situation, every pixel value is influenced by the four adjacent pixels (up, down, left, right).

2-D EXAMPLE

Say we are going to copy some pixels into the white areas of the above image (red pixels are boundary pixels), we doing that by using the equation as same as 1D situation:

Ax = b

in which A is a N*N matrix, N is the number of pixels we are going to copy, in this case, N=10.
We can create the 10*10 matrix as this:

01	for i=1:row number

02	for j=1:col number

03

if(i==j)

04	matrix(i, j)=-4

05	elif(adjacent(pixel(j), pixel(i)))

06	matrix(i, j)=1

07

else

08	matrix(i, j)=0

09

end

10

end

11

end

So the matrix in this case is like:


In the equation, b is a N-elements vector, and is created by:

b[i] = div ( G( Source(x,y) ) ) + Neighbor(target i) ;          i=1..N

div v can be calculated by the formula in front of this article, and  Neighbor means pixel i’s 4 neighbors which belong to the boundary, for example:
Neighbor(pixel 1) are pixels left, up, right;
Neighbor(pixel 8) are pixels left, up;
Neighbor(pixel 5) is pixels up.
and if we know matrix A and vector b, we can calculate vector x, which are the value of pixels in target image inside.
x = A \ b   (‘\’ means matrix left division)
or, x = A^-1 * b

MORE ABOUT POISSON SOLVER

Prof. Gilbert Strang in MIT tells more about FAST POISSON SOLVER, but this only works when the region we want to copy is a square, i.e. with same height and width.
In the above case, the matrix A (as Gilbert’s interpretation, matrix K2D) is in a very beautiful form that we can use the same method of FFT to deal with it, and it can save a lot of time. More details is here: