1.5. Scipy:高级科学计算

作者Adrien Chauve,Andre Espaze,Emmanuelle Gouillart,GaëlVaroquaux,Ralf Gommers

Scipy

scipy包包含专用于科学计算中常见问题的各种工具箱。其不同的子模块对应于不同的应用,如插值、积分、优化、图像处理、统计、特殊功能等。

scipy可以与其他标准科学计算库进行比较,例如GSL(GNU Scientific Library for C and C++)或Matlab的工具箱。scipy是Python中用于科学计算程序的核心包;它的意图之一是有效地操作numpy数组,以便numpy和scipy可以同步工作。

在实现一个例程之前,值得检查所需的数据处理是否尚未在Scipy中实现。作为非专业程序员,科学家通常倾向于重新创造轮子,这会导致错误、非最佳、难以共享和无法维护的代码。相比之下,Scipy的例程已优化和测试,因此应尽可能使用。

警告

本教程远不是数值计算的介绍。因为在scipy中枚举不同的子模块和函数将是非常无聊,所以我们集中精力在几个例子上给出如何使用scipy进行科学计算的一般想法。

scipy由任务特定的子模块组成:

scipy.cluster 矢量量化/Kmeans
scipy.constants 物理和数学常数
scipy.fftpack 傅里叶变换
scipy.integrate 集成例程
scipy.interpolate 插值
scipy.io 数据输入和输出
scipy.linalg 线性代数例程
scipy.ndimage n维图像包
scipy.odr 正交距离回归
scipy.optimize 优化
scipy.signal 信号处理
scipy.sparse 稀疏矩阵
scipy.spatial 空间数据结构和算法
scipy.special 任何特殊的数学函数
scipy.stats 统计

它们都依赖于numpy,但是大多数是彼此独立的。导入Numpy和这些Scipy模块的标准方法是:

>>> import numpy as np
>>> from scipy import stats # same for other sub-modules

scipy主命名空间包含的函数大部分真实为numpy函数(试一下scipy.cos is np.cos)。它们暴露出来只是历史原因;通常没有理由在代码中使用import scipy

1.5.1. 文件输入/输出: scipy.io

  • 加载和保存matlab文件:

    >>> from scipy import io as spio
    
    >>> a = np.ones((3, 3))
    >>> spio.savemat('file.mat', {'a': a}) # savemat expects a dictionary
    >>> data = spio.loadmat('file.mat', struct_as_record=True)
    >>> data['a']
    array([[ 1., 1., 1.],
    [ 1., 1., 1.],
    [ 1., 1., 1.]])
  • 读取图片:

    >>> from scipy import misc
    
    >>> misc.imread('fname.png')
    array(...)
    >>> # Matplotlib also has a similar function
    >>> import matplotlib.pyplot as plt
    >>> plt.imread('fname.png')
    array(...)

另见:

1.5.2. 特殊函数: scipy.special

特殊的函数是超越数函数。scipy.special模块的文档字符串编写得很好,因此我们不会在此列出所有函数。常用的有:

  • Bessel函数,如scipy.special.jn()(第n个整数阶贝塞尔函数)
  • 椭圆函数(scipy.special.ellipj()用于雅可比椭圆函数,...)
  • 伽马函数:scipy.special.gamma(),还要注意scipy.special.gammaln(),它将Gamma的log值赋予更高的数值精度。
  • Erf,高斯曲线下方的面积:scipy.special.erf()

1.5.3. 线性代数运算: scipy.linalg

