f:id:aizu-vr:20200810180146g:plain

始めに

こんにちは、なおしです。

積読していた「作って動かすALife」を読もうとした際にちょうどブログ当番が回ってきたので、復習がてらComputeShaderを使って先ほどの本で最初に紹介されるGray-Scottモデルを実装してみようと思います。

作って動かすALife ―実装を通した人工生命モデル理論入門

作者:岡瑞起,池上高志,ドミニク・チェン,青木竜太,丸山典宏
発売日: 2018/07/28
メディア: 単行本（ソフトカバー）

ここではComputeShaderについては解説しませんので、もしわからなかった場合はこれらのサイトを参考にしてください。

neareal.com

また、このプロジェクトはgithubに上げています。

github.com

Gray-Scottモデルについて

書籍の内容を出来るだけ噛み砕いて解説しようと思います。間違いがあったり、話が飛躍してしまうかもしれませんがご了承下さい。

自己組織化

例えば雪の結晶や雲の形状、地面のひび割れなど完成形がわかっていないにも関わらずある構造を形成することがあります。このことを自己組織化といいます。

反応拡散系

空間に分布された一種あるいは複数種の物質の濃度が、物質がお互いに変化し合うような局所的な化学反応と空間全体に物質が広がる拡散の二つのプロセスの影響によって変化する様子を数理モデル化したものである。（wikipediaより引用）

これにより、自然界の様々なパターンを自己組織的なロジックで作れることが示されています。

Gray-Scottモデル

Gray-Scottモデルは反応拡散系モデルの一つで、空間内で2種類の物質UとVを使用して化学反応と拡散過程を表したモデルです。

この物質Uは一定の補充率で追加され、逆に物質Vは一定の減量率で除かれます。

Gray-Scottモデルでは、1つの物質Uと2つの物質Vが反応して3つの物質Vになります。また、一定の減量率で除かれた物質Vは不活性生成物Pになり、これ以上の反応をせず変化しない物質になります。

$U + 2V \rightarrow 3V \\ V \rightarrow P$

2種類の物質UとVは異なる拡散係数を持ち、物質Uは物質Vに比べて素早く空間全体に拡散します。

数式・解説

先程解説した物質Uと物質Vの変化量は次の数式で表されます。

Gray-Scottモデルの反応拡散方程式 $\dfrac{\partial u}{\partial t} = Du \Delta u - uv^2 + f(1 - u) \\ \dfrac{\partial v}{\partial t} = Dv \Delta v +uv^2 - v(f + k)$

$\dfrac{\partial u}{\partial t}$ 、 $\dfrac{\partial v}{\partial t}$ は「時間に対するU, Vの変化量」という意味です。

$Du$ 、 $Dv$ は拡散定数でUとVの拡散の速さを表します。

次に $\Delta u$ 、 $\Delta v$ について解説します。

今回想定している空間は2次元なので次のように表せます。

$\Delta u = \dfrac{\partial^2}{\partial x^2}u + \dfrac{\partial^2}{\partial y^2}u$

$\dfrac{\partial ^2}{\partial ^2 x}u$ はどういう意味かというと2階微分なのでx軸に対する変化量の変化量となります。

変化量の変化量と言われてもぱっとしないと思うので軽く説明します。

ある地点 $x$ の物質Uの濃度を $u(x)$ とします。

右隣の地点 $x+1$ との濃度の変化量は $u(x + 1) - u(x)$ となります。

同様にして左隣の地点 $x - 1$ との変化量は $u(x) - u(x - 1)$ です。

この2つの変化量の差を求めることで変化量の変化量を求めることが出来ます。

$\dfrac{\partial ^2}{\partial ^2 x}u = (u(x + 1) - u(x)) - (u(x) - u(x - 1)) = u(x + 1) + u(x - 1) - 2u(x)$

同様にしてy軸も考慮して以下のようになります。

$\Delta u = \dfrac{\partial^2}{\partial x^2}u + \dfrac{\partial^2}{\partial y^2}u \\ = u(x + 1, y) + u(x - 1, y) + u(x, y + 1) +u(x, y - 1) - 4u(x, y)$

次に $uv ^2$ について解説します。

先ほど解説したこの式がこの部分に該当します。

