SciPy 1.12 中文文档（四）

原文：docs.scipy.org/doc/scipy-1.12.0/index.html

多维图像处理（scipy.ndimage）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/ndimage.html

简介

图像处理和分析通常被视为对值的二维数组进行操作。然而，在一些领域中，必须分析更高维度的图像。医学成像和生物成像是这方面的典型例子。由于其固有的多维特性，numpy非常适合于这类应用。scipy.ndimage包提供了许多通用的图像处理和分析函数，设计用于处理任意维度的数组。目前这些包括：线性和非线性滤波函数，二值形态学，B 样条插值以及对象测量。 ## 所有函数共享的属性

所有函数都具有一些共同的特性。特别是，所有函数都允许使用output参数来指定输出数组。通过这个参数，你可以指定一个数组，在操作中会就地改变这个数组并存储结果。在这种情况下，不会返回结果。通常情况下，使用output参数更加高效，因为可以利用现有数组来存储结果。

返回的数组类型取决于操作的类型，但大多数情况下等于输入的类型。然而，如果使用output参数，则结果的类型将等于指定的输出参数的类型。如果没有给出输出参数，则仍然可以指定输出结果的类型。这可以通过简单地将所需的numpy类型对象分配给输出参数来实现。例如：

>>> from scipy.ndimage import correlate
>>> import numpy as np
>>> correlate(np.arange(10), [1, 2.5])
array([ 0,  2,  6,  9, 13, 16, 20, 23, 27, 30])
>>> correlate(np.arange(10), [1, 2.5], output=np.float64)
array([  0\. ,   2.5,   6\. ,   9.5,  13\. ,  16.5,  20\. ,  23.5,  27\. ,  30.5]) 
```  ## 过滤函数

本节描述的函数都执行某种类型的空间过滤输入数组的操作：输出中的元素是相应输入元素邻域值的某种函数。我们称这些元素邻域为滤波核，通常是矩形的形状，但也可以具有任意的足迹。下面描述的许多函数允许通过*footprint*参数定义核的足迹。例如，可以如下定义十字形的核：

```py
>>> footprint = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
>>> footprint
array([[0, 1, 0],
 [1, 1, 1],
 [0, 1, 0]]) 

通常，核的原点位于核形状的维度的中心，通过将核形状的维度除以 2 来计算。例如，长度为 3 的一维核的原点位于第二个元素。例如，与长度为 3 的由 1 组成的滤波器的相关性：

>>> from scipy.ndimage import correlate1d
>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> correlate1d(a, [1, 1, 1])
array([0, 0, 1, 1, 1, 0, 0]) 

有时，为了选择内核的不同起点更为方便。因此，大多数函数支持 origin 参数，该参数指定了滤波器相对于其中心的起点。例如：

>>> a = [0, 0, 0, 1, 0, 0, 0]
>>> correlate1d(a, [1, 1, 1], origin = -1)
array([0, 1, 1, 1, 0, 0, 0]) 

这会使结果向左移动。这个特性通常不经常需要，但可能会很有用，特别是对于具有偶数大小的滤波器。一个很好的例子是计算后向和前向差分：

>>> a = [0, 0, 1, 1, 1, 0, 0]
>>> correlate1d(a, [-1, 1])               # backward difference
array([ 0,  0,  1,  0,  0, -1,  0])
>>> correlate1d(a, [-1, 1], origin = -1)  # forward difference
array([ 0,  1,  0,  0, -1,  0,  0]) 

我们也可以按以下方式计算前向差分：

>>> correlate1d(a, [0, -1, 1])
array([ 0,  1,  0,  0, -1,  0,  0]) 

然而，相比于较大的内核，使用 origin 参数更为高效。对于多维内核，origin 可以是一个数字，此时假定各轴上的起点相等，或者是一个序列，分别给出每个轴上的起点。

由于输出元素是输入元素邻域内元素的函数，所以数组的边界需要适当处理，即在边界外提供值。这是通过假设根据特定边界条件扩展数组来完成的。在下述描述的函数中，可以使用 mode 参数来选择边界条件，该参数必须是一个字符串，表示边界条件的名称。当前支持以下边界条件：

模式	描述	示例
“nearest”	使用边界处的值	[1 2 3]->[1 1 2 3 3]
“wrap”	周期性地复制数组	[1 2 3]->[3 1 2 3 1]
“reflect”	在边界处反射数组	[1 2 3]->[1 1 2 3 3]
“mirror”	在边界处镜像数组	[1 2 3]->[2 1 2 3 2]
“constant”	使用常量值，默认为 0.0	[1 2 3]->[0 1 2 3 0]

为了保持与插值例程的一致性，还支持以下同义词：

模式	描述
“grid-constant”	等同于 “constant”*
“grid-mirror”	等同于 “reflect”
“grid-wrap”	等同于 “wrap”

• “grid-constant” 和 “constant” 在滤波操作中等效，但在插值函数中具有不同的行为。为了 API 的一致性，滤波函数接受任一名称。

“constant” 模式是特殊的，因为它需要额外的参数来指定应使用的常量值。

注意，mirror 和 reflect 模式的区别仅在于边界处的样本是否在反射时重复。对于 mirror 模式，对称点正好位于最后一个样本，因此该值不重复。这种模式也被称为整样本对称，因为对称点位于最后一个样本上。类似地，reflect 常被称为半样本对称，因为对称点位于数组边界的半样本之外。

注意

实现这种边界条件的最简单方法是将数据复制到一个更大的数组中，并根据边界条件在边界处扩展数据。对于大数组和大滤波器核，这将非常消耗内存，因此下面描述的函数使用一种不需要分配大临时缓冲区的不同方法。

相关和卷积

• 函数correlate1d沿指定轴计算 1-D 相关。数组沿指定轴的行与给定的权重相关。权重参数必须是一个 1-D 数字序列。

• 函数correlate实现输入数组与给定核的多维相关。

• 函数convolve1d沿指定轴计算 1-D 卷积。数组沿指定轴的行与给定的权重进行卷积。权重参数必须是一个 1-D 数字序列。

• 函数convolve实现输入数组与给定核的多维卷积。

  注意

  卷积本质上是在镜像核之后进行相关。因此，origin参数的行为与相关的情况不同：结果向相反方向移动。

平滑滤波器

• 函数gaussian_filter1d实现 1-D 高斯滤波器。高斯滤波器的标准差通过参数sigma传递。将order = 0 设置为与高斯核的卷积。顺序为 1、2 或 3 对应于与高斯的一、二或三阶导数的卷积。高阶导数未实现。

• gaussian_filter 函数实现了多维高斯滤波器。高斯滤波器沿着每个轴的标准差通过参数sigma作为数字序列传递。如果sigma不是一个序列而是一个单一数字，则滤波器的标准差在所有方向上都相等。滤波器的顺序可以分别为每个轴指定。顺序为 0 对应于与高斯核的卷积。顺序为 1、2 或 3 对应于与高斯的一阶、二阶或三阶导数的卷积。高阶导数未实现。order参数必须是一个数字，以指定所有轴的相同顺序，或者是一个数字序列，以指定每个轴的不同顺序。下面的示例显示了在具有不同sigma值的测试数据上应用的滤波器。order参数保持为 0。

  

  注意

  多维滤波器被实现为一系列 1-D 高斯滤波器。中间数组以与输出相同的数据类型存储。因此，对于精度较低的输出类型，结果可能不精确，因为中间结果可能以不足的精度存储。可以通过指定更精确的输出类型来防止这种情况。

• uniform_filter1d 函数沿着给定轴计算给定size的 1-D 均匀滤波器。

• uniform_filter 实现了多维均匀滤波器。均匀滤波器的大小由size参数作为整数序列给出，每个轴一个。如果size不是一个序列，而是一个单一数字，则假定所有轴上的大小相等。

  注意

  多维滤波器被实现为一系列 1-D 均匀滤波器。中间数组以与输出相同的数据类型存储。因此，对于精度较低的输出类型，结果可能不精确，因为中间结果可能以不足的精度存储。可以通过指定更精确的输出类型来防止这种情况。

基于顺序统计的滤波器

• minimum_filter1d 函数沿着给定轴计算给定size的 1-D 最小值滤波器。

• maximum_filter1d 函数沿着给定轴计算给定size的 1-D 最大值滤波器。

• minimum_filter 函数用于计算多维最小值滤波器。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• maximum_filter 函数用于计算多维最大值滤波器。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• rank_filter 函数用于计算多维等级滤波器。rank 可以小于零，例如rank = -1 表示最大元素。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• percentile_filter 函数用于计算多维百分位数滤波器。percentile 可以小于零，例如percentile = -20 等同于 percentile = 80。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

• median_filter 函数用于计算多维中值滤波器。必须提供矩形核的尺寸或核的足迹。如果提供了size参数，必须是一系列尺寸或单个数字，此时滤波器沿每个轴的尺寸被假定相等。如果提供了footprint，必须是一个定义核形状的数组，其非零元素定义了核的形状。

导数

滤波器可以通过几种方法构造。函数gaussian_filter1d，描述于平滑滤波器，可用于沿指定轴使用order参数计算导数。其他导数滤波器包括 Prewitt 和 Sobel 滤波器：

• 函数prewitt计算沿给定轴的导数。

• 函数sobel计算沿给定轴的导数。

Laplace 滤波器通过所有轴上的二阶导数之和计算。因此，可以使用不同的二阶导数函数构造不同的 Laplace 滤波器。因此，我们提供了一个通用函数，该函数接受函数参数以计算给定方向上的二阶导数。

