声明:本博客为原创博客,未经允许,不得转载!小伙伴们如果是在别的地方看到的话,建议还是来csdn上看吧(原文链接地址为http://blog.youkuaiyun.com/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,float delta,float x,float y);
float**getExpandedH(