Matlab插值曲面及其曲率的计算|Jerkwin

2016-06-20 10:41:02

对未知曲面进行插值, 在处理数据中经常用到. 我们这里所谓的曲面, 指的是以 $x, y, z$ 坐标给出的一些点. 这些点决定了一个单值曲面 $z=f(x,y)$, 但我们并不知道 $f(x,y)$ 的具体形式. 否则的话, 我们就可以进行拟合, 而无须插值了.

给出的格点数据可分为两类. 一类是规则数据, 也称均匀数据, 就是按一定步长对X-Y平面进行分格, 给出 $N_x \times N_y$ 个 $x,y, z$ 数据. 对这种数据, 插值相对容易, 样条函数法, 卷积法效果都不错. 我也写过一段代码进行这种规则数据的三次卷积插值. 如果给出的 $x, y,z$ 数据是随意的, 并不遵循某种规则, 这种无规则数据也称为非均匀数据. 非均匀数据的插值相对较困难, 不同插值的方法可以给出大不相同的曲面.

对非均匀曲面的插值, 常用的方法有样条函数法, 径向基函数法, Matlab的v4方法, 还有高斯过程法(Kriging为其中一种). 具体的理论还是很复杂的, 这里就不细究了, 只关注如何使用Matlab对曲面进行插值.

Matlab自带的曲面插值函数为griddata, 但这个函数效果不好, 所以有人开发了新函数gridfit, 很多人使用这个函数. 又有人在gridfit基础之上开发了RegularizeData3D, 效果更好, 这是目前最好的曲面插值函数, 建议优先使用.

除了使用这两个函数进行曲面插值之外, 也还有其他的一些方法. 如有人提到

先用tri = delaunay(x,y)对区域进行三角剖分, 让点自行连接成一个个三角形, 然后使用trisurf(tri,x,y,z)生成曲面, 再用shading interp插值拟合. 注意, 如果你的曲面在xy平面的投影不是矩形, 要用inpolygon把不在区域内的点删掉.

这种方法可行, 但稍嫌麻烦, 效果也未必好, 仅供参考.

对于Kriging方法, 效果也不错, 但比较复杂, 控制参数很多, 要知晓很多知识才能熟练使用, 不太适合用于简单的插值. 如果需要使用, 请参阅下面参考资料中的相关链接.

下面是使用griddata以及RegularizeData3D函数进行曲面插值并计算插值曲面曲率的示例代码, 供参考. 点击这里下载相关的文件.

matlab
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81	%% 根据离散点坐标(x,y,z), 插值生成空间曲面 clc; clear; clear all; % 载入数据文件 % 数据文件每行三个数值分别为x, y, z坐标, 各数值之间以空格作为分隔符 % load xyz.dat % 或使用下面的测试数据，8个离散点的xyz坐标 xyz=[ 1 0.5 0.5 1 2 1 1 3.5 0.5 3 0.75 0.5 3 1.5 1 3 2.25 0.5 3 3 1 3 4 0.5 ]; % 获取xyz数据xyz坐标 x=xyz(:,1); y=xyz(:,2); z=xyz(:,3); % 生成规则网格坐标X和Y dX=0.1; dY=dX; % X和Y方向步长均取0.1 Xmin=min(x)-dX; Ymin=min(y)-dY; % X和Y方向范围 Xmax=max(x)+dX; Ymax=max(y)+dY; % 对网格进行插值生成相应的z坐标 % 注意: 不同插值方法得到的曲面不同 % 内置函数, 建议优先使用v4 [X,Y]=meshgrid(Xmin:dX:Xmax, Ymin:dY:Ymax); Zv4=griddata(x,y,z, X,Y,'v4'); % 外部函数, 优先使用bicubic方法, Bilinear得到的结果可能不正确 X=Xmin:dX:Xmax; Y=Ymin:dY:Ymax; Smoothness=0.0001; % 控制曲面光滑程度, 越小越接近数据点 Zreg=RegularizeData3D(x,y,z, X,Y, 'interp', 'bicubic', 'smoothness', Smoothness); % 绘制离散点及插值曲面 figure(1) plot3(x,y,z, 'r.','MarkerSize',30) hold on; grid on mesh(X,Y,Zv4, 'facealpha',0) surf(X,Y,Zreg, 'facealpha',.75, 'FaceColor','interp', 'EdgeColor','none') title({'\fontname{微软雅黑}不同插值方法得到的曲面'; ... '线框: griddata-v4 表面: RegularizeData3D-bicubic'}) %% 根据插值曲面计算曲面曲率 figure(2) [X,Y]=meshgrid(Xmin:dX:Xmax, Ymin:dY:Ymax); [K,H,P1,P2] = surfature(X,Y,Zreg); surf(X,Y,Zreg,H,'facecolor','interp'); title '\fontname{微软雅黑}RegularizeData3D-bicubic插值曲面的曲率' %% 径向基Multiquadric插值方法 % 结果与RegularizeData3D-bicubic相差不大, 可用于研究具体算法的实现 %{ [xj,xi]=meshgrid(x); [yj,yi]=meshgrid(y); dij2=(xi-xj).^2+(yi-yj).^2; idls=logical(tril(ones(length(x)),-1)); idx=dij2(idls)>0; if all(~idx) delta=1/618; else delta=1/max(sqrt(dij2(idx))); end Q=sqrt(dij2+delta^2); alpha=Q\z; [X,Y]=meshgrid(Xmin:dX:Xmax, Ymin:dY:Ymax); [Xj,Xi]=meshgrid(x,X); [Yj,Yi]=meshgrid(y,Y); Z=sqrt((Xi-Xj).^2+(Yi-Yj).^2+delta^2)*alpha; hold on surf(X,Y,reshape(Z,size(X))); %}

