かみのメモ

コンピュータビジョン・プログラムな話題中心の勉強メモ

Python scipyのRotationモジュールで三次元回転を扱う

このブログでは以前、回転ベクトル, 回転行列, クォータニオン四元数), オイラー角についてまとめる記事を書きました。

その後、実際のコーディングで利用できるライブラリはないものかと探していたところ、scipyにぴったりのモジュールを見つけました。 今回はこのモジュールについて軽くまとめてみます。

※ 2020/02/03 : コメントにて「*_euler()では軸の指定を大文字/小文字のどちらで行うかによってintrinsic/extrinsicが切り替わるよ」とご指摘いただき、該当箇所を修正しました。教えていただきありがとうございます!

スポンサーリンク

紹介するモジュール

今回紹介するのはscipyのRotationモジュールです。

公式ドキュメントがしっかりしてるので、英語アレルギーのない人はこの記事よりドキュメントを読むほうがいいかもしれません(笑)

使い方

このモジュールでできることは

  1. 回転ベクトル, 回転行列, クォータニオン, オイラー角の相互変換
  2. 三次元点への回転の適用
  3. 回転の合成

の3つです。

注意点として、このモジュールのAPIはすべて右手座標系向けに設計されています。 もし左手系の回転表現の計算に使うなら、回転表現と座標をすべて右手系に変換してからこのモジュールで計算しその結果を左手系に変換し直す必要があります。 右手系と左手系の変換方法については、以前の記事を参照してください。

まずは、回転ベクトル, 回転行列, クォータニオン, オイラー角のいずれかからRotationインスタンスを作成します。 これにはRotationクラスのfrom_*()メソッドを使います。

from scipy.spatial.transform import Rotation
import numpy as np

# 回転ベクトルから
rotvec = np.array([0, np.pi/2, np.pi/3])
rot = Rotation.from_rotvec(rotvec)

# 回転行列から
# - 与えられた行列に近い回転表現を最適化で計算しているので
# - 浮動小数点の誤差などで要素値が多少ブレてもエラーを出さずに動作するらしい
mat = np.array([[-0.31178007, -0.52705078,  0.79057616],
                [ 0.52705078,  0.59637536,  0.60543695],
                [-0.79057616,  0.60543695,  0.09184457]])
rot = Rotation.from_dcm(mat)

# クォータニオンから
# - x, y, z, w の順
quat = np.array([0.        , 0.67385289, 0.44923526, 0.58660887])
rot = Rotation.from_quat(quat)

# オイラー角(xyx intrinsic Euler angles)から(ラジアン)
euler = np.array([ 0.5880026 ,  1.88786223, -0.5880026])
rot = Rotation.from_euler('XYX', euler)
# オイラー角(xyx intrinsic Euler angles)から(弧度)
euler = np.array([ 33.69006753, 108.16653826, -33.69006753])
rot = Rotation.from_euler('XYX', euler, degrees=True)

# オイラー角(xyx exstrinsic Euler angles)から(ラジアン)
euler = np.array([-0.5880026 ,  1.88786223,  0.5880026 ])
rot = Rotation.from_euler('xyx', euler)
# オイラー角(xyx exstrinsic Euler angles)から(弧度)
euler = np.array([-33.69006732, 108.16653808,  33.69006732])
rot = Rotation.from_euler('xyx', euler, degrees=True)

オイラー角のintrinsic/extrinsicの違いについても、以前の記事を参照いただければと思います。

Rotationインスタンス内部でどのように回転の情報を保持しているのかは不明ですが、共役クォータニオンを区別して保持しているので、クォータニオンを使っているのではないかという気がします。

Rotationクラスが持つas_*()メソッドを呼び出せば、他の回転表現に変換できます。

from scipy.spatial.transform import Rotation
import numpy as np

rot = Rotation.from_rotvec(np.array([0, np.pi/2, np.pi/3]))

# 回転ベクトルに変換
rot.as_rotvec()
# > array([0.        , 1.57079633, 1.04719755])

# 回転行列に変換
rot.as_dcm()
# > array([[-0.31178007, -0.52705078,  0.79057616],
# >        [ 0.52705078,  0.59637536,  0.60543695],
# >        [-0.79057616,  0.60543695,  0.09184457]])

# クォータニオンに変換
# - x, y, z, w の順
rot.as_quat()
# > array([-0.        , -0.67385289, -0.44923526, -0.58660887])

# オイラー角(xyx intrinsic Euler angles)に変換
# - ジンバルロック(無数のオイラー角表現が存在しうる)時は、
# - warningを出した上で3番目の要素が0であるオイラー角を返す
rot.as_euler('XYX')
# > array([ 0.5880026 ,  1.88786223, -0.5880026 ])

# オイラー角(xyx extrinsic Euler angles)に変換
rot.as_euler('xyx')
# > array([-0.5880026 ,  1.88786223,  0.5880026 ])

Rotationインスタンスを使って三次元点を原点中心に回転させることもできます。

from scipy.spatial.transform import Rotation
import numpy as np

rot = Rotation.from_rotvec(np.array([0, np.pi/2, np.pi/3]))

# 三次元座標点(1,2,3)を原点を中心に回転させる
rot.apply(np.array([1, 2, 3]))
# > array([1.00584687, 3.53611236, 0.69583145])

また、Rotationインスタンス同士で回転の重ねがけ(合成)もできます。

from scipy.spatial.transform import Rotation
import numpy as np

rot1 = Rotation.from_rotvec([0, np.pi/2, np.pi/3])
rot2 = Rotation.from_rotvec([0, np.pi, np.pi])

# 回転を重ねがけする(rot1 -> rot2の順に回転を適用するとき)
rot3 = rot2 * rot1

rot3.apply(np.array([1, 2, 3]))
# > array([1.66807229, 1.05228478, 3.17965903])

rot2.apply(rot1.apply(np.array([1, 2, 3])))
# > array([1.66807229, 1.05228478, 3.17965903])

また、Rotationはndarrayっぽく複数の回転情報をまとめて保持することもできます。

from scipy.spatial.transform import Rotation
import numpy as np

# zyx extrinsicオイラー角を使って3つの回転の情報を生成する
rot = R.from_euler('zyx', [
    [90, 0, 0],
    [0, 45, 0],
    [45, 60, 30]], degrees=True)

# 3つの回転の情報をクォータニオン形式で表示させる
rot.as_quat()
# > array([[0.        , 0.        , 0.70710678, 0.70710678],
# >        [0.        , 0.38268343, 0.        , 0.92387953],
# >        [0.39190384, 0.36042341, 0.43967974, 0.72331741]])

以上、簡単ではありますがscipyのRotationモジュールを使って三次元回転を表現・変換する方法についてまとめてみました。