博文

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

已有 16767 次阅读 2016-6-21 00:21 |系统分类:科研笔记

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=[ 10.50.512113.50.530.750.531.5132.250.5331340.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)));endQ=sqrt(dij2+delta^2);alpha=Qz;[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 onsurf(X,Y,reshape(Z,size(X)));%}

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=[
	10.50.512113.50.530.750.531.5132.250.5331340.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)));endQ=sqrt(dij2+delta^2);alpha=Qz;[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 onsurf(X,Y,reshape(Z,size(X)));%}