一般而言, 二元函数或者曲面面的插值maltab已经有很强大的库函数了, 例如三次样条插值spline, 曲面插值interp, interp2, griddata等等. 但是对于一组三维曲线的拟合似乎并没有很好的函数可以解决这个问题, 本文简单介绍一下我的方法.

我们知道, 插值实际上就是通过某种插值方式, 例如拉格朗日插值, 牛顿插值, 三次样条插值等等方式对 邻近 的几个点插入一个多项式函数, 此即为插值. 对于二维曲线而言, 在数据集的基础上, 我们选择的插值的数据在区间[min(x), max(x)]即可, 但是对于三维曲线而言就没有简单[min(x), max(x)]或者[min(y), max(y)]. 我的思路就是, 从曲线的一个端点开始, 寻找其在三维空间中的最近邻, 作为其近邻位置, 然后在这两点间线性插值或者其他方法进行插值, 知道行走到曲线的另一端为止. 这就实现了三维曲线的空间插点的操作. 下面给出代码, 代码中会给出重要的注释.

  • 读入数据, 设置必要参数 (这里没有考虑到数据的稀疏或者密集来自动选取插值的个数, 后面还会更新的…)
clear
clc
data = load('XYZ'); %load data
index = 10; % index is  the posiiton of origin point
node = 50; % 任意两个点插入的点的个数
  • 求数据的距离矩阵, 以便接下来寻找数据点的近邻
d = zeros(size(data,1));
for ii = 1:size(data,1)
    for jj = ii+1:size(data,1)
        d(ii,jj) = norm(data(ii,:) - data(jj,:));
    end
end
d = d + d';
  • 利用上面的距离矩阵, 寻找近邻关系, 线性插入数据点

注意到这里需要做一个atom_list的表格, 用于记录已经插值过的数据点, 防止寻找数据近邻的时候重复寻找.

node = 50;
X = zeros(size(data,1)*node-node,1);
Y = zeros(size(data,1)*node-node,1);
[minimum, index] = min(data(:,3));
ind_G = index;
G = data(ind_G, 1:2);
atom_list = ind_G;

for ii = 1:size(data,1)-1
    [~, q] = sort(d(ind_G,:));
    [~,qq] = ismember(atom_list, q);q(qq) = [];
    M = data(q(1),1:2);
    G_M = M - G;
    split_point = repmat(G,node+1,1) +  bsxfun(@times, linspace(0,1,node+1)', repmat(G_M,node+1,1));
    X(node*(ii-1)+1:node*ii,1) = split_point(1:end-1,1);
    Y(node*(ii-1)+1:node*ii,1) = split_point(1:end-1,2);
    G = M;ind_G = q(1);
    atom_list = [atom_list ind_G];
end

  • 使用griddata进行插值, 并画图比较
Z = griddata(data(:,1), data(:,2), data(:,3), X, Y, 'natural');

figure
subplot(1,2,1)
scatter3(data(:,1),data(:,2),data(:,3),50,'filled');
subplot(1,2,2)
scatter3(X,Y,Z,50,'filled');

下图便是一个效果图: