指数関数的に向上した回転技術 (2022)

 (thenumb.at)
1 ポイント 投稿者 GN⁺ 2024-06-16 | 1件のコメント | WhatsAppで共有
  • 3D回転は表現方式ごとに強みが異なるため、点の変換には回転行列が便利だが、補間・合成・平均には別の道具が必要になる
  • Euler angles は人が扱いやすい一方で、gimbal lock、一定でない角速度、最短経路を外してしまう線形補間の問題が起こりうる
  • 単位 quaternion は slerp により一定速度の最短経路補間を提供するが、ベクトル空間ではないため、直接の作成・スカラー倍・平均が直感的ではない
  • 指数/対数マップ は axis/angle ベクトルと回転行列を結びつけ、R(t) = exp(t log(R1 R0^-1)) R0 という形の 2D・3D 最短経路補間を構成する
  • 複数の回転の平均では、単純な axis/angle 平均だけでは catastrophic cancellation が起こりうるが、Karcher mean は角距離の二乗和を最小化する回転を反復的に求め、より一貫した結果を生む

回転表現ごとの長所と短所

  • 3D回転には複数の表現があり、変換・作成・補間・平均のどれを行うかによって適した選択が変わる

  • 回転行列

    • 最も直接的な線形代数の表現は、正の行列式を持つ orthonormal 3x3 matrix である
    • 回転行列の3つの列は、回転後に x, y, z 軸がどこへ移動するかを表す
    • 点の変換は行列積で処理でき、他の線形変換とも行列積で合成できる
    • 描画時に回転行列を使うのは、world-space から screen へ点を移すのに1回の行列積だけで済むためである
    • 回転行列はベクトル空間ではないため、2つの回転行列を足しても再び回転行列にはならない
    • 2つの回転行列を線形補間すると、回転だけでなく scaling が混ざることがある

  • Euler angles

    • Euler angles は x, y, z 軸に関する3つの回転を指定し、pitch、yaw、roll とも呼ばれる
    • 3つの構成回転の適用順序は慣例によって異なり、例では x, y, z の順序を使う
    • 人には理解しやすく回転の作成によく使われるが、単純な補間は望ましくない結果を生むことがある
    • ある構成回転によって残り2つの回転軸が平行になる gimbal lock は singularity に相当する
    • singularity では、固定された2つの角度のどちらを変えても同じ出力回転を作れてしまう
    • 補間経路が singularity に到達すると、現在位置を表現する自由度が増え、任意の表現を選んでつなぐと出力補間が不連続になることがある
    • 各構成角は cyclic であるため、線形補間が常に2つの回転のあいだの最短経路を選べるとは限らない
    • 経路が singularity を通らなければ補間は滑らかであり、「真上」と「真下」を表現する必要がなければ制約を回避できる

  • Quaternions

    • 単位 quaternion は回転の合成と補間の標準的な道具として使われる
    • spherical linear interpolation、すなわち slerp は、2つの quaternion のあいだで一定速度の最短経路を選ぶ
    • 単位 quaternion もベクトル空間ではなく、人が直接作成するのは難しく、補間計算にもコストがかかることがある
    • スカラー倍や平均についての直感的な概念も乏しい
    • quaternion は回転空間を double-cover するため、場合によっては Q(1) が -Q1 に向かうことがある

  • Axis/angle

    • axis/angle 回転は実数の3Dベクトルで表される
    • ベクトルの向きは回転軸を、大きさはその軸まわりの回転角を指定する
    • θu と書き、u は単位ベクトル、θ は回転角である
    • 3Dベクトルなので ベクトル空間 をなし、加算、スケーリング、補間が可能である
    • 2つの axis/angle 回転を線形補間すると、滑らかで一定の角速度を与えられる
    • ただし、目標回転をどの axis/angle 表現で指定するかによって、線形補間が最短経路を選べないことがある
    • quaternion と同様に、axis/angle ベクトルも回転空間を double-cover する

