%	function y=smooth9(x,n)
%
%	nine point smoother, applied n times
%	x is a matrix
%

	function y=smooth9(x,n)

[N,M]=size(x);
xm=zeros(N,M);

if(n~=0)
  
  for l=1:n

x1=x(1:N-2,1:M-2);
x2=x(2:N-1,1:M-2);
x3=x(3:N,1:M-2);
x4=x(1:N-2,2:M-1);
x5=x(2:N-1,2:M-1);
x6=x(3:N,2:M-1);
x7=x(1:N-2,3:M);
x8=x(2:N-1,3:M);
x9=x(3:N,3:M);

xm(2:N-1,2:M-1)=1/9*(x1+x2+x3+x4+x5+x6+x7+x8+x9);
xm(1,1)=(x(1,1)+x(1,2)+x(2,1)+x(2,2))/4;
xm(N,1)=(x(N,1)+x(N-1,1)+x(N,2)+x(N-1,2))/4;
xm(1,M)=(x(1,M)+x(1,M-1)+x(2,M)+x(2,M-1))/4;
xm(N,M)=(x(N,M)+x(N-1,M)+x(N,M-1)+x(N-1,M-1))/4;

xm(2:N-1,1)=(x(1:N-2,1)+x(2:N-1,1)+x(3:N,1)+x(1:N-2,2)+x(2:N-1,2)+x(3:N,2))/6;
xm(2:N-1,M)=(x(1:N-2,M)+x(2:N-1,M)+x(3:N,M)+x(1:N-2,M-1)+x(2:N-1,M-1)+x(3:N,M-1))/6;
xm(1,2:M-1)=(x(1,1:M-2)+x(1,2:M-1)+x(1,3:M)+x(2,1:M-2)+x(2,2:M-1)+x(3,3:M))/6;
xm(N,2:M-1)=(x(N,1:M-2)+x(N,2:M-1)+x(N,3:M)+x(N-1,1:M-2)+x(N-1,2:M-1)+x(N-1,3:M))/6;

x=xm;

  end

y=xm;

else
  
y=x;

end
