如何在分配期间以矢量化方式移动输出变量的索引（以numpy为单位）

Question

动机：假设我有一个 RGB 图像 J，我想对 J 的像素应用变换 T（如旋转）。我将创建一个新的黑色图像 K，其像素为通过 K[x,y]=J[T[x,y]] 与 J 相关。现在的问题是 T[x,y] 必须在 J 内部，如果我想完全捕获 J 的转换图像，我可能必须处理一些 x 或 y 的负值或大于大小的值J 的。所以，首先我必须确定 K 的大小，然后将 K 的像素移动一个适当的向量以避免负值。

现在，假设我已经确定了合适的平移向量。我想做一个坐标转换，将（x，y）发送到（x+a，y+k）。

目标：使用for循环，我想做的事情如下：

for i in range(0,J.shape[0]):
    for j in range(0, J.shape[1]):
        K[i+a,j+b] = J[T[i,j]]

如何以矢量化方式执行此操作？任何帮助表示赞赏。

---

编辑：

img = face() # dummy RGB data

i,j = np.mgrid[:img.shape[0], :img.shape[1]] # 2d arrays each
i_min, i_max, j_min, j_max = func(*) # assume that these values have been found
i = i + i_min
j = j + j_min
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
inew,jnew = np.linalg.inv(T) @ [i.ravel(), j.ravel()] # 1d arrays each

inew = np.floor(inew).astype(int)
jnew = np.floor(jnew).astype(int)

out = np.zeros((i_max - i_min, j_max - j_min, 3), dtype=img.dtype)

for i in inew:
    for j in jnew:
        out[i-i_min,j-j_min, :] = img[i,j,:]

现在我想取消在数组中移动 i_min 和 j_min 的效果，就像我使用 for 循环编写的代码一样。

Answer 1

朴素版

据我了解您的问题：您有一个输入图像，您转换其像素位置，并希望将结果放入可以容纳它的更大数组中。以下是我的做法：

import numpy as np
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data
img = face() # dummy RGB data

# transform pixels by 45 degrees
i,j = np.mgrid[:img.shape[0], :img.shape[1]] # 2d arrays each
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
inew,jnew = T @ [i.ravel(), j.ravel()] # 1d arrays each

# new coordinates now range into negatives, shift back into positives
# and the non-integer pixel indices will be normalized with floor
inew = np.floor(inew - inew.min()).astype(int)
jnew = np.floor(jnew - jnew.min()).astype(int)

# now the new coordinates are all non-negative, this defines the size of the output
out = np.zeros((inew.max() + 1, jnew.max() + 1, 3), dtype=img.dtype)

# fill the necessary indices of out with pixels from img
# reshape the indices to 2d for matching broadcast
inew = inew.reshape(img.shape[:-1])
jnew = jnew.reshape(img.shape[:-1])
out[inew, jnew, :] = img
# OR, alternative with 1d index arrays:
#out[inew, jnew, :] = img.reshape(-1, 3)

# check what we've done
plt.imshow(out)
plt.show()

代码的要点是将旋转的像素坐标移回正数（这对应于您的[i+a, j+b] 移位），分配一个新的零数组以适应所有新索引，并且索引仅适用于右侧！这与您的代码不匹配，但我相信这是您真正想要做的：对于原始（未索引）图像中的每个像素，我们将其 RGB 值设置在 new 位置结果数组。

如您所见，图像中有很多黑色像素，这是因为非整数转换后的坐标是用floor 四舍五入的。这不好，所以如果我们追求这条路径，我们应该执行 2d 插值以消除这些伪影。请注意，这需要相当多的内存和 CPU 时间：

import numpy as np
import scipy.interpolate as interp
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data
img = face() # dummy RGB data

# transform pixels by 45 degrees
i,j = np.mgrid[:img.shape[0], :img.shape[1]] # 2d arrays each
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
inew,jnew = T @ [i.ravel(), j.ravel()] # 1d arrays each

