多维图像处理 (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]])

通常，核的原点位于通过将其核形状的维度除以二来计算的中心。例如，长度为三的一维核的原点在第二个元素处。以一维数组与由全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 函数计算沿给定轴的一维相关。沿给定轴的数组行与给定 weights 相关联。weights 参数必须是一维数字序列。

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

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

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

    注意

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

平滑滤波器#

  • gaussian_filter1d 函数实现一维高斯滤波器。高斯滤波器的标准差通过参数 sigma 传递。设置 order = 0 对应于与高斯核卷积。阶数为 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 函数计算多维最小滤波器。必须提供矩形核的大小或核的足迹。如果提供 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 函数计算沿给定轴的导数。

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

  • 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 参数可用于传递一个额外参数的元组和/或一个命名参数的字典，这些参数和命名参数在每次调用时都会传递给 derivative。例如，我们可以将滤波器的参数作为参数传递

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

    这里，指定了一个只包含两个元素的核足迹。因此，滤波函数接收到一个长度为二的缓冲区，该缓冲区乘以适当的权重，然后结果求和。

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

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

    >>> 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 个样本。实心圆圈表示在红色叉号位置处值的插值所涉及的采样位置。

" "

插值函数#

  • 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 参数可用于传递一个额外参数的元组和/或一个命名参数的字典，这些参数和命名参数在每次调用时都会传递给 derivative。例如，我们可以将示例中的移位作为参数传递

    >>> 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。返回的结构大小在每个方向上都等于三。如果从元素到中心的欧几里得距离的平方小于或等于 connectivity，则每个元素的值等于一。例如，二维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 生成一个连通性等于一的元素。border_value 参数给出边界外数组的值。腐蚀重复 iterations 次。如果 iterations 小于一，则腐蚀重复直到结果不再改变。如果给定 mask 数组，则在每次迭代时只修改那些在相应掩码元素处具有真值的元素。

  • binary_dilation 函数使用给定的结构元素实现任意秩数组的二值膨胀。origin 参数控制结构元素的放置，如滤波函数中所述。如果没有提供结构元素，则使用 generate_binary_structure 生成一个连通性等于一的元素。border_value 参数给出边界外数组的值。膨胀重复 iterations 次。如果 iterations 小于一，则膨胀重复直到结果不再改变。如果给定 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 生成一个连通性等于一的元素。iterations 参数给出执行腐蚀的次数，后接相同次数的膨胀。

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

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

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

灰度形态学#

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

与二值腐蚀和膨胀类似，灰度腐蚀和膨胀也有相应的操作

灰度开运算和闭运算可以与它们的二值对应物类似地定义

  • 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（参见 二值形态学）并使用连通性为一（在二维中是第一个示例的 4 连通结构）来生成一个。输入可以是任何类型，任何不等于零的值都被视为对象的一部分。这在您需要“重新标记”对象索引数组时很有用，例如在移除不需要的对象之后。只需将标签函数再次应用于索引数组即可。例如

    >>> 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 函数应用基于标记的分水岭算法，使用图像森林变换（Image Foresting Transform），如 [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（参见 二值形态学）并使用连通性为一（在二维中是 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 扩展，以及一个返回包含指向回调指针的 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 回调函数返回一个整数状态，成功时为一，否则为零。

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

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

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

  • 上面，我们使用 scipy.LowLevelCallable 来指定通过 ctypes 生成的 user_data。

    另一种方法是在胶囊上下文（可以通过 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 文档中找到更多信息。

参考文献#