[コードリーディング vol.3] Fire Storm Cube

## 概要

今回はShadertoyで見かけたこちらの作品のコードリーディングをしてみようと思います。

![Fire Storm Cube](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/82665e0e-f88e-eb85-6070-23d971e8f5f7.png)
- [Fire Storm Cube](https://www.shadertoy.com/view/ldyyWm)

ちなみに、これを解読してちょっとコード整理しつつ音に反応して見るようにしてみました↓
![Burn Sound Wave](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/5415b4a4-ba14-0e7d-517a-8ba14eb42493.gif)
- [Burn Sound Wave](https://www.shadertoy.com/view/WlB3zc)


## 解説

まずはコード全体を。そんなに長くありません。

```glsl:コード全体
float burn;

mat2 rot(float a)
{
    float s = sin(a);
    float c = cos(a);
    return mat2(s, c, -c, s);
}

float map(vec3 p)
{
    vec3 p1 = abs(p);
    
    // ========================================
    // Render the cube.
    float cubeSize = 0.5;
    float d1 = max(max(p1.x, p1.y), p1.z) - cubeSize;
    
    // ========================================
    // Rotate grids.
    mat2 rm = rot(-iTime / 3.0 + length(p));
    p.xy *= rm;
    p.zy *= rm;
    
    // PingPon on 0.5
    vec3 q = abs(p) - iTime;
    q = abs(q - round(q));
    
    rm = rot(iTime);
    q.xy *= rm;
    q.xz *= rm;
    
    // ========================================
    // Render grids space.
    float l1 = length(q.xy); // Plane XY
    float l2 = length(q.yz); // Plane YZ
    float l3 = length(q.xz); // Plane XZ
    
    float gridOffs = 0.01;
    float d2 = min(min(l1, l2), l3) + gridOffs;
    
    burn = pow(d2 - d1, 2.0);
    
    return min(d1, d2);
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (2.0 * fragCoord - iResolution.xy);
    vec3 rd = normalize(vec3(uv, iResolution.y));
    vec3 ro = vec3(0, 0, -2);
    
    // Rotate ray position and ray target.
    mat2 r1 = rot(iTime / 4.0);
    rd.xz *= r1;
    ro.xz *= r1;
    
    mat2 r2 = rot(iTime / 2.0);
    rd.yz *= r2;
    ro.yz *= r2;
    
    float t = 0.0;
    float i = 30.0 * (1.0 - exp(-0.2 * iTime - 0.1));

    for(; i-- > 0.0;)
    {
        t += map(ro + rd * t) / 2.0;
    }
    
    fragColor = vec4(1.0 - burn, exp(-t), 2.0 * exp(-t / 4.0 - 1.0), 1.0);
}
```

理解を促すために最初にポイントを書いておくと、

1. Cubeの（ちょっとトリッキーな）距離関数でCubeを描く
2. 軸に沿ったラインを求める距離関数を視点を回転させてエモくする
3. それを合成する

というような流れです。

では、まずは`mainImage`関数から順番に見ていきましょう。
（念のため補足しておくとShadertoyにおいてこのメソッドがエントリーポイントになります）



## メイン関数冒頭のrotは視点の回転

冒頭で以下のようにベクトルを回転させています。

```glsl:視点の回転
mat2 r1 = rot(iTime / 4.0);
rd.xz *= r1;
ro.xz *= r1;
    
mat2 r2 = rot(iTime / 2.0);
rd.yz *= r2;
ro.yz *= r2;
```

これは、レイの位置と注視点を回転させています。
これをコメントアウトしてもエフェクトのイメージはそこまで変わりません。
大きく変化が出るのは中央のCubeが回る部分ですね。
内容自体はシンプルです。



## レイマーチの回数は時間で急激に増加するが、一定値以上はほぼ横ばい

さて、次はレイマーチのループ部分です。
レイマーチの回数`i`の部分が、ちょっと複雑な計算によって行われています。

```glsl:マーチ回数の計算
float i = 30.0 * (1.0 - exp(-0.2 * iTime - 0.1));
```

一瞬、なにをやっているのか分かりづらいですが、これをグラフにすると見えてきます。
グラフ化したのが以下です↓
![マーチ回数のグラフ化](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/89d17ff5-c91a-c96b-c495-681fab2a8b8d.png)

これを見ると、`x`（つまり`iTime`）の増加によって、最初は急激に上昇していますが、一定値を超えるとほぼ横ばいとなっているのが分かります。
計算式の中で`exp`を使っていることからもなんとなくイメージできるかと思います。

要は、冒頭と後半で見た目を大きく変えている、と考えるといいかと思います。
なので、解読後に自分で作った方はループ回数を固定にしていますが、見た目の変化はほとんどありません。


## 色の決定

ループ内の`map`関数はあとに回すとして、最後の色を決定する部分を見てみましょう。

```glsl:色の決定
fragColor = vec4(1.0 - burn, exp(-t), 2.0 * exp(-t / 4.0 - 1.0), 1.0);
```

`burn`変数は関数外で定義されています。`map`関数内で更新される想定のものになっています。
`burn`が使われているのは`r`要素だけなので、つまりはレイマーチの距離に応じて赤みが変化する、という感じですね。
残りの`gb`要素はレイマーチの距離結果（`t`）をそのまま使い、`exp`関数によって値が変化するようになっています。

これをグラフにすると以下のようになります。
![expを使った値の変化](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/33f3d5d6-959c-c347-35c1-f0511f659723.png)
グラフを見ると分かりますが、`t`の値が小さいほど大きく、逆に大きくなるほど小さくなっていきます。
つまりは画面に近いあたりの部分は青くなる、というようなイメージですね。

ちなみに白黒だけにするとこんな感じ。これはこれでかっこいいですね。
![モノクロ版](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/40040f0c-62b1-2032-b059-769c86abfaa1.png)


## map関数の中を見る

さて、次からはいよいよ本丸、`map`関数の中を見ていきましょう。



## ちょっと特殊なBoxの距離関数

`map`関数の冒頭で行われている計算は一瞬難解な感じですがただのBoxの距離関数です。

```glsl:Planeを組み合わせたBoxの距離関数
float d = max(max(abs(p.x), abs(p.y)), abs(p.z)) - 0.5;
```

そして最後の`-0.5`がBoxのサイズですね。
これを少し分かりやすく整形すると以下のようになります。

```glsl:整形コード
vec3 p1 = abs(p);
float d = max(max(p.x, p.y), p.z) - 0.5;
```

`p`の絶対値を取り、さらに`xyz`要素の中で最大のものを採用します。
これはつまり、**原点**から各軸の距離を測っていることになります。
そして`-0.5`とすることで、各軸が`0.5`だけオフセットします。

理解を簡単にするためにXY平面に限って考えてみると分かりやすいです。
XY軸上で、原点からの距離を測り、さらにそこから`-0.5`すると、以下のようになりますね。
![平面との距離のイメージ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/12a091c4-b9e8-182d-2408-97c1a62dacec.jpeg)

つまり、XY平面上でいうと赤い線と緑の線が地平線となっていて、まだベクトルが小さいときはマイナスになります。（＝地平線の下にいる）
それがある程度進むと地平線の上に出て、あとは普通に平面との距離計算になる、というわけですね。

そしてそれを`max`による`AND`合成することでBoxが描かれる、というわけです。
イメージ的には上の図で書いた方向に法線を持った3平面を合成するとBoxが現れる、というわけです。
（`abs`によって折りたたんでいるため、結果、全方位を平面が囲んでいるような状態になります）

これを実際に利用したレイマーチングの結果は以下です↓
![Planeを合成したCube](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/3829ee9d-b8c7-ca9b-2037-a26b7fe54184.png)
- [Cube by planes combining](https://www.shadertoy.com/view/WtSGzc)



## 0.0～0.5～0.0で往復

これを見たときはちょっと驚きました。
だいぶスマートに`0.0～0.5～0.0`を往復させています。

コードは以下。

```glsl:0.0～0.5～0.0の範囲で往復
vec3 q = abs(p) - iTime;
q = abs(q - round(q));
```

こちらの計算。`round`は四捨五入をする関数です。
つまり、`q`の値が`0.4`までは`0.4 - 0`となって`0.4`となります。
一方、`0.5`を超えると`0.55 - 1.0`となって、`0.45`（`abs`なので絶対値）となり、値が大きくなるにつれて今度は`0`に近づいて行くことになります。（`0.9`のときは`0.9 - 1.0 = 0.1`）

また、`q`の値が極端に大きくても相殺されて常に`0.5`の範囲にとどまる、というのもポイントです。
この計算の意図は、レイの位置`p`を折りたたんでいる、と考えるといいと思います。



## XYZ軸それぞれへの距離の合成でグリッドを表現

次はグリッド表示の部分です。（実際にはこれを歪ませているのでグリッドには感じませんが）

次のコードを見てください。

```glsl:XYZ軸のライン描画（＝グリッドのライン描画）
// PingPon on 0.5
vec3 q = abs(p) - iTime;
q = abs(q - round(q));
    
rm = rot(iTime);
q.xy *= rm;
q.xz *= rm;
    
float l1 = length(q.xy); // Plane XY
float l2 = length(q.yz); // Plane YZ
float l3 = length(q.xz); // Plane XZ
    
float d2 = min(d1, min(min(l1, l2), l3) + .01);
```

冒頭の部分は前述の`0.5`で往復する、という前述の部分ですね。
これによってレイの位置を折りたたみます。

その下の回転は視点回転なので飛ばして、`l1〜l3`の変数部分に注目してください。

それぞれ、`XY`平面、`YZ`平面、`XZ`平面と３つの面に限定したベクトルの長さを測っています。
これは（おそらく）基底の軸との距離関数だと思われます。

試しにXY平面のベクトルについて考えてみましょう。
3次元ベクトルのXY要素のみ、ということはつまり、そのベクトルをXY平面に投影したもの、と考えることができます。
つまり、

```glsl:平面への投影のイメージ
vec3 p = vec3(1, 1, 1); // というベクトルがあったとして
p.z = 0.0; // と、Zの値を0にしているのに等しい
```

さて、ではこの「**XY平面のベクトル**」の長さはなにとの距離なのか。
それは「**Z軸**」との距離になります。
なぜなら、XY平面はZ軸に垂直だからです。なのでそのままベクトルの長さがZ軸との距離になる、というわけですね。

そしてそれ以外の平面についても考え方は同様です。
最後にそれを`min`によって合成しています。
ちなみに最後の`0.01`は線の太さを表しています。

なのでなんとなく線が引かれるイメージが湧いてきませんか？
加えて、評価しているベクトル`q`は`abs`により折り畳まれているのと、`round`によるトリックで`0〜0.5〜0`を行き来しています。
つまり「繰り返し」が現れているということですね。

ということで、この表現部分だけを抜き出して視覚化したのが以下のものです。

![グリッド表現](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/75f59a44-e082-bd24-8868-464deabc10f8.png)

繰り返しによって無限のグリッドが表示されているのが分かるかと思います。
非常にシンプルなのにこうした絵が出てくるのは本当に面白いですね。

ちなみにこちらも[Shadertoy](https://www.shadertoy.com/view/3sBXDy)にも上げてあります。


### rot関数によって視点を曲げる

グリッドが理解できたらあと少しです。
このグリッド上の表現を`rot`関数によって曲げることで、冒頭のようなカッコいいエフェクトに変化させている、というわけなんですね。

```glsl:rotによる視点回転
// ========================================
// Rotate grids.
mat2 rm = rot(-iTime / 3.0 + length(p));
p.xy *= rm;
p.zy *= rm;

// PingPon on 0.5
vec3 q = abs(p) - iTime;
q = abs(q - round(q));

rm = rot(iTime);
q.xy *= rm;
q.xz *= rm;
```

ちなみにこの回転をオフにしてみると以下のような絵になります。
![視点回転オフ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/a424d8f4-d44e-235f-f5e0-0d98be656504.png)
なんとなくグリッドがうっすらと見えているのが分かるかと思います。

同じタイミングで、視点の回転を入れると以下のような感じになります。
![視点回転オン](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/9752142c-cd16-c798-7be2-4a65f7f3f04b.png)

ぐにゃっと曲がっているのがなんとなく分かるかと思います。

この歪みを表現しているコツは以下の部分です。

```glsl:歪みのコツ
mat2 rm = rot(-iTime / 3.0 + length(p));
```

`rot`関数の引数に渡している値計算のうち、`length(p)`がそのポイントです。
これは、レイの位置に応じて変化する値になりますよね。
そして外に飛んでいけばいくほど値は大きくなります。（レイの進む距離が大きくなるから）
つまり、画面の際にいくほど大きく歪むように調整している、というわけです。


## まとめ

いかがだったでしょうか。
非常にシンプルなコードの中に学びがたくさんありました。
Shadertoyの作品はこの手のテクニックの宝庫なので、気に入った表現をこうして読み解いて行くととても学びが多いし、またそれを応用して自分自身でもカッコいい作品が作れるようになっていくのでオススメです。