Numpyで値のペアのための行列を生成する

Question

私は回転行列を使って各ベクトルを変換したい3D配列（ベクトルの2次元配列）を持っています。回転は、colsとrowsと呼ばれるラジアン角度値の2つの別々の2D配列にあります。

私は既にNumPyにPythonループを使わずに私の角度を計算させることができました。今私は、NumPyに回転行列を生成させる方法を探しています。うまくいけば、パフォーマンスが大幅に向上します。

size = img.shape[:2] 

# Create an array that assigns each pixel the percentage of 
# the correction (value between -1 and 1, distributed linearly). 
cols = np.array([np.arange(size[1]) for __ in range(size[0])]) /(size[1] - 1) * 2 - 1 
rows = np.array([np.arange(size[0]) for __ in range(size[1])]).T/(size[0] - 1) * 2 - 1 

# Atan distribution based on F-number and Sensor size. 
cols = np.arctan(sh * cols/(2 * f)) 
rows = np.arctan(sv * rows/(2 * f)) 

### This is the loop that I would like to remove and find a 
### clever way to make NumPy do the same operation natively. 
for i in range(size[0]): 
    for j in range(size[1]): 
    ah = cols[i,j] 
    av = rows[i,j] 

    # Y-rotation. 
    mat = np.matrix([ 
     [ np.cos(ah), 0, np.sin(ah)], 
     [0, 1, 0], 
     [-np.sin(ah), 0, np.cos(ah)] 
    ]) 

    # X-rotation. 
    mat *= np.matrix([ 
     [1, 0, 0], 
     [0, np.cos(av), -np.sin(av)], 
     [0, np.sin(av), np.cos(av)] 
    ]) 

    img[i,j] = img[i,j] * mat 

return img 

NumPy操作でループを書き直す方法はありますか？

Answer 1

（のは(a, b, 3)ことimgの形状を想定してみましょう。）

まず、colsとrowsは（あなたの代わりにcols[i,j]のcols[j]を書くことができる）完全(a, b)に拡張する必要はありません。そして、彼らは簡単にnp.linspaceを使用して生成することができます。

cols = np.linspace(-1, 1, size[1]) # shape: (b,) 
rows = np.linspace(-1, 1, size[0]) # shape: (a,) 

cols = np.arctan(sh * cols/(2*f)) 
rows = np.arctan(sv * rows/(2*f)) 

その後、我々は、行列の要素を事前計算得ます。

# shape: (b,) 
cos_ah = np.cos(cols) 
sin_ah = np.sin(cols) 
zeros_ah = np.zeros_like(cols) 
ones_ah = np.ones_like(cols) 

# shape: (a,) 
cos_av = np.cos(rows) 
sin_av = np.sin(rows) 
zeros_av = np.zeros_like(rows) 
ones_av = np.ones_like(rows) 

そして、回転行列を構築する：

# shape: (3, 3, b) 
y_mat = np.array([ 
    [cos_ah, zeros_ah, sin_ah], 
    [zeros_ah, ones_ah, zeros_ah], 
    [-sin_ah, zeros_ah, cos_ah], 
]) 

# shape: (3, 3, a) 
x_mat = np.array([ 
    [ones_av, zeros_av, zeros_av], 
    [zeros_av, cos_av, -sin_av], 
    [zeros_av, sin_av, cos_av], 
]) 

今見てみましょう。私たちは、マトリックス乗算アウト展開する場合、

for i in range(size[0]): 
    for j in range(size[1]): 
     img[i, j, :] = img[i, j, :] @ y_mat[:, :, j] @ x_mat[:, :, i] 

のか：：私たちはループを持っている場合は、私たちが書くでしょう。これはnp.einsumを使用してきれいに処理することができ

を（I、の点に注意してください。 j,k,m,nは上記の式と全く同じです。

img = np.einsum('ijk,kmj,mni->ijn', img, y_mat, x_mat) 

---

要約する：回転行列の

size = img.shape[:2] 

cols = np.linspace(-1, 1, size[1]) # shape: (b,) 
rows = np.linspace(-1, 1, size[0]) # shape: (a,) 

cols = np.arctan(sh * cols/(2*f)) 
rows = np.arctan(sv * rows/(2*f)) 

cos_ah = np.cos(cols) 
sin_ah = np.sin(cols) 
zeros_ah = np.zeros_like(cols) 
ones_ah = np.ones_like(cols) 

cos_av = np.cos(rows) 
sin_av = np.sin(rows) 
zeros_av = np.zeros_like(rows) 
ones_av = np.ones_like(rows) 

y_mat = np.array([ 
    [cos_ah, zeros_ah, sin_ah], 
    [zeros_ah, ones_ah, zeros_ah], 
    [-sin_ah, zeros_ah, cos_ah], 
]) 

x_mat = np.array([ 
    [ones_av, zeros_av, zeros_av], 
    [zeros_av, cos_av, -sin_av], 
    [zeros_av, sin_av, cos_av], 
]) 

return np.einsum('ijk,kmj,mni->ijn', img, y_mat, x_mat)