投影幾何学的代数PGAによる3次元計算が楽で面白い

作成日: 2023-09-05
youtubeで見かけたPGAのチュートリアルを見て感動しました。3次元の座標系を計算するとき、手を動かすには複雑な概念がわさわさ出てきて嫌な気持ちになりますよね。それらが全部おなじ基底であつかえるとなると俄然楽しそうです。動画を参考にPGAの概要を紹介します。
Video preview

PGAの基底は平面から始まる

初めにすべての3次元平面(ax+by+cz+d=0)(ax+by+cz+d=0)を表す基底を考えます。
ex=(x=0),  ey=(y=0),  ez=(z=0),  eo=(const=0)e_x=(x=0),\ \ e_y=(y=0),\ \ e_z=(z=0),\ \ e_o=(\rm{const}=0)
ex,ey,ez,eoe_x,e_y,e_z,e_oの4本の基底の線形結合で任意の平面を表すことができます。
aex+bey+cez+deo=(ax+by+cz+d=0)ae_x+be_y+ce_z+de_o=(ax+by+cz+d=0)

PGAのMeetで共通部分が求まる

PGAの演算にはいくつか種類があります。その中の1つがMeetです。Meetは基底を次のルールで計算します。
ea∧ea=0ea∧eb=eaebe_a \wedge e_a = 0\\ e_a\wedge e_b =e_a e_b
同じ基底があれば0。そうでなければ基底を並べて新しい基底にします。eaeb=eabe_a e_b = e_{ab}と略記します。
eab=eaeb=−ebea=−ebae_{ab} = e_a e_b = -e_be_a=-e_{ba}
並んだ基底を1箇所入れ替えると符号を反転させます。
Meetが賢いのは、これだけで平面の共通部分である直線を計算できることです。(x+y−1=0)(x+y-1=0)と(z−1=0)(z-1=0)の共通部分を求めてみましょう
(ex+ey−eo)∧(ez−eo)=ex∧ez−ex∧eo+ey∧ez−ey∧eo−eo∧ez+eo∧eo=exz−exo+eyz−eyo−eoz+0=−ezx+eox+eyz+eoy−eoz(e_x+e_y-e_o)\wedge(e_z-e_o)\\=e_x\wedge e_z-e_x\wedge e_o+e_y\wedge e_z-e_y\wedge e_o -e_o\wedge e_z+e_o\wedge e_o\\=e_{xz}-e_{xo}+e_{yz}-e_{yo}-e_{oz}+0\\=-e_{zx}+e_{ox}+e_{yz}+e_{oy}-e_{oz}
なんだかあまり見慣れない形になりましたが、これが直線の表現らしいです。そもそも三次元の直線の表現ってなんだと思いましたが、どうやらPlücker座標というらしいです。
Plücker coordinates - Wikipedia
https://en.wikipedia.org/wiki/Plücker_coordinates
3次元の直線は向きのベクトルddと直積のベクトルmmで表現できます。Plücker座標が次のように展開されています。
dxeyz+dyezx+dzexy+mxeox+myeoy+mzeozd_xe_{yz}+d_ye_{zx}+d_ze_{xy}+m_xe_{ox}+m_ye_{oy}+m_ze_{oz}
先程の共通部分では次のようになります。
d=(1,−1,0)  m=(1,1,−1)d=(1,-1,0)\ \ m=(1,1,-1)
点A(1, 0, 1)から点B(0, 1, 1)への直線なので、d=A−B,  m=A×Bd=A-B, \ \ m=A\times Bを満たしてますね！すごいです！

Meetでもっと共通部分が求まる

Meetは平面と平面の共通部分だけではありません。直線と平面の共通部分も求められます。先程の直線と平面(x−2=0)(x-2=0)の共通部分も求めてみましょう。
(−ezx+eox+eyz+eoy−eoz)∧(ex−2eo)=eyzx+eoyx−eozx+2ezxo−2eyzo=2eozy−eoxz+eoyx+exyz(-e_{zx}+e_{ox}+e_{yz}+e_{oy}-e_{oz})\wedge(e_x-2e_o)\\= e_{yzx}+e_{oyx}-e_{ozx}+2e_{zxo}-2e_{yzo}\\=2e_{ozy}-e_{oxz}+e_{oyx}+e_{xyz}
∧\wedgeの両側に同じ基底があるともりもり消えます。計算された点のPGAでの表現は次のようになります。
xeozy+yeoxz+zeoyx+wexyzxe_{ozy}+ye_{oxz}+ze_{oyx}+we_{xyz}
ただし、wが1でないときは全体をwで割って正規化する必要があります。
直線と平面の共通部分の計算結果は点(2, -1, 1)になりました。(x+y−1=0)(x+y-1=0)と(z−1=0)(z-1=0)と(x−2=0)(x-2=0)をすべて満たすので間違いないですね！平面∧\wedge平面∧\wedge平面で点が求まるのはとても美しいです。

Joinでつながる点と点

ここまでは共通部分を求める演算でしたが、逆につなげる演算もできます。Meetに対して双対なJoinという演算を考えます。Joinのために、基底を反転する演算を考えます。
!exz=eoy!e_{xz}=e_{oy}
詳細を省きますが、↑みたいな感じです。Joinは次のように定義できます。
a∨b=!(!a∧!b)a\vee b=!(!a\wedge !b)
これを使うと点と点をつなぐことができます。A(1,0,1) B(0,1,1)を繋いでみましょう。
(eozy+eoyx+exyz)∨(eoxz+eoyx+exyz)=−ezx+eox+eyz+eoy−eoz(e_{ozy}+e_{oyx}+e_{xyz})\vee(e_{oxz}+e_{oyx}+e_{xyz})\\=-e_{zx}+e_{ox}+e_{yz}+e_{oy}-e_{oz}
さっきと同じ直線の表現が無事得られました。ここまで来たら同様に平面までつなげることができます。偉いですね。

回転とか変換とか

ここまでは三次元のモノどうしの関係を扱いました。実は回転や変換などもここまでの基底で表すことができます。詳細を省きますが、直線ll(軸dd)に対する回転子mmは次のようになります。
m=cos⁡θ2+sin⁡θ2ldm = \cos \frac{\theta}{2} + \sin\frac{\theta}{2}l_d
回転子mmを使って点ppを回転させるには両側から挟んで積を取ります。
p′=mpm~p'=mp\widetilde{m}
長くなってきたので、PGAの積やm~\widetilde{m}については次のページを参照してください。チートシートや計算機もあり非常に便利なサイトです。
bivector.net
https://bivector.net/tools.html?p=3&q=0&r=1
つまり、3次元回転でおなじみのクオータニオンもPGAの中に含まれています。いたせりつくせりです。16個の基底を簡単にこねくり回すだけで、これだけの表現ができるので感動します。
ちなみに、Rustの実装も存在します。使ってみたいですね。

おまけ

ここまでの知識で螺旋をくるくるするコードを書いてみました。行列のことを知らないで実装するのは不思議で楽しいですね。
codepen.io
https://codepen.io/wass80/pen/QWzbXoo
メモ: クリフォード代数だと平面を基底にするのではなく、点を基底にして始めるほうが普通らしい？