# new coordinates now range into negatives, shift back into positives
# keep them non-integer for interpolation later
inew -= inew.min()
jnew -= jnew.min()
# (inew, jnew, img) contain the data from which the output should be interpolated


# now the new coordinates are all non-negative, this defines the size of the output
out = np.zeros((int(round(inew.max())) + 1, int(round(jnew.max())) + 1, 3), dtype=img.dtype)
i_interp,j_interp = np.mgrid[:out.shape[0], :out.shape[1]]

# interpolate for each channel
for channel in range(3):
    out[..., channel] = interp.griddata(np.array([inew.ravel(), jnew.ravel()]).T, img[..., channel].ravel(), (i_interp, j_interp), fill_value=0)

# check what we've done
plt.imshow(out)
plt.show()

至少结果看起来好多了：

scipy.ndimage: map_coordinates

直接按照您的想法的方法可以利用scipy.ndimage.map_coordinates 使用逆 变换执行插值。这应该比之前使用griddata 的尝试具有更好的性能，因为map_coordinates 可以利用输入数据在网格上定义的事实。事实证明，它确实使用更少的内存和更少的 CPU：

import numpy as np
import scipy.ndimage as ndi
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data

img = face() # dummy RGB data
n,m = img.shape[:-1]

# transform pixels by 45 degrees
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)

# find out the extent of the transformed pixels from the four corners
inew_tmp,jnew_tmp = T @ [[0, 0, n-1, n-1], [0, m-1, 0, m-1]] # 1d arrays each
imin,imax,jmin,jmax = inew_tmp.min(),inew_tmp.max(),jnew_tmp.min(),jnew_tmp.max()
imin,imax,jmin,jmax = (int(round(val)) for val in (imin,imax,jmin,jmax))

# so the pixels of the original map inside [imin, imax] x [jmin, jmax]
# we need an image of size (imax - imin + 1, jmax - jmin + 1) to house this
out = np.zeros((imax - imin + 1, jmax - jmin + 1, 3), dtype=img.dtype)
# indices have to be shifted by [imin, imax]

# compute the corresponding (non-integer) coordinates on the domain for interpolation
inew,jnew = np.mgrid[:out.shape[0], :out.shape[1]]
i_back,j_back = np.linalg.inv(T) @ [inew.ravel() + imin, jnew.ravel() + jmin]

# perform 2d interpolation for each colour channel separately
for channel in range(3):
    out[inew, jnew, channel] = ndi.map_coordinates(img[..., channel], [i_back, j_back]).reshape(inew.shape)

# check what we've done
plt.imshow(out)
plt.show()

结果还是不错的：

scipy.ndimage: 几何变换

最后，我意识到我们可以再上一层，直接使用scipy.ndimage.geometric_transform。对于旋转的浣熊案例，这似乎比使用map_coordinates 的手动版本要慢，但会导致代码更简洁：

import numpy as np
import scipy.ndimage as ndi
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data

img = face() # dummy RGB data
n,m = img.shape[:-1]

# transform pixels by 45 degrees
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)
Tinv = np.linalg.inv(T)

# find out the extent of the transformed pixels from the four corners
inew_tmp,jnew_tmp = T @ [[0, 0, n-1, n-1], [0, m-1, 0, m-1]] # 1d arrays each
imin,imax,jmin,jmax = inew_tmp.min(),inew_tmp.max(),jnew_tmp.min(),jnew_tmp.max()
imin,imax,jmin,jmax = (int(round(val)) for val in (imin,imax,jmin,jmax))

# so the pixels of the original map inside [imin, imax] x [jmin, jmax]
# we need an image of size (imax - imin + 1, jmax - jmin + 1) to house this

def transform_func(output_coords):
    """Inverse transform output coordinates back into input coordinates"""
    inew,jnew,channel = output_coords
    i,j = Tinv @ [inew + imin, jnew + jmin]
    return i,j,channel

out = ndi.geometric_transform(img, transform_func, output_shape = (imax - imin + 1, jmax - jmin + 1, 3))