$U + 2V \rightarrow 3V \\ V \rightarrow P$

1つの物質Uと2つの物質Vが反応することにより物質Uは減るので $\dfrac{\partial u}{\partial t}$ では減算、逆に物質Vは増えるので $\dfrac{\partial v}{\partial t}$ では加算します。

最後に $f(1 - u)$ と $v(f + k)$ は物質Uの補充率、物質Vの補充率を表します。

物質Uではuの値が1に近づくにつれて値は0になります。物質Vの場合はvの値が0より大きい場合、常に減り続けます。

実装

実装については作って動かすALifeを参考に作成しています。

C#プログラム

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEngine.Serialization;

public class Simulator : MonoBehaviour
{
    private RenderTexture m_uvTexture;

    [SerializeField] 
    private ComputeShader m_grayScottCalculator;

    [Header("モデルの各パラメータ")]
    [SerializeField]
    private float m_f = 0.04f;

    [SerializeField]
    private float m_k = 0.06f;

    [SerializeField]
    private float m_simulateSpeed = 1f;

    [SerializeField]
    private int m_step = 8;

    [SerializeField]
    private float m_gridSize = 0.01f;

    [SerializeField]
    private Texture2D m_noiseTexture;

    [Space(10)]
    [SerializeField]
    private float m_Du = 2e-5f;

    [SerializeField]
    private float m_Dv = 1e-5f;
    

    [SerializeField]
    [Header("初期化項目")] 
    private int m_initialQuadSize;
    

    [SerializeField] 
    private int m_textureSize;

    public int TextureSize { get { return m_textureSize; } }

    private int m_updateGrayScottKernel;
    private Vector3Int m_groupSize;


    // Start is called before the first frame update
    void Awake()
    {
        CreateTexture();

        InitializeTexture();

        InitializeUpdateGrayScott();

        // 生成したRenderTextureをマテリアルに設定
        GetComponent<Renderer>().sharedMaterial.SetTexture("_MainTex", m_uvTexture);
    }

    void Update() {

        for(int i = 0; i < m_step; ++i)
            UpdateGrayScott();

    }

    private void CreateTexture() {

        m_uvTexture = new RenderTexture(m_textureSize, m_textureSize, 0, RenderTextureFormat.RG32);
        m_uvTexture.wrapMode = TextureWrapMode.Clamp;
        m_uvTexture.enableRandomWrite = true;
        m_uvTexture.Create();

    }

    private void InitializeTexture() {
        
        int initializeTextureKernel = m_grayScottCalculator.FindKernel("InitializeTexture");

        // テクスチャの初期状態を設定
        m_grayScottCalculator.SetInt("InitialQuadSize", m_initialQuadSize);
        m_grayScottCalculator.SetInt("minThreadhold", m_textureSize / 2 - m_initialQuadSize / 2);
        m_grayScottCalculator.SetInt("maxThreadhold", m_textureSize / 2 + m_initialQuadSize / 2);
        m_grayScottCalculator.SetTexture(initializeTextureKernel, "Texture", m_uvTexture);
        m_grayScottCalculator.SetTexture(initializeTextureKernel, "NoiseTexture", m_noiseTexture);

        m_grayScottCalculator.GetKernelThreadGroupSizes(initializeTextureKernel, out uint x, out uint y, out uint z);

        m_grayScottCalculator.Dispatch(initializeTextureKernel,
            m_textureSize / (int)x,
            m_textureSize / (int)y,
            (int) z);

    }

    private void InitializeUpdateGrayScott() {

        m_updateGrayScottKernel = m_grayScottCalculator.FindKernel("UpdateGrayScotte");
        
        // 実行する際のグループサイズを求める
        m_grayScottCalculator.GetKernelThreadGroupSizes(m_updateGrayScottKernel, out uint x, out uint y, out uint z);
        m_groupSize = new Vector3Int(m_textureSize / (int)x, m_textureSize / (int)y, (int)z);

        // GrayScottのパラメータを設定
        m_grayScottCalculator.SetFloat("f", m_f);
        m_grayScottCalculator.SetFloat("k", m_k);
        m_grayScottCalculator.SetFloat("dt", m_simulateSpeed);
        m_grayScottCalculator.SetFloat("dx", m_gridSize);
        m_grayScottCalculator.SetFloat("Du", m_Du);
        m_grayScottCalculator.SetFloat("Dv", m_Dv);
        m_grayScottCalculator.SetInt("size", m_textureSize);

        m_grayScottCalculator.SetTexture(m_updateGrayScottKernel, "Texture", m_uvTexture);
        
    }

