matlab中二维插值中cubic方法的实现原理（个人见解）

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

        }

//判断点是否在凸包内

/// <summary>
        /// 判断点是否在凸包内
        /// 原理：将待定点‘按照顺序’分别与凸包点形成向量，按顺序两两叉积，最后一个点的向量与第一个点的向量形成叉积
        ///      如果点在凸包内，则叉积的符号是相同的，否则，叉积的符号不完全相同。
        /// </summary>
        /// <param name="boxPoints">离散的凸包点</param>
        /// <param name="x">待定点x坐标</param>
        /// <param name="y">待定点y坐标</param>
        /// <returns></returns>
        bool isInRegion(List<value> boxPoints,double x,double y)
        {
            double[] cj = new double[boxPoints.Count];
            for(int i=0;i<boxPoints.Count;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;
        }