• 函数generic_laplace使用通过derivative2传递的函数计算二阶导数来计算拉普拉斯滤波器。函数derivative2应具有以下签名

  derivative2(input, axis, output, mode, cval, *extra_arguments, **extra_keywords) 
  

  应计算沿尺寸axis的二阶导数。如果output不为None，则应将其用于输出并返回None，否则应返回结果。mode，cval具有通常的意义。

  extra_arguments和extra_keywords参数可用于传递每次调用derivative2时传递给其的额外参数元组和命名参数字典。

  例如

  >>> def d2(input, axis, output, mode, cval):
  ...     return correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)
  ...
  >>> a = np.zeros((5, 5))
  >>> a[2, 2] = 1
  >>> from scipy.ndimage import generic_laplace
  >>> generic_laplace(a, d2)
  array([[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]) 
  

  为了演示extra_arguments参数的使用，我们可以执行

  >>> def d2(input, axis, output, mode, cval, weights):
  ...     return correlate1d(input, weights, axis, output, mode, cval, 0,)
  ...
  >>> a = np.zeros((5, 5))
  >>> a[2, 2] = 1
  >>> generic_laplace(a, d2, extra_arguments = ([1, -2, 1],))
  array([[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]) 
  

  或者

  >>> generic_laplace(a, d2, extra_keywords = {'weights': [1, -2, 1]})
  array([[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]) 
  

通过提供适当的二阶导数函数，以下两个函数使用generic_laplace实现：

• 函数laplace使用离散差分计算二阶导数的拉普拉斯（即与[1, -2, 1]卷积）。

• 函数gaussian_laplace使用gaussian_filter计算二阶导数的拉普拉斯滤波器。通过参数sigma作为数字序列传递高斯滤波器沿每个轴的标准差。如果sigma不是序列而是单个数字，则滤波器的标准差在所有方向上是相等的。

梯度幅度定义为所有方向梯度平方和的平方根。与通用拉普拉斯函数类似，有一个generic_gradient_magnitude函数用于计算数组的梯度幅度。

• 函数generic_gradient_magnitude使用通过derivative计算一阶导数的函数来计算梯度幅度。函数derivative应具有以下签名

  derivative(input, axis, output, mode, cval, *extra_arguments, **extra_keywords) 
  

  它应计算维度axis上的导数。如果output不是None，则应用于输出并返回None，否则应返回结果。mode，cval具有通常的意义。

  extra_arguments和extra_keywords参数可用于传递额外参数的元组和传递给derivative的命名参数字典。

  例如，sobel函数符合所需的签名

  >>> a = np.zeros((5, 5))
  >>> a[2, 2] = 1
  >>> from scipy.ndimage import sobel, generic_gradient_magnitude
  >>> generic_gradient_magnitude(a, sobel)
  array([[ 0\.        ,  0\.        ,  0\.        ,  0\.        ,  0\.        ],
   [ 0\.        ,  1.41421356,  2\.        ,  1.41421356,  0\.        ],
   [ 0\.        ,  2\.        ,  0\.        ,  2\.        ,  0\.        ],
   [ 0\.        ,  1.41421356,  2\.        ,  1.41421356,  0\.        ],
   [ 0\.        ,  0\.        ,  0\.        ,  0\.        ,  0\.        ]]) 
  

  查看generic_laplace的文档以获取使用extra_arguments和extra_keywords参数的示例。

函数sobel和prewitt函数符合所需的签名，因此可以直接与generic_gradient_magnitude一起使用。

• 函数gaussian_gradient_magnitude使用gaussian_filter来计算梯度幅度。高斯滤波器沿每个轴的标准偏差通过参数sigma作为数字序列传递。如果sigma不是序列而是一个单一数字，则滤波器的标准偏差沿所有方向均相等。

通用滤波函数

要实现滤波函数，可以使用通用函数，它们接受一个实现滤波操作的可调用对象。这些通用函数处理输入和输出数组的迭代，以及边界条件的实现等详细信息。只需提供一个实现实际滤波工作的回调函数的可调用对象。回调函数也可以用 C 语言编写，并使用PyCapsule传递（有关更多信息，请参见在 C 中扩展 scipy.ndimage）。

• generic_filter1d函数实现了一个通用的一维滤波函数，其中实际的滤波操作必须作为 Python 函数（或其他可调用对象）提供。generic_filter1d函数迭代数组的行，并在每行调用function。传递给function的参数是numpy.float64类型的一维数组。第一个包含当前行的值。根据filter_size和origin参数，在开头和结尾进行扩展。第二个数组应该就地修改，以提供行的输出值。例如，考虑沿一个维度的相关性：

  >>> a = np.arange(12).reshape(3,4)
  >>> correlate1d(a, [1, 2, 3])
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

  使用generic_filter1d函数可以实现相同的操作，如下所示：

  >>> def fnc(iline, oline):
  ...     oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:]
  ...
  >>> from scipy.ndimage import generic_filter1d
  >>> generic_filter1d(a, fnc, 3)
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

  在这里，默认情况下，内核的起源被假定为长度为 3 的滤波器的中间。因此，在调用函数之前，每个输入行都在开始和结束时扩展了一个值。

  可以选择定义并传递额外参数给滤波器函数。extra_arguments和extra_keywords参数可用于传递一组额外参数的元组和/或传递给每次调用的命名参数的字典。例如，我们可以将我们的滤波器的参数作为一个参数传递

  >>> def fnc(iline, oline, a, b):
  ...     oline[...] = iline[:-2] + a * iline[1:-1] + b * iline[2:]
  ...
  >>> generic_filter1d(a, fnc, 3, extra_arguments = (2, 3))
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

  或

  >>> generic_filter1d(a, fnc, 3, extra_keywords = {'a':2, 'b':3})
  array([[ 3,  8, 14, 17],
   [27, 32, 38, 41],
   [51, 56, 62, 65]]) 
  

• generic_filter函数实现了一个通用的滤波器函数，其中实际的过滤操作必须作为 Python 函数（或其他可调用对象）提供。generic_filter函数迭代数组并在每个元素上调用function。function的参数是一个 1-D 数组，类型为numpy.float64，包含在滤波器足迹内的当前元素周围的值。函数应返回一个可以转换为双精度数的单个值。例如，考虑一个相关性：

  >>> a = np.arange(12).reshape(3,4)
  >>> correlate(a, [[1, 0], [0, 3]])
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

  可以使用generic_filter执行相同的操作，如下所示：

  >>> def fnc(buffer):
  ...     return (buffer * np.array([1, 3])).sum()
  ...
  >>> from scipy.ndimage import generic_filter
  >>> generic_filter(a, fnc, footprint = [[1, 0], [0, 1]])
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

  在这里，指定了一个仅包含两个元素的核足迹。因此，过滤器函数接收到一个长度等于两个的缓冲区，该缓冲区与适当的权重相乘并求和得到结果。

  当调用generic_filter函数时，必须提供矩形核的大小或核的足迹。如果提供了size参数，则必须是大小序列或单个数字，其中情况下假定过滤器沿每个轴的大小相等。如果提供了footprint参数，则必须是一个数组，通过其中的非零元素定义核的形状。

  可以选择定义并传递额外参数给滤波器函数。extra_arguments和extra_keywords参数可用于传递一组额外参数的元组和/或传递给每次调用的命名参数的字典。例如，我们可以将我们的滤波器的参数作为一个参数传递

  >>> def fnc(buffer, weights):
  ...     weights = np.asarray(weights)
  ...     return (buffer * weights).sum()
  ...
  >>> generic_filter(a, fnc, footprint = [[1, 0], [0, 1]], extra_arguments = ([1, 3],))
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

  或

  >>> generic_filter(a, fnc, footprint = [[1, 0], [0, 1]], extra_keywords= {'weights': [1, 3]})
  array([[ 0,  3,  7, 11],
   [12, 15, 19, 23],
   [28, 31, 35, 39]]) 
  

这些函数按照从最后一个轴开始的行或元素进行迭代，即最后一个索引变化最快。对于重要的情况，这种迭代顺序是有保证的，以便根据空间位置调整滤波器。以下是使用实现滤波器并在迭代时跟踪当前坐标的类的示例。它执行与上述generic_filter相同的滤波操作，但还打印当前坐标：

>>> a = np.arange(12).reshape(3,4)
>>>
>>> class fnc_class:
...     def __init__(self, shape):
...         # store the shape:
...         self.shape = shape
...         # initialize the coordinates:
...         self.coordinates = [0] * len(shape)
...
...     def filter(self, buffer):
...         result = (buffer * np.array([1, 3])).sum()
...         print(self.coordinates)
...         # calculate the next coordinates:
...         axes = list(range(len(self.shape)))
...         axes.reverse()
...         for jj in axes:
...             if self.coordinates[jj] < self.shape[jj] - 1:
...                 self.coordinates[jj] += 1
...                 break
...             else:
...                 self.coordinates[jj] = 0
...         return result
...
>>> fnc = fnc_class(shape = (3,4))
>>> generic_filter(a, fnc.filter, footprint = [[1, 0], [0, 1]])
[0, 0]
[0, 1]
[0, 2]
[0, 3]
[1, 0]
[1, 1]
[1, 2]
[1, 3]
[2, 0]
[2, 1]
[2, 2]
[2, 3]
array([[ 0,  3,  7, 11],
 [12, 15, 19, 23],
 [28, 31, 35, 39]]) 

对于generic_filter1d函数，相同的方法适用，只是这个函数不会迭代正在被过滤的轴。然后generic_filter1d的示例如下：

>>> a = np.arange(12).reshape(3,4)
>>>
>>> class fnc1d_class:
...     def __init__(self, shape, axis = -1):
...         # store the filter axis:
...         self.axis = axis
...         # store the shape:
...         self.shape = shape
...         # initialize the coordinates:
...         self.coordinates = [0] * len(shape)
...
...     def filter(self, iline, oline):
...         oline[...] = iline[:-2] + 2 * iline[1:-1] + 3 * iline[2:]
...         print(self.coordinates)
...         # calculate the next coordinates:
...         axes = list(range(len(self.shape)))
...         # skip the filter axis:
...         del axes[self.axis]
...         axes.reverse()
...         for jj in axes:
...             if self.coordinates[jj] < self.shape[jj] - 1:
...                 self.coordinates[jj] += 1
...                 break
...             else:
...                 self.coordinates[jj] = 0
...
>>> fnc = fnc1d_class(shape = (3,4))
>>> generic_filter1d(a, fnc.filter, 3)
[0, 0]
[1, 0]
[2, 0]
array([[ 3,  8, 14, 17],
 [27, 32, 38, 41],
 [51, 56, 62, 65]]) 

傅里叶域滤波器

本节描述的函数在傅里叶域执行滤波操作。因此，此类函数的输入数组应与逆傅里叶变换函数兼容，例如numpy.fft模块中的函数。因此，我们需要处理可能是实部或复数傅里叶变换结果的数组。在实傅里叶变换的情况下，仅存储对称复数变换的一半。此外，需要知道在实 fft 变换之前被转换的轴的长度。本节描述的函数提供了一个n参数，在实变换的情况下，必须等于变换之前的实变换轴的长度。如果此参数小于零，则假定输入数组是复数傅里叶变换的结果。axis参数可用于指示执行实变换的轴。

• fourier_shift函数将输入数组乘以给定移位操作的多维傅里叶变换。shift参数是每个维度的移位序列或所有维度的单个值。

• fourier_gaussian函数将输入数组乘以具有给定标准差sigma的高斯滤波器的多维傅里叶变换。sigma参数是每个维度的值序列或所有维度的单个值。

• fourier_uniform函数将输入数组乘以具有给定大小size的均匀滤波器的多维傅里叶变换。size参数是每个维度的值序列或所有维度的单个值。

• fourier_ellipsoid 函数将输入数组与给定尺寸size的椭球形滤波器的多维傅里叶变换相乘。size参数是每个维度的值序列或所有维度的单个值。此函数仅实现维度为 1、2 和 3 的情况。## 插值函数

本节描述基于 B 样条理论的各种插值函数。关于 B 样条的良好介绍可参见[1]，图像插值的详细算法见[5]。

样条预滤波器

使用大于 1 阶的样条进行插值需要预滤波步骤。在 section [Interpolation functions]中描述的插值函数通过调用spline_filter应用预滤波，但可以通过将prefilter关键字设置为 False 来指示不执行此操作。如果在同一数组上执行多次插值操作，则这很有用。在这种情况下，只需进行一次预滤波并使用预滤波后的数组作为插值函数的输入更有效。以下两个函数实现预滤波：

• spline_filter1d 函数沿给定轴计算 1-D 样条滤波器。可选择提供输出数组。样条的阶数必须大于 1 且小于 6。

• spline_filter 函数计算多维样条滤波器。

  注意

  多维滤波器是通过一系列 1-D 样条滤波器实现的。中间数组与输出相同数据类型存储。因此，如果要求有限精度的输出，由于中间结果可能存储不足精度，结果可能不准确。可以通过指定高精度的输出类型来避免这种情况。

插值边界处理

所有的插值函数都使用样条插值来实现输入数组的某种几何变换。这要求将输出坐标映射到输入坐标，因此，可能需要处理超出边界的输入值。对于多维滤波函数中描述的相同问题，解决方法与过滤函数相同。因此，这些函数都支持一个 mode 参数，用于确定如何处理边界，并支持一个 cval 参数，在使用“constant”模式时给出一个常数值。下面展示了所有模式的行为，包括在非整数位置的行为。请注意，不同模式下的边界处理方式不同；reflect（又称grid-mirror）和grid-wrap涉及关于距离图像样本之间一半位置的对称或重复（虚线垂直线），而mirror和wrap将图像视为其范围恰好结束于第一个和最后一个样本点，而不是过了 0.5 个样本点。

图像样本的坐标落在每个轴上从 0 到shape[i] - 1的整数采样位置范围内，其中 i 是轴的索引。下图展示了在形状为(7, 7)的图像中位置为(3.7, 3.3)处的插值。对于阶数为 n 的插值，每个轴上涉及 n + 1 个样本。填充的圆圈显示了插值中涉及到的采样位置，在红色 x 的位置插值的图中。

插值函数

• geometric_transform 函数将任意几何变换应用于输入。给定的 mapping 函数在输出中的每个点被调用，以找到对应的输入坐标。mapping 必须是一个可调用对象，接受一个长度等于输出数组秩的元组，并返回与输入数组秩相等的对应输入坐标的元组。输出形状和输出类型可以选择性地提供。如果没有给出，它们等于输入的形状和类型。

  例如：

  >>> a = np.arange(12).reshape(4,3).astype(np.float64)
  >>> def shift_func(output_coordinates):
  ...     return (output_coordinates[0] - 0.5, output_coordinates[1] - 0.5)
  ...
  >>> from scipy.ndimage import geometric_transform
  >>> geometric_transform(a, shift_func)
  array([[ 0\.    ,  0\.    ,  0\.    ],
   [ 0\.    ,  1.3625,  2.7375],
   [ 0\.    ,  4.8125,  6.1875],
   [ 0\.    ,  8.2625,  9.6375]]) 
  

  可选地，可以定义并传递额外的参数给过滤器函数。extra_arguments 和 extra_keywords 参数可用于传递额外参数的元组和/或传递给每次调用中的导数的命名参数的字典。例如，我们可以像以下示例中一样传递偏移量。

  >>> def shift_func(output_coordinates, s0, s1):
  ...     return (output_coordinates[0] - s0, output_coordinates[1] - s1)
  ...
  >>> geometric_transform(a, shift_func, extra_arguments = (0.5, 0.5))
  array([[ 0\.    ,  0\.    ,  0\.    ],
   [ 0\.    ,  1.3625,  2.7375],
   [ 0\.    ,  4.8125,  6.1875],
   [ 0\.    ,  8.2625,  9.6375]]) 
  

  或

  >>> geometric_transform(a, shift_func, extra_keywords = {'s0': 0.5, 's1': 0.5})
  array([[ 0\.    ,  0\.    ,  0\.    ],
   [ 0\.    ,  1.3625,  2.7375],
   [ 0\.    ,  4.8125,  6.1875],
   [ 0\.    ,  8.2625,  9.6375]]) 
  

  注意

  映射函数也可以用 C 语言编写，并使用scipy.LowLevelCallable进行传递。更多信息请参阅在 C 中扩展 scipy.ndimage。

• 函数map_coordinates使用给定的坐标数组应用任意坐标变换。输出的形状从坐标数组的形状派生，通过删除第一个轴。参数coordinates用于找到输出中每个点在输入中对应的坐标。coordinates沿第一个轴的值是输入数组中找到输出值的坐标。（另见 numarray coordinates函数。）由于坐标可以是非整数坐标，因此在这些坐标处的输入值由请求的阶数的样条插值确定。

  这里是一个示例，插值一个 2D 数组在(0.5, 0.5)和(1, 2)处：

  >>> a = np.arange(12).reshape(4,3).astype(np.float64)
  >>> a
  array([[  0.,   1.,   2.],
   [  3.,   4.,   5.],
   [  6.,   7.,   8.],
   [  9.,  10.,  11.]])
  >>> from scipy.ndimage import map_coordinates
  >>> map_coordinates(a, [[0.5, 2], [0.5, 1]])
  array([ 1.3625,  7.]) 
  

• 函数affine_transform将仿射变换应用于输入数组。给定的变换矩阵和偏移量用于查找输出中每个点在输入中对应的坐标。在计算的坐标处的输入值由请求的阶数的样条插值确定。变换矩阵必须是 2-D，或者也可以作为 1-D 序列或数组给出。在后一种情况下，假定矩阵是对角的。然后应用更有效的插值算法，利用问题的可分离性。输出形状和输出类型可以选择性地提供。如果未给出，则它们与输入的形状和类型相等。

• 函数shift返回输入的移位版本，使用请求的order的样条插值。

• 函数zoom返回输入的重新缩放版本，使用请求的order的样条插值。

• 函数rotate返回在由参数axes给出的两个轴定义的平面中旋转的输入数组，使用请求的order的样条插值。角度必须以度为单位给出。如果reshape为真，则输出数组的大小将适应包含旋转输入的内容。## 形态学

二进制形态学

• generate_binary_structure 函数生成一个用于二值形态学操作的二值结构元素。必须提供结构的rank。返回的结构大小在每个方向上都等于三。每个元素的值等于 1，如果从元素到中心的欧氏距离的平方小于或等于connectivity。例如，生成 2D 4 连接和 8 连接结构如下：

  >>> from scipy.ndimage import generate_binary_structure
  >>> generate_binary_structure(2, 1)
  array([[False,  True, False],
   [ True,  True,  True],
   [False,  True, False]], dtype=bool)
  >>> generate_binary_structure(2, 2)
  array([[ True,  True,  True],
   [ True,  True,  True],
   [ True,  True,  True]], dtype=bool) 
  

这是generate_binary_structure在 3D 中的视觉呈现：

大多数二值形态学函数可以用基本操作腐蚀和膨胀来表示，这些操作可以在这里看到：

• binary_erosion 函数使用给定的结构元素对任意秩的数组进行二值腐蚀。原点参数控制结构元素的放置，如过滤器函数中所述。如果没有提供结构元素，则使用generate_binary_structure生成一个连接度等于 1 的元素。border_value参数给出边界外数组的值。腐蚀重复iterations次。如果iterations小于 1，则重复腐蚀直到结果不再改变。如果给出mask数组，则仅在每次迭代时修改对应掩码元素处具有真值的元素。

• binary_dilation 函数实现了带有给定结构元素的任意秩数组的二值膨胀。原点参数控制结构元素的放置位置，如滤波函数中所述。如果没有提供结构元素，则使用 generate_binary_structure 生成连接等于一的元素。border_value 参数指定边界外数组的值。膨胀重复 iterations 次。如果 iterations 小于一，则膨胀重复直到结果不再改变。如果给定 mask 数组，则只有在每次迭代时对应掩模元素处为 true 的元素才会被修改。

这里是使用 binary_dilation 的示例，通过多次使用数据数组作为掩模，从边界膨胀空数组来找到所有接触边界的元素：

>>> struct = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
>>> a = np.array([[1,0,0,0,0], [1,1,0,1,0], [0,0,1,1,0], [0,0,0,0,0]])
>>> a
array([[1, 0, 0, 0, 0],
 [1, 1, 0, 1, 0],
 [0, 0, 1, 1, 0],
 [0, 0, 0, 0, 0]])
>>> from scipy.ndimage import binary_dilation
>>> binary_dilation(np.zeros(a.shape), struct, -1, a, border_value=1)
array([[ True, False, False, False, False],
 [ True,  True, False, False, False],
 [False, False, False, False, False],
 [False, False, False, False, False]], dtype=bool) 

binary_erosion 和 binary_dilation 函数都有一个 iterations 参数，允许腐蚀或膨胀重复多次。用给定结构 n 次腐蚀或膨胀等同于使用自身膨胀 n-1 次的结构进行腐蚀或膨胀。提供了一个函数来计算多次膨胀自身的结构：

• iterate_structure 函数通过将输入结构自身膨胀 迭代 - 1 次来返回一个结构。

  例如：

  >>> struct = generate_binary_structure(2, 1)
  >>> struct
  array([[False,  True, False],
   [ True,  True,  True],
   [False,  True, False]], dtype=bool)
  >>> from scipy.ndimage import iterate_structure
  >>> iterate_structure(struct, 2)
  array([[False, False,  True, False, False],
   [False,  True,  True,  True, False],
   [ True,  True,  True,  True,  True],
   [False,  True,  True,  True, False],
   [False, False,  True, False, False]], dtype=bool)
  
  If the origin of the original structure is equal to 0, then it is
  also equal to 0 for the iterated structure. If not, the origin
  must also be adapted if the equivalent of the *iterations*
  erosions or dilations must be achieved with the iterated
  structure. The adapted origin is simply obtained by multiplying
  with the number of iterations. For convenience, the
  :func:`iterate_structure` also returns the adapted origin if the
  *origin* parameter is not ``None``:
  
  .. code:: python
  
   >>> iterate_structure(struct, 2, -1)
   (array([[False, False,  True, False, False],
   [False,  True,  True,  True, False],
   [ True,  True,  True,  True,  True],
   [False,  True,  True,  True, False],
   [False, False,  True, False, False]], dtype=bool), [-2, -2]) 
  

其他形态学操作可以通过腐蚀和膨胀来定义。以下函数为方便起见提供了一些这些操作：

• binary_opening 函数实现了使用给定结构元素进行任意秩数组的二进制开运算。二进制开运算相当于使用相同结构元素进行二进制腐蚀，然后进行二进制膨胀。原点参数控制结构元素的放置，如滤波函数中描述的那样。如果没有提供结构元素，则使用generate_binary_structure生成具有连接性等于一的元素。iterations 参数指定进行的腐蚀次数，然后是相同次数的膨胀。

• binary_closing 函数实现了使用给定结构元素进行任意秩数组的二进制闭运算。二进制闭运算相当于使用相同结构元素进行二进制膨胀，然后进行二进制腐蚀。原点参数控制结构元素的放置，如滤波函数中描述的那样。如果没有提供结构元素，则使用generate_binary_structure生成具有连接性等于一的元素。iterations 参数指定进行的膨胀次数，然后是相同次数的腐蚀。

• binary_fill_holes 函数用于关闭二进制图像中物体的孔洞，其中结构定义了孔洞的连接性。原点参数控制结构元素的放置，如滤波函数中描述的那样。如果没有提供结构元素，则使用generate_binary_structure生成具有连接性等于一的元素。

• binary_hit_or_miss 函数实现了对具有给定结构元素的任意秩数组的二进制击中或未击中变换。击中或未击中变换通过对输入与第一个结构元素进行腐蚀，对输入的逻辑 not 与第二个结构元素进行腐蚀，然后对这两个腐蚀结果进行逻辑 and 来计算。原点参数控制结构元素的放置位置，如 Filter functions 中所述。如果 origin2 等于 None，则设为 origin1 参数的值。如果未提供第一个结构元素，则使用 generate_binary_structure 生成一个连接度等于一的结构元素。如果未提供 structure2，则设置为 structure1 的逻辑 not。### 灰度形态学

灰度形态学操作是在具有任意值的数组上操作的二进制形态学操作的等价物。下面我们描述了侵蚀、膨胀、开运算和闭运算的灰度等价物。这些操作的实现方式与 Filter functions 中描述的滤波器类似，并且我们参考此部分来描述滤波器核和足迹的处理以及数组边界的处理。灰度形态学操作可以选择接受一个 structure 参数，该参数给出结构元素的值。如果未提供此参数，则假定结构元素是一个值为零的平坦结构元素。结构的形状可以通过 footprint 参数进行可选定义。如果未提供此参数，则假定结构是矩形的，大小等于 structure 数组的尺寸，或者由 size 参数定义（如果未提供 structure）。仅当 structure 和 footprint 都未给出时，size 参数才会被使用，此时结构元素被假定为矩形和平坦的，其尺寸由 size 给出。如果提供了 size 参数，则必须是一个尺寸序列或一个单独的数字，此时滤波器在每个轴上的尺寸是相同的。如果提供了 footprint 参数，则必须是一个数组，通过其非零元素定义核的形状。

类似于二进制腐蚀和膨胀，灰度腐蚀和膨胀也有相应的操作：

• grey_erosion 函数计算多维灰度侵蚀。

• grey_dilation 函数计算多维灰度膨胀。

灰度开运算和闭运算操作可以类似于它们的二值化对应物定义：

• grey_opening 函数实现了任意秩数组的灰度开运算。灰度开运算等效于灰度腐蚀后再进行灰度膨胀。

• grey_closing 函数实现了任意秩数组的灰度闭运算。灰度闭运算等效于灰度膨胀后再进行灰度腐蚀。

• morphological_gradient 函数实现了任意秩数组的灰度形态梯度。灰度形态梯度等于灰度膨胀与灰度腐蚀的差。

• morphological_laplace 函数实现了任意秩数组的灰度形态拉普拉斯。灰度形态拉普拉斯等于灰度膨胀与灰度腐蚀的和减去两倍的输入。

• white_tophat 函数实现了任意秩数组的白顶帽滤波器。白顶帽等于输入与灰度开运算的差。

• black_tophat 函数实现了任意秩数组的黑顶帽滤波器。黑顶帽等于灰度闭运算与输入的差。

距离变换用于计算对象中每个元素到背景的最小距离。以下函数实现了三种不同距离度量的距离变换：欧几里得距离、城市街区距离和棋盘距离。

• 函数distance_transform_cdt使用 chamfer 类型算法计算输入的距离变换，通过将每个对象元素（定义为大于零的值）替换为到背景的最短距离（所有非对象元素）。结构确定所做 chamfering 的类型。如果结构等于“cityblock”，则使用generate_binary_structure生成一个距离平方等于 1 的结构。如果结构等于“chessboard”，则使用generate_binary_structure生成一个距离平方等于数组秩的结构。这些选择对应于二维中城市街区距离和棋盘距离度量的常见解释。

  除了距离变换外，还可以计算特征变换。在这种情况下，返回结果的第一个轴上最接近背景元素的索引。可以使用return_distances和return_indices标志来指示是否返回距离变换、特征变换或两者。

  可以使用distances和indices参数来提供必须是正确大小和类型（numpy.int32）的可选输出数组。用于实现此功能的算法的基本原理在[2]中描述。

• 函数distance_transform_edt通过将每个对象元素（定义为大于零的值）替换为到背景的最短欧几里得距离，计算输入的精确欧几里得距离变换。

  除了距离变换外，还可以计算特征变换。在这种情况下，返回结果的第一个轴上最接近背景元素的索引。可以使用return_distances和return_indices标志来指示是否返回距离变换、特征变换或两者。

  可选地，可以通过sampling参数给出每个轴向的采样，它应该是与输入秩相等的长度序列，或者是一个单一的数，在这种情况下，假定沿所有轴向的采样是相等的。

  可以使用distances和indices参数来提供必须是正确大小和类型（numpy.float64和numpy.int32）的可选输出数组。用于实现此功能的算法的基本原理在[3]中描述。

• 函数distance_transform_bf使用蛮力算法计算输入的距离变换，通过将每个对象元素（由大于零的值定义）替换为到背景（所有非对象元素）的最短距离。度量必须是“euclidean”、“cityblock”或“chessboard”之一。

  除了距离变换之外，还可以计算特征变换。在这种情况下，返回结果的第一个轴上最接近的背景元素的索引。return_distances和return_indices标志可以用于指示是否必须返回距离变换、特征变换或两者。

  可选地，沿每个轴的采样可以由sampling参数给出，该参数应该是与输入秩相等的长度序列，或者是一个单独的数字，其中假定沿所有轴的采样相等。这个参数仅在欧几里得距离变换的情况下使用。

  distances和indices参数可用于提供必须具有正确大小和类型(numpy.float64和numpy.int32)的可选输出数组。

  注意

  此函数使用缓慢的蛮力算法，函数distance_transform_cdt可用于更有效地计算城市块和棋盘距离变换。函数distance_transform_edt可用于更有效地计算精确的欧几里得距离变换。

分割和标记

分割是从背景中分离感兴趣的对象的过程。可能最简单的方法是强度阈值分割，可以轻松使用numpy函数完成：

>>> a = np.array([[1,2,2,1,1,0],
...               [0,2,3,1,2,0],
...               [1,1,1,3,3,2],
...               [1,1,1,1,2,1]])
>>> np.where(a > 1, 1, 0)
array([[0, 1, 1, 0, 0, 0],
 [0, 1, 1, 0, 1, 0],
 [0, 0, 0, 1, 1, 1],
 [0, 0, 0, 0, 1, 0]]) 

结果是一个二进制图像，其中仍需要识别和标记各个对象。函数label生成一个数组，其中每个对象被分配一个唯一的编号：

• 函数label生成一个数组，其中输入中的对象用整数索引标记。它返回一个元组，包括对象标签数组和找到的对象数，除非给出了output参数，在这种情况下仅返回对象数。对象的连接性由一个结构元素定义。例如，在 2D 中使用 4 连接结构元素：

  >>> a = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
  >>> s = [[0, 1, 0], [1,1,1], [0,1,0]]
  >>> from scipy.ndimage import label
  >>> label(a, s)
  (array([[0, 1, 1, 0, 0, 0],
   [0, 1, 1, 0, 2, 0],
   [0, 0, 0, 2, 2, 2],
   [0, 0, 0, 0, 2, 0]]), 2) 
  

  这两个对象之间没有连接，因为我们无法将结构元素放置在既与第一个对象又与第二个对象重叠的方式。然而，8 连接结构元素会导致只有一个对象：

  >>> a = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
  >>> s = [[1,1,1], [1,1,1], [1,1,1]]
  >>> label(a, s)[0]
  array([[0, 1, 1, 0, 0, 0],
   [0, 1, 1, 0, 1, 0],
   [0, 0, 0, 1, 1, 1],
   [0, 0, 0, 0, 1, 0]]) 
  

  如果未提供结构元素，则通过调用generate_binary_structure（参见二元形态学）生成一个结构元素，使用连接性为 1（在 2D 中是第一个示例的 4 连接结构）。输入可以是任何类型，任何不等于零的值被视为对象的一部分。如果需要在移除不需要的对象后重新标记对象索引数组，则这是有用的。只需再次将标签函数应用于索引数组。例如：

  >>> l, n = label([1, 0, 1, 0, 1])
  >>> l
  array([1, 0, 2, 0, 3])
  >>> l = np.where(l != 2, l, 0)
  >>> l
  array([1, 0, 0, 0, 3])
  >>> label(l)[0]
  array([1, 0, 0, 0, 2]) 
  

  注意

  label使用的结构元素被假定为对称的。

有大量其他分割方法，例如，可以通过导数滤波器估计边界来获得对象的边界。其中一种方法是分水岭分割。函数watershed_ift生成一个数组，其中每个对象被分配一个唯一的标签，从定位对象边界的数组生成，例如通过梯度幅度滤波器。它使用一个包含对象初始标记的数组：

• 函数watershed_ift应用从标记算法开始的分水岭算法，使用图像森林变换，如[4]中所述。

• 该函数的输入是应用变换的数组，以及一个由唯一标签指定对象的标记数组，其中任何非零值均为标记。例如：

  >>> input = np.array([[0, 0, 0, 0, 0, 0, 0],
  ...                   [0, 1, 1, 1, 1, 1, 0],
  ...                   [0, 1, 0, 0, 0, 1, 0],
  ...                   [0, 1, 0, 0, 0, 1, 0],
  ...                   [0, 1, 0, 0, 0, 1, 0],
  ...                   [0, 1, 1, 1, 1, 1, 0],
  ...                   [0, 0, 0, 0, 0, 0, 0]], np.uint8)
  >>> markers = np.array([[1, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 2, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0]], np.int8)
  >>> from scipy.ndimage import watershed_ift
  >>> watershed_ift(input, markers)
  array([[1, 1, 1, 1, 1, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 2, 2, 2, 2, 2, 1],
   [1, 2, 2, 2, 2, 2, 1],
   [1, 2, 2, 2, 2, 2, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 1, 1, 1, 1, 1]], dtype=int8) 
  

  这里，使用了两个标记来指定一个对象（marker = 2）和背景（marker = 1）。这些标记的处理顺序是任意的：将背景的标记移动到数组的右下角会产生不同的结果：

  >>> markers = np.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, 2, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 1]], np.int8)
  >>> watershed_ift(input, markers)
  array([[1, 1, 1, 1, 1, 1, 1],
   [1, 1, 1, 1, 1, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 2, 2, 2, 1, 1],
   [1, 1, 1, 1, 1, 1, 1],
   [1, 1, 1, 1, 1, 1, 1]], dtype=int8) 
  

  结果是对象（marker = 2）较小，因为第二个标记被先处理了。如果第一个标记本应指定背景对象，则可能不是期望的效果。因此，watershed_ift明确地将具有负值的标记作为背景标记并在常规标记之后处理它们。例如，将第一个标记替换为负标记会产生与第一个示例类似的结果：

  >>> markers = np.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, 2, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, 0],
  ...                     [0, 0, 0, 0, 0, 0, -1]], np.int8)
  >>> watershed_ift(input, markers)
  array([[-1, -1, -1, -1, -1, -1, -1],
   [-1, -1,  2,  2,  2, -1, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1, -1,  2,  2,  2, -1, -1],
   [-1, -1, -1, -1, -1, -1, -1]], dtype=int8) 
  

  对象的连通性由结构元素定义。如果未提供结构元素，则通过调用 generate_binary_structure 函数生成一个结构元素（参见二进制形态学），使用连通性为一（在 2D 中为 4 连通结构）。例如，使用 8 连通结构来处理上一个示例将得到一个不同的对象：

  >>> watershed_ift(input, markers,
  ...               structure = [[1,1,1], [1,1,1], [1,1,1]])
  array([[-1, -1, -1, -1, -1, -1, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1,  2,  2,  2,  2,  2, -1],
   [-1, -1, -1, -1, -1, -1, -1]], dtype=int8) 
  

  注意

  watershed_ift 的实现将输入数据类型限制为 numpy.uint8 和 numpy.uint16。

