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.integrate
1.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 ...>