Numpy坐标矩阵

Question

我想得到一个坐标数组的矩阵。这与numpy.meshgrid不同。例如，对于2x2大小，我想要2x2x2输出

[[[0,0],[0,1]],
 [[1,0],[1,1]]]

作为一个核阵列。这可能看起来和读取更干净的2x2元组矩阵：

[[(0,0),(0,1)],
 [(1,0),(1,1)]]

（除了我认为你不能在numpy数组中有元组，这不是重点）

这个简单的例子可以通过切换numpy meshgrid输出的轴来完成（具体来说，移动要最后一个的第一个轴）：

np.array(np.meshgrid([0,1],[0,1])).transpose([1,2,0])

这可以很容易地推广到任意维度，除了meshgrid的行为不像我对2个以上输入的期望。具体来说，返回的矩阵具有沿轴以奇数顺序变化的坐标值：

In [627]: np.meshgrid([0,1],[0,1],[0,1])
Out[627]:
[array([[[0, 0],
        [1, 1]],

       [[0, 0],
        [1, 1]]]),
 array([[[0, 0],
        [0, 0]],

       [[1, 1],
        [1, 1]]]),
 array([[[0, 1],
        [0, 1]],

       [[0, 1],
        [0, 1]]])]

请注意，此输出的元素分别沿轴1、0和2变化。这将建立一个不正确的坐标矩阵；我需要输出沿轴0、1和2按顺序变化。所以我可以

In [642]: np.array(np.meshgrid([0,1],[0,1],[0,1])).swapaxes(1,2)
Out[642]:
array([[[[0, 0],
         [0, 0]],

        [[1, 1],
         [1, 1]]],


       [[[0, 0],
         [1, 1]],

        [[0, 0],
         [1, 1]]],


       [[[0, 1],
         [0, 1]],

        [[0, 1],
         [0, 1]]]])

但这已经开始变得非常麻烦，我不知道我是否可以指望在更高维度的网格输出这个顺序。numpy.mgrid给出了正确的顺序，但似乎不允许任意值，这是我需要的。所以这可以归结为两个问题：

1）是否有一种更简洁的方法，可能是缺少numpy中的某个函数，它将生成一个坐标向量矩阵，如前所述？ 2） 这个奇怪的顺序真的是我们对meshgrid的期望吗？这一点有我可以信赖的规范吗？

[编辑]继Jaime的解决方案之后，这里有一个更通用的函数来为任何感兴趣的人更显式地构建它：[编辑2，修复了一个bug，可能是另一个bug，现在不能在这个上面花费太多时间，这真的需要一个更常见的函数…]

def build_coords(*vecs):
    coords = numpy.empty(map(len,vecs)+[len(vecs)])
    for ii in xrange(len(vecs)):
        s = np.hstack((len(vecs[ii]), np.ones(len(vecs)-ii-1)))
        v = vecs[ii].reshape(s)
        coords[...,ii] = v
    return coords

Answer 1

给定1D坐标：

rows = np.arange(2)
cols = np.arange(3)

我希望这能起到作用：

np.dstack((rows[:, None, None], cols[:, None]))

但显然dstack等需要完全匹配的维度，它们不会广播它们，我认为这是一个遗憾。

所以这个选择有点长，但是显式比隐式好，你可以把它都包装成一个小函数：

>>> coords = np.empty((len(rows), len(cols), 2), dtype=np.intp)
>>> coords[..., 0] = rows[:, None]
>>> coords[..., 1] = cols

>>> coords
array([[[0, 0],
        [0, 1],
        [0, 2]],

       [[1, 0],
        [1, 1],
        [1, 2]]])

Answer 2

numpy函数^{}也可以用于此效果，它的功能也可以从名称中清除。

>>> import numpy as np
>>> np.indices((2,3))
array([[[0, 0, 0],
        [1, 1, 1]],

       [[0, 1, 2],
        [0, 1, 2]]])

