Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

概要

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

一部引用の明記

 マンデルブロ集合の計算については、この記事を参考にさせていただきました(@notargsさん、ありがとうございます)。本コードではこれの一部を改変して、

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

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

ソースコード

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

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)$として表せるとき、

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

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

            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$について、

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

 という計算をn回おこなった時、$z_n$の長さが一定以上になる時を「発散」とします。規定の繰り返しまでに「発散」すればその回数の逆数(≒ 1をMAXとした発散しづらさ)を、しなければあきらめて1を返すような関数Mandelbrotを定義します。これによってマンデルブロ集合が描かれます。何回計算して(MaxIteration)何を閾値とするか(Threshold)はプロパティで調整しますが、それぞれ512、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;
            }

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

            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);
            }

 この部分で

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

 を行っています。

拡大・縮小

                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がいい感じなのを見つけたのでそれにしてますが、各自エモい組み合わせを見つけてください。

グラデーション

                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文の個別処理で黒色にしています。

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

hibit
数学とか3Dとか翻訳とか
https://deux-hibi.hatenablog.com/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away