指数マップと対数マップ

  • 複数の回転表現を目的に応じて行き来できれば、各表現の利点を組み合わせて使える
  • 最終的な変換には回転行列が必要なので、行列を canonical form として置く
  • exponential map は回転オブジェクトを受け取り、等価な回転行列を返す関数である
  • logarithmic map は回転行列を受け取り、回転オブジェクトへ戻す対応関数である
  • ここでは、回転行列と axis/angle ベクトルのあいだを行き来する exp と log マップを扱う

2D axis/angle から出発する直感

  • 2Dでは回転軸は平面の外を向く1つしかないため、axis/angle 回転は角度 θ 1つで表せる
  • 2Dの点 p を θ だけ回転させた pθ は次のように書ける
    • pθ = p cosθ + Jp sinθ
    • J は 2D ベクトルを90度回転させる行列である
  • J は [[0, -1], [1, 0]] で、J² = -I なので2回適用すると180度回転になる
  • この式を展開すると、標準的な 2D 回転行列 [[cosθ, -sinθ], [sinθ, cosθ]] が得られる

2D 指数マップと対数マップ

  • 複素数の Euler の公式 e^(iθ) = cosθ + i sinθ で i が quarter turn の役割を果たすように、2D 行列式では J が同じ役割を果たす
  • 指数関数の Taylor series に行列 A = θJ を入れると、行列の加算、積、スケーリングによって同じ計算を行える
  • 展開の結果、sinθ と cosθ の Taylor series が現れ、次の式を得る
    • e^(θJ) = [[cosθ, -sinθ], [sinθ, cosθ]]
  • したがって 2D 指数マップは、角度 θ を対応する回転行列へ変換する
  • 対数マップは指数マップの逆として定義される
    • R = exp(θJ) なら log(R) = θJ
    • θ = atan2(R21, R11) によって復元できる
  • 指数マップは injective ではない
    • exp(θJ) = exp((θ + 2π)J) なので、1周加えても同じ回転行列になる
    • 対数マップは、その回転行列に対応する 最小の角度 を返すように定義される
    • atan2 がこの定義を実装している

exp/log ベースの補間

  • 2つの 2D 回転角 θ0、θ1 を単純に線形補間してから回転行列を作ることもできる
  • しかし θ0 と θ1 の差が π より大きいと、角度の cyclic な性質を反映できず、遠回りの経路を選んでしまう
  • exp/log ベースの補間では、2つの回転行列 R0、R1 から直接移動回転を計算する
    • R1 R0^-1 は、まず R0 を打ち消してから R1 を適用する回転である
    • log(R1 R0^-1) は、R0 から R1 へ行く最小角を与える
    • この axis/angle 回転を t 倍し、exp で再び行列に戻す
  • 最終的な補間式は次の通りである
    • R(t) = exp(t log(R1 R0^-1)) R0
    • R(0) = R0, R(1) = R1
  • 2Dでは角度差を直接確認してもよいが、この方法は修正なしで 3D や任意次元へ一般化できる

3D axis/angle と skew-symmetric matrix

  • 3Dでも axis/angle θu を指数化して回転行列を作れる
  • 核心は、単位ベクトル u を軸とする quarter turn 変換を見つけることにある
  • cross product u × p は u と p が作る平面に垂直なベクトルとして定義されるが、p を u に垂直な平面へ射影した p⊥ の quarter turn としても解釈できる
  • u × p と同じ結果を与える行列 û を作ることができる
    • û = [[0, -uz, uy], [uz, 0, -ux], [-uy, ux, 0]]
    • ûp = u × p
  • ûᵀ = -û なので、û は skew-symmetric matrix である
  • 2D の J も skew-symmetric で 2D cross product を表すので、同じ構造が続いている
  • skew-symmetric matrix の和とスカラー倍も skew-symmetric であるため、axis/angle ベクトル空間の性質はこの行列表現でも保たれる
  • û^(k+2) = -û^k という恒等式は、cross product を3回適用すると p⊥ を3回 quarter turn させることになり、負の quarter turn と等しいという幾何学的解釈から導かれる

