三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)

三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)

样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法,三次样条又是其中用的较为广泛的一种。本篇介绍力求用容易理解的方式,介绍一下三次样条插值的原理,并附C语言的实现代码。

1. 三次样条曲线原理


假设有以下节点



1.1 定义


样条曲线 是一个分段定义的公式。给定n+1个数据点,共有n个区间,三次样条方程满足以下条件:

a. 在每个分段区间 (i = 0, 1, …, n-1,x递增), 都是一个三次多项式。

b. 满足 (i = 0, 1, …, n )

c. ,导数 ,二阶导数 在[a, b]区间都是连续的,即曲线是光滑的。

所以n个三次多项式分段可以写作:

i = 0, 1, …, n-1

其中ai, bi, ci, di代表4n个未知系数。

1.2 求解


已知:

a. n+1个数据点[xi, yi], i = 0, 1, …, n

b. 每一分段都是三次多项式函数曲线

c. 节点达到二阶连续

d. 左右两端点处特性(自然边界,固定边界,非节点边界)

根据定点,求出每段样条曲线方程中的系数,即可得到每段曲线的具体表达式。

插值和连续性:

, 其中 i = 0, 1, …, n-1

微分连续性:

, 其中 i = 0, 1, …, n-2

样条曲线的微分式:


将步长 带入样条曲线的条件:

a. 由 (i = 0, 1, …, n-1)推出

b. 由 (i = 0, 1, …, n-1)推出

c. 由 (i = 0, 1, …, n-2)推出



由此可得:

d. 由 (i = 0, 1, …, n-2)推出

设 ,则

a. 可写为:

,推出

b. 将ci, di带入 可得:

c. 将bi, ci, di带入 (i = 0, 1, …, n-2)可得:


端点条件

由i的取值范围可知,共有n-1个公式, 但却有n+1个未知量m 。要想求解该方程组,还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。 选择不是唯一的,3种比较常用的限制如下。

a. 自由边界(Natural)

首尾两端没有受到任何让它们弯曲的力,即 。具体表示为 和

则要求解的方程组可写为:


b. 固定边界(Clamped)

首尾两端点的微分值是被指定的,这里分别定为A和B。则可以推出





将上述两个公式带入方程组,新的方程组左侧为



c. 非节点边界(Not-A-Knot)

指定样条曲线的三次微分匹配,即

根据 和 ,则上述条件变为

新的方程组系数矩阵可写为:



下图可以看出不同的端点边界对样条曲线的影响:


1.3 算法总结


假定有n+1个数据节点

a. 计算步长 (i = 0, 1, …, n-1)

b. 将数据节点和指定的首位端点条件带入矩阵方程

c. 解矩阵方程,求得二次微分值。该矩阵为三对角矩阵,具体求法参见后一篇。

d. 计算样条曲线的系数:


其中i = 0, 1, …, n-1

e. 在每个子区间 中,创建方程

2. C语言实现

用C语言写了一个三次样条插值(自然边界)的S-Function,代码如下:

#define S_FUNCTION_NAME cubic
#define S_FUNCTION_LEVEL 2
#include "simstruc.h"
#include "malloc.h" //方便使用变量定义数组大小

static void mdlInitializeSizes(SimStruct *S)
{
/*参数只有一个,是n乘2的定点数组[xi, yi]:
* [ x1,y1;
* x2, y2;
* ..., ...;
* xn, yn;
*/

ssSetNumSFcnParams(S, 1);
if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;
ssSetNumContStates(S, 0);
ssSetNumDiscStates(S, 0);

if (!ssSetNumInputPorts(S, 1)) return; //输入是x
ssSetInputPortWidth(S, 0, 1);
ssSetInputPortRequiredContiguous(S, 0, true);
ssSetInputPortDirectFeedThrough(S, 0, 1);

if (!ssSetNumOutputPorts(S, 1)) return; //输出是S(x)
ssSetOutputPortWidth(S, 0, 1);
ssSetNumSampleTimes(S, 1);
ssSetNumRWork(S, 0);
ssSetNumIWork(S, 0);
ssSetNumPWork(S, 0);
ssSetNumModes(S, 0);
ssSetNumNonsampledZCs(S, 0);

ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE);

ssSetOptions(S, 0);
}


