はじめに

本記事ではUnityを用いて以下のシミュレーション映像を作るためのプロセスをメモついでに書いていきます．

https://yoshidev.tumblr.com/post/158231100497/oscillation-sphere-simulating-kuramoto-model

なお，記述に関して，不正確な部分が生じると思いますので指摘頂ければ幸いです．
また今回の実装にあたって，Keijiro TakahashiさんのOSSを参考している部分があります．

同期現象と蔵本モデルについて

まず，同期現象とは，例えば蛍の集団の点滅周期が森の中の蛍の集団で同期が見られたり，
人間の脳内の領野間でも脳活動の同期によって脳の可塑性を実現している説等，自然界に様々に見られる異なる振動子のリズムが揃っていく現象のことです．

そこで蔵本モデルとは

蔵本由紀によって提案された同期現象を記述する数学モデルである。特に、相互作用のある非線形振動子集団の振る舞いを記述するモデルである。
(https://ja.wikipedia.org/wiki/蔵本モデル)

本質的なことは蔵本先生の講義ノートとか著書を是非読んでみてください．
また同期現象に関してもSYNCという書籍が読み物としても面白いのでオススメです．

今回はこの蔵本モデルをUnityを用いて簡単なシミュレーションを行います．

実装

環境

・Unity 5.6.0b3
・MacBook Pro (15-inch, Late 2016)
・CPU：2.9 GHz Intel Core i7/　メモリ：16 GB
・Radeon Pro 460 4096 MB
・macOS Sierra 10.12.3

使用するスクリプト

・ModelSimulation.cs
・Floor.cs
・CameraController.cs
・RandomBoxMuller.cs (参考)
・Kernel.shader
・Surface.shader

モデルの更新処理の実装

まず最初に以下の数式で表されるモデルの更新処理の実装を行います．

文字の意味は
・文字の上のドットは時間微分
・R，Θ:集団の同期度を表す複素パラメタ
・K:それぞれの振動子の相互作用の結合係数
・N:振動子数
・i:振動子インデックス
・ω:振動子の自然振動数
・φ:振動子の位相

この式に従って実装していきます．

ModelSimulation.csとKernel.shaderの役割

今回の実装では多くの振動子集団に対して上記の式に従い更新処理を行うため計算時間の肥大化が考えられます．
そのため，上記の速度の式の変形で振動子間の相互作用は陽には現れない，
を用いることで個々の振動子の位相およびその速度の更新についてはShaderによって並列計算させることを考えます．
同期度を表す複素パラメタの算出については振動子集合の角位相が影響するため ModelSimulation.csにおいて計算を行います．

ModelSimulation.cs

Shaderでの処理のために，RenderTextureを位相，速度のバッファに，Texture2Dを自然振動数のバッファとして利用します．
振動子の数はRenderTextureのサイズ制限により4096までとしています．

フィールドおよびプロパティ

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEditor;

public class ModelSimulation : MonoBehaviour {

    [SerializeField, Range(1, 4096)]
    private int _pointNum; //振動子の数
    public int pointNum {
        get { return _pointNum;}
        private set { _pointNum = value;}
    }

    [SerializeField] 
    private float _baseFreq; //振動数分布の中央値
    public float baseFreq {
        get { return _baseFreq;}
        private set { _baseFreq = value;}
    }

    [SerializeField]
    private float _connectionCoefficient; //結合係数
    public float connectionCoefficient {
        get { return _connectionCoefficient;}
        private set { _connectionCoefficient = value;}
    }

    [SerializeField] Shader kernelShader; //位相，速度更新用のShader
    private Material kernelMat; 

    private Texture2D naturalFreqBuffer; //自然振動数のバッファ
    private RenderTexture phaseBuffer; //位相のバッファ
    private RenderTexture velocityBuffer; //位相の速度バッファ

}

マテリアル，バッファの初期化のための関数

ここは素直に．
フォーマットを浮動小数点にしておきます．

public class ModelSimulation : MonoBehaviour {
            ~~~~~略~~~~~
    Material CreateMaterial(Shader shader)
    {
        var material = new Material(shader);
        material.hideFlags = HideFlags.DontSave;
        return material;
    }

    RenderTexture createRTBuffer(int size){
        var format = RenderTextureFormat.ARGBFloat;
        var buffer = new RenderTexture(size, 1, 0, format);
        buffer.hideFlags = HideFlags.DontSave;
        buffer.filterMode = FilterMode.Point;
        buffer.wrapMode = TextureWrapMode.Clamp;
        return buffer;
    }

    Texture2D createT2Buffer(int size){

        var texture = new Texture2D(size, 1, TextureFormat.RGBAFloat, false);

        for(int i = 0; i < size; i++){
            texture.SetPixel (i, 1, new Color(0, 0, 0, 0));
        }
        texture.Apply ();
        return texture;
    }   
}

初期化

各バッファ，マテリアルの初期化を行います．
マテリアルに流し込むデータに関してはKernel.shaderの部分で記述します．

Graphics.Blit(Texture source, RenderTexture dest, Material mat, int pass = -1);

を用いることで，Shaderのパスを指定して演算結果をdestに格納します．
詳しくはhttps://docs.unity3d.com/jp/540/ScriptReference/Graphics.Blit.htmlを参考にしてください．

public class ModelSimulation : MonoBehaviour {
            ~~~~~略~~~~~
            
    void Start () {
        Initialize ();
    }   
        
    void Initialize(){

        naturalFreqBuffer = createT2Buffer (pointNum);
        phaseBuffer = createRTBuffer (pointNum);
        velocityBuffer = createRTBuffer (pointNum);

        InitializeMat ();　//マテリアル生成とパラメタ初期化
        InitializePosition ();  //Kernel.shaderによる位相の初期化
        InitializeNaturalFreqs (); //μ=0のがウス分布から自然周波数の初期化
    }

    void InitializeMat(){

        kernelMat = CreateMaterial (kernelShader);
        kernelMat.SetInt ("_PointNum", pointNum);
        kernelMat.SetFloat ("_K", connectionCoefficient);
        kernelMat.SetFloat ("_BaseFreq", baseFreq);
        kernelMat.SetFloat ("_DeltaTime", 0);
    }

    void InitializePosition(){
        
        Graphics.Blit (null, phaseBuffer, kernelMat, 0);
        kernelMat.SetTexture ("_PhaseTex", phaseBuffer);
    }

    void InitializeNaturalFreqs(){

        RandomBoxMuller random = new RandomBoxMuller();

        for(int i = 0; i < pointNum; i++){
            naturalFreqBuffer.SetPixel (1, i, new Color(((float)random.next(0, 2.0, true) - 0.5f) * baseFreq , 0, 0, 0)); //ガウス分布に従う自然周波数
        }
        naturalFreqBuffer.Apply ();

        kernelMat.SetTexture ("_NaturalFreqTex", naturalFreqBuffer);
        surfaceMat.SetTexture ("_NaturalFreqTex", naturalFreqBuffer);
    }
}

更新処理

初期化と同じようにGraphics.Blitを用いてKernel.shaderによる更新処理を行います．

public class ModelSimulation : MonoBehaviour {
            ~~~~~略~~~~~
    
    void Update () {
        
        kernelMat.SetFloat ("_DeltaTime", Time.deltaTime);
        
        UpdateVelocity (); //速度バッファの更新
        UpdatePhase (); //位相バッファの更新
    }
                    
    void UpdateVelocity(){
    
        Graphics.Blit (null, velocityBuffer, kernelMat, 1);

        kernelMat.SetTexture ("_VelocityTex", velocityBuffer);
    }

    void UpdatePhase(){
        
        Graphics.Blit (null, phaseBuffer, kernelMat, 2);

        float[] param = calcParams ();

        kernelMat.SetTexture ("_PhaseTex", phaseBuffer);
        kernelMat.SetFloat ("_ParamTheta", param[0]);
        kernelMat.SetFloat ("_ParamR", param[1]);
    }

    //同期度を表す複素パラメタの算出
    float[] calcParams(){

        Texture2D tex = new Texture2D(phaseBuffer.width, phaseBuffer.height, TextureFormat.RGBAFloat, false);

        RenderTexture.active = phaseBuffer;
        tex.ReadPixels(new Rect(0, 0, phaseBuffer.width, phaseBuffer.height), 0, 0);
        tex.Apply();

        float real = 0;
        float imag = 0;

        for(int i = 0; i < pointNum; i++){
            float phai = tex.GetPixel (i, 0).r;
            real += Mathf.Cos (phai);
            imag += Mathf.Sin (phai);
        }

        real /= (float)pointNum;
        imag /= (float)pointNum;

        float paramTheta = Mathf.Atan (imag/real);
        float paramR = Mathf.Sqrt (real * real + imag * imag);

        radius = paramR * 0.7f + 0.5f;


        return new float[2]{paramTheta, paramR};
    }
}

特に難しいことはしていないですがこれでModelSimulation.csnのモデルの更新処理の実装は終わりです．
次は実際に位相の初期化，位相，速度の更新を行うKernel.shaderの実装をします．

Kernel.shader

プロパティ

Shaderに流し込むデータをここで定義します．

Shader "Hidden/Kernel"
{
    Properties
    {
       _PhaseTex ("-", 2D) = ""{} //位相バッファ
       _VelocityTex ("-", 2D) = ""{} //速度バッファ
       _NaturalFreqTex("-", 2D) = ""{} //自然振動数バッファ
       _K           ("-", Float) = 0 //相互作用の結合係数
       _ParamTheta  ("-", Float) = 0 //同期度を表す複素パラメタのΘ
       _ParamR      ("-", Float) = 0 //同期度を表す複素パラメタのR
       _BaseFreq    ("-", Float) = 0 //自然振動数の分布の中央値
       _DeltaTime   ("-", Float) = 0 //dt
       _PointNum    ("-", int) = 0 //振動子の数
    }

    CGINCLUDE

    #include "UnityCG.cginc"
   
    sampler2D _PhaseTex;
    sampler2D _VelocityTex;
    sampler2D _NaturalFreqTex;

    float _BaseFreq;
    float _ParamTheta;
    float _ParamR;
    float _K;
    float _DeltaTime;

    int _PointNum;
    
    ENDCG
}

パスの定義

SubShader{}

の中でパスの定義を行います．
上から順にGraphics.Blit()で指定するpassになります．
VertexShaderはUnityCG.cgingで定義されているvert_imgを用いて，FragmentShaderのみ自前で実装します．

Shader "Hidden/Kernel"
{
            ~~~~~略~~~~~
            
    CGINCLUDE
            ~~~~~略~~~~~
    ENDCG
                
    SubShader
    {
        Pass //位相の初期化 pass=0
        {
            CGPROGRAM
            #pragma target 3.0
            #pragma vertex vert_img
            #pragma fragment frag_init_phase
            ENDCG
        }


        Pass //速度の更新 pass=1
        {
            CGPROGRAM
           #pragma target 3.0
           #pragma vertex vert_img
           #pragma fragment frag_update_velocity
            ENDCG
  
        }

        Pass //位相の更新 pass=2
        {
            CGPROGRAM
            #pragma target 3.0
            #pragma vertex vert_img
            #pragma fragment frag_update_phase
            ENDCG
        }
    }
    
}

FragmentShaderの実装

振動子の状態の初期化，更新を上記の式の通りに実装するだけです．

Shader "Hidden/Kernel"
{
            ~~~~~略~~~~~
    CGINCLUDE
            ~~~~~略~~~~~

    float nrand(float2 uv, float salt)
    {
        uv += float2(salt, 0);
        return frac(sin(dot(uv, float2(12.9898, 78.233))) * 43758.5453);
    }

    float4 frag_init_phase(v2f_img i) : SV_Target {     

        float v = nrand(i.uv, 10) * _BaseFreq * 2.0;

        return float4(v,0,0,0);
    }


    float4 frag_update_velocity(v2f_img i) : SV_Target {

        //バッファのColor.rにが格納されている位相，速度を取り出して計算する．
        float p = tex2D(_PhaseTex, i.uv).x; 
        float nf = tex2D(_NaturalFreqTex, i.uv).x;
        float v = nf + _K * _ParamR * sin(_ParamTheta - p);

        return float4(v, 0, 0, 0);
    }

    float4 frag_update_phase(v2f_img i) : SV_Target {

        float p = tex2D(_PhaseTex, i.uv).x;
        float v = tex2D(_VelocityTex, i.uv).x;

        v = p + v * _DeltaTime;

        return float4(v, 0, 0, 0);
    }       
            
    ENDCG
                        
            ~~~~~略~~~~~           
}

モデルの更新処理は以上になります．
この状態でGameObjectに貼り付けて実行してみるとRが大きくなっていくのがわかると思います．
以降では可視化のための実装を行います．

可視化の実装

可視化にはModelSimulation.cs，Surface.shaderに加え，Blenderで簡単に頂点数と三角形をいい感じにした球体を用いました．
この球体は別に他のものでもいいのでどこかにあげることはしない予定です．

ModelSimulation.cs

可視化とシミュレーションを同じスクリプトに書くのはどうかと思ったのですが，
ファイル数を減らしたかったのでそのままModelSimulation.csに可視化の実装も書いちゃいました．

フィールドおよびプロパティ

Shaderに流し込む経過時間，色とかの設定のためのフィールドとプロパティをつらつらと．

public class ModelSimulation : MonoBehaviour {

            ~~~~~略~~~~~
            
    [SerializeField] Shader surfaceShader;
    [SerializeField] private Mesh mesh;　//入力メッシュ
    private Material surfaceMat;
            
    [SerializeField] private Color _color1;
    public Color color1 {
        get{ return _color1;}
        set{ _color1 = value;}
    }

    [SerializeField] private Color _color2;
    public Color color2 {
        get{ return _color2;}
        set{ _color2 = value;}
    }

    [SerializeField, Range(0f, 1.0f)] private float _metallic = 0.5f;
    public float metallic {
        get{ return _metallic;}
        set{ _metallic = value;}
    }

    [SerializeField, Range(0f, 1.0f)] private float _smoothness = 0.5f;
    public float smoothness {
        get{ return _smoothness;}
        set{ _smoothness = value;}
    }

    [SerializeField, Range(0.5f, 2.0f)] private float _radius = 1.0f;
    public float radius {
        get{ return _radius;}
        set{ _radius = value;}
    }
        
    private float elapsedTime;
    
            ~~~~~略~~~~~
}

初期化

例によってここでマテリアルの初期化を行います．
また，それぞれの頂点に対して振動子を割り当てるためにUVの初期化を行います．
頂点数と振動子の数を合わせる必要はありませんがなるべく近い方がいいかもしれない．

public class ModelSimulation : MonoBehaviour {

            ~~~~~略~~~~~
            
    void Initialize(){

            ~~~~~略~~~~~
            
        elapsedTime = 0f;
        int vNum = mesh.vertexCount;

        Vector2[] uv = new Vector2[vNum];

        for(int i = 0; i < vNum; i++){
            uv [i] = new Vector2 (((float)i + 0.5f)/ (float)vNum, 0.5f);
        }

        mesh.uv = uv;
    }
    
    void InitializeMat(){

        surfaceMat = CreateMaterial (surfaceShader);
            ~~~~~略~~~~~   
    }
    
    
            ~~~~~略~~~~~
}

更新処理

そのままですね，毎フレームマテリアルにバッファ等を流し込んで

Graphics.DrawMesh(Mesh mesh, Matrix4x4 matrix, Material material, int layer)

で描画を行うだけです．

public class ModelSimulation : MonoBehaviour {

            ~~~~~略~~~~~
            
    void Update () {
        

            ~~~~~略~~~~~
        elapsedTime += Time.deltaTime;
        DrawMesh ();
    }   
    
            ~~~~~略~~~~~
            
    void DrawMesh(){
        
        surfaceMat.SetColor ("_Color1", color1);
        surfaceMat.SetColor ("_Color2", color2);
        surfaceMat.SetTexture ("_PhaseTex", phaseBuffer);
        surfaceMat.SetFloat ("_Metallic", metallic);
        surfaceMat.SetFloat ("_Smoothness", smoothness);
        surfaceMat.SetFloat ("_Radius", radius);
        surfaceMat.SetFloat ("_ElapsedTime", elapsedTime);
        surfaceMat.SetFloat ("_BaseFreq", baseFreq);

        Graphics.DrawMesh (mesh, transform.localToWorldMatrix, surfaceMat, gameObject.layer);
    }           
}

以上でModelSimulation.csは完成です．
次に描画のためのShaderを書いていきます．

Surface.shader

プロパティ

Shader "Custom/Surface" {
    Properties
    {
        _PhaseTex ("-", 2D) = ""{}
        _NaturalFreqTex("-", 2D) = ""{}

        _Color1 ("-", Color) = (1, 1, 1, 1)
        _Color2 ("-", Color) = (1, 1, 1, 1)

        _Metallic ("-", Float) = 0
        _Smoothness ("-", Float) = 0

        _Radius ("-", Float) = 0
        _PointNum ("-", Float) = 0
    }

    CGINCLUDE

    sampler2D _PhaseTex;
    sampler2D _NaturalFreqTex;
    float4 _PositionTex_TexelSize;
    float4 _NaturalFreqTex_TexelSize;

    half4 _Color1;
    half4 _Color2;

    half _Metallic;
    half _Smoothness;

    float _Radius;
    float _PointNum;

    float _ElapsedTime;
    float _BaseFreq;
    
    ENDCG

}

Vertex, Framgment Shaderの記述

振動子の位相によって頂点の座標を変換します．
特に他に書くことはなさそうなのでこれぐらいにしておきます.

Shader "Custom/Surface" {
            ~~~~~略~~~~~
    CGINCLUDE
            ~~~~~略~~~~~
            
    struct Input {
        half color;
    };

    float nrand(float2 uv, float salt)
    {
        uv += float2(salt, 0);
        return frac(sin(dot(uv, float2(12.9898, 78.233))) * 43758.5453);
    }

    void vert(inout appdata_base v, out Input data)
    {
        UNITY_INITIALIZE_OUTPUT(Input, data);

        float2 uv = v.texcoord;
        float p = tex2Dlod(_PhaseTex, float4(uv , 0, 0)).x;
        float theta = _BaseFreq * _ElapsedTime + p;

        v.vertex.xyz = v.vertex.xyz + v.normal * (sin(theta) + 0.3) * (length(v.vertex.xyz)) * _Radius;

        float ln = (sin(theta) + 1.0) / 2.0;
        data.color = ln;
    }
    ENDCG
    
     SubShader
    {
        Tags { "RenderType"="Opaque" }

        CGPROGRAM

        #pragma surface surf Standard vertex:vert nolightmap addshadow
        #pragma target 3.0

        void surf(Input IN, inout SurfaceOutputStandard o)
        {
            o.Albedo = lerp(_Color1, _Color2, IN.color);
            o.Metallic = _Metallic;
            o.Smoothness = _Smoothness;
        }

        ENDCG
    }                   
}

Surface.shaderの実装は以上です．
現時点でこんな感じにしてやれば振動子の振る舞いが可視化されると思います．
パラメータを変えて振る舞いを鑑賞します．笑

スクリーンショット 2017-03-11 18.43.53.png

次の項目はおまけです．

見た目を整える

背景

Quadを６枚内側に向けて以下の画像のように配置．（Walls）

スクリーンショット 2017-03-11 18.53.42.png

新しくマテリアルWallMat，Wall.csを作ってそれぞれのWallに以下のように貼り付ける．

Wall.csの中身は以下．

public class Wall : MonoBehaviour {

    Texture2D emmitionTex;
    int texSize = 1024;

    void Start () {
        emmitionTex = new Texture2D (texSize, texSize);

        for(int i = 0; i < texSize; i ++){
            for(int j = 0; j < texSize; j++){
                if (i % 64 == 0 || j % 64 == 0) {
                    emmitionTex.SetPixel (i, j, new Color (1.0f, 1.0f, 1.0f)); 
                } else {
                    emmitionTex.SetPixel (i, j, new Color (32f / 255f, 32f / 255f, 32f / 255f));    
                }
            }
        }

        emmitionTex.Apply ();

        Material fMat = gameObject.GetComponent<Renderer> ().material;
        fMat.SetTexture ("_EmissionMap", emmitionTex);
        fMat.SetColor ("_EmissionColor", new Color(1.0f, 1.0f, 1.0f, 1.0f));
    }
}

Cameraとライト

Cameraの動き

テキトーにぐるぐるさせます．

public class CameraController : MonoBehaviour {

    float radius = 3.5f;
    float time = 0;

    void Update () {
        time += Time.deltaTime;
        this.transform.position = new Vector3 (radius * Mathf.Cos (time / 2.0f), radius/ 1.5f * Mathf.Sin (time / 3.0f), (radius + 1.0f) * Mathf.Sin (time / 2.0f));
        this.transform.LookAt (Vector3.zero, Vector3.up);
    }
}

ImageEffectとライト

いい感じになるようにカメラにImageEffectをつける．
いい感じになるようにPointLightを配置する．

終わりに

長くなりましたが以上になります．
この手のモデルの可視化が好きなので試してほしいモデルなどありましたら紹介ください．

また，意見や質問，歓迎しますのでよろしくお願いします．
ありがとうございました．

2017.03.11

SnobbishNoob

プログラミングとか研究とか

RenderTextureとShaderを利用した集団同期現象のシミュレーションと視覚化

はじめに

同期現象と蔵本モデルについて

実装

環境

使用するスクリプト

モデルの更新処理の実装

ModelSimulation.csとKernel.shaderの役割

ModelSimulation.cs

フィールドおよびプロパティ

マテリアル，バッファの初期化のための関数

初期化

更新処理

Kernel.shader

プロパティ

パスの定義

FragmentShaderの実装

可視化の実装

ModelSimulation.cs

フィールドおよびプロパティ

初期化

更新処理

Surface.shader

プロパティ

Vertex, Framgment Shaderの記述

見た目を整える

背景

Cameraとライト

Cameraの動き

ImageEffectとライト

終わりに