可以认为是y坐标的2×3矩阵和x坐标的2×3矩阵（y,x = np.indices((2,3))）。它可以通过调换轴重铸成詹姆提议的形式：

>>> np.indices((2,3)).transpose((1,2,0))

它在功能上等同于使用indexing='ij'的the ^{} solution，但不要求您提供坐标数组，当您有许多维度时，这是一个好处。

>>> def f1(shape):
...     return np.array(np.meshgrid(*(np.arange(s) for s in shape), indexing='ij'))
...
>>> shape = (200, 31, 15, 4)
>>> np.all(f1(shape) == np.indices(shape))
True

从时间上看，当考虑到生成meshgrid操作的一维数组所需的时间时，这些解决方案是相似的，但是meshgrid返回一个（数组的）列表，而不是像indices这样的nd数组。通过添加对np.array的额外调用（如上面f1中所做的），与meshgrid相比，indices具有明显的优势：

In [14]: %timeit f1(shape)
100 loops, best of 3: 14 ms per loop

In [15]: %timeit np.indices(shape)
100 loops, best of 3: 5.77 ms per loop

如果没有对array的额外调用：

In [16]: def f2(shape):
    return np.meshgrid(*(np.arange(s) for s in shape), indexing='ij')
   .....: 

In [17]: %timeit f2(shape)
100 loops, best of 3: 5.78 ms per loop

不过，在解释计时时一定要小心。这可能不会成为你解决任何问题的瓶颈。

在任何情况下，meshgrid可以做比indices更多的事情，例如生成更通用的rectilinear grid而不是笛卡尔网格，因此在适当的时候使用它们。在本例中，我将使用更具描述性的命名方式indices。

Answer 3

尝试np.meshgrid([0, 1], [0, 1], [0, 1], indexing="ij")。与非默认的indexing="ij"相比，meshgrid文档实际上非常明确地说明了默认的indexing="xy"是如何生成有趣的轴顺序的，因此您可以查看更多详细信息。（他们对为什么会这样工作还不太清楚，唉…）

Answer 4

numpy函数^{}也可以用于此效果，它的功能也可以从名称中清除。

>>> import numpy as np
>>> np.indices((2,3))
array([[[0, 0, 0],
        [1, 1, 1]],

       [[0, 1, 2],
        [0, 1, 2]]])

可以认为是y坐标的2×3矩阵和x坐标的2×3矩阵（y,x = np.indices((2,3))）。它可以通过调换轴重铸成詹姆提议的形式：

>>> np.indices((2,3)).transpose((1,2,0))

它在功能上等同于使用indexing='ij'的the ^{} solution，但不要求您提供坐标数组，当您具有许多维度时，这将是一个好处。

>>> def f1(shape):
...     return np.array(np.meshgrid(*(np.arange(s) for s in shape), indexing='ij'))
...
>>> shape = (200, 31, 15, 4)
>>> np.all(f1(shape) == np.indices(shape))
True

从时间上看，当考虑到生成meshgrid操作的一维数组所需的时间时，这些解决方案是相似的，但是meshgrid返回一个（数组的）列表，而不是像indices这样的nd数组。通过添加对np.array的额外调用（如上面f1中所做的），与meshgrid相比，indices具有明显的优势：

In [14]: %timeit f1(shape)
100 loops, best of 3: 14 ms per loop

In [15]: %timeit np.indices(shape)
100 loops, best of 3: 5.77 ms per loop

如果没有对array的额外调用：

In [16]: def f2(shape):
    return np.meshgrid(*(np.arange(s) for s in shape), indexing='ij')
   .....: 

In [17]: %timeit f2(shape)
100 loops, best of 3: 5.78 ms per loop

不过，在解释计时时一定要小心。这可能不会成为你解决任何问题的瓶颈。

在任何情况下，meshgrid可以做比indices更多的事情，例如生成更通用的rectilinear grid而不是笛卡尔网格，因此在适当的时候使用它们。在本例中，我将使用更具描述性的命名方式indices。