从python中的(2，MN)矩阵中读取(M，N)图像的值

Question

假设我有一个(M，N(图像J(形状为(M，N，3((。我有一个(2，MN(矩阵K，像这样：

0 0 0 ... 0 1 1 1 ... 1 ............. M M M ... M
0 1 2 ... N 1 2 3 ... N ............. 1 2 3 ... N

现在我把上面的矩阵乘以一个2乘2的矩阵，得到一个与K大小相同的新矩阵T。

现在我想创建一个新图像，其中这个新图像中的像素(r，s(等于旧图像的像素(r，G，B(值，旧图像位于T的r*N+M列中。

如果可能的话，我想以矢量化的方式来做这件事。我不想使用for循环(我已经知道如何使用for循环来实现这一点，但速度非常慢(。事实上，我对这个问题感兴趣是因为我想以矢量化的方式应用单应性变换。

感谢您的帮助。这是一个让我想要的东西变得清晰的古怪版本：

for r in range(0,M):
for s in range(0,N):
x, y = T[:,r*N+s]
new_image[r,s] = J[x,y]

Answer 1

由于内部的索引布局很好，您几乎可以直接使用索引数组

import numpy as np
# set up dummy input
M,N = 300,400
J = np.random.rand(M, N, 3)
T = np.array([np.random.randint(0, M, M*N), np.random.randint(0, N, M*N)])
# original    
new_image = np.empty_like(J)
for r in range(0,M):
for s in range(0,N):
x, y = T[:,r*N+s]
new_image[r,s] = J[x,y]
# vectorized new
new_image_vect = J[tuple(T)].reshape(J.shape)

检查：

>>> np.array_equal(new_image, new_image_vect)
True

上面的工作方式并不完全是琐碎的，因为高级索引是一件变化无常的事情。我上面写的相当于

J[(T[0,...], T[1,...])].reshape(J.shape) -> J[T[0,...], T[1,...]].reshape(J.shape)

现在第一部分更清楚了：取T第一行中的每个元素，并将其用作J的第一索引，然后取T第二行中的相应元素，并使用其作为J的相应第二索引。这部分在某种程度上涵盖了循环版本中的J[x,y]。

但是，由于我们本质上是用长度为M*N的1d个数组对数组进行索引，因此所得数组的形状也将沿其第一维度(以及大小为3的尾部维度(具有形状M*N。因此，我们需要将结果reshape返回到(M,N,3)。

所有这些仅直接起作用，因为T中的索引是根据C连续存储器顺序存储的。如果不是这样的话，我们将不得不来回transpose我们的数组，以便生成具有正确布局的结果数组。