matlab

%% 根据离散点坐标(x,y,z), 插值生成空间曲面
clc; clear; clear all;

% 载入数据文件
% 数据文件每行三个数值分别为x, y, z坐标, 各数值之间以空格作为分隔符
% load xyz.dat

% 或使用下面的测试数据，8个离散点的xyz坐标
xyz=[
	1  0.5   0.5
	1  2     1
	1  3.5   0.5
	3  0.75  0.5
	3  1.5   1
	3  2.25  0.5
	3  3     1
	3  4     0.5
];

% 获取xyz数据xyz坐标
x=xyz(:,1); y=xyz(:,2); z=xyz(:,3);

% 生成规则网格坐标X和Y
dX=0.1;         dY=dX;          % X和Y方向步长均取0.1
Xmin=min(x)-dX; Ymin=min(y)-dY; % X和Y方向范围
Xmax=max(x)+dX; Ymax=max(y)+dY;

% 对网格进行插值生成相应的z坐标
% 注意: 不同插值方法得到的曲面不同

% 内置函数, 建议优先使用v4
[X,Y]=meshgrid(Xmin:dX:Xmax, Ymin:dY:Ymax);
Zv4=griddata(x,y,z, X,Y,'v4');

% 外部函数, 优先使用bicubic方法, Bilinear得到的结果可能不正确
X=Xmin:dX:Xmax;
Y=Ymin:dY:Ymax;
Smoothness=0.0001; % 控制曲面光滑程度, 越小越接近数据点
Zreg=RegularizeData3D(x,y,z, X,Y, 'interp', 'bicubic', 'smoothness', Smoothness);

% 绘制离散点及插值曲面
figure(1)
plot3(x,y,z, 'r.','MarkerSize',30)
hold on; grid on
mesh(X,Y,Zv4,  'facealpha',0)
surf(X,Y,Zreg, 'facealpha',.75, 'FaceColor','interp', 'EdgeColor','none')

title({'\fontname{微软雅黑}不同插值方法得到的曲面'; ...
    '线框: griddata-v4    表面: RegularizeData3D-bicubic'})

%% 根据插值曲面计算曲面曲率
figure(2)
[X,Y]=meshgrid(Xmin:dX:Xmax, Ymin:dY:Ymax);
[K,H,P1,P2] = surfature(X,Y,Zreg);
surf(X,Y,Zreg,H,'facecolor','interp');
title '\fontname{微软雅黑}RegularizeData3D-bicubic插值曲面的曲率'

%% 径向基Multiquadric插值方法
% 结果与RegularizeData3D-bicubic相差不大, 可用于研究具体算法的实现
%{
[xj,xi]=meshgrid(x);
[yj,yi]=meshgrid(y);
dij2=(xi-xj).^2+(yi-yj).^2;
idls=logical(tril(ones(length(x)),-1));
idx=dij2(idls)>0;
if all(~idx)
    delta=1/618;
else
    delta=1/max(sqrt(dij2(idx)));
end
Q=sqrt(dij2+delta^2);
alpha=Q\z;

[X,Y]=meshgrid(Xmin:dX:Xmax, Ymin:dY:Ymax);
[Xj,Xi]=meshgrid(x,X);
[Yj,Yi]=meshgrid(y,Y);
Z=sqrt((Xi-Xj).^2+(Yi-Yj).^2+delta^2)*alpha;

hold on
surf(X,Y,reshape(Z,size(X)));
%}

参考资料