static void mdlInitializeSampleTimes(SimStruct *S)
{
ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
ssSetOffsetTime(S, 0, 0.0);

}

#define MDL_INITIALIZE_CONDITIONS
#if defined(MDL_INITIALIZE_CONDITIONS)
static void mdlInitializeConditions(SimStruct *S)
{
}
#endif

#define MDL_START

#if defined(MDL_START)
static void mdlStart(SimStruct *S)
{
}
#endif /* MDL_START */

static void mdlOutputs(SimStruct *S, int_T tid)
{
const real_T *map = mxGetPr(ssGetSFcnParam(S,0)); //获取定点数据
const int_T *mapSize = mxGetDimensions(ssGetSFcnParam(S,0)); //定点数组维数
const real_T *x = (const real_T*) ssGetInputPortSignal(S,0); //输入x
real_T *y = ssGetOutputPortSignal(S,0); //输出y
int_T step = 0; //输入x在定点数中的位置
int_T i;
real_T yval;

for (i = 0; i < mapSize[0]; i++)
{
if (x[0] >= map[i] && x[0] < map[i + 1])
{
step = i;
break;
}
}

cubic_getval(&yval, mapSize, map, x[0], step);
y[0] = yval;
}

//自然边界的三次样条曲线函数
void cubic_getval(real_T* y, const int_T* size, const real_T* map, const real_T x, const int_T step)
{
int_T n = size[0];

//曲线系数
real_T* ai = (real_T*)malloc(sizeof(real_T) * (n-1));
real_T* bi = (real_T*)malloc(sizeof(real_T) * (n-1));
real_T* ci = (real_T*)malloc(sizeof(real_T) * (n-1));
real_T* di = (real_T*)malloc(sizeof(real_T) * (n-1));
real_T* h = (real_T*)malloc(sizeof(real_T) * (n-1)); //x的??

/* M矩阵的系数
*[B0, C0, ...
*[A1, B1, C1, ...
*[0, A2, B2, C2, ...
*[0, ... An-1, Bn-1]
*/

real_T* A = (real_T*)malloc(sizeof(real_T) * (n-2));
real_T* B = (real_T*)malloc(sizeof(real_T) * (n-2));
real_T* C = (real_T*)malloc(sizeof(real_T) * (n-2));
real_T* D = (real_T*)malloc(sizeof(real_T) * (n-2)); //等号右边的常数矩阵
real_T* E = (real_T*)malloc(sizeof(real_T) * (n-2)); //M矩阵
real_T* M = (real_T*)malloc(sizeof(real_T) * (n)); //包含端点的M矩阵
int_T i;

//计算x的步长
for ( i = 0; i < n -1; i++)
{
h[i] = map[i + 1] - map[i];
}

//指定系数
for( i = 0; i< n - 3; i++)
{
A[i] = h[i]; //忽略A[0]
B[i] = 2 * (h[i] + h[i+1]);
C[i] = h[i+1]; //忽略C(n-1)
}

//指定常数D
for (i = 0; i3; i++)
{
D[i] = 6 * ((map[n + i + 2] - map[n + i + 1]) / h[i + 1] - (map[n + i + 1] - map[n + i]) / h[i]);
}

//求解三对角矩阵,结果赋值给E
TDMA(E, n-3, A, B, C, D);
M[0] = 0; //自然边界的首端M为0
M[n-1] = 0; //自然边界的末端M为0

for(i=1; i<>1; i++)
{
M[i] = E[i-1]; //其它的M值
}

//?算三次样条曲线的系数
for( i = 0; i < n-1; i++)
{
ai[i] = map[n + i];
bi[i] = (map[n + i + 1] - map[n + i]) / h[i] - (2 * h[i] * M[i] + h[i] * M[i + 1]) / 6;
ci[i] = M[i] / 2;
di[i] = (M[i + 1] - M[i]) / (6 * h[i]);
}

*y = ai[step] + bi[step]*(x - map[step]) + ci[step] * (x - map[step]) * (x - map[step]) + di[step] * (x - map[step]) * (x - map[step]) * (x - map[step]);

free(h);
free(A);
free(B);
free(C);
free(D);
free(E);
free(M);
free(ai);
free(bi);
free(ci);
free(di);
}

