具有多维（或非标量）输出的 Scipy 过滤器答案

【问题标题】：Scipy filter with multi-dimensional (or non-scalar) output具有多维（或非标量）输出的 Scipy 过滤器
【发布时间】：2015-04-30 17:31:22
【问题描述】：

有没有类似ndimage的generic_filter的过滤器支持向量输出？我没有设法让scipy.ndimage.filters.generic_filter 返回超过一个标量。取消注释下面代码中的行以获取错误：TypeError: only length-1 arrays can be converted to Python scalars。

我正在寻找一个处理 2D 或 3D 数组并在每个点返回一个向量的通用过滤器。因此，输出将增加一个维度。对于下面的示例，我希望是这样的：

m.shape    # (10,10)
res.shape  # (10,10,2)

示例代码

import numpy as np
from scipy import ndimage

a = np.ones((10, 10)) * np.arange(10)

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

def myfunc(x):
    r = sum(x)
    #r = np.array([1,1])  # uncomment this
    return r

res = ndimage.generic_filter(a, myfunc, footprint=footprint)

【问题讨论】：

标签： python image-processing numpy scipy ndimage

【解决方案1】：

generic_filter 期望 myfunc 返回一个标量，而不是向量。但是，没有什么可以阻止myfunc也添加信息例如，一个作为额外参数传递给myfunc 的列表。

不使用generic_filter返回的数组，我们可以通过重塑这个列表来生成我们的向量值数组。

例如，

import numpy as np
from scipy import ndimage

a = np.ones((10, 10)) * np.arange(10)

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

ndim = 2
def myfunc(x, out):
    r = np.arange(ndim, dtype='float64')
    out.extend(r)
    return 0

result = []
ndimage.generic_filter(
    a, myfunc, footprint=footprint, extra_arguments=(result,))
result = np.array(result).reshape(a.shape+(ndim,))

【讨论】：

使用“可变默认值”确实是一个聪明的技巧，我还没有看到。我认为您的解决方案优雅简洁，并回答了我的问题。你能评论一下表演吗？看起来除了result 变量的操作之外没有任何开销。另外，r[1:][::-1] 是否等同于 r[:0:-1]？
generic_filter 的性能从来都不是很好，因为它需要为图像的每个像素调用一次 Python 函数。我会尝试仅将它用作最后的手段——在寻找一种方法来表达计算后，使用更快的内置 NumPy/SciPy 函数应用于整个数组，这些函数将大部分工作卸载到用 C 编写的底层函数/C++ 或 Fortran。
你是对的r[:0:-1] 相当于r[1:][::-1]。由于我很容易混淆，我选择了（在我看来）更简单的r[1:][::-1]。 r[:0:-1] 有点快，但要注意预优化。 myfunc 的真正瓶颈可能在于r 的计算，而不是逆向它。
如果您追求速度并愿意牺牲准确性，您可以使用另一种技巧：myfunc 必须返回一个标量，例如 np.float128。但是NumPy也支持np.float16，所以你可以将8个np.float16s打包成一个np.float128。因此，如果您想返回向量 np.arange(8, dtype=np.float16)，则可以改为返回单个标量 np.arange(8, dtype=np.float16).view(np.float128)，然后使用 res = res.view(np.float16) 恢复 8 个 np.float16 值。还有一个np.complex256dtype，可以让你打包16个np.float16s。
我找到了一种更简单的方法来达到预期的效果。只需将所有向量附加到列表中。 generic_filter 完成后，将列表转换为数组，然后对其进行整形。这也比我的第一个答案略快。

【解决方案2】：

我想我明白了你的要求，但我不完全确定 ndimage.generic_filter 是如何工作的（来源是多么深奥！）。

这里只是一个简单的包装函数。这个函数将接受一个数组，ndimage.generic_filter 需要的所有参数。函数返回一个数组，其中前一个数组的每个 element 现在由一个形状为 (2,) 的 array 表示，函数的结果存储为该数组的第二个元素数组。

def generic_expand_filter(inarr, func, **kwargs):
    shape = inarr.shape
    res = np.empty((  shape+(2,) ))
    temp = ndimage.generic_filter(inarr, func, **kwargs)
    for row in range(shape[0]):
        for val in range(shape[1]):
            res[row][val][0] = inarr[row][val]
            res[row][val][1] = temp[row][val]
    return res

输出，其中res 仅表示generic_filter，res2 表示generic_expand_filter，此函数为：

>>> a.shape #same as res.shape
(10, 10)
>>> res2.shape
(10, 10, 2)

>>> a[0]
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.])
>>> res[0]
array([  3.,   8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,  69.])
>>> print(*res2[0], sep=", ") #this is just to avoid the vertical default output
[ 0.  3.], [ 1.  8.], [  2.  16.], [  3.  24.], [  4.  32.], [  5.  40.], [  6.  48.], [  7.  56.], [  8.  64.], [  9.  69.]

>>> a[0][0]
0.0
>>> res[0][0]
3.0
>>> res2[0][0]
array([ 0.,  3.])

当然，您可能不想保存旧数组，而是将两个字段都作为新结果。除非我不知道您到底在想什么，如果您要存储的两个值不相关，只需添加一个 temp2 和 func2 并使用相同的 generic_filter 调用另一个 generic_filter 并将其存储为第一个值。

但是，如果您想要使用多个 inarr 元素计算的实际向量数量，这意味着两个新创建的字段不是独立的，您只需要编写那种函数，一个需要在数组中，idx，idy 索引并返回一个元组\列表\数组值，然后您可以将其解包并分配给结果。

【讨论】：