Numpy: 快速计算考虑数组内项的邻居及其位置

Question

我有4个2dnumpy数组，称为a, b, c, d，每个数组由n行和m列组成。我需要做的是给b和d的每个元素一个值，计算如下（伪代码）：

min_coords = min_of_neighbors_coords(x, y)
b[x,y] = a[x,y] * a[min_coords];
d[x,y] = c[min_coords];

其中min_of_neighbors_coords是一个函数，在给定数组元素的坐标后，返回值较小的“neighbor”元素的坐标。一、 e.考虑到阵列：

^{pr2}$

min_of_neighbors_coords(1, 1)将引用值为7的中心元素，并将返回元组(0, 0)：数字1的坐标。在

我设法用for循环（每个元素一个元素）来实现这一点，但是这个算法非常慢，我正在寻找一种方法来改进它，避免循环并要求计算达到numpy。在

有可能吗？在

Answer 1

编辑我将原始答案保留在底部。正如Paul在评论中指出的，最初的答案并没有真正回答OP的问题，使用ndimage过滤器可以更容易地实现。下面这个更麻烦的函数应该做正确的事情。它接受两个数组a和c，并返回a的加窗最小值和a中开窗最小值位置的值：

def neighbor_min(a, c):
    ac = np.concatenate((a[None], c[None]))
    rows, cols = ac.shape[1:]
    ret = np.empty_like(ac)

    # Fill in the center
    win_ac = as_strided(ac, shape=(2, rows-2, cols, 3),
                        strides=ac.strides+ac.strides[1:2])
    win_ac = win_ac[np.ogrid[:2, :rows-2, :cols] +
                    [np.argmin(win_ac[0], axis=2)]]
    win_ac = as_strided(win_ac, shape=(2, rows-2, cols-2, 3),
                        strides=win_ac.strides+win_ac.strides[2:3])
    ret[:, 1:-1, 1:-1] =  win_ac[np.ogrid[:2, :rows-2, :cols-2] +
                                 [np.argmin(win_ac[0], axis=2)]]

    # Fill the top, bottom, left and right borders
    win_ac = as_strided(ac[:, :2, :], shape=(2, 2, cols-2, 3),
                        strides=ac.strides+ac.strides[2:3])
    win_ac = win_ac[np.ogrid[:2, :2, :cols-2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, 0, 1:-1] = win_ac[:, np.argmin(win_ac[0], axis=0),
                             np.ogrid[:cols-2]]
    win_ac = as_strided(ac[:, -2:, :], shape=(2, 2, cols-2, 3),
                        strides=ac.strides+ac.strides[2:3])
    win_ac = win_ac[np.ogrid[:2, :2, :cols-2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, -1, 1:-1] = win_ac[:, np.argmin(win_ac[0], axis=0),
                             np.ogrid[:cols-2]]
    win_ac = as_strided(ac[:, :, :2], shape=(2, rows-2, 2, 3),
                        strides=ac.strides+ac.strides[1:2])
    win_ac = win_ac[np.ogrid[:2, :rows-2, :2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, 1:-1, 0] = win_ac[:, np.ogrid[:rows-2],
                             np.argmin(win_ac[0], axis=1)]
    win_ac = as_strided(ac[:, :, -2:], shape=(2, rows-2, 2, 3),
                        strides=ac.strides+ac.strides[1:2])
    win_ac = win_ac[np.ogrid[:2, :rows-2, :2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, 1:-1, -1] = win_ac[:, np.ogrid[:rows-2],
                             np.argmin(win_ac[0], axis=1)]
    # Fill the corners
    win_ac = ac[:, :2, :2]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, 0, 0] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
    win_ac = ac[:, :2, -2:]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, 0, -1] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
    win_ac = ac[:, -2:, -2:]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, -1, -1] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
    win_ac = ac[:, -2:, :2]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, -1, 0] = win_ac[:, np.argmin(win_ac[0], axis=-1)]

    return ret

返回的是一个(2, rows, cols)数组，可以将其解包到两个数组中：

^{pr2}$

然后，OP的案子可以解决为：

def bd_from_ac(a, c):
    b,d = neighbor_min(a, c)
    return a*b, d

虽然性能受到严重影响，但速度仍然相当快：

In [3]: a = np.random.rand(1000, 1000)

In [4]: c = np.random.rand(1000, 1000)

In [5]: %timeit bd_from_ac(a, c)
1 loops, best of 3: 570 ms per loop

---

除了获取元素之外，您并没有真正使用最小相邻元素的坐标，所以您可以跳过这一部分，创建一个min_neighbor函数。如果你不想借助cython来实现快速循环，你将不得不使用滚动窗口视图，比如在Paul的链接中概述的那样。这通常会将(m, n)数组转换为相同数据的(m-2, n-2, 3, 3)视图，然后在最后两个轴上应用np.min。在