scipy.linalg模块提供标准线性代数运算,依赖于底层有效的实现(BLAS,LAPACK)。

  • scipy.linalg.det()函数计算方阵的行列式:

    >>> from scipy import linalg
    
    >>> arr = np.array([[1, 2],
    ... [3, 4]])
    >>> linalg.det(arr)
    -2.0
    >>> arr = np.array([[3, 2],
    ... [6, 4]])
    >>> linalg.det(arr)
    0.0
    >>> linalg.det(np.ones((3, 4)))
    Traceback (most recent call last):
    ...
    ValueError: expected square matrix
  • scipy.linalg.inv()函数计算方阵的逆:

    >>> arr = np.array([[1, 2],
    
    ... [3, 4]])
    >>> iarr = linalg.inv(arr)
    >>> iarr
    array([[-2. , 1. ],
    [ 1.5, -0.5]])
    >>> np.allclose(np.dot(arr, iarr), np.eye(2))
    True

    最后计算奇异矩阵的逆矩阵(其行列式为零)将产生LinAlgError

    >>> arr = np.array([[3, 2],
    
    ... [6, 4]])
    >>> linalg.inv(arr)
    Traceback (most recent call last):
    ...
    ...LinAlgError: singular matrix
  • 还有更高级的操作,例如奇异值分解(SVD):

    >>> arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])
    
    >>> uarr, spec, vharr = linalg.svd(arr)

    得到的阵列谱为:

    >>> spec    
    
    array([ 14.88982544, 0.45294236, 0.29654967])

    可以通过svd的输出与np.dot的输出的矩阵乘法来重组原始矩阵:

    >>> sarr = np.diag(spec)
    
    >>> svd_mat = uarr.dot(sarr).dot(vharr)
    >>> np.allclose(svd_mat, arr)
    True

    SVD通常用于统计和信号处理。许多其他标准分解(QR,LU,Cholesky,Schur)以及线性系统的求解器在scipy.linalg中提供。

1.5.4. 快速傅里叶变换: scipy.fftpack

scipy.fftpack模块允许计算快速傅立叶变换。作为示例,(噪声)输入信号可以看起来像:

>>> time_step = 0.02
>>> period = 5.
>>> time_vec = np.arange(0, 20, time_step)
>>> sig = np.sin(2 * np.pi / period * time_vec) + \
... 0.5 * np.random.randn(time_vec.size)

观察者不知道信号频率,仅知道信号的采样时间步长sig该信号应来自实函数,因此傅里叶变换将是对称的。scipy.fftpack.fftfreq()函数将生成采样频率,scipy.fftpack.fft()将计算快速傅立叶变换:

>>> from scipy import fftpack
>>> sample_freq = fftpack.fftfreq(sig.size, d=time_step)
>>> sig_fft = fftpack.fft(sig)

因为所得到的功率是对称的,所以仅需要使用频谱的正部分来找到频率:

>>> pidxs = np.where(sample_freq > 0)
>>> freqs = sample_freq[pidxs]
>>> power = np.abs(sig_fft)[pidxs]

[源代码hires.pngpdf]

../_images/fftpack_frequency.png

信号频率可以通过以下公式找到:

>>> freq = freqs[power.argmax()]
>>> np.allclose(freq, 1./period) # check that correct freq is found
True

现在高频噪声将从傅立叶变换信号中去除:

>>> sig_fft[np.abs(sample_freq) > freq] = 0

可以通过scipy.fftpack.ifft()函数计算得到的滤波信号:

>>> main_sig = fftpack.ifft(sig_fft)

结果可以查看:

>>> import pylab as plt
>>> plt.figure()
<matplotlib.figure.Figure object at 0x...>
>>> plt.plot(time_vec, sig)
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.plot(time_vec, main_sig, linewidth=3)
[<matplotlib.lines.Line2D object at 0x...>]
>>> plt.xlabel('Time [s]')
<matplotlib.text.Text object at 0x...>
>>> plt.ylabel('Amplitude')
<matplotlib.text.Text object at 0x...>

[源代码hires.pngpdf]

../_images/fftpack_signals.png

numpy.fft

Numpy还具有FFT(numpy.fft)的实现。然而,一般来说,scipy应该是首选,因为它使用更有效的底层实现。

工作实例:粗周期性发现

[源代码hires.pngpdf]

../_images/periodicity_finder_00.png

[源代码hires.pngpdf]

../_images/periodicity_finder_01.png

工作示例:高斯图像模糊

卷积:

[源代码hires.pngpdf]

../_images/image_blur.png

练习:去噪月亮着陆图像

../_images/moonlanding.png
  1. 检查提供的图像moonlanding.png,它被周期性噪声严重污染。在本练习中,我们旨在使用快速傅立叶变换来清除噪声。
  2. 使用pylab.imread()加载图片。
  3. scipy.fftpack中找到并使用2-D FFT函数,并绘制图像的光谱(傅里叶变换)。你对光谱的可视化有什么麻烦吗?如果是,为什么?
  4. 频谱由高频和低频分量组成。噪声包含在频谱的高频部分,因此将某些分量设置为零(使用阵列切片)。
  5. 应用傅里叶逆变换以查看生成的图像。