void TDMA(real_T* X, const int_T n, real_T* A, real_T* B, real_T* C, real_T* D)
{
int_T i;
real_T tmp;

//上三角矩阵
C[0] = C[0] / B[0];
D[0] = D[0] / B[0];

for(i = 1; i)
{
tmp = (B[i] - A[i] * C[i-1]);
C[i] = C[i] / tmp;
D[i] = (D[i] - A[i] * D[i-1]) / tmp;
}

//直接求出X的最后一个值
X[n-1] = D[n-1];

//逆向迭代, 求出X
for(i = n-2; i>=0; i--)
{
X[i] = D[i] - C[i] * X[i+1];
}
}

#define MDL_UPDATE
#if defined(MDL_UPDATE)
static void mdlUpdate(SimStruct *S, int_T tid)
{
}
#endif

#define MDL_DERIVATIVES
#if defined(MDL_DERIVATIVES)
static void mdlDerivatives(SimStruct *S)
{
}
#endif

static void mdlTerminate(SimStruct *S)
{
}

#ifdef MATLAB_MEX_FILE
#include "simulink.c"
#else
#include "cg_sfun.h"
#endif

3. 例子


以y=sin(x)为例, x步长为1,x取值范围是[0,10]。对它使用三次样条插值,插值前后对比如下:


三对角矩阵(Tridiagonal Matrices)的求法:Thomas Algorithm(TDMA)

做三次样条曲线时,需要解三对角矩阵(Tridiagonal Matrices)。常用解法为Thomas Algorithm,又叫The tridiagonal matrix algorithm (TDMA)。它是一种基于高斯消元法的算法, 分为两个阶段:向前消元forward elimination和回代backward substitution。本文以一个6乘6矩阵为例,介绍一下使用TDMA的求解过程。

1.范例求解



步骤1: 将矩阵变为上三角矩阵

首先要把上面公式中的系数矩阵变为一个上三角矩阵。

第一行:



将上式除以b1:



可写作:



所以矩阵方程可写为:


第二行:



将变换后的第一行乘以a2,再与第二行相减,即可消去x1,得:



所以新的矩阵方程为:


同理可推,

第三行:



第四行:



第五行:



第六行:



最后得到新的上三角矩阵公式为:


步骤2:求解

x逆序可以求出,如下:




2. 一般性公式:










注意:

使用TDMA求解,系数矩阵需时diagonally dominant, 即:



3. 实现代码(C语言)


void tdma(float x[], const size_t N, const float a[], const float b[], float c[])
{
size_t n;
c[0] = c[0] / b[0];
x[0] = x[0] / b[0];
for (n = 1; n < N; n++) {
float m = 1.0f / (b[n] - a[n] * c[n - 1]);
c[n] = c[n] * m;
x[n] = (x[n] - a[n] * x[n - 1]) * m;
}

for (n = N - 1; n-- > 0; )
x[n] = x[n] - c[n] * x[n + 1];

}


双线性插值(Bilinear Interpolation)

最近用到插值算法,使用三次样条插值时仿真速度太慢,于是采用算法简单的线性插值。本篇主要介绍一下双线性插值的实现方法。

1. 线性插值

已知坐标 (x0, y0) 与 (x1, y1),要得到 [x0, x1] 区间内某一位置 x 在直线上的值。


由于 x 值已知,所以可以从公式得到 y 的值


已知 y 求 x 的过程与以上过程相同,只是 x 与 y 要进行交换。


2. 双线性插值(Bilinear Interpolation)

在数学上,双线性插值是有两个变量的插值函数的线性插值扩展,其核心思想是在两个方向分别进行一次线性插值。



图中:红色的数据点与待插值得到的绿色点

假如我们想得到未知函数 f 在点 P = (x, y) 的值,假设我们已知函数 f 在 Q11 = (x1, y1)、Q12 = (x1, y2), Q21 = (x2, y1) 以及 Q22 = (x2, y2) 四个点的值。

首先在 x 方向进行线性插值,得到



然后在 y 方向进行线性插值,得到
                           

这样就得到所要的结果 f(x, y),



双线性插值在三维空间的延伸是三线性插值。

相关推荐
©️2020 CSDN 皮肤主题: 编程工作室 设计师:CSDN官方博客 返回首页