多维图像处理 (scipy.ndimage)

简介

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

所有函数的共有属性

所有函数都具有一些共同属性。值得注意的是，所有函数都允许通过 output 参数指定输出数组。使用此参数，您可以指定一个数组，该数组将根据操作结果进行原地（in-place）修改。在这种情况下，函数不返回结果。通常，使用 output 参数效率更高，因为使用现有数组来存储结果。

返回数组的类型取决于操作类型，但在大多数情况下，它与输入的类型相同。但是，如果使用了 output 参数，则结果的类型与指定的输出参数类型相同。如果没有给出输出参数，仍然可以指定输出结果应有的类型。这只需通过将所需的 numpy 类型对象赋值给 output 参数即可实现。例如：

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

滤波函数

本节描述的函数均对输入数组执行某种类型的空间滤波：输出中的元素是对应输入元素邻域内值的某种函数。我们将这个元素邻域称为滤波核（filter kernel），它通常是矩形的，但也可能有任意形状。下面描述的许多函数允许您通过 footprint 参数传递一个掩码来定义核的形状。例如，交叉形状的核可以定义如下：

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

通常，核的原点位于中心，通过将核形状的尺寸除以二来计算。例如，长度为 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 参数选择边界条件，该参数必须是带有边界条件名称的字符串。目前支持以下边界条件：

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]

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

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

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

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

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

注意

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

相关和卷积

• correlate1d 函数沿给定轴计算一维互相关。数组沿给定轴的行与给定的 weights 进行互相关。weights 参数必须是一个一维数字序列。

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

• convolve1d 函数沿给定轴计算一维卷积。数组沿给定轴的行与给定的 weights 进行卷积。weights 参数必须是一个一维数字序列。

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

  注意

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

平滑滤波器

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

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

  

  注意

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

• uniform_filter1d 函数沿给定轴计算给定 size 的一维均匀滤波器（均值滤波）。

• uniform_filter 实现多维均匀滤波器。均匀滤波器的尺寸由 size 参数作为整数序列提供给每个轴。如果 size 不是序列而是一个数字，则假定沿所有轴的尺寸相等。

  注意

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

基于顺序统计量的滤波器

• minimum_filter1d 函数沿给定轴计算给定 size 的一维最小值滤波器。

• maximum_filter1d 函数沿给定轴计算给定 size 的一维最大值滤波器。

• minimum_filter 函数计算多维最小值滤波器。必须提供矩形核的尺寸或核的脚印（footprint）。如果提供 size 参数，它必须是一个尺寸序列或一个数字（在这种情况下，假定滤波器在每个轴上的尺寸相等）。如果提供 footprint，它必须是一个通过其非零元素定义核形状的数组。

• maximum_filter 函数计算多维最大值滤波器。必须提供矩形核的尺寸或核的脚印（footprint）。如果提供 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 函数计算沿给定轴的导数。

拉普拉斯滤波器由所有轴上的二阶导数之和计算得出。因此，可以使用不同的二阶导数函数构建不同的拉普拉斯滤波器。为此，我们提供了一个通用函数，它接受一个函数参数来计算沿给定方向的二阶导数。

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

  有关 extra_arguments 和 extra_keywords 参数的使用示例，请参阅 generic_laplace 的文档。

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 的参数是一个 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]])
  

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

  调用 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 阶的样条进行插值需要预滤波步骤。在 插值函数 节中描述的插值函数通过调用 spline_filter 来应用预滤波，但可以通过将 prefilter 关键字设置为 False 来指示它们不执行此操作。如果对同一个数组执行多个插值操作，这非常有用。在这种情况下，仅执行一次预滤波并将预滤波后的数组用作插值函数的输入会更有效。以下两个函数实现了预滤波：

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

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

  注意

  多维滤波器被实现为一系列一维样条滤波器。中间数组以与输出相同的数据类型存储。因此，如果请求精度有限的输出，结果可能不精确，因为中间结果可能以不足的精度存储。这可以通过指定高精度的输出类型来防止。

插值边界处理

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