1.5.5. 优化和拟合: scipy.optimize

优化是找到最小化或相等的数值解的问题。

scipy.optimize模块为函数最小化(标量或多维),曲线拟合和根查找提供了有用的算法。

>>> from scipy import optimize

查找标量函数的最小值

让我们定义以下函数:

>>> def f(x):
... return x**2 + 10*np.sin(x)

并绘制:

>>> x = np.arange(-10, 10, 0.1)
>>> plt.plot(x, f(x))
>>> plt.show()

[源代码hires.pngpdf]

../_images/scipy_optimize_example1.png

该函数的全局最小值约为-1.3,局部最小值约为3.8。

为该函数找到最小值的一般和有效的方法是从给定的初始点开始进行梯度下降。BFGS算法是一个很好的方法:

>>> optimize.fmin_bfgs(f, 0)
Optimization terminated successfully.
Current function value: -7.945823
Iterations: 5
Function evaluations: 24
Gradient evaluations: 8
array([-1.30644003])

这种方法的一个可能的问题是,如果函数具有局部最小值,则算法可以根据初始点找到这些局部最小值而不是全局最小值:

>>> optimize.fmin_bfgs(f, 3, disp=0)
array([ 3.83746663])

如果我们不知道全局最小值的邻域来选择初始点,我们需要诉诸于昂贵的全局优化。为了找到全局最小值,我们使用scipy.optimize.basinhopping()(它结合了局部优化器和本地优化器的起始点的随机抽样):

版本0.12.0中的新功能: basinhopping在版本0.12.0的Scipy中添加

>>> optimize.basinhopping(f, 0)  
nfev: 1725
minimization_failures: 0
fun: -7.9458233756152845
x: array([-1.30644001])
message: ['requested number of basinhopping iterations completed successfully']
njev: 575
nit: 100

另一个可用的(但效率低得多的)全局优化器是scipy.optimize.brute()(网格上的强力优化)。存在用于不同类别的全局优化问题的更高效的算法,但是这超出了scipy的范围。用于全局优化的一些有用的包是OpenOpt,IPOPTPyGMOPyEvolve

注意

scipy用于包含例程退火,其自SciPy 0.14.0起已被弃用,并在SciPy 0.16.0中被移除。

To find the local minimum, let’s constraint the variable to the interval (0, 10) using scipy.optimize.fminbound():

>>> xmin_local = optimize.fminbound(f, 0, 10)
>>> xmin_local
3.8374671...

注意

在高级章节中更详细地讨论了函数的最小值:Mathematical optimization: finding minima of functions

查找标量函数的根

为了找到根,即上面的函数ff(x) = 0点,我们可以使用例如scipy.optimize.fsolve()

>>> root = optimize.fsolve(f, 1)  # our initial guess is 1
>>> root
array([ 0.])

注意,只找到一个根。检查f的图表明有大约-2.5的第二根。我们通过调整我们的初始猜测找到它的确切值:

>>> root2 = optimize.fsolve(f, -2.5)
>>> root2
array([-2.47948183])

曲线拟合

假设我们从f采样了一些噪声的数据:

>>> xdata = np.linspace(-10, 10, num=20)
>>> ydata = f(xdata) + np.random.randn(xdata.size)

现在,如果我们知道绘制样本的函数的函数形式(在这种情况下x^2 + \sin(x)),而不是项的幅度,我们可以通过最小二乘曲线拟合找到那些。首先我们要定义函数来拟合:

>>> def f2(x, a, b):
... return a*x**2 + b*np.sin(x)

然后,我们可以使用scipy.optimize.curve_fit()来查找ab

>>> guess = [2, 2]
>>> params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)
>>> params
array([ 0.99667386, 10.17808313])

现在我们已经找到了f的最小值和根,并对它使用了曲线拟合,我们将所有的resuls放在一个图中:

[源代码hires.pngpdf]

../_images/scipy_optimize_example2.png

注意

In Scipy >= 0.11 unified interfaces to all minimization and root finding algorithms are available: scipy.optimize.minimize(), scipy.optimize.minimize_scalar() and scipy.optimize.root(). 它们允许通过method关键字轻松比较各种算法。

