matlab-grid-generator/bin/griddata_preprocess.m

function [del] = griddata_preprocess(x,y,xi,yi,method)
%GRIDDATA_PREPROCESS Pre-calculate Delaunay triangulation for use
%   with GRIDDATA_FAST.
%
%   DEL = GRIDDATA_PREPROCESS(X,Y,XI,YI)

%   Clay M. Thompson 8-21-95
%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 5.28 $  $Date: 2001/04/15 11:59:14 $

error(nargchk(4,5,nargin))

if prod(size(xi)) ~= prod(size(yi))
 [yi,xi]=ndgrid(yi,xi);
end

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


switch lower(method),
  case 'linear'
    del = linear(x,y,xi,yi);
% case 'cubic'
%   zi = cubic(x,y,z,xi,yi);
% case 'nearest'
%   zi = nearest(x,y,z,xi,yi);
% case {'invdist','v4'}
%   zi = gdatav4(x,y,z,xi,yi);
  otherwise
    error('Unknown method.');
end
  

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

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

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

% Triangularize the data
tri = delaunayn([x y]);
if isempty(tri),
  warning('Data cannot be triangulated.');
  return
end

% Find the nearest triangle (t)
t = tsearch(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);

delau.tri=tri;
delau.w=w;
delau.siz=siz;
delau.out=out;

%------------------------------------------------------------
1	function [del] = griddata_preprocess(x,y,xi,yi,method)
2	%GRIDDATA_PREPROCESS Pre-calculate Delaunay triangulation for use
3	% with GRIDDATA_FAST.
4	%
5	% DEL = GRIDDATA_PREPROCESS(X,Y,XI,YI)
6
7	% Clay M. Thompson 8-21-95
8	% Copyright 1984-2001 The MathWorks, Inc.
9	% $Revision: 5.28 $ $Date: 2001/04/15 11:59:14 $
10
11	error(nargchk(4,5,nargin))
12
13	if prod(size(xi)) ~= prod(size(yi))
14	[yi,xi]=ndgrid(yi,xi);
15	end
16
17	if nargin<6, method = 'linear'; end
18	if ~isstr(method),
19	error('METHOD must be one of ''linear'',''cubic'',''nearest'', or ''v4''.');
20	end
21
22
23	switch lower(method),
24	case 'linear'
25	del = linear(x,y,xi,yi);
26	% case 'cubic'
27	% zi = cubic(x,y,z,xi,yi);
28	% case 'nearest'
29	% zi = nearest(x,y,z,xi,yi);
30	% case {'invdist','v4'}
31	% zi = gdatav4(x,y,z,xi,yi);
32	otherwise
33	error('Unknown method.');
34	end
35
36
37
38	%------------------------------------------------------------
39	function delau = linear(x,y,xi,yi)
40	%LINEAR Triangle-based linear interpolation
41
42	% Reference: David F. Watson, "Contouring: A guide
43	% to the analysis and display of spacial data", Pergamon, 1994.
44
45	siz = size(xi);
46	xi = xi(:); yi = yi(:); % Treat these as columns
47	x = x(:); y = y(:); % Treat these as columns
48
49	% Triangularize the data
50	tri = delaunayn([x y]);
51	if isempty(tri),
52	warning('Data cannot be triangulated.');
53	return
54	end
55
56	% Find the nearest triangle (t)
57	t = tsearch(x,y,tri,xi,yi);
58
59	% Only keep the relevant triangles.
60	out = find(isnan(t));
61	if ~isempty(out), t(out) = ones(size(out)); end
62	tri = tri(t,:);
63
64	% Compute Barycentric coordinates (w). P. 78 in Watson.
65	del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
66	(x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
67	w(:,3) = ((x(tri(:,1))-xi).*(y(tri(:,2))-yi) - ...
68	(x(tri(:,2))-xi).*(y(tri(:,1))-yi)) ./ del;
69	w(:,2) = ((x(tri(:,3))-xi).*(y(tri(:,1))-yi) - ...
70	(x(tri(:,1))-xi).*(y(tri(:,3))-yi)) ./ del;
71	w(:,1) = ((x(tri(:,2))-xi).*(y(tri(:,3))-yi) - ...
72	(x(tri(:,3))-xi).*(y(tri(:,2))-yi)) ./ del;
73	w(out,:) = zeros(length(out),3);
74
75	delau.tri=tri;
76	delau.w=w;
77	delau.siz=siz;
78	delau.out=out;
79
80	%------------------------------------------------------------