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

Annotation 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 - (hide 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 enderton 1.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