function y = smooth_mat(x)
% smooth_mat: Smooths a matrix by taking weighted mean of neighbours
% function y = smooth_mat(x)
% Algorithm: Y(i, j) is equivalent to sum(sum(X(i-1:i+1, j-1:j+1), 1), 2),
% i.e., each value in Y is equivalent to the sum of the values at the 9
% points surrounding the same point in X. Only the non-NaN points are used
% and points outside of the boundaries are ignored. The algoritm is
% vectorized.
  
% $Id: smooth_mat.m,v 1.1 2000/08/01 01:49:53 mansbrid Exp $
% Copyright J. V. Mansbridge, CSIRO, Mon Jul 31 17:57:11 EST 2000

% Initialise
  
[m, n] = size(x);
mp = m + 1;
mpp = m + 2;
np = n + 1;
npp = n + 2;

% Make z, the same as x but with an outer layer of NaNs. Then make z_ok to
% have 0 where z has a NaN and 1 where it does not. Next replace each NaN in
% z with a 0.

z = repmat(NaN, mpp, npp);
z(2:mp, 2:np) = x;
z_ok = ~isnan(z);
ff = find(isnan(z));
z(ff) = repmat(0, size(ff));

% Make s, the sum of the 9 neighbouring points of z. Since s_ok is the
% equivalent sum of z_ok it records how many non-NaNs there were in each of
% the 9 element blocks. We can then find the mean values by division. We take
% into account those blocks with 9 NaNs in them, replacing their values with
% NaNs.

s = z(1:m, 1:n) + z(2:mp, 1:n) + z(3:mpp, 1:n) + ...
    z(1:m, 2:np) + z(2:mp, 2:np) + z(3:mpp, 2:np) + ...
    z(1:m, 3:npp) + z(2:mp, 3:npp) + z(3:mpp, 3:npp);
s_ok = z_ok(1:m, 1:n) + z_ok(2:mp, 1:n) + z_ok(3:mpp, 1:n) + ...
    z_ok(1:m, 2:np) + z_ok(2:mp, 2:np) + z_ok(3:mpp, 2:np) + ...
    z_ok(1:m, 3:npp) + z_ok(2:mp, 3:npp) + z_ok(3:mpp, 3:npp);
ff_good = find(s_ok > 0);
ff_bad = find(s_ok == 0);
s(ff_good) = s(ff_good)./s_ok(ff_good);
s(ff_bad) = repmat(NaN, size(ff_bad));
y = s;
