概説

Unityでオブジェクトファイルのメッシュ形状に応じた反応拡散シミュレーションを行いました。

全体の流れ

1.オブジェクトファイルのインポート
2.オブジェクトの頂点インデックスから隣接リストを生成する
3.各頂点における反応拡散を計算する
4.計算結果を頂点カラーに反映させる

I.オブジェクトファイルのインポート

最初に、適当なUnityプロジェクト(3D)を作成します。次に、オブジェクトファイル(.obj)をプロジェクトビューのAssetにいれます。今回はUtah teapotとStanford bunnyをobjファイルに変換してからインポートしています。
インポートしたオブジェクトファイルのセッティングは変更する必要があります。インスペクタビューの"Model"=>"Meshes"の"Read/Write Enabled"にチェックを入れておきスクリプトからメッシュが読み込めるようにしておきます。

また"Materials"の"Import Materials"のチェックを外しておきます。
設定を変更したオブジェクトをヒエラルキービューにドラッグ＆ドロップします。
新しいマテリアルを作成し、アタッチします。
Standard Surface shaderを新しく作成し、中身を以下のように書き換えマテリアルにひもづけます（マテリアル標準のシェーダーだと頂点カラーを反映してくれないらしいので）。

Shader "Custom/NewSurfaceShader"
{
    
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        LOD 200

        CGPROGRAM
        // Physically based Standard lighting model, and enable shadows on all light types
        #pragma surface surf Lambert vertex:vert

        // Use shader model 3.0 target, to get nicer looking lighting
        #pragma target 3.0

        struct Input
        {
            float4 vertColor;
        };

        void vert(inout appdata_full v, out Input o){
            UNITY_INITIALIZE_OUTPUT(Input,o);
            o.vertColor=v.color;
        }

        void surf (Input IN, inout SurfaceOutput o)
        {
            o.Albedo = IN.vertColor.rgb;
        }
        ENDCG
    }
    FallBack "Diffuse"
}

II.オブジェクトの頂点インデックスから隣接リストを生成する

私たちがこれからやりたいことは、各頂点における活性因子と抑制因子の濃度を計算し、状態を更新することです。濃度計算の中には拡散項があり、2頂点間の濃度勾配によって隣接する頂点に活性因子と抑制因子が拡散されます。

この拡散計算を行うためには「ある1つの頂点がどの頂点と隣接しているか」という情報が必要になってきます。
ここで、UnityのMeshが持つ情報を見てみましょう。

vertices（頂点の位置）
triangles（メッシュ内の全ての三角形を含む配列）
normals（メッシュの法線）
tangents（メッシュの接線）
uv（テクスチャ座標）
colors（メッシュの各頂点の色）
他にもいろいろ。

trianglesはint型の一次元配列で「どの頂点を結んで三角形ポリゴンをつくるか」という頂点インデックス情報を持っています。1つの三角形ポリゴンごとに3つの頂点が指定されます。時計回りが表です。

この情報を使えば、ある頂点が他のどの頂点につながっているのかが分かりそうです。ただし、この形式のままではなく隣接リスト形式を変える必要があります。隣接リストはある頂点をキーとしてその頂点が隣接している頂点のみをリストにしたものです。アルゴリズムとしては3の倍数ごとに値を読み、隣接する頂点をそれぞれのリストに加えます。
triangles=[0,1,3]ならば、

List[0]に1を加える＆List[1]に0を加える
List[1]に3を加える＆List[3]に1を加える
List[3]に0を加える＆List[0]に3を加える
という操作を行い隣接リストを作ります。

上記図のように頂点インデックスが

vertices=[0,1,3,
          0,3,2]

ならば、隣接リスト（リスト型の配列）は

0:[1,2,3]
1:[0,3]
2:[0,3]
3:[0,1,2]

です。
メッシュを読み込んでから隣接リストを作るまでのコードを以下に示します。（上に書いたやり方では2つの三角形ポリゴンが1つの辺を共有していたとき隣接リストの頂点が2回登録されて重複します。そのため実装では重複を避けるため「隣接リスト」という名前ですがListではなくHashedSetを使っています。）

//メッシュの取得
mesh=GetComponent<MeshFilter>().mesh;
//メッシュの頂点インデックス
int[] tris=mesh.triangles;
Color[] colors=new Color[mesh.vertexCount];
        
//隣接リストの宣言
adjacencyList=new HashSet<int>[mesh.vertexCount];
//隣接リストのインスタンス作成
for(int i=0;i<adjacencyList.Length;i++){
      adjacencyList[i]=new HashSet<int>();
}

//trianglesから隣接リストの作成
for(int i=0;i<tris.Length/3;i++){
    for(int j=0;j<3;j++){
        int a=i*3+j;
        int b=i*3+(j+1)%3;

        adjacencyList[tris[a]].Add(tris[b]);
        adjacencyList[tris[b]].Add(tris[a]);
    }
}

III.各頂点における反応拡散を計算する

反応拡散の例としてチューリングパターンを描いてみます。
活性因子と抑制因子の濃度をそれぞれ$p$,$q$と表します。$t$秒後の$i$番目の頂点における活性因子と抑制因子の濃度はそれぞれ$p(i,t)$,$q(i,t)$とします。

反応項

活性因子は、因子の産生を促進します。抑制因子は、因子の産生を抑制します。詳しくは以下の表のように定義されます。

	活性因子によって	抑制因子によって
活性因子の産生は	促進される(+0.6)	抑制される(-1.0)
抑制因子の産生は	促進される(+1.5)	抑制される(-2.0)

また、活性因子はある程度の濃度以上になると自身の濃度の3乗で産生が抑制されます。