    private void UpdateGrayScott() {

        m_grayScottCalculator.Dispatch(m_updateGrayScottKernel, m_groupSize.x, m_groupSize.y, m_groupSize.z);

    }

}

ComputeShader

#pragma kernel InitializeTexture
#pragma kernel UpdateGrayScotte


RWTexture2D<float2> Texture;
Texture2D<float4> NoiseTexture;

int InitialQuadSize;
int minThreadhold;
int maxThreadhold;

[numthreads(8, 8, 1)]
void InitializeTexture(uint3 id : SV_DispatchThreadID) {
    
    float2 initValue;

    if(minThreadhold <= id.x && id.x <= maxThreadhold &&
        minThreadhold <= id.y && id.y <= maxThreadhold)
        initValue = float2(0.5, 0.25);
    else
        initValue = float2(1, 0);

    initValue.x += NoiseTexture[id.xy].x * 0.1;
    initValue.y += NoiseTexture[id.xy].y * 0.1;

    Texture[id.xy] = initValue;

}

float dt;
float dx;
float Du;
float Dv;
float f;
float k;
int size;

[numthreads(8, 8, 1)]
void UpdateGrayScotte(uint3 id : SV_DispatchThreadID) {

    float u = Texture[id.xy].x;
    float v = Texture[id.xy].y;

    float2 laplacian =
        Texture[id.xy + uint2(-1, 0)] +
        Texture[id.xy + uint2(1, 0)] +
        Texture[id.xy + uint2(0, -1)] +
        Texture[id.xy + uint2(0, 1)] - 4 * Texture[id.xy];

    laplacian /= (dx*dx);
    
    float dudt = Du * laplacian.x - u * v * v + f * (1.0 - u);
    float dvdt = Dv * laplacian.y + u * v * v - (f + k) * v;
    
    Texture[id.xy] = float2(u + dt * dudt, v +dt * dvdt);
}

解説

初期化

平面上に2つの物質u、vが存在します。それらを管理するために、RenderTextureを作成します。フォーマットがRenderTextureFormat.RG32となっている点に注意してください。

また、m_uvTextureのR成分を物質U、G成分を物質Vとします。

        // Simulator.csの74行目から
    private void CreateTexture() {

        m_uvTexture = new RenderTexture(m_textureSize, m_textureSize, 0, RenderTextureFormat.RG32);
        m_uvTexture.wrapMode = TextureWrapMode.Clamp;
        m_uvTexture.enableRandomWrite = true;
        m_uvTexture.Create();

    }

次に、RenderTextureに初期状態を書き込みます。初期状態はテクスチャの中央に指定された幅で(u, v) = (0.5, 0.25)、それ以外は(u, v) = (1, 0)にします。下の画像が実行した結果となります。

