/[MITgcm]/MITgcm/utils/matlab/griddata_preprocess.m
ViewVC logotype

Contents of /MITgcm/utils/matlab/griddata_preprocess.m

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


Revision 1.2 - (show annotations) (download)
Sat Feb 17 23:49:43 2007 UTC (17 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58w_post, checkpoint60, checkpoint61, checkpoint58x_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint58y_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.1: +4 -1 lines
add $Header:  $ and $Name:  & (for CVS)

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

  ViewVC Help
Powered by ViewVC 1.1.22