大家好，又见面了，我是你们的朋友全栈君。
>>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
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/141324.html原文链接：https://javaforall.net