// GrayScottCalclator.computeの13行目から
void InitializeTexture(uint3 id : SV_DispatchThreadID) {
    
    float2 initValue;

    if(minThreadhold <= id.x && id.x <= maxThreadhold &&
        minThreadhold <= id.y && id.y <= maxThreadhold)
        initValue = float2(0.5, 0.25);
    else
        initValue = float2(1, 0);

f:id:aizu-vr:20200810174051j:plain

最後にテクスチャにランダム性を追加します。

今回はスクリプトからではなく、事前にランダムな値が書かれたノイズテクスチャを作成します。ノイズテクスチャの作成は以下のサイトのプログラムをお借りして作成しました。

【Unity】ノイズと創造をつなぐ - Qiita

f:id:aizu-vr:20200810174537p:plain

あとはこのノイズテクスチャの値を先程の初期状態に足し合わせます。結果の画像はあまり違いが分かりにくいのでここでは載せません。

        // GrayScottCalculator.computeの23行目から
    initValue.x += NoiseTexture[id.xy].x * 0.1;
    initValue.y += NoiseTexture[id.xy].y * 0.1;

    Texture[id.xy] = initValue;

また、状態の更新を行う際に必要となる各パラメータは実行前にComputeShaderに設定します。

        // Simulator.csの103行目から
    private void InitializeUpdateGrayScott() {

        m_updateGrayScottKernel = m_grayScottCalculator.FindKernel("UpdateGrayScotte");
        
        // 実行する際のグループサイズを求める
        m_grayScottCalculator.GetKernelThreadGroupSizes(m_updateGrayScottKernel, out uint x, out uint y, out uint z);
        m_groupSize = new Vector3Int(m_textureSize / (int)x, m_textureSize / (int)y, (int)z);

        // GrayScottのパラメータを設定
        m_grayScottCalculator.SetFloat("f", m_f);
        m_grayScottCalculator.SetFloat("k", m_k);
        m_grayScottCalculator.SetFloat("dt", m_simulateSpeed);
        m_grayScottCalculator.SetFloat("dx", m_gridSize);
        m_grayScottCalculator.SetFloat("Du", m_Du);
        m_grayScottCalculator.SetFloat("Dv", m_Dv);
        m_grayScottCalculator.SetInt("size", m_textureSize);

        m_grayScottCalculator.SetTexture(m_updateGrayScottKernel, "Texture", m_uvTexture);
        
    }

状態の更新

C#スクリプトの方ではとても単純でstep回計算を行います。

        // Simulator.csの67行目から
    void Update() {

        for(int i = 0; i < m_step; ++i)
            UpdateGrayScott();

    }

次にComputeShaderの解説します。

// GrayScottCalculator.computeの38行目から
void UpdateGrayScotte(uint3 id : SV_DispatchThreadID) {

    float u = Texture[id.xy].x;
    float v = Texture[id.xy].y;

    float2 laplacian =
        Texture[id.xy + uint2(-1, 0)] +
        Texture[id.xy + uint2(1, 0)] +
        Texture[id.xy + uint2(0, -1)] +
        Texture[id.xy + uint2(0, 1)] - 4 * Texture[id.xy];

    laplacian /= (dx*dx);
    
    float dudt = Du * laplacian.x - u * v * v + f * (1.0 - u);
    float dvdt = Dv * laplacian.y + u * v * v - (f + k) * v;
    
    Texture[id.xy] = float2(u + dt * dudt, v + dt * dvdt);
}

始めに座標 (x, y) に対するUとVの値を取り出します。

次に $\Delta u$ と $\Delta v$ を求めて変数laplacianに保存しています。また、計算する際にUとVを一度に計算しています。

求められたlaplacianを隣の成分との距離 $dx$ の2乗で割ります。

後は、前述した計算式通りに進めて最後に1回のシミュレーションの時間 $dt$ を掛けてテクスチャを更新します。

結果

物質Uが赤色、物質Vが緑色です。また、初期化項目のInital Quad Sizeは50、Texture Sizeは256で実行しています。

f:id:aizu-vr:20200811235922g:plain

このままだとおどろおどろしいので、シェーダーで色を変更します。

シェーダープログラム

Shader "Unlit/ColorShader"
{
    Properties
    {
        _MainTex ("Texture", 2D) = "white" {}
    }
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        LOD 100

        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            // make fog work
            #pragma multi_compile_fog

            #include "UnityCG.cginc"

            struct appdata
            {
                float4 vertex : POSITION;
                float2 uv : TEXCOORD0;
            };

            struct v2f
            {
                float2 uv : TEXCOORD0;
                UNITY_FOG_COORDS(1)
                float4 vertex : SV_POSITION;
            };

            sampler2D _MainTex;
            float4 _MainTex_ST;

            v2f vert (appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.uv = TRANSFORM_TEX(v.uv, _MainTex);
                UNITY_TRANSFER_FOG(o,o.vertex);
                return o;
            }

            fixed4 frag (v2f i) : SV_Target
            {
                // sample the texture
                fixed4 col = tex2D(_MainTex, i.uv);
                // apply fog
                UNITY_APPLY_FOG(i.fogCoord, col);
                return fixed4(col.r, col.r, col.r, 1);
            }
            ENDCG
        }
    }
}

やっていることは単純で、赤色を白色に変換しているだけです。こうすることで次のように変わります。

f:id:aizu-vr:20200812001455g:plain