通过查找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;

}