function w = griddata3(x,y,z,v,xi,yi,zi,method)
%GRIDDATA3 Data gridding and hyper-surface fitting for 3-dimensional data.
%   W = GRIDDATAN(X, Y, Z, V, XI, YI, ZI) fits a hyper-surface of the
%   form W = F(X,Y,Z) to the data in the (usually) nonuniformly-spaced
%   vectors (X, Y, Z, V).  GRIDDATA3 interpolates this hyper-surface at
%   the points specified by (XI,YI,ZI) to produce W.
%
%   (XI,YI,ZI) is usually a uniform grid (as produced by MESHGRID) and is
%   where GRIDDATA3 gets its name. 
%
%   [...] = GRIDDATA3(...,'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 triangulation of
%   the data via qhull. 
%
%   See also GRIDDATA, GRIDDATAN, QHULL, DELAUNAYN, MESHGRID

%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 1.5 $  $Date: 2000/08/28 15:53:33 $

if nargin < 7
  error('Need 8 arguments.');
end
if nargin < 8, method = 'linear'; end
if ~isequal(method(1),'l') & ~isequal(method(1),'n')
  error('METHOD must be one of ''linear'', or ''nearest''.');
end

x = x(:); y=y(:); z=z(:); v = v(:);
m = length(x);
if m < 3, error('Not enough points'); end
if m ~= length(y) | m ~= length(z) | m ~= length(v)
  error('X,Y,Z,V must all have the same size.');
end

X = [x y z];

% Sort (x,y,z) so duplicate points can be averaged before passing to delaunay

[X, ind] = sortrows(X);
v = v(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
    v(fe(i)) = mean(v(fs(i):fe(i)));
  end
  X = X(~ind(2:end),:);
  v = v(~ind(2:end));
end

switch lower(method(1)),
  case 'l'
    w = linear(X,v,[xi(:) yi(:) zi(:)]);
  case 'n'
    w = nearest(X,v,[xi(:) yi(:) zi(:)]);
  otherwise
    error('Unknown method.');
end
w = reshape(w,size(xi));


%------------------------------------------------------------
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 triangulated.');
  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) = y(tri(t(i),:),:)'*p(:,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 triangulated.');
  zi = repmat(NaN,size(xi));
  return
end

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

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

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