您可以在scipy.optimize中找到多维问题具有相同功能的算法。

练习:温度数据的曲线拟合

从1月开始,每月的阿拉斯加气温极值由(摄氏度)给出:

max:  17,  19,  21,  28,  33,  38, 37,  37,  31,  23,  19,  18
min: -62, -59, -56, -46, -32, -18, -9, -13, -25, -46, -52, -58
  1. 绘制这些温度极值。
  2. 定义一个可以描述最小和最大温度的函数。提示:此功能必须有1年的期限。提示:包括时间偏移。
  3. 使用scipy.optimize.curve_fit()将此函数适用于数据。
  4. 绘制结果。适合是否合理?如果不是,为什么?
  5. 在拟合精度内,最小和最大温度的时间偏移是否相同?

练习:2-D最小化

[源代码hires.pngpdf]

../_images/scipy_optimize_sixhump.png

六驼峰回弹功能

有多个全局和局部最小值。找到此函数的全局最小值。

提示:

有多少个全局最小值,这些点的函数值是什么?(x, y) = (0, 0)的初始猜测会发生什么?

请参阅Non linear least squares curve fitting: application to point extraction in topographical lidar data中的总结练习。

1.5.6. 统计资料和随机数字: scipy.stats

模块scipy.stats包含统计工具和随机过程的概率描述。用于各种随机过程的随机数发生器可以在numpy.random中找到。

1.5.6.1. Histogram and probability density function

给定随机过程的观察,它们的直方图是随机过程的PDF(概率密度函数)的估计:

>>> a = np.random.normal(size=1000)
>>> bins = np.arange(-4, 5)
>>> bins
array([-4, -3, -2, -1, 0, 1, 2, 3, 4])
>>> histogram = np.histogram(a, bins=bins, normed=True)[0]
>>> bins = 0.5*(bins[1:] + bins[:-1])
>>> bins
array([-3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5])
>>> from scipy import stats
>>> b = stats.norm.pdf(bins) # norm is a distribution
>>> plt.plot(bins, histogram)
[<matplotlib.lines.Line2D object at ...>]
>>> plt.plot(bins, b)
[<matplotlib.lines.Line2D object at ...>]

[源代码hires.pngpdf]

../_images/normal_distribution.png

如果我们知道随机过程属于一个给定的随机过程家族,例如正常过程,我们可以做最大似然拟合的观察值来估计基础分布的参数。这里我们对观察数据拟合正态过程:

>>> loc, std = stats.norm.fit(a)
>>> loc
0.0314345570...
>>> std
0.9778613090...

练习:概率分布

从形状参数为1的伽马分布生成1000个随机变量,然后从这些样本绘制直方图。你可以在顶部绘制pdf(它应该匹配)?

额外:分布有很多有用的方法。通过阅读docstring或使用IPython选项卡完成来探索它们。你可以使用随机变量上的fit方法找到1的形状参数吗?

1.5.6.2. Percentiles

中值是具有下面一半观察值的值,并且是上面的一半:

>>> np.median(a)     
0.04041769593...

它也称为百分位数50,因为50%的观察值低于它:

>>> stats.scoreatpercentile(a, 50)     
0.0404176959...

同样,我们可以计算百分位数90:

>>> stats.scoreatpercentile(a, 90)     
1.3185699120...

百分位数是CDF的累积分布函数的估计量。

1.5.6.3. Statistical tests

统计测试是决策指标。例如,如果我们有两组观察值,我们假设从高斯过程产生,我们可以使用T-test来决定两组观察值是否有明显的不同:

>>> a = np.random.normal(0, 1, size=100)
>>> b = np.random.normal(1, 1, size=10)
>>> stats.ttest_ind(a, b)
(array(-3.177574054...), 0.0019370639...)

结果输出包括:

  • T统计值:它是一个数,其符号与两个随机过程之间的差成比例,并且该量值与该差的显着性相关。
  • p值:两个进程的概率相同。如果接近1,这两个过程几乎肯定是相同的。越接近零,过程具有不同手段的可能性越大。

另见

有关statistics的章节介绍了更多精细的工具,用于统计测试和统计数据加载以及在scipy之外的可视化。

1.5.7. 插值: scipy.interpolate

