Warning, /utils/matlab/griddata_preprocess.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4ec37fd8 on 2023-10-03 22:00:32 UTC
a33039f58e Alis*0001 function [del] = griddata_preprocess(x,y,xi,yi,method)
0002 %GRIDDATA_PREPROCESS Pre-calculate Delaunay triangulation for use
0003 % with GRIDDATA_FAST.
0004 %
0005 % DEL = GRIDDATA_PREPROCESS(X,Y,XI,YI)
0006
0007 % Based on
0008 % Clay M. Thompson 8-21-95
0009 % Copyright 1984-2001 The MathWorks, Inc.
0010
4ec37fd829 Jean*0011 narginchk(4,5)
a33039f58e Alis*0012
dc4631a0e8 Jean*0013 %- To interpolate to a section (with position vector xi,yi), skip the
0014 % ndgrid call hereafter. However, since this is poorly documented,
0015 % leave the unconditional call to ndgrid to avoid falling into this trap
0016 % when doing 2-D interpolation on to a square grid (xi & yi of equal size)
0017 %if prod(size(xi)) ~= prod(size(yi))
a33039f58e Alis*0018 [yi,xi]=ndgrid(yi,xi);
dc4631a0e8 Jean*0019 %end
a33039f58e Alis*0020
0021 if nargin<6, method = 'linear'; end
0022 if ~isstr(method),
0023 error('METHOD must be one of ''linear'',''cubic'',''nearest'', or ''v4''.');
0024 end
0025
0026
0027 switch lower(method),
0028 case 'linear'
0029 del = linear(x,y,xi,yi);
0030 % case 'cubic'
0031 % zi = cubic(x,y,z,xi,yi);
0032 % case 'nearest'
0033 % zi = nearest(x,y,z,xi,yi);
0034 % case {'invdist','v4'}
0035 % zi = gdatav4(x,y,z,xi,yi);
0036 otherwise
0037 error('Unknown method.');
0038 end
0039
0040
0041
0042 %------------------------------------------------------------
0043 function delau = linear(x,y,xi,yi)
0044 %LINEAR Triangle-based linear interpolation
0045
0046 % Reference: David F. Watson, "Contouring: A guide
0047 % to the analysis and display of spacial data", Pergamon, 1994.
0048
0049 siz = size(xi);
0050 xi = xi(:); yi = yi(:); % Treat these as columns
0051 x = x(:); y = y(:); % Treat these as columns
0052
4f46ae74ad Mart*0053 % older version of matlab do not have DelaunayTri, later versions do not
0054 % have tsearch
0055 msg=which('DelaunayTri');
0056 if length(msg) > 0 % version is 2012 or newer
0057 % Triangularize the data
0058 tri=DelaunayTri(x,y);
0059 if isempty(tri),
0060 warning('Data cannot be triangulated.');
0061 return
0062 end
0063
0064 % Find the nearest triangle (t)
0065 [t, bcs] = pointLocation(tri, xi, yi);
0066
0067 else % use delaunay and tsearch (and hope it is available)
0068
0069 % Triangularize the data
0070 tri = delaunayn([x y]);
0071 if isempty(tri),
0072 warning('Data cannot be triangulated.');
0073 return
0074 end
0075
0076 % Find the nearest triangle (t)
0077 t = tsearch(x,y,tri,xi,yi);
a33039f58e Alis*0078
4f46ae74ad Mart*0079 end % end of selecting version
a33039f58e Alis*0080
0081 % Only keep the relevant triangles.
0082 out = find(isnan(t));
0083 if ~isempty(out), t(out) = ones(size(out)); end
0084 tri = tri(t,:);
0085
0086 % Compute Barycentric coordinates (w). P. 78 in Watson.
0087 del = (x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
0088 (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1)));
0089 w(:,3) = ((x(tri(:,1))-xi).*(y(tri(:,2))-yi) - ...
0090 (x(tri(:,2))-xi).*(y(tri(:,1))-yi)) ./ del;
0091 w(:,2) = ((x(tri(:,3))-xi).*(y(tri(:,1))-yi) - ...
0092 (x(tri(:,1))-xi).*(y(tri(:,3))-yi)) ./ del;
0093 w(:,1) = ((x(tri(:,2))-xi).*(y(tri(:,3))-yi) - ...
0094 (x(tri(:,3))-xi).*(y(tri(:,2))-yi)) ./ del;
0095 w(out,:) = zeros(length(out),3);
0096
0097 delau.tri=tri;
0098 delau.w=w;
0099 delau.siz=siz;
0100 delau.out=out;
0101
0102 %------------------------------------------------------------
