四元数の回転操作と空間群の表現(Julia)

今回は、Juliaで四元数、二重四元数を実装して、それを用いて空間群の対象操作を行ってみたいと思います。

四元数

四元数(Quaternion)とは、実数$s, u, v, w$と虚数単位$i, j, k$を用いて以下のように表せる数体系のことです。
$$q = s + ui + vj + wk$$
虚数単位は以下の条件を満たします。
$$i^2 = j^2 = k^2 = -1$$
$$ij = -ji = k, jk = -kj = i, ki = -ik = j$$
$q$を実数部$s$とベクトル部$\boldsymbol{v}$を用いて$q = (s, \boldsymbol{v})$と書くと便利です。
四元数は三次元空間での任意の回転操作を表現できることが知られています。

四元数の演算

和、差

$$q_1 \pm q_2 = (s_1 \pm s_2,\boldsymbol{v_1} \pm \boldsymbol{v_2})$$
スカラー積

$$q_1 \cdot q_2 = s_1s_2 + \boldsymbol{v_1}\boldsymbol{v_2}$$
ベクトル積

$$q_{1} \times q_{2} = (0, s_1\boldsymbol{v_2} + s_2\boldsymbol{v_1}+\boldsymbol{v_1}\times\boldsymbol{v_2})$$
クオータニオン積

$$q_1 \otimes q_2 = (s_1s_2 - \boldsymbol{v_1}\boldsymbol{v_2}, s_1\boldsymbol{v_1} + s_2\boldsymbol{v_1} + \boldsymbol{v_1}\times\boldsymbol{v_2}) = ((q_1\cdot q_2)^{*}, (q_1\times q_2))$$
共役

$$q^{*} = (s, -\boldsymbol{v})$$
ノルム

$$|q| = \sqrt{qq^{*}}$$

四元数による回転

導出は今回の本題ではないので、結果だけ示します。
行いたい回転の回転角$\theta$, 回転軸$\boldsymbol{\xi}$、対象となる点の座標を$\boldsymbol{r}$とする。

四元数

$$p = (0, \boldsymbol{r})$$
$$q = (\cos\frac{\theta}{2}, \xi\sin\frac{\theta}{2})$$

と共役$q^{*}$を用いて、

$$p' = qpq^{*}$$

とすると回転後の座標は$p'$のベクトル部として得られます。

例. x軸周りの4回回転操作(4回操作を行うと元の位置に戻る、90°回転)

$\theta = \frac{\pi}{2}$, $\xi = (1, 0, 0)$

$q = (\cos\frac{\pi}{4}, \xi\sin\frac{\pi}{4}) = (\frac{1}{\sqrt{2}}, \xi\frac{1}{\sqrt{2}}) = \frac{1}{\sqrt{2}} + i\frac{1}{\sqrt{2}}$

$q^{*} = \frac{1}{\sqrt{2}} - i\frac{1}{\sqrt{2}}$

$r = (0, 1, 0)$とすると、$p = (0, r) = 0 + 0i + 1j + 0*k = j$

$p' = qpq^{*} = (\frac{1}{\sqrt{2}}+i\frac{1}{\sqrt{2}})j(\frac{1}{\sqrt{2}}-i\frac{1}{\sqrt{2}}) = (\frac{1}{\sqrt{2}}+i\frac{1}{\sqrt{2}})(j\frac{1}{\sqrt{2}}+k\frac{1}{\sqrt{2}}) = k$

よって$r' = (0, 0, 1)$です。

Juliaでの実装

四元数による回転操作をJuliaで実装してみましょう。
まずは四元数の構造体を作ります。実数と実3次元ベクトルがメンバです。

mutable struct Quaternion <: Number
    s::Real
    v::Vector{Real}
end #Quaternion

<: Numberの部分で上位の型としてNumberを指定しています。
続いて四元数の演算を定義します。
Juliaでは多重ディスパッチを用いることで、既存の四則演算を自作の構造体に拡張できます。

import Base

Base.:+(a::Quaternion, b::Quaternion) = Quaternion(a.s + b.s, a.v + b.v)
Base.:-(a::Quaternion, b::Quaternion) = Quaternion(a.s - b.s, a.v - b.v)
Base.:-(a::Quaternion) = Quaternion(-a.s, -a.v)
Base.:*(a::Quaternion, b::Quaternion) = Quaternion(a.s*b.s - a.v⋅b.v, (a × b).v)

Juliaの四則演算はBaseで定義されているため、Baseをimportする必要があります(四則演算を使うだけならわざわざimportする必要はありません。書き換えるときのみ)。

スカラー積、ベクトル積も四元数に拡張します。

using LinearAlgebra

function LinearAlgebra.:dot(a::Quaternion, b::Quaternion)
    s = a.s * b.s
    v = dot(a.v, b.v)

    result = s + v

    return result
end #quaternion_dot

