UnityでPhysarum Simulationをやってみた話

はじめに

個人的な興味と事情からUnity上でPhysarum Simulationをやってみました。
アニメーションの表現の幅が広がったらいいなというのが主な動機です。

この@entagmaのチュートリアルを元にPhysarum SimulationUnityで試作してみた#generativeart #Unity #madewithunity #creativecoding https://t.co/CuVuFmsy6m pic.twitter.com/Y2OPgJjzZy

— 北村　圭 (@KeiKitamura_JP) June 10, 2022

Physarum Simulationとは

粘菌の挙動をモデル化したアルゴリズムです。
手付けアニメーションや簡単な数式では出せないような独特な動きになるので気になっていました。
参考: https://cargocollective.com/sagejenson/physarum

アルゴリズムの概要を大まかに説明すると以下の通りです。
・シミュレーションはパーティクルと2次元データで構成
・２次元データの平面上にパーティクルを大量にばら撒く
・各パーティクルの位置・角度で２次元データをサンプリングした結果に基づきパーティクルを動かす
・パーティクルが動いた結果を２次元データに記録し、減衰させる

実装手順

まずはHoudiniのチュートリアルを見ながら作ってみて仕組みを理解するところから始めました。
理解するのに手を付けたチュートリアルは以下の２点です。
　・https://entagma.com/physarum-slime-mold/
　・https://www.youtube.com/watch?v=1JSX4VcgrIM

その後、上記のHoudiniでの実装に基づきUnityに移植を試みます。
シミュレーション部分はCompute Shaderで計算し、
２次元データはRender Textureに書き込んで結果を確認します。

実装内容

以下のURLがスクリプト・Compute Shaderの全文です。

・PhysarumSimulationHandler.cs (Compute Shaderを動かすスクリプト)
https://gist.github.com/kayk5654/6a647df906f93292b9b12af2e35f6009

・PhysarumSimulation.compute (シミュレーション本体)
https://gist.github.com/kayk5654/513dcc833d3e6e9f1f1342702d645916

シミュレーション実装（抜粋）

// パーティクルの挙動のシミュレーション
[numthreads(8,1,1)]
void Simulate(uint3 id : SV_DispatchThreadID)
{
    // インスペクタ上で設定した各種パラメータをシミュレーション用変数に代入
    // マスク用のテクスチャを使用する場合は、そのテクスチャのサンプリング結果に応じてパラメータを分ける
    float sensorAngle;
    float turnAngle;
    float sensorDistance;
    float moveDistance;
    float depositValue;
    if (_useInputTexture && _inputTexture[(uint2)_particles[id.x].position].r > 0.5)
    {
        sensorAngle = _sensorAngle1;
        turnAngle = _turnAngle1;
        sensorDistance = _sensorDistance1;
        moveDistance = _moveDistance1;
        depositValue = _depositValue1;
    }
    else 
    {
        sensorAngle = _sensorAngle0;
        turnAngle = _turnAngle0;
        sensorDistance = _sensorDistance0;
        moveDistance = _moveDistance0;
        depositValue = _depositValue0;
    }
    
    // 各パーティクルに付随するセンサー（テクスチャをサンプリングする部分）の位置を求める
    float2 sensorCPos = _particles[id.x].position + _particles[id.x].direction * sensorDistance;
    float2 sensorLPos = _particles[id.x].position + mul(_sensorLRot, _particles[id.x].direction) * sensorDistance;
    float2 sensorRPos = _particles[id.x].position + mul(_sensorRRot, _particles[id.x].direction) * sensorDistance;

    // センサーの位置がテクスチャの領域をはみ出さないように区切る
    sensorCPos = WrapPosition(sensorCPos);
    sensorLPos = WrapPosition(sensorLPos);
    sensorRPos = WrapPosition(sensorRPos);

    // テクスチャのピクセルの値を読む
    float centerValue = _resultTexture[uint2(sensorCPos)].r;
    float leftValue = _resultTexture[uint2(sensorLPos)].r;
    float rightValue = _resultTexture[uint2(sensorRPos)].r;

    // 回転行列の初期化
    _sensorRRot = GetRotMatrix(sensorAngle);
    _sensorLRot = GetRotMatrix(-sensorAngle);
    _turnRRot = GetRotMatrix(turnAngle);
    _turnLRot = GetRotMatrix(-turnAngle);
    
    // 回転角度を計算
    float2 newDir = _particles[id.x].direction;
    if (rightValue > centerValue && leftValue > centerValue)
    {
        newDir = random(float3(_particles[id.x].position, 0), float3(_particles[id.x].direction, 0)) < 0.5 ?
            mul(_turnRRot, _particles[id.x].direction) : mul(_turnLRot, _particles[id.x].direction);
    }
    else if (leftValue > centerValue && leftValue > rightValue)
    {
        // 左のセンサーの値の方が大きい場合
        newDir = mul(_turnLRot, _particles[id.x].direction);
    }
    else if (rightValue > leftValue && rightValue > centerValue)
    {
        // 右のセンサーの値の方が大きい場合
        newDir = mul(_turnRRot, _particles[id.x].direction);
    }
    else
    {
        newDir = mul(_turnRRot, _particles[id.x].direction);
    }
    _particles[id.x].direction = newDir;

    // パーティクルの位置を更新
    float2 newPos = _particles[id.x].position + _particles[id.x].direction * moveDistance;
    newPos = WrapPosition(newPos);
    _particles[id.x].position = newPos;

    // テクスチャにパーティクルの現在位置を記録
    uint2 depositCoord = (uint2)newPos;
    _resultTexture[depositCoord] += depositValue;
    _resultTexture[depositCoord] = clamp(_resultTexture[depositCoord], 0, 1);
}