图像样本的坐标落在沿每个轴 i 的 0 到 shape[i] - 1 范围内的整数采样位置上。下图说明了在形状为 (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 函数。）由于坐标可以是非整数坐标，输入在这些坐标处的值由请求阶数的样条插值确定。

  这是一个在 (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 函数对输入数组应用仿射变换。给定的变换 matrix（矩阵）和 offset（偏移）用于为输出中的每个点查找输入中对应的坐标。计算出的坐标处输入的值由请求阶数的样条插值确定。变换 matrix 必须是二维的，也可以作为一维序列或数组给出。在后一种情况下，假定矩阵是对角矩阵。然后会应用一种更高效的插值算法，该算法利用了问题的可分离性。输出形状和输出类型可以可选地提供。如果未给出，它们与输入形状和类型相等。

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

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

• rotate 函数使用请求 order 的样条插值，返回在由 axes 参数给定的两个轴定义的平面中旋转后的输入数组。角度必须以度为单位给出。如果 reshape 为 true，则调整输出数组的大小以包含旋转后的输入。

形态学

二值形态学

• generate_binary_structure 函数生成一个二值结构元素，用于二值形态学操作。必须提供结构的 rank（秩/维度）。返回的结构在每个方向上的大小等于 3。如果元素到中心的欧几里得距离平方小于或等于 connectivity（连通性），则每个元素的值等于 1。例如，二维 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 函数使用给定的结构元素实现任意秩数组的二值腐蚀。origin 参数控制结构元素的放置位置，如 滤波函数 中所述。如果未提供结构元素，则使用 generate_binary_structure 生成连通性等于 1 的元素。border_value 参数给出数组边界外的值。腐蚀重复 iterations 次。如果 iterations 小于 1，则重复腐蚀直到结果不再改变。如果给出了 mask 数组，则在每次迭代时仅修改对应掩码元素为 true 的元素。

• binary_dilation 函数使用给定的结构元素实现任意秩数组的二值膨胀。origin 参数控制结构元素的放置位置，如 滤波函数 中所述。如果未提供结构元素，则使用 generate_binary_structure 生成连通性等于 1 的元素。border_value 参数给出数组边界外的值。膨胀重复 iterations 次。如果 iterations 小于 1，则重复膨胀直到结果不再改变。如果给出了 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 函数通过将输入结构与其自身膨胀 iteration - 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)
  

  如果原始结构的原点等于 0，则迭代后的结构的原点也等于 0。如果不是，如果要通过迭代结构实现相当于 iterations 次腐蚀或膨胀的效果，则原点也必须进行调整。调整后的原点只需通过乘以迭代次数即可获得。为了方便起见，如果 origin 参数不是 None，iterate_structure 还会返回调整后的原点：

  >>> 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 函数使用给定的结构元素实现任意秩数组的二值开运算。二值开运算相当于二值腐蚀后跟一个使用相同结构元素的二值膨胀。origin 参数控制结构元素的放置位置。如果未提供结构元素，则生成连通性等于 1 的元素。iterations 参数给出了执行腐蚀的次数，随后是相同次数的膨胀。

• binary_closing 函数使用给定的结构元素实现任意秩数组的二值闭运算。二值闭运算相当于二值膨胀后跟一个使用相同结构元素的二值腐蚀。origin 参数控制结构元素的放置位置。如果未提供结构元素，则生成连通性等于 1 的元素。iterations 参数给出了执行膨胀的次数，随后是相同次数的腐蚀。

• binary_fill_holes 函数用于填充二值图像中对象的孔洞，其中结构定义了孔洞的连通性。origin 参数控制结构元素的放置。如果未提供结构元素，则生成连通性等于 1 的元素。

• binary_hit_or_miss 函数使用给定的结构元素实现任意秩数组的二值击中-击不中变换（hit-or-miss transform）。击中-击不中变换的计算方法是：输入与第一个结构的腐蚀，输入逻辑“非”与第二个结构的腐蚀，然后对这两个腐蚀结果进行逻辑“与”。origin 参数控制结构元素的放置。如果 origin2 等于 None，则将其设置为等于 origin1。如果未提供第一个结构元素，则生成连通性等于 1 的结构。如果未提供 structure2，则将其设置为 structure1 的逻辑“非”。

灰度形态学

灰度形态学操作是操作任意值数组的二值形态学操作的等效形式。下面，我们描述了腐蚀、膨胀、开运算和闭运算的灰度等效操作。这些操作的实现方式与 滤波函数 中描述的滤波器类似，有关滤波核、脚印以及数组边界处理的描述请参考该节。灰度形态学操作可选地接受一个 structure 参数，给出结构元素的值。如果未给出此参数，则假定结构元素是值为 0 的扁平结构。结构的形状可以可选地由 footprint 参数定义。如果未给出此参数，则假定结构是矩形的，其大小等于 structure 数组的维度，或者如果未给出 structure，则由 size 参数定义。仅在 structure 和 footprint 均未给出的情况下才使用 size 参数，在这种情况下，假定结构元素是矩形且扁平的，其维度由 size 给出。如果提供 size 参数，它必须是尺寸序列或单个数字。如果提供 footprint 参数，它必须是一个通过其非零元素定义核形状的数组。

类似于二值腐蚀和膨胀，也有灰度腐蚀和膨胀的操作：

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

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

灰度开运算和闭运算可以类似于其二值对应操作进行定义：

• grey_opening 函数实现任意秩数组的灰度开运算。灰度开运算相当于灰度腐蚀后跟一个灰度膨胀。

• grey_closing 函数实现任意秩数组的灰度闭运算。灰度闭运算相当于灰度膨胀后跟一个灰度腐蚀。

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

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

• white_tophat 函数实现任意秩数组的白顶帽（white top-hat）滤波。白顶帽等于输入与灰度开运算的差值。

• black_tophat 函数实现任意秩数组的黑顶帽（black top-hat）滤波。黑顶帽等于灰度闭运算与输入的差值。

距离变换

距离变换用于计算从对象的每个元素到背景的最小距离。以下函数实现了三种不同距离度量的距离变换：欧几里得距离、曼哈顿距离（city block）和切比雪夫距离（chessboard）。

• distance_transform_cdt 函数使用倒角（chamfer）类型算法计算输入的距离变换，方法是将每个对象元素（定义为大于零的值）替换为到背景（所有非对象元素）的最短距离。结构决定了执行的倒角类型。如果结构等于 'cityblock'，则使用 generate_binary_structure 生成平方距离等于 1 的结构。如果结构等于 'chessboard'，则生成平方距离等于数组秩的结构。这些选择对应于二维中曼哈顿距离和切比雪夫距离度量的常用解释。

  除了距离变换之外，还可以计算特征变换。在这种情况下，最近的背景元素的索引将沿结果的第一轴返回。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]], dtype=int32), 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]], dtype=int32)
  

  如果没有提供结构元素，则会通过调用 generate_binary_structure（参见 二值形态学）并使用连通度为 1 的设置来生成一个（在 2D 中即为第一个示例中的 4 连通结构）。输入可以是任何类型，任何不等于零的值都被视为对象的一部分。如果你需要对对象索引数组进行“重新标记”（例如在移除不需要的对象后），这非常有用。只需再次对索引数组应用 label 函数即可。例如：

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

  注意

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