对象测量

给定一个带标签对象的数组，可以测量各个对象的属性。find_objects 函数可用于生成一个切片列表，每个对象对应一个完全包含该对象的最小子数组：

• find_objects 函数查找标记数组中的所有对象，并返回一个切片列表，对应于包含对象的数组中最小的区域。

  例如：

  >>> a = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
  >>> l, n = label(a)
  >>> from scipy.ndimage import find_objects
  >>> f = find_objects(l)
  >>> a[f[0]]
  array([[1, 1],
   [1, 1]])
  >>> a[f[1]]
  array([[0, 1, 0],
   [1, 1, 1],
   [0, 1, 0]]) 
  

  函数 find_objects 返回所有对象的切片，除非 max_label 参数大于零，此时仅返回前 max_label 个对象。如果 label 数组中缺少索引，则返回 None 而不是切片。例如：

  >>> from scipy.ndimage import find_objects
  >>> find_objects([1, 0, 3, 4], max_label = 3)
  [(slice(0, 1, None),), None, (slice(2, 3, None),)] 
  

由 find_objects 生成的切片列表有助于找到数组中对象的位置和尺寸，也可用于对各个对象进行测量。例如，我们想要找到图像中一个对象强度的总和：

>>> image = np.arange(4 * 6).reshape(4, 6)
>>> mask = np.array([[0,1,1,0,0,0],[0,1,1,0,1,0],[0,0,0,1,1,1],[0,0,0,0,1,0]])
>>> labels = label(mask)[0]
>>> slices = find_objects(labels) 

