レイマーチングでメンガーのスポンジを描く

前に書きかけだったやつが、いちおうそれなりに形になっていたので公開しちゃいますｗ

Menger sponge（メンガーのスポンジ）を描く

Wikipediaから画像を引用させてもらうとこんな感じの絵。

また、文章では以下のように説明されています。

メンガーのスポンジとは自己相似なフラクタル図形の一種であり、立方体に穴をあけたものである。そのフラクタル次元（ハウスドルフ次元、相似次元）は $\frac {\log 20}{\log 3}(=2.7268\ldots )$ $\frac{\log 20}{\log 3}(=2.7268\ldots )$次元である。メンガーのスポンジの面は同じくフラクタル図形のシェルピンスキーのカーペットでできている。

それをレイマーチングで描き出し、かつきれいな絵にしている作品がこちら。

- Visualization Menger sponge

ちなみに前回は[GLSL] Shadertoyのシェーダ芸人になるためのTips集という記事も書いたのでよかったら見てみてください。

さて、上記作品のメンガーのスポンジはどうやって描いているのか。それを紐解いてみようと思います。

大枠としての処理はおそらく以下のようなフローで行っていると考えています。

1. 指定サイズ空間での繰り返し（modなどを使うリピートパターン）
2. その指定サイズ空間内のXY平面、YZ平面、XZ平面それぞれで一番遠い位置を3点探し、かつそれらの中で一番近い位置を探す（距離関数）
3. (2)で計算したものよりやや小さい空間で同様の計算を行う（ここは複雑さに応じて複数回実行）
4. 最後に、(2)と(3)の計算結果をmax関数による合成を行う

という手順でフラクタル構造を作っていると考えています。

これを踏まえて、上記フローの(2)と(3)の部分を時間によってモーフィングしてみると以下のような絵が得られます。
小さい空間と大きい空間が複雑にモーフィングしているのが分かるかと思います。

イメージ的には、各空間ごとにmax関数で距離を採択しているために、ひとつ分であれば通路のような見た目になるものの、それが複数組み合わさることで距離の算出に差異が出てこうした図が描かれる、という感じでしょうか。（ちょっとまだ頭の中に正確なイメージが湧いていない・・・）

ちなみにこれを、モーフィングではなくmax関数によって合成すると以下のような絵が得られます。

max関数での合成なので、重なり合っている部分だけ描かれることになります。

（詳細についてはwgld.orgで解説されているのでそちらを見るといいでしょう→ オブジェクトの重なりを考慮した描画 | wgld.org）

モーフィングの最初と最後だけだとこの形は想像できないですが、おそらく、遮蔽されてしまっている「壁の向こう」ではもう少し複雑な形状になっていることに気づいていないだけかな、と思っています。

これの距離関数のコードは以下です。

float s = 10.0;
vec3 p = abs(fract(q / s) * s - s * 0.5); // Repeat space.

// Get maximum distance of surface of cube?
float a1 = max(p.x, p.y); // XY平面で一番遠い位置
float b1 = max(p.y, p.z); // YZ平面で一番遠い位置
float c1 = max(p.x, p.z); // XZ平面で一番遠い位置

float size = 1.0;

// さらにその中で一番近い位置を計算。そこから`size`を引くことで、Void Cubeの距離関数として機能する
float d1 = min(a1, min(b1, c1)) - size + 0.05; // Void Cube.

// --------------------------------------------------------

// もうひとつ、少し小さめのVoid Cubeを計算し、それを最後合成する
s = 1.0;
// 上記のfractを使った計算とやっていることは同じ。2パターン明示。
p = abs(mod(q, s) - s * 0.5); // Repeat space.

// Get maximum distance of surface of cube?
float a2 = max(p.x, p.y);
float b2 = max(p.y, p.z);
float c2 = max(p.x, p.z);

size = 1.0 / 3.0;

float d2 = min(a2, min(b2, c2)) - size + 0.05;

float d = max(d1, d2);

float m = abs(sin(iTime * 0.1));
float r = mix(d1, d2, m);

return r;