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

scipy可以与其他标准科学计算库进行比较，例如 GSL（用于 C 和 C++ 的 GNU 科学库）或 Matlab 的工具箱。scipy是 Python 中科学计算的核心包；它旨在有效地在numpy 数组上运行，以便 numpy 和 scipy 共同使用。

作为非专业程序员，科学家往往倾向于重新发明轮子，这会导致错误、非最优、难以共享和不可维护的代码。相比之下，Scipy经过优化和测试，因此应尽可能使用。

目录

• 文件输入/输出： scipy.io
• 特殊函数： scipy.special
• 线性代数运算： scipy.linalg
• 插值： scipy.interpolate
• 优化和拟合： scipy.optimize
• 统计和随机数： scipy.stats
• 数值积分： scipy.integrate
• 快速傅里叶变换： scipy.fftpack
• 信号处理： scipy.signal
• 图像处理： scipy.ndimage

scipy 非常丰富，这里我们只介绍一些重点部分，帮助我们了解如何将scipy用于科学计算。

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

它们都依赖于numpy，但大多是相互独立的。导入 Numpy 和这些 Scipy 模块的标准方法是：

>>> import numpy as np
>>> from scipy import stats  

主scipy命名空间大部分包含了真正的numpy函数(scipy.cos ， np.cos)。这些都是由于历史原因造成的;没有理由在代码中使用import scipy。

.文件输入/输出：scipy.io

Matlab 文件：加载和保存：

>>>

>>> from scipy import io as spio
>>> a = np.ones((3, 3))
>>> spio.savemat('file.mat', {'a': a}) 
>>> data = spio.loadmat('file.mat')
>>> data['a']
array([[1.,  1.,  1.],
       [1.,  1.,  1.],
       [1.,  1.,  1.]])

Python/Matlab 不匹配，例如

>>>

>>> a = np.ones(3)
>>> a
array([1.,  1.,  1.])
>>> spio.savemat('file.mat', {'a': a})
>>> spio.loadmat('file.mat')['a']
array([[1.,  1.,  1.]])

图像文件：读取图像：

>>>

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

此外

• 加载文本文件：numpy.loadtxt()/numpy.savetxt()
• 巧妙加载文本/csv文件： numpy.genfromtxt()/numpy.recfromcsv()
• 快速高效，但特定于 numpy 的二进制格式： numpy.save()/numpy.load()
• scikit-image 中更高级的图像输入/输出： skimage.io

.特殊函数：scipy.special

特殊函数是超越函数。scipy.special模块的文档写得很好，所以我们不会在这里列出所有的函数。常用的有：

贝塞尔函数，如scipy.special.jn()（nth integer order Bessel function）

椭圆函数（scipy.special.ellipj()对于雅可比椭圆函数，...）

Gamma function: scipy.special.gamma()，还要注意 scipy.special.gammaln()这将使 Gamma 的对数具有更高的数值精度。

Erf，高斯曲线下的面积： scipy.special.erf()

.线性代数运算：scipy.linalg

scipy.linalg模块提供标准的线性代数运算。

from scipy import linalg
>>> arr = np.array([[1, 2],
...                 [3, 4]])
>>> linalg.det(arr)
-
>>> arr = np.array([[3, 2],
...                 [6, 4]])
>>> linalg.det(arr) 

>>> 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. ],
       [ , -]])
>>> np.allclose(np.dot(arr, iarr), np.eye(2))
True

• 可以使用更高级的操作，例如奇异值分解 (SVD)：

arr = np.array([[3, 2], ... [6, 4]])

>>> linalg.inv(arr) Traceback (most recent call last): ... ...

LinAlgError: singular matrix

• 得到的数组谱为：

>>> spec array([, , ])

• 原始矩阵可以通过svdwith的输出的矩阵乘法重新组合

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.

.插值：scipy.interpolate

scipy.interpolate对于从实验数据拟合函数非常有用，因此可以对不存在测度的点进行评估。该模块基于FITPACK Fortran。

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

>>>

>>> measured_time = np.linspace(0, 1, )
>>> noise = (np.random.random()*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)

然后可以在感兴趣的点估算结果：

>>>

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

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

>>>

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

scipy.interpolate.interp2d类似于
scipy.interpolate.interp1d，但适用于二维数组。请注意，对于interp系列，插值点必须保持在给定数据点的范围内。

优化和拟合：scipy.optimize

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

该scipy.optimize模块提供函数最小化（标量或多维）、曲线拟合和求根的算法。

>>>

>>> from scipy import optimize

。曲线拟合

假设我们有一个正弦波数据，带有一些噪声：

>>>

>>> x_data = np.linspace(-5, 5, num=)
>>> y_data =  * np.sin( * x_data) + np.random.normal(size=)

