二维三次卷积插值算法

      声明：本博客为原创博客，未经允许，不得转载！小伙伴们如果是在别的地方看到的话，建议还是来csdn上看吧（原文链接地址为http://blog.csdn.net/bettarwang/article/details/11376241)，毕竟在csdn上看代码和提问、讨论都更方便。

     生产或生活中，我们往往只能得到有限的关于某个物理量的离散点（比如一张图片，我们只能得到其原大小内每个坐标处的颜色信息），但是很多时候，我们还需要不在已有栅格点的物理量信息（比如图像缩放时）。此时，就需要某种估计方法，用于得到该坐标处的物理量值，此时就需要通过插值得到。常用的插值算法有很多，从简单的线性插值，Lagrange插值，Newton插值，到较复杂的Herminate插值，二维三次卷积插值，分段多项式插值以及样条插值。

     那么在使用时，该如何选择呢？

    一般来说，有两个准则：第一个是要满足精度要求、连续性要求（比如曲面过渡往往要求一阶或者二阶导数连续);第二个是要满足具体使用时的算法复杂度要求。算法本无优劣，须要结合具体情况灵活选择。

    这里介绍一种在图像处理和DEM(Digital Elevation Model)中应用很广的一种插值算法----二维三次卷积插值算法。

    DEM数据是二维数据，对于二维数据，经常使用的插值方法有双线性插值和双立方插值，它们都是在最邻近的四个网点之间进行插值。双线性插值算法简单，但是却无法得到网点处的导数.双立方Herminate插值是四个邻近网点高度的加权和，且方程连接可微。但是该方法的缺点是网点处的微分被强迫为零，虽保证了插值方程的连续性，但是地形在网点处为平面时和直觉矛盾。

    综合考虑算法效率和插值精度，二维三次卷积插值是比较理想的插值方法。

   在处理地形时， 二维三次插值算法是以待插值点周围16个点的高程值来得到它自己的高程值。如下图所示。

     

 

由于该插值方法需要使用未知点周围16个点的数据，当数字地形风格为有限大小时，如果未知点处于风格国家级，就无法得到风格以下的网点数据，因此，需要研究该算法的边界条件。

     Robert G K在《Cubic convolution interpolation for digital image processing》一文中给出了一种较好的边界条件。设

网点在x方向由0到M，在y方向由0到N，则边界条件为

下面是用C++写的二维三次卷积插值的源代码。

     首先是point.h文件

#ifndef POINT_H
#define POINT_H
typedef struct 
{
    float x;
    float y;
    float z;
}Point;
 
#endif

然后是interpolation.h文件

#ifndef INTERPOLATION_h
#define INTERPOLATION_H
#include "Point.h"
#include <string>
#include <Math.h>
#include <iostream>
#include <fstream>
using namespace std;
float getInterpolationHeight(float**,int,int,float,float);
float funA(float);
float funB(float);
float funC(float);
float funD(float);
float getDEMInterpolationHeight(Point**p,float**expandedH,int rowNum,int colNum,fl