scipy.interpolate适用于从实验数据拟合函数,从而评估不存在度量的点。该模块基于netlib项目中的FITPACK Fortran子程序

通过想象接近正弦函数的实验数据:

>>> measured_time = np.linspace(0, 1, 10)
>>> noise = (np.random.random(10)*2 - 1) * 1e-1
>>> measures = np.sin(2 * np.pi * measured_time) + noise

scipy.interpolate.interp1d类可以构建一个线性插值函数:

>>> from scipy.interpolate import interp1d
>>> linear_interp = interp1d(measured_time, measures)

然后,scipy.interpolate.linear_interp实例需要在感兴趣的时间进行评估:

>>> computed_time = np.linspace(0, 1, 50)
>>> linear_results = linear_interp(computed_time)

也可以通过提供kind可选关键字参数来选择三次插值:

>>> cubic_interp = interp1d(measured_time, measures, kind='cubic')
>>> cubic_results = cubic_interp(computed_time)

结果现在收集在以下Matplotlib图:

[源代码hires.pngpdf]

../_images/scipy_interpolation.png

scipy.interpolate.interp2d类似于scipy.interpolate.interp1d,但适用于2-D数组。请注意,对于interp系列,计算的时间必须保持在测量的时间范围内。有关更为先进的样条插值示例,请参见Maximum wind speed prediction at the Sprogø station上的总结练习。

1.5.8. 数值积分: scipy.integrate

最通用的集成例程是scipy.integrate.quad()

>>> from scipy.integrate import quad
>>> res, err = quad(np.sin, 0, np.pi/2)
>>> np.allclose(res, 1)
True
>>> np.allclose(err, 1 - res)
True

其他集成方案可以使用fixed_quadquadratureromberg

scipy.integrate还具有用于集成普通微分方程(ODE)的例程。特别地,scipy.integrate.odeint()是使用LSODA的通用积分器(Livermore Solver for Ordinary Differential equation with Automatic method switching for stiff and non-rigid problems),参见ODEPACK Fortran库了解详情。

odeint解决以下形式的一阶ODE系统:

dy/dt = rhs(y1, y2, .., t0,...)

作为介绍,让我们用初始条件y(t=0) = 1解决t = 0 \dots 4之间的ODE \frac{dy}{dt} = -2 y首先,需要定义计算位置的导数的函数:

>>> def calc_derivative(ypos, time, counter_arr):
... counter_arr += 1
... return -2 * ypos
...

添加了一个额外的参数counter_arr以说明该函数可以针对单个时间步长被调用多次,直到求解器收敛。计数器阵列定义为:

>>> counter = np.zeros((1,), dtype=np.uint16)

现在将计算轨迹:

>>> from scipy.integrate import odeint
>>> time_vec = np.linspace(0, 4, 40)
>>> yvec, info = odeint(calc_derivative, 1, time_vec,
... args=(counter,), full_output=True)

因此,导数函数已被调用超过40次(这是时间步长数):

>>> counter
array([129], dtype=uint16)

并且10个第一时间步长中的每一个的迭代的累积数目可以通过以下获得:

>>> info['nfe'][:10]
array([31, 35, 43, 49, 53, 57, 59, 63, 65, 69], dtype=int32)

注意,求解器对于第一时间步骤需要更多迭代。现在可以绘制轨迹的解:yvec

scipy.integrate.odeint()的另一个示例是阻尼弹簧质量振荡器(2阶振荡器)。The position of a mass attached to a spring obeys the 2nd order ODE y'' + 2 \varepsilon \omega_0 y' + \omega_0^2 y = 0 with \omega_0^2 = k/m with k the spring constant, m the mass and \varepsilon = c/(2 m \omega_0) with c the damping coefficient. 对于这个例子,我们选择参数:

>>> mass = 0.5  # kg
>>> kspring = 4 # N/m
>>> cviscous = 0.4 # N s/m

因此系统将处于欠阻尼,因为:

>>> eps = cviscous / (2 * mass * np.sqrt(kspring/mass))
>>> eps < 1
True

