一维插值
线性插值
- 线性插值就是将相邻两点用直线连接起来
- 用线性插值进行近似计算,当插值区间小时,近似程度较高。
%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)))处也这样做。
那么:
从特殊到一般
1.定义插值基函数:
基函数性质
(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.拉格朗日插值多项式:
#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个方程:
2.第一个端点和最后一个端点必须过第一个和最后一个方程:这里一共是2个方程
3.节点处的一阶导数的值必须相等。这里为两个方程。
4.在这里假设第一个方程的二阶导数为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.端点处的二阶导数为零,这里是两个方程。
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实现
4.三次样条插值法