多维图像处理 (scipy.ndimage)#

简介#

图像处理和分析通常被视为对二维数组的值进行的操作。但是，在许多领域，必须分析更高维度的图像。医学成像和生物成像就是很好的例子。 numpy 非常适合这种应用，因为它具有固有的多维特性。 scipy.ndimage 包提供了许多通用图像处理和分析函数，这些函数旨在对任意维度的数组进行操作。该软件包目前包括：线性与非线性滤波、二值形态学、B 样条插值和对象测量函数。

所有函数的共有属性#

所有函数都共享一些共同的属性。值得注意的是，所有函数都允许使用 output 参数指定输出数组。使用此参数，您可以指定一个将被结果更改的数组。在这种情况下，不会返回结果。通常，使用 output 参数更有效率，因为使用现有数组来存储结果。

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

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

滤波函数#

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

>>> 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 参数选择边界条件，该参数必须是包含边界条件名称的字符串。目前支持以下边界条件

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 等同于用高斯内核进行卷积。1、2 或 3 的阶数对应于用高斯的第一个、第二个或第三个导数进行卷积。未实现更高阶导数。

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

    " "

    注意

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

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

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

    注意

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

基于顺序统计的滤波器#

  • minimum_filter1d 函数计算给定轴上给定大小的一维最小滤波器。

  • maximum_filter1d 函数计算给定轴上给定大小的一维最大滤波器。

  • 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 实现了一个通用的 1-D 滤波函数，其中实际的滤波操作必须作为 Python 函数（或其他可调用对象）提供。函数 generic_filter1d 遍历数组的行，并在每行调用 function。传递给 function 的参数是 numpy.float64 类型的 1-D 数组。第一个包含当前行的值。它根据 *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 类型的 1-D 数组，它包含当前元素周围在滤波器足迹内的值。该函数应该返回一个可以转换为双精度数字的单个值。例如，考虑一个相关性

    >>> 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-D 样条滤波器。可以选择提供输出数组。样条的阶数必须大于 1 小于 6。

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

    注意

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

插值边界处理#

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

" "

图像样本的坐标落在每个轴上从 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参数可以用来传递一个额外的参数元组和/或一个命名的参数字典，这些参数和/或字典在每次调用时都会传递给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返回由参数axes给出的两个轴定义的平面内旋转的输入数组，使用请求的order的样条插值。角度必须以度为单位给出。如果reshape为真，则输出数组的大小将调整为包含旋转后的输入。

形态学#

二值形态学#

  • 函数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数组，则在每次迭代中只修改那些在相应掩码元素处具有真值的元素。

  • 函数binary_dilation使用给定的结构元素实现任意秩数组的二值膨胀。origin 参数控制结构元素的位置，如过滤器函数中所述。如果没有提供结构元素，则使用generate_binary_structure生成一个连接性等于 1 的元素。border_value参数给出边界之外数组的值。膨胀重复iterations次。如果iterations小于 1，则膨胀重复，直到结果不再改变。如果给出了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 数组的尺寸，或者由 size 参数定义（如果未提供 structure）。 size 参数仅在未提供 structure 和 footprint 时使用，在这种情况下，假定结构元素是矩形的且平坦的，其尺寸由 size 给出。 如果提供 size 参数，则它必须是一个尺寸序列或一个单个数字，在这种情况下，假定滤波器的尺寸沿每个轴都相等。 如果提供 footprint 参数，则它必须是一个数组，该数组通过其非零元素来定义内核的形状。

与二进制腐蚀和膨胀类似，存在灰度腐蚀和膨胀运算。

灰度开运算和闭运算可以类似于它们的二进制对应物来定义。

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

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

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

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

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

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

