Luna Flow
日本語

core の設計

設計目標

core は geometry3d のほかのすべてのパッケージが前提とする部分です。メッシュとは何か、点が変換でどう動くか、面のどちら側が外側か、光のもとで面がどれだけ明るいか。一度に読み通せるほど小さく、規約について正確で、カメラ、画面、出力デバイスに属するものを一切含まないことが求められます。また、Luna-Flow/linear-algebra の密な型だけで 3D パイプラインを支えられることを示す例でもあります。

数学的背景

同次座標

R3\mathbb{R}^3 のアフィン写像 x↦Ax+tx \mapsto A x + t は線形ではないので、3×3 の行列を持ちません。R3\mathbb{R}^3 を R4\mathbb{R}^4 に埋め込むとこれが解決します。点 pp を (p,1)(p, 1)、方向 dd を (d,0)(d, 0) と書くと、アフィン写像は次のブロック行列になります。

M=(At0T1),M(p1)=(Ap+t1),M(d0)=(Ad0).M = \begin{pmatrix} A & t \\ 0^\mathsf{T} & 1 \end{pmatrix},\qquad M \begin{pmatrix} p \\ 1 \end{pmatrix} = \begin{pmatrix} A p + t \\ 1 \end{pmatrix},\qquad M \begin{pmatrix} d \\ 0 \end{pmatrix} = \begin{pmatrix} A d \\ 0 \end{pmatrix}.

最後の座標はベクトルがどの種類の対象かを記録します。方向は 2 点の差 (p,1)−(q,1)=(p−q,0)(p, 1) - (q, 1) = (p - q, 0) であり、上の計算から AA だけで動かされることがわかります。両方の点を平行移動しても差は変わりません。これがまさに Transform3::apply_point（w=1w = 1 を付ける）と Transform3::apply_direction（w=0w = 0 を付ける）の区別です。

一般の 4×4 行列は最終行が (hT,k)≠(0,0,0,1)(h^\mathsf{T}, k) \ne (0, 0, 0, 1) のこともあります。そのとき (p,1)(p, 1) は w′=h⋅p+kw' = h \cdot p + k として (x′,y′,z′,w′)(x', y', z', w') に送られ、それが表す点は同次除算 (x′/w′,y′/w′,z′/w′)(x'/w', y'/w', z'/w') で求まります。任意の λ≠0\lambda \ne 0 について (x′,y′,z′,w′)(x', y', z', w') と λ(x′,y′,z′,w′)\lambda (x', y', z', w') は同じ点を表すので、これらの行列は射影空間に作用します。view が使う透視投影もこれに含まれます。apply_point は常に除算するので両方の場合を扱え、∣w′∣≤ε|w'| \le \varepsilon = DEPTH_EPSILON のとき、つまり無限遠平面（またはその近く）に送られた点では除算を省きます。

合成

2 つの変換 M1M_1（先に適用）と M2M_2（後に適用）について、行列積の結合法則により合成は列ベクトルに M2(M1v)=(M2M1)vM_2 (M_1 v) = (M_2 M_1) v として作用します。Transform3::compose は M2M1M_2 M_1 を格納します。

t1.compose(t2).matrix=M2M1.\texttt{t1.compose(t2)}.\mathrm{matrix} = M_2 M_1 .

メソッドは適用順に読み、行列積は右から左に読みます。代数から 2 つのことが従います。合成は結合的なので、連鎖は自由にまとめられます。可換ではありません。TT を tt による平行移動、RR を回転とすると、

TR(p1)=(Rp+t1)≠RT(p1)=(Rp+Rt1)T R \begin{pmatrix} p \\ 1 \end{pmatrix} = \begin{pmatrix} R p + t \\ 1 \end{pmatrix} \quad\ne\quad R T \begin{pmatrix} p \\ 1 \end{pmatrix} = \begin{pmatrix} R p + R t \\ 1 \end{pmatrix}

であり、Rt=tR t = t でない限り両者は異なります。そのため通常のモデル変換は、拡大縮小、回転、平行移動の順に行います。s.compose(r).compose(t) で、行列は TRST R S です。

