多维图像处理 (scipy.ndimage)

简介

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

所有函数共享的属性

所有函数都共享一些共同的属性。值得注意的是，所有函数都允许使用 output 参数指定输出数组。使用此参数，您可以指定一个数组，该数组将通过操作就地更改为结果。在这种情况下，不会返回结果。通常，使用 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])

滤波函数

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

>>> 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 的一维内核的原点位于第二个元素。例如，取一维数组与由 1 组成的长度为 3 的滤波器的相关性

>>> 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 相关性。沿给定轴的数组线与给定的 weights 相关。weights 参数必须是一个一维的数字序列。

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

• convolve1d 函数计算沿给定轴的 1-D 卷积。沿给定轴的数组线与给定的 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。

  

  注意

  多维滤波器实现为一系列 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 传递的函数来计算 Laplace 滤波器，以计算二阶导数。函数 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 使用离散微分来计算 Laplace（即，与 [1, -2, 1] 的卷积）。

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

梯度幅度定义为所有方向上梯度平方和的平方根。与通用 Laplace 函数类似，有一个 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函数实现了一个通用的 1D 滤波器函数，其中实际的滤波操作必须作为 Python 函数（或其他可调用对象）提供。generic_filter1d函数迭代数组的每一行，并在每一行调用function。传递给function的参数是numpy.float64类型的 1D 数组。第一个数组包含当前行的值。根据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类型的 1D 数组，其中包含当前元素周围的、在滤波器覆盖区内的值。该函数应返回可以转换为双精度数的单个值。例如，考虑一个相关性

  >>> 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函数计算沿给定轴的 1D 样条滤波器。可以可选地提供输出数组。样条的阶数必须大于 1 且小于 6。

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

  注意

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

插值边界处理

插值函数都使用样条插值来实现输入数组的某种几何变换。这需要将输出坐标映射到输入坐标，因此，可能会需要边界之外的输入值。这个问题与 滤波器函数 中描述的多维滤波器函数采用相同的方式解决。因此，这些函数都支持一个 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) 处对 2D 数组进行插值的示例

  >>> 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 必须是 2-D 的，也可以作为 1-D 序列或数组给出。在后一种情况下，假定矩阵是对角矩阵。然后应用更有效的插值算法，该算法利用问题的可分离性。输出形状和输出类型可以可选地提供。如果未给出，则它们等于输入形状和类型。

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

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

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

形态学

二值形态学

• generate_binary_structure 函数生成一个二值结构元素，用于二值形态学操作。必须提供结构的 rank。返回的结构的大小在每个方向上都等于 3。如果从元素到中心的欧几里德距离的平方小于或等于 connectivity，则每个元素的值等于 1。例如，2-D 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 数组，则每次迭代时仅修改在对应的 mask 元素处具有真值的元素。

• binary_dilation 函数使用给定的结构元素实现任意秩的数组的二值膨胀。origin 参数控制结构元素的放置，如 滤波器函数 中所述。如果未提供结构元素，则使用 generate_binary_structure 生成一个连通性等于 1 的元素。border_value 参数给出边界之外数组的值。膨胀重复 iterations 次。如果 iterations 小于 1，则重复膨胀，直到结果不再变化。如果给定了 mask 数组，则每次迭代时仅修改在对应的 mask 元素处具有真值的元素。

这是一个使用 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)
  
  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 函数使用给定的结构元素实现任意秩数组的二值开运算。二值开运算等效于使用相同结构元素进行二值腐蚀，然后进行二值膨胀。origin 参数控制结构元素的位置，如 滤波函数 中所述。如果没有提供结构元素，则使用 generate_binary_structure 生成一个连通性等于 1 的元素。iterations 参数指定执行的腐蚀次数，然后执行相同次数的膨胀。

• binary_closing 函数使用给定的结构元素实现任意秩数组的二值闭运算。二值闭运算等效于使用相同结构元素进行二值膨胀，然后进行二值腐蚀。origin 参数控制结构元素的位置，如 滤波函数 中所述。如果没有提供结构元素，则使用 generate_binary_structure 生成一个连通性等于 1 的元素。iterations 参数指定执行的膨胀次数，然后执行相同次数的腐蚀。

• binary_fill_holes 函数用于闭合二值图像中对象的孔洞，其中结构定义了孔洞的连通性。origin 参数控制结构元素的位置，如 滤波函数 中所述。如果没有提供结构元素，则使用 generate_binary_structure 生成一个连通性等于 1 的元素。

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

灰度形态学

灰度形态学运算是作用于任意值数组的二值形态学运算的等效运算。下面，我们将描述腐蚀、膨胀、开运算和闭运算的灰度等效运算。这些运算的实现方式与 滤波函数 中描述的滤波器类似，我们参考该部分了解滤波器核和足迹的描述，以及数组边界的处理。灰度形态学运算可以选择性地使用一个 structure 参数，该参数给出结构元素的值。如果没有给出此参数，则假定结构元素是平面的，值为零。结构的形状可以选择性地由 footprint 参数定义。如果未给出此参数，则假定结构是矩形的，大小等于 structure 数组的维度，或者如果未给出 structure 则由 size 参数给出。size 参数仅在未给出 structure 和 footprint 时使用，在这种情况下，假定结构元素是矩形的且是平面的，其尺寸由 size 给出。size 参数（如果提供）必须是一个大小序列或一个数字，在这种情况下，假定滤波器的大小在每个轴上都相等。footprint 参数（如果提供）必须是一个数组，该数组通过其非零元素定义核的形状。

与二值腐蚀和膨胀类似，也有用于灰度腐蚀和膨胀的运算

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

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

灰度开运算和闭运算可以类似于它们的二值运算来定义

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

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

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

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

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

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

距离变换

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

• 函数 distance_transform_cdt 使用倒角类型算法来计算输入的距离变换，方法是将每个对象元素（由大于零的值定义）替换为到背景（所有非对象元素）的最短距离。结构确定执行的倒角类型。如果结构等于 ‘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]], 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 生成一个数组，其中每个对象都从一个定位对象边界的数组（例如，通过梯度幅度滤波器生成）分配一个唯一的标签。它使用一个包含对象初始标记的数组

• 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 （请参见 二进制形态学）使用连通性 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]])
  

  函数 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 函数使用 labels 数组作为对象标签，计算由 index 给出的标签的对象元素的总和。如果 index 为 None，则所有标签值非零的元素都被视为单个对象。如果 label 为 None，则计算中使用 input 的所有元素。

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

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

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

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

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

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

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

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

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

• histogram 函数使用 labels 数组作为对象标签，计算由 index 给出的标签的对象的直方图。如果 index 为 None，则所有标签值非零的元素都被视为单个对象。如果 label 为 None，则计算中使用 input 的所有元素。直方图由其最小值（min）、最大值（max）和 bin 的数量（bins）定义。它们作为类型为 numpy.int32 的 1-D 数组返回。

在 C 中扩展 scipy.ndimage

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

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

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 context）中提供数据，这可以通过 PyCapsule_SetContext 设置，并且在 scipy.LowLevelCallable 中省略指定 user_data。但是，在这种方法中，我们需要处理数据的分配/释放——在胶囊被销毁后释放数据可以通过在 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