通过查找matlab的帮助程序,对离散数据格网化采用的方法有如下5种:
griddata(..., METHOD) where METHOD is one of
'nearest' - Nearest neighbor interpolation
'linear' - Linear interpolation (default)
'natural' - Natural neighbor interpolation
'cubic' - Cubic interpolation (2D only)
'v4' - MATLAB 4 griddata method (2D only)
本文主要对cubic的实现原理进行介绍。
根据matlab的帮助说明:
defines the interpolation method. The 'nearest' and 'linear' methods
have discontinuities in the zero-th and first derivatives respectively,
while the 'cubic' and 'v4' methods produce smooth surfaces. All the
methods except 'v4' are based on a Delaunay triangulation of the data.
意思是说除了‘v4’方法,其他的方法都是基于Delaunary三角形的,所以,cubic方法也不例外。因此matlab中cubic方法就是基于Delaunary的三次方程插值。
步骤:
1、根据离散坐标,参考凸包算法,获得离散点的凸包;
2、构件Delaunary三角网;
3、根据格网化方法,获取格网点的坐标:
可以根据离散点控制的范围,然后纵横坐标的格网间隔进行格网化,也可以根据离散点的控制范围,然后纵横坐标的节点数进行格网化,根据需要来确定。
4、因为要求待插值点在Delaunary三角形内,所以,遍历插值区域内所有的格网点,根据插值点所在的三角形内的三个点,使用三次方程内插出待定点坐标上的值。
注:三次方程插值,并没有找到具体的介绍什么叫三次方程插值,本文使用的方法是根据不在同一条直线上的三个点拟合出一个平面,然后及可以求出平面内任意坐标上点的值。(或者这就叫三次方程插值)
通过和matlab结果进行对比,统计结果如下:
可见插值结果和matlab的结果具有非常高的一致性,统计特征表明,其差值的中误差为0.6mm,精度还是非常高的。
控制区域的对比:
matlab:
本方法:
附:
平面拟合函数:
//点法式方程求平面方程,内插待定点坐标。
void CubicInterpolation(value[] triScatters,ref value target)
{
double x1 = triScatters[0].lon;
double y1 = triScatters[0].lat;
double z1 = triScatters[0].v;
double x2 = triScatters[1].lon;
double y2 = triScatters[1].lat;
double z2 = triScatters[1].v;
double x3 = triScatters[2].lon;
double y3 = triScatters[2].lat;
double z3 = triScatters[2].v;
double x = target.lon;
double y = target.lat;
target.v = (-((y3 - y1) * (z2 - z1) - (z3 - z1) * (y2 - y1)) * (x - x1) - ((z3 - z1) * (x2 - x1) - (x3 - x1) * (z2 - z1)) * (y - y1)) / ((x3 - x1) * (y2 - y1) - (y3 - y1) * (x2 - x1)) + z1;
}
//判断点是否在凸包内
///
/// 判断点是否在凸包内
/// 原理:将待定点‘按照顺序’分别与凸包点形成向量,按顺序两两叉积,最后一个点的向量与第一个点的向量形成叉积
/// 如果点在凸包内,则叉积的符号是相同的,否则,叉积的符号不完全相同。
///
/// 离散的凸包点
/// 待定点x坐标
/// 待定点y坐标
///
bool isInRegion(List boxPoints,double x,double y)
{
double[] cj = new double[boxPoints.Count];
for(int i=0;i
{
double x1=0, y1=0, x2=0, y2=0;
//终点和起点连接
if(i==boxPoints.Count-1)
{
x1 = boxPoints[i].lon - x;
y1 = boxPoints[i].lat - y;
x2 = boxPoints[0].lon - x;
y2 = boxPoints[0].lat - y;
}
else
{
x1 = boxPoints[i].lon - x;
y1 = boxPoints[i].lat - y;
x2 = boxPoints[i + 1].lon - x;
y2 = boxPoints[i + 1].lat - y;
}
double chaji = x1 * y2 - x2 * y1;
cj[i ] = chaji;
//在凸包外
if (i>0&&(cj[i] * cj[i - 1] < 0))
{
return false;
}
}
return true;
}