图形图像处理算法(3) ---- 多项式插值法

1. 多项式的定义

多项式（polynomial）是代数学的基本概念，简单的说是由称为未知数的变量和称为系数的常数通过有限次的加减法，乘法，以及自然数幂次的乘方运算得到的代数表达式。多项式是整式的一种。未知数只有一个的多项式称为一元多项式，下面就是一个一元二次多项式

下面所示是一个三元三次多项式：

一个多项式有几次取决于最高次幂的幂指数，有几元取决于未知数的个数，如果某一项不含未知数，称为常数项，多项式在数学的很多分支乃至自然科学中以及工程学中都有很重要的应用。

 2. 多项式插值法

有很多问题都可以按照函数的方式进行描述，然而这个函数通常是未知的，我们只能通过少量的已知的点来推断函数的大致模型，为了实现这个目的，在已知点之间做插值处理

多项式插值的根本点是建立一个特殊形式的多项式，称为插值多项式，多项式是具有如下形式的函数：

 这里的 .....  是系数，当  是非零整数时，这种形式的多项式称为 n 阶多项式，这个是多项式的指数形式，在某些特定的环境中其他形式的多项式则更为简便，牛顿多项式的形式如下：

 这里的 a0..... an 是系数，而 c0..... cn 是中值，注意到   c0..... cn全为0时，牛顿插值多项式就退化为前面的定义的  n 阶多项，对于一个含有 n 个已知点的插值情景，就需要构建 n 次的的多项式，多项式的项数是 n + 1个

 3. 牛顿插值法

牛顿插值法的多项式定义是：

 其中的  是已知点的  坐标， 称为差商，差商的计算方法如下：

也就是说  的差商要由  和  计算出来

 3.1 差商的计算

根据差商的定义，差商需要重复迭代的方法计算出来，下图是 的差商计算过程：

 3.2 C语言实现

牛顿插值法的C语言代码实现如下：

//  x 已知的点的x轴坐标   fx 已知的点y轴左边  n已知点的个数
//  z 插值点的x 轴坐标    pz 输出值 插值点的 y轴坐标  m 插值点的个数
static int interpol(const double *x, const double *fx, int n, double *z, double *pz, int m) {

	double term, *table, *coeff;
	int i, j, k;

	if ((table = (double*)malloc(sizeof(double)*n)) == NULL) {
		return -1;
	}

	if ((coeff = (double*)malloc(sizeof(double)*n)) == NULL) {
		free(table);
		return -1;
	}

	memcpy(table, fx, sizeof(double)*n);

	coeff[0] = table[0];


	//  计算差商的值
	for (k = 1; k < n; k++) {
        //  k = 1 计算第一阶差商的值
        //  k = 2 计算第二阶差商的值 table 被更新
        //  k = n-1  计算最后一阶差商的值,只有 table[0] 被更新
		for (i = 0; i < n - k; i++) {
			j = i + k;
			table[i] = (table[i + 1] - table[i]) / (x[j] - x[i]);
		}

		coeff[k] = table[0];
	}
	free(table);

	//  使用厂商构建插值多项式，以5阶多项式为例
	//  fx = coeff[0] (z-x0)(z-x1)(z-x2)(z-x3)(z-x4) +  coeff[1] (z-x0)(z-x1)(z-x2)(z-x3) + coeff[2](z-x0)(z-x1)(z-x2) +  coeff[3](z-x0)(z-x1) +  coeff[4](z-x0) + pz[0]

	for (k = 0; k < m; k++) {
		pz[k] = coeff[0];
		for (j = 1; j < n; j++) {
				term = coeff[j];
				for (i = 0; i < j; i++) {
					term = term*(z[k] - x[i]);
				}
				pz[k] = pz[k] + term;
		}
	}

	free(coeff);
	return 0;
}

#define INTERPOLSIZE 5

void testinterpol(int argc, char* argv[]) {
    
    //  以正弦为例，1000sinx 的插值
	double inputx[INTERPOLSIZE] = {0.0f, 45.0f, 90.0f, 135.0f, 180.0f};
	double fx[INTERPOLSIZE] = {0.0f, 707.0f, 1000.0f, 707.0f, 0.0f};

	double z[4] = {20.0f, 40.0f, 60.0f,80.0f};

	double pz[4] = { 0.0f };

	interpol(inputx, fx, INTERPOLSIZE, z, pz, 4);

	printf("output value: ");
	for (int i = 0; i < 4; i++) {
		printf("%f ", pz[i]);
	}
}

4. 牛顿插值法的验证

以  为例，插值得到的结果和实际验证的结果误差只有 1 左右：

已知 通过插值得到 的值