matlab griddata,MATLAB 中有关griddata的计算原理，看不懂，求详解！ – 程序语言 – 小木虫 – 学术 科研 互动社区…

>>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