我有 4 个 2D numpy 数组，称为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];
wheremin_of_neighbors_coords是一个函数，给定数组元素的坐标，返回具有较低值的“邻居”元素的坐标。即，考虑数组：
1, 2, 5
3, 7, 2
2, 3, 6
min_of_neighbors_coords(1, 1)将引用值为 的中心元素7，并返回元组(0, 0)：数字的坐标1。
我设法使用 for 循环（每个元素的元素）来做到这一点，但算法非常慢，我正在寻找一种改进它的方法，避免循环并要求计算 numpy。
是否可以？
编辑我将原始答案保留在底部。正如保罗在评论中指出的那样，原始答案并没有真正回答 OP 的问题，并且可以使用 ndimage 过滤器更容易地实现。以下更繁琐的功能应该做正确的事情。它需要两个数组a和c，并返回 的窗口最小值和中窗口最小值位置处的a值：ca
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)数组，可以解包成两个数组：
>>> a = np.random.randint(100, size=(5,5))
>>> c = np.random.randint(100, size=(5,5))
>>> a
array([[42, 54, 18, 88, 26],
[80, 65, 83, 31, 4],
[51, 52, 18, 88, 52],
[ 1, 70, 5, 0, 89],
[47, 34, 27, 67, 68]])
>>> c
array([[94, 94, 29, 6, 76],
[81, 47, 67, 21, 26],
[44, 92, 20, 32, 90],
[81, 25, 32, 68, 25],
[49, 43, 71, 79, 77]])
>>> neighbor_min(a, c)
array([[[42, 18, 18, 4, 4],
[42, 18, 18, 4, 4],
[ 1, 1, 0, 0, 0],
[ 1, 1, 0, 0, 0],
[ 1, 1, 0, 0, 0]],
[[94, 29, 29, 26, 26],
[94, 29, 29, 26, 26],
[81, 81, 68, 68, 68],
[81, 81, 68, 68, 68],
[81, 81, 68, 68, 68]]])
然后可以将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)数据的副本。幸运的是，您可以分两步计算最小值，首先沿一个轴开窗和最小化，然后沿另一个轴，并获得相同的结果。所以最多你会有输入大小的中间存储。如果需要，您甚至可以重用输出数组作为中间存储并避免内存分配，但这留作练习......
下面的函数就是这样做的。它有点冗长，因为它不仅要处理中心区域，还要处理四个边缘和四个角的特殊情况。除此之外，它是一个非常紧凑的实现：
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
查找a[min_coords]是滚动窗口操作。我们在这篇文章中概述了几个巧妙的解决方案。您需要使 c[min_coords] 数组的创建成为您选择的任何解决方案的副作用。
我希望这有帮助。稍后我有时间可以发布一些示例代码。
我有兴趣帮助你，我相信在你的问题范围之外可能有更好的解决方案，但是为了把我自己的时间用于编写代码，我必须得到你的一些反馈，因为我不是 100% 确定我了解你需要什么。
需要考虑的一件事：如果您是 C# 开发人员，也许 C# 的“蛮力”实现可以胜过 Numpy 的巧妙实现，因此您可以考虑至少测试您在 C# 中实现的相当简单的操作。Geotiff（我想你正在阅读）有一个相对友好的规范，我想可能有 .NET GeoTiff 库。
但是假设你想尝试一下 Numpy（我相信你应该），让我们看看你想要实现的目标：
min_coords(array)在数组的每个元素中运行aand c，则可以考虑使用numpy.dstack()and将同一数组的九个副本“堆叠”，每个副本滚动一些偏移量numpy.roll()。然后，您应用numpy.argmin(stacked_array, axis=2)并获得一个包含 0 到 8 之间值的数组，其中每个值都映射到包含偏移索引的元组。然后，使用此原理，您的min_coords()函数将被矢量化，一次在整个数组中操作，并返回一个数组，该数组为您提供一个偏移量，该偏移量将是包含偏移量的查找表的索引。
如果您有兴趣详细说明这一点，请发表评论。
希望这可以帮助！