function new_vals = find_weighted_mean_neighbours(index_pts, array_2d, ...
				weight_2d, weight_x, weight_y, weight_xy)
% FIND_WEIGHTED_MEAN_NEIGHBOURS: find weighted mean values of neighbours of
% some points in matrix
%
%  function new_weighted_vals = find_weighted_mean_neighbours(index_pts, ...
%                       array_2d, weight_2d, weight_x, weight_y, weight_xy)
%
% The user passes 6 identifiers. ARRAY_2D is a 2-d array of dimension m X n.
% INDEX_PTS is a vector of indices into ARRAY_2D. (If INDEX_PTS is actually
% the string 'all' then all of the indices into ARRAY_2D will be used, i.e,
% index_pts = 1:m*n.) WEIGHT_2D is also of dimension m X n, while weight_x and
% weight_y are scalars. If WEIGHT_2D is empty then it will be set to be all
% ones; if it is a scalar then it will be treated as an array of scalars.
%
% The function returns a vector NEW_VALS where the I'th element is a weighted
% mean of values in ARRAY_2D centred in the 9 points forming a square around
% the I'th point specified in INDEX_PTS. (If INDEX_PTS is actually the string
% 'all' then NEW_VALS will be dimensioned m X n.)
%
%    ALGORITHM for calculating the weighted mean
%
% Construct a 9 point square from WEIGHT_2D centred on the the point of
% interest (as specified by the I'th element of INDEX_PTS). Diagramatically
% this is multiplied by : C A C
%                         B 1 B
%                         C A C
% to produce a new 9 element square, say W. Here A = WEIGHT_Y, B = WEIGHT_X,
% C = WEIGHT_XY.
% 
% The weighted mean is found by multiplying W by the corresponding
% elements in ARRAY_2D and then taking the sum, ignoring NaNs. This is then
% normalised by dividing by the sum of the elements in the weighting square
% (counting only those terms that did not multiply by a NaN).
%
%    USAGE:
% Typically we could use FIND_WEIGHTED_MEAN_NEIGHBOURS to replace the NaNs in
% an array as follows:
%
%  >> index_pts = find(isnan(array_2d));
%  >> new_vals = find_weighted_mean_neighbours(index_pts, array_2d, 1, 1, 1, 1);
%  >> array_2d(index_pts) = new_vals;
%
% Even better (see note 3) is:
%  
%  >> index_pts = find(isnan(array_2d));
%  >> while ~isempty(index_pts)
%  >>   new_vals = find_weighted_mean_neighbours(index_pts, array_2d, 1, 1, 1, 1);
%  >>   array_2d(index_pts) = new_vals;
%  >>   index_pts = find(isnan(array_2d));
%  >> end
%
% The other use is to smooth data in the presence of NaNs. You could do:
%  >> array_2d = find_weighted_mean_neighbours('all', array_2d, 1, 1, 1, 1);
%
% 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_WEIGHTED_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
%    at least one point will have a NaN for its mean value. However,
%    successive calls to FIND_WEIGHTED_MEAN_NEIGHBOURS along with replacement
%    must converge.
% 4) Using the weights enables some smoothing of data. However, think
%    carefully before making iterative calls to
%    find_weighted_mean_neighbours. While the the elements of WEIGHT_2D,
%    WEIGHT_X and WEIGHT_Y may be suitable for the first call (they are like
%    covariance values) they may not be correct for later calls as
%    information propagates and NaNs are eliminated.
  
% $Id: find_weighted_mean_neighbours.m,v 1.3 2004/10/05 00:29:40 man133 Exp $
% Copyright J. V. Mansbridge, CSIRO, Wed Apr 26 12:31:02 EST 2000

  [m, n] = size(array_2d);
  mpp = m + 2;
  npp = n + 2;
  get_full = 0;
  if ischar(index_pts)
    if strcmp(lower(index_pts), 'all')
      index_pts = (1:m*n)';
      get_full = 1;
    else
      error(['passed index_points = ' index_pts])
    end
  end
  if isempty(weight_2d)
    weight_2d = ones(m, n);
  end
  if prod(size(weight_2d)) == 1
    weight_2d = repmat(weight_2d, [m n]);
  end
  
  % 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;
  ww = repmat(NaN, mpp, npp);
  ww(2:mpp-1, 2:npp-1) = weight_2d;
  ff = find(isnan(xx));
  ww(ff) = NaN;
  wwxx = ww.*xx;
  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.
  
  % weight_xy = 1/sqrt(1/weight_x.^2 + 1/weight_y.^2);
  len_ff = length(ff);
  tmp = repmat(NaN, len_ff, 9);
  tmp(:, 1) = weight_xy*wwxx(ff - mpp - 1);
  tmp(:, 2) = weight_x*wwxx(ff - mpp);
  tmp(:, 3) = weight_xy*wwxx(ff - mpp + 1);
  tmp(:, 4) = weight_y*wwxx(ff - 1);
  tmp(:, 5) = wwxx(ff);
  tmp(:, 6) = weight_y*wwxx(ff + 1);
  tmp(:, 7) = weight_xy*wwxx(ff + mpp - 1);
  tmp(:, 8) = weight_x*wwxx(ff + mpp);
  tmp(:, 9) = weight_xy*wwxx(ff + mpp + 1);
  
  new_vals = nansum(tmp, 2);
  
  tmp = repmat(NaN, len_ff, 9);
  tmp(:, 1) = weight_xy*ww(ff - mpp - 1);
  tmp(:, 2) = weight_x*ww(ff - mpp);
  tmp(:, 3) = weight_xy*ww(ff - mpp + 1);
  tmp(:, 4) = weight_y*ww(ff - 1);
  tmp(:, 5) = ww(ff);
  tmp(:, 6) = weight_y*ww(ff + 1);
  tmp(:, 7) = weight_xy*ww(ff + mpp - 1);
  tmp(:, 8) = weight_x*ww(ff + mpp);
  tmp(:, 9) = weight_xy*ww(ff + mpp + 1);
  
  new_vals =  new_vals./ nansum(tmp, 2);
  if get_full
    new_vals = reshape(new_vals, m, n);
  end