如果我们知道数据位于正弦波上，而不知道振幅或周期，我们可以通过最小二乘曲线拟合找到它们。首先，我们必须定义要拟合的测试函数，这里有一个振幅和周期未知的正弦：

>>>

>>> def test_func(x, a, b):
...     return a * np.sin(b * x)

然后我们使用scipy.optimize.curve_fit()来寻找 a 和b：

>>>

>>> params, params_covariance = optimize.curve_fit(test_func, x_data, y_data, p0=[2, 2])
>>> print(params)
[  ]

。求标量函数的最小值

让我们定义以下函数：

>>>

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

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

这个函数的全局最小值大约为 - ，局部最小值大约为。

搜索最小值可以用 scipy.optimize.minimize() 来做，给定起点 x0，找到最小值的位置：

scipy.optimize.minimize()是一个复合对象，包含所有关于收敛的信息

>>>

>>> result = optimize.minimize(f, x0=0)
>>> result 
      fun: -...
 hess_inv: array([[...]])
      jac: array([-...e-])
  message: 'Optimization terminated successfully.'
     nfev: 
      nit: 5
     njev: 6
   status: 0
  success: True
        x: array([-...])
>>> result.x # The coordinate of the minimum  
array([-...])

方法：由于该函数是一个平滑函数，因此基于梯度下降的方法是不错的选择。一般lBFGS算法是一个不错的选择：

>>>

>>> optimize.minimize(f, x0=0, method="L-BFGS-B")  
      fun: array([-])
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([-1.42108547e-])
  message: ...'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 
      nit: 5
   status: 0
  success: True
        x: array([-])

上面只需要个函数就可以为最小值找到一个合适的值。

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

>>>

>>> res = optimize.minimize(f, x0=3, method="L-BFGS-B")
>>> res.x
array([])

如果我们不知道全局最小值的邻域来选择初始点，我们需要求助于更复杂的全局优化。为了找到全局最小值，我们使用
scipy.optimize.basinhopping() ，它将局部优化器与起点采样相结合：

>>>

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

约束：我们可以使用“bounds”参数将变量约束到区间 ：(0, )

边界列表

由于minimize()通常适用于x个多维，“bounds”参数是每个维度上的一个边界列表。

>>>

>>> res = optimize.minimize(f, x0=1,
...                         bounds=((0, ), ))
>>> res.x    
array([0.])

发生了什么？为什么我们会找到 0，这不是函数的最小值。

最小化多个变量的函数

要最小化多个变量，诀窍是将它们转换为多维变量(向量)的函数。

scipy.optimize.minimize_scalar() 是一种具有专用方法的函数，可以最小化只有一个变量的函数。

。寻找标量函数的根

要找到 f(x) = 0 的根，

我们可以使用scipy.optimize.root()：

>>>

>>> root = optimize.root(f, x0=1)  # our initial guess is 1
>>> root    # The full result
    fjac: array([[-1.]])
     fun: array([0.])
 message: 'The solution converged.'
    nfev: 
     qtf: array([1.33310463e-])
       r: array([-.])
  status: 1
 success: True
       x: array([0.])
>>> root.x  # Only the root found
array([0.])

请注意，仅找到一个根。图像显示在 - 附近有第二个根。我们通过调整我们的初始猜测来找到它的确切值：

>>>

>>> root2 = optimize.root(f, x0=-)
>>> root2.x
array([-])

scipy. optimization .root()也提供了各种算法，通过“method”参数进行设置。

现在我们已经找到了最小值和根并对其进行了曲线拟合，我们将所有这些结果放在一个图中。

统计和随机数：scipy.stats

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

分布：直方图和概率密度函数

对给定随机过程的观察，它们的直方图是随机过程 PDF（概率密度函数）的估计量：

>>>

>>> samples = np.random.normal(size=)
>>> bins = np.arange(-4, 5)
>>> bins
array([-4, -3, -2, -1,  0,  1,  2,  3,  4])
>>> histogram = np.histogram(samples, bins=bins, density=True)[0]
>>> bins = *(bins[1:] + bins[:-1])
>>> bins
array([-, -, -, -,  ,  ,  ,  ])
>>> from scipy import stats
>>> pdf = stats.norm.pdf(bins)  # norm is a distribution object

>>> plt.plot(bins, histogram) 
[]
>>> plt.plot(bins, pdf) 
[]

分布对象

scipy.stats.norm是一个分布对象：每个分布都表示为一个对象。norm 是正态分布，包含有 PDF、CDF 等等。

如果我们知道随机过程属于给定的随机过程族，例如正态过程，我们可以对观测值进行最大似然拟合以估计基础分布的参数。在这里，我们将正态过程拟合到观察到的数据：

>>>

loc, std = stats.norm.fit(samples)
print(loc,std)
-0. 

平均值、中位数和百分位数

均值

