基于径向基函数数据插值

径向基函数的定义：一个取值仅仅依赖于离支撑点距离的函数

径向基函数插值具有很多优点：比如径向基函数形式简单，插值只需用到每个网格节点的坐标，与网格之间是怎么连接的无关，数据结构简单  ，编程较容易实现。

径向基函数插值：利用已知有限个离散点坐标及已知的插值数据信息（比如：网格变形的变形量，数据传递的物理量） 构造一个只与这些离散点距离相关的插值系统。

径向基函数插值系统表达式

R表示支撑半径，插值系统里，未知量只有各个点的权重αi

而权重可以根据已知点的坐标及已知的插值数据信息代入插值系统里面，最后将待插值的坐标再一次代入插值系统，由于此时权重已经可以计算得到，即可求得待插值点的数据信息。

下面以网格变形为例：

#include <iostream>
#include <vector>
#include <algorithm>
#include<string>
#include<sstream>
#include<cmath>
#include <D:\\eigen-3.4.0\\Eigen/LU> 

#define PI acos(-1)

using namespace std;
using namespace Eigen;
int _max_number;//网格节点数
double R;//物面最大距离

class Point
{
public: double x;
public: double y;
public: double z;
};

//径向基函数表达式:I(r)=权重(w)*系数矩阵(M)

Eigen::Vector<double, Dynamic>rbf_values;//变形量，rbf函数右边的值，相当于I（r)
Eigen::Matrix<double, Dynamic, Dynamic> rbf_system_mat;//rbf里面的系数矩阵，相当于M
Eigen::Matrix<double, Dynamic, 3> rbf_weights;//权重，相当于w



//计算距离
double distance(Point P1, Point P2)
{
	return sqrt((P1.x - P2.x) * (P1.x - P2.x) + (P1.y - P2.y) * (P1.y - P2.y) + (P1.z - P2.z) * (P1.z - P2.z));
}





//定义径向基函数
double fun(Point P1, Point P2)
{
	double f;
	double d = distance(P1, P2)/R;
	double dd = (1 - d);
	dd *= dd;
	dd *= dd;
	f = dd * (4 * d + 1);
	return f;
}


std::vector<Point>rbf_point;//保存点的坐标
void RBF_mesh_deformation()
{
	rbf_values.resize(_max_number);
	rbf_system_mat.resize(_max_number, _max_number);
	double original_coord_x = 0.0;double original_coord_y = 0.0;double original_coord_z = 0.0;//原网格的坐标
	double rbf_x = 0.0;double rbf_y = 0.0;double rbf_z = 0.0;//插值得到的坐标
	double final_x = 0.0;double final_y = 0.0;double final_z = 0.0;//插值后的坐标
	//求系数矩阵
	for (int id = 0; id < _max_number; id++)
	{
		for (int jd = 0; jd < _max_number; jd++)
		{
			rbf_system_mat(id,jd)= fun(rbf_point[id], rbf_point[jd]);
		}
	}
	//给定变形量
	for (int id = 0; id < _max_number; id++)
	{
		rbf_values[id] = sin(2 * PI * rbf_point[id].x);//给定正弦函数的变形量
	}
	//求权重
	rbf_weights = rbf_system_mat.lu().solve(rbf_values);



	//rbf插值
	double def_coord_x = 0.0;double def_coord_y = 0.0;double def_coord_z = 0.0;
	double f2 = 0.0;
	for (int id = 0; id < _max_number; id++)
	{
		def_coord_x = 0;
		def_coord_y = 0;
		def_coord_z = 0;
		for (int jd = 0; jd < _max_number; jd++)
		{
			f2 = fun(rbf_point[id], rbf_point[jd]);
			def_coord_x += rbf_weights(jd, 0) * f2;//x坐标
			def_coord_y += rbf_weights(jd, 1) * f2;//y坐标
			def_coord_z += rbf_weights(jd, 2) * f2;//z坐标
		}
		original_coord_x = rbf_point[id].x;
		original_coord_y = rbf_point[id].y;
		original_coord_z = rbf_point[id].z;
		final_x = original_coord_x + def_coord_x;
		final_y = original_coord_y + def_coord_y;
		final_z = original_coord_z + def_coord_z;
		rbf_point[id].x = final_x;//更新坐标
		rbf_point[id].y = final_y;//更新坐标
		rbf_point[id].z = final_z;//更新坐标
	}
}

上面代码中求解线性方程组时调用了Eigen库