Gradient Domain Imaging


Submitted in the scope of ‘Computer Graphics Seminar’ course in Winter-Semester 2011 at the TU Berlin, department of CG.




Implement a seamless image blending algorithm (benefiting from the principles of Poisson Image Editing), which enables copying regions from source image to target image, not simply as pixels themselves, but gradients instead, which achieves better results by adjusting copied pixels to the conditions of the new target image !!

Mathematical Overview

Suppose that

g    known function of source
f*    known function in target over domain S
f     unknown function in target over domain Ω
∂Ω boundary between known and unknown regions

Solution idea is to find f whose gradient is the closest to G, gradient of source g.
This can be formulated into :

min ∫∫ ||▼f – g||2

whose solution is the solution of Poisson equation

∆ f = div G = ∆g    (*)

under the condition of Dirichlit-boundary : known pixels at boundaries:

f(x,y) = f0(x,y) ∀(x,y) ∈ ∂Ω

this condition is important, because the known boundary works as a reference to reconstruct gradients into the resulting image. And it also plays a big role in constructing the linear equation system.
The condition is already fulfilled, because pixels at boundaries are read from target image.

From (*): div G is the Divergence Property of gradient field

div G = ∂Gx/∂x + ∂Gy/∂y

and ∆  is the Laplacian operator

|0  1  0|
|1 -4  1|
|0  1  0|

Applying these to (*) equation

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

We turn into a classical linear equation

Ax = b

With number of  Unknowns = N pixels
A: NxN matrix , b: N-elements vector
The problem now is to build both b and A !

Vector b

Suppose Pi denotes the unknown pixel (i) ; i = 1 .. N
And S(x,y) denotes the corresponding pixel (i) in source image

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

Where Neighbor(Pi) is the subset of the 4 neighbors of Pi (top, bottom, right, left) that contains only neighbors belonging to the boundary. Consider the following image for example

The function Neighbor(Pi) for the stated pixel Pi returns the top and right neighbors, that they belong to the boundary at target image.

Matrix A

The first thing to do is to build an index i.e. numerating pixels under the mask

Then we build an NxN matrix A as follows, for each pixel denoted by A’s row number :

  • Put -4 in its column
  • Put 1 in columns referring to its neighbors
  • Put 0 elsewhere

So, we get a symmetric matrix with diagonal = -4But A is huge: usually number of pixels to be processed is > 30,000
=> A : 30,000 x 30,000 = 900,000,000 elements
Problem : How to store and deal with it ?!
Solution : Notice that A is a sparse matrix i.e. populated primarily with zeros, so we benefit from sparse notation in Matlab!
Sparse matrix is stored using Coordinated-list compression method, where a sparse matrix (A) is represented with 3 vectors referring to non-zero elements
  1. Row : row index
  2. Column : column index
  3. Value : value at the index A(row, column)

An example considering the last constructed matrix A, at index (1,1) there’s the value (-4), and at index (1,2) there’s the value (1) .. etc.
Row = [ 1, 1, 1, 2, 2, ..]
Col = [ 1, 2, 4, 1, 2, ..]
Val = [-4, 1, 1, 1,-4, ..]
Number of non-zero elements = O(5xN); 5 = pixel itself + 4 neighbors => store 3x5xN instead of NxN

Solving Equation

Once A and b are there, X can be easily computed

Xk = A \ bk ; k = R, G and B
where \ is the mathematical “Matrix Left Division” operation, which is equivalent to x = b * A-1. But Matlab internally applies the “Direct Solver” approach.Finally, combine the 3 channels and re-arrange the resulting vector X into the mask-form.


For easier usability we created a GUI. It allows the user to specify the location of the source and target images.
The next step is initiated by pressing the “Start Computation” button.
Trough lasso selection the area is selected wich is to be copied into the target image. The lasso selection has auto completion and thereby simplifies the selection process.

In the next step the position of the source image in the target image needs to be specified. The crosshairs determine the position of the upper left corner of the source image.
Left-click to the desired position shifts the source image’s upper left corner to the position of the crosshairs. Subsequently by right-clicking the position is locked and the computation process initiated

When finished the result is shown.

In order to show the possibilities of the algorithm let’s have a look at some examples that work and those that don’t.

Both images have a similar setting and are taken with similar ambient light conditions. That makes it easier for the algorithm to compute a good result.
Also the overall smoothness of the images allows the algorithm to create the illusion of seamless blending.

In this example the background of both the source and the target picture have a very different structure and extremely different color. Additionally the ambient light conditions differ quiet a bit.
The source image is taken under the glaring flood light in a football stadium. While the target image is taken under the very soft and warm artificial ambient light of a restaurant. All these factors lead to a very grainy and discolored result.
The edges of the source image are still visible and not at all blended with the background.

The images share similarities though color and ambient lighting don’t match. Nonetheless the result is seamlessly blended and appears natural.
The structures of the stones around the tracks and the leafs on the forrest floor are alike thus allowing for easy blending.

