matlab中曲面的实际问题,科学网—Matlab插值曲面及其曲率的计算

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

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;elsedelta=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)));%}