法治建设网站模块,win10运行wordpress,湖南云网站建设,服务佳的广州网站建设基于MATLAB的径向基函数插值#xff08;RBF插值#xff09;#xff08;一维、二维、三维#xff09; 0 前言1 RBF思路2 1维RBF函数2.1 参数说明2.1.1 核函数选择2.1.2 作用半径2.1.3 多项式拟合2.1.4 误差项#xff08;光滑项#xff09; 3 2维RBF函数4 3维RBF函数 惯例声… 基于MATLAB的径向基函数插值RBF插值一维、二维、三维 0 前言1 RBF思路2 1维RBF函数2.1 参数说明2.1.1 核函数选择2.1.2 作用半径2.1.3 多项式拟合2.1.4 误差项光滑项 3 2维RBF函数4 3维RBF函数 惯例声明本人没有相关的工程应用经验只是纯粹对相关算法感兴趣才写此博客。所以如果有错误欢迎在评论区指正不胜感激。本文主要关注于算法的实现对于实际应用等问题本人没有任何经验所以也不再涉及。
0 前言
插值是一个工程中非常常见的扩展数据方法。通常数据测量数量永远是已知的数据的储存空间也是有限的但工程中的数据需求却永远是无限的。
如何用较少位置处的数据点来推广到任意位置处的数据点是工程中常见的问题。其中径向基函数RBF插值具有不依赖数据网格的特点省去了传统插值的网格剖分是一种基于拟合的插值。
本文主要注重于具有几何意义上的RBF插值所以只列举了一维二维和三维插值更高维插值数据可以稍加改写代码就可以实现本文也不再涉及。
本文的参考文献如下 [1]Meshfree Approximation Methods with MATLAB.Gregory E. Fasshauer.
1 RBF思路
径向基函数的大概原理是利用一系列函数叠加对原函数进行拟合。也就是: F ( x ) ∑ w i ∗ f i ( x 0 , x ) F(x)\sum w_i*f_i(x_0,x) F(x)∑wi∗fi(x0,x) 其中f(x0,x)为基函数是中心对称函数函数中心点在x0上。w为每个函数的权重。
因此只需要求出权重w就可以利用基函数在各个点的值计算出整个域的函数。而求解权重也是一个简单的线性代数问题直接用线性方程组求逆的方式就可以得到。
因此RBF方法也具有原理简单编程容易的特点。
下面用一个简单的例子来解释。
首先问题假设如下我有下面6个点的数据xi,yi想插值出蓝色的曲线结果该如何处理 那么第一步我们构造核函数。这里选用高斯函数作为核函数 总共构造6个核函数f(x,xi)每个核函数中心点在xi上。 f i ( x ) e x p ( − ( x − x i ) 2 / 2 / s 2 ) f_i(x)exp(-(x-x_i)^2/2/s^2) fi(x)exp(−(x−xi)2/2/s2) 每个函数乘以权重可以得到当前基函数对应各个点的数值。以第一个基函数为例 w 1 ∗ [ f 1 ( x 1 ) f 1 ( x 2 ) f 1 ( x 3 ) f 1 ( x 4 ) f 1 ( x 5 ) f 1 ( x 6 ) ] w_1* \begin{bmatrix} f_1(x_1)\\ f_1(x_2)\\ f_1(x_3)\\ f_1(x_4)\\ f_1(x_5)\\ f_1(x_6)\\ \end{bmatrix} w1∗ f1(x1)f1(x2)f1(x3)f1(x4)f1(x5)f1(x6) 则原函数F(x)可以计算为 F ( x ) ∑ w i ∗ f i ( x ) F(x)\sum w_i*f_i(x) F(x)∑wi∗fi(x) ∑ w i ∗ [ f i ( x 1 ) f i ( x 2 ) f i ( x 3 ) f i ( x 4 ) f i ( x 5 ) f i ( x 6 ) ] \sum w_i*\begin{bmatrix} f_i(x_1)\\ f_i(x_2)\\ f_i(x_3)\\ f_i(x_4)\\ f_i(x_5)\\ f_i(x_6)\\ \end{bmatrix} ∑wi∗ fi(x1)fi(x2)fi(x3)fi(x4)fi(x5)fi(x6) [ f 1 ( x 1 ) f 2 ( x 1 ) f 3 ( x 1 ) f 4 ( x 1 ) f 5 ( x 1 ) f 6 ( x 1 ) f 1 ( x 2 ) f 2 ( x 2 ) f 3 ( x 2 ) f 4 ( x 2 ) f 5 ( x 2 ) f 6 ( x 2 ) f 1 ( x 3 ) f 2 ( x 3 ) f 3 ( x 3 ) f 4 ( x 3 ) f 5 ( x 3 ) f 6 ( x 3 ) f 1 ( x 4 ) f 2 ( x 4 ) f 3 ( x 4 ) f 4 ( x 4 ) f 5 ( x 4 ) f 6 ( x 4 ) f 1 ( x 5 ) f 2 ( x 5 ) f 3 ( x 5 ) f 4 ( x 5 ) f 5 ( x 5 ) f 6 ( x 5 ) f 1 ( x 6 ) f 2 ( x 6 ) f 3 ( x 6 ) f 4 ( x 6 ) f 5 ( x 6 ) f 6 ( x 6 ) ] ∗ [ w 1 w 2 w 3 w 4 w 5 w 6 ] \begin{bmatrix} f_1(x_1)f_2(x_1)f_3(x_1)f_4(x_1)f_5(x_1)f_6(x_1)\\ f_1(x_2)f_2(x_2)f_3(x_2)f_4(x_2)f_5(x_2)f_6(x_2)\\ f_1(x_3)f_2(x_3)f_3(x_3)f_4(x_3)f_5(x_3)f_6(x_3)\\ f_1(x_4)f_2(x_4)f_3(x_4)f_4(x_4)f_5(x_4)f_6(x_4)\\ f_1(x_5)f_2(x_5)f_3(x_5)f_4(x_5)f_5(x_5)f_6(x_5)\\ f_1(x_6)f_2(x_6)f_3(x_6)f_4(x_6)f_5(x_6)f_6(x_6)\\ \end{bmatrix}* \begin{bmatrix} w_1\\ w_2\\ w_3\\ w_4\\ w_5\\ w_6\\ \end{bmatrix} f1(x1)f1(x2)f1(x3)f1(x4)f1(x5)f1(x6)f2(x1)f2(x2)f2(x3)f2(x4)f2(x5)f2(x6)f3(x1)f3(x2)f3(x3)f3(x4)f3(x5)f3(x6)f4(x1)f4(x2)f4(x3)f4(x4)f4(x5)f4(x6)f5(x1)f5(x2)f5(x3)f5(x4)f5(x5)f5(x6)f6(x1)f6(x2)f6(x3)f6(x4)f6(x5)f6(x6) ∗ w1w2w3w4w5w6 然后利用线性代数方式就可以得到系数w。 则原函数F可以利用这些个基函数叠加得到 叠加后的函数和原函数对比见下图 可以看到计算结果还可以曲线过渡也比较光滑。
如果已知的插值点更多还可以拟合的更好。下图为10个已知点进行插值的结果和预想曲线几乎重合。
上面代码如下
%RBF基本原理
clear
clc
close all%插值点
x00.2:0.5:3;
y0sin(2*pi*0.5*x0)cos(2*pi*0.8*x0)2;
%原函数
x20:0.02:3;
y2sin(2*pi*0.5*x2)cos(2*pi*0.8*x2)2;figure()
hold on
plot(x2,y2)
plot(x0,y0,o,MarkerSize,8)
hold off%1构造函数
N_klength(x0);%构造函数的数量
RBF_Kernelcell(N_k,1);
for k1:N_kx_kx0(k);%中心位置插值节点RBF_Kernel_k(x) exp(-(x-x_k).^2/2/0.3^2);%正态函数RBF_Kernel{k}RBF_Kernel_k;%储存每个函数
end%2计算出插值矩阵
InterpMatzeros(N_k,N_k);
for k1:N_kRBF_phi_kRBF_Kernel{k}(x0);%计算出当前正态函数InterpMat(:,k)RBF_phi_k(:);%插值矩阵每一列储存的都是对应的正态函数
end%3利用插值矩阵求解出每个高斯函数的权重
%∑φ(k)*w(k)InterpMat*wy0可以线性求逆直接得到系数w
wInterpMat\y0;%绘制出每个正态函数
figure()
hold on
mcpcolormap(lines);
for k1:N_kRBF_phi2_kRBF_Kernel{k}(x2)*w(k);plot(x2,RBF_phi2_k,Color,mcp(k,:))plot([x0(k),x0(k)],[0,RBF_Kernel{k}(x0(k))*w(k)],Color,mcp(k,:),LineStyle,--)
end
plot(x2,y2,Color,k,LineWidth,2)
hold off%绘制出最终相加的结果
y_RBFzeros(size(y2));
for k1:N_ky_RBFy_RBFRBF_Kernel{k}(x2)*w(k);%叠加每一个正态函数
end
figure()
hold on
plot(x2,y2)
plot(x2,y_RBF,--)
plot(x0,y0,o,MarkerSize,8)
legend({原函数,RBF插值结果},Location,best)2 1维RBF函数
下面为1维RBF函数插值的代码基本思路和上面第一章节的一样。但是对于计算速度和输入输出等方面简单做了一些优化。
计算结果如下
代码如下
clear
clc
close all%初始已知点
x00:0.1:3;
y0sin(2*pi*0.5*x0)cos(2*pi*0.8*x0)20.02*rand(size(x0));x(-1:0.01:4);
yRBF1(x0,y0,x,gaussian,0.3,0,0.001);figure()
hold on
plot(x0,y0,o)
plot(x,y)
hold offfunction yRBF1(x0,y0,x,Method,Rs,Npoly,Error)
%一维RBF插值输入x0散点y0值。x是输出插值得到的点。
%Method方法默认linear。
%linear,|R|
%gaussian,exp(-(R/Rs)^2)
%thin_plate,R^2*log(R)
%linearEpsR,|R|。分段线性插值要求s更大一些
%cubic,|R^3|%Rs插值核作用半径。Rs对于linearcubicthin_plate无影响。Rs大概和点和点之间的距离差不多就行。
%Npoly多项式拟合。默认是1只拟合1次项。
%Error误差。在(-∞,∞)区间。默认是0无误差。表示可以在部分误差范围内去插值。
%Error一般在0~1之内。如果只是为了收敛可以比较小在0.1以内就可以有很好效果如果想实现平滑可适当增大大于1。%示例x00:0.3:pi;y0sin(2*x0)0.2*x0;x-1:0.01:4;yRBF1(x0,y0,x,linear,0.2,1,0.0);
%示例x00:0.1:3;y00.25*x0.^2-0.5*x00.1*rand(size(x0));x-1:0.01:4;yRBF1(x0,y0,x,gaussian,0.3,2,0.01);Nlength(x0);%点的数量也是RBF核的数量
%整理为列向量
x0x0(:);
y0y0(:);
xx(:);%函数默认信息
if nargin4 || isempty(Method)Methodlinear;%默认线性核函数
elseif nargin5 || isempty(Rs)Rs1.1*(max(x0)-min(x0))/N;%假设空间均匀分布
elseif nargin6 || isempty(Npoly)Npoly1;%默认拟合1次项
elseif nargin7 || isempty(Error)Error0;%默认输入数据没有误差
end%选择核函数
K_method0;
switch Methodcase linearK_method1;fun(RMat) Kernel_Linear(RMat,Rs);case gaussianK_method2;fun(RMat) Kernel_Gaussian(RMat,Rs);Error-Error;%因为零点不是零远点是0所以这里误差为负case thin_plateK_method3;fun(RMat) Kernel_Thin_plate(RMat,Rs);case linearEpsRK_method4;fun(RMat) Kernel_LinearEpsR(RMat,Rs);Error-Error;%因为零点不是零远点是0所以这里误差为负
end%将原始数据分离出多项式项Npoly
if Npoly0Cones(N,1);%常数项wCmean(y0);
elseif Npoly1C[ones(N,1),x0];%常数项一次项wCC\y0;
elseif Npoly2C[ones(N,1),x0,x0.^2];%常数项一次项二次项wCC\y0;
elseerror(只支持Npoly0,1,2)
end
yCC*wC;%多项式项
y1y0-yC;%计算每个基函数在各个点的值
K2zeros(N);%每一列对应一个函数
%计算距离矩阵
DisMatsquareform(pdist(x0));
K3fun(DisMat);%每一个核函数的中心点在节点上
K3K3-Error*eye(N);%把误差函数加入wK3\y1;%计算权重%根据权重计算插值
yzeros(size(x));
for k1:N %计算每个核的贡献然后叠加R_kabs(x-x0(k));ytw(k)*fun(R_k );%plot(t,yt);yyyt;
endNOutlength(y);
%再还原回多项式项
if Npoly0Cones(NOut,1);%常数项
elseif Npoly1C[ones(NOut,1),x];%常数项一次项
elseif Npoly2C[ones(NOut,1),x,x.^2];%常数项一次项二次项
end
yyC*wC;%再加上多项式项%检查结果
if max(abs(w))/max(abs(y1))1e3warning(结果未收敛建议调整误差Error);
elseif (max(y)-min(y))5*(max(y0)-min(y0))warning(结果未收敛建议增大区间Eps或者调整误差Error);
endendfunction zKernel_Linear(R,s)
%R距离
zabs(R);%线性函数
endfunction zKernel_Gaussian(R,s)
%R距离
%s大概和采样点间距差不多就行可以略大更胖
zexp(-R.^2/2/s^2);%正态函数
endfunction zKernel_Thin_plate(R,s)
%R距离
zR.^2.*log(R);%thin_plate
endfunction zKernel_LinearEpsR(R,s)
%R距离
%s大概和采样点间距差不多就行可以略大更胖
s2*s;%这里要更大
zs-R;%线性函数
z(z0)0;
end2.1 参数说明
下面说一下上面函数RBF1的后面各项参数结果
yRBF1(x0,y0,x,Method,Rs,Npoly,Error)2.1.1 核函数选择
RBF函数的核函数有非常多种比如高斯函数、线性函数、薄板函数等。
下面以分别选择高斯函数和线性函数作为核函数的RBF进行对比
clear
clc
close all%初始已知点
x00:0.1:3;
y0sin(2*pi*0.5*x0)cos(2*pi*0.8*x0)20.02*rand(size(x0));x(-1:0.01:4);
yRBF1(x0,y0,x,gaussian,0.3,0,0.001);figure()
hold on
plot(x0,y0,o)
plot(x,y)
hold off可以看到线性核函数的效果相当于线性插值而高斯核函数对峰值的拟合预测很好各有侧重点。
事实上核函数的选择对RBF插值有重要影响不同核函数对应着不同的场景。
2.1.2 作用半径
对于某些核函数还需要确定其作用半径。比如对于高斯函数半径越大核越宽。
对于高斯函数作用半径通常要大于两个散点之间距离但最好不要太大。如果太小会导致每个点都是一个孤立的峰。如果太大会导致求系数时方程求解困难。
下面代码展示了作用半径为0.05、0.3和0.8半径的插值效果。
%初始节点
x00:0.3:3;
y0sin(2*pi*0.5*x0)cos(2*pi*0.8*x0)20.02*rand(size(x0));x(-1:0.01:4);
yRBF1(x0,y0,x,gaussian,0.3,0,0.01);
y3RBF1(x0,y0,x,gaussian,0.05,0,0.0001);
y4RBF1(x0,y0,x,gaussian,0.8,0,0.0001);
figure()
hold on
plot(x,y3)
plot(x,y)
plot(x,y4)
plot(x0,y0,o)
hold off
legend({Rs0.05,Rs0.3,Rs0.7},location,best)
ylim([-0.5,4.5])可以看到上面每个点的间距为0.2。当Rs很小的时候每个点只有一个孤立的高斯峰。当Rs大于间距时比如Rs0.3和0.7都可以求解出结果。
但当Rs比较大比如Rs0.7时会出现计算很奇怪的解比如系数w特别大。这种如果适当的加一小点误差Error有一点就行1e-3量级的已经足够了放宽求解精度可以避免这种系数w很大不收敛的情况。
2.1.3 多项式拟合
本RBF插值函数默认带一个常数项。因为像高斯函数等函数无穷大处是0对于本身数值在0附近的函数插值效果还可以但是函数距离x轴很远时再用这些基函数去拟合效果就会非常差。
因此在RBF插值之前需要先求出常数项C然后将所有散点的yi减去这个常数项C。之后计算完系数再插值后再把这个常数项C加上。
一次项也是同样的原理就是先拟合出y’kxb然后把所有散点yi-y’得到计算系数用的y值。之后再把这个一次项加上即可。
通常常数项一般就够除非数据有很明显的线性度。
下面代码展示了常数项和一次项对比插值的效果
x00:0.3:3;
x(-1:0.01:4);
y0sin(2*pi*0.5*x0)cos(2*pi*0.8*x0)12*x00.02*rand(size(x0));
y1RBF1(x0,y0,x,gaussian,0.3,0,0.0001);
y2RBF1(x0,y0,x,gaussian,0.3,1,0.0001);
figure()
hold on
plot(x,y1)
plot(x,y2)
plot(x0,y0,o)
hold off
legend({常数项,一次项},location,best)2.1.4 误差项光滑项
当数据采集具有一定的误差或者计算不需要完全贴合每一个点的数据或者拟合出的曲线太曲折需要光滑或者计算出的系数w过大甚至不收敛都可以适当的添加误差项来解决。
其中系数w不收敛的问题可能是由于核半径Rs过大或者原始数据本身不太好导致的。适当添加一点误差1e-3甚至更小就可以解决这个问题。
而如果想要消除误差光滑曲线则需要较大的数比如0.2甚至超过1。这个根据具体效果来看。
这个的原理是直接在矩阵K的对角线上减去一个数。这样可以略微放宽求解的约束使得节点处的数值不一定是控制点数值。
下面展示了添加了随机项后大误差和小误差项的对比代码和结果
x00:0.1:3;
y0sin(2*pi*0.5*x0)cos(2*pi*0.8*x0)10.4*randn(size(x0));
x(-0.0:0.01:3.0);
y1RBF1(x0,y0,x,gaussian,0.15,0,0.0001);
y2RBF1(x0,y0,x,gaussian,0.15,0,0.3);figure()
hold on
plot(x,y1)
plot(x,y2)
plot(x0,y0,o)
hold off
legend({小误差项,大误差项},location,best)可以看到小误差项曲线经过了每一个散点。大误差项曲线只是大致经过了散点。
3 2维RBF函数
二维RBF函数的原理和一维相同。
zRBF2(x0,y0,z0,x,y,Method,Rs,Npoly,ExtrapMethod,Error)其中x0、y0是已知散点z0是对应的值。然后要插值出x、y对应的值。
二维RBF函数额外考虑了外插的方法。通常并不建议数据外插因为缺乏足够的信息。
本函数总共考虑了三种外插方法 第一个是’rbf’其实就是RBF方法算出来多少就是多少。 第二个是’nearest’就是不外插超出包线的点的数值直接等于相邻数据点数值。 第三个是’ploy’是多项式外插超出包线的点按照多项式拟合的值计算得到适用于波动不大的数据。 推荐’rbf’方法因为不需要计算包线如果不怎么考虑外插的话非常节省时间。
下面展示了二维RBF插值计算的结果 可以看到RBF插值可以较好的拟合出预期的效果。
展示的代码如下
clear
clc
close all%初始节点
N60;
x0rand(N,1)*6-3;%-3~3
y0rand(N,1)*6-3;%-3~3
z0peaks(x0,y0)0.05*x0;[x,y]meshgrid(-5:0.02:5);
%二维RBF插值
z3RBF2(x0,y0,z0,x,y,gaussian,0.5,1,rbf,0.001);zzeros(size(x));z(:)z3(:);figure()
hold on
surf(x,y,z,EdgeColor,none);
scatter3(x0,y0,z0)
hold offfunction zRBF2(x0,y0,z0,x,y,Method,Rs,Npoly,ExtrapMethod,Error)
%二维RBF插值输入x0y0散点z0值。xy是输出插值得到的点。
%Method方法默认linear。
%linear,|R|
%gaussian,exp(-(R/Rs)^2)
%thin_plate,R^2*log(R)
%linearEpsR,|R|。分段线性插值要求s更大一些
%cubic,|R^3|%Rs插值核作用半径。Rs对于linearcubicthin_plate无影响。Rs大概和点和点之间的距离差不多就行。
%Npoly多项式拟合。默认是1只拟合1次项。
%ExtrapMethod,外插方法。某个数全部赋值为这个数。nearest按照临近值外插。ploy多项式外插。rbf默认用RBF插值得到的值。
%Error误差。在(-∞,∞)区间。默认是0无误差。表示可以在部分误差范围内去插值。
%Error一般在0~1之内。如果只是为了收敛可以比较小在0.1以内就可以有很好效果如果想实现平滑可适当增大大于1。%示例Nnumel(x0);%点的数量也是RBF核的数量
%整理为列向量
x0x0(:);
y0y0(:);
z0z0(:);
xx(:);
yy(:);
%计算包线
[k_ch,V]convhull(x0,y0);
if V0warning(输入散点过于集中);
end
%函数默认信息
narginchk(5,10)
if nargin7 || isempty(Method)Methodlinear;%默认线性核函数
elseif nargin8 || isempty(Rs)Rs1.1*(max(x0)-min(x0))/N;%假设空间均匀分布
elseif nargin9 || isempty(Npoly)Npoly1;%默认拟合1次项
elseif nargin10 || isempty(Npoly)ExtrapMethodrbf;%默认不外插
elseif nargin11 || isempty(Error)Error0;%默认输入数据没有误差
end%选择核函数
switch Methodcase linearfun(RMat) Kernel_Linear(RMat,Rs);case gaussianfun(RMat) Kernel_Gaussian(RMat,Rs);Error-Error;%因为零点不是零远点是0所以这里误差为负case thin_platefun(RMat) Kernel_Thin_plate(RMat,Rs);case linearEpsRfun(RMat) Kernel_LinearEpsR(RMat,Rs);Error-Error;%因为零点不是零远点是0所以这里误差为负
end%将原始数据分离出多项式项Npoly
if Npoly0Cones(N,1);%常数项wCmean(z0);
elseif Npoly1if N2;warning(点数太少建议Npoly等于0);endC[ones(N,1),x0,y0];%常数项一次项wCC\z0;
elseif Npoly2if N5;warning(点数太少建议Npoly等于1);endC[ones(N,1),x0,y0,x0.^2,y0.^2,x0.*y0];%常数项一次项二次项wCC\z0;
elseerror(只支持Npoly0,1,2)
end
zCC*wC;%多项式项
z1z0-zC;%计算距离矩阵
DisMatsquareform(pdist([x0,y0]));
K3fun(DisMat);%每一个核函数的中心点在节点上
K3K3-Error*eye(N);%把误差函数加入wK3\z1;%计算权重%根据权重计算插值
zzeros(size(x));
for k1:N %计算每个核的贡献然后叠加R_ksqrt((x-x0(k)).^2(y-y0(k)).^2);ztw(k)*fun(R_k );zzzt;
endNOutlength(z);
%再还原回多项式项
if Npoly0Cones(NOut,1);%常数项
elseif Npoly1C[ones(NOut,1),x,y];%常数项一次项
elseif Npoly2C[ones(NOut,1),x,y,x.^2,y.^2,x.*y];%常数项一次项二次项
end
zC2C*wC;%再加上多项式项%判断外插方法
Inpolyinpolygon(x,y,x0(k_ch),y0(k_ch));%多边形内
indxInpolyfind(Inpoly);
if isnumeric(ExtrapMethod)z(Inpoly)z(Inpoly)zC1(Inpoly);%内部的正常z(~Inpoly)ExtrapMethod;%外部的固定值
elseif ischar(ExtrapMethod)switch ExtrapMethodcase rbfzzzC2;%再加上多项式项case ployz(Inpoly)z(Inpoly)zC2(Inpoly);%内部的正常z(~Inpoly)zC2(~Inpoly);%外部的直接用多项式值case nearestz(Inpoly)z(Inpoly)zC2(Inpoly);%内部的正常%找到最近的点F scatteredInterpolant(x(Inpoly),y(Inpoly),z(Inpoly),nearest,nearest);%外部的直接插值z(~Inpoly)F(x(~Inpoly),y(~Inpoly));end
elseerror(未识别ExtrapMethod)
end%检查结果
if max(abs(w))/max(abs(z1))1e3warning(结果未收敛建议调整误差Error);
elseif (max(z)-min(z))5*(max(z0)-min(z0))warning(结果未收敛建议增大区间Eps或者调整误差Error);
endendfunction zKernel_Linear(R,s)
%R距离
zabs(R);%线性函数
endfunction zKernel_Gaussian(R,s)
%R距离
%s大概和采样点间距差不多就行可以略大更胖
zexp(-R.^2/2/s^2);%正态函数
endfunction zKernel_Thin_plate(R,s)
%R距离
zR.^2.*log(R);%thin_plate
endfunction zKernel_LinearEpsR(R,s)
%R距离
%s大概和采样点间距差不多就行可以略大更胖
s2*s;%这里要更大
zs-R;%线性函数
z(z0)0;
end4 3维RBF函数
计算方法和二维的方法一致。其实更高维的数据计算方法也相差不多。 对于实际工程的几何插值使用最高维度也就是三维。对于其它应用可能还是需要高维的插值本文暂不涉及由于RBF代码较为简单因此稍加改动即可更改为高维RBF插值结果。
三维RBF函数的原理和二维也相同。
PRBF3(x0,y0,z0,P0,x,y,z,Method,Rs,Npoly,ExtrapMethod,Error)其中x0、y0、z0是已知散点P0是对应的值。然后要插值出x、y、z对应的值。
下面展示了三维RBF插值计算的结果和对应代码 clear
clc
close all%初始节点
N1000;
x0rand(N,1)*6-3;%-3~3
y0rand(N,1)*6-3;%-3~3
z0rand(N,1)*6-3;%-3~3
P0sin(1*pi*x0)cos(1.5*pi*y0)sin(1*pi*z0)1;[x,y,z]meshgrid(-5:0.2:5);P3RBF3(x0,y0,z0,P0,x,y,z,linear,0.5,1,nearest,0.01);%RBF插值Pzeros(size(x));P(:)P3(:);%实际插值结果
figure()
hold on
sslice(x,y,z,P,[0],[0],[0]);
set(s,LineStyle,none);
scatter3(x0,y0,z0)
hold off
xlabel(x);ylabel(y);zlabel(z);
view(3)%理论目标结果
figure()
P3sin(1*pi*x)cos(1.5*pi*y)sin(1*pi*z)1;
P3(x3)0;P3(y3)0;P3(z3)0;P3(x-3)0;P3(y-3)0;P3(z-3)0;
sslice(x,y,z,P3,[0],[0],[0]);
xlabel(x);ylabel(y);zlabel(z);
view(3)function PRBF3(x0,y0,z0,P0,x,y,z,Method,Rs,Npoly,ExtrapMethod,Error)
%二维RBF插值输入x0y0散点z0值。xy是输出插值得到的点。
%Method方法默认linear。
%linear,|R|
%gaussian,exp(-(R/Rs)^2)
%thin_plate,R^2*log(R)
%linearEpsR,|R|。分段线性插值要求s更大一些
%cubic,|R^3|%Rs插值核作用半径。Rs对于linearcubicthin_plate无影响。Rs大概和点和点之间的距离差不多就行。
%Npoly多项式拟合。默认是1只拟合1次项。
%ExtrapMethod,外插方法。某个数全部赋值为这个数。nearest按照临近值外插。rbf默认用RBF插值得到的值。三维外插判据复杂不建议用。
%Error误差。在(-∞,∞)区间。默认是0无误差。表示可以在部分误差范围内去插值。
%Error一般在0~1之内。如果只是为了收敛可以比较小在0.1以内就可以有很好效果如果想实现平滑可适当增大大于1。
Nnumel(x0);%点的数量也是RBF核的数量
%整理为列向量
x0x0(:);
y0y0(:);
z0z0(:);
P0P0(:);
xx(:);
yy(:);
zz(:);
%计算包线
[k_ch,V]convhull(x0,y0,z0);
%函数默认信息
narginchk(7,12)
if nargin7 || isempty(Method)Methodlinear;%默认线性核函数
elseif nargin8 || isempty(Rs)Rs1.1*(max(x0)-min(x0))/N;%假设空间均匀分布
elseif nargin9 || isempty(Npoly)Npoly1;%默认拟合1次项
elseif nargin10 || isempty(Npoly)ExtrapMethodrbf;%默认不外插
elseif nargin11 || isempty(Error)Error0;%默认输入数据没有误差
end%选择外插包线
if strcmp(ExtrapMethod,rbf)Inpoly[];
elsek_chunique(k_ch);InpolyIsPointInShape3([x0(k_ch),y0(k_ch),z0(k_ch)],[x,y,z]);
end
if V0warning(输入散点过于集中);
end
%选择核函数
switch Methodcase linearfun(RMat) Kernel_Linear(RMat,Rs);case gaussianfun(RMat) Kernel_Gaussian(RMat,Rs);Error-Error;%因为零点不是零远点是0所以这里误差为负case thin_platefun(RMat) Kernel_Thin_plate(RMat,Rs);case linearEpsRfun(RMat) Kernel_LinearEpsR(RMat,Rs);Error-Error;%因为零点不是零远点是0所以这里误差为负
end
%将原始数据分离出多项式项Npoly
if Npoly0Cones(N,1);%常数项wCmean(P0);
elseif Npoly1if N3;warning(点数太少建议Npoly等于0);endC[ones(N,1),x0,y0,z0];%常数项一次项wCC\P0;
elseif Npoly2if N9;warning(点数太少建议Npoly等于1);endC[ones(N,1),x0,y0,z0,x0.^2,y0.^2,z0.^2,x0.*y0,x0.*z0,z0.*y0];%常数项一次项二次项wCC\P0;
elseerror(只支持Npoly0,1,2)
end
PCC*wC;%多项式项
P1P0-PC;
%计算距离矩阵
DisMatsquareform(pdist([x0,y0,z0]));
K3fun(DisMat);%每一个核函数的中心点在节点上
K3K3-Error*eye(N);%把误差函数加入wK3\P1;%计算权重
%根据权重计算插值
Pzeros(size(x));
for k1:N %计算每个核的贡献然后叠加R_ksqrt((x-x0(k)).^2(y-y0(k)).^2(z-z0(k)).^2);Ptw(k)*fun(R_k );PPPt;
end
NOutlength(P);
%再还原回多项式项
if Npoly0Cones(NOut,1);%常数项
elseif Npoly1C[ones(NOut,1),x,y,z];%常数项一次项
elseif Npoly2C[ones(NOut,1),x,y,z,x.^2,y.^2,z.^2,x.*y,x.*z,z.*y];%常数项一次项二次项
end
PC2C*wC;%再加上多项式项%判断外插方法
if isnumeric(ExtrapMethod)P(Inpoly)P(Inpoly)zC1(Inpoly);%内部的正常P(~Inpoly)ExtrapMethod;%外部的固定值
elseif ischar(ExtrapMethod)switch ExtrapMethodcase rbfPPPC2;%再加上多项式项case ployP(Inpoly)P(Inpoly)PC2(Inpoly);%内部的正常P(~Inpoly)PC2(~Inpoly);%外部的直接用多项式值case nearestP(Inpoly)P(Inpoly)PC2(Inpoly);%内部的正常%找到最近的点F scatteredInterpolant(x(Inpoly),y(Inpoly),z(Inpoly),P(Inpoly),nearest,nearest);%外部的直接插值P(~Inpoly)F(x(~Inpoly),y(~Inpoly),z(~Inpoly));end
elseerror(未识别ExtrapMethod)
end
%检查结果
if max(abs(w))/max(abs(P1))1e3warning(结果未收敛建议调整误差Error);
elseif (max(P)-min(P))5*(max(P0)-min(P0))warning(结果未收敛建议增大区间Eps或者调整误差Error);
end
endfunction zKernel_Linear(R,s)
%R距离
zabs(R);%线性函数
endfunction zKernel_Gaussian(R,s)
%R距离
%s大概和采样点间距差不多就行可以略大更胖
zexp(-R.^2/2/s^2);%正态函数
endfunction zKernel_Thin_plate(R,s)
%R距离
zR.^2.*log(R);%thin_plate
endfunction zKernel_LinearEpsR(R,s)
%R距离
%s大概和采样点间距差不多就行可以略大更胖
s2*s;%这里要更大
zs-R;%线性函数
z(z0)0;
endfunction IsInShapeIsPointInShape3(A,points)
%判断点是否在某个三维凸多面体内
%A凸多面体的顶点。pointsN行3列。
%例子A[0,0,0;0,0,1;0,1,0;0,1,1;1,0,0;1,0,1;1,1,0;1,1,1];pointsrand(1000,3)*3-1;
%IsInShapeIsPointInShape3(A,points)NPointssize(points,1);%计算有多少个点
IsInShapefalse(NPoints,1);
NConvhullsize(A,1);%计算凸包上有多少个点
BD_A[min(A,[],1);max(A,[],1)];
for kp1:NPointsP_kpoints(kp,:);if any(P_kBD_A(1,:)) || any(P_kBD_A(2,:)) %如果超过A的矩形边界说明肯定不在A内continueendNewA[A;P_k];%把这个点加入新凸包计算NewPmax(unique(convhull(NewA,Simplify,false)));%查看新凸包这个点会不会涉及到新点IsInShape(kp)(NewPNConvhull);
end
end
