Unityで超立方体の投影・回転・切断をしてみた話

# はじめに
四次元図形、いいですよね。
何だそれは、という人は「超立方体」「超球」なんかでググると未知の世界が拓けると思います。
私が初めて四次元に触れたのは[この動画シリーズ](https://www.nicovideo.jp/mylist/11674226)でした。

四次元図形を使った面白い表現をUnityでできないかな、投影とか切断とかどうやってるんだろう、と思って諸々調べて実装したので、まとめておこうと思います。
扱う図形は[四次元の超立方体（正八胞体）](https://ja.wikipedia.org/wiki/%E6%AD%A3%E5%85%AB%E8%83%9E%E4%BD%93)のみ。理論の方は深入りせず、実装メインで。

↓こんなのを作りました。
<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr">四次元遊びまとめ<br>黒球：超立方体の頂点（透視投影）<br>水色線：超立方体の辺（透視投影）<br>水膜：超立方体の面（透視投影）<br>赤線：超立方体の切断面（u=0）<br>QL, QRは６軸の回転速度 <a href="https://t.co/Jej4D3J28J">pic.twitter.com/Jej4D3J28J</a></p>&mdash; りじちょー (@Rijicho_nl) <a href="https://twitter.com/Rijicho_nl/status/1098316500732043264?ref_src=twsrc%5Etfw">2019年2月20日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

※シェーダは[NextGen Sprites](https://assetstore.unity.com/packages/tools/sprite-management/nextgen-sprites-47482)を使用。無料Asset。

# 四次元座標型を作る
Unityで普段使っている座標系は三次元です。しかしこれから扱う図形は四次元なので、四次元座標を表す型を自分で用意する必要があります。
Vector4を使う手もありますが、どうやら四次元座標の回転処理のときにクォータニオン（四元数）を使うと都合がいいっぽい。
しかしUnityのQuaternionは勝手に余計な処理を挟んでくるみたいなので、ここではオレオレQuaternionを作ることにします。


## 定義とか
いつもの$x,y,z$軸に第四次元の$u$軸を加えた四次元座標$(u,x,y,z)$を、クォータニオン $q=u+xi+yj+zk$という形で表すことにしましょう。


```cs:MyQuaternion.cs
//くぉーたにおん！
public struct MyQuaternion
{
    //q = u + xi + yj + zk
    public float u;
    public float x;
    public float y;
    public float z;

    //コンストラクタ
    public MyQuaternion(float u, float x, float y, float z)
        => (this.u, this.x, this.y, this.z) = (u, x, y, z);
}
```

なお、回転以外の文脈では単なる四次元ベクトルとしても扱います。

## オイラー角からの変換
UnityでQuaternionを作る場合、クォータニオンを完全に理解した~~変態~~プロ以外は四引数コンストラクタは使わないでしょう。
割とよく使うであろうQuaternion.Euler(x, y, z)に相当するコンストラクタも作っておきます。

```cs:MyQuaternion.cs
    //Quaternion.Euler(x,y,z)のやつ
    public MyQuaternion(float x, float y, float z)
    {
        var rx = x * Mathf.Deg2Rad;
        var ry = y * Mathf.Deg2Rad;
        var rz = z * Mathf.Deg2Rad;
        var cx = Mathf.Cos(rx / 2);
        var sx = Mathf.Sin(rx / 2);
        var cy = Mathf.Cos(ry / 2);
        var sy = Mathf.Sin(ry / 2);
        var cz = Mathf.Cos(rz / 2);
        var sz = Mathf.Sin(rz / 2);

        this.u = cx * cy * cz + sx * sy * sz;
        this.x = sx * cy * cz - cx * sy * sz;
        this.y = cx * sy * cz + sx * cy * sz;
        this.z = cx * cy * sz - sx * sy * cz;
    }

    public MyQuaternion(Vector3 v) : this(v.x, v.y, v.z) { }
```

## 演算子オーバーライド
クォータニオン同士やクォータニオンとスカラーの間の演算を定義します。
加減算やスカラー倍なんかはやるだけなので省略して、ちょっとややこしいクォータニオン間の積だけ。
適当に調べれば出てきますが、↓こんな感じ。

```cs:MyQuaternion.cs
    public static MyQuaternion operator *(MyQuaternion l, MyQuaternion r)
    {
        return new MyQuaternion
        (
            l.u * r.u - l.x * r.x - l.y * r.y - l.z * r.z,
            l.u * r.x + l.x * r.u + l.y * r.z - l.z * r.y,
            l.u * r.y - l.x * r.z + l.y * r.u + l.z * r.x,
            l.u * r.z + l.x * r.y - l.y * r.x + l.z * r.u
        );
    }
```

# 投影
我々は四次元人ではないので、四次元の物体をそのまま知覚することはできません。どうにかして三次元に落とし込む必要がありますが、その手段の一つとして**投影**があります。三次元の物体に光を当てると二次元の影が得られるように、四次元の物体に光を当てると三次元の影を得ることができます。ここでは見栄えや$u$軸方向の情報保持も考えて**透視投影**をすることにしましょう。

簡単のため、光源は$u$軸上の点$\ {\bf f}=(f,0,0,0)$、投影先のスクリーン（我々が知覚できる空間）は$\ u=0$で固定します。
ある四次元空間上の点 ${\bf v}=(v_u,v_x,v_y,v_z)\ $の投影先の点を ${\bf v'}=(v'_u,v'_x,v'_y,v'_z)$とすると、ある実数$t$があって

```math
{\bf v'}=t{\bf v}+(1-t){\bf f}\ 
```

ですが、$v'_u=0$ なので、

```math
t=\frac{f}{f-v_u}\ 
```
が得られ、${\bf v'}$ が求まります。三次元→二次元の場合とやってることは変わりません。
MyQuaternion構造体のインスタンスメソッドにしておきましょう。

```cs:MyQuaternion.cs
    //光源(f,0,0,0)からスクリーンu=0への透視投影
    public Vector3 Project(float f)    
        => Vector3.Lerp(Vector3.zero, new Vector3(x, y, z), f / (f - u));
    
    //平行投影はu座標を消すだけ
    public Vector3 ProjectParallel() => new Vector3(x, y, z);
```

## 例：超立方体の透視投影
三次元空間において、原点を中心とした立方体は、$(\pm1,\pm1,\pm1)$に8頂点を置いて最近傍頂点間を結べば作成できます。
四次元空間において、原点を中心とした超立方体は、$(\pm1,\pm1,\pm1,\pm1)$に16頂点を置いて最近傍頂点間を結べば作成できます。距離は普通にユークリッドノルムです。簡単ですね。

$(\pm1,\pm1,\pm1,\pm1)$をそれぞれ透視投影 $(f=2)$ した位置にSphereを置き、頂点の接続をLineRendererで示したのが以下の図です。[Wikipediaの図](https://ja.wikipedia.org/wiki/%E6%AD%A3%E5%85%AB%E8%83%9E%E4%BD%93)とも一致しますね。

![image.png](https://qiita-image-store.s3.amazonaws.com/0/84499/ad97831b-c4b1-58f3-cafe-dfd1a40fb689.png)

# 回転
三次元空間における（原点を中心とする）回転は、$xy\ $平面、$yz\ $平面、$zx\ $平面に沿った回転に分解することができます。
四次元空間においてはここに$ux\ $平面、$uy\ $平面、$uz\ $平面が加わるので、６種類の回転を組み合わせて一つの回転を表現することになります。

詳しいことは[Wikipedia](https://en.wikipedia.org/wiki/Rotations_in_4-dimensional_Euclidean_space)あたりに譲るとして、結論だけ。

- 四次元空間の（原点を中心とする）回転はデュアルクォータニオン$(q_L, q_R)$で定義される。
- 現在座標がクォータニオン$q_0$で表されているとき、回転後の座標は $q=q_Lq_0q_R$ で表される。

簡単ですね。MyQuaternion構造体のインスタンスメソッドにしておきます。

```cs:MyQuaternion.cs
    //回転（原点中心）
    public MyQuaternion Rotate(MyQuaternion QL, MyQuaternion QR) => QL * this * QR;
```

ではどう６軸の回転角を指定するかですが、ここで先程定義しておいたオイラー角からの変換が使えます。

```cs:MyQuaternion.cs
    //６つのオイラー角からデュアルクォータニオンを作って適用
    public MyQuaternion Rotate(float x1, float y1, float z1, float x2, float y2, float z2)
        => new MyQuaternion(x1, y1, z1) * this * new MyQuaternion(x2, y2, z2);
```

これを使って先程の超立方体を回転させてみたのがこちら。
<blockquote class="twitter-tweet" data-conversation="none" data-lang="ja"><p lang="ja" dir="ltr">基本的な回転 <a href="https://t.co/FnXlrEVVuK">pic.twitter.com/FnXlrEVVuK</a></p>&mdash; りじちょー (@Rijicho_nl) <a href="https://twitter.com/Rijicho_nl/status/1098325177367113728?ref_src=twsrc%5Etfw">2019年2月20日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

こんな感じで基本的な回転は比較的わかりやすいですが、複数の軸を組み合わせると冒頭に示したようなわけのわからないものになります。

# 切断
四次元のオブジェクトを三次元に落とし込む別の方法として、**超平面による切断**があります。
三次元のオブジェクトを平面で切断すると二次元の断面図が得られるのと同じように、四次元のオブジェクトを超平面で切断すると三次元の断面図が得られるわけです。

三次元の物体をある平面で切断する場合、物体の各面と切断面との交線で囲まれた領域（二次元）が断面となります。
四次元の物体をある超平面で切断する場合、物体の各胞と切断面との交面で囲まれた領域（三次元）が断面（断胞？）となります。
このとき、交面は物体の各面と切断面の交線で囲まれた領域（二次元）なので、結局やることは三次元のときと同じで、**「（超）平面と物体の各面との交線を求める」**です。

ここでは三角メッシュを総なめしていく手法で実装します（もっといい実装あったら教えてください）。

1. 切断面を$({\bf v}-{\bf c})\cdot {\bf n}=0\ $（${\bf c}$は面上の一点、${\bf n}$は法線）で表される平面とします。
2. ３つの頂点座標${\bf p}[i]\ (i=0,1,2)$で定義された三角メッシュについて、$({\bf p}[i]-{\bf c})\cdot {\bf n}\ $の正負を調べることで、各頂点が平面のどちら側にあるかを判定します。
3. 全ての頂点が同じ側にある場合、その三角メッシュは切断されません。
4. それ以外の場合、頂点は１個と２個に分けられます。１個の方を ${\bf p_0}$、２個の方を ${\bf p_1}, {\bf p_2}\ $とします。
5. ${\bf p_0}, {\bf p_1}\ $及び${\bf p_0}, {\bf p_2}\ $それぞれの組について、切断面上の内分点${\bf p_{01}}, {\bf p_{02}}\ $を求めます。
6. ${\bf p_{01}}, {\bf p_{02}}\ $を結ぶ線分が、切断面と三角メッシュの交線です。
7. 全ての三角メッシュについて2～6を繰り返して交線全体を得ます。

ここで、太字で示している変数はベクトルですが、三次元でも四次元でも全く同じアルゴリズムを適用できます。
結果として現れる線分の集合は、三次元であれば二次元領域を、四次元であれば三次元領域を囲っているはずです。

切断面を$u=0\ $で固定した場合、頂点のグループ判定は$u$成分の正負で決まります。内分点の方は

```math
{\bf p_{01}}=t_1{\bf p_0}+(1-t_1){\bf p_1}\\
{\bf p_{02}}=t_2{\bf p_0}+(1-t_2){\bf p_2}
```

ですが、${\bf p_{01}},{\bf p_{02}}$の$u$成分は0なので、投影のときと同じように$t_1, t_2$を求めることができます。

というわけで実装したのがこちら。計算された交線を全てLineRendererで出力しているだけですが、常に破綻のない三次元領域を示しているのが分かると思います。

<blockquote class="twitter-tweet" data-conversation="none" data-lang="ja"><p lang="ja" dir="ltr">回転する超立方体の切断 <a href="https://t.co/cUAFl66YFk">pic.twitter.com/cUAFl66YFk</a></p>&mdash; りじちょー (@Rijicho_nl) <a href="https://twitter.com/Rijicho_nl/status/1098339221230649344?ref_src=twsrc%5Etfw">2019年2月20日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>


ナイーブ実装ですが一応コードも。

```cs
    MyQuaternion[] CutAtU0(MyQuaternion A, MyQuaternion B, MyQuaternion C)
    {
        //超立方体の面を構成する各三角形ABCについて、
        //超平面u=0との交線を求める。
        
        //すべての頂点が同じ側にある場合は交わらない。
        if (A.u < 0 && B.u < 0 && C.u < 0 || A.u > 0 && B.u > 0 && C.u > 0)
        {
            return null;
        }
        //AとBが同じ側にある場合、AC間とBC間がu=0で切断される。
        else if (A.u * B.u >= 0)
        {
            //AC上、u=0の点Pを求める。
            var t = -C.u / (A.u - C.u);
            var P = t * A + (1 - t) * C;
            //BC上、u=0の点Qを求める。
            var s = -C.u / (B.u - C.u);
            var Q = s * B + (1 - s) * C;
            //PQは求めたい交線である。
            return new[] { P, Q };
        }
        //他の場合も同様
        else if (A.u * C.u >= 0)
        {
            var t = -B.u / (A.u - B.u);
            var P = t * A + (1 - t) * B;
            var s = -B.u / (C.u - B.u);
            var Q = s * C + (1 - s) * B;
            return new[] { P, Q };
        }
        else
        {
            var t = -A.u / (B.u - A.u);
            var P = t * B + (1 - t) * A;
            var s = -A.u / (C.u - A.u);
            var Q = s * C + (1 - s) * A;
            return new[] { P, Q };
        }
    }
```
# おわりに
で、これが何に使えるのかという話ですが。
今サークルで作っているゲームのグラフィックに取り入れてみようかなと考えてたり。