| 查看: 1090 | 回复: 3 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
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 |
回复此楼» 猜你喜欢
面上项目filecode邪修
已经有6人回复
-
应该是93bebmhtak前后十一个字符比较关键
已经有14人回复
-
2026国自然放榜时间
已经有4人回复
-
关于代码变化问题，想知道的进来
已经有22人回复
-
这样的filecode谁见过
已经有21人回复
-
filecode与中标关系的预测
已经有5人回复
-
fileCode有新解读？
已经有18人回复
-
关于filecode，很负责任的告诉大家
已经有7人回复
-
【2027博士申请】纳米药物递送方向
已经有5人回复
-
关于豆爷回答的JTJC与%2F数量
已经有7人回复
3楼2014-12-26 10:07:04
acalephs
捐助贵宾 (著名写手)
- 应助: 26 (小学生)
- 金币: 11614.6
- 散金: 181
- 红花: 13
- 帖子: 2562
- 在线: 330.3小时
- 虫号: 1939629
- 注册: 2012-08-15
- 专业: 光学
4楼2015-01-13 20:01:38
