在numpy数组上用二进制结构迭代得到单元和

Question

在包scipy中，有一个定义二进制结构的函数（例如taxicab（2,1）或chessboard（2,2））。在

import numpy
from scipy import ndimage
a = numpy.zeros((6,6), dtype=numpy.int) 
a[1:5, 1:5] = 1;a[3,3] = 0 ; a[2,2] = 2
s = ndimage.generate_binary_structure(2,2) # Binary structure
#.... Calculate Sum of 
result_array = numpy.zeros_like(a)

我想要的是用给定的结构s遍历这个数组的所有单元格，然后在空数组中索引的当前单元格值上附加一个函数（例如sum函数），它使用二进制结构中所有单元格的值。在

例如：

^{pr2}$

#数组a。单元格1,2中的值当前为1。给定结构s和sum等示例函数，结果数组（result_array）中的值变为7（如果当前单元格值被排除，则为6）。在

有人有主意吗？在

Answer 1

对于求和的特殊情况，可以使用ndimage.convolve：

In [42]: import numpy as np

In [43]: a = np.zeros((6,6), dtype=np.int) 
a[1:5, 1:5] = 1;
a[3,3] = 0;
a[2,2] = 2

In [48]: s = ndimage.generate_binary_structure(2,2) # Binary structure

In [49]: ndimage.convolve(a,s)
Out[49]: 
array([[1, 2, 3, 3, 2, 1],
       [2, 5, 7, 7, 4, 2],
       [3, 7, 9, 9, 5, 3],
       [3, 7, 9, 9, 5, 3],
       [2, 4, 5, 5, 3, 2],
       [1, 2, 3, 3, 2, 1]])

对于产品的特定情况，可以使用log(a*b) = log(a)+log(b)将问题转换回涉及和的问题。例如，如果我们想要“乘积卷积”b：

^{pr2}$

我们可以计算：

print(np.exp(ndimage.convolve(np.log(b), s, mode = 'constant')))
# [[ 2.  2.  2.  1.]
#  [ 2.  0.  0.  0.]
#  [ 2.  0.  0.  0.]
#  [ 1.  0.  0.  0.]]

如果b包含负值，情况会变得更复杂：

b[0,1] = -1
print(b)
# [[ 1 -1  1  1]
#  [ 1  2  1  1]
#  [ 1  1  0  1]
#  [ 1  1  1  1]]

但并非不可能：

logb = np.log(b.astype('complex'))
real, imag = logb.real, logb.imag
print(np.real_if_close(
    np.exp(
        sum(j * ndimage.convolve(x, s, mode = 'constant')
            for x,j in zip((real, imag),(1,1j))))))
# [[-2. -2. -2.  1.]
#  [-2. -0. -0.  0.]
#  [ 2.  0.  0.  0.]
#  [ 1.  0.  0.  0.]]

Answer 2

如果使用2深的零墙，则更容易：

In [11]: a0
Out[11]: 
array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  2.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

In [12]: b0 = zeros_like(a0)

In [13]: for i in range(1,len(a0)-1):
   ....:     for j in range(1,len(a0)-1):
   ....:         b0[i,j] = sum(a0[i-1:i+2, j-1:j+2] * s)

这样可以根据需要将两个子矩阵相乘并求和。（你也可以在这里做一些更详细的事情…）

^{pr2}$