不幸的是，您必须一次应用一个轴，因此您必须创建数据的(m-2, n-2, 3)副本。幸运的是，您可以通过两个步骤计算最小值，首先沿一个轴开窗并最小化，然后沿另一个轴最小化，并获得相同的结果。所以你最多会有一个输入大小的中间存储器。如果需要，您甚至可以将输出数组重用为中间存储并避免内存分配，但这是left as exercise。。。在

下面的函数可以做到这一点。它有点冗长，因为它不仅要处理中心区域，还要处理四个边和四个角的特殊情况。除此之外，它是一个非常紧凑的实现：

def neighbor_min(a):
    rows, cols = a.shape
    ret = np.empty_like(a)

    # Fill in the center
    win_a = as_strided(a, shape=(m-2, n, 3),
                       strides=a.strides+a.strides[:1])
    win_a = win_a.min(axis=2)
    win_a = as_strided(win_a, shape=(m-2, n-2, 3),
                       strides=win_a.strides+win_a.strides[1:])
    ret[1:-1, 1:-1] = win_a.min(axis=2)

    # Fill the top, bottom, left and right borders
    win_a = as_strided(a[:2, :], shape=(2, cols-2, 3),
                       strides=a.strides+a.strides[1:])
    ret[0, 1:-1] = win_a.min(axis=2).min(axis=0)
    win_a = as_strided(a[-2:, :], shape=(2, cols-2, 3),
                       strides=a.strides+a.strides[1:])
    ret[-1, 1:-1] = win_a.min(axis=2).min(axis=0)
    win_a = as_strided(a[:, :2], shape=(rows-2, 2, 3),
                       strides=a.strides+a.strides[:1])
    ret[1:-1, 0] = win_a.min(axis=2).min(axis=1)
    win_a = as_strided(a[:, -2:], shape=(rows-2, 2, 3),
                       strides=a.strides+a.strides[:1])
    ret[1:-1, -1] = win_a.min(axis=2).min(axis=1)

    # Fill the corners
    ret[0, 0] = a[:2, :2].min()
    ret[0, -1] = a[:2, -2:].min()
    ret[-1, -1] = a[-2:, -2:].min()
    ret[-1, 0] = a[-2:, :2].min()

    return ret

现在您可以执行以下操作：

>>> a = np.random.randint(10, size=(5, 5))
>>> a
array([[0, 3, 1, 8, 9],
       [7, 2, 7, 5, 7],
       [4, 2, 6, 1, 9],
       [2, 8, 1, 2, 3],
       [7, 7, 6, 8, 0]])
>>> neighbor_min(a)
array([[0, 0, 1, 1, 5],
       [0, 0, 1, 1, 1],
       [2, 1, 1, 1, 1],
       [2, 1, 1, 0, 0],
       [2, 1, 1, 0, 0]])

你最初的问题可以解决为：

def bd_from_ac(a, c):
    return a*neighbor_min(a), neighbor_min(c)

作为绩效基准：

In [2]: m, n = 1000, 1000

In [3]: a = np.random.rand(m, n)

In [4]: c = np.random.rand(m, n)

In [5]: %timeit bd_from_ac(a, c)
1 loops, best of 3: 123 ms per loop

Answer 2

查找a[min_coords]是一个滚动窗口操作。我们在this post中概述了几个聪明的解决方案。您将希望创建c[min_coords]数组成为您选择的任何解决方案的副作用。在

我希望这有帮助。我可以稍后有时间发布一些示例代码。在

Answer 3

我有兴趣帮助你，我相信在你的问题范围之外可能有更好的解决方案，但是为了把我自己的时间花在编写代码上，我必须得到你的一些反馈，因为我不能百分之百地确定我理解你需要什么。在

需要考虑的一件事是：如果您是一名C#开发人员，那么C#的“强力”实现可能会比聪明的Numpy实现更好，因此您至少可以考虑测试用C#实现的相当简单的操作。Geotiff（我想您正在阅读）有一个相对友好的规范，我想可能有.NET Geotiff库。在

但是假设你想尝试一下Numpy（我相信你应该这样做），让我们看看你想要达到的目标：

1. 如果要在数组a和c的每个元素中运行min_coords(array)，那么可以考虑使用numpy.dstack()和numpy.roll()来“堆叠”同一数组的九个副本，每个副本滚动一定的偏移量。然后，应用numpy.argmin(stacked_array, axis=2)，得到一个包含0到8之间的值的数组，其中每个值都映射到一个包含偏移索引的元组。在

然后，使用这个原理，min_coords()函数将被矢量化，一次在整个数组中操作，并返回一个数组，该数组给出一个偏移量，该偏移量是包含偏移量的查找表的索引。在

如果你有兴趣详细说明这一点，请留下评论。在

希望这有帮助！在