在numpy中，将值与相邻的元素进行比较

Question

假设我有一个numpy数组

    a b c
A = i j k
    u v w

我希望能够将中心元素的值与其八个相邻元素（沿轴线或对角线）进行比较。除了嵌套for循环之外，是否还有更快的方法（对于大矩阵而言，它太慢了）？

更具体地说，我想要做的是比较元素的值与其相邻元素，并赋予新值。

例如：

    if (j == 1):
      if (j>i) & (j>k):
        j = 999
      else:
        j = 0
    if (j == 2):
      if (j>c) & (j>u):
        j = 999
      else:
        j = 0
   ...

类似于这样的内容。

- Arthur

---

你能解释一下在矩阵尺寸较大的情况下，你希望程序如何运行吗？例如，你是否希望将每个元素与所有邻居进行比较？ - user20160

例如，将每个元素与其所有邻居进行比较，如果它是最大的，则将其分配给一个值，否则将其分配为零。 - Arthur

只有当中心值为1或2时，您才需要执行操作。这些值在您的矩阵中相对密集还是稀疏出现？ - user20160

这只是一个例子，在实际情况中它可能是密集的或稀疏的。谢谢。 - Arthur

Answer 1

你的操作包含许多条件语句，因此在一般情况下（任何类型的条件语句，任何类型的操作）最有效的方法是使用循环。可以使用numba或cython来高效地完成这项工作。在特殊情况下，你可以使用numpy/scipy中的更高级别函数进行实现。我将展示一个特定示例的解决方案，希望你可以从中推广。

首先使用一些虚假数据：

A = np.asarray([
    [1, 1, 1, 2, 0],
    [1, 0, 2, 2, 2],
    [0, 2, 0, 1, 0],
    [1, 2, 2, 1, 0],
    [2, 1, 1, 1, 2]
])

我们将在A中找到符合各种条件的位置。

• 1a）该值为1
• 1b）该值大于其水平邻居
• 2a）该值为2
• 2b）该值大于其对角线邻居

查找A中出现指定值的位置：

cond1a = A == 1
cond2a = A == 2

这将生成一个布尔值矩阵，其大小与A相同。当条件成立时，值为true，否则为false。

查找A中每个元素与其相邻元素之间的指定关系的位置：

# condition 1b: value greater than horizontal neighbors
f1 = np.asarray([[1, 0, 1]])
cond1b = A > scipy.ndimage.maximum_filter(
    A, footprint=f1, mode='constant', cval=-np.inf)

# condition 2b: value greater than diagonal neighbors
f2 = np.asarray([
    [0, 0, 1],
    [0, 0, 0],
    [1, 0, 0]
])
cond2b = A > scipy.ndimage.maximum_filter(
    A, footprint=f2, mode='constant', cval=-np.inf)

与之前类似，这会给出一个由布尔值矩阵组成的结果，指示条件的真假。此代码使用 scipy.ndimage.maximum_filter()。该函数迭代地将一个“足印”移动到A的每个元素中心。该位置返回值为足印元素值等于1的所有元素的最大值。参数mode指定如何处理矩阵边缘外的隐式值，即足印超出矩阵范围时的情况。在这里，我们将它们视为负无穷大，这与忽略它们相同（因为我们使用了最大操作）。

根据条件设置结果的值。如果条件1a和1b都成立，或者如果条件2a和2b都成立，则该值为999。否则，该值为0。

result = np.zeros(A.shape)
result[(cond1a & cond1b) | (cond2a & cond2b)] = 999

结果如下：

[
    [  0,   0,   0,   0,   0],
    [999,   0,   0, 999, 999],
    [  0,   0,   0, 999,   0],
    [  0,   0, 999,   0,   0],
    [  0,   0,   0,   0, 999]
]

您可以通过更改过滤器印记来将此方法推广到其他邻居模式。您可以使用其他类型的过滤器（请参见 scipy.ndimage）将其推广到其他操作（最小值、中位数、百分位数等）。对于可以表示为加权和的操作，请使用二维交叉相关。

这种方法应该比在Python中循环要快得多。但是，它确实执行了不必要的计算（例如，仅当值为1或2时才需要计算最大值，但我们正在为所有元素进行计算）。手动循环可以让您避免这些计算。在Python中循环可能比这里的代码要慢得多。但是，在numba或cython中实现它可能会更快，因为这些工具生成编译代码。

- user20160

Answer 2

我使用了numpy的：

• concatenate 来填充零
• dstack 和 roll 来正确对齐

在不同的维度上应用custom_roll两次并减去原始值。

import numpy as np

def custom_roll(a, axis=0):
    n = 3
    a = a.T if axis==1 else a

    pad = np.zeros((n-1, a.shape[1]))
    a = np.concatenate([a, pad], axis=0)
    ad = np.dstack([np.roll(a, i, axis=0) for i in range(n)])
    a = ad.sum(2)[1:-1, :]

    a = a.T if axis==1 else a
    return a

考虑以下的 ndarray：

A = np.arange(25).reshape(5, 5)
A

array([[ 0,  1,  2,  3,  4],
       [ 5,  6,  7,  8,  9],
       [10, 11, 12, 13, 14],
       [15, 16, 17, 18, 19],
       [20, 21, 22, 23, 24]])


sum_of_eight_around_me = custom_roll(custom_roll(A), axis=1) - A

sum_of_eight_around_me

array([[  12.,   20.,   25.,   30.,   20.],
       [  28.,   48.,   56.,   64.,   42.],
       [  53.,   88.,   96.,  104.,   67.],
       [  78.,  128.,  136.,  144.,   92.],
       [  52.,   90.,   95.,  100.,   60.]])

- piRSquared

---

有没有一种只使用numpy和scipy而不需要其他包的方法来完成它？ - Arthur

仅更新到NumPy的帖子。 - piRSquared

除了求和之外，有没有其他逻辑运算的方式呢？例如：如果 j 比它的邻居大？ - Arthur

是的，但是你的问题不够明确。我用求和的具体示例说明了如何做到这一点。如果你想要更符合你需求的东西，请提出更具体的问题。 - piRSquared

抱歉，我更新了问题。感谢您的帮助。 - Arthur