3D 指数マップ: Rodrigues’ formula

  • axis/angle 回転 θu から θû を作って指数化すると、3D 回転行列が得られる
  • Taylor series と û^(k+2) = -û^k を使うと、次の式が導かれる
    • e^(θû) = I + sin(θ)û + (1 - cos(θ))û²
  • この式は Rodrigues’ formula として知られている
  • θ = 0 なら e^(0û)p = p となり、点はそのまま保たれる
  • θ = π/2 なら u × p + p∥ となり、quarter rotation になる
  • θ = π なら -p⊥ + p∥ となり、half rotation になる
  • この行列は orthonormal である
    • AᵀA = I の条件は ûᵀ = -û と û^(k+2) = -û^k から確認できる
  • 行列式は θ = 0 で 1 であり、行列式が 0 になる場合はなく、exp は θ と û に関して連続なので負になることもない
  • したがって exp(θû) は 3D 回転行列である

3D 対数マップ

  • 3D 指数マップも injective ではないため、3D 対数マップは、与えられた行列に対応する 最小の大きさ の axis/angle 回転を返すように定義される
  • R = exp(θû) = I + sin(θ)û + (1 - cos(θ))û² で trace を取ると回転角を求められる
    • trace は対角成分の和である
    • tr(I) = 3
    • û は skew-symmetric なので対角成分の和は 0 である
    • tr(û²) = -2 である
    • したがって tr(R) = 1 + 2cosθ
    • θ = arccos((tr(R) - 1) / 2)
  • 回転軸は R を antisymmetrize することで復元できる
    • R - Rᵀ = 2 sin(θ)û
    • û = (R - Rᵀ) / (2 sinθ)
    • u = 1/(2 sinθ) [R32 - R23, R13 - R31, R21 - R12]ᵀ
  • これで 3D 回転行列から axis/angle に戻る完全な対数マップが完成する

3D 補間の結果

  • 3D でも 2D と同じ補間式がそのまま適用できる
    • R(t) = exp(t log(R1 R0^-1)) R0
  • この補間は axis/angle 回転の利点を保ちながら、常に 最短経路 を選ぶ
  • Euler angles では同じ例が滑らかに見えないことがある

複数の回転の平均

  • quaternion を使えば、exp/log 行列の数学を使わなくても良い補間を得られるため、補間問題自体は解決できる
  • axis/angle 回転でより簡単にできる作業の1つが、複数の回転行列の 平均 である
  • 最も単純な方法は、各行列を axis/angle に変換し、ベクトルを平均してから戻すことである
  • この方法は有効ではあるが、直感に反する挙動を生むことがある
  • とくに axis/angle ベクトルを足し合わせると catastrophic cancellation が起こりうる
    • 例として、[π, 0, 0] と [-π, 0, 0] を平均すると 0 になる場合がある
    • 2つの値は等価な回転だが、平均結果の 0 はその2つの回転を代表していない

Karcher mean

  • 平面上の点の平均は、すべての点までの二乗距離の総和を最小化する点とみなせる
  • それを反復最適化で求める手順は次の通りである
    • 初期推定値 x̄ ∈ R² を選ぶ
    • 各点から推定値までの translation ui = xi - x̄ を計算する
    • ベクトル平均 u = (1/n) Σ ui を求める
    • x̄ = x̄ + τu として平均方向へ移動する
    • |u| > ε のあいだ繰り返す
  • 同じ考え方を回転 R0, ..., Rn に適用できる
    • 初期推定回転 R̄ ∈ R^(3×3) を選ぶ
    • 各行列について、推定値からその回転までの axis/angle ui = log(Ri R̄^-1) を計算する
    • ベクトル平均 u = (1/n) Σ ui を求める
    • R̄ = exp(τu) R̄ として平均回転方向へ移動する
    • |u| > ε のあいだ繰り返す
  • このアルゴリズムの結果が Karcher mean である
  • Karcher mean は、他のすべての回転までの二乗角距離を最小化する回転である
  • catastrophic cancellation の影響を受けず、常に 0 ではない中間的な回転へ収束する
  • 単純な axis/angle 平均と Karcher mean の結果はしばしば似ているが、Karcher mean のほうがより一貫した挙動を示す