This example shows the pasting of three source images into one target image. Because of similar ambient lighting and settings as well as homogeneity of the backgrounds a neary perfect result is achieved. All edges seem to blend seamlessly and no discolorations can be found.

Source Code
You can copy and reuse our source code, but referencing us would be a nice gesture 🙂

Amit Agrawal and Ramesh Raskar, Gradient Domain Manipulation, Cambridge, USA, 2007
Patrick Pèrez, Michel Gangnet and Andrew Blake, Poisson Image Editing, Microsoft Research UK, 2003

I welcome any comment or question 😀

24 Responses

[…] Gradient Domain Imaging by Merza Klaghstan […]

  • I’m trying to implement poisson image editing and came across your blog. It is very useful. But I have the following doubts,

    1. In vector b calculation, S(x,y) is the corresponding pixel value of the source image. It’s value is going to be a number between 0 to 255. How can you apply divergence to a number?

    2. In the same vector b calculation, the function Neighbour(pi) returns the neighbours belonging to the boundary. According to your image example given in your blog, the function Neighbour(pi) on pixels 8, 9, 10 and 11 will return NULL as their neighbouring pixels doesnot lie in the boundary?

    3. Why are you using values -4, 1 and 0 for the matrix?

    It would be really helpful if you can answer these questions.

    • Thanks for your contact. Concerning your questions:
      1. Divergence is not calculated from the value of the pixel, but from the gradient at the that pixel. Gradient itself is calculated around the pixel.
      2. You are absolutely right. We are only concerned with pixels near to boundary.
      3. These are the values of Laplacian operator.

      I hope this answers your doubts !

  • Mr. Merza,

    Thank you for your reply. I’m clear with doubt 2 and 3. But still I’m not clear with doubt 1. So, this is how I’m suppose to calculate the gradient for a pixel,

    Horizontal_gradient[row][col] = pixel[row][col+1] – pixel[row][col-1];

    Vertical_gradient[row][col] = pixel[row+1][col] – pixel[row-1][col];

    These 2 values are going to be in the range -255 to 255. So, scaling these values back to 0 to 255.

    gradient[row][col] = Horizontal_gradient[row][col] + Vertical_gradient[row][col];

    Is this correct? If yes, I don’t know what to do next.

    Is it possible to post some kind of Pseudo code for these steps?

    • please refer to the downloadable source code. It’s written in simple matlab, along with enough detailed comments. Calculating (b) vector takes place in “cloneImagePoisson.m”

  • Hello Merza,

    It seems your program is still containing a bug. I try to lasso a small area (for example a ball), and I choose the target close to corner (for example bottom / top left corner). The result, the integration is not on the area what I expected . Would you fix it, please?

    • Hey,

      thanks for your comment.
      Unfortunately this project is 2 years old, and not open for further developments !
      However the code is well commented, and free to be reused 🙂

    • Hi Jelas,
      Did you ever debug the code? I’m having the same problem with the integration not being in the area I expected.

  • Hi,
    I am grateful for you blog about gradient domain imaging, but the source code can not be downloaded.
    May you fix it?

  • Hi Merza, I have three questions.
    Question 1:
    you say that “when we build an NxN matrix A as follows, for each pixel denoted by A’s row number :
    o Put -4 in its column
    o Put 1 in columns referring to its neighbors
    o Put 0 elsewhere
    However, in your code:
    A(count, colIndex) = -1; A(count, count) = 4;
    That is, you put 4, not -4 in its column, and put -1, not 1 in its neighbor. Why?

    Question 2:
    When constructing the laplacian image, your code is:
    imLaplacianR = conv2(imSrcR, -laplacian, ‘same’);
    Why you negative the laplacian?

    Quesion 3:
    your code create the boundary condition b as:
    b = zeros(3, n);
    and solve for X as:
    xR = A\( b(1,:)’ );
    I think if we define b as:
    b = zeros(n, 3); b(count, 1) = ; b(count, 2) = ; b(count, 3) = ;
    then xR = A\(b(1, :);
    Thus we can avoid the transpose operation.

  • i have a prob with this source code…i run thi GUI.m file and select both target and source images but when start computation figure source open nd i select the portion but further its not copying and dnt knw what to do

    • Hi Aisha,

      have you followed the steps accurately ?!
      select a region from source code -> paste on destination using left-click -> right click to fix !

      • yes i have but any click doesnt take it to next step mean after selection nothing works…iam using MATLAB R2008a and in command prompt error occurs of undefined method createmask

        • well, it might take me a while until I can reply you on this !
          but please keep trying till then, and the question is however open for other viewers, that might be able to contribute to this 🙂

  • hi
    please tell me what is offset value ? i am confused ,

  • hi
    i used your cloneimagepoisson code but when i run this program they cannot give me result some time offset value problem and some time “Too many input arguments” ans. i put all the requirement but still ans is ” Too many input arguments.”
    PLEASE tell me what i do ??? .

  • Hi. can you reupdate the file link?
    It only shows text {filelink=1}

  • Leave a Reply