function yi = griddatan(x, y, xi, method)
%GRIDDATAN Data gridding and hyper-surface fitting (dimension >= 2).
%   YI = GRIDDATAN(X, Y, XI) fits a hyper-surface of the form Y = F(X) to the
%   data in the (usually) nonuniformly-spaced vectors (X, Y).  GRIDDATAN
%   interpolates this hyper-surface at the points specified by XI to
%   produce YI. XI can be nonuniform.
%
%   X is of dimension m-by-n, representing m points in n-D space. Y is
%   of dimension m-by-1, representing m values of the hyper-surface
%   F(X). XI is a vector of size p-by-n, representing p points in the n-D
%   space whose surface value is to be fitted. YI is a vector of length
%   p approximating the values F(XI).  The hyper-surface always
%   goes through the data points (X,Y).  XI is usually a uniform grid
%   (as produced by MESHGRID).
%
%   [...] = GRIDDATAN(...,'method') where 'method' is one of
%       'linear'    - Tessellation-based linear interpolation (default).
%       'nearest'   - Nearest neighbor interpolation.
%   defines the type of surface fit to the data. 
%   All the methods are based on a Delaunay tessellation of the data.
%
%   See also GRIDDATA, GRIDDATA3, QHULL, DELAUNAYN, MESHGRID.

%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 1.10 $  $Date: 2000/09/01 16:21:38 $

[m,n] = size(x);
if m < n+1, error('Not enough points'); end
if m ~= size(y,1)
  error('Y must have the same number of points as X.');
end

if any(isinf(x(:)) | isnan(x(:)))
  error('Data containing Inf or NaN cannot be tessellated.');
end

if nargin < 4, method = 'linear'; end
if ~isstr(method), 
  error('METHOD must be one of ''linear'',''cubic'',''nearest'', or ''v4''.');
end


% Sort x and y so duplicate points can be averaged before passing to delaunay

[x, ind] = sortrows(x);
y = y(ind);
ind = all(diff(x)'==0);
if any(ind)
  warning('Duplicate x data points detected: using average of the y values');
  ind = [0 ind];
  ind1 = diff(ind);
  fs = find(ind1==1);
  fe = find(ind1==-1);
  if fs(end) == length(ind1) % add an extra term if the last one start at end
     fe = [fe fs(end)+1];
  end
  
  for i = 1 : length(fs)
    % averaging z values
    y(fe(i)) = mean(y(fs(i):fe(i)));
  end
  x = x(~[ind(2:end) 0],:);
  y = y(~[ind(2:end) 0]);
end

switch lower(method),
  case 'linear'
    yi = linear(x,y,xi);
  case 'nearest'
    yi = nearest(x,y,xi);
  otherwise
    error('Unknown method.');
end

%------------------------------------------------------------
function zi = linear(x,y,xi)
%LINEAR Triangle-based linear interpolation

%   Reference: David F. Watson, "Contouring: A guide
%   to the analysis and display of spacial data", Pergamon, 1994.

% Triangularize the data
tri = delaunayn(x);
if isempty(tri),
  warning('Data cannot be tessellated.');
  zi = NaN*zeros(size(xi));
  return
end

% Find the nearest triangle (t)
[t,p] = tsearchn(x,tri,xi);

m1 = size(xi,1);
onev = ones(1,size(x,2)+1);
zi = NaN*zeros(m1,1);

for i = 1:m1
  if ~isnan(t(i))
%    zi(i) = y(tri(t(i),:),:)'*([onev; x(tri(t(i),:),:)']\[1; xi(i,:)']);
     zi(i) = p(i,:)*y(tri(t(i),:));
  end
end

%------------------------------------------------------------

%------------------------------------------------------------
function zi = nearest(x,y,xi)
%NEAREST Triangle-based nearest neightbor interpolation

%   Reference: David F. Watson, "Contouring: A guide
%   to the analysis and display of spacial data", Pergamon, 1994.

% Triangularize the data
tri = delaunayn(x);
if isempty(tri), 
  warning('Data cannot be tessellated.');
  zi = repmat(NaN,size(xi));
  return
end

% Find the nearest vertex
k = dsearchn(x,tri,xi);

zi = k;
d = find(isfinite(k));
zi(d) = y(k(d));

%----------------------------------------------------------