分割还有许多其他方法，例如通过导数滤波器获得的物体边界估计。其中一种方法是分水岭分割。函数 watershed_ift 会根据定位物体边界的数组（例如由梯度幅度滤波器生成）生成一个数组，其中每个物体都被分配一个唯一的标记。它使用一个包含物体初始标记（markers）的数组。

• watershed_ift 函数应用了基于图像森林变换（Image Foresting Transform）的标记分水岭算法，如 [4] 中所述。

• 该函数的输入是应用变换的数组，以及一个通过唯一标记指定物体的 markers 数组，其中任何非零值都是一个标记。例如：

  >>> 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（参见 二值形态学）并使用连通度为 1 的设置来生成一个（在 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]])
  

  除非 max_label 参数大于零（此时仅返回前 max_label 个对象），否则 find_objects 函数返回所有对象的切片。如果 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 类型的一维数组形式返回。

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

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

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 函数进行交互。首先，我们需要编写共享库，我们用 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 时，C 代码和 so/DLL 的编译与上述 cffi 相同。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 文档中找到更多信息。

参考文献

[1]

M. Unser, “Splines: A Perfect Fit for Signal and Image Processing,” IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, November 1999.

[2]

G. Borgefors, “Distance transformations in arbitrary dimensions.”, Computer Vision, Graphics, and Image Processing, 27:321-345, 1984.

[3]

C. R. Maurer, Jr., R. Qi, and V. Raghavan, “A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions.” IEEE Trans. PAMI 25, 265-270, 2003.

[4]

A. X. Falcão, J. Stolfi, and R. A. Lotufo. “The image foresting transform: Theory, algorithms, and applications.” IEEE Trans. PAMI 26, 19-29. 2004.

[5]

T. Briand and P. Monasse, “Theory and Practice of Image B-Spline Interpolation”, Image Processing On Line, 8, pp. 99–141, 2018. https://doi.org/10.5201/ipol.2018.221