Quaternion と exp/log の関係

  • この部分は quaternion の知識を前提とする
  • 複素数の指数化が skew-symmetric な 2D 行列の指数化と等価だったのと同様に、quaternion の指数化は skew-symmetric な 3D 行列の指数化と等価である
  • 2D では、axis/angle 回転 θ から pure-imaginary complex number iθ を作って指数化する
    • e^(iθ) = cosθ + i sinθ
    • 結果は点に掛けると θ だけ回転させる複素数である
    • 常に norm が 1 なので、2D 回転は unit-norm complex number で表現できる
  • 3D では、axis/angle 回転ベクトル u から pure-imaginary quaternion q = ux i + uy j + uz k を作れる
  • quaternion の乗算規則を使うと q² = -||q||² = -θ² となり、これは skew-symmetric matrix で使った恒等式に似ている
  • 指数化の結果は次の通りである
    • e^q = cosθ + (q/θ) sinθ
    • 2D の式とほとんど同じだが、imaginary axis が1つではなく3つある
  • 3D axis/angle 回転は unit-norm quaternion へ変換される
  • rotation matrix が不要なら、quaternion 指数マップは計算しやすい選択肢である
  • quaternion 対数マップも単純である
    • θ = arccos(Re(q))
    • u = Im(q) / sinθ
  • quaternion q で点 p を回転させるには、conjugation q p q^-1 を計算する
    • 点は p = px i + py j + pz k という pure-imaginary quaternion として表す
    • conjugation は厳密には u 軸まわりに 2θ 回転させるので、最初に |u| = θ/2 としておけばよい

さらに読む

関連記事

