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.