基本回転

線形写像は基底ベクトルの像で決まり、それが行列の列になります。xyxy 平面を zz 軸まわりに θ\theta 回転させると、次のように写ります。

ex↦(cos⁡θ,sin⁡θ,0),ey↦(−sin⁡θ,cos⁡θ,0),ez↦ez,e_x \mapsto (\cos\theta, \sin\theta, 0),\qquad e_y \mapsto (-\sin\theta, \cos\theta, 0),\qquad e_z \mapsto e_z ,

これが rotation_z の行列 Rz(θ)R_z(\theta) を与えます。xx と yy まわりの回転は軸を x→y→z→xx \to y \to z \to x と巡回させて得られます。xx まわりでは (y,z)(y, z) が、yy まわりでは (z,x)(z, x) が (x,y)(x, y) の役割を果たします。2 つ目の置き換えのため、RyR_y の負号は対角線の下にあります。

Ry(θ):ez↦(sin⁡θ,0,cos⁡θ),ex↦(cos⁡θ,0,−sin⁡θ).R_y(\theta):\quad e_z \mapsto (\sin\theta, 0, \cos\theta),\qquad e_x \mapsto (\cos\theta, 0, -\sin\theta).

各 RR は行列式 +1+1 の直交行列です。列は正規直交で、det⁡Rz(θ)=cos⁡2θ+sin⁡2θ=1\det R_z(\theta) = \cos^2\theta + \sin^2\theta = 1 です。したがって R−1=RT=R(−θ)R^{-1} = R^\mathsf{T} = R(-\theta) であり、回転は長さ、角度、座標系の向きを保ちます（右手系の 3 つ組を右手系の 3 つ組に写します）。

オイラー角

rotation_matrix(α, β, γ) は基本回転を R=Rz(γ)Ry(β)Rx(α)R = R_z(\gamma) R_y(\beta) R_x(\alpha) と合成します。対応する角度の c∙=cos⁡c_\bullet = \cos、s∙=sin⁡s_\bullet = \sin と書いて掛け合わせると、

Ry(β)Rx(α)=(cβsβsαsβcα0cα−sα−sβcβsαcβcα),R=(cγcβcγsβsα−sγcαcγsβcα+sγsαsγcβsγsβsα+cγcαsγsβcα−cγsα−sβcβsαcβcα).\begin{aligned} R_y(\beta) R_x(\alpha) &= \begin{pmatrix} c_\beta & s_\beta s_\alpha & s_\beta c_\alpha \\ 0 & c_\alpha & -s_\alpha \\ -s_\beta & c_\beta s_\alpha & c_\beta c_\alpha \end{pmatrix},\\[4pt] R &= \begin{pmatrix} c_\gamma c_\beta & c_\gamma s_\beta s_\alpha - s_\gamma c_\alpha & c_\gamma s_\beta c_\alpha + s_\gamma s_\alpha \\ s_\gamma c_\beta & s_\gamma s_\beta s_\alpha + c_\gamma c_\alpha & s_\gamma s_\beta c_\alpha - c_\gamma s_\alpha \\ -s_\beta & c_\beta s_\alpha & c_\beta c_\alpha \end{pmatrix}. \end{aligned}

列ベクトルに作用させると、xx まわり、yy、zz の順に、それぞれ固定されたワールド座標軸について回転します（外因的な xx-yy-zz、同値な内因的 zz-y′y'-x′′x''）。回転の積は回転なので、オイラーの回転定理により RR はある軸まわりのある角度 θ\theta の単一の回転で、tr⁡R=1+2cos⁡θ\operatorname{tr} R = 1 + 2\cos\theta です。11 トレースは基底の変換で不変であり、回転軸を第 3 の基底ベクトルとする基底では行列は Rz(θ)R_z(\theta) で、そのトレースは 1+2cos⁡θ1 + 2\cos\theta です。geometry3d は軸と角度を取り出しませんが、この恒等式はテストで合成した回転を確かめるのに役立ちます。

