三次三角插值样条曲线的数控插补计算方法
三次样条插值的方法和思路
三次样条插值的方法和思路摘要：1.三次样条插值的基本概念2.三次样条插值的数学原理3.三次样条插值的实现步骤4.三次样条插值的优缺点5.三次样条插值在实际应用中的案例正文：在日常的科学研究和工程应用中，我们经常会遇到需要对一组数据进行插值的问题。
插值方法有很多，其中三次样条插值是一种常见且有效的方法。
本文将从基本概念、数学原理、实现步骤、优缺点以及实际应用案例等方面，全面介绍三次样条插值的方法和思路。
一、三次样条插值的基本概念三次样条插值（Cubic Spline Interpolation）是一种基于分段多项式的插值方法。
它通过在各个节点上构建一条三次多项式曲线，使得这条曲线在节点之间满足插值条件，从而达到拟合数据的目的。
二、三次样条插值的数学原理三次样条插值的数学原理可以分为两个部分：一是分段三次多项式的构建，二是插值条件的满足。
1.分段三次多项式的构建假设有一组数据点序列为（x0，y0），（x1，y1），（x2，y2），（x3，y3），我们可以将这些数据点连接起来，构建一条分段三次多项式曲线。
分段三次多项式在每个子区间上都是一个三次多项式，它们之间通过节点值进行连接。
2.插值条件的满足为了使分段三次多项式在节点之间满足插值条件，我们需要在每个子区间上满足以下四个条件：（1）端点条件：三次多项式在区间的端点上分别等于节点值；（2）二阶导数条件：三次多项式在区间内的二阶导数等于节点间的斜率；（3）三阶导数条件：三次多项式在区间内的三阶导数等于节点间的曲率；（4）内部点条件：三次多项式在区间内部满足插值函数的连续性。
通过求解这四个条件，我们可以得到分段三次多项式的系数，从而实现插值。
三、三次样条插值的实现步骤1.确定插值节点：根据数据点的位置，选取合适的节点；2.构建分段三次多项式：根据节点值和插值条件，求解分段三次多项式的系数；3.计算插值结果：将待插值点的横坐标代入分段三次多项式，得到插值结果。
三次样条插值函数求解例题
三次样条插值函数求解例题三次样条插值函数是一种常用的插值方法，用于在给定的一组数据点上构建一个连续的曲线。
下面我将通过一个例题来解释三次样条插值函数的求解过程。
假设我们有一组数据点{(x0, y0), (x1, y1), ..., (xn, yn)}，其中x0 < x1 < ... < xn。
我们的目标是构建一个连续的曲线，使得曲线经过这些数据点。
首先，我们需要确定每个数据点之间的插值多项式。
在三次样条插值中，每个插值多项式的形式为：Si(x) = ai + bi(x xi) + ci(x xi)^2 + di(x xi)^3。
其中，ai、bi、ci、di是待求的系数，Si(x)是第i段插值多项式。
接下来，我们需要确定每个插值多项式的系数。
为了满足插值条件，我们需要确定每个数据点处的函数值和导数值。
具体而言，我们需要满足以下条件：1. 函数值条件，Si(xi) = yi，即插值多项式通过每个数据点。
2. 导数值条件，Si'(xi) = Si-1'(xi)，即相邻插值多项式在数据点处的导数值相等。
通过这些条件，我们可以得到一系列的线性方程组，其中未知数为插值多项式的系数。
解这个线性方程组即可得到每个插值多项式的系数。
最后，我们可以将每个插值多项式的系数代入到对应的插值多项式中，得到最终的三次样条插值函数。
需要注意的是，在边界处，我们需要额外的条件来确定插值多项式的系数。
常见的边界条件有自然边界条件和固定边界条件。
自然边界条件要求插值函数的二阶导数在边界处为零，而固定边界条件要求插值函数在边界处通过给定的导数值。
综上所述，三次样条插值函数的求解过程包括确定插值多项式的系数和边界条件的确定。
通过解线性方程组，我们可以得到每个插值多项式的系数，从而构建出连续的三次样条插值函数。
希望以上回答能够满足你的要求。
如果你有任何其他问题，请随时提出。
详细讲解三次样条插值法及其实现方法
样条函数的定义 定义4.1 设区间[a,b]上给定一个节点划分
a=x0<x1<……<xn-1<xn=b 如果存在正整数k使得[a,b]上的分段函数s(x)满足 如下两条： (1)在[a,b]上有直到k-1阶连续导数。 (2)在每个小区间[xi,xi+1]上是次数不大于k的多项式。 则称分段函数s(x)是以(2.6)为节点集的k次样条函数。
x xi i i 1 hi
) mi1hi
( ) xi1x
1 hi
x [xi , xi1], hi xi1 xi , i 0,1,, n 1
(x) (2x 1)( x 1)2,1(x) x(x 1)2 13
对Si (x)求二阶导数 ,并整理后得
Si( x)
6( xi
xi 1 hi3
2x)
因为分段三次Hermite插值多项式已经至少是一阶连续 可导了，为了让它成为三次样条函数只需确定节点处 的一阶导数使这些节点处的二阶导数连续即可！
S(xi 0) S(xi 0), i 1,, n 1
S(x)
y ( xxi i 0 hi
)
y ( ) m h ( xi1x i1 0 hi
( yn
yn 1 )
2 hn1
(mn1
2mn )
立即可得下式：
21
其中：
nm1 nmn1 2mn gn
n
h0
h0 hn1
, n
hn1 h0 hn1
1 n
gn
3 n
y1 y0 h0
n
yn
yn1 hn1
联合基本方程得一个广义三对角或周期三对角方程组：
2 1
1
1
2
样条函数及三次样条插值PPT课件
(x)
lim
x xk
Sk 1( x)
lim
x
x
k
Sk (x)
lim
x
x
k
Sk1( x)
k 1,2,,n 1
------(4)
lim
x
x
k
Sk( x)
lim
x
x
k
Sk1( x)
共4n 2个条件
5
Sk (x)是[xk , xk 1 ]上的三次样条插值多项式,应有4个待定的系数 即要确定S(x)必须确定4n个待定的系数 少两个条件 并且我们不能只对插值函数在中间节点的状态进行限制 也要对插值多项式在两端点的状态加以要求 也就是所谓的边界条件:
例. 使用不同的插值方法于函数
y
1
1 x2
x [5,5]
最后，介绍一个有用的结论
定理 . 设f (x) C 2[a,b], S(x)是以xk (k 0,1,, n)
为节点, 满足任意边界条件的三次样条插值函数，
设hi
xi 1
xi
,
h
max
0in1
hi
,
min
0in1
hi
,
则当 h
c 时
S(x)和S(x)在[a,b]上一致收敛到f (x)和f (x)
------(6)
13
由(11)式,可知
S0( x0
)
6( x0
x1 h03
2 x0
) ( y1
y0 )
6 x0
2 x0 h02
4 x1
m0
6 x0
4 x0 h02
2 x1
m1
6 h02
(
数值分析——样条函数及三次样条插值
S k ( x )是[ xk , xk + 1 ]上的(两点)三次样条插值多项式, 满足
Sk ( x j ) = y j
x → xk
k = 0,1,2, ⋯ , n − 1; j = k , k + 1
lim S k ( x ) = lim S k − 1 ( x ) + −
一、三次样条插值函数
定义1. 定义
a ≤ x0 , x1 ,⋯ , xn ≤ b为区间[ a , b ]的一个分割
如果函数S ( x )在区间[ a , b ]上满足条件 :
( 1) S ( x ), S ′( x ), S ′′( x )都在区间[ a , b ]上连续 ,即 S ( x ) ∈ C 2 [ a , b ]
f ( x j ) = y j , j = 0 ,1,⋯ , n 如果S ( x )是f ( x )的三次样条插值函数, 则其必满足
S ( x j ) = y j , j = 0 ,1,⋯ , n lim S ( x ) = S ( x j ) = y j , j = 1,⋯ , n − 1 x→ x xlim S ′( x ) = S ′( x j ) = m j , j = 1,⋯ , n − 1 →x lim S ′′( x ) = S ′′( x j ), j = 1 ,⋯ , n − 1
+ x → xk
Sk ( x j ) = y j
k
k = 0,1, ⋯ , n − 1; j = k , k + 1
− x → xk
k −1
k = 1 , 2 ,⋯ , n − 1 k = 1,2 ,⋯ , n − 1 ------(8) k = 1, 2 ,⋯ , n − 1
三次样条插值算法详解
三次样条插值算法要求数据点数量较多，且在某些情况下可能存在数值不稳定性，如数据 点过多或数据点分布不均等情况。此外，该算法对于离散数据点的拟合效果可能不如其他 插值方法。
对未来研究的展望
01
02
03
改进算法稳定性
针对数值不稳定性问题， 未来研究可以探索改进算 法的数值稳定性，提高算 法的鲁棒性。
3
数据转换
对数据进行必要的转换，如标准化、归一化等， 以适应算法需求。
构建插值函数
确定插值节点
根据数据点确定插值节点，确保插值函数在节点处连续且光滑。
构造插值多项式
根据节点和数据点，构造三次多项式作为插值函数。
确定边界条件
根据实际情况确定插值函数的边界条件，如周期性、对称性等。
求解插值函数
求解线性方程组
06
结论
三次样条插值算法总结
适用性
三次样条插值算法适用于各种连续、光滑、可微的分段函数插值问题，尤其在处理具有复 杂变化趋势的数据时表现出色。
优点
该算法能够保证插值函数在分段连接处连续且具有二阶导数，从而在插值过程中保持数据 的平滑性和连续性。此外，三次样条插值算法具有简单、易实现的特点，且计算效率较高 。
根据数据点的数量和分布，合理分段，确保 拟合的精度和连续性。
求解线性方程组
使用高效的方法求解线性方程组，如高斯消 元法或迭代法。
结果输出
输出拟合得到的插值函数，以及相关的误差 分析和图表。
03
三次样条插值算法步骤
数据准备
1 2
数据收集
收集需要插值的原始数据点，确保数据准确可靠。
数据清洗
对数据进行预处理，如去除异常值、缺失值处理 等。
常见插值算法--拉格朗日插值、三次卷积插值、三次样条插值、兰克索斯插值
常见插值算法--拉格朗⽇插值、三次卷积插值、三次样条插值、兰克索斯插值写在前⾯本⽂简单介绍了⼏种常见的插值算法并附带了相应的python代码，本⽂公式使⽤latex编写，如有错误欢迎评论指出，如果谁知道如何修改latex字号也欢迎留⾔关于⼀维、⼆维和多维插值三次卷积插值、拉格朗⽇两点插值（线性插值）、兰克索斯插值在⼆维插值时改变x和y⽅向的计算顺序不影响最终结果，这三个也是图像缩放插值时常⽤的插值算法，⽽其他插值在改变计算顺序时会产⽣明显差异，多维的情况笔者没有尝试，读者可以⾃⾏尝试或推导最近邻插值法（Nearest Neighbour Interpolation）在待求像素的四邻像素中，将距离待求像素最近的像素值赋给待求像素p_{11}p_{12}pp_{21}p_{22}python代码1def NN_interpolation(srcImg, dstH, dstW):2 scrH, scrW, _ = srcImg.shape3 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)4for i in range(dstH - 1):5for j in range(dstW - 1):6 scrX = round(i * (scrH / dstH))7 scrY = round(j * (scrW / dstW))8 dstImg[i, j] = srcImg[scrX, scrY]9return dstImg拉格朗⽇插值（Lagrange Interpolation）拉格朗⽇插值法需要找到k个p_i(x)函数，使得每个函数分别在在x_i处取值为1，其余点取值为0，则y_ip_i(x)可以保证在x_i处取值为y_i，在其余点取值为0，因此L_k(x)能恰好经过所有点，这样的多项式被称为拉格朗⽇插值多项式，记为L_k(x)=\sum_{i=1}^ky_ip_i(x)p_i(x)=\prod_{j \neq i}^{1 \leq j \leq k}\frac{x-x_j}{x_i-x_j}以四点即三次图像插值为例，因为横坐标间隔为1，则设四个点横坐标为-1、0、1和2，可得p_1(x)、p_2(x)、p_3(x)和p_4(x)假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得拉格朗⽇函数如下图所⽰，待插值点横坐标范围为[0,1]在K=2时在k=2时，也被称为线性插值通⽤公式p_1=\frac{x-x_2}{x_1-x_2}p_2=\frac{x-x_1}{x_2-x_1}\begin{align} L_2x &= p_1y_1+p_2y_2 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}pp_{21}p_{22}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_2x &= \frac{x-x_2}{x_1-x_2}y_1 + \frac{x-x_1}{x_2-x_1}y_2 \nonumber \\ &= (x_2-x)y_1+(x-x_1)y_2 \nonumber \\ &= (1-dx)y_1+dxy_2 \nonumber \\ &= (y_2-y_1)dx+y_1 \nonumber \end{align}L_2'x=y_2-y_1在K=3时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\begin{align} L_3x &= p_1y_1+p_2y_2+p_3y_3 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1+\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2+\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \end{align}图像插值像素分布如图所⽰p_{11}p_{12}p_{13}p_{21}p_{22}p_{23}pp_{31}p_{32}p_{33}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_3x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}y_3 \nonumber \\ &= \frac{-dx(1-dx)}{(-1)\cdot(-2)}y_1 + \frac{-(1+dx)(1-dx)}{1\cdot(-1)}y_2 + \frac{(1+dx)dx}{2\cdot 1}y_3 \nonumber \\ &= (\frac{1}{2}d^2x-\frac{1}{2}dx)y_1 - (d^2x-1)y_2 + (\frac{1}{2}d^2x+\frac{1}{2}dx)y_3 \nonumber \\ &= d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{2}y_1+\frac{1}{2}y_3)+y_2 \nonumber \end{align}L_3'x=dx(y_1-2y_2+y_3)+(\frac{1}{2}y_3-\frac{1}{2}y_1)在K=4时通⽤公式p_1=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}p_2=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}p_3=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}p_4=\frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}\begin{align} L_4x &= p_1y_1+p_2y_2+p_3y_3+p_4y_4 \nonumber \\ &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3} {x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4\nonumber \end{align}图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}即当x_{i+1}=x_i+1时，设p与p_{11}的横纵坐标差分别为dx和dy\begin{align} L_4x &= \frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}\frac{x-x_4}{x_1-x_4}y_1 + \frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}\frac{x-x_4}{x_2-x_4}y_2 + \frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}\frac{x-x_4}{x_3-x_4}y_3 + \frac{x-x_1}{x_4-x_1}\frac{x-x_2}{x_4-x_2}\frac{x-x_3}{x_4-x_3}y_4 \nonumber \\ &= \frac{dx[-(1-dx)][-(2-dx)]}{(-1)\cdot(-2)\cdot(-3)}y_1 + \frac{(1+dx)[-(1-dx)][-(2-dx)]}{1\cdot(-1)\cdot(-2)}y_2 + \frac{(1+dx)dx[-(2-dx)]}{2\cdot 1\cdot(-1)}y_3 + \frac{(1+dx)dx[-(1-dx)]}{3\cdot 2\cdot 1}y_4 \nonumber \\ &= \frac{d^3x-3d^2x+2dx}{-6}y1 + \frac{d^3x-2d^2x-dx+2}{2}y_2 + \frac{d^3x-d^2x-2dx}{-2}y_3 + \frac{d^3x-dx}{6}y_4 \nonumber \\ &= d^3x(-\frac{1}{6}y_1+\frac{1}{2}y_2-\frac{1} {2}y_3+\frac{1}{6}y_4)+d^2x(\frac{1}{2}y_1-y_2+\frac{1}{2}y_3)+dx(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4)+y_2 \nonumber \end{align}\begin{align} L_4'x &= d^2x(-\frac{1}{2}y_1+\frac{3}{2}y_2-\frac{3}{2}y_3+\frac{1}{2}y_4)+dx(y_1-2y_2+y_3)+(-\frac{1}{3}y_1-\frac{1}{2}y_2+y_3-\frac{1}{6}y_4) \nonumber \\ &= -[\frac{1}{2}d^2x(y_1-3y_2+3y_3-y_4)-dx(y_1-2y_2+y_3)+\frac{1}{6}(2y_1+3y_2-6y_3+y_4)] \nonumber \end{align}python代码插值核计算的时候乘法和加减法计算的顺序不同可能会导致结果存在细微的差异，读者可以⾃⾏研究⼀下1class BiLagrangeInterpolation:2 @staticmethod3def LagrangeInterpolation2(x, y1, y2):4 f1 = 1 - x5 f2 = x6 result = y1 * f1 + y2 * f27return result89 @staticmethod10def LagrangeInterpolation3(x, y1, y2, y3):11 f1 = (x ** 2 - x) / 2.012 f2 = 1 - x ** 213 f3 = (x ** 2 + x) / 2.014 result = y1 * f1 + y2 * f2 + y3 * f315return result1617 @staticmethod18def LagrangeInterpolation4(x, y1, y2, y3, y4):19 f1 = - (x ** 3 - 3 * x ** 2 + 2 * x) / 6.020 f2 = (x ** 3 - 2 * x ** 2 - x + 2) / 2.021 f3 = - (x ** 3 - x ** 2 - 2 * x) / 2.022 f4 = (x ** 3 - x) / 6.023 result = y1 * f1 + y2 * f2 + y3 * f3 + y4 * f424return result2526def biLag2_2(self, srcImg, dstH, dstW):27 dstH, dstW = int(dstH), int(dstW)28 srcH, srcW, _ = srcImg.shape29 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')30 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)31for dstY in range(dstH):32for dstX in range(dstW):33for channel in [0, 1, 2]:34# p11 p1235# p36# p21 p2237# 储存为 p(y, x)38 p = [dstY * srcH / dstH, dstX * srcW / dstW]39 p11 = [math.floor(p[0]), math.floor(p[1])]40 p12 = [p11[0], p11[1] + 1]4142 p21 = [p11[0] + 1, p11[1]]43 p22 = [p21[0], p12[1]]4445 diff_y, diff_x = p[0] - p11[0], p[1] - p11[1]46 r1 = grangeInterpolation2(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel])47 r2 = grangeInterpolation2(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel])4849 c = grangeInterpolation2(diff_y, r1, r2)5051 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)52return dstImg5354def biLag3_3(self, srcImg, dstH, dstW):55 dstH, dstW = int(dstH), int(dstW)56 srcH, srcW, _ = srcImg.shape57 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')58 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)59for dstY in range(dstH):60for dstX in range(dstW):61for channel in [0, 1, 2]:62# p11 p12 p1363#64# p21 p22 p2365# p66# p31 p32 p3367# 储存为 p(y, x)68 p = [dstY * srcH / dstH, dstX * srcW / dstW]69 p22 = [math.floor(p[0]), math.floor(p[1])]70 p21 = [p22[0], p22[1] - 1]71 p23 = [p22[0], p22[1] + 1]7273 p11 = [p21[0] - 1, p21[1]]74 p12 = [p11[0], p22[1]]75 p13 = [p11[0], p23[1]]7677 p31 = [p21[0] + 1, p21[1]]78 p32 = [p31[0], p22[1]]79 p33 = [p31[0], p23[1]]8081 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]82 r1 = grangeInterpolation3(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel])83 r2 = grangeInterpolation3(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel])84 r3 = grangeInterpolation3(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel]) 8586 c = grangeInterpolation3(diff_y, r1, r2, r3)8788 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)89return dstImg9091def biLag4_4(self, srcImg, dstH, dstW):92 dstH, dstW = int(dstH), int(dstW)93 srcH, srcW, _ = srcImg.shape94 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')95 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)96for dstY in range(dstH):97for dstX in range(dstW):98for channel in [0, 1, 2]:99# p11 p12 p13 p14100#101# p21 p22 p23 p24102# p103# p31 p32 p33 p34104#105# p41 p42 p43 p44106# 储存为 p(y, x)107 p = [dstY * srcH / dstH, dstX * srcW / dstW]108 p22 = [math.floor(p[0]), math.floor(p[1])]109 p21 = [p22[0], p22[1] - 1]110 p23 = [p22[0], p22[1] + 1]111 p24 = [p22[0], p22[1] + 2]112113 p11 = [p21[0] - 1, p21[1]]114 p12 = [p11[0], p22[1]]115 p13 = [p11[0], p23[1]]116 p14 = [p11[0], p24[1]]117118 p31 = [p21[0] + 1, p21[1]]119 p32 = [p31[0], p22[1]]120 p33 = [p31[0], p23[1]]121 p34 = [p31[0], p24[1]]122123 p41 = [p21[0] + 2, p21[1]]124 p42 = [p41[0], p22[1]]125 p43 = [p41[0], p23[1]]126 p44 = [p41[0], p24[1]]127128 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]129 r1 = grangeInterpolation4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 130 r2 = grangeInterpolation4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 131 r3 = grangeInterpolation4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 132 r4 = grangeInterpolation4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 133134 c = grangeInterpolation4(diff_y, r1, r2, r3, r4)135136 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)137return dstImg三次卷积插值法（Cubic Convolution Interpolation）使⽤上图中的卷积核进⾏加权平均计算，卷积核为u(s)，四个等距（距离为1）的采样点记为x_0、x_1、x_2和x_3，采样数值记为y_0、y_1、y_2和y_3，且保证四个点均在[-2,2]区间上，计算得到g(x)，假设y_1、y_2、y_3和y_4分别为1、2、-1、4，则可得三次卷积插值函数如下图所⽰，待插值点横坐标范围为[0,1]公式推导设u(s)=\begin{cases} A_1|s|^3+B_1|s|^2+C_1|s|+D_1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\because函数在s=0,1,2处连续\therefore\begin{cases} 1=u(0^+)=D_1 \\ 0=u(1^-)=A_1+B_1+C_1+D_1 \\ 0=u(1^+)=A_2+B_2+C_2+D_2 \\ 0=u(2^-)=8A_2+4B_2+2C_2+D_2 \end{cases} (1)\because函数在s=0,1,2处导函数连续\therefore\begin{cases} u'(0^-)=u'(0+) \\ u'(1^-)=u'(1+) \\ u'(2^-)=u'(2+)\end{cases} \Rightarrow \begin{cases} -C_1=C_1 \\ 3A_1+2B_1+C_1=3A_2+2B_2+C_2\\ 12A_2+4B_2+C+2=0 \end{cases} ~~~~ (2)联⽴⽅程组(1)(2)，设A_2=a，解得\begin{cases} A_1=a+2 \\ B_1=-(a+3) \\ C_1=0 \\ D_1=1 \\ A_2=a \\ B_2=-5a \\ C_2=8a \\ D_2=-4a \end{cases}\Rightarrow u(s)=\begin{cases} (a+2)|s|^3-(a+3)|s|^2+1, &0<|s|<1 \\ A_2|s|^3+B_2|s|^2+C_2|s|+D_2, &1<|s|<2\\ 1, &s=0 \\ 0, &otherwise \end{cases}\because g(x)=\sum_kC_ku(s+j-k), ~~~~k=j-1,j, j+1,j+2且0<s<1⼜\because \begin{cases}\begin{align} u(s+1)&=as^3-2as^2+as \nonumber \\ u(s)&=(a+2)s^3-(a+3)s^2+1 \nonumber \\ u(s-1)&=-(a+2)s^3+(2a+3)s^2-as \nonumber \\ u(s-2)&=-as^3+as^2 \nonumber \end{align}\end{cases}\begin{align} \therefore g(x) &= C_{j-1}u(s+1)+C_{j}u(s)+C_{j+1}u(s-1)+C_{j+2}u(s-2) \nonumber \\ &= C_{j-1}(as^3-2as^2+as)+C_j[(a+2)s^3-(a+3)s^2+1]+C_{j+1}[-(a+2)s^3+ (2a+3)s^2-as]+C_{j+2}[-a^3+as^2] \nonumber \\ &= s^3[aC_{j-1}+(a+2)C_j-(a+2)C_{j+1}-aC_{j+2}]+s^2[-2aC_{j-1}-(a+3)C_j+(2a+3)C_{j+1}+aC_{j+2}]+s[aC_{j-1}-aC_{j+1}]+C_j \nonumber \end{align} ~~(3)f在x_j处泰勒展开得到f(x)=f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+\cdots\therefore \begin{cases} f(x_{j+1})=f(x_j)+f'(x_j)(x_{j+1}-x_j)+\frac{1}{2}f''(x_j)(x_{j+1}-x_j)^2+\cdots \\ f(x_{j+2})=f(x_j)+f'(x_j)(x_{j+2}-x_j)+\frac{1}{2}f''(x_j)(x_{j+2}-x_j)^2+\cdots \\ f(x_{j-1})=f(x_j)+f'(x_j)(x_{j-1}-x_j)+\frac{1}{2}f''(x_j)(x_{j-1}-x_j)^2+\cdots \end{cases}令x_{j+1}-x_j=h\because x_{i+1}=x_i+1\therefore x_{j+2}-x_j=2h，x_{j-1}-x_j=-h\therefore \begin{cases} f(x_{j+2})=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+\cdots \\ f(x_{j+1})=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \\ f(x_{j-1})=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+\cdots \end{cases}\therefore \begin{cases} c_{j-1}=f(x_j)-f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3) \\ c_j=f(x_j) \\ c_{j+1}=f(x_j)+f'(x_j)h+\frac{1}{2}f''(x_j)h^2+o(h^3)\\ c_{j+2}=f(x_j)+2f'(x_j)h+2f''(x_j)h^2+o(h^3) \end{cases} ~~ (4)将(4)代⼊(3)，得g(x)=-(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3+[(6a+3)hf'(x_j)+\frac{4a+3}{2}h^2f''(x_j)]s^2-2ahf'(x_j)s+f(x_j)+o(h^3)\because h=1,s=x-x_J\therefore sh=x-x_j\begin{align}\therefore f(x)&= f(x_j)+f'(x_j)(x-x_j)+\frac{1}{2}f''(x_j)(x-x_j)^2+o(h^3) \nonumber \\ &= f(x_j)+f'(x_j)sh+\frac{1}{2}f''(x_j)s^2h^2+o(h^3) \nonumber \end{align}\therefore f(x)-g(x)=(2a+1)[2hf'(x_j)+h^2f''(x_j)]s^3-(2a+1)[3hf'(x_j)+h^2f''(x_j)]s^2+[(2a+1)hf'(x_j)]s+o(h^3)\because 期望f(x)-g(x)趋于0\therefore 2a+1=0 \Rightarrow a=-\frac{1}{2}\therefore u(s)=\begin{cases} \frac{3}{2}|s|^3-\frac{5}{2}|s|^2+1, &0<|s|<1 \\ -\frac{1}{2}|s|^3+\frac{5}{2}|s|^2-4|s|+2, &1<|s|<2 \\ 1, &s=0 \\ 0, &otherwise \end{cases}\therefore g(s)=s^3[-\frac{1}{2}c_{j-1}+\frac{3}{2}c_j-\frac{3}{2}c_{j+1}+\frac{1}{2}c_{j+2}]+s^2[c_{j-1}-\frac{5}{2}c_j+2c_{j+1}-\frac{1}{2}c_{j+2}]+s[-\frac{1}{2}c_{j-1}+\frac{1} {2}c_{j+1}]+c_j图像插值p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}python代码1class BiCubicConvInterpolation:2 @staticmethod3def CubicConvInterpolation1(p0, p1, p2, p3, s):4# ⽤g（s）公式计算，已经将四个u（s）计算完毕并整理5# as^3 + bs^2 + cs + d6 a = 0.5 * (-p0 + 3.0 * p1 - 3.0 * p2 + p3)7 b = 0.5 * (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3)8 c = 0.5 * (-p0 + p2)9 d = p110return d + s * (c + s * (b + s * a))1112 @staticmethod13def CubicConvInterpolation2(s):14# ⽤u（s）公式计算15 s = abs(s)16if s <= 1:17return 1.5 * s ** 3 - 2.5 * s ** 2 + 118elif s <= 2:19return -0.5 * s ** 3 + 2.5 * s ** 2 - 4 * s + 220else:21return 02223def biCubic1(self, srcImg, dstH, dstW):24# p11 p12 p13 p1425#26# p21 p22 p23 p2427# p28# p31 p32 p33 p3429#30# p41 p42 p43 p4431 dstH, dstW = int(dstH), int(dstW)32 scrH, scrW, _ = srcImg.shape33 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')34 dstImg = np.zeros((dstH, dstW, 1), dtype=np.uint8)35for dstY in range(dstH):36for dstX in range(dstW):37for channel in [0]:38 y = dstY * scrH / dstH39 x = dstX * scrW / dstW40 y1 = math.floor(y)41 x1 = math.floor(x)4243 array = []44for i in [-1, 0, 1, 2]:45 temp = self.CubicConvInterpolation1(srcImg[y1 + i, x1 - 1, channel],46 srcImg[y1 + i, x1, channel],47 srcImg[y1 + i, x1 + 1, channel],48 srcImg[y1 + i, x1 + 2, channel],49 x - x1)50 array.append(temp)5152 temp = self.CubicConvInterpolation1(array[0], array[1], array[2], array[3], y - y1)53 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)5455return dstImg5657def biCubic2(self, srcImg, dstH, dstW):58# p11 p12 p13 p1459#60# p21 p22 p23 p2461# p62# p31 p32 p33 p3463#64# p41 p42 p43 p4465 dstH, dstW = int(dstH), int(dstW)66 scrH, scrW, _ = srcImg.shape67 srcImg = np.pad(srcImg, ((1, 1), (1, 1), (0, 0)), 'edge')68 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)69for dstY in range(dstH):70for dstX in range(dstW):71for channel in [0, 1, 2]:72 y = dstY * scrH / dstH73 x = dstX * scrW / dstW74 y1 = math.floor(y)75 x1 = math.floor(x)7677 array = []78for i in [-1, 0, 1, 2]:79 temp = 080for j in [-1, 0, 1, 2]:81 temp += srcImg[y1 + i, x1 + j, channel] * self.CubicConvInterpolation2(x - (x1 + j))82 array.append(temp)8384 temp = 085for i in [-1, 0, 1, 2]:86 temp += array[i + 1] * self.CubicConvInterpolation2(y - (y1 + i))87 dstImg[dstY, dstX, channel] = np.clip(temp, 0, 255)8889return dstImg三次样条插值在n-1个区间上寻找n-1个三次曲线，使其满⾜相邻曲线在端点处值相等、⼀阶导数相等，⼆阶导数相等，在加以边界条件后可得每个曲线的⽅程，然后沿x轴依次偏移对应的距离即可得到插值结果，如仅需要特定范围内的结果，则可以⼤幅减少计算量公式推导设S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3, ~~~~i=0,1,...,n-1则 \begin{cases} S_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2\\ S_i''(x)=2c_i+6d_i(x-x_i)\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1设h_i(x)=x_{i+1}-x_i，可得\begin{cases} S_i(x)=a_i+b_ih_i+c_ih_i^2+d_ih_i^3\\ S_i'(x)=b_i+2c_ih_i+3d_ih_i^2\\ S_i''(x)=2c_i+6d_ih_i\\ S_i'''(x)=6d_i\\ \end{cases} ~~~~i=0,1,...,n-1\because S_i(x)过点(x_i,y_i)\therefore S_i(x)=a_i=y+i, ~~~~i=0,1,...,n-1 ~~~~~~(1)\because S_i(x)与S_{i+1}(x)在X_{i+1}处相等\therefore S_i(x_{i+1})=S_{i+1}(x_{i+1})\Rightarrow a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1}, ~~~~i=0,1,...,n-2~~~~~~(2)\because S_i'(x)与S_{i+1}'(x)在X_{i+1}处相等\therefore S_i'(x)-S_{i+1}'(x)=0\Rightarrow b_i+2c_ih_i+3d_ih_i^2-b_{i+1}=0~~~~~~(3)\because S_i''(x)与S_{i+1}''(x)在X_{i+1}处相等\therefore S_i''(x)-S_{i+1}''(x)=0\Rightarrow 2c_i+6d_ih_i-2c_{i+1}=0, ~~~~i=0,1,...,n-2~~~~~~(4)设m_i=S_i(x_i)=2c_i，即c_i=\frac{1}{2}m_i, ~~~~i=0,1,...,n-1~~~~~~(5)将(5)代⼊(4)，得2c_i+6d_ih_i-2c_{i+1}=0\Rightarrow m_i+6h_id_i-m_{i+1}=0\Rightarrow d_i=\frac{m_{i+1}-m_i}{6h_i}, ~~~~i=0,1,...,n-2~~~~~~(6)将(1)(5)(6)代⼊(2)，得\begin{align} &a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1} \nonumber \\ \Rightarrow&y_i+b_ih_i+\frac{1}{2}m_ih_i^2+\frac{m_{i+1}-m_i}{6h_i}h_i^3=y_{i+1} \nonumber \\\Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{2}m_ih_i-\frac{1}{6}(m_{i+1}-m_i)h_i \nonumber \\ \Rightarrow&b_i=\frac{y_{i+1}-y_i}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i, ~~~~i=0,1,...,n-2~~~~~~(7) \nonumber \end{align}将(5)(6)(7)代⼊(3)，得\begin{align} &\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+2\cdot\frac{1}{2}m_ih_i+3\frac{m_{i+1}-m_i}{6h_i}h_i^2-(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{1}{3}m_{i+1}h_{i+1}-\frac{1}{6}m_{i+2}h_{i+1})=0 \nonumber \\ \Rightarrow&\frac{y_{i+1}-y{i}}{h_i}-\frac{1}{3}m_ih_i-\frac{1}{6}m_{i+1}h_i+m_ih_i+\frac{1}{2}(m_{i+1}-m_i)h_i-\frac{y_{i+2}-y_{i+1}}{h_{i+1}}+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=0 \nonumber \\ \Rightarrow&m_ih_i(-\frac{1}{3}+1-\frac{1}{2})+m_{i+1}h_i(-\frac{1}{6}+\frac{1} {2})+\frac{1}{3}m_{i+1}h_{i+1}+\frac{1}{6}m_{i+2}h_{i+1}=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&\frac{1}{6}(m_ih_i+2m_{i+1}h_i+2m_{i+1}h_{i+1}+m_{i+2}h_{i+1})=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}} \nonumber \\ \Rightarrow&m_ih_i+2m_{i+1}(h_i+h_{i+1})+m_{i+2}h_{i+1}=6(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}}), ~~~~i=0,1,...,n-2~~~~~~(8) \nonumber \end{align}由(8)可知i=0,1,...,n-2，则有m_0,m_1,...,m_n，需要两个额外条件⽅程组才有解⾃然边界(Natural)m_0=0，m_n=0\begin{bmatrix} \tiny 1 & 0 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}固定边界(Clamped)\begin{align} &\begin{cases} S_0'(x_0)=A\\ S_{n-1}'(x_n)=B \end{cases} \nonumber \\ \Rightarrow&\begin{cases} b_0=A\\ b_{n-1}+2c_{n-1}h_{n-1}+3d_{n-1}h_{n-1}^2=B\end{cases} \nonumber \\ \Rightarrow&\begin{cases} A=\frac{y_1-y_0}{h_0}-\frac{h_0}{2}m_0-\frac{h_0}{6}(m_1-m_0)\\ B=\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{1}{3}m_{n-1}h_{n-1}+m_{n-1}h_{n-1}+\frac{1}{2}m_nh_{n-1}-\frac{1}{2}m_{n-1}h_{n-1} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 2h_0m_0+h_0m_1=6(\frac{y_1-y_0}{h_0}-A)\\ h_{n-1}m_{n-1}+2h_{n-1}m_{n}=6(B-\frac{y_n-y_{n-1}}{h_{n-1}}) \end{cases} \nonumber \\ \end{align}\begin{bmatrix} 2 & 1 & 0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 & 2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & 0 & 1 & 2 \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\\frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ B-\frac{y_n-y_{n-1}}{h_{n-1}} \end{bmatrix}⾮节点边界(Not-A-Knot)\begin{align} &\begin{cases} S_0'''(x_1)=S_1'''(x_1)\\ S_{n-2}'''(x_{n-1})=S_{n-1}'''(x_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} 6\cdot\frac{m_1-m_0}{6h_0}=6\cdot\frac{m_2-m_1}{6h_1}\\ 6\cdot\frac{m_{n-1}-m_{n-2}}{6h_{n-2}}=6\cdot\frac{m_n-m_{n-1}}{6h_{n-1}} \end{cases} \nonumber \\ \Rightarrow&\begin{cases} h_1(m_1-m_0)=h_0(m_2-m_1)\\ h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_n-m_{n-1}) \end{cases} \nonumber \\ \Rightarrow&\begin{cases} -h_1m_0+(h_1+h_0)m_1-h_0m_2=0\\ -h_{n-1}m_{n-2}+(h_{n-1}+h_{n-2})m_{n-1}-h_{n-2}m_n=0 \end{cases} \nonumber \\ \end{align}\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 & 0 & \cdots & 0\\ h_0 & 2(h_0+h_1) & h_1 & 0 & 0 & \cdots & 0\\ 0 & h_1 & 2(h_1+h_2) & h_2 & 0 & \cdots & 0\\ 0 & 0 & h_2 &2(h_2+h_3) & h_3 & \cdots & 0\\ \vdots& & & \ddots & \ddots & \ddots & \vdots\\ 0 & \cdots & & & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1}\\ 0 & \cdots & & & -h_{n-1} & h_{n-1}+h_{n-2} & -h_{n-2} \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3\\\vdots\\m_{n-1}\\m_n \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ \frac{y_4-y_3}{h_3}-\frac{y_3-y_2}{h_2}\\ \vdots\\ \frac{y_n-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\\ 0 \end{bmatrix}在n=4时通⽤公式⾃然边界\begin{bmatrix} 1 & 0 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}固定边界\begin{bmatrix} 2 & 1 & 0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} \frac{y_1-y_0}{h_0}-A\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ B-\frac{y_3-y_2}{h_2} \end{bmatrix}⾮节点边界\begin{bmatrix} -h_1 & h_0+h_1 & -h_0 & 0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 \\ 0 & -h_2 & h_1+h_2 & -h_1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0}\\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1}\\ 0 \end{bmatrix}图像插值x_{i+1}-x_i=1 \Rightarrow h_i(x)=1在n=4时，即四个点时如下所⽰p_{11}p_{12}p_{13}p_{14}p_{21}p_{22}p_{23}p_{24}pp_{31}p_{32}p_{33}p_{34}p_{41}p_{42}p_{43}p_{44}⾃然边界（可⽤TDMA或化简计算）\begin{bmatrix} 1 & 0 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}固定边界（只能⽤TDMA计算）\begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} y_1-y_0-A\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ y_2-y_3+B \end{bmatrix}⾮节点边界（只能化简计算）\begin{bmatrix} -1 & 2 & -1 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & -1 & 2 & -1 \\ \end{bmatrix}\begin{bmatrix} m_0\\m_1\\m_2\\m_3 \end{bmatrix}=6\begin{bmatrix} 0\\ y_0+y_2-2y_1\\ y_1+y_3-2y_2\\ 0 \end{bmatrix}python代码1class BiSplineInterpolation:2 @staticmethod3 def TDMA(a, b, c, d):4 n = len(d)56 c[0] = c[0] / b[0]7 d[0] = d[0] / b[0]89for i in range(1, n):10 coef = 1.0 / (b[i] - a[i] * c[i - 1])11 c[i] = coef * c[i]12 d[i] = coef * (d[i] - a[i] * d[i - 1])1314for i in range(n - 2, -1, -1):15 d[i] = d[i] - c[i] * d[i + 1]1617return d1819 @staticmethod20 def Simplified_Natural4(y1, y2, y3, y4):21 # 四点⾃然边界化简公式22 d1 = y1 + y3 - 2 * y223 d2 = y2 + y4 - 2 * y32425 k0 = 026 k1 = (4 * d1 - d2) * 0.427 k2 = (4 * d2 - d1) * 0.428 k3 = 02930return [k0, k1, k2, k3]3132 @staticmethod33 def Simplified_Not_A_Knot4(y1, y2, y3, y4):34 # 四点⾮节点边界化简公式35 d1 = y1 + y3 - 2 * y236 d2 = y2 + y4 - 2 * y33738 k0 = 2 * d1 - d239 k1 = d140 k2 = d241 k3 = 2 * d2 - d14243return [k0, k1, k2, k3]4445 # TDMA矩阵说明46 # a0 和 c3 没有实际意义，占位⽤47 # a0 [b0 c0 00 ] [x0] [d0]48 # [a1 b1 c1 0 ] [x1] = [d1]49 # [0 a2 b2 c2] [x2] [d2]50 # [00 a3 b3] c3 [x3] [d3]5152 def SplineInterpolationNatural4(self, x, y1, y2, y3, y4):53 # ⽤TDMA计算54 # matrix_a = [0, 1, 1, 0]55 # matrix_b = [1, 4, 4, 1]56 # matrix_c = [0, 1, 1, 0]57 # matrix_d = [0, 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 0]58 # matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)5960 # 化简计算61 matrix_x = self.Simplified_Natural4(y1, y2, y3, y4)6263 a = y264 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.065 c = matrix_x[1] / 2.066 d = (matrix_x[2] - matrix_x[1]) / 6.06768 s = a + b * x + c * x * x + d * x * x * x69return s7071 def SplineInterpolationClamped4(self, x, y1, y2, y3, y4):72 # 仅有TDMA计算，⽆法化简73 A, B = 1, 17475 matrix_a = [0, 1, 1, 1]76 matrix_b = [2, 4, 4, 2]77 matrix_c = [1, 1, 1, 0]78 matrix_d = [6 * (y2 - y1 - A), 6 * (y1 + y3 - 2 * y2), 6 * (y2 + y4 - 2 * y3), 6 * (B - y4 + y3)]79 matrix_x = self.TDMA(matrix_a, matrix_b, matrix_c, matrix_d)8081 a = y282 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.083 c = matrix_x[1] / 2.084 d = (matrix_x[2] - matrix_x[1]) / 6.08586 s = a + b * x + c * x * x + d * x * x * x87return s8889 def SplineInterpolationNotAKnot4(self, x, y1, y2, y3, y4):90 # ⽆法使⽤TDMA计算91 matrix_x = self.Simplified_Not_A_Knot4(y1, y2, y3, y4)9293 a = y294 b = y3 - y2 - matrix_x[1] / 3.0 - matrix_x[2] / 6.095 c = matrix_x[1] / 2.096 d = (matrix_x[2] - matrix_x[1]) / 6.09798 s = a + b * x + c * x * x + d * x * x * x99return s100101 def biSpline4(self, srcImg, dstH, dstW):102 dstH, dstW = int(dstH), int(dstW)103 srcH, srcW, _ = srcImg.shape104 srcImg = np.pad(srcImg, ((1, 2), (1, 2), (0, 0)), 'edge')105 dstImg = np.zeros((dstH, dstW, 3), dtype=np.uint8)106for dstY in range(dstH):107for dstX in range(dstW):108for channel in [0, 1, 2]:109 # p11 p12 p13 p14110 #111 # p21 p22 p23 p24112 # p113 # p31 p32 p33 p34114 #115 # p41 p42 p43 p44116 # 储存为 p(y, x)117 p = [dstY * srcH / dstH, dstX * srcW / dstW]118 p22 = [math.floor(p[0]), math.floor(p[1])]119 p21 = [p22[0], p22[1] - 1]120 p23 = [p22[0], p22[1] + 1]121 p24 = [p22[0], p22[1] + 2]122123 p11 = [p21[0] - 1, p21[1]]124 p12 = [p11[0], p22[1]]125 p13 = [p11[0], p23[1]]126 p14 = [p11[0], p24[1]]127128 p31 = [p21[0] + 1, p21[1]]129 p32 = [p31[0], p22[1]]130 p33 = [p31[0], p23[1]]131 p34 = [p31[0], p24[1]]132133 p41 = [p21[0] + 2, p21[1]]134 p42 = [p41[0], p22[1]]135 p43 = [p41[0], p23[1]]136 p44 = [p41[0], p24[1]]137138 diff_y, diff_x = p[0] - p22[0], p[1] - p22[1]139 r1 = self.SplineInterpolationNatural4(diff_x, srcImg[p11[0], p11[1], channel], srcImg[p12[0], p12[1], channel], srcImg[p13[0], p13[1], channel], srcImg[p14[0], p14[1], channel]) 140 r2 = self.SplineInterpolationNatural4(diff_x, srcImg[p21[0], p21[1], channel], srcImg[p22[0], p22[1], channel], srcImg[p23[0], p23[1], channel], srcImg[p24[0], p24[1], channel]) 141 r3 = self.SplineInterpolationNatural4(diff_x, srcImg[p31[0], p31[1], channel], srcImg[p32[0], p32[1], channel], srcImg[p33[0], p33[1], channel], srcImg[p34[0], p34[1], channel]) 142 r4 = self.SplineInterpolationNatural4(diff_x, srcImg[p41[0], p41[1], channel], srcImg[p42[0], p42[1], channel], srcImg[p43[0], p43[1], channel], srcImg[p44[0], p44[1], channel]) 143144 c = self.SplineInterpolationNatural4(diff_y, r1, r2, r3, r4)145146 dstImg[dstY, dstX, channel] = np.clip(c, 0, 255)。
计算方法大作业——三次样条插值
计算方法上机报告
此完成所有数据的输入。继续按 Enter 键会出现提示“选择封闭方程组的边界条件: 第 一类边界条件输入 1，第二类边界条件输入 2，第三类边界条件输入 3。 ”根据已知情况 选择相应的边界条件，若为自然三次样条插值，则选 1，并将插值区间两端点的二阶导 数值设置为 0。输入完成之后按 Enter 开始求解，程序运行结束后命令窗口会显示要求 的三次样条插值函数，同时会出现该插值函数以及插值节点的图像，便于直接观察。 2.3 算例及计算结果 (1) 《数值分析》课本第 137 页的例题 4.6.1，已知函数 y=f(x)的数值如下表，求它 的自然三次样条插值函数。 xi yi -3 7 -1 11 0 26 3 56 4 29
(2) 给定函数 f ( x)
3 x 1 1 x 0 0 x3 3 x 4
1 (1 x 1) 。取等距节点，构造牛顿插值多项式 N5(x) 1 25x 2 和 N10(x)及三次样条插值函数 S10(x)。分别将三种插值多项式与 f(x)的曲线画在同一个
N10 x
22757 10 5444 8 20216 6 17147 4 3725 2 x x x x x 1 103 11 53 139 221
将牛顿插值多项式 N5(x)和 N10(x)及三次样条插值函数 S10(x)分别与 f(x)的曲线画在 同一个坐标系上进行比较，如图 12。可以看出三次样条函数与原函数符合的非常好， 对于低次的牛顿插值多项式，与原函数的大致趋势相同，而高次的牛顿插值多项式由 于龙格现象的出现，与原函数之间相差比较大。
S ( xi ) S ( xi ), ( xi ) S ( xi ), S S ( x ) S ( x ), i i i 1, 2, , n 1
三次样条插值
三次样条插值
三次样条插值是插值运算的一种，它具有计算精度高、收敛性好以及曲线拟合准确等特点，是插值运算中最常用的插值方法之
三次样条插值是以曲线为基本元素，把离散点数据连接成一个曲线，并能够在曲线上求出任意点的函数值。
它通过拟合所有离散数据点，来求出一个连续曲线，从而解决了插值法的局限性。
三次样条插值的基本原理是：在离散点的两端，曲线的曲率是零，由此可以计算出曲线的系数，从而得到曲线的表达式，这样就可以得到曲线上任意点的函数值。
三次样条插值的优点在于计算精度高、收敛性好，可以很好地拟合离散数据，并且经过插值后得到的曲线更加平滑，其结果更加可靠。
由于它的优点，三次样条插值得到了广泛的应用，如在统计分析中，用于拟合离散数据；在机械工程中，用于优化加工轨迹；在号处理中，用于滤波等。
总之，三次样条插值是插值运算的一种，它的准确性高，拟合性好，广泛应用于各种领域，是科学研究中的一种重要方法。
三次样条插值的方法和思路 -回复
三次样条插值的方法和思路-回复三次样条插值是一种常用的插值方法，它可以在已知的离散数据点上构造出一条光滑的曲线。
这种方法被广泛应用在曲线拟合、图像处理、数据分析等领域。
本文将介绍三次样条插值的方法和思路，并详细阐述每个步骤。
第一步是确定插值段数。
在进行三次样条插值时，首先需要将已知数据点划分成若干个插值段。
插值段越多，插值曲线越接近原始数据，但也会使插值算法复杂度增加。
因此，在确定插值段数时需要权衡精度和计算效率。
第二步是计算每个插值段的系数。
对于每个插值段，我们需要计算出一个三次曲线，该曲线会通过该段的两个端点。
具体的计算方法是，假设有n 个插值点，则有n-1个插值段，每个插值段的系数需要通过以下步骤计算：1. 计算边界条件：这是三次样条插值的关键一步。
我们需要根据已知数据点的性质，来确定边界条件是自然边界、固定边界还是其他类型的边界。
自然边界要求二阶导数在两个端点处为0，即S''(x_0) = S''(x_n) = 0。
固定边界要求插值曲线通过端点的给定导数值，即S'(x_0) = d_0、S'(x_n) = d_n。
2. 构建三对角矩阵：三次样条插值的求解过程可以转化为解线性方程组的问题。
为了解这个方程组，我们需要构建一个三对角矩阵。
其中的对角线元素是2，上下对角线元素是1。
3. 计算方程组的右侧：方程组的右侧是一个n-1维的向量，每个元素对应插值段的边界条件。
对于自然边界，右侧元素都是0；对于固定边界，则通过求解给定的导数值得到。
4. 解线性方程组：将三对角矩阵与右侧向量相乘，即可得到每个插值段的系数。
第三步是构造插值曲线。
在前两步中，我们计算得到了每个插值段的系数。
现在，我们需要将这些系数整合起来，构造出整个插值曲线。
具体的构造方法为，对于第i个插值段，其插值函数可表示为：S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3其中x_i和x_{i+1}为插值段的端点，a_i、b_i、c_i、d_i为第i个插值段的系数。
三次样条曲线插补改进算法
三次样条曲线插补改进算法朱宁【摘要】重点研究了三次样条曲线的插补算法及其改进,解决了现有样条曲线插补算法存在的弓高难控制、加工存在安全隐患的问题.改进后的算法简单,适合现代有强大计算能力的计算机系统,低成本情况下能在低档数控机床上直接加工样条曲线.【期刊名称】《通信电源技术》【年(卷),期】2018(035)004【总页数】3页(P29-30,32)【关键词】三次曲线;插补;样条曲线【作者】朱宁【作者单位】徐州机电技师学院,江苏徐州 221000【正文语种】中文0 引言三次样条曲线是在生产实践中产生和发展起来的。
在CAD/CAM技术还没有得到广泛应用时，技术人员绘制飞机、船舶和汽车上的复杂轮廓曲线都是借助于样条通过手工来完成[1]。
绘制样条曲线时，先选定好支点位置，然后在另一端放上重物或压铁使其作自由弹性弯曲，获得的曲线即所需要的样条曲线。
1 三次样条曲线函数的一般插补算法1.1 三次样条曲线函数的一般定义已知n个点P1(x1,y1)，P2(x2,y2)，…，Pn(xn,yn)，且x1＜x2…＜xn，若函数S(x)满足条件：（1）曲线通过所有的型值点，即S(xi)=yi(i=1,2,…n)；（2）S(x)在[xi,xn]区间上有连续的一阶和二阶导数；（3）S(x)在每一个子区间[xi,xi+1]上都是三次多项式，即每一个子区间内有Si(x)=Ai+Bi(x-xi)+Ci(xxi)2+Di(x-xi)3,(i=1,2，…n-1)。
则称S(x)为[xi,xn]上以xi(i=1,2,…n)为结点的三次样条函数。
1.2 三次样条曲线的常见插补运算令t为弦长参数，x=x(t)，y=y(t)。
可以看出：对应于n个型值点(xi,yi)(i=1,2,…n)有n个弦长参数t i(i=1,2,…n)，如图1所示。
图1 弦长参数示意图令t1=0，t2=[(x2-x1)2+(y2-y1)2]1/2，…，tn=[(xn-x1)2+(yn-y1)2]1/2。
三次样条插值计算算法
/* 三次样条插值计算算法*/#include "math.h "#include "stdio.h "#include "stdlib.h "/*N：已知节点数N+1R：欲求插值点数R+1x，y为给定函数f(x)的节点值{x(i)} (x(i) <x(i+1)) ,以及相应的函数值{f(i)} 0 <=i <=NP0=f(x0)的二阶导数；Pn=f(xn)的二阶导数u:存插值点{u(i)} 0 <=i <=R求得的结果s(ui)放入s[R+1] 0 <=i <=R返回0表示成功，1表示失败*/int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[]){/*声明局部变量*/double *h; /*存放步长:{hi} 0 <=i <=N-1 */double *a; /*存放系数矩阵{ai} 1 <=i <=N ；分量0没有利用*/ double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0 <=i <=N-1 */double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0 <=i <=N */double *af; /*存放系数矩阵{a(f)i} 1 <=i <=N ；*/double *ba; /*存放中间结果0 <=i <=N-1*/double *m; /*存放方程组的解{m(i)} 0 <=i <=N ；*/int i,k;double p1,p2,p3,p4;/*分配空间*/if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);/*第一步：计算方程组的系数*/for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++)a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++)c[k]=1-a[k];for(k=1;k <N;k++)g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]); c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;/*第二步：用追赶法解方程组求{m(i)} */ba[0]=c[0]/2;g[0]=g[0]/2;for(i=1;i <N;i++){af[i]=2-a[i]*ba[i-1];g[i]=(g[i]-a[i]*g[i-1])/af[i];ba[i]=c[i]/af[i];}af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N]; /*P110 公式：6.32*/ for(i=N-1;i> =0;i--)m[i]=g[i]-ba[i]*m[i+1];/*第三步：求值*/for(i=0;i <=R;i++){/*判断u(i)属于哪一个子区间，即确定k */if(u[i] <x[0] || u[i]> x[N]){/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 1;}k=0;while(u[i]> x[k+1])k++;//p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); //p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p1=(h[k]+2*(u[i]-x[k]))*pow((u[i]-x[k+1]),2)*y[k]/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1]))*pow((u[i]-x[k]),2)*y[k+1]/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 0;}void main(){int N,R;double *x,*y,*u,*s;double P0,Pn;int i;/*验证算法：*/N=7;R=6;/*分配空间*/if(!(x=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(y=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(u=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(s=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0 .9463;u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;P0=-0.4794;Pn=-0.9463;if(!SPL( N, R, x, y, P0, Pn, u, s)){/*打印结果*/printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\ns= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf( "\nsin= ");for(i=0;i <=R;i++)printf( "%9.5f ",sin(u[i]));}/*释放空间*/free(x);free(y);free(u);free(s);}/* 测试数据来自课本55页例5 《数值分析》清华大学出版社第四版*/ //输入327.7 4.128 4.329 4.130 3.013.0 -4.0//输出输出三次样条插值函数：1: [27.7 , 28]13.07*(x - 28)^3 + 0.22*(x - 27.7)^3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: [28 , 29]0.066*(29 - x)^3 + 0.1383*(x - 28)^3+ 4.234*(29 - x) + 3.962*(x - 28)3: [29 , 30]0.1383*(30 - x)^3 - 1.519*(x - 29)^3+ 3.962*(30 - x) + 4.519*(x - 29)//三次样条插值函数#include<iostream>#include<iomanip>using namespace std;const int MAX = 50;float x[MAX], y[MAX], h[MAX];float c[MAX], a[MAX], fxym[MAX];float f(int x1, int x2, int x3){float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);return (a - b)/(x[x3] - x[x1]);} //求差分void cal_m(int n){ //用追赶法求解出弯矩向量M……float B[MAX];B[0] = c[0] / 2;for(int i = 1; i < n; i++)B[i] = c[i] / (2 - a[i]*B[i-1]);fxym[0] = fxym[0] / 2;for(i = 1; i <= n; i++)fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);for(i = n-1; i >= 0; i--)fxym[i] = fxym[i] - B[i]*fxym[i+1];}void printout(int n);int main(){int n,i; char ch;do{cout<<"Please put in the number of the dots:";cin>>n;for(i = 0; i <= n; i++){cout<<"Please put in X"<<i<<':';cin>>x[i]; //cout<<endl;cout<<"Please put in Y"<<i<<':';cin>>y[i]; //cout<<endl;}for(i = 0; i < n; i++) //求步长h[i] = x[i+1] - x[i];cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";int t;float f0, f1;cin>>t;switch(t){case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";cin>>f0>>f1;c[0] = 1; a[n] = 1;fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];break;case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";cin>>f0>>f1;c[0] = a[n] = 0;fxym[0] = 2*f0; fxym[n] = 2*f1;break;default:cout<<"不可用\n";//待定};//switchfor(i = 1; i < n; i++)fxym[i] = 6 * f(i-1, i, i+1);for(i = 1; i < n; i++){a[i] = h[i-1] / (h[i] + h[i-1]);c[i] = 1 - a[i];}a[n] = h[n-1] / (h[n-1] + h[n]);cal_m(n);cout<<"\n输出三次样条插值函数：\n";printout(n);cout<<"Do you to have anther try ? y/n :";cin>>ch;}while(ch == 'y' || ch == 'Y');return 0;}void printout(int n){cout<<setprecision(6);for(int i = 0; i < n; i++){cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";/*cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";cout<<endl;*/float t = fxym[i]/(6*h[i]);if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";else cout<<-t<<"*(x - "<<x[i+1]<<")^3";t = fxym[i+1]/(6*h[i]);if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";cout<<"\n\t";t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";cout<<endl<<endl;}cout<<endl;}。
详细讲解三次样条插值法及其实现方法61页PPT
39、没有不老的誓言，没有不变的承 诺，踏 上旅途 ，义无 反顾。 40、对时间的价值没有没有深切认识 的人， 决不会 坚韧勤 勉。
66、节制使快乐增加并使享受加强。 ——德 谟克利 特 67、今天应做的事没有做，明天再早也 是耽误 了。——裴斯 泰洛齐 68、决定一个人的一生，以及整个命运 的，只 是一瞬 之间。 ——歌 德 69、懒人无法享受休息之乐。——拉布 克 70、浪费时间是一桩大罪过。——卢梭
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
关键词
数据采样插补，三次三角插值样条曲线，刀具运动轨迹，轨迹空间
1. 引言
在数控加工中，经常遇到如飞机的的机翼、汽车流线型覆盖件成形模具型腔、汽轮机叶片等许多复 杂外形型面的零件。面对这类曲线曲面的描述通常采用的方法是 NURBS 方法，它既能表示自由曲线曲 面，又能表示一些规则的二次曲线曲面。但 NURBS 曲线的参数表达式比较复杂，插补运算量大，且由 型值点逆求控制点的反求运算较难实现，不易于实现直接插补。为此，本课题组在文献[1]中提出了一种 新型样条曲线——可调形三次三角插值样条曲线，这种样条曲线的插值样条表达式中含有三角函数多项 式，通过选择合适的控制点和形状可调参数，该插值样条可以精确表示直线、圆弧、椭圆以及自由曲线 等常见工程曲线。由于该使用的三次三角插值样条避免了使用有理的数学表达形式，故其表达式比较简 洁，所需的计算量也相对较小，并且能够实现实时插补，这很好的为插值曲线及曲面的设计，提供了一 种全新的理论方法。 本文针对该三次三角样条曲线提出了一种新的插补计算方法，该法有效地提高了加工精度和加工过 程的可靠性。其插补思想是：使用若干顺序连接的微小直线段逼近用户程序给定的插值样条曲线轮廓。 将直线段映射到参变量空间，得到参变量的增量值，再通过控制参变量 u，求取对应的轨迹空间中的映 射点，最终得到插补点的坐标值。以此不断重复直到插补终点，可得到整个离散化插补轨迹。
Mechanical Engineering and Technology 机械工程与技术, 2015, 4, 151-158 Published Online June 2015 in Hans. /journal/met /10.12677/met.2015.42016
扬州大学机械工程学院，江苏 扬州 Email: ydzhangyihan@
151
三次三角插值样条曲线的数控插补计算方法
收稿日期：2015年6月3日；录用日期：2015年6月22日；发布日期：2015年6月25日
摘
Hale Waihona Puke 要论文提出一种新的样条曲线——可调形三次三角插值样条曲线，该样条曲线实现了对直线、圆弧、椭圆 以及自由曲线等常见数控运动轨迹曲线的统一精确表达，并根据参数曲线的数据采样插补原理，对三次 三角插值样条曲线的插补算法进行了研究。在保证自由曲线数控插补运动高速同时兼顾运动的平稳性， 采用一系列首尾相连的微小直线段逼近给定的插值样条曲线， 再利用轨迹空间和参变量空间的对应关系， 控制运动加速度，得到整个离散化插补轨迹。此算法提高了自由曲线的插补速度和运动平稳性，保证了 曲面数控加工的表面质量。
(4)
2 2 2 2 π k π k −1 1− π k − π k − 2 k 1 − 2 k 2 k − 1 2 k π π π π 2 2 0 1− k 0 k π π M = 2 2 k 0 1− k 0 π π 2 2 0 0 k − k π π 2 2 k 0 − k 0 π π 0 0 0 0
Keywords
Data Sampling Interpolation Theory, Cubic Trigonometric Interpolation Spline, Tool Motion Path Curves, Track Space
三次三角插值样条曲线的数控插补计算方法
张益汉，宋爱平，刘祖奇，邱 林
j =0
3
qi −1 B0,3 ( u ) + qi B1,3 ( u ) + qi +1 B2,3 ( u ) + qi + 2 B3,3 ( u )
(2)
= i 1, 2, , n − 2 )，为由每四个顺序控制点 qi−1，qi，qi+1，qi+2 所构造的分 ( u ∈ [ 0,1] ，i 为曲线段号，
p1 , p2 , , pi , 的坐标值。3) 将直线映射到参变量空间，得到与其对应的参变量空间中的微小直线段即
参变量的增量值，求取与参变量空间中下一点插补点的坐标。
4. 三次三角插值样条曲线的插补的实现方法
4.1. 插补预处理
为了减少三次三插值样条曲线在插补中的计算量，可将插补算法中一些可以一次性计算的任务放在 实时插补前的预处理阶段完成。这些任务主要包括轨迹表达系数的求取和导数表达式系数的求取[2]。 轨迹表达式 P ( u ) 系数的求取是根据 CAD/CAM 给出的样条曲线描述信息， 在实时插补前一次性求解
3.2. 三次三角插值样条曲线插补实现原理和步骤
三次三角样条曲线的插补过程的直接控制量是参变量 u，而最终的被控量是插补点的坐标值以及插 补点所沿插补轨迹的进给速度，先利用轨迹空间到参变量空间的映射和参变量，然后再通过参变量空间 到轨迹空间的映射这两大步骤的相互转换来完成。具体步骤如下：1) 根据三次三角插值样条曲线现有的 几何信息，确定轨迹计算公式中的相关系数。2) 按照用户程序的进给速度、加减速要求以及允许误差， 在 插 补 周 期 中 产 生 空 间 小 直 线 段 ∆l1 , ∆l2 , , ∆li , 去 逼 近 被 插 补 曲 线 ， 逐 步 求 得 各 插 补 直 线 端 点
段三次三角样条曲线，所有曲线段 Pi(u)有序组成的曲线 P(u)称为三次三角样条曲线。
3. 三次三角插值样条曲线的插补计算方法
数据采样插补实质上是使用一系列首尾相连的微小直线段来逼近给定曲线。由于这些线段是按加工 时间分割的，因此数据采样插补又称为时间分割法插补，近代数控系统轨迹插补在数据采样原理的基础 上产生了很多新的插补算法。本文就数据采样插补原理产生的三次三角插值样条曲线的插补算法进行研 究。
(1)
152
三次三角插值样条曲线的数控插补计算方法
其中 S := sin ( πu 2 ) ， C := cos ( πu 2 ) 。 设 qt ( t = 0,1, , n ) 为 R2 或 R3 中给定的一组控制点(n ≥ 3)，定义参数曲线段
Pi ( u )=
∑ qi −1+ j B j ,3 ( u )=
P ( u ) = aS 3 + bC 3 + cS 2 + dC 2 + eS + fC
式中：a，b，c， d，e，f 为常数矢量， u ∈ [ 0,1] ， S = sin ( πu 2 ) ， C = cos ( πu 2 )
(3)
各项系数将随插值点位置的变化而变化，即取不同的值，故三次三角插值样条曲线是一个变系数三 次参数曲线插补问题。
rd nd th
Abstract
The paper proposes a new spline curve called cubic trigonometric cardinal interpolation spline with adjustable shape. The spline curve realized the unification to accurately represent some common NC tool motion path curves such as straight lines, arcs, ellipses and free curves, etc. An interpolation algorithm of cubic trigonometric cardinal interpolation spline based on data sampling interpolation theory is researched. This curve interpolation method ensures the high-speed of spline curve CNC interpolation while taking into account the smooth movement, using a series of conterminous short straight-line segments to approach the given interpolation spline, and then utilizing the relation between track space and parameter space to control motion acceleration and get the whole discrete interpolation trajectory. This algorithm improves the speed and motion stability of free curve interpolation and guarantees the surface quality of surface numerical control machining.
An NC Interpolation Calculation Method of Cubic Trigonometric Cardinal Interpolation Spline
Yihan Zhang, Aiping Song, Zuqi Liu, Lin Qiu
College of Mechanical Engineering, Yangzhou University, Yangzhou Jiangsu Email: ydzhangyihan@ Received: Jun. 3 , 2015; accepted: Jun. 22 , 2015; published: Jun. 25 , 2015 Copyright © 2015 by authors and Hans Publishers Inc. This work is licensed under the Creative Commons Attribution International License (CC BY). /licenses/by/4.0/