# Implementing discrete Poisson equation wtih Neumann boundary condition

I am trying to reconstruct an image from gradients in an arbitrary shaped region of an image. I understand how to implement a discrete 2D poisson solution with Dirchlet boundary conditions. Using http://en.wikipedia.org/wiki/Discrete_Poisson_equation#On_a_two-dimensional_rectangular_grid , you just replace any of the u_i,j that are on the source side of the boundary with the boundary values that you have. However, I don't understand how to implement the Neumann boundary condition. If you want to "match" the external derivative, wouldn't you have to just set the nodes that are on the target side of the boundary to a value that makes the derivative between the source boundary pixel and the target boundary pixel equal the derivative outside the boundary? If you do this, haven't you just simply "filled in" one node around the boundary and provided a slightly smaller version of the Dirchlet problem?

Can anyone explain this?

For image in painting and blending a pixel its new value is dependent on its neighboring pixels. These neighbors may or may not lie on a border. The aim is to interpolate and create smooth transitions. It is thus not only dependent on a border condition, nor does a value on the border get replicated exactly to inside pixels. Should the area between boundaries be exactly 1 pixel wide, this pixel will most likely be the average of the boundary values.

A simple blender algorithm will work like this. Given:

• A background image $\textbf{b}$
• A foreground image $\textbf{f}$
• A mask $\textbf{m}$

where we will put the foreground image $\textbf{f}$ onto the background $\textbf{g}$, where the mask $\textbf{m}$ is true.

To avoid a visible seam we will make a correction image $\textbf{c}$. The run-of-the-mill patching to form new image $\textbf{g}$ then is:

$g(x, y) = \left \{\begin{array}{cc} b(x, y) & \textrm{if } m(x, y) = 0 \\ f(x, y) + c(x, y) & \textrm{if } m(x, y) = 1 \end{array} \right .$

This image $\textbf{c}$ is obviously some form of interpolation, but informally we want two things from $\textbf{c}$:

• $\textbf{f}$+$\textbf{c}$ and $\textbf{b}$ should match on the edges
• It should be 'smooth' as to not create disturbances

Smooth transitions can be attained by simply stipulating that a pixel its new value is the weighted average of its four neighbors (up, down, left and right):

$c(x,y) = \frac{c(x-1, y)+c(x+1, y)+c(x, y-1)+c(x-1, y+1)}{4}$

And for boundary pixels the corresponding values of $\textbf{c}$ are replaced by $\textbf{b}-\textbf{f}$. If for instance for $c(x,y)$ the right neighbour is on de boundary:

$c(x,y) = \frac{c(x-1, y)+c(x+1, y)+c(x, y-1)+(b(x-1, y+1)-f(x-1, y+1))}{4}$

This means that $\textbf{every}$ pixel of $\textbf{f}$ under mask $\textbf{m}$ gives a linear equation, and each pixel in the patch is also a variable.

The result is a square linear system: $\textbf{A}\textbf{c} = \textbf{d}$, which can be solved for $\textbf{c}$.

Note that

• For boundary pixels, the values of $\textbf{b}$ and $\textbf{f}$ are $\textbf{constants}$ in this system, not variables. They thus go `to the right' of the equals sign.
• Pixels outside the masked region, as defined by $\textbf{m}$ are not variables.

• As I understand that, in the inpainting context, you can get superior blending results by taking gradients into account, for which humans are sensitive. For example by taking $\textbf{b} + \nabla_\textbf{f}$ instead of $\textbf{b}-\textbf{f}$. An other alternative is to take $\textbf{b} + \nabla_\textbf{h}$ where $\textbf{h}$ is a mixture of the dominant gradients (up, down..) of $\textbf{f}$ and $\textbf{g}$. And otherwise, I don't know.. – Maurits Mar 12 '12 at 20:05