# Solving poisson equation on image with neumann boundary condition

I am trying to solve a standard Poisson equation on image with Neumann boundary condition. The Poisson equation has the form:

$Δp=b$

where $b$ is the divergence of the input image,which is as shown:

I tried to write it in Matlab with centred finite difference method and Gauss-Seidel iterative method but the result turns out to be very weird, which is as shown:

Can any kind soul help to spot errors in my code below:

{%Specifying parameters
nx=rows;                           %Number of steps in space(x)(vertical)
ny=columns;                        %Number of steps in space(y)(horizontal)
dx=1;
dy=1;
h=dx;

%Solving (Dxx+Dyy)p(i,j)=b(i,j) (eqn 3.16) using centred FDM (5-points difference)
%Gauss-Seidel iterative method
pn=zeros(nx,ny);                 %Preallocating pn
p=zeros(nx,ny);                  %Preallocating p
b=zeros(nx,ny);
ax=0;                            %counter

for j=2:ny-1
for i=2:nx-1
b(i,j)=0.5*(grayImage1(i+1,j)-grayImage1(i-1,j)+grayImage(i,j+1)-grayImage(i,j-1))./h;

end
end

i=2:nx-1;
j=2:ny-1;
for it=1:200000
p(1,1)=0;
pn=p;
p(i,j)=0.25*(((pn(i+1,j)+p(i-1,j)))+((pn(i,j+1)+p(i,j-1)))-h^2*(b(i,j)));

%%boundary conditions
%The four corners
%p(1,1)=1/2*pn(2,1)+1/2*pn(1,2)-1/4*dx^2*b(1,1);
p(1,ny)=1/2*pn(2,ny)+1/2*pn(1,ny-1)-1/4*dx^2*b(1,ny);
p(nx,1)=1/2*pn(nx,2)+1/2*pn(nx-1,1)-1/4*dx^2*b(nx,1);
p(nx,ny)=1/2*pn(nx-1,ny)+1/2*pn(nx,ny-1)-1/4*dx^2*b(nx,ny);

%The four sides
p(1,j)=1/4*(2*pn(2,j)+pn(1,j+1)+pn(1,j-1)-dx^2*b(1,j));
p(nx,j)=1/4*(2*pn(nx-1,j)+pn(nx,j+1)+pn(nx,j-1)-dx^2*b(nx,j));
p(i,1)=1/4*(2*pn(i,2)+pn(i+1,1)+pn(i-1,1)-dy^2*b(i,1));
p(i,ny)=1/4*(2*pn(i,ny-1)+pn(i+1,ny)+pn(i-1,ny)-dy^2*b(i,ny));

if abs(p(i,j)-pn(i,j))<=zeros(nx-2,ny-2)+0.0001
break
end
ax=ax+1
end }


Thanks very much!!

p/s: I do notice that pure Neumann boundary condition will give rise to non-unique solutions. Hence I followed one of the advice from other answers to set one of the corner to zero. Unfortunately it seems to have no effect on the results.

• What exactly do you mean by "the result turns out to be very weird"? Knowing how a solution is wrong often tells a lot about where the bug could be. – Wolfgang Bangerth Jan 5 '17 at 16:39
• As a side note: The proper spelling of the name is "Seidel", not "Seidler". – Wolfgang Bangerth Jan 5 '17 at 16:39

I cannot recommend fixing a corner to zero, that changes the problem you are trying to solve, instead just add before you iterate a line with:

b = b - mean(b);


This is called the Discrete Compatibility Criteria. May I refer to my own answer to the following: Solving a linear equation system with pure Neumann condition

In an effort to make this an answer, the following code will generate the matrix version of your code and let you use e.g. conjugate gradient as a solver. I grabbed it out of some of the code on: http://math.mit.edu/classes/18.086/2008/ It combines two 1d matrices for the all-Neumann Poisson problem into a large matrix for the 2d problem using the Kronecker Product:

maxit = 1000;
tol = 1e-6;

e = ones(nx-4,1);
A1= spdiags([1 -1 1; e*[1 -2 1]; 1 -1 1],-1:1,nx-2,nx-2); %/dx^2;
e = ones(ny-4,1);
A2= spdiags([1 -1 1; e*[1 -2 1]; 1 -1 1],-1:1,ny-2,ny-2); %/dy^2;
A = -( kron(speye(ny-2),A1)+kron(A2,speye(nx-2)) ); clear A1 A2
b = reshape(b,(nx-2)*(ny-2),1);
b = b - mean(b);

x = pcg(A,b,tol,maxit);

x = reshape(x,nx-2,ny-2);