插值法

一维插值

线性插值

• 线性插值就是将相邻两点用直线连接起来
• 用线性插值进行近似计算，当插值区间小时，近似程度较高。

%pylab
%matplotlib inline
from scipy import interpolate 
#解决中文显示问题
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)

x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y)

plt.title(u"线性插值")

plt.plot(x, y, "o")
plt.plot(x_new, f_linear(x_new))
plt.show()

多项式插值

用多项式$p(x) = a_0 + a_1 * x + a_2 * x^2 + ... + a_n * x^n $拟合

%pylab
%matplotlib inline
from scipy import interpolate 

x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)

x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y, kind = 'quadratic')

plt.title(u"多项式插值")

plt.plot(x, y, "o")
plt.plot(x_new, f_linear(x_new))
plt.show()

Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib

拉格朗日插值

拉格朗日插值是用多条特殊的多项式函数，确定相应的系数，加权得到多项式函数

一个例子

已知下面这几个点，我想找到一根穿过它们的曲线：

使用多项式是可以的， 但是大神的脑回路总是令人捉摸不透，拉格朗日是这么想的，第一，这肯定是多项式曲线，第二，他认为可以用多条曲线叠加得到，具体是这样的。
第一根曲线 (f_1(x)) ，在(x_1)点处，取值为1，其余两点取值为0：

同理在 ((x_2，f_2(x))) ，((x_3，f_3(x)))处也这样做。

那么：

[f(x) = y_1*f_1(x) + y_2*f_2x) + y_3*f_3(x) ]

从特殊到一般

1.定义插值基函数：

[l_k(x) = frac{(x-x_0)...(x-x_k-1)(x-x_k+1)...(x-x_n)}{(x_k-x_0)...(x_k-x_k-1)(x_k-x_k+1)...(x_k-x_n)} = prod_{j = 0, j eq k}^{n} frac{x - x_j}{x_k - x_j} ]

基函数性质

(1)
( l_k(x_i) = egin{Bmatrix} 1, i=k \ 0, i eq k end{Bmatrix} (k = 0, 1 ... n) )

(2)(l_k(x_i))是惟一确定的n次多项式

(3)拉格朗日插值所包含基函数个数与插值节点个数相同。

2.拉格朗日插值多项式：

[L_n(x) = sum_{k=0}^{n}y_k l_k(x) = sum_{k=0}^{n}y_k prod_{j = 0, j eq k}^{n} frac{x - x_j}{x_k - x_j} ]

#coding=utf-8
from matplotlib import pyplot as plt

def Lg(data,testdata):
    predict=0
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    if testdata in data_x:
        #print "testdata is already known"
        return data_y[data_x.index(testdata)]
    for i in range(len(data_x)):
        af=1
        for j in range(len(data_x)):
            if j!=i:
                af*=(1.0*(testdata-data_x[j])/(data_x[i]-data_x[j]))
        predict+=data_y[i]*af
    return predict

def plot(data,nums):
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    Area=[min(data_x),max(data_x)]
    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]
    X[len(X)-1]=Area[1]
    Y=[Lg(data,x) for x in X]
    plt.title(u"拉格朗日插值")
    plt.plot(X,Y)
    for i in range(len(data_x)):
        plt.plot(data_x[i],data_y[i],'ro',label="point")
    plt.show()

data=[[0,0],[1,2],[2,3],[3,8],[4,2]]
print (Lg(data,1.5))
plot(data,100)

1.84375

牛顿(Newton)插值

牛顿Newton插值基本思想是将待求的n次插值多项式(P_n（x）)改写为具有承袭性的形式，然后利用插值条件⑴确定(P_n（x）)的待定系数，以求出所要的插值函数。

一个例子

用牛顿向前插值多项式计算

简单的推导

一阶均差：

二阶均差是一阶均差的均差：

三阶均差就是二阶均差的均差，以此类推，我们得到牛顿插值法为：

联系泰勒公式理解吧，大名鼎鼎的泰勒公式就是这么来的。

#coding=utf-8
from matplotlib import pyplot as plt

def calF(data):
    #差商计算  n个数据 0-(n-1)阶个差商 n个数据
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    F= [1 for i in range(len(data))]   
    FM=[]
    for i in range(len(data)):
        FME=[]
        if i==0:
            FME=data_y
        else:
            for j in range(len(FM[len(FM)-1])-1):
                delta=data_x[i+j]-data_x[j]
                value=1.0*(FM[len(FM)-1][j+1]-FM[len(FM)-1][j])/delta
                FME.append(value)
        FM.append(FME)
    F=[fme[0] for fme in FM]
    print (FM)
    return F