# check what we've done
plt.imshow(out)
plt.show()

结果：

最终修复：仅 numpy

我主要关心图像质量，因此上述所有解决方案都以一种或另一种方式使用插值。正如您在 cmets 中解释的那样，这不是您最关心的问题。如果是这种情况，我们可以使用map_coordinates 修改版本并自己计算近似（舍入整数）索引并执行矢量化赋值：

import numpy as np
import matplotlib.pyplot as plt # for plotting the result
from scipy.misc import face # for dummy data

img = face() # dummy RGB data
n,m = img.shape[:-1]

# transform pixels by 45 degrees
T = np.array([[1, -1],[1, 1]])/np.sqrt(2)

# find out the extent of the transformed pixels from the four corners
inew_tmp,jnew_tmp = T @ [[0, 0, n-1, n-1], [0, m-1, 0, m-1]] # 1d arrays each
imin,imax,jmin,jmax = inew_tmp.min(),inew_tmp.max(),jnew_tmp.min(),jnew_tmp.max()
imin,imax,jmin,jmax = (int(round(val)) for val in (imin,imax,jmin,jmax))

# so the pixels of the original map inside [imin, imax] x [jmin, jmax]
# we need an image of size (imax - imin + 1, jmax - jmin + 1) to house this
out = np.zeros((imax - imin + 1, jmax - jmin + 1, 3), dtype=img.dtype)

# compute the corresponding coordinates on the domain for matching
inew,jnew = np.mgrid[:out.shape[0], :out.shape[1]]
inew = inew.ravel() # 1d array, indices of output array
jnew = jnew.ravel() # 1d array, indices of output array
i_back,j_back = np.linalg.inv(T) @ [inew + imin, jnew + jmin]

# create a mask to grab only those rounded (i_back,j_back) indices which make sense
i_back = i_back.round().astype(int)
j_back = j_back.round().astype(int)
inds = (0 <= i_back) & (i_back < n) & (0 <= j_back) & (j_back < m)
# (i_back[inds], j_back[inds]) maps to (inew[inds], jnew[inds])
# the rest stays black

out[inew[inds], jnew[inds], :] = img[i_back[inds], j_back[inds], :]

# check what we've done
plt.imshow(out)
plt.show()

结果虽然充满了单像素的不准确，但看起来已经足够好了：

Answer 2

你可以使用地图功能

for i in range(0,J.shape[0]):
    for j in range(0, J.shape[1]):
        K[i+a,j+b] = J[T[i,j]]

例如，您可以生成矩阵的所有索引元组

indexes = [ (i,j) for i in range(J.shape[0]) for j in range(J.shape[1]) ]

然后使用 lambda 函数应用地图

f = lambda coords:  J[T[coords[0],coords[1]]]
resp = list(map(f, indexes))

此时 resp 包含 f 对索引的所有应用的列表。现在你必须把它重塑成好的形状。对于K

所以这里有两种可能，你可以将范围列表设为 K 的大小，然后在 lambda 函数中需要时返回零

旧答案...

这里的问题是你必须事先知道输出图像的大小。 所以有两种可能性，要么计算它，要么假设它不会大于某个估计值。

因此，如果您计算它，要走的路取决于您要应用的转换。 例如，转置表示 X 轴和 Y 轴长度的交换。 对于旋转，结果的大小取决于形状和角度。

所以

如果你想保持非常非常简单 但不一定对内存友好。假设您的转换不会输出大于 X 和 Y 长度最大值的三倍的图像。

这样做，您可以轻松处理偏移量

如果您的图像是NxM 和N > M，则转换的画布将为3*Nx3*N

现在假设输出图像将在此画布中居中。 在这种情况下，您必须计算您在问题中描述的 a 和 b 偏移量

沿垂直轴的变换图像的中心应与原始图像的中心匹配。

if i=N/2 then a+i=3*N/2 这意味着a=N

同样适用于水平轴，在这种情况下

if j=M/2 then b+j=3*N/2 这意味着b=(3*N - M)/2

我希望清楚