// テクスチャへの書き込み結果を減衰させる
[numthreads(8, 8, 1)]
void DiffuseDecay(uint3 id : SV_DispatchThreadID)
{
    // パラメータをマスク用テクスチャのサンプリング結果に応じて分ける
    float decayFactor;
    if (_useInputTexture && _inputTexture[id.xy].r > 0.5)
    {
        decayFactor = _decayFactor1;
    }
    else
    {
        decayFactor = _decayFactor0;
    }
    
    // テクスチャの書き込み結果を拡散させる
    // 拡散の計算に使用する範囲を定義
    int2 sampleRage = int2(-5, 5);

    // 色の合計値
    float4 colorSum = float4(0, 0, 0, 0);
    
    // 加算するピクセルの数
    int count = 0;

    // ピクセルの色をサンプリング
    for (int x = sampleRage.x; x <= sampleRage.y; x++)
    {
        for (int y = sampleRage.x; y < sampleRage.y; y++)
        {
            uint2 coord = uint2((int)id.x + x, (int)id.y + y);

            // execute in the range of the resolution of _resultTexture
            //if (coord.x < 0 || coord.x >_resultTexRes.x || coord.y < 0 || coord.y >_resultTexRes.y) { continue; }

            colorSum += _resultTexture[coord];
            count++;
        }
    }
    
    colorSum /= (float)count;

    // 減衰を適用
    colorSum /= decayFactor;
    _resultTexture[id.xy] = float4(colorSum.xyz, 1);
}

動作結果

まずはそのまま組み込んでみます。
動画中のSimulation ParamsのElement 0の値を変えることでパターンが変化します。

　

次に画面内の表情に変化を付ける為、２次元データ上に適用するパラメータのセットを２種類用意します。
適用するパラメータの振り分けには画像をマスクとして使用します。
マスクにはは以下のような画像を使用しました。

マスク画像の黒い領域が動画中のElement0、白い領域がElement 1の割り当てに使用されます。

　

今度はマスクのテクスチャとして動的にペイントしているRenderTextureを使用します。
マウスやペンタブレット等でペイントされた線は徐々に減衰して消えていきます。
ペイントした線が残っている部分にElement 1を適用します。

参考

・jasonwebb/morphogenesis-resources
　https://github.com/jasonwebb/morphogenesis-resources#physarum

・physarum - Sage Jenson
　https://cargocollective.com/sagejenson/physarum

・Physarum Slime Mold - Entagma
　https://entagma.com/physarum-slime-mold/

・[Houdini Tutorial] #0052 Constrained Physarum Simulation - Junichiro Horikawa
　https://www.youtube.com/watch?v=1JSX4VcgrIM