然后我们可以计算第二个对象中元素的总和：

>>> np.where(labels[slices[1]] == 2, image[slices[1]], 0).sum()
80 

然而，这并不特别高效，对于其他类型的测量可能也更加复杂。因此，定义了一些测量函数，这些函数接受对象标签数组和要测量的对象的索引。例如，可以通过以下方式计算强度的总和：

>>> from scipy.ndimage import sum as ndi_sum
>>> ndi_sum(image, labels, 2)
80 

对于大数组和小对象，调用数组切片后再调用测量函数效率更高：

>>> ndi_sum(image[slices[1]], labels[slices[1]], 2)
80 

或者，我们可以通过单个函数调用对多个标签进行测量，返回一个结果列表。例如，在我们的示例中，要测量背景和第二个对象的值之和，我们提供一个标签列表：

>>> ndi_sum(image, labels, [0, 2])
array([178.0, 80.0]) 

下面描述的测量函数都支持index参数，以指示应测量哪些对象。index的默认值为None。这表示所有标签大于零的元素都应视为单个对象并进行测量。因此，在这种情况下，labels数组被视为由大于零的元素定义的掩码。如果index是一个数字或数字序列，则给出被测量对象的标签。如果index是一个序列，则返回结果的列表。如果返回多个结果的函数在index为单个数字时返回其结果作为元组，或在index为序列时返回其结果作为列表的元组。

