python - 如何在分配期间以矢量化方式移动输出变量的索引(以numpy为单位)
问题描述
动机:假设我有一个 RGB 图像 J,我想对 J 的像素应用一个变换 T(如旋转)。我将创建一个新的黑色图像 K,它的像素与 J 的关系为 K[x,y] =J[T[x,y]]。现在的问题是 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 循环编写的代码一样。
解决方案
天真的版本
据我了解您的问题:您有一个输入图像,您转换其像素位置,并希望将结果放入一个可以容纳它的更大数组中。这是我的做法:
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 值设置在结果数组的新位置。
如您所见,由于非整数变换后的坐标被四舍五入,因此图像中有很多黑色像素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:地图坐标
直接按照您的想法的方法可以利用逆变换scipy.ndimage.map_coordinates
执行插值。这应该比之前的尝试具有更好的性能,因为可以利用输入数据在网格上定义的事实。事实证明,它确实使用了更少的内存和更少的 CPU:griddata
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)
# 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
我主要关心图像质量,因此上述所有解决方案都以一种或另一种方式使用插值。正如您在评论中解释的那样,这不是您最关心的问题。如果是这种情况,我们可以使用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()
结果虽然充满了单像素的不准确,但看起来已经足够好了:
推荐阅读
- php - 如何在视图 laravel 上显示用户个人资料图片
- c# - Asp.Net Core 3.1 Cookie 身份验证循环到登录页面
- javascript - 尝试检查 firebase firestore 文档中的某个字段,但在其背后的逻辑上苦苦挣扎
- r - 有没有办法可以绘制这个密度图
- django - Django Rest Framework 中的购物车和购物车项目关系
- wordpress - WP中的链接没有转到正确的页面
- apache-kafka - Kafka 节点宕机时 Kafka Connect 挂起
- asp.net-mvc - 如何使用自定义数据库实现角色授权?
- c# - GroupPrincipal 是否应提供所有 AD 文件?
- promise - 抛出错误时如何在循环中重试获取