55869


2つの長さが同じ3次元実数値ベクトル（空間上の位置）xとyが与えられた時、  
xからy（あるいはyからx）へ回転する行列Rを求めたいです。  
numpyやscipyを使って、どのようなコードでこれを実現できるでしょうか？


CC BY-SA 4.0
2019-06-17T02:22:21.547
55865
2019-06-17T10:41:49.117
2019-06-17T04:22:29.777
19110
12874
question
python
2つの3次元実数値ベクトルxとyが与えられた時、xからyへ回転する行列を求めたい


まず、Mathematics Stack Exchange
上の[同様の質問](https://math.stackexchange.com/q/180418/371701)に対する [Nico
Schlömer](https://math.stackexchange.com/users/36678/nico-schl%C3%B6mer)
さんの[回答](https://math.stackexchange.com/a/2672702/371701)を和訳します：

>
> [ロドリゲスの回転公式](https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula)を使うと、ベクトル
> `a` を回転軸 `k` まわりに角 `θ` だけ回転したときのベクトルが得られます。ここで単位ベクトル `a` を `b` に回転するには `a` を
> `(a + b) / 2` まわりに `π` だけ回転すれば良いことを考えると、この公式から以下の美しい式を得ます：
>
> [![R = 2 * \(a+b\)\(a+b\)^T / \(\(a+b\)^T \(a+b\)\) -
> I](https://i.stack.imgur.com/js4gK.png)](https://i.stack.imgur.com/js4gK.png)

この式は `a`, `b` が単位ベクトルのときに限らず、`a` と `b` の長さが等しく `a ≠ -b`
であれば使えます。このことから、上の式をそのまま実装してみました。

[code]

    import numpy as np
    
    def rotation_matrix_3d(frm, to):
        if np.array_equal(frm, to):
            return np.identity(3)
        if np.array_equal(frm, -to):
            return -np.identity(3)
    
        s = frm + to
        return 2.0 * np.outer(s, s) / np.dot(s, s) - np.identity(3)
    
[/code]

以下は動作テストの様子です。

[code]

    def generate(r):
        theta = np.random.uniform(0, 2.0 * np.pi)
        cosphi = np.random.uniform(-1.0, 1.0)
        sinphi = np.sqrt(1.0 - cosphi * cosphi)
        v = r * np.array([np.cos(theta) * sinphi, np.sin(theta) * sinphi, cosphi])
        return v
    
    if __name__ == '__main__':
        r = np.random.random()
        a = generate(r)
        b = generate(r)
        r = rotation_matrix_3d(a, b)
    
        # test
        print(" R =", r)
        print()
        print(" a =", a)
        print("Ra =", np.dot(r, a))
        print(" b =", b)
    
[/code]

[code]

    $ python test.py
     R = [[-0.33930792  0.7584119   0.55649036]
     [ 0.7584119  -0.129415    0.6387982 ]
     [ 0.55649036  0.6387982  -0.53127708]]
    
     a = [0.83254887 0.9255484  0.65935197]
    Ra = [0.78637951 0.93282798 0.70424549]
     b = [0.78637951 0.93282798 0.70424549]
    
[/code]


CC BY-SA 4.0
2019-06-17T04:19:23.050
55869
2019-06-17T10:41:49.117
2019-06-17T10:41:49.117
19110
19110
55865
answer