• sum 函数计算由index给定标签的对象元素的总和，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• mean 函数计算由index给定标签的对象元素的平均值，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• variance 函数计算由index给定标签的对象元素的方差，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• standard_deviation 函数计算由index给定标签的对象元素的标准差，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• minimum 函数计算由index给定标签的对象元素的最小值，使用labels数组作为对象标签。如果index为None，则所有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• maximum 函数计算由index给定标签的对象元素的最大值，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• minimum_position 函数计算由index给定标签的对象元素的最小位置，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• maximum_position 函数计算由index给定标签的对象元素的最大位置，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• extrema 函数计算由index给定标签的对象元素的最小值、最大值及其位置，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。结果是一个元组，包含最小值、最大值、最小位置和最大位置。其结果与上述minimum、maximum、minimum_position 和 maximum_position 函数的结果组成的元组相同。

• center_of_mass 函数计算由index给定标签的对象的质心，使用labels数组作为对象标签。如果index为None，则所有具有非零标签值的元素被视为单个对象。如果label为None，则计算中使用input的所有元素。

• histogram 函数使用 index 指定的标签的对象计算直方图，并使用 labels 数组作为对象标签。如果 index 为 None，则所有具有非零标签值的元素将被视为单个对象。如果 label 为 None，则计算时使用 input 的所有元素。直方图由其最小值 (min)、最大值 (max) 和箱数 (bins) 定义。它们作为 numpy.int32 类型的 1-D 数组返回。## 在 C 中扩展 scipy.ndimage

