/[MITgcm]/MITgcm_contrib/high_res_cube/matlab-grid-generator/bin/griddata_preprocess.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/matlab-grid-generator/bin/griddata_preprocess.m

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


Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Tue Nov 11 18:08:08 2003 UTC (21 years, 8 months ago) by cnh
Branch: MAIN, initial
CVS Tags: baseline, HEAD
Changes since 1.1: +0 -0 lines
Checking in work done with Dimitri on high-resolution cube gridding and parallel 
communications. 
   o code is in a contrib experiment for now so we can continue collaborating
     on it. However most code is general and will be moved into main branch once 
     it is fully hardened.
   o There are README files in the contrib root and in the subdirectories that
     explain the contents

1 cnh 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     % 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     %------------------------------------------------------------

  ViewVC Help
Powered by ViewVC 1.1.22