PythonでQuaternionを使う ~ numpy-quaternion ~

はじめに

Pythonでクォータニオンを扱うライブラリはpyquaternionとnumpy-quaternionが世界でのトップ2のようですが，日本ではpyquaternionの参考ページを作った人が最初にいたからか，巷に溢れているPythonでのクォータニオン計算はpyquaternionばっか（しかない？）です．

しかし，numpy-quaternionのほうが計算コストが低そうです．
参考: https://www.theoj.org/joss-papers/joss.00787/10.21105.joss.00787.pdf

また，githubのスターの数を比べてもnumpy-quaternionのほうが多いです（2019/12/13時点でnumpy-quaternionが276，pyquaternionが138）．

numpyとの親和性も一括変換の扱いやすさなどでpyquaternionよりも良い気がします．使ってみませんか？

githubのリポジトリ: https://github.com/moble/quaternion
ドキュメント: https://quaternion.readthedocs.io/en/latest/
PyPI: https://pypi.org/project/numpy-quaternion/

インストール

pip install numpy-quaternion

で簡単に入ります．

また，

import quaternion

でライブラリがインポートできます．

クォータニオンの表記

quat = (w, x, y, z) = \left(\underbrace{\cos\frac{\theta}{2}}_{実数部},\ \underbrace{\lambda_x\sin\frac{\theta}{2},\ \lambda_y\sin\frac{\theta}{2},\ \lambda_z\sin\frac{\theta}{2}}_{虚数(ベクトル)部}\right)

使い方

以下，

import numpy as np
import quaternion

とします．

クォータニオンを作る

quat = np.quaternion(w, x, y, z)

numpyの属性にquaternionが追加されます．
例えば，

print(np.quaternion(1,0,0,0))
# -> quaternion(1, 0, 0, 0)

type(np.quaternion(1,0,0,0))
# -> quaternion.quaternion

となります．

クォータニオンの合成

q1 = np.quaternion(1, 2, 3, 4)
q2 = np.quaternion(4, 5, 6, 7)
print(q1*q2)
# -> quaternion(-52, 10, 24, 20)

掛け算が定義されています．一応他の四則演算も定義されていますが，和とか使うことありますかね…C++のEigenでは定義されていません．

np.quaternionのメソッド

一部のみ取り上げます．球面関数関係のものや，線形補間などは取り上げていません．
参考: https://quaternion.readthedocs.io/en/latest/_autosummary/quaternion.html

メンバ変数

メンバ変数	機能
w	実数部の要素
x	虚数部の最初の要素
y	虚数部の2番目の要素
z	虚数部の3番目の要素
components	(w,x,y,z)がnumpy.arrayで返ってくる
imag	虚数部(x,y,z)がnumpy.arrayで返ってくる
vec	虚数部(x,y,z)がnumpy.arrayで返ってくる
real	実数部(w)が返ってくる

メンバ関数

メンバ関数	機能
abs()	クォータニオン（ユークリッド距離）の絶対値
absolute()	クォータニオンの絶対値
angle()	回転角度
conj()	クォータニオンの複素共役を返す
conjugate()	クォータニオンの複素共役を返す
equal(quat)	引数の中身のクォータニオンと等しいか
exp	exponetialを返す（$e^q$）
inverse()	逆クォータニオンを返す
isfinite()	全ての要素が有限か
ininf()	1つでもinfの要素が存在するか
innan()	1つでもnanの要素が存在するか
log()	クォータニオンのログを返す
nonzero()	全ての要素が0か
norm()	クォータニオンのCayley norm（絶対値の二乗のルート）
normalized()	正規化したクォータニオンを返す
notequal(quat)	引数の中身のクォータニオンと等しくないか
sqrt()	クォータニオンのsquare-root($quat=q*q$を満たす$q$)を返す
square()	クォータニオンの2乗$(quat*quat)$を返す

quaternionのメソッド

この辺のメソッドは1次元のものではなく，多次元配列にも利用できるのがポイント．最後の次元の大きさがクォータニオンを要求するものなら4，3次元のものなら3次元，3x3なら最後の2次元が3x3にするなどはしないといけない．

メンバ関数	機能
quaternion.as_quat_array(a)	numpy.arrayをquaternionに変換．aの最後の次元のサイズは4でないといけない
quaternion.as_float_array(a)	numpy.quaternionをnumpy.arrayに変換．出力の次元は入力より1大きい．
quaternion.from_float_array(a)	as_quat_arrayと同じ
quaternion.as_rotation_matrix(q)	numpy.quaternionを3x3の回転行列に変換．
quaternion.from_rotation_matrix(rot, nonorthogonal=True)	3x3の回転行列をnumpy.quaternionに変換
quaternion.as_rotation_vector(q)	クォータニオンから回転軸を求める．出力の最後の次元の大きさは3．
quaternion.from_rotation_vector(rot)	サイズ3の回転軸からクォータニオンに変換する．
quaternion.as_euler_angles(q)	クォータニオンからオイラー角に変換．オイラー角の変換順などはドキュメント参照
quaternion.from_euler_angles(alpha_beta_gamma, beta=None, gamma=None)	オイラー角からクォータニオンに変換．オイラー角の変換順などはドキュメント参照
quaternion.rotate_vectors(R, v, axis=-1)	Rはクォータニオン，vはベクトル．クォータニオンに応じてベクトルを回転させる．
quaternion.allclose(a, b, rtol=8.881784197001252e-16, atol=0.0, equal_nan=False, verbose=False)	2つのクォータニオンを比較する
quaternion.integrate_angular_velocity(Omega, t0, t1, R0=None, tolerance=1e-12)	角速度に応じて回転させる

公式ドキュメント

ちなみにこのライブラリの作者はオイラー角が大ッキライみたいです．使ってる人見たらやめろ！って言ってあげてその場から立ち去ってお母さんにちくっちゃえ！とか書いてあります

最後に

Pythonでクォータニオンを扱う需要ってどこにあるんでしょうか？UnityはC#だしそもそもUnityの機能使えって話だし，ロボットではPythonは遅すぎて論外だし…