GLSLでIFSシェルピンスキーのカーペット

最近、GLSLでIFSフラクタルを描く勉強しています。

手始めにMneger Spongeと思ったのですが、3Dよりも2Dのほうが簡単だろうということでMenger Spongeの2D版といった感じのシェルペンスキーのカーペットを考えてみました。

<img width="421" alt="sherpenski-carpet.png" src="https://qiita-image-store.s3.amazonaws.com/0/61569/b8d0e365-0a5b-28d5-7411-161b2671866f.png">

実際にやろうするとどのように空間を操作するかが難しく、半日ぐらい悩んだので備忘録として記事にしました。

自分で考えた方法とiq氏が[こちらの記事](ttps://www.iquilezles.org/www/articles/menger/menger.htm)でMenger Spongeについて解説している方法を参考にしたものの二通りを紹介します。
どちらを使用しても結果は同じになるはずです(たぶん)。

## 自分で考えた方法

ブーリアン演算で繰り返し差を取ることでシェルピンスキーのカーペットを作成します。
実際のコードと動くものは[こちら](http://glslsandbox.com/e#51001.3)で見られます。


```glsl
float sierpinskiCarpet(vec2 p) {
	float s = 1.0;
	float f  = 1.0 / 3.0;
	float d = sdBox(p, vec2(s));
	for (int i = 0; i < 4; i++) {
		d = max(d, -sdBox(p, vec2(f)) / s);
		p = abs(p);
		p -= f;
		p = abs(p);
		p -= f;
		p *= 3.0;
		s *= 3.0;	
	}
	return d;	
}
```

まず、`float d = sdBox(p, vec2(s));`を正方形を描きます。
そこから`d = max(d, -sdBox(p, vec2(f)) / s);`で論理差を取り、真ん中に穴を開けます。

![method1_01.png](https://qiita-image-store.s3.amazonaws.com/0/61569/42ace6ee-ca2b-a865-cb8b-6415dd48904a.png)

次に周囲8近傍についても穴を開ける必要がありますが、X軸、Y軸それぞれについて対称なので`abs(p);`で空間を折りたたみます。これにより$x \geqq 0$かつ$y \geqq 0$の座標だけを考えればよくなります。

![method1_02.png](https://qiita-image-store.s3.amazonaws.com/0/61569/183e2d0c-d55c-4a2c-11d6-65e7f9c3d094.png)

折りたたんだものを見ると、X軸、Y軸それぞれについて$(x, y)=(0.333..., 0.333...)$で（左下を無視すれば）対称になので、`p -= f;`で$(x, y)=(0.333..., 0.333...)$を原点に移動させてから`p = abs(p);`で空間をもう一度折りたたみます。

![method1_03.png](https://qiita-image-store.s3.amazonaws.com/0/61569/c4b75ecf-f880-d590-aef8-18d95187b828.png)

折りたたんだものを`p -= f;`で穴が真ん中になるように座標を移動させます。

![method1_04.png](https://qiita-image-store.s3.amazonaws.com/0/61569/eebe291e-b316-b9b4-8c3c-a67eefd9e0df.png)

その結果、以下のようになり、`p *= 3.0;`することで最初の画像と同じ状態に戻るので、同じ操作を繰り返し穴を開けていくことができます。

![method1_05.png](https://qiita-image-store.s3.amazonaws.com/0/61569/9f30aed1-8645-9e26-c247-c5b352960e31.png)



## iq氏のMenger Spongeを参考にした方法

ブーリアン演算で繰り返し積を取ることでシェルピンスキーのカーペットを作成します。
コードと実際に動くものは[こちら](http://glslsandbox.com/e#51038.0)でみられます。

```glsl

float sdCross(vec2 p) {
	p = abs(p);
	return min(p.x, p.y) - 1.0;
}

float sierpinskiCarpet(vec2 p) {
	float d = sdBox(p, vec2(1.0));
	
	float s = 1.0;
	for (int i = 0; i < 4; i++) {
		vec2 a = mod(p * s, 2.0) - 1.0;
		vec2 r = 1.0 - 3.0 * abs(a);
		s *= 3.0;
		float c = sdCross(r) / s;
		d = max(d, c);
	}
	return d;
}
```

最初に`float d = sdBox(p, vec2(1.0));`で四角形を描くのは先に述べた手法と同じです。

`vec2 a = mod(p * s, 2.0) - 1.0;`で座標を操作をすると赤で示したようになります。

![method2_01.png](https://qiita-image-store.s3.amazonaws.com/0/61569/27871c5e-b159-21d1-4ddd-6b137ff85a2b.png)

さらに`vec2 r = 1.0 - 3.0 * abs(a);`で座標は次のようになります。

![method2_02.png](https://qiita-image-store.s3.amazonaws.com/0/61569/b7e00cb5-323d-776d-396c-15bb42845f58.png)

`float c = sdCross(r) / s;`は以下のような形となっており、` d = max(d, c);`で論理積を取ることで穴以外の部分を残すことができます。

![method2_03.png](https://qiita-image-store.s3.amazonaws.com/0/61569/6db4e114-8239-b0b9-82de-61411c92249a.png)

##　余談

[iq氏のMenger Spongeの記事](https://www.iquilezles.org/www/articles/menger/menger.htm)は大変参考になるのですが、最初のコードでは論理差を取るのに、繰り返しが導入されてからは論理積をとるようになるので記事を読み進める上で注意が必要です。

http://glslsandbox.com/e#51032.0

<img width="513" alt="menger-sponge.png" src="https://qiita-image-store.s3.amazonaws.com/0/61569/8a52413a-28fe-2a65-b603-10074be4a96d.png">