function LinearAlgebra.:cross(a::Quaternion, b::Quaternion)
    v = a.s*b.v + b.s*a.v + cross(a.v, b.v)

    result = Quaternion(0, v)

    return result
end #quaternion_cross

Juliaのスカラー積、ベクトル積はLinearAlgebraで定義されているのでusingします(通常、LinearAlgebraはusingしないと使えません)。LinearAlgebraで定義したdot、crossはコード中で　"$\cdot$"　、"$\times$"で呼び出せます。便利ですね。

ちなみに、筆者はusingとimportの違いがいまいちわかっていません。

共役も定義します。

function Base.:conj(a::Quaternion)
    result = Quaternion(a.s, -a.v)

    return result
end #Quaternion

最後に、四元数と他のNumberの接続を定義します。

Base.promote_rule(::Type{<:Number}, ::Type{Quaternion}) = Quaternion
Base.convert(::Type{Quaternion}, x::T) where {T<:Number} = Quaternion(x, [0, 0, 0])
Base.convert(::Type{Quaternion}, x::Quaternion) = x

四元数ができましたので、回転操作を行なってみます。先の例をJuliaでやってみましょう。
回転操作の四元数を作ります。

r = [0, 1, 0]
p = Quaternion(0, r)

θ = π/2
ξ = [1, 0, 0]

s = cos(θ/2)
v = ξ*sin(θ/2)
q = Quaternion(s, v)
q̄ = conj(q)

$p' = qpq^{*}$を計算し、ベクトル部を取り出します。

pp = q*p*q̄
result = pp.v

結果を確認してみましょう。

Real[0.0, 2.220446049250313e-16, 1.0]

元がr = (0, 1, 0)なので、ちゃんとできてそうです。
行列計算を行うことなく三次元空間での回転ができました。四元数による回転操作の利点ですね。

空間群

話は変わって空間群の話です。空間群とは、対称操作が元となっている群です。対称操作とは、原子や分子に変換操作を行ったとき、変換後の位置に同種の原子、分子がある操作のことです。

空間群は結晶点群(恒等操作、回転操作、鏡映操作、反転操作、回映操作、回反操作)と並進操作からなります。ヘルマン・モーガン記号、シェーンフリース記号という二種類の表現法があり、ヘルマン・モーガン記号は主に結晶の対称性を使う物性物理の分野で、シェーンフリース記号は分子の対称性を使う化学の分野でよく用いられます(正直どっちもわかりづらいので何使っても一緒です)。

以後は著者にとってどちらかといえば親しみのあるヘルマン・モーガン記号で説明します。

空間群のうち、結晶点群の操作について四元数で表現してみましょう。

回転操作

回転操作は先に示した方法そのままです。例のx軸周りの4回回転操作は回転操作のひとつで、4xと表されます。ヘルマン・モーガン記号では軸を略して次数だけを書くことも多いです(4xなら4だけとか)。

鏡映操作

鏡映操作とは、文字通り鏡写しです。鏡写しをするとき、(仮想的に)置かれる鏡の面を鏡映面といい、鏡映面で鏡映操作は定義されます。ヘルマン・モーガン記号ではmで表されます。

例　xy面を鏡映面とした鏡映操作

r = [0, 0, 1]
p = Quaternion(0, r)

q = Quaternion(0, [0, 0, 1])
q̄ = q

pp = q*p*q̄
print(pp.v)

結果を確認しましょう。

Real[0, 0, -1]

鏡映操作は回転操作と反転操作を組み合わせた、回反操作の特殊な場合と言えます。

反転操作

反転操作もそのまま、対象を反転させたもの(座標の符号を反転したもの)です。この操作については四元数を使うまでもありません。鏡映操作に近いですが、面の指定がないのでより自由度が高い操作です。

回映操作、回反操作はそれぞれ回転操作と鏡映操作、回転操作と反転操作の組み合わせであるため省略します。

これで結晶点群については四元数で表現できました。空間群には結晶点群に加えて並進操作があります。上記で述べた方法では並進操作を表現できませんので、二重数というものを導入します。

二重数

二重数(dual number)とは、$a, b \in \mathbb{R}$と$\epsilon \neq 0, \epsilon^2 = 0$となる$\epsilon$を用いて$z = a + b\epsilon$で表す非実数です。虚数単位を$\epsilon$に変更した複素数のようなものです。

二重四元数

二重数の実数a, bを四元数に置き換えたものです(八元数とか十六元数とかではありません、それらはほかにあります)。

四元数$q_p , q_d, \epsilon \neq 0, \epsilon^2 = 0$を用いて$\hat{q} = q_p + q_d\epsilon$と定義します。

二重四元数の演算

和、差

\hat{q}_1 \pm \hat{q}_2 = q_{1, p} \pm q_{2, p} + (q_{1, d} \pm q_{2, d})\epsilon

スカラー積

\hat{q}_1 \cdot \hat{q}_2 = q_{1, p}q_{2, p} + q_{1, d}q_{2, d}

ベクトル積

