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

类别:    标签: 数理 matlab   阅读次数:   版权: (CC) BY-NC-SA

对未知曲面进行插值, 在处理数据中经常用到. 我们这里所谓的曲面, 指的是以 $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)));
%}

参考资料

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 【译】蛋白质折叠问题50年的历程
后一篇: 分子尺寸大小的计算

访问人次(2015年7月 9日起): | 最后更新: 2017-08-15 19:57:07 UTC | 版权所有 © 2008 - 2017 Jerkwin