function y = extrapolate(x);
% EXTRAPOLATE: this does a very crude extrapolation/interpolation through
% missing values (as marked by NaNs). The array x may be either 2 or 3
% dimensional. If it is 2 dimensional then missing values are replaced by
% the nanmean of the of the 8 neighbours. If the point is on a boundary
% (i.e., has index 1 or end) then there is no replacement. For the 3-d case
% it is treated as a series of 2-d slices (this makes sense if the 3rd
% dimesnion is time).
%
% CALL: function y = extrapolate(x);
%
% Tested by t_extrapolate.m.

% $Id:$
% Copyright J. V. Mansbridge, CSIRO, Mon Feb 23 16:09:33 EST 1998

si = size(x);
num_dims = length(si);

if num_dims < 2
  error('x can only be 2 or 3 dimensional')
elseif num_dims == 2
  yy = repmat(NaN, [si 9]);
  yy(:, :, 1) = x;
  yy(2:end-1, 2:end-1, 2) = x(1:end-2, 1:end-2);
  yy(2:end-1, 2:end-1, 3) = x(1:end-2, 2:end-1);
  yy(2:end-1, 2:end-1, 4) = x(1:end-2, 3:end);
  yy(2:end-1, 2:end-1, 5) = x(2:end-1, 1:end-2);
  yy(2:end-1, 2:end-1, 6) = x(2:end-1, 3:end);
  yy(2:end-1, 2:end-1, 7) = x(3:end, 1:end-2);
  yy(2:end-1, 2:end-1, 8) = x(3:end, 2:end-1);
  yy(2:end-1, 2:end-1, 9) = x(3:end, 3:end);
  
  y = x;
  xx = nanmean(yy, 3);
  ff = find(isnan(y));
  y(ff) = xx(ff);
elseif num_dims == 3
  % Deal with an array that is reshaped to be 3 dimensional
  yy = repmat(NaN, [si 9]);
  yy(:, :, :, 1) = x;
  yy(2:end-1, 2:end-1, :, 2) = x(1:end-2, 1:end-2, :);
  yy(2:end-1, 2:end-1, :, 3) = x(1:end-2, 2:end-1, :);
  yy(2:end-1, 2:end-1, :, 4) = x(1:end-2, 3:end, :);
  yy(2:end-1, 2:end-1, :, 5) = x(2:end-1, 1:end-2, :);
  yy(2:end-1, 2:end-1, :, 6) = x(2:end-1, 3:end, :);
  yy(2:end-1, 2:end-1, :, 7) = x(3:end, 1:end-2, :);
  yy(2:end-1, 2:end-1, :, 8) = x(3:end, 2:end-1, :);
  yy(2:end-1, 2:end-1, :, 9) = x(3:end, 3:end, :);
  
  y = x;
  xx = nanmean(yy, 4);
  ff = find(isnan(y));
  y(ff) = xx(ff);
elseif num_dims > 3
  error('x can only be 2 or 3 dimensional')
end
