24小时热门版块排行榜

返回列表

雨人007

银虫 (小有名气)

应助: 12 (小学生)
金币: 247.5
帖子: 212
在线: 43.7小时
虫号: 2672857
注册: 2013-09-22
性别: GG
专业: 资源化工

[求助] 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

回复此楼