オイラー角には特異点があります。β=π/2\beta = \pi/2 では行列は次のようになります。

R=(0sin⁡(α−γ)cos⁡(α−γ)0cos⁡(α−γ)−sin⁡(α−γ)−100),R = \begin{pmatrix} 0 & \sin(\alpha - \gamma) & \cos(\alpha - \gamma) \\ 0 & \cos(\alpha - \gamma) & -\sin(\alpha - \gamma) \\ -1 & 0 & 0 \end{pmatrix},

これは α−γ\alpha - \gamma だけに依存します。α\alpha と γ\gamma を同じ量だけ変えても同じ回転になります。自由度が 1 つ失われ（ジンバルロック）、β=±π/2\beta = \pm\pi/2 の近くでは姿勢の小さな変化に角度の大きな変化が必要になります。デモでは 3 つの角度をそれぞれ一定の速さで独立に動かすので、これは問題になりません。

面、法線、向き

面 (a,b,c,d)(a, b, c, d) は 4 頂点が同一平面上にあるとき平面上にあります。法線は最初の 3 頂点から計算します。

n=(vb−va)×(vc−va),n^=n/∥n∥.n = (v_b - v_a) \times (v_c - v_a),\qquad \hat n = n / \lVert n \rVert .

外積は反対称なので、2 頂点を入れ替えると nn は逆向きになります。頂点の順序こそが向きです。生成関数はすべての面を n^\hat n が立体の外を向くように並べます。z=−sz = -s にある立方体の面 (0,3,2,1)(0, 3, 2, 1) では、

v3−v0=(0,2s,0),v2−v0=(2s,2s,0),(0,2s,0)×(2s,2s,0)=(2s⋅0−0⋅2s, 0⋅2s−0⋅0, 0⋅2s−2s⋅2s)=(0,0,−4s2),\begin{aligned} v_3 - v_0 &= (0, 2s, 0), \qquad v_2 - v_0 = (2s, 2s, 0),\\ (0, 2s, 0) \times (2s, 2s, 0) &= (2s \cdot 0 - 0 \cdot 2s,\ 0 \cdot 2s - 0 \cdot 0,\ 0 \cdot 2s - 2s \cdot 2s) = (0, 0, -4s^2), \end{aligned}

となり、立方体から離れる −z-z を向きます。トーラス P(u,v)=((R+rcos⁡v)cos⁡u, rsin⁡v, (R+rcos⁡v)sin⁡u)P(u, v) = ((R + r\cos v)\cos u,\ r \sin v,\ (R + r\cos v)\sin u) では、(u,v)(u, v) の面は頂点 P(u,v)P(u, v)、P(u,v+Δv)P(u, v + \Delta v)、P(u+Δu,v+Δv)P(u + \Delta u, v + \Delta v) を使います。1 次の近似で 2 つの辺は PvΔvP_v \Delta v と PvΔv+PuΔuP_v \Delta v + P_u \Delta u なので、

n≈(Pv×Pu) Δu Δv,n \approx (P_v \times P_u)\, \Delta u\, \Delta v ,

であり、u=v=0u = v = 0 では Pv=(0,r,0)P_v = (0, r, 0) と Pu=(0,0,R+r)P_u = (0, 0, R + r) から Pv×Pu=(r(R+r),0,0)P_v \times P_u = (r (R + r), 0, 0) となります。法線は軸から離れ、管の外を向きます。背面カリングが内壁を表示するように退行しないよう、リポジトリのテストでこれを確かめています。

生成関数が出力する四角形はすべて平面なので、法線はきちんと定まります。球とトーラスでは、パラメータ u,u′u, u' の間の面は、yy 軸を含み角度 (u+u′)/2(u + u')/2 にある平面についての鏡映で対称です。この鏡映は P(u,v)↔P(u′,v)P(u, v) \leftrightarrow P(u', v) と P(u,v′)↔P(u′,v′)P(u, v') \leftrightarrow P(u', v') を入れ替えます。線分 P(u,v)P(u′,v)P(u, v)P(u', v) と P(u,v′)P(u′,v′)P(u, v')P(u', v') はどちらもその鏡映面に垂直なので平行であり、2 本の平行な線分は 1 つの平面を張ります。立方体、円柱の側面、退化した四角形は作り方から平面です。

