Solving poisson equation on image with neumann boundary condition

Question

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.

Answer 1

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);