/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagUtility/griddata_preprocess.m
ViewVC logotype

Contents of /MITgcm_contrib/enderton/Diagnostics/DiagUtility/griddata_preprocess.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Mon Jan 31 15:43:28 2005 UTC (20 years, 5 months ago) by enderton
Branch: MAIN
CVS Tags: HEAD
 o Initial check in.

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

  ViewVC Help
Powered by ViewVC 1.1.22