>>type griddata
function [xi,yi,zi] = gdatav4(x,y,z,xi,yi)
%GDATAV4 MATLAB 4 GRIDDATA interpolation
% Reference: David T. Sandwell, Biharmonic spline
% interpolation of GEOS-3 and SEASAT altimeter
% data, Geophysical Research Letters, 2, 139-142,
% 1987. Describes interpolation using value or
% gradient of value in any dimension.
xy = x( + y(*sqrt(-1);
% Determine distances between points
d = xy(:,ones(1,length(xy)));
d = abs(d – d.’);
n = size(d,1);
% Replace zeros along diagonal with ones (so these don’t show up in the
% find below or in the Green’s function calculation).
d(1:n+1:numel(d)) = ones(1,n);
non = find(d == 0, 1);
if ~isempty(non),
% If we’ve made it to here, then some points aren’t distinct. Remove
% the non-distinct points by averaging.
[r,c] = find(d == 0);
k = find(r < c);
r = r(k); c = c(k); % Extract unique (row,col) pairs
v = (z(r) + z(c))/2; % Average non-distinct pairs
rep = find(diff(c)==0);
if ~isempty(rep), % More than two points need to be averaged.
runs = find(diff(diff(c)==0)==1)+1;
for i=1:length(runs),
k = (c==c(runs(i))); % All the points in a run
v(runs(i)) = mean(z([r(k);c(runs(i))])); % Average (again)
end
end
z(r) = v;
if ~isempty(rep),
z(r(runs)) = v(runs); % Make sure average is in the dataset
end
% Now remove the extra points.
z(c) = [];
xy(c, = [];
xy(:,c) = [];
d(c, = [];
d(:,c) = [];
% Determine the non distinct points
ndp = sort([r;c]);
ndp(ndp(1:length(ndp)-1)==ndp(2:length(ndp))) = [];
warning(‘MATLAB:griddata:NonDistinctPoints’,[‘Averaged %d non-distinct ‘ …
‘points.\n Indices are: %s.’],length(ndp),num2str(ndp’))
end
% Determine weights for interpolation
g = (d.^2) .* (log(d)-1); % Green’s function.
% Fixup value of Green’s function along diagonal
g(1:size(d,1)+1:numel(d)) = zeros(size(d,1),1);
weights = g \ z(;
[m,n] = size(xi);
zi = zeros(size(xi));
jay = sqrt(-1);
xy = xy.’;
% Evaluate at requested points (xi,yi). Loop to save memory.
for i=1:m
for j=1:n
d = abs(xi(i,j)+jay*yi(i,j) – xy);
mask = find(d == 0);
if ~isempty(mask), d(mask) = ones(length(mask),1); end
g = (d.^2) .* (log(d)-1); % Green’s function.
% Value of Green’s function at zero
if ~isempty(mask), g(mask) = zeros(length(mask),1); end
zi(i,j) = g * weights;
end
end
if nargout<=1,
xi = zi;
end
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:http://bianchenghao.cn/35914.html