1件のコメント

 
GN⁺ 2024-06-16
Hacker News のコメント

  • リー群／リー代数対応は、学校で教わりたかった最高にクールな概念の一つ。記事で言っている指数写像と対数写像のことだが、はるかに再利用しやすい形で現れる
    3D 回転のように扱いたい抽象対象を、座標の細部に巻き込まれずに捉えるとそれがリー群で、そこからうまく機能する座標表現を導くと対応するリー代数になる
    すると、座標と抽象対象の間を行き来する方法や合成する方法などが、ほとんどただで手に入り、エンジニアリングでよく出会うケースでは補間や平均もかなり妥当に扱える
    問題をリー群の組み合わせとして表現できれば、それぞれの代数が何かを見つけることで、自力でやると時間のかかる作業をかなり節約できる
    ここでの対象には滑らかな変化という概念と少しの追加構造が必要で、行き来する過程では連結成分の問題が生じることもあるが、それは既知の結果を借りて使いやすい理由でもある

  • 長い一週間がほぼ終わりかけているところで、スライダーで牛を回転させるのはまさに必要としていた休憩だった

    • 数字が大量に見えた瞬間に目がかすんだが、牛は本当にかわいかった
    • 低労力の iOS キャッシュカウゲームにしたらちょうどよさそうな感じ

  • 多くの 3D ソフトウェアが回転に Arcball インターフェースを使っていないことは、昔からの不満
    Autodesk 製品の 3DSmax と Maya は使っているが、Blender と OpenSCAD は使っておらず、Roblox で働いていたときも、既存方式でもユーザーは耐えているという理由で PM を説得できなかった
    Arcball は四元数ベースで、補間に指数関数を使い、1 回のドラッグで任意の回転が可能で、ジンバルロックがなく、閉じたループを描いてドラッグすると開始位置に戻るという性質を持つ
    単位四元数が SO(3) の二重被覆になるという事実から数学的に証明でき、円周上の回転を単位複素数が正確に表すのと似ている
    直接触れる四元数／Arcball の参考実装: https://romankogan.net/math/arcball_js/index.html
    コードはコメントの多い Java で、Processing ライブラリを ProcessingJS により JavaScript 上で動かす形
    手と体は脳より先に四元数を理解できるので、3D ソフトウェアを作るならこの方式を使ってほしい
    Arcball: http://courses.cms.caltech.edu/cs171/assignments/hw3/hw3-not...
    回転用四元数: https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotati...
    Processing の Arcball: https://romankogan.net/math/arcball_js/index.html

    • 中央ではなく右側のような地点でドラッグすると、なぜずっと悪く感じるのか気になる。何か引っかかる感じがして、ときどき跳ねるのだが、回転を計算する基準点が変わるなら、その意味は画面上で視覚的に示されるべきだと思う
    • モバイルでも動くようにできるのか気になる
      今モバイルで見ているが、https://asliceofrendering.com/camera/2019/11/30/ArcballCamer... が Arcball を理解する助けになった

  • この文脈で、なぜ四元数がそこまで好まれるのかよく分からない。行列が四元数より直感的でないとは言いにくい
    行列はベクトルに作用し、回転もベクトルに作用するのだから、回転を行列として見ること以上に自然なものが何かあるのかと思う
    行列指数も常微分方程式と結び付ければ直感的だ。dx/dt = Ax の解は exp(t A) であり、A が反対称なら x の変化は常に x に直交するため、長さを変えない回転になる
    リー群／リー代数はこれを大きく一般化するが、核心は直交する変化を連続的に作って回転を生成し、指数写像がその過程を説明するという点にある。この図式のほうが、はるかに幾何学的で直感的に感じる

    • 四元数自体も Pauli 行列のように行列で表現されることが多く、これは量子スピンのモデリングで広く使われている
      四元数の利点は、同じ行列乗算を手で計算するより、紙とペンだけで扱いやすい点にあると思う
      個人的には複素数と似た意味で直感的だ。最初は奇妙だったが、今では知っている代替手段よりも使いやすく、推論もしやすいと感じる
    • 記事でも、四元数が行列に対して持つ大きな利点を説明していた。四元数は補間がうまくいき、行列はそうではない
      この性質は、アニメーションや 3D スプライン曲線に沿ってフレームを計算するグラフィックス作業で非常に重要
    • 計算の観点での四元数の利点は、3x3 行列の 9 個の数の代わりに4 個の数だけを使い、回転の適用も同様に演算量が減る点にある
    • 四元数の核心的な利点は回転の合成
      2D 回転を複素数で扱うのと同じように、2 つの複素数を掛けると回転が合成され、2D では偏角が足し合わされる形になる。同様に、2 つの四元数を掛けると 3D 回転を合成でき、3x3 行列の積よりはるかに効率的
      直感のために言うと、四元数は軸角表現と密接で、これはリー代数 so(3) と同じ
      ベクトルに作用するという観点では、複数の回転パラメータ化を、同じ抽象 Rotation トレイトの実装として見ればよい。内部実装が行列、四元数、オイラーベクトル、オイラー角、Gibbs ベクトルのいずれであっても、回転はベクトルに作用し、合成される方法は同じ

  • 単位四元数はリー群であり、任意に加算できるものが欲しいなら、回転速度を表す全四元数のリー代数を見るべき。これは軸角の回転速度を表すのと同じ
    単位四元数と軸角を比較するのはカテゴリが少しずれていて、単位四元数は回転行列と、全四元数は軸角と比較するほうがより適切
    四元数を使うと指数写像を簡単に計算できる利点があるが、四元数を使う場合、回転行列はほとんど必要ない。記事のように pqp^-1 で回転を計算できる
    四元数を理解する最も簡単な道は、幾何代数を読むことだと思う。四元数の発明には数百年かかったが、驚くほど単純な幾何代数を理解すれば、数分で四元数を再発明できる
    数年前によい入門として読んだ記事: https://crypto.stanford.edu/~blynn/haskell/ga.html

    • 回転を考えるときは、今でも 指数写像 が最も堅牢な方法だと思う。SO(3) のリー群を、座標変換、微分、接空間の扱いまで最も直接的に扱えるようにしてくれるから
      幾何代数のさまざまな定式化を経ても、結局 SO(3)/SE(3) 空間を表現するためにローターとモーターを使うことになり、これらはそれぞれ四元数と双対四元数に同型
      ただ、その目的なら、指数写像を伴う 3x3 回転行列と 4x4 変換行列のほうが、今でもはるかに有用だと思う。四元数は保存領域が少なく、互いの乗算も速いが、点を変換するときは行列のほうが速く、全体の効率は状況によって変わる
    • ここで「全四元数」と言っているのが 純虚四元数 のことなのか気になる

  • 大学で学んだすごいことの一つは、+ 演算子を行列とベクトル空間の変化量に合わせて、- 演算子を2つの行列に合わせて再定義すれば、Kalman フィルタの状態 の中に回転行列をそのまま入れられるという点だった
    こうすれば、ジンバルロックを心配せずに回転を推定できる
    https://openslam-org.github.io/MTK

    • SO(2) の 接空間 を状態に持つという意味なのか気になる
    • 接空間で + 演算子を実装するのは、Kalman フィルタだけでなく非線形最適化全般でもかなり一般的。Ceres ライブラリもそのための LocalParameterization をサポートしている

  • 本当に良かったし、一部だけが良かったわけではない
    特に、これらの方法が結局 標準の回転行列 を計算するという点が良かった。100万個のベクトルを回転させる必要があるなら、面白い計算は一度だけ行い、その後は高度に最適化された行列乗算パイプラインを回せばいい

  • すばらしいブログだったが、著者プロフィールをクリックして他の記事を見ていたら、「2010年ごろ、9歳のときに初めてプログラミングに触れた」という文を見た
    私は2010年には13歳で、中等教育の数学と科学を頭に詰め込もうと必死だった
    すばらしいコンピュータグラフィックス関連の記事を見るたびに、自分より若く、はるかに才能のある人が書いたように思えて 強い劣等感 を覚える

    • まだ遅くない。ここ数年で独学でニッチな技術を身につけ、ある大企業がそれにお金を払って使いたがっているようだ。年齢は35歳くらい
      ただ、どうすればいいかという助言はあまりできない。単にそれが本当に好きで続けてきただけだから。それでも、構造化された練習のようなものも効く可能性はありそう

  • 複数の回転の平均を求める方法を探していて、https://mathweb.ucsd.edu/~sbuss/ResearchWeb/spheremean/paper... を見つけた
    この記事の方法は、少なくとも私の数学レベルでは、その論文よりずっと易しそうに見える

    • 複数の回転の平均を求めようとしている 文脈 が気になる
      平均は和でできる操作であり、和では合成順序は重要ではない。しかし回転は交換法則に従わないので、私たちが普通に理解している平均の概念はそのまま適用できない
      スマートフォンを持って画面を180°回して向こう側に向けたあと、地面基準で時計回りに90°回すと、カメラは左を向く。逆に同じ2つの回転を逆順に行うと、カメラは右を向く
      この2つの回転の「平均」でカメラがどこを向くべきかには単一の答えはなく、平均に求める性質によって変わる
    • 数日前、平面回転に関して少し脇道のコメントを書いたことがある。平面ではずっと単純になり、プログラミング環境が複素数演算を第一級の機能としてサポートしていると特に楽になる
      https://news.ycombinator.com/item?id=40333541
      核心となる直観は、加算は平行移動で、乗算は回転だということ。だから平均平行移動には算術平均を使い、平均回転には 幾何平均 を使える

  • 数学でも、ソフトウェア工学で抽象化を思い浮かべるのと似た形で 抽象化 を作り出しているのだと気づくまで、長い時間がかかった
    子どものころは、なぜ虚数を作るのか、行列はいったい何の意味があるのか混乱していた
    後になってようやく、こうした表現は設計されたものなのだと分かった。虚数というものを作るとある種の計算が簡単になり、線形方程式を行列で書くと、全部を展開して書くよりはるかに推論しやすくなる
    当たり前に見えるが、誰もこうは教えてくれなかった

    • もっと抽象的な数学を学んでも、この見方はずっと当てはまる。たとえば 群 はインターフェースであり、「ある群」は必要な3つの性質/メソッドを実装する演算を持つ型だと見なせる
      ベクトル空間、環、距離空間、圏も同じで、数学は設計パターンとしてのインターフェースで満ちている
      ただし数学のインターフェースは、プログラミングの継承よりも型クラスに近い。同じ集合/型でも複数の方法で群になり得るから