\hat{q}_1 \times \hat{q}_2 = q_{1, p} \times q_{2, p} + (q_{1, p}q_{2, d} +　q_{2, p}q_{1, d})\epsilon

共役(二重数)

\hat{p}　' = q_p - q_d\epsilon

共役(四元数)

\hat{q}^{\ast} = q_p^{\ast} + q_d^{\ast}\epsilon

共役(二重四元数)

\bar{\hat{q}} = q^{*}_p - q^{*}_d\epsilon

二重四元数を用いた対称操作では、並進操作と回転操作を別の四元数で表現します。

回転操作の四元数を$q_R$、並進操作の四元数を$q_I$、回転角$\theta$、回転軸$\xi$、並行移動$t$として、二重四元数を定義する。

q_R = (\cos\frac{\theta}{2}, \xi\sin\frac{\theta}{2}) \\
\hat{q} = q_r + \frac{1}{2}q_rq_I\epsilon

この二重四元数を使って、四元数の回転操作と同様に、
$p' = \hat{q}p\bar{\hat{q}}$
とします。
ただし、$p = 1 + p_d\epsilon$, $p_d = (0, \boldsymbol{r})$です。

これをJuliaで実装しましょう。
まず、二重四元数の構造体を作り、演算を定義します。

mutable struct DualNumber <: Number
    primary
    dual
end #DualNumber

Base.promote_rule(::Type{<:Number}, ::Type{DualNumber}) = DualNumber
Base.convert(::Type{DualNumber}, x::T) where {T <: Number} = DualNumber(x, 0)
Base.convert(::Type{DualNumber}, x::DualNumber) = x

Base.:+(a::DualNumber, b::DualNumber) = DualNumber(a.primary + b.primary, a.dual + b.dual)
Base.:-(a::DualNumber, b::DualNumber) = DualNumber(a.primary - b.primary, a.dual - b.dual)
Base.:-(a::DualNumber) = DualNumber(-a.primary, -a.dual)
Base.:*(a::DualNumber, b::DualNumber) = DualNumber(a.primary * b.primary, a.primary * b.dual + a.dual * b.primary)

共役も必要なので、関数として定義します。四元数のようにBase.conj()に上書きしてもいいのですが、四元数、二重数、二重四元数とそれぞれ別の共役が考えられるので、ここでは二重四元数の共役として別にします。

function dual_quater_conj(q::DualNumber)
    if (typeof(q.primary) == Quaternion && typeof(q.dual) == Quaternion)
        p = conj(q.primary)
        d = conj(q.dual)

        result = DualNumber(p, -d)
        return result
    end
end #dual_quater_conj

回転、並進操作の二重四元数を作ります。[0, 0, 1]方向に1/2(ここでは単位はないものとします)並進して、180°回転する(xz面で鏡映する)操作を例にします。

r = [0, 1, 0]
pd = Quaternion(0, r)
p = DualNumber(1, pd)

θ = π
ξ = [0, 0, 1]
t = [0, 0, 1/2]

qr = rotation(θ, ξ)
qi = Quaternion(0, t)

d = 1/2*qr*qi
q = DualNumber(qr, d)

q̄ = dual_quater_conj(q)

$q$とその共役を使って$p' = qp\bar{\hat{q}}$を計算します。

pp = q*p*q̄

操作後の座標は$p' = (1, 0) + \epsilon(0, \boldsymbol{r'})$の$\boldsymbol{r'}$として得られます。

result = pp.dual.v

Plotして確認してみましょう(わかりやすさのために二度映進操作してます)。

ちょっと注釈の方が目立っちゃってますが、(0, 1, 0)が一度の映進操作で(0, -1.0, 0.5)に、もう一度映進操作して(0, 1, 1)ができています。

ある方向に並進して180°回転する(並進した方向に並行な面に対して鏡映する)操作は空間群において映進と呼ばれる操作です。映進は並進と並行する鏡映面で表され、上の例はc映進面による映進操作と言われます。

まとめ

四元数を用いることで、行列計算を行わず3次元の回転操作を表現できました。Juliaは行列計算がそこそこ高速なので、たくさんのオブジェクトを回転させるような場合を除けばオイラー角での回転操作でも十分です。しかし、空間群の元と回転の演算子が1対1で対応するほうがすっきりするので、今回は四元数を採用しました。

余談ですが、これはJuliaLang Advent Calender 2022用に書いたものなんですが、同じアドベントカレンダーの数日前にQuaternion.jlについての記事が投稿されています。そんな便利なモジュールがあることすら知らなかった……。

Juliaで四元数を扱うならこの記事みたく構造体を作るより、Quaternion.jlを使った方が楽だと思います。

参考文献、参考記事

・クォータニオン｜CGのための数学
四元数の回転操作について詳しく書かれています。この記事では省いた導出もしっかりされていますのでとても参考になります。

・Quaternions.jlをメンテナンスしてる話
余談で述べたQuaternion.jlについての記事です。