对于scipy.integrate.odeint()求解器,二阶方程需要在矢量Y = (y, y')的两个一阶方程的系统中变换。定义\nu = 2 \varepsilon * \omega_0 = c / m\Omega = \omega_0^2 = k/m将很方便:

>>> nu_coef = cviscous / mass  # nu
>>> om_coef = kspring / mass # Omega

因此,该函数将通过以下公式计算速度和加速度:

>>> def calc_deri(yvec, time, nu, om):
... return (yvec[1], -nu * yvec[1] - om * yvec[0])
...
>>> time_vec = np.linspace(0, 10, 100)
>>> yinit = (1, 0)
>>> yarr = odeint(calc_deri, yinit, time_vec, args=(nu_coef, om_coef))

最终位置和速度显示在下面的Matplotlib图上:

[源代码hires.pngpdf]

../_images/odeint_damped_spring_mass.png

这两个例子只是普通微分方程(ODE)。然而,在Scipy中没有部分微分方程(PDE)求解器。用于解决PDE的一些Python包可用,例如fipySfePy

1.5.9. 信号处理: scipy.signal

>>> from scipy import signal

1.5.10. 图片处理: scipy.ndimage

在scipy中专用于图像处理的子模块是scipy.ndimage

>>> from scipy import ndimage

图像处理例程可以根据它们执行的处理的类别来排序。

1.5.10.1. Geometrical transformations on images

更改方向,分辨率,..

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> shifted_face = ndimage.shift(face, (50, 50))
>>> shifted_face2 = ndimage.shift(face, (50, 50), mode='nearest')
>>> rotated_face = ndimage.rotate(face, 30)
>>> cropped_face = face[50:-50, 50:-50]
>>> zoomed_face = ndimage.zoom(face, 2)
>>> zoomed_face.shape
(1536, 2048)
../_images/face_transforms.png
>>> plt.subplot(151)    
<matplotlib.axes._subplots.AxesSubplot object at 0x...>
>>> plt.imshow(shifted_face, cmap=plt.cm.gray)
<matplotlib.image.AxesImage object at 0x...>
>>> plt.axis('off')
(-0.5, 1023.5, 767.5, -0.5)
>>> # etc.

1.5.10.2. Image filtering

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> face = face[:512, -512:] # crop out square on right
>>> import numpy as np
>>> noisy_face = np.copy(face).astype(np.float)
>>> noisy_face += face.std() * 0.5 * np.random.standard_normal(face.shape)
>>> blurred_face = ndimage.gaussian_filter(noisy_face, sigma=3)
>>> median_face = ndimage.median_filter(noisy_face, size=5)
>>> from scipy import signal
>>> wiener_face = signal.wiener(noisy_face, (5, 5))
../_images/filtered_face.png

scipy.ndimage.filtersscipy.signal中的许多其他过滤器可应用于图像。

行使

比较不同滤波图像的直方图。

1.5.10.3. Mathematical morphology

数学形态是源于集合理论的数学理论。它表征和变换几何结构。特别地,可以使用该理论来转换二进制(黑白)图像:要转换的集合是相邻非零值像素的集合。该理论也扩展到灰度图像。

../_images/morpho_mat1.png

基本数学形态运算使用结构化元素以修改其他几何结构。

让我们首先生成一个结构化元素

>>> el = ndimage.generate_binary_structure(2, 1)
>>> el
array([[False, True, False],
[...True, True, True],
[False, True, False]], dtype=bool)
>>> el.astype(np.int)
array([[0, 1, 0],
[1, 1, 1],
[0, 1, 0]])
  • 侵蚀

    >>> a = np.zeros((7, 7), dtype=np.int)
    
    >>> a[1:6, 2:5] = 1
    >>> a
    array([[0, 0, 0, 0, 0, 0, 0],
    [0, 0, 1, 1, 1, 0, 0],
    [0, 0, 1, 1, 1, 0, 0],
    [0, 0, 1, 1, 1, 0, 0],
    [0, 0, 1, 1, 1, 0, 0],
    [0, 0, 1, 1, 1, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]])
    >>> ndimage.binary_erosion(a).astype(a.dtype)
    array([[0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]])
    >>> #Erosion removes objects smaller than the structure
    >>> ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)
    array([[0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0]])
  • Dilation

    >>> a = np.zeros((5, 5))
    
    >>> a[2, 2] = 1
    >>> a
    array([[ 0., 0., 0., 0., 0.],
    [ 0., 0., 0., 0., 0.],
    [ 0., 0., 1., 0., 0.],
    [ 0., 0., 0., 0., 0.],
    [ 0., 0., 0., 0., 0.]])
    >>> ndimage.binary_dilation(a).astype(a.dtype)
    array([[ 0., 0., 0., 0., 0.],
    [ 0., 0., 1., 0., 0.],
    [ 0., 1., 1., 1., 0.],
    [ 0., 0., 1., 0., 0.],
    [ 0., 0., 0., 0., 0.]])
  • 打开

    >>> a = np.zeros((5, 5), dtype=np.int)
    
    >>> a[1:4, 1:4] = 1
    >>> a[4, 4] = 1
    >>> a
    array([[0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 0, 0, 1]])
    >>> # Opening removes small objects
    >>> ndimage.binary_opening(a, structure=np.ones((3, 3))).astype(np.int)
    array([[0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 0, 0, 0]])
    >>> # Opening can also smooth corners
    >>> ndimage.binary_opening(a).astype(np.int)
    array([[0, 0, 0, 0, 0],
    [0, 0, 1, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 1, 0, 0],
    [0, 0, 0, 0, 0]])
  • 关闭: ndimage.binary_closing