def NT(data,testdata,F):
    #差商之类的计算
    predict=0
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    if testdata in data_x:
        return data_y[data_x.index(testdata)]
    else:
        for i in range(len(data_x)):
            Eq=1
            if i!=0:
                for j in range(i):
                    Eq=Eq*(testdata-data_x[j])
            predict+=(F[i]*Eq)
        return predict

def plot(data,nums):
    data_x=[data[i][0] for i in range(len(data))]
    data_y=[data[i][1] for i in range(len(data))]
    Area=[min(data_x),max(data_x)]
    X=[Area[0]+1.0*i*(Area[1]-Area[0])/nums for i in range(nums)]
    X[len(X)-1]=Area[1]
    F=calF(data)
    Y=[NT(data,x,F) for x in X]
    plt.title(u"牛顿(Newton)插值")
    plt.plot(X,Y,label='result')
    for i in range(len(data_x)):
        plt.plot(data_x[i],data_y[i],'ro',label="point")
    plt.show()

data=[[0,0],[1,2],[2,3],[3,8],[4,2]]

plot(data,100)

[[0, 2, 3, 8, 2], [2.0, 1.0, 5.0, -6.0], [-0.5, 2.0, -5.5], [0.8333333333333334, -2.5], [-0.8333333333333334]]

分段低次插值

二次插条

我们用二次函数：(aX^2+bx+c)来描述曲线，下图中一共有4个点，可以分成3个区间。每一个区间都需要一个二次函数来描述，一共需要9个未知数。下面的任务就是找出9个方程。

如下图所示：一共有(x_0,x_1,x_2,x_3)四个点，三个区间，每个区间上都有一个方程。

1.曲线方程在节点处的值必须相等，即函数在(x_1,x_2)两个点处的值必须符合两个方程，这里一共是4个方程：

[a_1*x_1^2+b_1*x_1+c_1=f(x_1) \ a_2*x_1^2+b_2*x_1+c_2=f(x_1) \ a_2*x_2^2+b_2*x_2+c_2=f(x_2) \ a_3*x_2^2+b_3*x_2+c_3=f(x_2) \ ]

2.第一个端点和最后一个端点必须过第一个和最后一个方程：这里一共是2个方程

3.节点处的一阶导数的值必须相等。这里为两个方程。

[2*a_1*x_1+b_1 = 2*a_2*x_1+b_2 \ 2*a_2*x_1+b_2 = 2*a_3*x_2+b_3 \ ]

4.在这里假设第一个方程的二阶导数为0：这里为一个方程：

[a_1=0 ]

import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
二次样条实现：
函数x的自变量为:3,   4.5, 7,    9
      因变量为：2.5, 1   2.5,  0.5
"""
x = [3, 4.5, 7, 9]
y = [2.5, 1, 2.5, 0.5]
 
"""一共有三个区间，用二次样条求解，需要有9个方程"""
 
 
"""
功能：完后对二次样条函数求解方程参数的输入
参数：要进行二次样条曲线计算的自变量
返回值：方程的参数
"""
def calculateEquationParameters(x):
    #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
    parameter = []
    sizeOfInterval=len(x)-1;
    i = 1
    #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
    while i < len(x)-1:
        data = init(sizeOfInterval*3)
        data[(i-1)*3]=x[i]*x[i]
        data[(i-1)*3+1]=x[i]
        data[(i-1)*3+2]=1
        data1 =init(sizeOfInterval*3)
        data1[i * 3] = x[i] * x[i]
        data1[i * 3 + 1] = x[i]
        data1[i * 3 + 2] = 1
        temp=data[1:]
        parameter.append(temp)
        temp=data1[1:]
        parameter.append(temp)
        i += 1
    #输入端点处的函数值。为两个方程,加上前面的2n-2个方程，一共2n个方程
    data = init(sizeOfInterval*3-1)
    data[0] = x[0]
    data[1] = 1
    parameter.append(data)
    data = init(sizeOfInterval *3)
    data[(sizeOfInterval-1)*3+0] = x[-1] * x[-1]
    data[(sizeOfInterval-1)*3+1] = x[-1]
    data[(sizeOfInterval-1)*3+2] = 1
    temp=data[1:]
    parameter.append(temp)
    #端点函数值相等为n-1个方程。加上前面的方程为3n-1个方程,最后一个方程为a1=0总共为3n个方程
    i=1
    while i < len(x) - 1:
        data = init(sizeOfInterval * 3)
        data[(i - 1) * 3] =2*x[i]
        data[(i - 1) * 3 + 1] =1
        data[i*3]=-2*x[i]
        data[i*3+1]=-1
        temp=data[1:]
        parameter.append(temp)
        i += 1
    return parameter
 
"""
对一个size大小的元组初始化为0
"""
def init(size):
    j = 0;
    data = []
    while j < size:
        data.append(0)
        j += 1
    return data
 
 
"""
功能：计算样条函数的系数。
参数：parametes为方程的系数，y为要插值函数的因变量。
返回值：二次插值函数的系数。
"""
 
def solutionOfEquation(parametes,y):
    sizeOfInterval = len(x) - 1;
    result = init(sizeOfInterval*3-1)
    i=1
    while i<sizeOfInterval:
        result[(i-1)*2]=y[i]
        result[(i-1)*2+1]=y[i]
        i+=1
    result[(sizeOfInterval-1)*2]=y[0]
    result[(sizeOfInterval-1)*2+1]=y[-1]
    a = np.array(calculateEquationParameters(x))
    b = np.array(result)
    return np.linalg.solve(a,b)
 
"""
功能：根据所给参数，计算二次函数的函数值：
参数:parameters为二次函数的系数，x为自变量
返回值：为函数的因变量
"""
def calculate(paremeters,x):
    result=[]
    for data_x in x:
        result.append(paremeters[0]*data_x*data_x+paremeters[1]*data_x+paremeters[2])
    return  result
 
 
"""
功能：将函数绘制成图像
参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
返回值：空
"""
def  Draw(data_x,data_y,new_data_x,new_data_y):
        plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")
        plt.scatter(data_x,data_y, label="离散数据",color="red")
        mpl.rcParams['font.sans-serif'] = ['SimHei']
        mpl.rcParams['axes.unicode_minus'] = False
        plt.title("二次样条函数")
        plt.legend(loc="upper left")
        plt.show()
 
result=solutionOfEquation(calculateEquationParameters(x),y)
new_data_x1=np.arange(3, 4.5, 0.1)
new_data_y1=calculate([0,result[0],result[1]],new_data_x1)
new_data_x2=np.arange(4.5, 7, 0.1)
new_data_y2=calculate([result[2],result[3],result[4]],new_data_x2)
new_data_x3=np.arange(7, 9.5, 0.1)
new_data_y3=calculate([result[5],result[6],result[7]],new_data_x3)
new_data_x=[]
new_data_y=[]
new_data_x.extend(new_data_x1)
new_data_x.extend(new_data_x2)
new_data_x.extend(new_data_x3)
new_data_y.extend(new_data_y1)
new_data_y.extend(new_data_y2)
new_data_y.extend(new_data_y3)
Draw(x,y,new_data_x,new_data_y)

三次插条

三次样条的原理和二次样条的原理相同，我们用函数(aX^3+bX^2+cX+d)这个函数来进行操作，这里一共是4个点，分为3个区间，每个区间一个三次样条函数的话，一共是12个方程，只要我们找出这12个方程，这个问题就算解决了。

1.内部节点处的函数值应该相等，这里一共是4个方程。

2.函数的第一个端点和最后一个端点，应该分别在第一个方程和最后一个方程中。这里是2个方程。

3.两个函数在节点处的一阶导数应该相等。这里是两个方程。

4.两个函数在节点处的二阶导数应该相等，这里是两个方程。

5.端点处的二阶导数为零，这里是两个方程。

[a_1=0  \ b_1=0 \ ]

import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
三次样条实现：
函数x的自变量为:3,   4.5, 7,    9
      因变量为：2.5, 1   2.5,  0.5
"""
x = [3, 4.5, 7, 9]
y = [2.5, 1, 2.5, 0.5]
 
 
"""
功能：完后对三次样条函数求解方程参数的输入
参数：要进行三次样条曲线计算的自变量
返回值：方程的参数
"""
def calculateEquationParameters(x):
    #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
    parameter = []
    sizeOfInterval=len(x)-1;
    i = 1
    #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
    while i < len(x)-1:
        data = init(sizeOfInterval*4)
        data[(i-1)*4] = x[i]*x[i]*x[i]
        data[(i-1)*4+1] = x[i]*x[i]
        data[(i-1)*4+2] = x[i]
        data[(i-1)*4+3] = 1
        data1 =init(sizeOfInterval*4)
        data1[i*4] =x[i]*x[i]*x[i]
        data1[i*4+1] =x[i]*x[i]
        data1[i*4+2] =x[i]
        data1[i*4+3] = 1
        temp = data[2:]
        parameter.append(temp)
        temp = data1[2:]
        parameter.append(temp)
        i += 1
   # 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程，一共2n个方程
    data = init(sizeOfInterval * 4 - 2)
    data[0] = x[0]
    data[1] = 1
    parameter.append(data)
    data = init(sizeOfInterval * 4)
    data[(sizeOfInterval - 1) * 4 ] = x[-1] * x[-1] * x[-1]
    data[(sizeOfInterval - 1) * 4 + 1] = x[-1] * x[-1]
    data[(sizeOfInterval - 1) * 4 + 2] = x[-1]
    data[(sizeOfInterval - 1) * 4 + 3] = 1
    temp = data[2:]
    parameter.append(temp)
    # 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。
    i=1
    while i < sizeOfInterval:
        data = init(sizeOfInterval * 4)
        data[(i - 1) * 4] = 3 * x[i] * x[i]
        data[(i - 1) * 4 + 1] = 2 * x[i]
        data[(i - 1) * 4 + 2] = 1
        data[i * 4] = -3 * x[i] * x[i]
        data[i * 4 + 1] = -2 * x[i]
        data[i * 4 + 2] = -1
        temp = data[2:]
        parameter.append(temp)
        i += 1
    # 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零，为两个方程。总共为4n个方程。
    i = 1
    while i < len(x) - 1:
        data = init(sizeOfInterval * 4)
        data[(i - 1) * 4] = 6 * x[i]
        data[(i - 1) * 4 + 1] = 2
        data[i * 4] = -6 * x[i]
        data[i * 4 + 1] = -2
        temp = data[2:]
        parameter.append(temp)
        i += 1
    return parameter
 
 
 
"""
对一个size大小的元组初始化为0
"""
def init(size):
    j = 0;
    data = []
    while j < size:
        data.append(0)
        j += 1
    return data
 
"""
功能：计算样条函数的系数。
参数：parametes为方程的系数，y为要插值函数的因变量。
返回值：三次插值函数的系数。
"""
 
def solutionOfEquation(parametes,y):
    sizeOfInterval = len(x) - 1;
    result = init(sizeOfInterval*4-2)
    i=1
    while i<sizeOfInterval:
        result[(i-1)*2]=y[i]
        result[(i-1)*2+1]=y[i]
        i+=1
    result[(sizeOfInterval-1)*2]=y[0]
    result[(sizeOfInterval-1)*2+1]=y[-1]
    a = np.array(calculateEquationParameters(x))
    b = np.array(result)
    for data_x in b:
        print(data_x)
    return np.linalg.solve(a,b)
 
"""
功能：根据所给参数，计算三次函数的函数值：
参数:parameters为二次函数的系数，x为自变量
返回值：为函数的因变量
"""
def calculate(paremeters,x):
    result=[]
    for data_x in x:
        result.append(paremeters[0]*data_x*data_x*data_x+paremeters[1]*data_x*data_x+paremeters[2]*data_x+paremeters[3])
    return  result
 
 
"""
功能：将函数绘制成图像
参数：data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
返回值：空
"""
def  Draw(data_x,data_y,new_data_x,new_data_y):
        plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")
        plt.scatter(data_x,data_y, label="离散数据",color="red")
        mpl.rcParams['font.sans-serif'] = ['SimHei']
        mpl.rcParams['axes.unicode_minus'] = False
        plt.title("三次样条函数")
        plt.legend(loc="upper left")
        plt.show()
 
 
result=solutionOfEquation(calculateEquationParameters(x),y)
new_data_x1=np.arange(3, 4.5, 0.1)
new_data_y1=calculate([0,0,result[0],result[1]],new_data_x1)
new_data_x2=np.arange(4.5, 7, 0.1)
new_data_y2=calculate([result[2],result[3],result[4],result[5]],new_data_x2)
new_data_x3=np.arange(7, 9.5, 0.1)
new_data_y3=calculate([result[6],result[7],result[8],result[9]],new_data_x3)
new_data_x=[]
new_data_y=[]
new_data_x.extend(new_data_x1)
new_data_x.extend(new_data_x2)
new_data_x.extend(new_data_x3)
new_data_y.extend(new_data_y1)
new_data_y.extend(new_data_y2)
new_data_y.extend(new_data_y3)
Draw(x,y,new_data_x,new_data_y)

