function new_vals = find_mean_neighbours(index_pts, array_2d)
% FIND_MEAN_NEIGHBOURS: find mean values of neighbours of some points in matrix
%
%  function new_vals = find_mean_neighbours(index_pts, array_2d)
%
%  The user passes a 2-d array (ARRAY_2D) and a vector of indices into that
%  array (INDEX_PTS). The function returns a vector where the I'th element is
%  the nanmean value of the points in the square around the I'th point
%  specified in INDEX_PTS.
%
%    USAGE:
% Typically we could use FIND_MEAN_NEIGHBOURS to replace the NaNs in an array
% as follows:
%
%   index_pts = find(isnan(array_2d));
%   new_vals = find_mean_neighbours(index_pts, array_2d);
%   new_array(index_pts) = new_vals;
%
% Even better (see note 3) is:
%  
%   index_pts = find(isnan(array_2d));
%   while ~isempty(index_pts)
%     new_vals = find_mean_neighbours(index_pts, array_2d);
%     array_2d(index_pts) = new_vals;
%     index_pts = find(isnan(array_2d));
%   end
%
% The above is the basic code of FILL_GAPS_2D.
%
%    NOTES:
% 1) For interior values of ARRAY_2D there will be 9 values (including the
%    point itself whose nanmean will be taken. ARRAY_2D is padded with extra
%    rows and columns of NaNs so that edge values will be handled correctly.
% 2) The padding mentioned above and the 9 point mean can cause a lot of
%    local memory usage - typically twice the size of ARRAY_2D and 9 times
%    the size of INDEX_PTS. If this is too large then the user can choose the
%    call FIND_MEAN_NEIGHBOURS with only a submatrix of ARRAY_2D and/or
%    INDEX_PTS.
% 3) If there is a region at least 3X3 in size that contains only NaNs then
%    the at least one point will have a NaN for its mean value. However,
%    successive calls to FIND_MEAN_NEIGHBOURS along with replacement must
%    converge.
  
% $Id: template $
% Copyright J. V. Mansbridge, CSIRO, Wed Apr 26 12:31:02 EST 2000

  [m, n] = size(array_2d);
  mpp = m + 2;
  npp = n + 2;
  
  % Create bigger version of array_2d with NaNs around the edges so that
  % there will be no edge problems with the calculations.
  
  xx = repmat(NaN, mpp, npp);
  xx(2:mpp-1, 2:npp-1) = array_2d;
  j_index_pts = ceil(index_pts/m); % column indices in array_2d
  i_index_pts = index_pts - (j_index_pts - 1)*m; % row indices in array_2d
  j_index_xx = j_index_pts + 1;
  i_index_xx = i_index_pts + 1;
  ff = i_index_xx + (j_index_xx - 1)*mpp;
  
  % Make up a matrix of the 'neighbouring' points to each point in the vector
  % index_pts and then take the nanmean of this to find the new
  % interpolated/extrapolated point.
  
  len_ff = length(ff);
  tmp = repmat(NaN, len_ff, 9);
  tmp(:, 1) = xx(ff - mpp - 1);
  tmp(:, 2) = xx(ff - mpp);
  tmp(:, 3) = xx(ff - mpp + 1);
  tmp(:, 4) = xx(ff - 1);
  tmp(:, 5) = xx(ff);
  tmp(:, 6) = xx(ff + 1);
  tmp(:, 7) = xx(ff + mpp - 1);
  tmp(:, 8) = xx(ff + mpp);
  tmp(:, 9) = xx(ff + mpp + 1);
  new_vals = nanmean(tmp, 2);
  