>>> np.mean(samples)     
-0.

中位数

>>> np.median(samples)     
 

与平均值不同，中位数对分布的尾部不敏感。

中位数也是百分位数 ，因为 % 的观察值低于它：

>>>

>>> stats.scoreatpercentile(samples, )     
 

同样，我们可以计算百分位数 ：

>>> stats.scoreatpercentile(samples, )     

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

统计检验

统计检验是一个决策指标。例如，如果我们有两组观察值，我们假设它们是从高斯过程生成的，我们可以使用 T 检验来确定两组观察值的均值是否显着不同：

>>>

a = np.random.normal(0, 1, size=)
b = np.random.normal(1, 1, size=)
stats.ttest_ind(a, b)
Out[]: Ttest_indResult(statistic=-, pvalue=0.)

输出结果由以下部分组成:

• T统计值: 它是一个数字，其符号与两个随机过程的差值成正比，其大小与这个差值的显著性有关。
• p值:两个过程相同的概率。如果它接近于1，那么这两个过程几乎肯定是相同的。它越接近于零，过程就越有可能有不同的方法。

数值积分：scipy.integrate

函数积分

最通用的是scipy.integrate.quad(). 计算

>>>

>>> from scipy.integrate import quad
>>> res, err = quad(np.sin, 0, np.pi/2)
>>> np.allclose(res, 1)   # Res是结果，应该接近1
True
>>> np.allclose(err, 1 - res)  # Err是误差
True

其它：
scipy.integrate.fixed_quad(),
scipy.integrate.quadrature(), scipy.integrate.romberg()...

积分微分方程

scipy.integrate还具有用于常微分方程 (ODE) 的例程。特别是scipy.integrate.odeint()求解以下形式的 ODE：

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

作为介绍，让我们计算在t=之间的常微分方程 dy/dt = -2y，其初始条件y(t=0) = 1。首先需要定义计算位置导数的函数:

>>>

>>> def calc_derivative(ypos, time):
...     return -2 * ypos

然后，计算y为时间的函数：

>>>

>>> from scipy.integrate import odeint
>>> time_vec = np.linspace(0, 4, )
>>> y = odeint(calc_derivative, y0=1, t=time_vec)

让我们用模块去计算一个更复杂的 ODE：阻尼弹簧谐振子。连接到弹簧上的质点的位置服从二阶ODE

k 为 弹簧常数，m 质量和

c 为阻尼系数。我们设置：

>>>

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

因此：

>>>

>>> eps = cviscous / (2 * mass * np.sqrt(kspring/mass))
>>> omega = np.sqrt(kspring / mass)

系统是欠阻尼的，因为：

>>> eps < 1
True

对于odeint()，二阶方程需要转换两个一阶方程组

函数计算速度和加速度：

>>>

>>> def calc_deri(yvec, time, eps, omega):
...     return (yvec[1], - * eps * omega * yvec[1] - omega **2 * yvec[0])

结果如下：

>>>

>>> time_vec = np.linspace(0, , )
>>> yinit = (1, 0)
>>> yarr = odeint(calc_deri, yinit, time_vec, args=(eps, omega))

odeint()使用LSODA(常微分方程的Livermore求解器，具有刚性和非刚性问题的自动方法切换)，请参阅ODEPACK Fortran库了解更多细节。

快速傅里叶变换：scipy.fftpack

scipy.fftpack模块计算快速傅里叶变换 (FFT) 并提供处理它们的实用程序。主要功能有：

• scipy.fftpack.fft() 计算 FFT
• scipy.fftpack.fftfreq() 生成采样频率
• scipy.fftpack.ifft() 计算从频率空间到信号空间的逆 FFT

例如，一个（噪声）输入信号 ( sig) 及其 FFT：

>>>

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

由于信号来自实函数，因此傅里叶变换是对称的。

峰值信号频率可以用 freqs[power.argmax()]

将高于该频率的傅里叶分量设置为零，并用 反 FFT scipy.fftpack.ifft()，得到一个滤波信号。

numpy.fft

Numpy 也有 FFT ( numpy.fft)的实现。但是，应该首选 scipy，因为它使用更有效的底层实现。

信号处理：scipy.signal

scipy.signal 用于典型的信号处理：一维、规则采样信号。

重采样 scipy.signal.resample()： 使用 FFT将信号重采样到n个点。

>>>

>>> t = np.linspace(0, 5, )
>>> x = np.sin(t)

>>> from scipy import signal
>>> x_resampled = signal.resample(x, )

>>> plt.plot(t, x) 
[]
>>> plt.plot(t[::4], x_resampled, 'ko') 
[]

请注意，在窗口的一侧，重新采样的准确性会降低，并且会产生涟漪效应。

这种重采样不同于提供的插值，scipy.interpolate因为它仅适用于定期采样的数据。