scipy.ndimage 中的几个函数接受回调参数。这可以是一个 Python 函数或一个包含指向 C 函数指针的 scipy.LowLevelCallable。使用 C 函数通常更有效率，因为它避免了在数组的许多元素上调用 Python 函数的开销。要使用 C 函数，您必须编写一个 C 扩展，其中包含回调函数和一个返回 scipy.LowLevelCallable 的 Python 函数，该函数包含指向回调函数的指针。

支持回调的一个函数示例是 geometric_transform，它接受定义从所有输出坐标到输入数组中对应坐标的映射的回调函数。考虑以下 Python 示例，它使用 geometric_transform 实现了一个移位函数。

from scipy import ndimage

def transform(output_coordinates, shift):
    input_coordinates = output_coordinates[0] - shift, output_coordinates[1] - shift
    return input_coordinates

im = np.arange(12).reshape(4, 3).astype(np.float64)
shift = 0.5
print(ndimage.geometric_transform(im, transform, extra_arguments=(shift,))) 

我们还可以使用以下 C 代码来实现回调函数：

/* example.c */

#include  <Python.h>
#include  <numpy/npy_common.h>

static  int
_transform(npy_intp  *output_coordinates,  double  *input_coordinates,
  int  output_rank,  int  input_rank,  void  *user_data)
{
  npy_intp  i;
  double  shift  =  *(double  *)user_data;

  for  (i  =  0;  i  <  input_rank;  i++)  {
  input_coordinates[i]  =  output_coordinates[i]  -  shift;
  }
  return  1;
}

static  char  *transform_signature  =  "int (npy_intp *, double *, int, int, void *)";

static  PyObject  *
py_get_transform(PyObject  *obj,  PyObject  *args)
{
  if  (!PyArg_ParseTuple(args,  ""))  return  NULL;
  return  PyCapsule_New(_transform,  transform_signature,  NULL);
}

static  PyMethodDef  ExampleMethods[]  =  {
  {"get_transform",  (PyCFunction)py_get_transform,  METH_VARARGS,  ""},
  {NULL,  NULL,  0,  NULL}
};

/* Initialize the module */
static  struct  PyModuleDef  example  =  {
  PyModuleDef_HEAD_INIT,
  "example",
  NULL,
  -1,
  ExampleMethods,
  NULL,
  NULL,
  NULL,
  NULL
};

PyMODINIT_FUNC
PyInit_example(void)
{
  return  PyModule_Create(&example);
} 

关于编写 Python 扩展模块的更多信息可以在这里找到。如果 C 代码在文件 example.c 中，则可以在将其添加到 meson.build（查看 meson.build 文件中的示例）之后进行编译并遵循其中的步骤。完成后，运行脚本：

import ctypes
import numpy as np
from scipy import ndimage, LowLevelCallable

from example import get_transform

shift = 0.5

user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)
callback = LowLevelCallable(get_transform(), ptr)
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

产生与原始的 Python 脚本相同的结果。

在 C 版本中，_transform 是回调函数，参数 output_coordinates 和 input_coordinates 的作用与 Python 版本相同，而 output_rank 和 input_rank 则提供了 len(output_coordinates) 和 len(input_coordinates) 的等效值。变量 shift 通过 user_data 而不是 extra_arguments 传递。最后，C 回调函数返回一个整数状态，成功时为 1，否则为 0。

函数 py_transform 将回调函数封装在 PyCapsule 中。主要步骤如下：

• 初始化一个 PyCapsule。第一个参数是指向回调函数的指针。

• 第二个参数是函数签名，必须与 ndimage 预期的完全匹配。

• 上述示例中，我们使用了 scipy.LowLevelCallable 来指定由 ctypes 生成的 user_data。

  另一种方法是在胶囊上下文中提供数据，可以通过 PyCapsule_SetContext 设置，且在 scipy.LowLevelCallable 中省略指定 user_data。然而，在此方法中，我们需要处理数据的分配/释放 —— 在销毁胶囊后释放数据可以通过在 PyCapsule_New 的第三个参数中指定非空回调函数来完成。

用于 ndimage 的 C 回调函数都遵循此模式。下一节列出了接受 C 回调函数并给出函数原型的 ndimage 函数。

另见

支持低级回调参数的函数有：

generic_filter，generic_filter1d，geometric_transform

下面，我们展示了使用 Numba、Cython、ctypes 或 cffi 写代码的替代方式，而不是在 C 中编写包装器代码。

Numba

Numba 提供了在 Python 中轻松编写低级函数的方式。我们可以使用 Numba 来编写上述代码：

# example.py
import numpy as np
import ctypes
from scipy import ndimage, LowLevelCallable
from numba import cfunc, types, carray

@cfunc(types.intc(types.CPointer(types.intp),
                  types.CPointer(types.double),
                  types.intc,
                  types.intc,
                  types.voidptr))
def transform(output_coordinates_ptr, input_coordinates_ptr,
              output_rank, input_rank, user_data):
    input_coordinates = carray(input_coordinates_ptr, (input_rank,))
    output_coordinates = carray(output_coordinates_ptr, (output_rank,))
    shift = carray(user_data, (1,), types.double)[0]

    for i in range(input_rank):
        input_coordinates[i] = output_coordinates[i] - shift

    return 1

shift = 0.5

# Then call the function
user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)
callback = LowLevelCallable(transform.ctypes, ptr)

im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

Cython

从功能上看，与上述代码相同的代码可以在 Cython 中写成，而且更少繁文缛节，如下所示：

# example.pyx

from numpy cimport npy_intp as intp

cdef api int transform(intp *output_coordinates, double *input_coordinates,
                       int output_rank, int input_rank, void *user_data):
    cdef intp i
    cdef double shift = (<double *>user_data)[0]

    for i in range(input_rank):
        input_coordinates[i] = output_coordinates[i] - shift
    return 1 

# script.py

import ctypes
import numpy as np
from scipy import ndimage, LowLevelCallable

import example

shift = 0.5

user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)
callback = LowLevelCallable.from_cython(example, "transform", ptr)
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

cffi

使用 cffi，您可以与驻留在共享库（DLL）中的 C 函数进行接口。首先，我们需要编写共享库，在 Linux/OSX 上如下：

/*
 example.c
 Needs to be compiled with "gcc -std=c99 -shared -fPIC -o example.so example.c"
 or similar
 */

#include  <stdint.h>

int
_transform(intptr_t  *output_coordinates,  double  *input_coordinates,
  int  output_rank,  int  input_rank,  void  *user_data)
{
  int  i;
  double  shift  =  *(double  *)user_data;

  for  (i  =  0;  i  <  input_rank;  i++)  {
  input_coordinates[i]  =  output_coordinates[i]  -  shift;
  }
  return  1;
} 

调用库的 Python 代码如下：

import os
import numpy as np
from scipy import ndimage, LowLevelCallable
import cffi

# Construct the FFI object, and copypaste the function declaration
ffi = cffi.FFI()
ffi.cdef("""
int _transform(intptr_t *output_coordinates, double *input_coordinates,
 int output_rank, int input_rank, void *user_data);
""")

# Open library
lib = ffi.dlopen(os.path.abspath("example.so"))

# Do the function call
user_data = ffi.new('double *', 0.5)
callback = LowLevelCallable(lib._transform, user_data)
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

您可以在 cffi 文档中找到更多信息。

ctypes

使用 ctypes，与上述 cffi 的情况一样，C 代码和 so/DLL 的编译是一样的。Python 代码则有所不同：

# script.py

import os
import ctypes
import numpy as np
from scipy import ndimage, LowLevelCallable

lib = ctypes.CDLL(os.path.abspath('example.so'))

shift = 0.5

user_data = ctypes.c_double(shift)
ptr = ctypes.cast(ctypes.pointer(user_data), ctypes.c_void_p)

# Ctypes has no built-in intptr type, so override the signature
# instead of trying to get it via ctypes
callback = LowLevelCallable(lib._transform, ptr,
    "int _transform(intptr_t *, double *, int, int, void *)")

# Perform the call
im = np.arange(12).reshape(4, 3).astype(np.float64)
print(ndimage.geometric_transform(im, callback)) 

你可以在 ctypes 文档中找到更多信息。

参考资料

文件 IO（scipy.io）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/io.html

另请参见

NumPy IO routines

MATLAB 文件

loadmat(file_name[, mdict, appendmat])	Load MATLAB file.
savemat(file_name, mdict[, appendmat, ...])	保存字典的名称和数组到 MATLAB 风格的.mat 文件中。
whosmat(file_name[, appendmat])	列出 MATLAB 文件中的变量。

基本功能

我们将从导入scipy.io开始，并为方便起见称其为sio：

>>> import scipy.io as sio 

如果您正在使用 IPython，请尝试在sio上进行制表符完成。在众多选项中，您会找到：

sio.loadmat
sio.savemat
sio.whosmat 

这些是您在处理 MATLAB 文件时最可能使用的高级功能。您还会发现：

sio.matlab 

这是导入loadmat、savemat和whosmat的包。在sio.matlab中，您会找到mio模块。该模块包含loadmat和savemat使用的机制。偶尔您可能会发现自己重新使用此机制。

我该如何开始？

您可能有一个.mat文件，想要将其读入 SciPy。或者，您想要从 SciPy / NumPy 传递一些变量到 MATLAB。

为了避免使用 MATLAB 许可证，让我们从Octave开始。Octave 具有与 MATLAB 兼容的保存和加载功能。在命令行上启动 Octave（对我来说是octave）：

octave:1> a = 1:12
a =

   1   2   3   4   5   6   7   8   9  10  11  12

octave:2> a = reshape(a, [1 3 4])
a =

ans(:,:,1) =

   1   2   3

ans(:,:,2) =

   4   5   6

ans(:,:,3) =

   7   8   9

