function [new, varargout] = regrid2(data, x, y, varargin)
%  1.  NAME
%
%  regrid2 - Regrid two dimensional data
%
%  2.  SYNOPSIS
%
%  [new, grid] = regrid2(data, x, y, grid, method )
%
%  3.  DESCRIPTION
%
%  This Matlab function removes lfirst ines that contain bad data along
%  the specified x and y axes. It generates then x- andy-grid points from
%  the input data, and averages duplicate values at identical grid
%  points. If grid is not specified, x- and y-grids are generated from
%  the x- and y-values in the input data. Data are then interpolated on
%  the new grid points.
%
%  Input data do not need to be sorted and may contain bad values, or
%  duplicate values at identical points.
%
%  4.  PARAMETERS
%
%     data
%        Array containing two dimensional data in xyz-.
%
%     x  Column numbers specifying the x-axis.
%
%     y  Column numbers specifying the y-axis.
%
%     grid optional
%        A cell array of two column vectors specifying the new grid
%        points,where x=grid{1},and y=grid{2}. If not present,a default
%        regular grid is build,starting and ending at the same points as
%        the given axes,and with the same number of points.
%
%     method optional
%        The interpolation method to use for regridding data (default is
%        linear). See Matlab help on function griddata for more
%        information.
%
%  5.  RETURNS
%
%     new
%        An Y*X*V mesh where the values of index v correspond to the
%        columns in the input array.
%
%     grid optional
%        A cell array with two column vectors containing the values of
%        the new grid points,where  x=grid{1} and y=grid{2}.
%
%  6.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

%--- Default values:
grid = [];
method = 'linear';

%--- Arguments:
if (nargin < 3) | (nargin > 6),
	error([ 'Wrong number of arguments' ]);
end
if nargin > 3,
	grid = varargin{1};
end
if nargin > 4,
	method = varargin{2};
end

%--- Check that grid columns are defined:
if (x > size(data,2)) | (y > size(data,2)) | (x < 1) | (y < 1),
	error([ 'Grid columns must be within available columns' ]);
end

%--- Initializations
[ nl, nc ] = size(data);
vv = find( ((1:nc) ~= x) & ((1:nc) ~= y) );

%--- Get the grid and remove bad x- and y-points:
[ grd, data ] = ascii_grd(data, [ x, y ]);
nx = length(grd{1});
ny = length(grd{2});
nl = size(data, 1);

%--- Rebuild the data array si that it has the same size as that defined by
%    the x and y input grids:
data = sortrows(data, [ x, y ]);
% Build mask of duplicate lines:
ii = 2:nl;
ii = (data(ii,x) == data(ii-1, x)) & (data(ii,y) == data(ii-1,y));
ii = cat(1, 0, ii );
% Build mask of first duplicate grid lines:
i0 = 2:nl;
i0 = (ii(i0) & ~ii(i0-1));
i0 = cat(1, 0, i0);
i0 = find(i0);
% Build index of last duplicate grid lines:
i1 = 1:(nl-1);
i1 = (ii(i1) & ~(ii(i1+1)));
i1 = cat(1, i1, ii(nl));
i1 = find(i1);
% This is a control:
if length(i0) ~= length(i1),
	error([ 'Cannot find duplicate ranges on specified grid' ]);
end
% Loop through duplicat ranges:
for i=1:length(i0),
  ii = i0(i):i1(i);
  % Average duplicate grid points:
  for j=vv,
    data(i0(i),j) = mean_nan(data(ii,j));
	end
	% Flag bad rest of positions in this range:
	ii = (i0(i)+1):i1(i);
	data(ii,[x,y]) = NaN;
end
% Filter out bad pad positions:
igood = find(~isnan(data(:,x)) | ~isnan(data(:,y)));
if isempty(igood),
	error([ 'No good data' ]);
end
data = data(igood, :);

%--- Eventually, build a new regular grid:
if isempty(grid),
	for c=1:2,
		grid{c} = new_grid( min(grd{c}), max(grd{c}), length(grd{c}));
	end
end
nx = length(grid{1}); 
ny = length(grid{2}); 

%--- Regrid data:
new = NaN * ones(ny, nx, nc);
[xi, yi] = meshgrid(grid{1}, grid{2});
new(:,:,x) = xi;
new(:,:,y) = yi;
for v=vv,
  igood = find(~isnan(data(:,v)));
  if ~isempty(igood),
    new(:,:,v) = griddata(data(igood,x), data(igood,y), data(igood,v), ...
        xi, yi, method);
  end
end

%--- Final exit:
if nargout > 0,
	varargout{1} = grid;
end

return;





