消解趋势 scipy.signal.detrend()：从信号中去除线性趋势：

>>>

>>> t = np.linspace(0, 5, )
>>> x = t + np.random.normal(size=)

>>> from scipy import signal
>>> x_detrended = signal.detrend(x)

>>> plt.plot(t, x) 
[]
>>> plt.plot(t, x_detrended) 
[]

过滤：对于非线性过滤，scipy.signal有过滤（中值过滤scipy.signal.medfilt()，Wiener scipy.signal.wiener()）。

scipy.signal 还有一整套用于设计线性滤波器（有限和无限响应滤波器）的工具。

频谱分析： scipy.signal.spectrogram()计算频谱图——连续时间窗口上的频谱——同时计算 scipy.signal.welch()功率谱密度（PSD）。

图像处理：scipy.ndimage

scipy.ndimage 提供将 n 维数组作为图像进行操作。

。图像的几何变换

改变方向，分辨率，..

>>>

>>> from scipy import misc  # Load an image
>>> face = misc.face(gray=True)

>>> from scipy import ndimage # Shift, roate and zoom it
>>> shifted_face = ndimage.shift(face, (, ))
>>> shifted_face2 = ndimage.shift(face, (, ), mode='nearest')
>>> rotated_face = ndimage.rotate(face, )
>>> cropped_face = face[:-, :-]
>>> zoomed_face = ndimage.zoom(face, 2)
>>> zoomed_face.shape
(, )

>>>

>>> plt.subplot)    


>>> plt.imshow(shifted_face, cmap=plt.cm.gray)    


>>> plt.axis('off')
(-, , , -)

>>> # etc.

图像过滤

生成一张嘈杂的脸：

>>>

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> face = face[:, -:]  # crop out square on right
>>> import numpy as np
>>> noisy_face = np.copy(face).astype(np.float)
>>> noisy_face += face.std() *  * 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))

在其它过滤器scipy.ndimage.filters和scipy.signal 可应用于图像。

数学形态学

数学形态学源于集合论。它表征和转换几何结构。尤其是二进制（黑白）图像，可以使用该理论进行转换：要转换的集合是相邻非零值像素的集合。该理论还扩展到灰度图像。

数学形态学操作使用结构元素 来修改几何结构。

首先生成一个结构元素：

>>>

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

• Erosion scipy.ndimage.binary_erosion()

a = np.zeros((7, 7), dtype=np.int)
a[1:6, 2:5] = 1

a
Out[4]: 
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)
Out[7]: 
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]])

ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)
Out[8]: 
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 scipy.ndimage.binary_dilation()

a = np.zeros((5, 5))

a[2, 2] = 1

a
Out[]: 
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)
Out[]: 
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.]])

• Opening scipy.ndimage.binary_opening()

a = np.zeros((5, 5), dtype=np.int)

a[1:4, 1:4] = 1

a[4, 4] = 1

a
Out[]: 
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]])

ndimage.binary_opening(a, structure=np.ones((3, 3))).astype(np.int)
Out[]: 
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]])

ndimage.binary_opening(a).astype(np.int)
Out[]: 
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]])

• Closing: scipy.ndimage.binary_closing()

打开操作会移除小结构，而关闭操作会填充小孔。因此，此类操作可用于“清理”图像。

>>>

>>> a = np.zeros((, ))
>>> a[:-, :-] = 1
>>> a +=  * np.random.standard_normal(a.shape)
>>> mask = a>=
>>> opened_mask = ndimage.binary_opening(mask)
>>> closed_mask = ndimage.binary_closing(opened_mask)

对于灰度值图像，腐蚀（或膨胀）相当于用以感兴趣像素为中心的结构元素覆盖的像素中的最小（或最大）值替换像素。

>>>

>>> 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]])

连接组件和图像测量

我们首先生成一个漂亮的合成二进制图像。

>>>

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

scipy.ndimage.label() 为每个连接的组件分配不同的标签：

>>>

>>> labels, nb = ndimage.label(mask)
>>> nb
8

现在计算每个连接组件的测量值：

>>>

>>> areas = ndimage.sum(mask, labels, range(1, labels.max()+1))
>>> areas   # The number of pixels in each connected component
array([.,   .,  .,  .,  .,  .,  .,  .])
>>> maxima = ndimage.maximum(sig, labels, range(1, labels.max()+1))
>>> maxima  # The maximum signal in each connected component
array([ ,   ,   ,   , ,
        ,  ,   ])

提取第 4 个连通分量，并裁剪它周围的数组：

>>>

>>> ndimage.find_objects(labels==4) 
[(slice(30L, 48L, None), slice(30L, 48L, None))]
>>> sl = ndimage.find_objects(labels==4)
>>> from matplotlib import pyplot as plt
>>> plt.imshow(sig[sl[0]])