ans(:,:,4) =

   10   11   12

octave:3> save -6 octave_a.mat a % MATLAB 6 compatible
octave:4> ls octave_a.mat
octave_a.mat 

现在，到 Python：

>>> mat_contents = sio.loadmat('octave_a.mat')
>>> mat_contents
{'a': array([[[  1.,   4.,   7.,  10.],
 [  2.,   5.,   8.,  11.],
 [  3.,   6.,   9.,  12.]]]),
 '__version__': '1.0',
 '__header__': 'MATLAB 5.0 MAT-file, written by
 Octave 3.6.3, 2013-02-17 21:02:11 UTC',
 '__globals__': []}
>>> oct_a = mat_contents['a']
>>> oct_a
array([[[  1.,   4.,   7.,  10.],
 [  2.,   5.,   8.,  11.],
 [  3.,   6.,   9.,  12.]]])
>>> oct_a.shape
(1, 3, 4) 

现在让我们试着换个角度：

>>> import numpy as np
>>> vect = np.arange(10)
>>> vect.shape
(10,)
>>> sio.savemat('np_vector.mat', {'vect':vect}) 

然后回到 Octave：

octave:8> load np_vector.mat
octave:9> vect
vect =

  0  1  2  3  4  5  6  7  8  9

octave:10> size(vect)
ans =

    1   10 

如果要检查 MATLAB 文件的内容而不将数据读入内存，请使用whosmat命令：

>>> sio.whosmat('octave_a.mat')
[('a', (1, 3, 4), 'double')] 

whosmat返回一个元组列表，每个文件中的数组（或其他对象）都有一个。每个元组包含数组的名称、形状和数据类型。

MATLAB 结构

MATLAB 结构有点像 Python 字典，但字段名称必须是字符串。任何 MATLAB 对象都可以是字段的值。与 MATLAB 中的所有对象一样，结构实际上是结构数组，其中单个结构是形状为（1，1）的数组。

octave:11> my_struct = struct('field1', 1, 'field2', 2)
my_struct =
{
  field1 =  1
  field2 =  2
}

octave:12> save -6 octave_struct.mat my_struct 

我们可以在 Python 中加载它：

>>> mat_contents = sio.loadmat('octave_struct.mat')
>>> mat_contents
{'my_struct': array([[([[1.0]], [[2.0]])]],
 dtype=[('field1', 'O'), ('field2', 'O')]), '__version__': '1.0', '__header__': 'MATLAB 5.0 MAT-file, written by Octave 3.6.3, 2013-02-17 21:23:14 UTC', '__globals__': []}
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct.shape
(1, 1)
>>> val = oct_struct[0,0]
>>> val
([[1.0]], [[2.0]])
>>> val['field1']
array([[ 1.]])
>>> val['field2']
array([[ 2.]])
>>> val.dtype
dtype([('field1', 'O'), ('field2', 'O')]) 

在 SciPy 版本从 0.12.0 开始，MATLAB 结构返回为 NumPy 结构化数组，其字段命名为结构字段。您可以在上面的dtype输出中看到字段名称。还要注意：

>>> val = oct_struct[0,0] 

和：

octave:13> size(my_struct)
ans =

   1   1 

因此，在 MATLAB 中，结构数组必须至少是 2 维的，并且我们在读入 SciPy 时复制了这一点。如果您希望将所有长度为 1 的维度挤出，请尝试这样做：

>>> mat_contents = sio.loadmat('octave_struct.mat', squeeze_me=True)
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct.shape
() 

有时，将 MATLAB 结构加载为 Python 对象而不是 NumPy 结构化数组更方便 - 这可以使 Python 中的访问语法与 MATLAB 中的语法更加相似。为此，请使用struct_as_record=False参数设置为loadmat。

>>> mat_contents = sio.loadmat('octave_struct.mat', struct_as_record=False)
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct[0,0].field1
array([[ 1.]]) 

struct_as_record=False 与 squeeze_me 配合使用效果很好：

>>> mat_contents = sio.loadmat('octave_struct.mat', struct_as_record=False, squeeze_me=True)
>>> oct_struct = mat_contents['my_struct']
>>> oct_struct.shape # but no - it's a scalar
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'mat_struct' object has no attribute 'shape'
>>> type(oct_struct)
<class 'scipy.io.matlab.mio5_params.mat_struct'>
>>> oct_struct.field1
1.0 

可以以多种方式保存结构数组。一种简单的方法是使用字典：

>>> a_dict = {'field1': 0.5, 'field2': 'a string'}
>>> sio.savemat('saved_struct.mat', {'a_dict': a_dict}) 

被加载为：

octave:21> load saved_struct
octave:22> a_dict
a_dict =

  scalar structure containing the fields:

    field2 = a string
    field1 =  0.50000 

您还可以像这样将结构体再次保存回 MATLAB（或者在我们的情况下是 Octave）：

>>> dt = [('f1', 'f8'), ('f2', 'S10')]
>>> arr = np.zeros((2,), dtype=dt)
>>> arr
array([(0.0, ''), (0.0, '')],
 dtype=[('f1', '<f8'), ('f2', 'S10')])
>>> arr[0]['f1'] = 0.5
>>> arr[0]['f2'] = 'python'
>>> arr[1]['f1'] = 99
>>> arr[1]['f2'] = 'not perl'
>>> sio.savemat('np_struct_arr.mat', {'arr': arr}) 

MATLAB 单元数组

MATLAB 中的单元数组与 Python 列表相似，数组中的元素可以包含任何类型的 MATLAB 对象。事实上，它们最类似于 NumPy 对象数组，这就是我们如何将它们加载到 NumPy 中的方式。

octave:14> my_cells = {1, [2, 3]}
my_cells =
{
  [1,1] =  1
  [1,2] =

     2   3

}

octave:15> save -6 octave_cells.mat my_cells 

回到 Python：

>>> mat_contents = sio.loadmat('octave_cells.mat')
>>> oct_cells = mat_contents['my_cells']
>>> print(oct_cells.dtype)
object
>>> val = oct_cells[0,0]
>>> val
array([[ 1.]])
>>> print(val.dtype)
float64 

保存到 MATLAB 单元数组只需创建一个 NumPy 对象数组：

>>> obj_arr = np.zeros((2,), dtype=np.object)
>>> obj_arr[0] = 1
>>> obj_arr[1] = 'a string'
>>> obj_arr
array([1, 'a string'], dtype=object)
>>> sio.savemat('np_cells.mat', {'obj_arr':obj_arr}) 

octave:16> load np_cells.mat
octave:17> obj_arr
obj_arr =
{
  [1,1] = 1
  [2,1] = a string
} 

IDL 文件

readsav(文件名[, idict, python_dict, ...])	读取 IDL 的 .sav 文件。

Matrix Market 文件

mminfo(源)	从类似于 Matrix Market 文件的 '源' 返回大小和存储参数。
mmread(源)	从类似于 Matrix Market 的 '源' 中读取内容到矩阵中。
mmwrite(目标, a[, 注释, 字段, ...])	将稀疏或密集数组 a 写入类似于 Matrix Market 的 '目标' 文件。

Wav 声音文件（scipy.io.wavfile）

read(文件名[, mmap])	打开 WAV 文件。
write(文件名, rate, 数据)	将 NumPy 数组写入 WAV 文件。

Arff 文件（scipy.io.arff）

loadarff(f)	读取 arff 文件。

Netcdf

netcdf_file(文件名[, 模式, mmap, 版本, ...])	用于 NetCDF 数据的文件对象。

允许读取 NetCDF 文件（使用 pupynere 包的版本）

可执行教程

插值过渡指南

原文：docs.scipy.org/doc/scipy-1.12.0/notebooks/interp_transition_guide.html

本笔记本包含三组演示：

• 用于遗留 bug 兼容scipy.interpolate.interp2d的 FITPACK 低级替代品；

• 建议用于新代码中的scipy.interpolate.interp2d替代品；

• 展示了基于 2D FITPACK 的线性插值失败模式及推荐的替代方案。

注意： 由于本笔记本展示了interp2d的用法（标记为已弃用），我们将简单起见静默处理弃用警告：

import warnings
warnings.filterwarnings('ignore') 

1. 如何过渡到不再使用interp2d

interp2d在 2D 常规网格上和插值 2D 分散数据时会静默切换。这种切换基于（拉直后的）x、y和z数组的长度。简而言之，对于常规网格，请使用scipy.interpolate.RectBivariateSpline；对于分散插值，请使用bisprep/bisplev组合。下面我们提供了逐点转换的文字示例，这应该完全保留interp2d的结果。

1.1 interp2d在常规网格上的应用

我们从（稍作修改的）文档字符串示例开始。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d, RectBivariateSpline

x = np.arange(-5.01, 5.01, 0.25)
y = np.arange(-5.01, 7.51, 0.25)
xx, yy = np.meshgrid(x, y)
z = np.sin(xx**2 + 2.*yy**2)
f = interp2d(x, y, z, kind='cubic') 

这是“常规网格”代码路径的示例，因为

z.size == len(x) * len(y) 

True 

还要注意x.size != y.size：

x.size, y.size 

(41, 51) 

现在，让我们构建一个方便的函数来构造插值器并绘制它。

def plot(f, xnew, ynew):
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
    znew = f(xnew, ynew)

    ax1.plot(x, z[0, :], 'ro-', xnew, znew[0, :], 'b-')

    im = ax2.imshow(znew)
    plt.colorbar(im, ax=ax2)

    plt.show()
    return znew 

绘图：

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 7.51, 1e-2)
znew_i = plot(f, xnew, ynew) 

替代方案：使用RectBivariateSpline，结果完全相同

注意转置：首先在构造函数中，其次在评估结果时需要转置。这是为了撤消interp2d的转置操作。

r = RectBivariateSpline(x, y, z.T)

rt = lambda xnew, ynew: r(xnew, ynew).T
znew_r = plot(rt, xnew, ynew) 