行使

检查开口量是否侵蚀,然后扩张。

打开操作移除小结构,而关闭操作填充小孔。因此,这样的操作可以用于“清洁”图像。

>>> a = np.zeros((50, 50))
>>> a[10:-10, 10:-10] = 1
>>> a += 0.25 * np.random.standard_normal(a.shape)
>>> mask = a>=0.5
>>> opened_mask = ndimage.binary_opening(mask)
>>> closed_mask = ndimage.binary_closing(opened_mask)
../_images/morpho.png

行使

检查重建的正方形的面积是否小于初始正方形的面积。(如果在开启之前执行关闭步骤,则会发生相反的情况)。

对于灰度值图像,扩张)等于用最小(相应地,最大)值,其中该像素由以感兴趣像素为中心的结构化元素覆盖。

>>> a = np.zeros((7, 7), dtype=np.int)
>>> a[1:6, 1:6] = 3
>>> a[4, 4] = 2; a[2, 3] = 1
>>> a
array([[0, 0, 0, 0, 0, 0, 0],
[0, 3, 3, 3, 3, 3, 0],
[0, 3, 3, 1, 3, 3, 0],
[0, 3, 3, 3, 3, 3, 0],
[0, 3, 3, 3, 2, 3, 0],
[0, 3, 3, 3, 3, 3, 0],
[0, 0, 0, 0, 0, 0, 0]])
>>> ndimage.grey_erosion(a, size=(3, 3))
array([[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 3, 2, 2, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]])

1.5.10.4. Measurements on images

让我们首先生成一个不错的合成二进制图像。

>>> x, y = np.indices((100, 100))
>>> sig = np.sin(2*np.pi*x/50.) * np.sin(2*np.pi*y/50.) * (1+x*y/50.**2)**2
>>> mask = sig > 1

现在我们寻找关于图像中的对象的各种信息:

>>> labels, nb = ndimage.label(mask)
>>> nb
8
>>> areas = ndimage.sum(mask, labels, range(1, labels.max()+1))
>>> areas
array([ 190., 45., 424., 278., 459., 190., 549., 424.])
>>> maxima = ndimage.maximum(sig, labels, range(1, labels.max()+1))
>>> maxima
array([ 1.80238238, 1.13527605, 5.51954079, 2.49611818,
6.71673619, 1.80238238, 16.76547217, 5.51954079])
>>> ndimage.find_objects(labels==4)
[(slice(30L, 48L, None), slice(30L, 48L, None))]
>>> sl = ndimage.find_objects(labels==4)
>>> import pylab as pl
>>> pl.imshow(sig[sl[0]])
<matplotlib.image.AxesImage object at ...>
../_images/measures.png

有关更高级的示例,请参阅Image processing application: counting bubbles and unmolten grains中的摘要练习。

1.5.11. 科学计算概要练习

总结练习主要使用Numpy,Scipy和Matplotlib。他们用Python提供了一些科学计算的现实例子。现在介绍了使用Numpy和Scipy的基础知识,感兴趣的用户被邀请尝试这些练习。

练习:

建议的解决方案: