4.4 插值interpolate模块

插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。与拟合不同的是，要求曲线通过所有的已知数据。SciPy的interpolate模块提供了许多对数据进行插值运算的函数，范围涵盖简单的一维插值到复杂多维插值求解。当样本数据变化归因于一个独立的变量时，就使用一维插值；反之样本数据归因于多个独立变量时，使用多维插值。
计算插值有两种基本的方法，1、对一个完整的数据集去拟合一个函数；2、对数据集的不同部分拟合出不同的函数，而函数之间的曲线平滑对接。第二种方法又叫做仿样内插法，当数据拟合函数形式非常复杂时，这是一种非常强大的工具。我们首先介绍怎样对简单函数进行一维插值运算，然后进一步深入比较复杂的多维插值运算。

4.4.1 一维插值

一维数据的插值运算可以通过函数interp1d()完成。其调用形式如下，它实际上不是函数而是一个类：

interp1d(x, y, kind='linear', ...)

其中，x和y参数是一系列已知的数据点，kind参数是插值类型，可以是字符串或整数，它给出插值的B样条曲线的阶数，候选值及作用下表所示：

候选值	作用
‘zero’ 、'nearest'	阶梯插值，相当于0阶B样条曲线
‘slinear’ 、'linear'	线性插值，用一条直线连接所有的取样点，相当于一阶B样条曲线
‘quadratic’ 、'cubic'	二阶和三阶B样条曲线，更高阶的曲线可以直接使用整数值指定

下面的程序演示了通过不同的 kind参数（linear和quadratic），对一个正弦函数进行插值运算。示例代码：

import numpy as np
from scipy.interpolate import interp1d

#创建待插值的数据
x = np.linspace(0, 10*np.pi, 20)
y = np.cos(x)

# 分别用linear和quadratic插值
fl = interp1d(x, y, kind='linear')
fq = interp1d(x, y, kind='quadratic')

#设置x的最大值和最小值以防止插值数据越界
xint = np.linspace(x.min(), x.max(), 1000)
yintl = fl(xint)
yintq = fq(xint)

通过绘图的方式查看插值拟合效果，示例代码：

import pylab as pl
pl.plot(xint,fl(xint), color="green", label = "Linear")
pl.plot(xint,fq(xint), color="yellow", label ="Quadratic")
pl.legend(loc = "best")
pl.show()

绿色线表示线性插值曲线，黄色线表示二阶插值曲线，如图所示：

线性插值和二阶插值效果

4.4.2 噪声数据插值

我们可以通过interpolate模块中UnivariateSpline()函数对含有噪声的数据进行插值运算，示例代码：

import numpy as np
from scipy.interpolate import UnivariateSpline

# 通过人工方式添加噪声数据
sample = 30
x = np.linspace(1, 10*np.pi, sample)
y = np.cos(x) + np.log10(x) + np.random.randn(sample)/10

# 插值，参数s为smoothing factor 
f = UnivariateSpline(x, y, s=1)
xint = np.linspace(x.min(), x.max(), 1000)
yint = f(xint)

在UnivariateSpline()函数中，参数s是平滑向量参数，被用来拟合还有噪声的数据。如果参数s=0，将忽略噪声对所有点进行插值运算。通过绘图的方式查看插值效果，示例代码：

import pylab as pl
pl.plot(xint,f(xint), color="green", label = "Interpolation")
pl.plot(x, y, color="yellow", label ="Original")
pl.legend(loc = "best")
pl.show()

绿色线表示插值曲线，黄色线表示原始曲线，如图所示：

噪声数据插值拟合

4.4.3 多维插值

多维插值主要用于重构图片，interpolate模块中的griddata()函数有很强大的处理多维散列取样点进行插值运算的能力。其调用形式如下：

griddata(points, values, xi, method='linear', fill_value=nan)

其中points表示K维空间中的坐标，它可以是形状为(N,k)的数组，也可以是一个有k个数组的序列，N为数据的点数。values是points中每个点对应的值。xi是需要进行插值运算的坐标，其形状为(M,k)。method参数有三个选项：'nearest'、 ‘linear’、 'cubic'，分别对应0阶、1阶以及3阶插值。
以下示例利用1000个随机散列点对1000x1000像素的图片进行重构，示例代码：

import numpy as np
from scipy.interpolate import griddata

#定义一个函数
ripple = lambda x,y:np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

# 生成grid数据，复数定义了生成grid数据的step，若无该复数则step为5 
grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j] 

# 生成待插值的样本数据 
xy = np.random.rand(1000, 2) 
sample = ripple(xy[:,0] * 5 , xy[:,1] * 5) 

# 用cubic方法插值
grid_z0 = griddata(xy * 5, sample, (grid_x, grid_y), method='cubic')

我们还可以使用interpolate模块的SmoothBivariateSpline类进行多元仿样插值运算，对图片进行重构。示例代码：

import numpy as np
from scipy.interpolate import SmoothBivariateSpline as SBS

# 定义函数
ripple = lambda x,y:np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)

# 生成插值样本
xy = np.random.rand(1000, 2)
x, y = xy[:, 0], xy[:, 1]
sample = ripple(xy[:, 0] * 5, xy[:, 1] * 5)

# 插值运算
fit = SBS(x * 5, y * 5, sample, s=0.01, kx=4, ky=4)
interp = fit(np.linspace(0, 5, 1000), np.linspace(0, 5, 1000))

我们得到了一个与上个示例同样的结果。左图显示的是随机样本点的原始图像，右图是插值重构图像。除了图右上角部分，整体上SmoothBivariateSpline函数的表现略好于griddata函数。
通过反复测试，尽管SmoothBivariateSpline表现略好，但其对给定的样本数据非常敏感，就可能导致忽略一些显著特征。而griddata函数有很强的鲁棒性，不管给定的数据样本，能够合理的进行插值运算。