function [y, num_NaNs] = fill_gaps_1d(x)
% fill_gaps_1d: replace the NaNs in a vector with best guesses.
%   [y, num_NaNs] = fill_gaps_1d(x)
% The missing values are found by linear interpolation. NaNs at the end of
% the vector are found from the nearest non-NaN. x may be a row or column
% vector. If x is a matrix or multi-dimensional array the interpolation is
% done along the first index (down the column in a matrix).

% $Id: fill_gaps_1d.m,v 1.4 2001/11/20 04:39:57 man133 Exp $
% Copyright J. V. Mansbridge, CSIRO, Thu Nov 15 10:34:22 EST 2001

% determine array_type:
%  1: row vector
%  2: column vector
%  3: matrix
%  4: multi-dimensional array

si = size(x);
if length(si) == 2
  if si(1) == 1 % turn row vector into a column vector
    x = x(:);
    array_type = 1;
    len_vec = si(2);
  elseif si(2) == 1
    array_type = 2;
    len_vec = si(1);
  else
    array_type = 3;
    dim_big = si(2);
  end
else % turn multi-dimensional x into a matrix
  dim_big = prod(si(2:end));
  x = reshape(x, si(1), dim_big);
  array_type = 4;
end

if (array_type == 1) | (array_type == 2) % vector calcs
  ff = find(~isnan(x));
  len_ff = length(ff);
  if (len_ff == 0) | (len_ff == len_vec)
    y = reshape(x, si);
    num_NaNs = len_vec - len_ff;
    return
  end
   % Apply linear interpolation to interior values then nearest neighbour fit
   % to any missing end-points.
  ind = 1:len_vec;
  y = interp1(ind(ff), x(ff), ind);
  ff = find(~isnan(y));
  y = interp1(ind(ff), y(ff), ind, 'nearest', 'extrap');
  y = reshape(y, si);
  num_NaNs = len_vec - len_ff;
else % matrix calcs
  
  % Make a bigger matrix with a row of NaNs at the top. Then turn this matrix
  % into a column vector with a NaN separating each of the original columns. The
  % original NaNs are then replaced by linear interpolation. The extra NaNs then
  % ensure that the first pass with a linear interpolation puts in a NaN. The
  % NaNs at the ends of the columns can then be gotten rid of with the
  % subsequent nearest neighbour interpolation/extrapolation.  Note that this is
  % done in one big column so that there are no loops in the call to interp1
 
  x_ex = [repmat(0, [1 dim_big]); x];
  ff_bad = find(isnan(x_ex));
  ff_good = find(~isnan(x_ex));
  x_ex(1, :) = NaN;
  x_ex = x_ex(:);
  ind = 1:length(x_ex);
  rep = interp1(ind(ff_good), x_ex(ff_good), ind(ff_bad));
  y_int = x_ex;
  y_int(ff_bad) = rep;
  y_int = reshape(y_int, [(si(1) + 1) dim_big]);
  y_int = y_int(2:end, :);
  
  y_int = y_int(:);
  ind = reshape((1:length(y_int)), si) + repmat(si(1)*(1:dim_big), si(1), 1);
  ind = ind(:);
  ff_bad = find(isnan(y_int));
  ff_good = find(~isnan(y_int));
  rep = interp1(ind(ff_good), y_int(ff_good), ind(ff_bad), 'nearest', 'extrap');
  y = y_int;
  y(ff_bad) = rep;
  
  % Check for a whole column of NaNs and make sure that it is not replaced
  % with interpolated values.
  
  y = reshape(y, [si(1) dim_big]);
  count = sum(isnan(x));
  ff = find(count == si(1));
  y(:, ff) = NaN;
  y = reshape(y, si);
  num_NaNs = length(find(isnan(x)));
end