距离变换#

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

  • 函数 distance_transform_cdt 使用一种 chamfer 型算法来计算输入的距离变换，方法是将每个对象元素（由大于零的值定义）替换为到背景（所有非对象元素）的最短距离。 结构决定了进行的 chamfering 类型。 如果结构等于 ‘cityblock’，则使用 generate_binary_structure 生成一个结构，其平方距离等于 1。 如果结构等于 ‘chessboard’，则使用 generate_binary_structure 生成一个结构，其平方距离等于数组的秩。 这些选择对应于二维空间中城市街区距离和棋盘距离的常见解释。

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

    distances 和 indices 参数可用于给出必须具有正确大小和类型（两者均为 numpy.int32）的可选输出数组。 用于实现此函数的算法的基本原理在 [2] 中进行了描述。

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

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

    可以选择由 sampling 参数给出沿每个轴的采样，该参数应该是长度等于输入秩的序列，或者是一个单个数字，在这种情况下，假定沿所有轴的采样都相等。

    distances 和 indices 参数可用于给出必须具有正确大小和类型（numpy.float64 和 numpy.int32）的可选输出数组。 用于实现此函数的算法在 [3] 中进行了描述。

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

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

    可选地，可以使用sampling参数指定每个轴的采样方式，该参数应为与输入秩相等的长度的序列，或者是一个单一数字，在这种情况下，假定所有轴的采样相同。此参数仅在使用欧几里德距离变换的情况下使用。

    distances和indices参数可用于提供可选的输出数组，这些数组必须具有正确的尺寸和类型（numpy.float64 和 numpy.int32）。

    注意

    此函数使用慢速的暴力算法，函数distance_transform_cdt可用于更有效地计算城市街区和棋盘距离变换。函数distance_transform_edt可用于更有效地计算精确的欧几里德距离变换。

分割和标记#

分割是指将感兴趣的目标从背景中分离出来。最简单的方法可能是强度阈值化，这可以使用numpy函数轻松完成

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

结果是二值图像，其中单个目标仍然需要被识别和标记。函数label生成一个数组，其中每个目标都分配了一个唯一的数字

  • 函数label生成一个数组，其中输入中的目标用整数索引标记。它返回一个包含目标标签数组和找到的目标数量的元组，除非提供了output参数，在这种情况下只返回目标数量。目标的连通性由结构元素定义。例如，在二维空间中使用4连通结构元素会导致

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

    这两个目标没有连接，因为我们无法放置结构元素，使其与两个目标重叠。但是，8连通结构元素会导致只有一个目标

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

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

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

    注意

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

还有许多其他分割方法，例如，通过估计目标的边界可以获得的目标边界，这些边界可以通过导数滤波器获得。一种这样的方法是分水岭分割。函数watershed_ift生成一个数组，其中每个目标都分配了一个唯一的标签，该标签来自一个定位目标边界的数组，例如由梯度幅度滤波器生成。它使用一个包含目标初始标记的数组

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

  • 此函数的输入是应用变换的数组，以及一个标记数组，该数组用唯一的标签指定目标，其中任何非零值都是标记。例如

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

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

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

    结果是目标（marker = 2）更小，因为第二个标记更早被处理。如果第一个标记应该指定背景目标，这可能不是预期效果。因此，watershed_ift将负值的标记显式地视为背景标记，并在处理普通标记之后处理它们。例如，将第一个标记替换为负标记会产生与第一个示例类似的结果

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

    目标的连通性由结构元素定义。如果没有提供结构元素，则通过调用generate_binary_structure（参见二值形态学）使用连通性1（在二维空间中是4连通结构）来生成结构元素。例如，在最后一个示例中使用8连通结构会产生不同的目标

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

    注意

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

目标测量#

给定一个标记目标的数组，可以测量单个目标的属性。函数find_objects可用于生成一个切片列表，该列表对于每个目标，给出完全包含该目标的最小子数组

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

    例如

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

    函数find_objects返回所有目标的切片，除非max_label参数大于零，在这种情况下只返回前max_label个目标。如果label数组中缺少索引，则返回None，而不是切片。例如

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

函数find_objects生成的切片列表有助于查找数组中目标的位置和尺寸，但也可以用于对单个目标进行测量。假设我们要查找图像中目标的强度之和

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

然后，我们可以计算第二个目标中元素的总和

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

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

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

对于大型数组和小目标，在对数组进行切片后调用测量函数效率更高

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

或者，我们可以对多个标签执行测量，并使用单个函数调用返回结果列表。例如，要测量我们示例中背景和第二个目标的值之和，我们给出标签列表

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

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

  • 函数sum计算index给出的标签的目标元素的总和，使用labels数组作为目标标签。如果index为None，则所有标签值不为零的元素都被视为单个目标。如果label为None，则使用input的所有元素进行计算。

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

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

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

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

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

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

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

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

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

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

用 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 文档中找到更多信息。

参考文献#

本页内容