from numpy.testing import assert_allclose
assert_allclose(znew_i, znew_r, atol=1e-14) 

1.2. 使用点的完整坐标进行interp2d（分散插值）

在这里，我们展示了前一个练习中的网格平铺以说明功能。

xxr = xx.ravel()
yyr = yy.ravel()
zzr = z.ravel()

f = interp2d(xxr, yyr, zzr, kind='cubic') 

注意这是“非常规网格”代码路径，用于分散数据，其中len(x) == len(y) == len(z)。

len(xxr) == len(yyr) == len(zzr) 

True 

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 7.51, 1e-2)
znew_i = plot(f, xnew, ynew) 

替换：直接使用scipy.interpolate.bisplrep / scipy.interpolate.bisplev

from scipy.interpolate import bisplrep, bisplev
tck = bisplrep(xxr, yyr, zzr, kx=3, ky=3, s=0)
# convenience: make up a callable from bisplev
ff = lambda xnew, ynew: bisplev(xnew, ynew, tck).T   # Note the transpose, to mimic what interp2d does

znew_b = plot(ff, xnew, ynew) 

assert_allclose(znew_i, znew_b, atol=1e-15) 

2. 替代interp2d：正则网格

对于新代码，推荐的替代方案是RegularGridInterpolator。这是一个独立的实现，不基于 FITPACK。支持最近邻、线性插值和奇次张量积样条。

样条结节保证与数据点重合。

注意，这里：

1. 元组参数是(x, y)

2. z数组需要转置

3. 关键字名称是method，而不是kind

4. bounds_error参数默认为True。

from scipy.interpolate import RegularGridInterpolator as RGI

r = RGI((x, y), z.T, method='linear', bounds_error=False) 

评估：创建一个 2D 网格。使用indexing='ij'和sparse=True以节省一些内存：

xxnew, yynew = np.meshgrid(xnew, ynew, indexing='ij', sparse=True) 

评估时，请注意元组参数：

znew_reggrid = r((xxnew, yynew)) 

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))

# Again, note the transpose to undo the `interp2d` convention
znew_reggrid_t = znew_reggrid.T

ax1.plot(x, z[0, :], 'ro-', xnew, znew_reggrid_t[0, :], 'b-')

im = ax2.imshow(znew_reggrid_t)
plt.colorbar(im, ax=ax2) 

<matplotlib.colorbar.Colorbar at 0x7fa55ccd5f10> 

3. 散点 2D 线性插值：优先使用LinearNDInterpolator而不是SmoothBivariateSpline或bisplrep

对于 2D 散点线性插值，SmoothBivariateSpline和biplrep可能会发出警告，或者无法插值数据，或者产生带有结节远离数据点的样条。“相反，建议使用LinearNDInterpolator，它基于通过QHull对数据进行三角剖分。

# TestSmoothBivariateSpline::test_integral
from scipy.interpolate import SmoothBivariateSpline, LinearNDInterpolator

x = np.array([1,1,1,2,2,2,4,4,4])
y = np.array([1,2,3,1,2,3,1,2,3])
z = np.array([0,7,8,3,4,7,1,3,4]) 

现在，使用基于 Qhull 的数据三角剖分进行线性插值：

xy = np.c_[x, y]   # or just list(zip(x, y))
lut2 = LinearNDInterpolator(xy, z)

X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y) 

结果易于理解和解释：

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

ax.plot_wireframe(X, Y, lut2(X, Y))
ax.scatter(x, y, z,  'o', color='k', s=48) 

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa55cbd8250> 

注意，bisplrep做了一些不同的事情！它可能会将样条结节放在数据之外。

作为说明，考虑前面示例中的相同数据：

tck = bisplrep(x, y, z, kx=1, ky=1, s=0)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

xx = np.linspace(min(x), max(x))
yy = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xx, yy)
Z = bisplev(xx, yy, tck)
Z = Z.reshape(*X.shape).T

ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2)
ax.scatter(x, y, z,  'o', color='k', s=48) 

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa55cc26310> 

此外，SmoothBivariateSpline无法插值数据。再次使用前面示例中的相同数据。

lut = SmoothBivariateSpline(x, y, z, kx=1, ky=1, s=0)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

xx = np.linspace(min(x), max(x))
yy = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xx, yy)

ax.plot_wireframe(X, Y, lut(xx, yy).T, rstride=4, cstride=4)
ax.scatter(x, y, z,  'o', color='k', s=48) 

<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fa55cc6ebe0> 

注意，SmoothBivariateSpline和bisplrep的结果都存在缺陷，不像LinearNDInterpolator那样。此处所示的问题是针对线性插值报告的，然而 FITPACK 的结节选择机制并不保证对高阶（如三次）样条曲面避免这些问题。

SciPy API

原文：docs.scipy.org/doc/scipy-1.12.0/reference/index.html

从 SciPy 导入

在 Python 中，库的公共 API 和私有实现细节之间的区分并不总是清晰的。与 Java 等其他语言不同，Python 中可以访问“私有”函数或对象。偶尔这可能很方便，但请注意，如果这样做，您的代码在未来版本中可能会无预警地中断。一些广泛认可的 Python 公共 API 规则包括：

• 方法/函数/类和模块属性名称以下划线开头的是私有的。

• 如果类名以下划线开头，则其所有成员都不是公共的，无论它们是否以下划线开头。

• 如果包中的模块名称以下划线开头，则其所有成员都不是公共的，无论它们是否以下划线开头。

• 如果模块或包定义了__all__，则这是官方定义的公共接口。

• 如果模块或包没有定义__all__，则所有不以下划线开头的名称都是公共的。

注意

阅读上述指南，可以得出结论，每个私有模块或对象都以下划线开头。但事实并非如此；下划线的存在确实标志着某些内容为私有，但缺少下划线并不意味着某些内容为公共的。

在 SciPy 中，有些模块的名称不以下划线开头，但应视为私有。为了澄清这些模块是哪些，我们在下面定义了 SciPy 的公共 API，并提供了一些关于如何从 SciPy 导入模块/函数/对象的建议。

从 SciPy 导入函数的指南

SciPy 子模块的命名空间中的所有内容都是公共的。通常在 Python 中建议使用命名空间。例如，函数 curve_fit（在 scipy/optimize/_minpack_py.py 中定义）应该这样导入：

import scipy
result = scipy.optimize.curve_fit(...) 

或者，您可以像这样使用子模块作为命名空间：

from scipy import optimize
result = optimize.curve_fit(...) 

注意

对于 scipy.io，推荐使用 import scipy，因为 io 也是 Python 标准库中的模块名称。

在某些情况下，公共 API 是更深层次的。例如，scipy.sparse.linalg 模块是公共的，它包含的函数在 scipy.sparse 命名空间中不可用。如果选择第二种形式，则代码更容易理解，例如，以下代码立即清楚地表明 lomax 是一个分布：

# first form
from scipy import stats
stats.lomax(...)

# second form
from scipy.stats import distributions
distributions.lomax(...) 

在这种情况下，如果文档中指明该子模块是公共的，则可以选择第二种形式。当然，您仍然可以使用：

import scipy
scipy.stats.lomax(...)
# or
scipy.stats.distributions.lomax(...) 

注意

SciPy 使用延迟加载机制，这意味着只有在首次尝试访问模块时才会将其加载到内存中。

注意

scipy 命名空间本身还包含从 numpy 导入的函数。这些函数仍然存在以保持向后兼容性，但应直接从 numpy 导入。

API 定义

下面列出的每个子模块都是公共的。这意味着这些子模块不太可能被重命名或以不兼容的方式进行更改，如果必须更改，将在 SciPy 的一个版本中引发弃用警告。

• scipy

• scipy.cluster

  • scipy.cluster.vq

  • scipy.cluster.hierarchy

• scipy.constants

• scipy.datasets

• scipy.fft

• scipy.fftpack

• scipy.integrate

• scipy.interpolate

• scipy.io

  • scipy.io.arff

  • scipy.io.matlab

  • scipy.io.wavfile

• scipy.linalg

  • scipy.linalg.blas

  • scipy.linalg.cython_blas

  • scipy.linalg.lapack

  • scipy.linalg.cython_lapack

  • scipy.linalg.interpolative

• scipy.misc

• scipy.ndimage

• scipy.odr

• scipy.optimize

  • scipy.optimize.cython_optimize

• scipy.signal

  • scipy.signal.windows

• scipy.sparse

  • scipy.sparse.linalg

  • scipy.sparse.csgraph

• scipy.spatial

  • scipy.spatial.distance

  • scipy.spatial.transform

• scipy.special

• scipy.stats

  • scipy.stats.contingency

  • scipy.stats.distributions

  • scipy.stats.mstats

  • scipy.stats.qmc

  • scipy.stats.sampling

SciPy 结构

所有 SciPy 模块应遵循以下约定。在此处，SciPy 模块 定义为位于 scipy/ 目录中的 Python 包，比如 yyy。

• 理想情况下，每个 SciPy 模块应尽可能自包含。即应最小化对其他包或模块的依赖。当然，假定对 NumPy 的依赖。

• 目录 yyy/ 包含：

  • 文件 meson.build 包含子模块的构建配置。

  • 目录 tests/ 包含文件 test_<name>.py，对应模块 yyy/<name>{.py,.so,/}。

• 私有模块应以下划线 _ 前缀，例如 yyy/_somemodule.py。

• 用户可见的函数应遵循 NumPy 文档风格。

• 模块的 __init__.py 应包含其主要参考文档，位于其 docstring 中。这与 Sphinx 文档在 doc/ 下的连接通过 Sphinx 的 automodule 指令。

  参考文档应首先使用 autosummary:: 指令列出模块内容的分类列表，随后解释了解模块使用的重要点。

  教程风格的文档应详细示例，需单独放置于 doc/source/tutorial/。

参见现有的 SciPy 子模块以获取指导。