function zi = triterp(tri, x, y, z, xi, yi)

% triterp -- Linear interpolation on triangles.
%  triterp(tri, x, y, z, xi, yi) interpolates z(x, y)
%   at points (xi, yi), using the triangulation
%   tri.  If tri is empty, the Matlab "delaunay"
%   routine is called to compute it.
%  triterp(x, y, z, xi, yi) calls triterp([], x, y, z, xi, yi).

%  This routine is cloned directly from the "linear"
%  sub-function found in "griddata.m".

% Note by Rich Signell: If this routine is returning NaN where 
% it should be returning good interpolated values, it is probably
% the fault of "tsearch" (which doesn't work for non Delauney triangles
% or "tsearchsafe" (which is supposed to work for non-Delauney triangles,
% but doesn't always.  We've had cases where for a given mesh, "tsearch"
% fails in some regions, and "tsearchsafe" fails in others.

% Modified blocks below are marked by "ZYDECO":
% Copyright (C) 1999 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
% Version of 25-Jun-1999 09:29:57.
% Updated    23-Jul-1999 14:19:43.

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

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

if nargin < 5 | nargin > 6   % ZYDECO.
	help(mfilename)
	zi = [];
	return
end

if nargin < 5 | nargin > 6   % ZYDECO.
	error(' ## Requires exactly five or six input arguments.')
end

% Prepare for unavailable tri.

if nargin == 5   % ZYDECO.
	yi = xi;
	xi = z;
	z = y;
	y = x;
	x = tri;
	tri = [];
	zi = triterp(tri, x, y, z, xi, yi);
	return
end

% Pad z with NaNs if shorter than the (x, y) data.

if length(z) < length(x)   % ZYDECO.
	z(length(x)) = 0;
	z(length(x)+1:end) = NaN;
end

siz = size(xi);
xi = xi(:); yi = yi(:); % Treat these as columns
x = x(:); y = y(:); % Treat these as columns

% Triangularize the data

if isempty(tri)   % ZYDECO.
	tri = delaunay(x,y,'sorted');
end

if isempty(tri),
  warning('Data cannot be triangulated.');
  zi = repmat(NaN,size(xi));
  return
end

% Find the nearest triangle (t)
t = tsearch(x,y,tri,xi,yi);
%t = tsearchsafe(x,y,tri,xi,yi);

% Only keep the relevant triangles.
out = find(isnan(t));
if ~isempty(out), t(out) = ones(size(out)); end
tri = tri(t,:);

% Compute Barycentric coordinates (w).  P. 78 in Watson.
del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
      (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
w(:,3) = ((x(tri(:,1))-xi).*(y(tri(:,2))-yi) - ...
          (x(tri(:,2))-xi).*(y(tri(:,1))-yi)) ./ del;
w(:,2) = ((x(tri(:,3))-xi).*(y(tri(:,1))-yi) - ...
          (x(tri(:,1))-xi).*(y(tri(:,3))-yi)) ./ del;
w(:,1) = ((x(tri(:,2))-xi).*(y(tri(:,3))-yi) - ...
          (x(tri(:,3))-xi).*(y(tri(:,2))-yi)) ./ del;
w(out,:) = zeros(length(out),3);

z = z(:).'; % Treat z as a row so that code below involving
            % z(tri) works even when tri is 1-by-3.
zi = sum(z(tri) .* w,2);

zi = reshape(zi,siz);

if ~isempty(out), zi(out) = NaN; end