Unityシェーダで時間変化するマンデルブロ集合を描く

#概要

　先日、VRChatで以下のようなワールドを作成しPublic化しました。シェーダのソースコードを解説します。

<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr">「Fractal Room」Public化しました！　こんな感じで、拡大・縮小を続ける巨大なマンデルブロ集合を眺められるワールドです。全部シェーダで書いてるので、ワールドの容量は正真正銘の0MB！遊びに来てね～　<a href="https://twitter.com/hashtag/VRC?src=hash&amp;ref_src=twsrc%5Etfw">#VRC</a> <a href="https://twitter.com/hashtag/VRChat_world%E7%B4%B9%E4%BB%8B?src=hash&amp;ref_src=twsrc%5Etfw">#VRChat_world紹介</a> <a href="https://t.co/DGL4FSsmIM">pic.twitter.com/DGL4FSsmIM</a></p>&mdash; Hibit (@hibit_at) <a href="https://twitter.com/hibit_at/status/1043301724105011200?ref_src=twsrc%5Etfw">2018年9月22日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

#一部引用の明記

　マンデルブロ集合の計算については、[この記事](http://wordpress.notargs.com/blog/blog/2015/11/07/unity%E3%83%9E%E3%83%B3%E3%83%87%E3%83%AB%E3%83%96%E3%83%AD%E9%9B%86%E5%90%88%E3%82%92%E6%8F%8F%E3%81%8F%E3%82%B7%E3%82%A7%E3%83%BC%E3%83%80%E3%82%92%E4%BD%9C%E3%82%8B/)を参考にさせていただきました（@notargsさん、ありがとうございます）。本コードではこれの一部を改変して、

* 時間変化によって任意の座標を中心に拡大・縮小する
* 発散速度によってカラフルなグラデーションをつける

　という要素を付け加えています。

#ソースコード

　とりあえずソースを見せろという方へ。

```GLSL
Shader "Custom/Mandelbrot"
{
	Properties{
		_MaxIteration("MaxIteration", Range(1, 2048)) = 512
		_Threshold("Threshold", Range(1, 100)) = 2
		_PosX("PosX",Range(-2,2)) = -0.685
		_PosY("PosY",Range(-2,2)) = -0.3
	}
		SubShader
	{
		Tags { "RenderType" = "Opaque" }
		LOD 100
 
		Pass
		{
			CGPROGRAM
			#pragma vertex vert
			#pragma fragment frag
 
			float _Threshold;
			int _MaxIteration;
			float _PosX;
			float _PosY;

			float2 cpow(float2 v)
			{
				float a = 3.14;
				if (v.x != 0){
					a = atan2(v.y, v.x) * 2;
					}
				return float2(cos(a), sin(a)) * pow(length(v), 2);
			}
 
			half mandelbrot(half2 c)
			{
				half2 z = half2(0, 0);
				for (int i = 0; i < _MaxIteration; i++)
				{
					z = cpow(z);
					z += c;
					if (length(z) > _Threshold){
					return (float)i / _MaxIteration;
					}
				}
				return 1.0;
			}

			struct appdata
			{
				float4 vertex : POSITION;
				float2 uv : TEXCOORD0;
			};
 
			struct v2f
			{
				float2 uv : TEXCOORD0;
				float4 vertex : SV_POSITION;
			};
 
			v2f vert (appdata v)
			{
				v2f o;
				o.vertex = UnityObjectToClipPos(v.vertex);
				o.uv = v.uv;
				return o;
			}
			
			fixed4 frag (v2f i) : SV_Target
			{
				float scale = 2.1+pow(sin(_Time.x),5)*2;
				scale = pow(10,scale)-1;
				i.uv.x = i.uv.x + _PosX * scale;
				i.uv.y = i.uv.y + _PosY * scale;
				fixed col = mandelbrot((i.uv - 0.5)/scale);
				fixed col_r = 0.5+sin(col*6.28)/2;
				fixed col_g = 0.5+sin(col*6.28+2.09)/2;
				fixed col_b = 0.5+sin(col*6.28+4.19)/2;
				if (col == 1){
					col_r = 0;
					col_g = 0;
					col_b = 0;
				}
				return fixed4(col_r,col_g,col_b,0);
			}
			ENDCG
		}
	}
	Fallback "Diffuse"
}
```

#複素数の2乗を定義

　マンデルブロ集合の計算には、複素数の2乗を計算する必要がありますが、その計算関数を自前で用意します。ある複素数$z = x + yi$が複素数平面の極座標で$z(r,\theta)$として表せるとき、

```math
z(r,\theta)^2 = z(r^2,2\theta)
```

　が成り立ちます。$r$は$\sqrt{x^2+y^2}^{※1}$、$\theta$は$tan^{-1}(\frac{y}{x})^{※2}$で求められるので、コードとしては以下のようになります。

```GLSL
			float2 cpow(float2 v)
			{
				float a = 3.14;
				if (v.x != 0){
					a = atan2(v.y, v.x) * 2;
					}
				return float2(cos(a), sin(a)) * pow(length(v), 2);
			}
```

※1　…　この計算はlength関数としてデフォルトで用意されています。
※2　…　x = 0の場合は0除算でエラーになってしまうのでif文で個別処理します。

#マンデルブロ集合を計算

　ある複素数$c$について、

```math
z_0 = 0\\
z_{n+1} = z_n^2 + c
```

　という計算をn回おこなった時、$z_n$の長さが一定以上になる時を「発散」とします。規定の繰り返しまでに「発散」すればその回数の逆数（≒　1をMAXとした発散しづらさ）を、しなければあきらめて1を返すような関数Mandelbrotを定義します。これによってマンデルブロ集合が描かれます。何回計算して（MaxIteration）何を閾値とするか（Threshold）はプロパティで調整しますが、それぞれ512、2ぐらいがちょうど良いかなと思います。

```GLSL
			half mandelbrot(half2 c)
			{
				half2 z = half2(0, 0);
				for (int i = 0; i < _MaxIteration; i++)
				{
					z = cpow(z);
					z += c;
					if (length(z) > _Threshold){
					return (float)i / _MaxIteration;
					}
				}
				return 1.0;
			}
```

#時間によって拡大・縮小 + 色付け

```GLSL
			fixed4 frag (v2f i) : SV_Target
			{
				float scale = 2.1+pow(sin(_Time.x),5)*2;
				scale = pow(10,scale)-1;
				i.uv.x = i.uv.x + _PosX * scale;
				i.uv.y = i.uv.y + _PosY * scale;
				fixed col = mandelbrot((i.uv - 0.5)/scale);
				fixed col_r = 0.5+sin(col*6.28)/2;
				fixed col_g = 0.5+sin(col*6.28+2.09)/2;
				fixed col_b = 0.5+sin(col*6.28+4.19)/2;
				if (col == 1){
					col_r = 0;
					col_g = 0;
					col_b = 0;
				}
				return fixed4(col_r,col_g,col_b,0);
			}
```

　この部分で

* 時間変化による、任意の座標を中心とした拡大・縮小
* 発散速度を元にしたグラデーション色付け

　を行っています。
## 拡大・縮小

```GLSL
				float scale = 2.1+pow(sin(_Time.y),5)*2;
				scale = pow(10,scale)-1;
				i.uv.x = i.uv.x + _PosX * scale;
				i.uv.y = i.uv.y + _PosY * scale;
				fixed col = mandelbrot((i.uv - 0.5)/scale);
```

　倍率（scale）は、1倍から10,000倍ぐらいまでが周期変化するようにしたいのですが、線形的に変化させると前半が速くなりすぎる&後半が遅くなりすぎるので、指数的に変化するようにします。$10^n$の右肩にある指数部を、三角関数を用いて0.1～4.1で周期変化するようにしています。なぜ+0.1しているのかについては、0にすると$10^0-1 = 0$になってエラーが起こるからです。

　拡大・縮小の中心となる座標については、プロパティ（PosXとPosY)で調整します。縮小した時はともかく、拡大した時は0.01動くだけですごく動いてしまうので、ミリ単位での調整が必要になってきます。私はPosX = -0.685,PosY = -0.3がいい感じなのを見つけたのでそれにしてますが、各自エモい組み合わせを見つけてください。

## グラデーション

```GLSL
				fixed col_r = 0.5+sin(col*6.28)/2;
				fixed col_g = 0.5+sin(col*6.28+2.09)/2;
				fixed col_b = 0.5+sin(col*6.28+4.19)/2;
				if (col == 1){
					col_r = 0;
					col_g = 0;
					col_b = 0;
				}
				return fixed4(col_r,col_g,col_b,0);
```

　先程、colという変数に各座標におけるMandelbrot関数を適用した後の値をいれました。これはその座標の「発散しづらさ」を示しています。0だとすぐに発散して、何回計算しても発散しなかったら1になります。

　この0～1の間の数値にうまく色をあてはめてやればグラデーションができます。このコードではRGBの3色について三角関数を(2/3)πずつずらして足し合わせることによりそれを表現しています。また、発散しなかった場合はその領域をハッキリと明示したいので（巷でよくみるマンデルブロ集合の絵もよくそうなっています）そこだけはif文の個別処理で黒色にしています。

　ここの味付けも各自調節して、最も美しいと感じる組み合わせを見つけるのもいいと思います。本記事ではソースコードまで貼りませんが、グラデーション用のテクスチャ画像して、それを読み込むようにすればより自由にグラデーションを定義できます。