背面の可視判定

点 ee から平面の面の表側が見えるのは、その点が法線の向く開半空間にあるときに限ります。

n^⋅(e−q)>0for a point q of the plane.\hat n \cdot (e - q) > 0 \quad\text{for a point } q \text{ of the plane.}

符号は qq の選び方に依存しません。平面上の 2 点 q,q′q, q' について n^⋅(q−q′)=0\hat n \cdot (q - q') = 0 なので、n^⋅(e−q)=n^⋅(e−q′)\hat n \cdot (e - q) = \hat n \cdot (e - q') です。face_is_visible は qq = face_center を使います。閉じた立体では、視点に背を向けた面は同じ立体の前面に隠されるので、それを除いても見える面が失われることはなく、ラスタライザの仕事はおよそ半分になります。異なる前面どうしの遮蔽は解決しません。それは深度バッファの仕事です。

ランバートシェーディング

単位法線 n^\hat n を持つ面積 AA の小さな面の一部が、単位方向 ℓ\ell と逆向きに進む平行光で照らされると、光線に垂直な面積 Acos⁡θA \cos\theta を通る光束を受け取ります。ここで cos⁡θ=n^⋅ℓ\cos\theta = \hat n \cdot \ell です。完全拡散（ランバート）面は光をすべての方向に等しく反射するので、見かけの明るさは cos⁡θ\cos\theta に比例し、光源が面の裏側にあるときは 0 です。

I=max⁡(0, n^⋅ℓ)∈[0,1].I = \max(0,\ \hat n \cdot \ell) \in [0, 1] .

face_intensity はまさにこれを返します。だから light は光源を向く単位ベクトルでなければなりません。面ごとに 1 つの値を使うとフラット（ファセット状）なシェーディングになります。

設計上の判断

標準の位相としての四角形

課題は、すべての生成関数、カリングの判定、ラスタライザが理解できる 1 種類の面の型を決めることです。選択肢は三角形のみ、四角形のみ、一般の多角形でした。四角形を選んだのは、生成関数のパラメトリック曲面（球、円柱、トーラス）が (u,v)(u, v) のグリッドで、そのセルが四角形だからです。1 つの四角形は 1 つの法線と 1 つのシェーディング値を持ち、フラットシェーディングの面の数が半分になります。三角形は d=ad = a の退化した四角形として埋め込むので、角錐、円錐の側面、球の極の部分も同じ型に収まります。triangulate_quad はラスタライザが必要とするところでだけ四角形を 2 つの三角形にし、退化した四角形の面積 0 の 2 つ目の三角形は面積の判定で落とされます。

点と方向を別の操作にする

Transform3 は 4×4 行列を 1 つ格納し、呼び出し側に ww を組み立てさせる代わりに apply_point と apply_direction を提供します。これにより同次座標の規約はパッケージ内に閉じ、「方向は平行移動を無視する」という規則を呼び出し側で間違えることがなくなります。

法線は変換せずに計算し直す

一様でない拡大縮小のもとでは、法線は方向と同じようには変換されません。法線の正しい行列は (A−1)T(A^{-1})^\mathsf{T} です。面の接ベクトル tt は n⋅t=0n \cdot t = 0 を満たし、写像 t′=Att' = A t の後でもベクトル n′=A−Tnn' = A^{-\mathsf{T}} n は n′⋅t′=nTA−1At=0n' \cdot t' = n^\mathsf{T} A^{-1} A t = 0 を保つからです。core は法線を保持してこの行列で変換する代わりに、毎回変換後の頂点から n^\hat n を計算し直します（apply_mesh の後に face_normal）。コストは面ごとに外積 1 回で、利点は一様でない拡大縮小を含め det⁡A>0\det A > 0 のどんなアフィン写像でも法線が常に正しいことです。det⁡A<0\det A < 0 の写像（鏡映）は巡回順を逆にし、すべての面を裏返します。

軸と角度や四元数ではなくオイラー角

デモで必要なのは物体を 3 つの軸まわりに独立した速さで回すことだけで、それはまさにオイラー角が表すものであり、行列はそれぞれ 3 行で書けます。軸と角度からの構成（ロドリゲスの公式）や四元数の合成は実装していません。それらは Luna-Flow/quaternion の役割です。ほかで作った回転行列は Transform3::from_matrix で渡せます。

共有の許容誤差

DEPTH_EPSILON = 10−910^{-9} は、パイプライン内のあらゆる「これは 0 か」という判断に使われます。正規化、同次除算、三角形の面積、そして厳密な深度テスト depth + ε < stored です。定数が 1 つなので、退化の境界での振る舞いがパッケージ間で一貫します。絶対的な許容誤差なので、シーンの座標が 11 から 10310^3 程度であることを前提とします。デモはおよそ 1 程度の単位を使います。

正しさと不変条件

  • アフィン変換に対する apply_point は浮動小数点の丸めを除いて正確です。w′=1w' = 1 なので除算による誤差は生じません。
  • compose はアフィン変換について t1.compose(t2).apply_point(p)=t2.apply_point(t1.apply_point(p))\texttt{t1.compose(t2).apply\_point}(p) = \texttt{t2.apply\_point}(\texttt{t1.apply\_point}(p)) を満たし、射影変換については同次スケールを除いて満たします。平行移動の後に拡大縮小を行うテストで順序を確かめています。
  • 回転行列は行列式 11 の直交行列で、その積もまた回転です。
  • 生成されたメッシュは閉じており、原点を中心とし、面は平面で法線は外向きです。face_is_visible と face_intensity はこれを前提にしています。
  • normalize_vec と face_normal は DEPTH_EPSILON より小さい数で割ることはありません。退化した入力は零ベクトルになり、その面は不可視で照らされません。
  • どの操作も新しいベクトルかメッシュを割り当て、引数を変更する関数はありません。apply_mesh は入力と面配列を共有します。

生成関数と apply_mesh（頂点数と面数に比例）を除き、各関数は定数時間です。

採用しなかった案

  • 専用の小さなベクトル型（フィールド x、y、z を持つ Vec3 構造体）は速く、次元について型安全ですが、linear-algebra と重複します。このリポジトリの目的は Luna-Flow の基盤の上に作ることなので、ベクトルは @la.Vector[Double] で、次元は規約にとどめます。
  • ジェネリックなスカラー（任意の Field 上の Mesh[T]）は採用しませんでした。三角関数、平方根、許容誤差はすべて Double を前提とし、ほかの型を使えるバックエンドもないからです。
  • メッシュに法線を保持することは上で述べた理由で採用しませんでした。計算し直すのは安価で、変換後も常に正しいからです。
  • 三角形メッシュを主な型にすることは、フラットシェーディングされたグリッドの各面を、同じ法線を持つ 2 つの面にしてしまいます。

境界

core は次のことをしません。

  • カメラ、投影、ビューポート、ターミナル、色、DOM について知ること。
  • パッケージ外のコードに任意のメッシュを作らせること。Mesh と QuadFace のフィールドは外からは読み取り専用なので、メッシュの供給元は生成関数だけです。
  • メッシュの読み込みや保存、シーングラフ、マテリアル、テクスチャ、物理、空間インデックスの提供。
  • スムーズな（頂点ごとの）法線の計算、幾何のクリッピング、交差判定。
  • 逆変換、軸と角度による回転、四元数による回転の提供。

Footnotes

  1. トレースは基底の変換で不変であり、回転軸を第 3 の基底ベクトルとする基底では行列は Rz(θ)R_z(\theta) で、そのトレースは 1+2cos⁡θ1 + 2\cos\theta です。geometry3d は軸と角度を取り出しませんが、この恒等式はテストで合成した回転を確かめるのに役立ちます。 ↩