1.0
1.0
2.5
2.5
2.5
0.5
0.0
0.0
0.0
0.0

二维插值

方法与一维数据插值类似

演示二维插值

# -*- coding: utf-8 -*-
"""
演示二维插值。
"""
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl

def func(x, y):
    return (x+y)*np.exp(-5.0*(x**2 + y**2))

# X-Y轴分为15*15的网格
y,x= np.mgrid[-1:1:15j, -1:1:15j]

fvals = func(x,y) # 计算每个网格点上的函数值  15*15的值
print (len(fvals[0]))

#三次样条二维插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')

# 计算100*100的网格上的插值
xnew = np.linspace(-1,1,100)#x
ynew = np.linspace(-1,1,100)#y
fnew = newfunc(xnew, ynew)#仅仅是y值   100*100的值

# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
plt.figure(figsize = (20, 8))
pl.subplot(121)
im1=pl.imshow(fvals, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower")#pl.cm.jet
#extent=[-1,1,-1,1]为x,y范围  favals为
pl.colorbar(im1)

pl.subplot(122)
im2=pl.imshow(fnew, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()

# -*- coding: utf-8 -*-
"""
演示二维插值。
"""
# -*- coding: utf-8 -*-
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
from scipy import interpolate
import matplotlib.cm as cm
import matplotlib.pyplot as plt

def func(x, y):
    return (x+y)*np.exp(-5.0*(x**2 + y**2))

# X-Y轴分为20*20的网格
x = np.linspace(-1, 1, 20)
y = np.linspace(-1,1,20)
x, y = np.meshgrid(x, y)#20*20的网格数据

fvals = func(x,y) # 计算每个网格点上的函数值  15*15的值

fig = plt.figure(figsize=(20, 8))
#Draw sub-graph1
ax=plt.subplot(1, 2, 1,projection = '3d')
surf = ax.plot_surface(x, y, fvals, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
plt.colorbar(surf, shrink=0.5, aspect=5)#标注

#二维插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')#newfunc为一个函数

# 计算100*100的网格上的插值
xnew = np.linspace(-1,1,100)#x
ynew = np.linspace(-1,1,100)#y
fnew = newfunc(xnew, ynew)#仅仅是y值   100*100的值  np.shape(fnew) is 100*100
xnew, ynew = np.meshgrid(xnew, ynew)
ax2=plt.subplot(1, 2, 2,projection = '3d')
surf2 = ax2.plot_surface(xnew, ynew, fnew, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True)
ax2.set_xlabel('xnew')
ax2.set_ylabel('ynew')
ax2.set_zlabel('fnew(x, y)')
plt.colorbar(surf2, shrink=0.5, aspect=5)#标注

plt.show()

图像处理的应用

import cv2
import numpy as np
from matplotlib import pyplot as plt


img = np.uint8(np.random.randint(0,255,size=(5,5)))
height,width= img.shape


# 声明新的维度
new_dimension = (1000, 1000)

plt.subplot(231)
plt.title("SRC Image")
plt.imshow(img,cmap='seismic')

plt.subplot(232)
resized = cv2.resize(img, new_dimension, interpolation = cv2.INTER_NEAREST)
plt.title("INTER_NEAREST")
plt.imshow(resized,cmap='seismic')


plt.subplot(233)
resized = cv2.resize(img, new_dimension, interpolation = cv2.INTER_LINEAR)
plt.title("INTER_LINEAR")
plt.imshow(resized,cmap='seismic')


plt.subplot(234)
resized = cv2.resize(img, new_dimension, interpolation = cv2.INTER_AREA)
plt.title("INTER_AREA")
plt.imshow(resized,cmap='seismic')


plt.subplot(235)
resized = cv2.resize(img, new_dimension, interpolation = cv2.INTER_CUBIC)
plt.title("INTER_CUBIC")
plt.imshow(resized,cmap='seismic')


plt.subplot(236)
resized = cv2.resize(img, new_dimension, interpolation = cv2.INTER_LANCZOS4)
plt.title("INTER_LANCZOS4")
plt.imshow(resized,cmap='seismic')

plt.show()

参考

1.Lagrange、Newton、分段插值法及Python实现

2.如何直观地理解拉格朗日插值法？

3.牛顿插值的几何解释是怎么样的？

4.三次样条插值法

5.（数值分析）各种插值法的python实现