Python数据分析教程（用Python做科学计算工具篇）

Python数据分析教程（用Python做科学计算工具篇）scipy.constants矢量量化/Kmeansscipy 非常丰富，这里我们只介绍一些重点部分，帮助我们了解如何将scipy用于科学计算。scipy 由特定于任务的子模块组成：scipy.cluster

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 由特定于任务的子模块组成：

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

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

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

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

1.1.文件输入/输出：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

1.2.特殊函数：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()

1.3.线性代数运算：scipy.linalg

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

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

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

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

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

LinAlgError: singular matrix

• 得到的数组谱为：

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

• 原始矩阵可以通过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.

1.4.插值：scipy.interpolate

scipy.interpolate对于从实验数据拟合函数非常有用，因此可以对不存在测度的点进行评估。该模块基于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)

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

>>>

>>> interpolation_time = np.linspace(0 1 50) >>> 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系列，插值点必须保持在给定数据点的范围内。

1.5 优化和拟合：scipy.optimize

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

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

>>>

>>> from scipy import optimize 1.5.1。曲线拟合

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

>>>

>>> x_data = np.linspace(-5 5 num=50) >>> y_data = 2.9 * np.sin(1.5 * x_data) np.random.normal(size=50)

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

>>>

>>> 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) [3.05931973 1.45754553] 1.5.2。求标量函数的最小值

让我们定义以下函数：

>>>

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

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

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

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

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

>>>

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

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

>>>

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

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

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

>>>

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

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

>>>

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

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

边界列表

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

>>>

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

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

最小化多个变量的函数

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

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

1.5.3。寻找标量函数的根

要找到 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: 10 qtf: array([1.33310463e-32]) r: array([-10.]) status: 1 success: True x: array([0.]) >>> root.x # Only the root found array([0.])

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

>>>

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

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

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

1.6 统计和随机数：scipy.stats

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

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

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

>>>

>>> samples = np.random.normal(size=1000) >>> 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 = 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 >>> pdf = stats.norm.pdf(bins) # norm is a distribution object >>> plt.plot(bins histogram) [<matplotlib.lines.Line2D object at ...>] >>> plt.plot(bins pdf) [<matplotlib.lines.Line2D object at ...>]

分布对象

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

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

>>>

loc std = stats.norm.fit(samples) print(loc std) -0.004822316383745015 1.0050446891199836 1.6.2 平均值、中位数和百分位数

均值

>>> np.mean(samples) -0.004822316383745015

中位数

>>> np.median(samples) 0.03361659877914626

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

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

>>>

>>> stats.scoreatpercentile(samples 50) 0.03361659877914626

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

>>> stats.scoreatpercentile(samples 90) 1.2565894469998924

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

1.6.3 统计检验

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

>>>

a = np.random.normal(0 1 size=100) b = np.random.normal(1 1 size=10) stats.ttest_ind(a b) Out[14]: Ttest_indResult(statistic=-3.6062755503726986 pvalue=0.0004718656448315962)

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

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

1.7 数值积分：scipy.integrate1.7.1 函数积分

最通用的是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()...

1.7.2 积分微分方程

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

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

作为介绍，让我们计算在t=0-4之间的常微分方程 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 40) >>> y = odeint(calc_derivative y0=1 t=time_vec)

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

k 为 弹簧常数，m 质量和

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

>>>

>>> mass = 0.5 # kg >>> kspring = 4 # N/m >>> cviscous = 0.4 # 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] -2.0 * eps * omega * yvec[1] - omega **2 * yvec[0])

结果如下：

>>>

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

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

1.8 快速傅里叶变换：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，因为它使用更有效的底层实现。

1.9 信号处理：scipy.signal

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

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

>>>

>>> t = np.linspace(0 5 100) >>> x = np.sin(t) >>> from scipy import signal >>> x_resampled = signal.resample(x 25) >>> plt.plot(t x) [<matplotlib.lines.Line2D object at ...>] >>> plt.plot(t[::4] x_resampled 'ko') [<matplotlib.lines.Line2D object at ...>]

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

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

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

>>>

>>> t = np.linspace(0 5 100) >>> x = t np.random.normal(size=100) >>> from scipy import signal >>> x_detrended = signal.detrend(x) >>> plt.plot(t x) [<matplotlib.lines.Line2D object at ...>] >>> plt.plot(t x_detrended) [<matplotlib.lines.Line2D object at ...>]

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

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

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

1.10 图像处理：scipy.ndimage

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

1.10.1。图像的几何变换

改变方向，分辨率，..

>>>

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

>>>

>>> 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.10.2 图像过滤

生成一张嘈杂的脸：

>>>

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

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

1.10.3 数学形态学

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

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

首先生成一个结构元素：

>>>

>>> 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[11]: 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[12]: 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[16]: 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[17]: 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[18]: 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((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)

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

>>>

>>> 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.10.4 连接组件和图像测量

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

>>>

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

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([190. 45. 424. 278. 459. 190. 549. 424.]) >>> maxima = ndimage.maximum(sig labels range(1 labels.max() 1)) >>> maxima # The maximum signal in each connected component array([ 1.80238238 1.13527605 5.51954079 2.49611818 6.71673619 1.80238238 16.76547217 5.51954079])

提取第 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]]) <matplotlib.image.AxesImage object at ...>