これら産生量の変化は以下の関数で示されます。これを反応項と呼びます。

\begin{align}
f(p,q)&=+0.6p-1.0q-p^3\\
g(p,q)&=+1.5p-2.0q
\end{align}

拡散項

ある頂点にある活性/抑制因子は周辺の隣接した頂点に拡散します。どれくらい移動するかはそれぞれの頂点間の濃度勾配と拡散係数$d_p$,$d_q$に依存します。頂点間の距離はすべて一定の長さ$dx$とします（ほんとうは2頂点間の長さはばらばらです）。

計算式

$t+dt$秒の濃度は$t$秒における濃度に反応項と拡散項の結果を足したものとなります。拡散項では隣接リストを使って隣接する頂点との拡散を計算しています。

p(i,t+dt)=p(i,t)+\Biggl(f(p(i,t),q(i,t))+\frac{d_p}{dx^2}\sum_{j=neighbor}\Bigl(p(j,t)-p(i,t)\Bigl)dt\Biggr)\\
q(i,t+dt)=q(i,t)+\Biggl(g(p(i,t),q(i,t))+\frac{d_q}{dx^2}\sum_{j=neighbor}\Bigl(q(j,t)-q(i,t)\Bigl)dt\Biggr)\\

チューリングパターンの詳細については以前の記事：2次元Turingパターンの生成も参照してください。

IV.計算結果を頂点カラーに反映させる

計算した各頂点の活性因子濃度を色情報としてMesh.Colorに格納します。活性因子と抑制因子の濃度は-1.0~+1.0の値をとりうるので2で割って0.5を足すことで0~1に正規化します。

最終的なコード

全体をまとめたC#スクリプトを目的のオブジェクトにアタッチすることで冒頭のgif動画に描画されます。

using System.Linq;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class CalculateTuring : MonoBehaviour
{
    //宣言
    Mesh mesh;
    HashSet<int>[] adjacencyList;//隣接リスト

    float[] pconc,qconc;//活性因子の濃度、抑制因子の濃度
    float[] temppconc,tempqconc;//前世代の濃度を保存する
    float dp=0.0002f;
    float dq=0.01f;
    float dx=0.2f;
    float dt=0.05f;

  
    void Start()
    {
        //メッシュの取得
        mesh=GetComponent<MeshFilter>().mesh;
        //メッシュの頂点インデックス
        int[] tris=mesh.triangles;
        //頂点の色情報
        Color[] colors=new Color[mesh.vertexCount];

        adjacencyList=new HashSet<int>[mesh.vertexCount];
        //隣接リストのインスタンス作成
        for(int i=0;i<adjacencyList.Length;i++){
            adjacencyList[i]=new HashSet<int>();
        }

        //trianglesから隣接リストの作成
        for(int i=0;i<tris.Length/3;i++){
           for(int j=0;j<3;j++){
              int a=i*3+j;
              int b=i*3+(j+1)%3;

              adjacencyList[tris[a]].Add(tris[b]);
              adjacencyList[tris[b]].Add(tris[a]);
           }
        }
        //濃度の初期化
        pconc=new float[mesh.vertexCount];
        qconc=new float[mesh.vertexCount];

        temppconc=new float[mesh.vertexCount];
        tempqconc=new float[mesh.vertexCount];

        Random.InitState(2019);//乱数のシード値固定

        //初期濃度には小さなノイズを与えておく
        for(int i=0;i<mesh.vertexCount;i++){
            pconc[i]=Random.value*0.01f;
            qconc[i]=Random.value*0.01f;
            temppconc[i]=0.0f;
            tempqconc[i]=0.0f;
        }
        //Debug.Log(pconc[8]);
    }


    void Update()
    {
        pconc.CopyTo(temppconc,0);
        qconc.CopyTo(tempqconc,0);

        for(int i=0;i<pconc.Length;i++){
            pconc[i]=temppconc[i]+dt*(f(temppconc[i],tempqconc[i])+dp/dx/dx*DiffisionPQ(temppconc,i));
            qconc[i]=tempqconc[i]+dt*(g(temppconc[i],tempqconc[i])+dq/dx/dx*DiffisionPQ(tempqconc,i));
        }
        ChangeMeshColor(pconc,qconc);
    }

   
    //拡散項：1つの頂点について近隣頂点との濃度差からdt秒後の増減分を計算する
    float DiffisionPQ(float[] conc,int index){
        float ans=0.0f;
        foreach (int neighbor in adjacencyList[index])
        {
            ans+=conc[neighbor]-conc[index];
        }
        return ans;
    }
    //反応項
    float f(float p,float q){
        return 0.6f*p-1.0f*q-p*p*p;
    }
    float g(float p,float q){
        return 1.5f*p-2.0f*q;
    }

    //頂点の色情報を変える
    void ChangeMeshColor(float[] p,float[] q){
        Color[] colors=new Color[mesh.vertexCount];
        for(int i=0;i<colors.Length;i++){
            colors[i]=new Color(p[i]*10.0f/2.0f+0.5f,p[i]*10.0f/2.0f+0.5f,p[i]*10.0f/2.0f+0.5f);
        }
        mesh.colors=colors;
       
        mesh.RecalculateNormals();
        mesh.RecalculateBounds();
    }
}

参照

おもちゃラボ、【Unityシェーダ入門】頂点カラーを表示するシェーダを作る(http://nn-hokuson.hatenablog.com/entry/2017/03/29/195041)
三浦岳、発生の数理、京都大学学術出版会、2015、6.6
jonki, 【Unity】オリジナルMeshの作成とVertex ColoredなShader(https://www.jonki.net/entry/20140305/1394024333)