e.blog

主にUnity/UE周りのことについてまとめていきます

ComputeShaderで動かす様々な形に変化するパーティクルシステム

概要

今回は1/7〜1/10にかけて開催されたCES2020でMESONが展示した『PORTAL』の中で、メインの演出となったパーティクルシステムについて書きたいと思います。

このパーティクルシステムの特徴は、指定したモデルグループの頂点位置にパーティクルをまとわせるというものです。パーティクルの色は対象モデルの頂点位置と同じ色になるように調整されています。

以下の動画で、どんな動きをするかを見ることができます。

youtu.be

今回の記事の内容はGitHubにもアップしてあるので実際の動作を見たい方はそちらもご覧ください。

github.com

サンプルプロジェクトを実行すると以下のようなシーンを見ることができます。
(※ 動画内のモデルはUnity Asset Storeのものなのでリポジトリ内には含まれません。モデルは各自でご用意ください


Transform particle system demo

常にシーンに存在しているパーティクルが設定に応じて様々な形に変化します。
メッシュの頂点に張り付いたり、任意の形状になったり、テクスチャの絵の形になったり。

今回はこれを達成するために工夫した点についてまとめます。


Table of Contents


処理フロー

まずはざっくり概観を。

  1. 対象となるモデル(メッシュ)の情報を集める
  2. それをグルーピングする(すべての情報を配列にまとめる)
  3. 初期化用データのためのインデックスバッファを用意する
  4. データをComputeShaderに送る
  5. 必要に応じてComputeShaderでターゲット位置を更新
  6. GPUインスタンシングでまとめて描画する

方針

変化しないデータをひとつにまとめる

最初、モデル(メッシュ)のデータをターゲットごとに都度取り出して処理しようとしたところ単純に負荷が高すぎてうまくいきませんでした。
特に、プロジェクトでは1モデルあたり10,000〜20,000頂点くらいあり、さらにそれを5、6体分グルーピングする必要があったためそれなりの負荷になってしまいました。

ゲーミングPCのような高性能なPCでさえもプチフリーズをするくらいには負荷が高かったです。

そこでメモリアクセスを効率よくするため、またコピーを最適化するため頂点情報をまとめてひとつの配列にしSystem.Array.Copyを利用することで対応しました。(同様にUVなどの値も別途配列にまとめる)

つまり、メモリ効率を意識してすべての情報を配列にまとめるということをしました。
テクスチャもTexture2DArrayという仕組みで配列にまとめてシェーダに送ります。

頂点データも変化しないデータとしてまとめてしまっているのは対象モデルがSkinnedMeshRendererを持っておらずアニメーションしない前提のためです。

アニメーションする場合は毎フレーム頂点位置を更新してやる必要があるため負荷が高いかもしれません。(時間があったらSkinnedMeshRendererにも対応してみたいと思います)

初期化用データを用意してパーティクルのターゲット変更をCompute Shaderで行う

各パーティクルそれぞれは常になにかしらのターゲット位置を持っています。(ランタイムで位置を計算する場合もあります)
またターゲットは、冒頭の動画のように途中で切り替えることができるようになっています。

しかしながら、パーティクル数は膨大になることが多くそれをCPUで切り替えるには重すぎる可能性があります。
そこで本パーティクルシステムでは初期化(ターゲットの変更処理)もGPUで行っています。

イメージ的には、ターゲットが変わった際にグループ単位でごっそり内容を入れ替えてしまうイメージです。

Indexバッファを用いて初期化用データ配列にアクセスする

前述のターゲット変更の仕組みですが、ここにも少し工夫があります。
パーティクルの更新処理としてはスレッドIDを配列の添字にしてアクセスすることが多いと思いますが、今回のシステムではシーンに存在しているパーティクル数とターゲットの数が異なる場合がほとんどです。

つまりパーティクルの数と初期化用データの数が異なるということです。
当然、配列の数が異なるということは同一の添字を利用することはできません。

例えば冒頭の動画を見てもらうと分かるように、移動対象がメッシュのモデルの頂点だけではなく、テクスチャのピクセルなどにもなりえます。

そこで採用したのが『インデックスバッファ』です。

イメージ的には通常のレンダリングパイプラインで利用されるインデックスバッファと同様です。

各パーティクルを初期化する際、InitData構造体の配列を用いてパーティクルの情報を書き換えるようにしているのですが、どのInitData構造体にアクセスすればいいのかはスレッドIDからでは判断できません。
そこでインデックスバッファの出番、というわけです。

IndexBuffer配列はパーティクル数と同じになっていますが、格納されているのは対象となるInitData配列のインデックスです。
このインデックス番号を利用して初期化対象のデータを特定し、初期化を行うというフローになっています。

初期化時のカーネル関数を見てみるとイメージが掴めるかと思います。

[numthreads(THREAD_NUM, 1, 1)]
void SetupParticlesImmediately(uint id : SV_DispatchThreadID)
{
    TransformParticle p = _Particles[id];

    uint idx = _IndexBuffer[id];

    float4x4 mat = _MatrixData[_InitDataList[idx].targetId];

    p.isActive = _InitDataList[idx].isActive;
    p.position = mul(mat, float4(_InitDataList[idx].targetPosition, 1.0)).xyz;
    p.targetPosition = p.position;
    p.uv = _InitDataList[idx].uv;
    p.targetId = _InitDataList[idx].targetId;
    p.scale = _InitDataList[idx].scale;
    p.horizontal = _InitDataList[idx].horizontal;
    p.velocity = _InitDataList[idx].velocity;

    _Particles[id] = p;
}

上記のuint idx = _IndexBuffer[id];の部分がインデックスバッファからインデックスを取得している箇所ですね。

解説

全体を概観できたところで、コードを交えながらどう実装したかを解説していきたいと思います。

Compute Shaderの利用

今回のパーティクルシステムではCompute Shaderを利用しています。

主なCompute Shaderのカーネルは2種類

まずはCompute Shaderについてです。
そもそもCompute Shaderってなんぞ? って人は前に書いた記事も参考にしてみてください。

edom18.hateblo.jp

edom18.hateblo.jp

今回、パーティクルの計算のために用意したのは主に2種類のカーネルです。
後述しますが、それぞれに複数のカーネルがあります。

種類は以下。

  1. パーティクルをターゲットに動かすための更新用カーネル
  2. パーティクルのターゲット自体を変更するための初期化カーネル

の2つです。

構造体

カーネルの説明に入る前に、今回のパーティクルシステムで利用している構造体を整理しておきます。
利用している構造体は以下の2つ。

ひとつめがパーティクル用でふたつめが初期化用です。

public struct TransformParticle
{
    public int isActive;
    public int targetId;
    public Vector2 uv;

    public Vector3 targetPosition;

    public float speed;
    public Vector3 position;

    public int useTexture;
    public float scale;

    public Vector4 velocity;

    public Vector3 horizontal;
}

public struct InitData
{
    public int isActive;
    public Vector3 targetPosition;

    public int targetId;
    public float scale;

    public Vector4 velocity;

    public Vector2 uv;
    public Vector3 horizontal;
}

初期化用のほうは、初期化したい内容だけを定義した構造体になっています。
ターゲット位置とターゲットのID(モデル行列やテクスチャアクセスに利用)、UV値、スケール、そしてアクティブ状態の内容を保持します。

初期化用カーネル

初期化用カーネルでは前述の通り、すでにバインド済みのパーティクルバッファを初期化用構造体によって初期化します。
また前述のように、対象のInitData配列にアクセスするためのインデックスバッファも利用します。

なお、初期化用のカーネルはふたつあります。

ひとつは、ターゲットを更新し、そのターゲットにスムーズに遷移するための初期化(SetupParticles)。
そしてもうひとつは『即座に』パーティクルの状態を変更するための初期化(SetupParticlesImmediately)です。

ひとつめがメインの初期化処理ですが、表現によっては即座にパーティクルを指定の位置に移動させたい場合があります。
その場合に利用するのが2番目のカーネルというわけです。

具体的なコードは以下。

[numthreads(THREAD_NUM, 1, 1)]
void SetupParticles(uint id : SV_DispatchThreadID)
{
    TransformParticle p = _Particles[id];

    uint idx = _IndexBuffer[id];

    float4x4 mat = _MatrixData[_InitDataList[idx].targetId];

    p.isActive = _InitDataList[idx].isActive;
    p.targetPosition = mul(mat, float4(_InitDataList[idx].targetPosition, 1.0)).xyz;
    p.uv = _InitDataList[idx].uv;
    p.targetId = _InitDataList[idx].targetId;
    p.scale = _InitDataList[idx].scale;
    p.horizontal = _InitDataList[idx].horizontal;
    p.velocity = _InitDataList[idx].velocity;

    _Particles[id] = p;
}

[numthreads(THREAD_NUM, 1, 1)]
void SetupParticlesImmediately(uint id : SV_DispatchThreadID)
{
    TransformParticle p = _Particles[id];

    uint idx = _IndexBuffer[id];

    float4x4 mat = _MatrixData[_InitDataList[idx].targetId];

    p.isActive = _InitDataList[idx].isActive;
    p.position = mul(mat, float4(_InitDataList[idx].targetPosition, 1.0)).xyz;
    p.targetPosition = p.position;
    p.uv = _InitDataList[idx].uv;
    p.targetId = _InitDataList[idx].targetId;
    p.scale = _InitDataList[idx].scale;
    p.horizontal = _InitDataList[idx].horizontal;
    p.velocity = _InitDataList[idx].velocity;

    _Particles[id] = p;
}

実は違いは一箇所だけです。p.positionを設定するかしないかです。
位置更新用のカーネルでは、呼び出されるごとに位置を更新するわけですが、ターゲットが存在するパーティクルの場合は今の位置から徐々にp.targetPositionに移動する、という実装になっています。

つまりは、p.positionが即座に変更されることで、見た目上即時、位置が更新されたように見えている、というわけですね。

さて肝心の処理ですが、基本的にはInitData構造体によって渡されたデータをコピーしているだけです。

初期化処理の中でのポイントは_MatrixData部分です。
targetIdによって計算中のパーティクルがどのモデルに属しているのかを判別しており、該当のモデル行列を取り出してターゲット位置に掛けることで実際のモデルの頂点位置を計算しています。

ちなみに感のいい方なら気づいたかもしれませんが、テクスチャのピクセルに変形する場合はモデル行列は必要ありません。そのため、その場合にはただのMaterix4x4.identityを利用しています。

また同様に、レンダリング用シェーダでも対象テクスチャを判別するためにtargetIdが必要になるため、パーティクルのデータとして設定しています。
なお、なぜこれが必要になるのかはレンダリングの詳細のときに解説します。

更新用カーネル

これは特にむずかしいことはしていません。
コードを見てもらうと分かりますが、ターゲット位置と現在位置の差分距離をDeltaTimeによって徐々にターゲット位置に近づけるように計算しているだけです。
(ただ、パーティクルごとにアニメーションの差を生ませるために『速度』の項目がありますが、本質的にはターゲットに近づけていることには変わりありません)

[numthreads(THREAD_NUM, 1, 1)]
void UpdateAsTarget(uint id : SV_DispatchThreadID)
{
    TransformParticle p = _Particles[id];

    float3 delta = p.targetPosition - p.position;
    float3 pos = (delta + p.velocity.xyz * 0.2) * _DeltaTime * p.speed;

    const float k = 5.5;
    p.velocity.xyz -= k * p.velocity.xyz * _DeltaTime;

    p.position += pos;
    p.useTexture = 1;

    _Particles[id] = p;
}

特定の形状に移動させるカーネル

基本的には設定されたモデルの頂点にパーティクルが張り付く動作をしますが、パーティクルの移動先をランタイムに計算することで任意の形状に散らすこともできます。

サンプルの動画では拡散したりぐるぐる周ったりしているのが分かるかと思います。

該当のコードは以下の通り。

[numthreads(THREAD_NUM, 1, 1)]
void UpdateAsExplosion(uint id : SV_DispatchThreadID)
{
    TransformParticle p = _Particles[id];

    float3 pos = (p.velocity.xyz) * p.velocity.w * _DeltaTime;

    float s = sin(rand(id) + _Time) * 0.00003;
    p.velocity.xyz += s;
    float k = 2.0;
    p.velocity.xyz -= k * p.velocity.xyz * _DeltaTime;

    p.position += pos;

    p.useTexture = 0;

    _Particles[id] = p;
}

パーティクルシステムでは更新方法のカーネルを変えることでパーティクルの位置計算の仕方を変更しています。
つまり、カーネルを増やせば様々な形状の変化を与えることができるわけです。

なので例では四散するパーティクルの更新用カーネルを示しましたが、サンプルの動画ではぐるぐると周る演出もあります。
これは更新用カーネルを変えることによって実現しています。

この更新用カーネルを増やすことでさらに様々な表現を作ることが可能になっています。

C#コード

次にC#側のコードを見ていきましょう。

今回のシステムで使う主なクラスは以下の3つです。

  • TransformParticleSystem
  • ParticleTarget
  • ParticleTargetGroup

各クラスの概要は以下の通りです。

TransformParticleSystem

TransformParticleSystemが今回のコア部分です。このクラスでCompute Shaderへのデータ設定や各種データの取得、整理などを行っています。

ParticleTarget

ParticleTargetはターゲットとなるMeshの各種情報を取得するためのラッパークラスです。
各種情報とは、頂点やUV、テクスチャ情報などパーティクルシステムで利用する情報です。

この派生型のクラスとして、テクスチャのピクセル位置に移動させるようなクラスも用意されています。
ただ基本概念は一緒なので説明は割愛します。

ParticleTargetGroup

ParticleTargetGroupParticleTargetをグループ化するためのクラスです。
今回のプロジェクトの要件では複数のモデルをひとつにまとめて扱う必要があったためこのグループクラスを定義しています。

このグループクラスでは、子どもに持つParticleTargetから頂点などの情報を取得しそれらをひとつにまとめる処理を担当しています。
最終的にこのグループ単位でパーティクルシステムにデータが渡され利用されます。

C#側では主にバッファの設定とCompute Shaderへの計算命令の実行(Dispatch)を行います。

バッファ周りの仕組みや概要については以前の記事を参考にしてください。

edom18.hateblo.jp

データの設定

今回のポイントはデータの変更処理にあります。

特にマトリクスの変更とターゲット位置の変更がメインです。
ということでマトリクスの変更処理部分を見てみます。

private void UpdateMatrices(ParticleTargetGroup group)
{
    group.UpdateMatrices();

    System.Array.Copy(group.MatrixData, 0, _matrixData, 0, group.MatrixData.Length);

    _matrixBuffer.SetData(_matrixData);
}

// ParticleTargetGroup.UpdateMatrices処理
public void UpdateMatrices()
{
    for (int i = 0; i < _targets.Length; i++)
    {
        _matrixData[i] = _targets[i].WorldMatrix; // localToWorldMatrixを返すプロパティ
    }
}

グルーピングされているターゲットのlocalToWorldMatrixを配列に詰めているだけです。

続いてターゲット位置の変更。

public void UpdateInitData(InitData[] updateData)
{
    if (updateData.Length > _initDataList.Length)
    {
        Debug.LogError("Init data list size is not enough to use.");
    }

    int len = updateData.Length > _initDataList.Length ? _initDataList.Length : updateData.Length;

    System.Array.Copy(updateData, _initDataList, len);

    _initDataListBuffer.SetData(_initDataList);
}

このInitDataが更新用データ構造体で、パーティクルひとつの更新情報になります。
このInitData配列を、各グループクラスが起動時にまとめあげて保持しているわけです。

そして変更のタイミングに応じて、すでに生成済みのInitData配列を渡すことで高速にターゲットを変更しているわけです。
なのでここをランタイムで生成して渡すことで完全に任意の位置にパーティクルを飛ばすことも可能になっています。

そしてこのあとに初期化用カーネルを呼び出してパーティクルの更新を行っています。

呼び出し部分は以下。

private void UpdateAllBuffers(int kernelId)
{
    SetBuffer(kernelId, _propertyDef.ParticleBufferID, _particleBuffer);
    SetBuffer(kernelId, _propertyDef.InitDataListID, _initDataListBuffer);
    SetBuffer(kernelId, _propertyDef.MatrixDataID, _matrixBuffer);
    SetBuffer(kernelId, _propertyDef.IndexBufferID, _indexBuffer);
}

private void Dispatch(int kernelId)
{
    _computeShader.Dispatch(kernelId, _maxCount / THREAD_NUM, 1, 1);
}

それぞれ必要なバッファを対象カーネルにセットし、そのあとでカーネルを起動します。
呼び出すカーネルは前述した初期化用カーネルです。

これを呼び出すことでパーティクルデータのバッファの内容が書き換わり、以後の位置更新ループによって別のターゲットへパーティクルが近づいていくことになります。

パーティクルのレンダリングについて

本パーティクルシステムの大まかな仕組みは当初から変わっていません。
しかしモック段階では問題にならなかったものも、実際のコンテンツに組み込むことで問題が出てくることが多々あります。

今回も同じように問題が発生しました。

問題というのは、モック段階ではパーティクルの移動のみに焦点を当てて実装を行っていたのですが、いざ組み込みを開始すると、パーティクルの移動先の対象モデルと同じ色にしないとなりません。
(でないとモデル形状になるだけの単色のパーティクルの塊になってしまう)

当然、対象モデルはテクスチャを持っており、頂点ごとに色など持っていません。
つまり、対象頂点位置の色をテクスチャから引っ張ってくる必要があるわけです。

しかし、パーティクルに利用するシェーダは当然ひとつだけで、GPUに転送できるデータにも限りがあります。
CPUのように柔軟にテクスチャ情報を取り出すわけにはいきません。
なにより、転送できるテクスチャの数にも限りがある上に、複数のテクスチャを配列以外でアクセスするには問題が多すぎました。

そこで、最初の『方針』のところでも書いたように、対象モデルのテクスチャをTexture2DArrayにまとめることで配列にし、かつ『どのテクスチャのどこをフェッチするか』という情報をパーティクル自体に持たせる必要があります。

ということを念頭に置いて、どうしているかのコード断片を抜き出すと以下のようになります。

// 頂点シェーダ
TransformParticle p = _Particles[id];

v2f o;

// 中略

o.uv.xy = p.uv.xy;
o.uv.z = p.targetId;

// フラグメントシェーダ
col = UNITY_SAMPLE_TEX2DARRAY(_Textures, i.uv);

Texture2DArrayでは最初の2要素(x, y)が通常のUV値として利用され、3要素目(z)が配列の添字として利用されます。
なのでそれを頂点シェーダからフラグメントシェーダへ転送し、UNITY_SAMPLE_TEX2DARRAYマクロを利用して対象の色をフェッチしている、というわけです。

正直、これが想定通りうまく行ったときは内心飛び跳ねていました。
移動自体は問題なく動いたものの、色味をつけるところで躓いていたので。

この仕組み自体はWindows, Mac OSX, iOS, Android, MagicLeapすべてで問題なく動いているので汎用的に使える仕組みだと思います。

なお、パーティクルを画面に表示するフローについては以前書いた以下の記事を参照ください。このパーティクルシステムを実装したあとに備忘録として書いた記事なのでほぼ内容はそのままです。

edom18.hateblo.jp

プロジェクトでのハマりポイント

最初のモック段階では問題なく動いていたものの、実際のプロジェクトでは画面になにも表示されないという問題が発生しました。

おそらく理由はメモリ不足。
ターゲットとなるモデルを複数指定していたのですが、同じモデルなのにも関わらず複数のグループに分割したために必要となるデータの重複が発生し、特にTexture2DArrayにまとめたテクスチャのデータ量が増えすぎたためだと思います。

なので参考にアップしてあるGitHubのコードではグループ化の最適化を行っています。
具体的には、サブグループという概念を取り入れてグループデータはひとつにまとめ、サブグループ単位でパーティクル位置を指定できるようにしました。

ただ今回の解説からは脱線するので説明は割愛します。

逆引きUniRx - 使用例から見る使い方

概要

UniRx。最近よく目にする気がしますが基本的な概念は『Reactive Extensions』です。
似た用語として『関数型リアクティブプログラミング:Functional Reactive Programming(FRP』がありますが概念としては違うようです。

ざっくりとRx系ライブラリを説明すると、連続した値を『ストリーム』と捉え、それをどう扱うかに焦点を当てたものです。
ストリームのイメージは、『なにか』がパイプの中を流れていくイメージです。
この『なにか』はデータであるかもしれないし、時間かもしれません。とにかく抽象化されたものがパイプ内を流れ、それに加工を加えながら処理するもの、という感じです。

今回はUniRxの使い方の説明は割愛します。使い方や基礎的なところは@toRisouPさんがとても詳しく記事を書いてくださっているのでそちらを見るといいでしょう。

qiita.com


今回書くのは、ある程度基礎が分かっていて概念も把握したものの『で、実務でどう使ったらいいんだ?』となった人向けに、具体的な事例を交えて逆引き的に利用できるようにまとめたものです。
(というか、完全に自分のためのメモです・・・w なので、今後使いやすい・思い出しやすいように逆引きで書いてるというわけです)

なので多分、随時更新されていくと思います。

ドラッグ処理

最近実装したものでドラッグ処理です。
例で使われている_raycasterはレイが当たっているかどうかを判定するためのクラスで、_inputControllerはコントローラのトリガー状態を提供してくれるものです。

例がだいぶ偏ったものになっていますが、VRでコントローラからレイを飛ばしてなにかをドラッグして加速度を計算、対象オブジェクトを移動させる、みたいなシーンを思い浮かべてください。

var startStream = this.UpdateAsObservable()
                      .Where(_ => _raycaster.IsHit && _inputController.IsTriggerDown);

var stopStream = this.UpdateAsObservable()
                     .Where(_ => !_raycaster.IsHit || _inputController.IsTriggerUp);

startStream
    .SelectMany(x => this.UpdateAsObservable())
    .TakeUntil(stopStream)
    .Select(_ => _raycaster.ResultRaycastHit.point)
    .Pairwise()
    .RepeatUntilDestroy(this)
    .Subscribe(Dragging)
    .AddTo(this);


private void Dragging(Pair<Vector3> points)
{
    Vector3 cur = _anyObj.transform.worldToLocalMatrix.MultiplyPoint3x4(points.Current);
    Vector3 prev = _anyObj.transform.worldToLocalMatrix.MultiplyPoint3x4(points.Previous);

    float velocity = (cur.x - prev.x) / Time.deltaTime;
    _acceleration += velocity / Time.deltaTime;
    _acceleration *= _attenuateOfAcceleration;

    float relativeVelocity = (_acceleration * Time.deltaTime * _coefOfVelocity) - _anyObj.Velocity;

    _anyObj.AddVelocity(relativeVelocity);
}

キモは以下の部分です。

startStream
    .SelectMany(x => this.UpdateAsObservable())
    .TakeUntil(stopStream)
    .Select(_ => _raycaster.ResultRaycastHit.point)
    .Pairwise()
    .RepeatUntilDestroy(this)
    .Subscribe(Dragging)
    .AddTo(this);

ここで行っていることは、『なにがしかのスタートタイミング(startStream)』から開始され、指定の終了ストリーム(stopStream)に値が流れてくるまで継続する、というものです。

そしてストリームが継続している間はレイのヒット位置をストリームに流し(Select)、それをペアにし(Pairwise)、オブジェクトが破棄されるまで継続する、というものです。

値が閾値を越えたら処理をする

例えば、速度が一定速度以上になり、またそれが一定速度以下になった、という閾値またぎを検知したい場合があるかと思います。
そんなときに利用できるのがDistinctUntilChangedフィルタです。

これは、同じ値が連続している間はその値を流さない、という動作をします。
つまり、特定の値(今回の例では速度)が閾値をまたいだときにtrue/falseを返すようにしておき、それをフィルタすることで閾値をまたいだことを検知することができます。

ちなみに以下の例でSkip(1)が入っているのは初期化時など最初に発生してしまうイベントを無視するために入れています。

this.UpdateAsObservable()
    .Select(_ => Mathf.Abs(Velocity) <= _stopLimit)
    .DistinctUntilChanged()
    .Where(x => x)
    .Skip(1)
    .Subscribe(_ => DoAnything())
    .AddTo(this);

一定時間経過したら無効化する

次はボタンなど『開始イベント』と『終了イベント』があり、かつ制限時間を設けたい場合の処理です。

以下のサンプルではまず、開始イベントにスペースキーのDown、終了イベントにスペースキーのUpを設定しています。
そしてさらに、制限時間(例では3秒)が経過した場合も終了するようになっています。

var startStream = this.UpdateAsObservable().Where(_ => Input.GetKeyDown(KeyCode.Space));
var stopStream = this.UpdateAsObservable().Where(_ => Input.GetKeyUp(KeyCode.Space));
var timeOut = Observable.Timer(System.TimeSpan.FromMilliseconds(3000)).Select(_ => false);

startStream
    .SelectMany(stopStream.Select(_ => true).Amb(timeOut ).First())
    .Where(x => x)
    .Subscribe(_ => DoHoge())
    .AddTo(this);

ここでの大事な点はAmbです。
Ambはストリームを合成し、どちらかのストリームに流れたものをそのままひとつのストリームとして流してくれるオペレータです。

なのでここでは『スペースキーUp』と『制限時間経過』を『終了イベント』として捉え、そのどちらかが流れたら終了するようにしています。

ボタン長押しを検知

次はシンプルな『ボタン長押し』の処理です。

1秒後に発火

まず最初は指定時間押し続けていたら発火するもの。

var clickDownStream = this.UpdateAsObservable().Where(_ => Input.GetKeyDown(KeyCode.Space));
var clickUpStream = this.UpdateAsObservable().Where(_ => Input.GetKeyUp(KeyCode.Space));

clickDownStream
    .SelectMany(_ => Observable.Interval(System.TimeSpan.FromSeconds(1)))
    .TakeUntil(clickUpStream)
    .DoOnCompleted(() =>
    {
        Debug.Log("Completed!");
    })
    .RepeatUntilDestroy(this)
    .Subscribe(_ =>
    {
        Debug.Log("pressing...");
    });

押している間、押下を検知

上記は『長押し』だけを検知するものでしたが、こちらは『押している間』のイベントも受け取れるようにしたものです。
使用想定としては、押している間だけ常に何か処理をする、みたいなケースです。

var clickDownStream = this.UpdateAsObservable().Where(_ => Input.GetKeyDown(KeyCode.Space));
var clickUpStream = this.UpdateAsObservable().Where(_ => Input.GetKeyUp(KeyCode.Space));

clickDownStream
    .SelectMany(_ => this.UpdateAsObservable())
    .TakeUntil(clickUpStream)
    .DoOnCompleted(() =>
    {
        Debug.Log("Completed!");
    })
    .RepeatUntilDestroy(this)
    .Subscribe(_ =>
    {
        Debug.Log("pressing...");
    });

一定時間経つ前にボタンが離された場合も処理

こちらはUnity開発者ギルドの質問チャンネルで質問した際に教えていただいた方法です。
上記では『押している間』のイベントを捉えることができましたが、『終了判定』は取れませんでした。
『終了判定』を加えたものが以下のものです。

var startStream = this.UpdateAsObservable().Where(_ => Input.GetKeyDown(KeyCode.Space));
var stopStream = this.UpdateAsObservable().Where(_ => Input.GetKeyUp(KeyCode.Space));
var timeOut = Observable.Timer(System.TimeSpan.FromSeconds(5)).AsUnitObservable();

startStream
    .SelectMany(stopStream.Amb(timeOut))
    .First()
    .RepeatUntilDestroy(this)
    .Subscribe(_ =>
    {
        Debug.Log("hoge");
    });

値の変化監視(true / false)

値の変化監視はよくあるニーズだと思います。
例えば、前回はfalseだったものがtrueになったときだけ処理したい、などです。

using UniRx.Triggers; // ←UpdateAsObservableを使うにはこれが必要

this.UpdateAsObservable()
    .Select(_ => IsHoge())
    .DistinctUntilChanged()
    .Where(x => x)
    .Subscribe(_ =>
    {
        Debug.Log("Is Hoge");
    });

ObserveEveryValueChangedを使ったほうがもっとシンプルに書ける

this.ObserveEveryValueChanged(x => x.IsHoge())
    .Where(x => x)
    .Subscribe(_ =>
    {
        Debug.Log("Is Hoge");
    });

UniRxでコルーチンを使ったアニメーション

こちらはコルーチンを交えてUniRxでアニメーションを行う例です。
Observable.FromCoroutine<T>を使うことでコルーチンをストリームに変え、かつコルーチン内で計算した結果を受け取ることができます。

使い方は、`に渡すラムダの引数にObservableが渡ってくるのでそれをコルーチンに渡し、そのObservableを介してOnNext`を呼んでやることで実現しています。

private void DoAnimation()
{
    Observable.FromCoroutine<float>(o => Animation(o, duration))
        .SubscribeWithState3(transform, transform.position, position,
        (t, trans, start, goal) =>
        {
            trans.position = Vector3.Lerp(start, goal, t);
        })
        .AddTo(this);
}

private IEnumerator Animation(IObserver<float> observer, float duration)
{
    float timer = duration;

    while (timer >= 0)
    {
        timer -= Time.deltaTime;

        float t = 1f - (timer / duration);

        observer.OnNext(t);

        yield return null;
    }
}

参考記事

以下は参考にした記事のリンク集です。

qiita.com

noriok.hatenadiary.jp

developer.aiming-inc.com

qiita.com

qiita.com

www.slideshare.net

rxmarbles.com

Generalized Perspective ProjectionをUnityで実装する

はじめに

元記事は、以前自分が英語で書いたものの日本語訳版になります。(自分で書いた英語を翻訳するという初体験w)

medium.com

概要

今回のこの記事は、Generalized Perspective ProjectionをUnityで実装するための解説記事です。
これが意味するのは、任意の視点からのGeneralized Perspective Projection用のマトリクスを生成する方法を示します。

本実装を行うにあたって以下の論文を参考にさせていただきました。

drive.google.com


以下の動画を見ると今回の記事でなにができるかが分かるかと思います。

まず、動画ではカメラに対して垂直になるように回転とプロジェクションをプレーンに適用するところを示しています。

www.youtube.com

次の動画は任意の視点から、例えば『窓の外』を映し出すようなプレーンの様子です。

www.youtube.com

今回の動画の実装は、GitHub上にアップしてあるので、実際の挙動を見たい方は以下からcloneしてご覧ください。

github.com

イデア

今回のメインアイデアは、カメラ(ビュー)から垂直に位置するプレーンに対してプロジェクションする、という方法をベースにしています。

それを実現するためにいくつかの点(ポイント)と直行するベクトルを求める必要があります。

コーナーポイントの定義

最初に、スクリーンとなるプレーンのコーナーポイントを定義します。

以下の画像は、各ポイントがどう配置されるかを示しています。

f:id:edo_m18:20200621154717p:plain

Pa, Pb, Pcの3点を定義しています。

平面上の直行ベクトルを計算する

次に、平面上の直行するベクトルを計算します。
計算自体はとてもシンプルです。すでにコーナーポイントがあるので、これを用いてベクトルを求めます。

計算は以下のようになります。

Vu = (pc - pa).normalized;
Vr = (pb - pa).normalized;
Vn = Vector3.Cross(vu, vr).normalized;

これを図示すると以下になります。

f:id:edo_m18:20200621155143p:plain

中心から外れたポイントを計算する

まず、通常の視錐台の状態を確認しておきましょう。

f:id:edo_m18:20200621155249p:plain

これは、通常の視錐台の注視点がどこにあるかを図示しています。
見て分かる通り、平面の中心に位置していることが分かります。

しかし今回は任意のビューからの視錐台を作る必要があります。
つまり、中心点が動かされた点(Peの正面にある点)を求める必要があります。
図にすると以下のようになります。

f:id:edo_m18:20200621155456p:plain

当然、この点を求める必要があります。が、すでに計算に必要な情報は揃っているので、それを使って求めます。

まず、カメラ位置から各コーナーポイントへのベクトルを求めます。

Vector3 va = pa - pe; 
Vector3 vb = pb - pe;
Vector3 vc = pc - pe;

これを図示すると以下のようになります。

f:id:edo_m18:20200621155812p:plain

視錐台の各値を求める

以上までで必要な情報が揃いました。
これらを用いて視錐台の情報(l,r,t,b)を求めます。

各変数の意味は以下の図の通りです。

f:id:edo_m18:20200621155957p:plain

これらの値は以下のように求めることができます。

// Find the distance from the eye to screen plane.
float d = -Vector3.Dot(va, vn);

// Find the extent of the perpendicular projection.
float nd = n / d;

float l = Vector3.Dot(vr, va) * nd;
float r = Vector3.Dot(vr, vb) * nd;
float b = Vector3.Dot(vu, va) * nd;
float t = Vector3.Dot(vu, vc) * nd;

さて、ここで登場したdがなんだろうと思ったと思います。
このdはスケールです。near planeまでの距離nに比べてどれくらい差があるか、を示しています。

図にしたほうが分かりやすいと思うので、図示すると以下のようになります。

f:id:edo_m18:20200621160205j:plain

図には2つの三角形が描かれているのが分かるかと思います。
これらは相似になっています。相似ということはひとつの比率から、その他の辺の比率も求められる、ということです。

今ほしいのは視錐台を生成するためのパラメータです。
そして前述のようにl,r,t,bを求めたいわけです。ただこれはnear plane想定の値である必要があるため、カメラと平面との距離を測り、その比率を掛けたものが、最終的に求めたい値となります。

なので、前述のコードではそれぞれの値にnd(スケール)を掛けていた、というわけです。

さて、これで視錐台を生成するための準備が整いました。

コード解説

前述までの結果を元に視錐台を生成します。まずはコード全文を見てください。

private void UpdateFrustrum()
{
    float n = _camera.nearClipPlane;
    float f = _camera.farClipPlane;
    Vector3 pa = _pa.position;
    Vector3 pb = _pb.position;
    Vector3 pc = _pc.position;
    Vector3 pe = _pe.position;
    // Compute an orthonormal basis for the screen.
    Vector3 vr = (pb - pa).normalized;
    Vector3 vu = (pc - pa).normalized;
    Vector3 vn = Vector3.Cross(vu, vr).normalized;
    // Compute the screen corner vectors.
    Vector3 va = pa - pe;
    Vector3 vb = pb - pe;
    Vector3 vc = pc - pe;
    // Find the distance from the eye to screen plane.
    float d = -Vector3.Dot(va, vn);
    // Find the extent of the perpendicular projection.
    float nd = n / d;
    float l = Vector3.Dot(vr, va) * nd;
    float r = Vector3.Dot(vr, vb) * nd;
    float b = Vector3.Dot(vu, va) * nd;
    float t = Vector3.Dot(vu, vc) * nd;
    // Load the perpendicular projection.
    _camera.projectionMatrix = Matrix4x4.Frustum(l, r, b, t, n, f);
    _camera.transform.rotation = Quaternion.LookRotation(-vn, vu);
}

※ 論文と比べると外積を求める方法が若干異なりますが、これはUnityが採用している座標系によるものです。

論文からの変更点

さて、実は今回の実装は、論文のものと比べて少し変更が入っています。
論文では最終的なマトリクス計算まで説明されています。(P, M, Tマトリクス)

論文でのTマトリクスは『カメラがどこにいるか』を示しています。
が、UnityはCGのレンダリングエンジンを持っているので、カメラの位置自体がすでに位置情報を持っています。
そのため、Unityでの実装ではTマトリクスの計算は必要ありません。

次に、Mマトリクスはすべての頂点をビューの前に移動させる行列です。言い換えると、すべての頂点はビューの正面に置かれます。
これが意味するところは、その行列はカメラを回転させることに他なりません。

ということで、これを達成するためにはUnity APIQuaternion.LookRotation)を使ってカメラを回転します。

Quaternion.LookRotationはふたつの引数を取ります。ひとつめはカメラのフォワードベクトル(つまり-vn)です。
そしてふたつめは上方向ベクトル(つまりvn)です。結果として、プレーンに垂直になるためのカメラの回転を求めることができます。

視錐台に入ったオブジェクトの検知

ここからは余談ですが、今回作成した視錐台を用いることで、この歪んだ視錐台でも、その視錐台内に入ったオブジェクトを検知することができます。

それには、Unityが提供してくれているAPIを持ちて以下のように達成することができます。

private void CheckTargets()
{
    Plane[] planes = GeometryUtility.CalculateFrustumPlanes(_camera);
    foreach (var col in _colliders)
    {
        col.GetComponent<Renderer>().material.color = Color.white;
        
        if (GeometryUtility.TestPlanesAABB(planes, col.bounds))
        {
            col.GetComponent<Renderer>().material.color =  Color.blue;
        }
    }
}

結果は以下の通りです。

www.youtube.com

上記コードで重要なのはGeometryUtility.CalculateFrustumPlanesGeometryUtility.TestPlanesAABBメソッドです。

これらは、視錐台内に入ったオブジェクトの検知に使われます。
上記動画を観てもらうと分かる通り、歪んだ視錐台でも正常に動作しているのが分かるかと思います。

ズームも行える

n/dがスケールを意味することは前述の通りです。そしてこれに特定のスケールを掛けることで視錐台をスケールさせズーム表現に使うこともできます。

www.youtube.com

最後に

自分の英語記事を翻訳するという新しい経験をしましたw

今回、この記事を日本語化しようと思った理由は、今のプロジェクトで使ってみて意外と有用だなーと感じたためです。
しかも歪ませた視錐台でも、オブジェクト検知などにそのまま使えるのは面白いですね。

例えば、プレイヤー視点から見える窓の外にオブジェクトが見えたら、みたいなことにも使えるかもしれません。

ただこれ、VRで利用するとカメラの設定が異なるせいなのか分かりませんが、若干見え方がずれることがあります・・。
このあたりについてはどなたか詳細ご存知の人がいたら教えてもらえるとうれしいです;

Unityで流体シミュレーション ~粒子法編~

はじめに

以前に、格子法を用いた2次元の流体シミュレーションについて概念編実装編に分けて2つの記事を書きました。

しかし格子法は、その定義領域の範囲を越えてしまうと計算ができなくなってしまうという問題があります。特に今回は3Dに拡張し、AR/VR内で利用することを目論んでいるのでうまく行かない可能性があります。
そこで今回は別の手法である『粒子法』を用いて流体シミュレーションを行う方法について調べてみたのでそれをまとめていきます。

今回の実装は特に、Unity Graphics Programming vol.1に掲載されているSPH法の解説とGitHubにアップされているサンプルを元にしました。

ここで書くことはあくまで自分の理解のためのメモであり、物理的・数学的見地からの詳細を解説しているものではありません。
また、多くの部分を書籍からの引用に頼っています。予めご了承ください。


Unity Graphics Programmingシリーズは無料でダウンロードできます

なお、上記書籍のシリーズは無料でPDFにて公開されています。
vol.4まであるので読み応えがあります。元の情報を参照されたい方は以下からダウンロードしてください。

github.com

これを元に実装してみたのが以下です。


Table of contents

粒子法の概要

以前書いた格子法と今回の粒子法では『視点の違い』があります。

格子法の視点は格子で、計算対象を格子に区切って計算を行います。
一方粒子法の視点は粒子で、計算対象を粒子そのものにして計算を行います。

ざっくりとしたイメージで言うとそれぞれの違いは以下です。

  • 格子法は観測地点(各格子)に定点カメラを設置してその『場』の動きを計算します。
  • 粒子法では粒子を追跡するカメラを用意して計算します。

そしてこれらの違いをそれぞれ、オイラー的視点ラグランジュ的視点と呼びます。

オイラー的視点が格子法、ラグランジュ的視点が粒子法です。

参考にさせていただいた書籍から図を引用させていただくと以下のようなイメージの違いになります。

f:id:edo_m18:20200523194416p:plain

左の図が示している丸は粒子ではなく、各格子の現在の速度場ベクトルを表しています。(なので各点間の距離が均一になっています)
一方、右の図は粒子そのものを観測しているため各点自体が現在の粒子の位置を表しています。(なので各点間の距離がバラバラになっています)

オイラー的視点とラグランジュ的視点での微分の違い

オイラー的視点とラグランジュ的視点では微分の演算の仕方が異なるようです。
なぜ突然微分の話なんだ、と思うかもしれませんが、流体の方程式は一発で粒子の位置を求めることはできず、微分して得られた速度を元に徐々に更新していくことで求めます。
そのため現在と過去の状態から微分によって速度を求め、次の位置を求める必要があります。なので微分の話なのです。

そして違いをとてもざっくり書いてしまうと、物理量に対して以下の式の違いがあります。

オイラー的視点での微分


\phi = \phi(\vec{x}, t)

これは、時刻tにおける位置 \vec{x}の物理量を表しています。これを時間微分すると以下の式が得られます。


\frac{\partial \phi}{\partial t}

なんてことはない、偏微分の式ですね。

ラグランジュ的視点での微分

一方、ラグランジュ的視点では以下の式が物理量を表します。


\phi = \phi(\vec{x}(\vec{x}_0, t), t)

これを微分すると、最終的に以下の式が得られます。


\begin{array}{c}
\lim _{\Delta t \rightarrow 0} \frac{\phi\left(\vec{x}\left(\vec{x}_{0}, t+\Delta t\right), t+\Delta t\right)-\phi\left(\vec{x}\left(\vec{x}_{0}, t\right), t\right)}{\Delta t} \\
=\sum_{i} \frac{\partial \phi}{\partial x_{i}} \frac{\partial x_{i}}{\partial t}+\frac{\partial \phi}{\partial t} \\
=\left(\left(\begin{array}{c}
u_{1} \\
u_{2} \\
u_{3}
\end{array}\right) \cdot\left(\begin{array}{c}
\frac{\partial}{\partial x_{1}} \\
\frac{\partial}{\partial x_{2}} \\
\frac{\partial}{\partial x_{3}}
\end{array}\right)+\frac{\partial}{\partial t}\right) \phi
\end{array}

※ この式が出てくる理由は、おそらく全微分によるものだと思います。全微分のざっくりとした説明はこちら


余談

微分について上記記事から引用させてもらうと、

接線の傾き(つまり、関数を1次関数で近似したときの傾き)が微分です。すなわち、
 f(x_0+h)=f(x_0)+ah+o(h)

とあります。
つまり、とある関数に対して x_0付近で1次関数で近似したものということですね。(接線の傾きということは直線=1次関数)

そして全微分の係数(a)は偏微分で求まります。(同様に、多数変数の場合はbcと変数分、係数が増え、それぞれの変数の偏微分で求まります)
(例えば a = \frac{\partial f}{\partial x}

これを用いて2変数の場合は


df = \frac{\partial f}{\partial x}dx + \frac{\partial f}{\partial y}dy

となります。このことから、偏微分の値を \sumで合計している理由が分かりますね。
おそらく \frac{\partial x_i}{\partial t}は、時間による偏微分なのでとある軸方向の微小量を求めているのだと思います。
なので上記式で言うところのdxっていうことですね。


閑話休題

さて、両者の違いが分かりますか?

オイラー的視点の式が意味するのは『位置が固定である』という事実です。なので式の中の位置ベクトル( \vec{x})はただの変数になっていますね。

一方、ラグランジュ的視点の式が意味するのは『位置は変動する』という事実です。なので位置ベクトルは時間の関数( \vec{x}(\vec{x}_0, t))になっています。

少しだけ補足すると、とある時間の位置ベクトル( \vec{x})は初期位置( \vec{x}_0)と時間tを引数に取る関数となっています。

粒子法のナビエ・ストークス方程式

流体シミュレーションをしようとすると必ず登場するナビエ・ストークス方程式
これを解くことで流体をシミュレーションすることができるようになります。つまり、流体シミュレーションにおいて一番重要な部分です。

粒子法のナビエ・ストークス方程式は以下で与えられます。
実際に式を見てもらうと分かりますが、微分の嵐ですねw
このことからも微分が重要なことが分かります。


\frac{D \vec{u}}{D t}=-\frac{1}{\rho} \nabla p+\nu \nabla^2 \vec{u}+\vec{f}

一方、格子法は以下でした。


\frac{\partial \vec{u}}{\partial t}=-(\vec{u} \cdot \nabla) \vec{u} -\frac{1}{\rho} \nabla p +\nu \nabla^{2} \vec{u}+\vec{f}

前の記事で解説した格子法によるナビエ・ストークス方程式と見比べてみると移流項 -(\vec{u} \cdot \nabla) \vec{u})がまるまるなくなっているのが分かるかと思います。
これは、格子法では観測点を『』に固定して観測している一方で、粒子法は『粒子自体』を観測しているため、そもそも『移流』が発生しないことによります。

ここの導出に関してはちょっとまだしっかり理解しきれていません。詳細については参考にした書籍をご覧ください。

そこでの解説を少しだけ引用させていただくと以下のように記載されています。

前章の格子法で出てきたNS方程式とは少し形が異なりますね。移流項がまるまる抜けてしまっていますが、先程のオイラー微分ラグランジュ微分の関係を見てみると、うまくこの形に変形できることが分かります。粒子法では観測点を流れに沿って移動させますから、NS方程式計算時に移流項を考慮する必要がありません。移流の計算はNS方程式で算出した加速度をもとに粒子位置を直接更新することで済ませる事ができます。

ここ、ちゃんと理解しきれていないのですが、微分の式を見てみるとちょうど -(\vec{u} \cdot \nabla) \vec{u}と同じ形の式が出来ていることが分かります。もしかしたらこれが相殺されてなくなる、ということなのかもしれません。

流体の非圧縮性

実はナビエ・ストークス方程式だけでは解を求めることはできません。これについては前回の「Unityで流体シミュレーションを実装する 〜 概念編 〜」の未知数と式の数のところでも触れていますが、未知数の数に対して式の数が合わないためです。

そのときの文を引用すると、

実は冒頭の式だけではこれを解くことができません。 理由は、未知数が4つに対して式が3つしかないためです。式がひとつ足りないんですね。

そこで、一般的な流体の持つ非圧縮性という性質を利用してもうひとつ式を作ります。

ということで式がもうひとつ必要になるわけです。
そこで出てくるのが見出しの『流体の非圧縮性』です。

日常観測できる範囲の液体は非圧縮性を考慮しても問題ないようです。
参考にした書籍では『音速よりも十分に小さい場合』と記載があり、やはり相当特殊な環境下でない限りは非圧縮性を持つものとして考えていいと思います。

非圧縮性は以下の式で表されます。


\nabla \cdot \vec{u} = 0

これは『発散』といい、流体内で湧き出しや消失がない(= 0)ことを意味します。なので非圧縮性。
要は、液体が押されてとある場所に移動した場合、その量と同じ分だけそこにあった液体が別のところに押し出される、という感じでしょうか。
(もしそうならないなら圧縮されていることに他ならないからです)

粒子法の計算の概念

ナビエ・ストークスの導出などは自分の知識ではまだまだできませんので、細かいことについては独自で調べてみてください。
以前に書いた格子法の概念編では色々参考にしたリンクを掲載してるのでそちらも見てみると理解が深まるかもしれません。

ここでは粒子法の計算概念について簡単に書いておきたいと思います。

いったん理論は忘れて、実際の流体に目を向けてみましょう。
みなさんご存知の通り、実際の流体も『分子』によって形成されています。つまり人の目には見えないだけで、実際は細かい粒子が集まってできているというわけですね。

つまり実際の流体も、極論で言えば超超超精細なパーティクルシステムと見なすことができます。
が、現在の科学ではそれらを真面目に計算することなど到底できません。
なので粒子法ではこれを、現在のコンピュータを用いて現実的な時間で計算できるくらいの数に減らすことを考えます。

簡単に言えば、ある程度大きな『粒』として分子をまとめて塊として扱う、ということです。
そして『人間の目に知覚できるくらいの粒子』として見た場合、いわゆる剛体の運動方程式 m\vec{a} = \vec{f}と同じことを各粒子に課すことで粒子を移動させ流体を表現できる、というわけです。(と理解しています)

この『ある程度の塊』として見た場合の粒子では、それぞれ、質量 m、位置ベクトル \vec{x}、速度ベクトル \vec{v}、体積 Vを持つと考えることができます。

流体とは、この粒子が様々な力を受けて相互干渉し、移動することで実現されます。
この様々な力は以下に挙げられるものがあります。

  • 外力(重力など)
  • 圧力
  • 粘性力

個別に見ていきましょう。

外力

これが一番分かりやすい力でしょう。
重力は言わずもがな、なにか剛体などがぶつかって力を加えたり、あるいは遮蔽物による反発力などもあるでしょう。そうした『外から与えられる力』がこの外力です。

圧力

圧力は『気圧』『水圧』などの単語からも分かるように、流体同士が及ぼし合う力です。
流体は圧力の高いところから低いところに向かって流れます。

それ以外でも、ある程度安定した状態の流体であっても、粒子間が押し合うため常に圧力は発生しています。
人が水の中に入ると流れがなくても『圧』を感じますよね。

粘性力

圧力はわりとイメージしやすいかと思いますが、粘性力とはなんでしょうか。
粘性力を一言で言うと『周りの粒子への影響力』と言えると思います。

粘性力が高い流体としては、はちみつや溶けたチョコレートなどがイメージしやすいでしょう。
『粘性力が高い』とは、流体の流れを平均化させ(周りに強く影響し)変形しづらくさせる、ということです。

書籍から引用させていただくと

周囲の平均をとるという演算は、ラプラシアンを用いて行うことができます。

とあります。
波動方程式でも周りの平均を取って周りに力を伝搬させていくのにラプラシアンを使いますね。

ちなみにそのあたりについては格子法のナビエ・ストークス方程式の解説のところで少し触れているので参考にしてみてください。

SPH実装のフロー

さて、ざーっと概要を書いてきましたが、ここからは実際の実装について書いていきます。

SPH法では以下のフローを繰り返し行ってシミュレーションしていきます。

  1. 重み関数の係数を事前計算(*)
  2. 密度を計算
  3. 圧力項を計算
  4. 粘性項を計算
  5. 衝突判定
  6. 位置更新

(*) ... (1)は事前計算なので初回のみ計算するだけでOKです。

以下から、ひとつずつ詳しく見ていきます。

重み関数

SPH法を実装していくにあたって、この『重み関数』が重要になります。
重み関数をざっくり一言で言うと『周りの粒子から受ける力の重みを求める関数』です。

前述の通り、各粒子は周りの粒子からの影響を受けて移動します。そしてその影響の受ける度合いを決めるのがこの重み関数です。

そしてナビエ・ストークス方程式の各項を求める際にこの重み関数を利用します。
各項のところで詳細は説明しますが、各項に使われる重み関数はそれぞれ異なったものを利用します。

重み関数の係数を事前計算

重み関数を見ていく・・・前に、まずは全体を通して必要となる「物理量の離散化」を取り上げます。

物理量の離散化

普通の数学では連続した値をそのまま扱って数式を解いていきますがコンピュータではそのままでは計算できません。

コンピュータで計算ができるようにするためには『離散化』が必要となります。
Wikipediaから引用させてもらうと以下のように記述があります。

数学において、離散化 (discretization) 連続関数、モデル、変数、方程式を離散的な対応する物へ移す過程のこと。この過程は普通、それらをデジタルコンピュータ上での数値評価および実装に適したものにするために最初に行われるステップである。

離散化、というとむずかしく聞こえますが、プログラミングで考えるととても簡単な概念です。
つまりコンピュータ上で取り扱える式(=四則演算のみ)に変換することを言います。

プログラムによる例

Unityを触っている人であれば以下のような計算は一度は行ったことがあると思います。

Vector3 velocity = (currentPosition - previousPosition) / Time.deltaTime;

この処理の意味は『今の位置から前の位置を引いてそれを \Delta tで除算する=微分』ということですね。

逆にここで求めた速度を使って計算を行うと積分となります。

transform.position += velocity * Time.deltaTime;

コンピュータではこのTime.deltaTimeを使って『1フレーム』を表現し、その『飛び飛びの時間(離散した時間)』で計算を行っているわけです。このことから『離散化』のイメージが湧くかと思います。

そして上で『四則演算のみで計算を行う』と書いたように、上の式はまさに『四則演算だけ』で成り立っているのが分かるかと思います。これが離散化の考え方です。

特に今回の記事での離散化とは、とある点ごと(Unityではフレームごと)に観測し計算していくこと、と考えていいと思います。

重み関数を使った物理量の離散化

SPH法では粒子ひとつひとつが影響範囲を持っていて、他の粒子との距離が近いほど強く影響を受ける、という挙動をします。

前述の書籍から図を引用させていただくと以下のようなグラフになります。

f:id:edo_m18:20200518143342p:plain

中央が盛り上がっているグラフなのが分かるかと思います。つまり、近くの粒子の距離が近づくほど影響が強く、離れるほど影響が小さくなっていくことを図示しています。(そして一定距離以上離れている粒子からの影響は最終的にゼロになります)

SPH法における物理量をΦとすると、重み関数を用いて以下のように離散化されます。


\phi(\vec{x})=\sum_{j \in N} m_{j} \frac{\phi_{j}}{\rho_{j}} W(\overrightarrow{x_{j}}-\vec{x}, h)

ここで、 N, m, ρ, hはそれぞれ以下の意味になります。

  •  N ... 近傍粒子の集合
  •  m ... 粒子の質量
  •  ρ ... 粒子の密度
  •  W ... カーネル関数(重み関数)
  •  h ... 重み関数の影響半径

ちなみにこちらの記事(SPH法の理論 - SPH法による離散化式)から引用させていただくと以下のように記載があります。

ここで,Nは近傍パーティクルの集合,mはパーティクル質量,sph_eq_rho.gifはパーティクル密度, Wはカーネル関数である. カーネル関数は有効半径hを持ち,有効半径外では0となる(コンパクトサポート). また,その積分が1となる( \int Wd\vec{r} = 1)ように設定される

上式のようにSPHでは物理量を周囲のパーティクルの重み付き和で近似する. そして,物理量の勾配 \nabla \phiカーネル関数の導関数を用いて表す.

SPH法ではいくつかの異なる重み関数を利用して計算を行っていきます。
その際、各重み関数には決まった係数があり、また毎フレーム同じ値が利用できるので初回に求めてそれを使い回すことができます。

重み関数および係数に関してはこちらの記事に詳しく掲載されています。

www.slis.tsukuba.ac.jp

各重み関数については各項の説明のときに詳細を説明しますが、ひとつ例として上げておきます。

Poly6関数

Poly6と呼ばれる重み関数です。

 W_{poly6}(\vec{r}, h) = \alpha\begin{cases}
(h^2 - r^2)^3 & (0 \leq r \leq h) \\\
0  & (otherwise)
\end{cases} \\
\alpha \begin{cases}
\frac{4}{\pi h^8} & ... 2D \\
\frac{315}{64\pi h^9} & ... 3D
\end{cases}

前述のように、 hはひとつの粒子に注目した際に、その粒子に影響を与える他粒子を見つけるための影響半径です。
要は、その半径の円内にある粒子からのみ影響を受ける、ということを言っているわけです。

そして影響を受ける粒子の場合には(h2 - r2)3という計算で求められる値を重みとする、というわけです。
ちなみにαはさらに場合分けされていますが、これは2Dと3Dで異なる値を利用することを意味しています。

密度の計算

ここからは各フローを詳細に見ていくことにしましょう。

まずは密度を計算できるように離散化します。密度の離散化は以下の式で与えられます。


\rho(\vec{x}) = \sum_{j \in N} m_j W_{poly6}(\vec{x}_j - \vec{x}, h)

そして密度計算で利用する重み関数は以下のPoly6関数です。

Poly6関数

Poly6は上で説明したものと同じです。これは密度計算にのみ使われます。再掲すると以下で示される関数です。

 W_{poly6}(\vec{r}, h) = \alpha\begin{cases}
(h^2 - r^2)^3 & (0 \leq r \leq h) \\\
0  & (otherwise)
\end{cases} \\
\alpha \begin{cases}
\frac{4}{\pi h^8} & ... 2D \\
\frac{315}{64\pi h^9} & ... 3D
\end{cases}

今回は3D版として実装したので3D版の係数を使います。

圧力項の計算

次は圧力項を計算します。
圧力項も密度と同様、離散化して計算を行います。

なお、こちらの記事から引用させていただくと以下のように説明されています。

圧力項は圧力値の勾配で構成される. 勾配値用の離散化式を用いることで圧力項を離散化することができるが,粘性拡散項と同様にasymmetricの問題が出るので,
// 中略。引用元では離散化の式が書かれている
と圧力の平均値を用いている.

ということで、圧力の平均値を利用するようです。

圧力項の離散化の式は以下で与えられます。


f_{i}^{\text {press}}=-\frac{1}{\rho_{i}} \nabla p_{i}=-\frac{1}{\rho_{i}} \sum_{j \in N} m_{j} \frac{p_{j}-p_{i}}{2 \rho_{j}} \nabla W_{s p i k y}(\vec{x_{j}}-\vec{x}, h)

そして圧力計算でも専用の重み関数を使います。また式から分かるように、重み関数の勾配をとったものを利用します。

∇Spiky

Spikyと呼ばれる重み関数に∇を作用させたもの。

 \nabla W_{spiky}(\vec{r}, h) = \alpha\begin{cases}
(h - r)^2 \frac{\vec{r}}{r} & (0 \leq r \leq h) \\\
0  & (otherwise)
\end{cases} \\
\alpha \begin{cases}
-\frac{30}{\pi h^5} & ... 2D \\
-\frac{45}{\pi h^6} & ... 3D
\end{cases}

これも2D版と3D版があり、3D版の係数を利用します。

Tait方程式による粒子圧力の事前計算

書籍によると各粒子単体の圧力を『Tait方程式』という方程式を用いて事前に計算しておくようです。

またこちらの記事では以下のように説明されています。

非圧縮性を強制するもっともポピュラーな方法はポアソン方程式 \nabla^{2} p=\rho \frac{\nabla u}{\Delta t}を解くことであるが,そのためには巨大な線形システムを解く必要があり,非常に計算時間がかかる(MPS法ではこれをやっている).

[Becker2007]ではTait方程式を用いることである程度の非圧縮性を確保している.

要は、処理負荷が高いから簡単な式に置き換えちゃいましょう、ということ。

そしてTait方程式は以下で与えられます。


p=B\left(\left(\frac{\rho}{\rho_{0}}\right)^{\gamma}-1\right)

ここで \gamma = 7です。
またB圧力定数で、以下の式で求められます。


B = \frac{\rho_0 c_s^2}{\gamma}

ここで c_{s}^2は音速です。

これらの詳細については以下の記事を参考にしてください。

www.slis.tsukuba.ac.jp

粘性項の計算

次に粘性項を計算します。

圧力同様、こちらの記事から引用させていただくと以下のように説明されています。

粘性項は流体の速度により構成される. 速度 \vec{u}ラプラシアンをそのまま離散化すると(離散化式の \phi \vec{u}で置き換える), あるパーティクルiからjに与える力とjからiへの力が異なる(asymmetric). このままでは速度の速いパーティクルは遅いパーティクルに影響を与えるが,その逆の影響力が小さくパーティクル速度が発散してしまう.

粘性力は流体の速度差に依存すると考えて,symmetricに作用するようにしたのが以下の式である.

ということで、得られる離散化の式は以下です。


f_{i}^{v i s c}=\mu \nabla^{2} \vec{u}_{i}=\mu \sum_{j \in N} m_{j} \frac{\vec{u}_{j}-\vec{u}_{i}}{\rho_{j}} \nabla^{2} W_{v i s c}(\vec{x_{j}}-\vec{x}, h)

△Viscosity

Viscosity重み関数にラプラシアンを作用させたもの。

 \nabla W_{viscosity}(\vec{r}, h) = \alpha\begin{cases}
h - r & (0 \leq r \leq h) \\\
0  & (otherwise)
\end{cases} \\
\alpha \begin{cases}
\frac{20}{3\pi h^5} & ... 2D \\
\frac{45}{\pi h^6} & ... 3D
\end{cases}

外力項の計算

最後は外力項です。
ただ外力項はその名の通り、粒子間の影響ではなく完全に外部からの力になります。
代表的なのは重力でしょう。なにもしていなくても下に引きつけられる力がかかります。

それ以外では、例えば剛体にぶつかった場合などの反発力などもあるでしょう。

概念としては粒子の外になるので、ここでの説明は割愛します。
次回は実装編(3D版編)を書く予定なのでそちらで書きたいと思います。

実装に向けて

さて、以上まででSPH法によるナビエ・ストークス方程式の各項を求める準備が整いました。
それぞれの項をしっかり解こうとすると一筋縄では行きませんが、(正直、今の自分ではまったくできません)コンピュータでは離散化された(四則演算だけで計算可能な)シンプルな式に置き換え、それを毎フレーム足し込んでいくだけで結果を表示することができます。

コンピュータによる計算のフロー

あとはここまで解説してきた各項の離散化された式をまとめ上げて計算してやればいいことになります。
コンピュータによる計算のフローは以下になります。

  1. 各項の重み関数の係数を事前計算
  2. 密度を更新
  3. 粒子単体の圧力の計算
  4. 圧力項の更新
  5. 粘性項の更新
  6. 圧力・粘性・密度を用いて加速度を計算
  7. 求まった加速度を元に速度を計算
  8. 求まった速度から粒子の位置を更新

ナビエ・ストークス方程式での計算は1~5までで、それ以降は求まった速度から粒子の位置を更新してやるだけです。

実際のコードを抜粋すると以下のように位置を計算しています。

velocity += _TimeStep * acceleration;
position += _TimeStep * velocity;

// Update particles buffer.
_ParticlesBufferWrite[P_ID].position = position;
_ParticlesBufferWrite[P_ID].velocity = velocity;

普通の剛体の位置計算と変わらないのが分かるかと思います。
この『速度』を求めるために色々やっていたのがナビエ・ストークス方程式、というわけです。(まぁそもそも式を見たらそう書いてありますしね)

ナビエ・ストークス方程式を再掲すると以下です。


\frac{D \vec{u}}{D t}=-\frac{1}{\rho} \nabla p+\nu \nabla^2 \vec{u}+\vec{f}

まさに『速度』を求める式ですね。

あとはこれを毎フレーム実行してやれば冒頭の動画のようになる、というわけです。

まとめ

いかがでしょうか。粒子法の概要だけでも伝わってくれれば幸いです。
そして書籍を見た方はお気づきかと思いますが、ベースは書籍で書かれていることほとんどそのままです。

そこに、自分なりのメモや理解、そして別途参考にした記事などのリンクを引用させてもらいました。
(基本的に自分の理解を助けるために頭の中の整理を記事にしたものなのでそのあたりはご了承ください)

次は実装編、と行きたいところですが、実装に関しては参考にさせていただいた書籍を見たほうがいいでしょう。
あるいは、親切にも、動作するサンプルをGitHubにアップしてくださっているので、動くものを見ながら理解していくのがいいと思います。(自分もそうしました)

なので次回はちょっと視点を変えて、2D版としてアップされているサンプルを3D化した部分についてメモを添えつつ、コードをざっくり解説しようかなと思っています。

近傍探索について

さてまとめたあとになんですが、今回サンプルとして上がっているものは最適化を行わず理解のしやすさを重視しているため、そのままだとだいぶFPSがきびしい感じになります。

なぜかと言うと、粒子法では粒子は自由に動き回ることができるため、影響のある近傍の粒子を集めて計算を行う必要があります。

今回の実装では書籍に習って愚直に全部の粒子に対して距離計算を行っています。(というか、サンプルに手を加えていないとも言う)
そのため、実際に公開しているサンプルを実行してもらうと分かりますが、8,000パーティクルも出すとだいぶFPSが落ちてしまいます。

自分のPCではGTX1080を搭載していますが、それでも60FPSを下回るくらいにFPSが落ちてしまいます。

通常、こうした『全粒子に対する計算』のオーダーは[tex: O(n2)]となり、粒子数の2乗に比例して計算負荷が高まっていきます。
つまり今回の8,000パーティクル程度でもその2乗値、つまり64,000,000回の計算となってしまいます。
いくらGPUによって並列計算していると言えどもかなりの高負荷になることが分かるかと思います。

そのため、『近傍探索』を最適化し、計算回数を抑える工夫をしなければなりません。
この近傍探索についても、重み関数で紹介したサイトに解説があります。
また、今回の書籍を執筆された@kodai100さんがQiitaで近傍探索の記事を公開されています。

次はこれらを参考にして近傍探索を最適化したものを実装する予定です。(そして願わくば、自分で作っているARコンテンツに入れ込みたい)

もし興味がある方は以下の記事を参考にしてみてください。

www.slis.tsukuba.ac.jp

qiita.com

その他参考にした記事

soysoftware.sakura.ne.jp

www.slis.tsukuba.ac.jp

epa.scitec.kobe-u.ac.jp

www.epii.jp

ドメインワープをテクスチャフェッチで実装する

概要

以前、フラクタルブラウン運動とドメインワープという記事を書きました。
ここで紹介した方法は、ドメインワープの計算のためにランタイムでノイズを計算する方法でした。

ちなみにドメインワープによるアニメーションはこんな感じのものになります↓

f:id:edo_m18:20200530154711g:plain

が、これはそこそこ処理負荷が高く、これをOculus Questなどの非力なマシン上で実行すると目に見えてFPSが落ちてしまいます。
そこで今回は、このノイズ計算をテクスチャフェッチに置き換えて軽量化した話を書きます。

ちなみになんとなーくイメージが湧いているだけで、実際にそれぞれの処理が具体的に(数学的に)どういう意味があるか、という点については理解しきれていません。
あくまでOculus Questなどの非力なマシンでもそれっぽく動いた、というのをまとめたものになります。

上の動画の動作版はShadertoyにアップしてあるので実際に動くものを見ることができます。

今回実装したものをUnity上で再生したものは以下のような感じになります。
ドメインワープをテクスチャフェッチ版に置き換え、さらにそれが上に流れていくようにして炎っぽいエフェクトにしてみました。

必要なテクスチャは下で解説する「パーリンノイズで作ったテクスチャ1枚だけ」なので、結構汎用的に使えるのではないかなと思います。

解説

ドメインワープ自体の仕組みについては以前書いたこちらの記事(フラクタルブラウン運動とドメインワープ)を参照ください。

ドメインワープの主役は、ノイズ関数を、異なるスケールで足し合わせるfBM(フラクタルブラウン運動)です。
これを実行して1枚の画像を作ると以下のような画像を作る(計算する)ことができます。

f:id:edo_m18:20200530155439j:plain

冒頭のShadertoyの例では、これをランタイムで計算することにより実現していました。
ただ、ここの処理を、各フラグメントごとに実行するのは相当負荷が高いので、ここを省略したい、というのが今回の記事を書くに至った理由です。

大まかな動作原理

今回達成したいことは処理負荷の軽減です。
そのほとんどがノイズ関数の重ね合わせの計算です。

そして上の図を見てもらうと分かりますが、これをわざわざランタイムで計算しなくても良いわけですね。

ランタイムで計算するメリットは、無限に広がるノイズ空間を使える点でしょうか。
一方、テクスチャを利用する場合はある程度の精度の問題が出てくることがあるかもしれません。

が、そこらへんはスケール調整などでごまかせる範囲かと思います。

今回はこの「テクスチャフェッチ」を利用して計算負荷を下げるアプローチを取ります。

ドメインワープをもう一度

さて、ではどうテクスチャフェッチしたらいいでしょうか。
ということで、ドメインワープのやっていることを少し深堀りしてみましょう。

まずはシェーダのコードを見てみましょう。

vec2 q = vec2(0.0);
q.x = fbm(st + vec2(0.0));
q.y = fbm(st + vec2(1.0));

// These numbers(such as 1.7, 9.2, etc.) are not special meaning.
vec2 r = vec2(0.0);
r.x = fbm(st + (4.0 * q) + vec2(1.7, 9.2) + (0.15 * iTime));
r.y = fbm(st + (4.0 * q) + vec2(8.3, 2.8) + (0.12 * iTime));

float f = fbm(st + 4.0 * r);

fbm関数の呼び出しが計5回行われているのが分かります。
そしてその引数に注目すると、そのどれもが『ひとつ前のfbmの計算結果を利用して』います。

イメージ的には積分を繰り返していくような感じでしょうか。
最初が加速度で、その加速度から速度を算出し、そして最終的な位置を計算しているような。

要は、『ノイズ関数によって得られた位置を使って次の位置を決定することを繰り返す』というのが本質です。

なので、この考え方をテクスチャフェッチの位置にも適用してやることを考えてみます。

ということでまず、今回実装したUnityのシェーダの全文を載せます。(そんなに長くないので)

Shader "Unlit/FireEffect"
{
    Properties
    {
        _Color1("Color1", Color) = (1, 1, 1, 1)
        _Color2("Color2", Color) = (1, 1, 1, 1)
        _Color3("Color3", Color) = (1, 1, 1, 1)
        _Scale ("Scale", Float) = 0.01
        _NoiseScale("Noise Scale", Float) = 10.0
        _Speed("Speed", Float) = 1.0
        _PerlinNoise("Perlin Noise", 2D) = "white" {}
    }
        SubShader
    {
        Tags { "RenderType" = "Transparent" "Queue" = "Transparent" }
        LOD 100

        Pass
        {
            Blend SrcAlpha OneMinusSrcAlpha
            ZWrite Off

            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            #include "UnityCG.cginc"

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

            struct v2f
            {
                float2 uv : TEXCOORD0;
                float4 vertex : SV_POSITION;
            };

            fixed _Scale;
            fixed _NoiseScale;
            fixed _Speed;
            fixed4 _Color1;
            fixed4 _Color2;
            fixed4 _Color3;
            sampler2D _PerlinNoise;

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

            fixed4 frag(v2f i) : SV_Target
            {
                float k = _Scale;

                float2 st = i.uv * _NoiseScale;

                float3 q = tex2D(_PerlinNoise, st * k).rgr;

                float2 r;
                r.xy = tex2D(_PerlinNoise, st * k + q + (0.15 * _Time.x)).rg;

                float2 puv = (st + r) * 0.05;
                puv.y -= _Time.x * _Speed;

                float f = tex2D(_PerlinNoise, puv);

                float coef = (f * f * f + (0.6 * f * f) + (0.5 * f));

                fixed4 color = fixed4(0, 0, 0, 1);
                color = lerp(_Color1, _Color2, length(q.xy));
                color = lerp(color, _Color2, length(r.xy));
                color *= coef;

                // for alpha
                float2 uv = abs(i.uv * 2.0 - 1.0);
                float2 auv = smoothstep(0, 1, 1.0 - uv);

                float a = auv.x * auv.y;

                float3 luminance = float3(0.3, 0.59, 0.11);
                float l = dot(luminance, coef.xxx);

                l = saturate(pow(1.0 - l, 5.0));

                color.rgb = lerp(color.rgb, _Color3, l);

                color.a = a;

                return color;
            }
            ENDCG
        }
    }
}

大事な点を抜粋すると以下です。

float k = _Scale;

float2 st = i.uv * _NoiseScale;

float3 q = tex2D(_PerlinNoise, st * k).rgr;

float2 r;
r.xy = tex2D(_PerlinNoise, st * k + q + (0.15 * _Time.x)).rg;

float2 puv = (st + r) * 0.05;
puv.y -= _Time.x * _Speed;

float f = tex2D(_PerlinNoise, puv);

float coef = (f * f * f + (0.6 * f * f) + (0.5 * f));

まず最初でUVにスケールを掛けています。これによって最初にフェッチするテクスチャの位置を調整しているわけですね。
そしてその引数をさらに0.01倍しています(_Scaleのデフォルト値は0.01)。このあたりはわりと適当な数値で、自分が望んだ形になるように調整するためのパラメータです。
他でも同様の値を使うために変数にしているに過ぎません。(もちろん、他の部分で使用しているkを別の数値にしても大丈夫です)

次の部分では、フェッチした値(q)をstに加算しています。
若干値を加工していますが、これは、時間経過を利用することで徐々に変化を生むためのものです。

最終的にrが最後のフェッチする値となり、得られた値をf3+0.6f2+0.5fの式でなめらかにしたものを係数として利用しています。
ちなみにこの式をグラフ化すると以下のようななめらかに変化するグラフになります。

f:id:edo_m18:20200530172049j:plain

最後に

水のような、ちょっと神秘的な流体表現を手軽に行えるので、ちょっとしたアクセントに使うといいかもしれません。
今回の調整によってOculus Questでも、複数配置しても問題ないレベルで動作しているので利用範囲は多いかなと思っています。

Unityで流体シミュレーションを実装する 〜 実装編 〜

概要

前回は流体シミュレーションの概念編を書きました。

edom18.hateblo.jp

今回は前回説明した内容を元に実際に実装したプログラムに関して解説していきたいと思います。
ちなみに実行したデモ動画はこちら↓

また今回実装したものはGitHubにアップしてあります。

github.com



おさらい

まずは簡単におさらいから始めましょう。
流体シミュレーションはナビエ・ストークス方程式というのを解くことで実現するのでした。

ナビエ・ストークス方程式は以下の形をした方程式です。

f:id:edo_m18:20200211103647j:plain

実装ではまさにこの右辺を求めたいわけです。毎フレーム、パーティクルがどう動くのかが知りたいわけなので。
つまりプログラムで達成すべきは、各ピクセルに保存された速度を更新して次にパーティクルをどう動かすべきかを計算することです。

今回はその計算にCompute Shaderを利用しました。

Compute Shaderについては過去に2つほど記事を書いているので使い方については以下の記事を参照ください。今回は使い方に関しては割愛します。

edom18.hateblo.jp

edom18.hateblo.jp

実装のフロー

さて、まずはどう実装するかのフローを確認しておきましょう。

なお、本実装にあたり以下の「Nikkei Development Book」を大いに参考にさせていただきました。ありがとうございます。

note.com

実装についてはCompute Shaderを用いて以下のステップで各ピクセルを更新していきます。

  1. 前フレームの速度を用いて移流を計算
  2. マウスによる外力を計算
  3. (1), (2)によって計算された速度を元に発散を計算
  4. (3)の発散を元に圧力を計算
  5. (4)の圧力を元に速度を調整
  6. 求まった速度を利用してピクセルの色を更新

という流れになります。
ナビエ・ストークス方程式で言うと(1)が移流項(2)が外力項(3)と(4)が圧力項となります。

前回の概念編では書きましたが、今回は粘性項は考慮していません。毎フレーム0.99などを掛けて徐々に減速されるようになっています。

実際の実装

ここからは実際の実装を元に解説していきたいと思います。

実装自体はそれほど複雑ではなく、実際に実装したCompute Shaderのソースコードはそこまで長くありません。(ソースコードはこちら
見てもらうと分かりますが140行足らずしかありません。(C#などプロジェクト全体はGitHubを参考にしてください)

カーネルの定義

定義されているカーネルは以下。

#pragma kernel UpdateAdvection
#pragma kernel InteractionForce
#pragma kernel UpdateDivergence
#pragma kernel UpdatePressure
#pragma kernel UpdateVelocity
#pragma kernel UpdateTexture

実装フロー分のカーネルが定義されているのが分かるかと思います。

上から順に、

  • 移流の計算
  • 外力の計算
  • 発散の計算
  • 圧力の計算
  • 速度の計算
  • ピクセル位置の計算

となっています。

念の為カーネルについて補足しておくと、Compute Shaderの「核」となる関数で、C#から呼び出す(Dispatch)ことができる関数です。つまりこれらカーネルを利用してテクスチャを更新していく、というわけですね。

テクセルへのアクセス


テクセルとはテクスチャのピクセルの意味です。 Wikipediaから引用すると、

テクセル(英: texel)は、コンピュータグラフィックスで使用するテクスチャ空間の基本単位である。

と記載されています。


実装の詳細に入る前に各計算で共通している部分について先に説明しておきます。
具体的には以下の部分です。

float w = _Width;
float h = _Height;

float3 px = float3(1.0/w, 1.0/h, 0.0);
float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

テクセルの基本単位の計算

_Width_Heightはuniform変数としてC#から渡された値です。つまりはテクスチャのサイズです。

そして最後の2行が、「1テクセル分の幅(基本単位)」と「サンプリングするUV値」を計算によって求めています。

シェーダでは大体の場合において値を0.0 ~ 1.0で扱います。そのため1.0を幅および高さで割ることで1テクセル分の幅を求めることができます。

例えば128 x 128のテクスチャの場合、1.0 / 128 = 0.0078125という値が1テクセル分の幅ということになります。

UVの計算

そして最後のUVの計算ですが、(前知識として)カーネルでは計算しているスレッドのIDをuint型の値で受け取ることができます。例えば以下のように受け取ります。

void UpdateTexture(uint2 id : SV_DispatchThreadID) { ... }

上の例ではuint2型で受け取っているので、つまりはテクスチャのxyの2次元配列の添字として利用することができるわけです。(Compute Shaderのこのあたりの細かい内容については以前のスレッド編の記事をご覧ください)

なので普通にテクスチャのテクセルにアクセスするだけなら以下のように指定することができます。

// uint2なのでそのまま指定すれば該当する位置のテクセルが取得できる
float4 col = texture[id];

しかし今回の例では以下のように計算してUV値を求めています。

float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

idを幅と高さで割ることで対象とするUVの位置の左下の値が計算できます。

ここであえて左下と書いたのにはわけがあります。

どういうことかと言うと、テクスチャ空間では実はピクセルの単位ごとにアクセスするのではなくさらに細かい位置までアクセスすることが可能となっています。

データとしてはピクセル単位で保存されているわけですが、その間の値を指定した場合はその間の値を補間した値が返されます。

この補間はフィルタと呼ばれ、バイリニアフィルタなどは聞いたことがある人もいるかと思います。
このあたりの内容については以下の記事を参照ください。

karanokan.info

このことから、1ピクセルの中だけを見ても細かな位置に意味があることが分かります。
そのために「左下」と書いたわけです。

ピクセルの中心にアクセスするためにオフセットさせる

計算ではさらに続けてpx.xy * 0.5を足しています。これは「基本単位の半分だけオフセットさせる」という意味になります。

前述の通り、オフセット前の値はテクセルの「左下」の値でした。
そしてテクセル単位の半分だけオフセットさせれば無事、テクセルの中心へのアクセスとなります。

そしてそれを加味した状態でのテクセルへのアクセスは以下のようになります。

f:id:edo_m18:20200223100124p:plain

全体で利用しているSwapBufferについて

これから説明していく中で「バッファ」という言葉を使いますが、実態はRenderTexureです。
今回の実装ではこれをラップするSwapBufferというクラスを実装しました。

これは単純にふたつのRenderTextureを内部に保持し、内容を更新する際に交換(スワップ)して簡単にCompute Shaderへテクスチャを渡すことができるようにした便利クラスです。

実装は以下のようにシンプルになっています。

public class SwapBuffer
{
    private RenderTexture[] _buffers = new RenderTexture[2];

    public RenderTexture Current => _buffers[0];
    public RenderTexture Other => _buffers[1];

    private int _width = 0;
    private int _height = 0;

    public int Width => _width;
    public int Height => _height;

    public SwapBuffer(int width, int height)
    {
        _width = width;
        _height = height;

        for (int i = 0; i < _buffers.Length; i++)
        {
            _buffers[i] = new RenderTexture(width, height, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Linear);
            _buffers[i].enableRandomWrite = true;
            _buffers[i].Create();
        }
    }

    public void Swap()
    {
        RenderTexture temp = _buffers[0];
        _buffers[0] = _buffers[1];
        _buffers[1] = temp;
    }

    public void Release()
    {
        foreach (var buf in _buffers)
        {
            buf.Release();
        }
    }
}

内部にRenderTextureをふたつ保持し、スワップを簡単に行えるようにするのと、Currentとそれ以外を指定して簡単にアクセスできるようにするだけのシンプルなものです。

移流項の計算

まずは移流項から見ていきましょう。移流項はとてもシンプルです。
移流項の計算は以下になります。

[numthreads(8,8,1)]
void UpdateAdvection(uint2 id : SV_DispatchThreadID)
{
    float w = _Width;
    float h = _Height;

    float3 px = float3(1.0/w, 1.0/h, 0.0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float2 velocity = _UpdateVelocity.SampleLevel(_LinearClamp, uv, 0).xy;
    float2 result = _SourceVelocity.SampleLevel(_LinearClamp, uv - velocity * _DeltaTime, 0).xy;

    _ResultVelocity[id] = float4(result, 0.0, 1.0);
}

ここでやっていることは移流、つまり流体が流れてきた際の速度の変化を計算しています。
上記コードではふたつの速度テクスチャを利用して計算を行っています。

C#側の呼び出しも合わせて見てみましょう。

private void UpdateAdvection()
{
    _shader.SetTexture(_kernelDef.UpdateAdvectionID, _propertyDef.SourceVelocityID, _velocityBuffer.Current);
    _shader.SetTexture(_kernelDef.UpdateAdvectionID, _propertyDef.UpdateVelocityID, _velocityBuffer.Current);
    _shader.SetTexture(_kernelDef.UpdateAdvectionID, _propertyDef.ResultVelocityID, _velocityBuffer.Other);

    _shader.Dispatch(_kernelDef.UpdateAdvectionID, _velocityBuffer.Width / 8, _velocityBuffer.Height / 8, 1);

    _velocityBuffer.Swap();
}

Source VelocityUpdate Velocityは同じバッファを指定しています。
Result Velocityは結果を保持するため別のバッファを渡しています。

イメージ的には、現在の速度場の各ピクセルの速度を利用して次のフレームのための速度に更新する、ということをやっています。

現在の速度場テクスチャの値を取得しそれをvelocityとし、そのマイナス方向のテクセルの値を自身の値とする、という計算を行っています。

テクセルのフェッチ位置の計算は以下の部分ですね。

uv - velocity * _DeltaTime

現在のテクセルの位置に速度を上乗せして、その位置のテクセルをフェッチしているということです。
また取得したのが速度なので_DeltaTimeを掛けて単位を秒単位にしてから計算している点に注意が必要です。

外力の計算

次は外力の計算です。
外力は(今回は)マウスの動きを利用しています。
もちろんマウス以外から外力を受けるようにしても問題ありません。

要はこの項は流体がなにによってどう動かされているかを計算します。
なので一定時間外力を加えたあとに外力を計算しなくしても(常にゼロにしても)、流体自体に残されている速度によって徐々に拡散が起こります。(ただ速度の拡散によって徐々に動きが穏やかになっていきます。無風のときの湖面が鏡のようになっているようなイメージですね)

[numthreads(8,8,1)]
void InteractionForce(uint2 id : SV_DispatchThreadID)
{
    float w = _Width;
    float h = _Height;

    float2 px = float2(1.0 / w, 1.0 / h);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float3 vec = _SourceVelocity.SampleLevel(_LinearClamp, uv, 0).xyz;

    float dist = distance(_Cursor * px * _Scale, uv);

    if (dist < 0.005)
    {
        vec.xy += _Velocity.xy * px;
    }

    _ResultVelocity[id] = float4(vec, 1.0);
}

計算に利用している_Cursorはマウスの位置です。これとpxを乗算することで正規化されたテクスチャの位置を知ることができます。(_Scaleは単にサイズ調整のためのものです)

そしてその値と現在のUV位置との距離を測ることで、マウスの影響下にあるテクセルなのかどうかを判断しています。
その閾値として0.005をハードコードしていますが、この値を大きくすると影響を受ける範囲が大きくなります。

あとは、マウスの影響範囲だと判断された場合に、C#側から渡されたマウス速度(_Velocity)を加算することで速度を更新しています。

注意点としてはpxを掛けてテクスチャ空間での単位に変換してやる必要があります。
これは、C#側で計算された速度がスクリーン座標系のため値が数十〜数百という値になっており、これをそのまま利用するととんでもない速度になってしまうためです。

発散の計算

続いて発散の計算です。
この前の段階までで、前フレームから今フレームの移流の計算、およびマウス操作による外力が加わり、速度場が変化した状態になっています。

ただ、このままだと流体としての振る舞いにはならないため、流体らしく動くよう全体の速度を調整する必要があります。

そのための項が圧力項なのですが、全体の圧力を求めるためにこの「発散」の値が必要になります。

ということで次は発散を求めていきます。

[numthreads(8,8,1)]
void UpdateDivergence(uint2 id : SV_DispatchThreadID)
{
    float w = _Width;
    float h = _Height;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float x0 = _SourceVelocity.SampleLevel(_LinearClamp, uv - px.xz, 0).x;
    float x1 = _SourceVelocity.SampleLevel(_LinearClamp, uv + px.xz, 0).x;
    float y0 = _SourceVelocity.SampleLevel(_LinearClamp, uv - px.zy, 0).y;
    float y1 = _SourceVelocity.SampleLevel(_LinearClamp, uv + px.zy, 0).y;

    float divergence = (x1 - x0 + y1 - y0);

    _ResultDivergence[id] = float4(divergence.xx, 0.0, 1.0);
}

ここでなぜ場の発散を求めているのかというと、\(\nabla \cdot v = 0\)(発散がゼロ)を前提としているからです。
非圧縮流体の場合どの地点においても発散がゼロになるように計算をする必要があります。

そのため、次の圧縮項の計算に利用するための発散をここで求めています。

発散の求め方

発散の求め方は以下の記事を参考にしてください。

hooktail.sub.jp

上記記事から引用すると発散は以下のように定義されます。

 \begin{eqnarray} div A = \nabla \cdot A \\
 = \frac{\partial A_1}{\partial x_1} + \frac{\partial A_2}{\partial x_2} + \frac{\partial A_3}{\partial x_3} \end{eqnarray}

これは3次元による表現なので今回の場合は2次元として考えます。(つまり\(x\)と\(y\)だけについて考えます)

上記式の意味は、それぞれの要素の偏微分の値を足し合わせるという意味です。
それを踏まえてコードを見てみるとまさにそうなっているのが分かるかと思います。

float x0 = _SourceVelocity.SampleLevel(_LinearClamp, uv - px.xz, 0).x;
float x1 = _SourceVelocity.SampleLevel(_LinearClamp, uv + px.xz, 0).x;
float y0 = _SourceVelocity.SampleLevel(_LinearClamp, uv - px.zy, 0).y;
float y1 = _SourceVelocity.SampleLevel(_LinearClamp, uv + px.zy, 0).y;

float divergence = (x1 - x0 + y1 - y0);

上下左右のテクセルを取得してその差分(e.g. x1 - x0)を取っています。これはまさに微分の計算ですね。
そしてxyをそれぞれ微分したものを足し合わせているので、まさに発散の計算となっているわけです。

圧力の計算

次は圧力の計算です。圧力の計算では「ヤコビ法」と呼ばれる、反復による計算で解(の近似)を求める方法を使います。
そのためこの部分だけはC#側から複数回(今回の実装では20回)呼び出して値を更新しています。

C#側のコードを見ると単純にforループで複数回カーネルを呼び出しています。

private void UpdatePressure()
{
    for (int i = 0; i < _numCalcPressure; i++)
    {
        _shader.SetTexture(_kernelDef.UpdatePressureID, _propertyDef.SourcePressureID, _pressureBuffer.Current);
        _shader.SetTexture(_kernelDef.UpdatePressureID, _propertyDef.ResultPressureID, _pressureBuffer.Other);
        _shader.SetTexture(_kernelDef.UpdatePressureID, _propertyDef.ResultDivergenceID, _divergenceTexture);

        _shader.Dispatch(_kernelDef.UpdatePressureID, _pressureBuffer.Width / 8, _pressureBuffer.Height / 8, 1);

        _pressureBuffer.Swap();
    }
}

また圧力自体の解は「ポアソン方程式」を用いて解いています。
つまり「ポアソン方程式をヤコビ法を用いて解く」というのが今回の実装ということになります。

ポアソン方程式やヤコビ法については前回の概念編の記事を参照ください。

以下がその実装です。
このカーネル関数を複数回呼び出すことで各ピクセルにおける圧力の近似解を得ます。

[numthreads(8,8,1)]
void UpdatePressure(uint2 id : SV_DispatchThreadID)
{
    float w = _Width;
    float h = _Height;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float x0 = _SourcePressure.SampleLevel(_LinearClamp, uv - px.xz, 0).r;
    float x1 = _SourcePressure.SampleLevel(_LinearClamp, uv + px.xz, 0).r;
    float y0 = _SourcePressure.SampleLevel(_LinearClamp, uv - px.zy, 0).r;
    float y1 = _SourcePressure.SampleLevel(_LinearClamp, uv + px.zy, 0).r;

    float d = _ResultDivergence[id].r;
    float relaxed = (x0 + x1 + y0 + y1 - d) * 0.25;

    _ResultPressure[id] = float4(relaxed.xx, 0.0, 1.0);
}

ポアソン方程式は注目している点の圧力は近傍の4点の平均から求めるものです。
そのためすべての要素を足して4で割っていますね。(1/4 = 0.25

概念編の記事で参考にした記事から引用すると、

電荷が無い場合のPoissonの方程式(Laplaceの方程式)の意味するところは,「ある場所の電位は,近傍の電位の平均値に等しい」と言っているに過ぎないことに気がついただろうか.

というわけです。

ここではさらに、前段で求めたdivergenceの値を引いています。
これは非圧縮流体を過程しているため発散がゼロになるように計算する必要があります。

発散がゼロということは、その分近傍の動きからそのまま圧力を受けることと同義ですよね。そのためその圧力としての発散(divergence)を計算に加えているわけです。

速度の計算

さて、いよいよ最後の計算です。(実際にはレンダリングのためのテクスチャ更新もありますが「流体のための計算」という意味ではここが最後です)

前段までで移流項、圧力項、外力項を求めました。
その値を利用して、実際に次のフレームで流体がどうなるのかを計算していきます。

[numthreads(8,8,1)]
void UpdateVelocity(uint2 id : SV_DispatchThreadID)
{
    float w = _Width;
    float h = _Height;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float x0 = _SourcePressure.SampleLevel(_LinearClamp, uv - px.xz, 0).r;
    float x1 = _SourcePressure.SampleLevel(_LinearClamp, uv + px.xz, 0).r;
    float y0 = _SourcePressure.SampleLevel(_LinearClamp, uv - px.zy, 0).r;
    float y1 = _SourcePressure.SampleLevel(_LinearClamp, uv + px.zy, 0).r;

    float2 v = _SourceVelocity.SampleLevel(_LinearClamp, uv, 0).xy;
    float4 v2 = float4((v - (float2(x1, y1) - float2(x0, y0)) * 0.5), 0.0, 1.0);
    v2 *= _Attenuation;

    _ResultVelocity[id] = v2;
}

速度の計算は以下のように行います。

  • 現在の速度のx成分 += (左の圧力 - 右の圧力) * 0.5
  • 現在の速度のy成分 += (上の圧力 – 下の圧力) * 0.5

これをまとめると(v - (float2(x1, y1) - float2(x0, y0)) * 0.5となるわけですね。

そして今回の実装では「粘性項」を求めていません。
近似的に減衰率を速度に掛けることで調整しています。

ちなみに最後の_Attenuationが「減衰率」です。
これはC#側から渡して自由に変更できるようにしています。

ピクセルの計算

ここは流体としての計算ではなく、流体っぽい動きをするための「表現・描画」のための計算となります。
なのでここの部分を別の計算に置き換えて(例えばパーティクルの動きに転化するなど)利用することも可能です。

[numthreads(8,8,1)]
void UpdateTexture(uint2 id : SV_DispatchThreadID)
{
    float w = _Width;
    float h = _Height;

    float3 px = float3(1.0 / w, 1.0 / h, 0);
    float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

    float2 vel = _SourceVelocity.SampleLevel(_LinearClamp, uv, 0).xy;

    float4 col = _SourceTexture.SampleLevel(_LinearClamp, uv - vel * _DeltaTime, 0);

    _ResultTexture[id] = col;
}

計算自体は実は移流の計算と同じです。
つまり、速度テクスチャに保存されている速度を元に、速度分移動した位置のピクセルを持ってきて自身の色とする、ということをやっているわけですね。

移流(つまり速度)では圧力項により徐々に速度が変化していきますが、色はその速度を利用して自身の色を変化させるだけなのでこれだけです。

最後に

最後に完成動画をもう一度。

実はこれを3次元に拡張して煙っぽい動きをするパーティクルの実装として使えないかなーなんて思っています。

おそらくですが、速度テクスチャを3次元に拡張し、それらを速度の影響を与える点として利用してやればいけるのではないかなと。

もしそれっぽく動いたらそれもブログで取り上げます。

ともあれ、かねてから流体表現を実装したいと思っていたことを達成できて感無量です。
とはいえもっと応用しながら理解を深めていきたいと思います。

Unityで流体シミュレーションを実装する 〜 概念編 〜

はじめに

今回の実装をするにあたってナビエ・ストークス方程式を解いて実装していく必要があるのですが、正直、数学・物理としての正確な意味の理解や方程式のしっかりとした理解には至ってません。

式の内容やどういうことを行っているのかというのは理解して書いていますが、本当の意味での理解をしたい場合はこの記事は参考にならないかもしれません。

あくまで「プログラムで動かし、実際にインタラクションさせる」という点においての記事になります。


概要

今回は流体力学を使った流体シミュレーションを実装したのでその解説をしたいと思います。

内容が濃くなりそうなので概念編実装編を分けて書きたいと思います。今回は概念編です。

ちなみに実際に動かした動画はこちら↓

だいぶ自然な感じで動いているのが分かるかと思います。

ナビエ・ストークス方程式

さて、流体シミュレーションや流体力学で検索をかけるとこのナビエ・ストークス方程式というのが出てきます。式は以下です。


\frac{\partial \vec{u}}{\partial t} = -(\vec{u} \cdot \nabla)\vec{u} - \frac{1}{ρ}\nabla p + ν\nabla^2 \vec{u} + \vec{F}

いきなり見ると「うっ・・・」ってなりますよね。(自分はなりました。でも右辺第一項がちょっとかわいく見えるw)

が、これでもかなり簡略化された記述です。これを解き、プログラムで動かすにはどうするか、がこの記事(と次の実装編)の目的となります。

ちなみにそれぞれの項の意味は以下です。

f:id:edo_m18:20200211103647j:plain

左辺は流体の速度の時間微分を表しています。つまり、流体の速度変化を求める式が右辺ということになります。

各項ごとに説明をつけましたが、これはつまり「流体の時間微分(=速度変化)」は「移流」「圧縮」「粘性」「外力」の合成として求めることができる、というのがこの式の言いたいことです。

動画による解説

上記の各項について、そもそもなぜ流体の速度変化がこの式で求まるのでしょうか。
その答えは、つまり各項の導出は以下の動画を見ると分かりやすいと思います。図と式を用いてとても分かりやすく解説してくれています。
自分で説明できるほど理解してないから丸投げ

www.youtube.com

www.youtube.com

www.youtube.com

www.youtube.com

動画の冒頭でだいたい小さいボケをするんですが、個人的にはとても好きですw

ナビエ・ストークス方程式の理解

導出や各項の意味については上記動画をご覧ください。
ここでは、理解を深めるために方程式の各項を別の視点から見てみたいと思います。

実はナビエ・ストークス方程式流体力学における運動方程式です。
形は複雑になっていますが、古典力学同様 F = maの形になっているのです。

こちらの記事ではまた少し異なる形の式が示されています。


ρ\biggl(\frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla)\vec{u}\biggr) = -\nabla p + μ\nabla\^2 \vec{u} + ρ\vec{F}

ここで、それぞれをニュートン力学に当てはめると「 ρ」が「 m」、「 \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla)\vec{u}」が「 a」、「 -\nabla p + μ\nabla^2 \vec{u} + ρ\vec{F}」が「 f」と考えることができます。

こう見ると少しだけ簡単な気がしてきますね。(錯覚

未知数と式の数

さて、実は冒頭の式だけではこれを解くことができません。
理由は、未知数が4つに対して式が3つしかないためです。式がひとつ足りないんですね。

そこで、一般的な流体の持つ非圧縮性という性質を利用してもうひとつ式を作ります。

非圧縮性を利用する

非圧縮性とは、とある流体が場所によらず一定で圧縮されないという意味です。
イレギュラーっぽい感じがしますが、実は水や空気など、身の回りにあるものはほぼ非圧縮性があるそうです。
動画でそのあたりにも言及されていました。

さてこの非圧縮性を導入するとなぜ解けるのか。
その理由は、非圧縮性を表す式として以下の式を定義することができるためです。


\nabla \cdot \vec{v} = 0

これが意味するところは流体の速度の発散が0になることを意味しています。

「発散」がどういう意味かというと、観測のとある地点において「プラスであれば湧き出し」、「マイナスであれば吸収」となります。
しかし水などの通常の流体は突然その状態が変わることはなく、とある地点から別の地点に移動した分の質量は、別の地点にあった同量分、さらに別の場所に押し出されるため常に密度は一定に保たれます。これが「非圧縮性」ってことなんですね。

つまり、流体の各地点において発散がゼロ=流入も流出も釣り合っている状態=圧縮されていない状態=非圧縮性、ということを意味しているわけです。

この式を加えることで「未知数が4つ」「式が4つ」になるのでこれを連立して解くことができるようになります。

なお、冒頭の動画で話されていますが、現時点ではこの方程式の一般解は説かれておらず、解ければ100万ドルがもらえる懸賞問題となっているようです。(そもそも一般解が存在するかも分かっていないらしい)

実際には数値解析によって流体を表現します。

移流項の計算

移流項は、プログラムでは該当位置の速度分後ろに進んだ位置の速度を自身の速度として利用することで達成します。詳細は後編の実装編を確認ください。

ちなみに式は -(\vec{u} \cdot \nabla)\vec{u}となっていることから、今の速度の微小量分前(マイナス)の値を使うというのがなんとなく想像できるかと思います。

粘性項

実は今回の実装では粘性項は計算していません。ざっくりと毎フレーム0.99などの値を乗算して辻褄を合わせています。

こちらの記事では粘性項についても触れられているので参考になるかもしれません。

blog.livedoor.jp

外力項

外力項は任意の値です。というのも、通常流体はなにもしなければ動きません。
静かな水面を見ているとほぼ動いていないのが分かるかと思います。

つまりこの外力項は「流体をどう動かしたいか」を表します。
なのでなにもしなければ常にゼロです。

今回はマウスをドラッグした際にそれを速度として扱い、対象の点(マウスの位置の点)にマウスの速度を加えて計算しています。

圧力項の計算 - ポアソン方程式

さて最後が圧力項です。これを最後に持ってきたのは、この計算が一番要になるからです。

色々資料を漁っていると、圧力に対する解法としてポアソン方程式が出てきます。
ポアソン方程式については以下の記事が分かりやすかったです。

teamcoil.sp.u-tokai.ac.jp

記事から引用させてもらうと、

電荷が無い場合のPoissonの方程式(Laplaceの方程式)の意味するところは,「ある場所の電位は,近傍の電位の平均値に等しい」と言っているに過ぎないことに気がついただろうか.

と書かれています。
そして記事に書かれている式を引用させてもらうと以下の式が最後に紹介されている式になります。


\phi (i, j) = \frac{1}{4} \biggl\{ \frac{\delta^2 \rho(i, j)}{\epsilon_0} + \phi (i + 1, j) + \phi (i - 1, j) + \phi (i, j + 1) + \phi (i, j - 1) \biggr\}

上記記事が解説している式をプログラム的に解くと、以下のように、計算したいピクセル位置の近傍、上下左右をサンプリングして計算を行うことが出来ます。
(後編で紹介する実際に実装したコードから抜粋)

float3 px = float3(1.0 / w, 1.0 / h, 0);
float2 uv = float2(id.x / w, id.y / h) + px.xy * 0.5;

float x0 = _SourcePressure.SampleLevel(_LinearClamp, uv - px.xz, 0).r;
float x1 = _SourcePressure.SampleLevel(_LinearClamp, uv + px.xz, 0).r;
float y0 = _SourcePressure.SampleLevel(_LinearClamp, uv - px.zy, 0).r;
float y1 = _SourcePressure.SampleLevel(_LinearClamp, uv + px.zy, 0).r;

float d = _ResultDivergence[id].r;
float relaxed = (x0 + x1 + y0 + y1 - d) * 0.25;

x0 + x1 + y0 + y1の部分が上下左右の近傍の値を足しているところですね。そして上記式の最初の項がdに該当するのだと思います。
(すみません、このあたりは理解が浅いです)
そして最後に0.25(=  \frac{1}{4})を掛けています。

ポアソン方程式をヤコビ法を使って解く

実際のプログラムではこれを複数回実行して近似値を求めています。(これをヤコビ法と言います)

ヤコビ法についてはこちらのPDFが分かりやすかったです。

例を引用させてもらうと以下のようになります。

1.1 ヤコビ法
例 1. 方程式


0 = 3x + y + z  \\\
4 = x + 3y + z  \\\
6 = x + y + 3z

が与えられたとき, 各方程式を各々の変数について解き


x = \frac{1}{3}(0 − y − z)  \\\
y = \frac{1}{3}(4 − x − z) \\\
z = \frac{1}{3}(6 − x − y)

 x^{(0)}, y^{(0)}, z^{(0)} を適当に決めて(例えば   (x^{(0)}, y^{(0)}, z^{(0)}) := (0, 0, 0)  x^{(1)}, y^{(1)}, z^{(1)}


x^{(1)} =\frac{1}{3}(0 − y^{(0)} − z^{(0)}) = 0 \\\
y^{(1)} = \frac{1}{3}(4 − x^{(0)} − z^{(0)}) = \frac{4}{3} \\\
z^{(1)} = \frac{1}{3}(6 − x^{(0)} − y^{(0)}) = 2

とし, 以下同様に


x^{(k)} = \frac{1}{3}(0 − y^{(k − 1)} − z^{(k − 1)}) \\\
y^{(k)} = \frac{1}{3}(4 − x^{(k − 1)} − z^{(k − 1)}) \\\
z^{(k)} = \frac{1}{3}(6 − x^{(k − 1)} − y^{(k − 1)})

で順次 x^{(k)}, y^{(k)}, z^{(k)} を計算していけば解に近付く事がある. 実際, 今の場合を計算すれば

[ k] x y z
[ 0] 0.000 0.000 0.000
[ 1] 0.000 1.333 2.000
[ 2] -1.111 0.667 1.556
[ 3] -0.741 1.185 2.148
[ 4] -1.111 0.864 1.852
[ 5] -0.905 1.086 2.082
[ 6] -1.056 0.941 1.940
[ 7] -0.960 1.039 2.038
[ 8] -1.026 0.974 1.974
[ 9] -0.983 1.017 2.017
[10] -1.012 0.988 1.988
[11] -0.992 1.008 2.008
[12] -1.005 0.995 1.995
[13] -0.997 1.003 2.003
[14] -1.002 0.998 1.998
[15] -0.998 1.002 2.002

となり, 解 (x, y, z) = (−1, 1, 2) に近付いている事が見える.

上記のように、繰り返せば繰り返すだけ解に近づいて行っているのが分かるかと思います。
これを用いて、本来は複雑な方程式を解いていくことになります。


ざーっと書いてきましたが、プログラムで扱うための前知識として各項について解説しました。

次回はこれを実際にプログラムに落とし込んで、今回実装した内容を解説しようと思います。

気づいた点や混乱した点のメモ

解説自体は以上なんですが、この実装をする前に色々と調べたものをだーっとメモしておこうと思います。(完全なる備忘録)
なのでここから先は完全に未来の自分宛てへのメモなので興味がある方だけ読んでください。


数式について調べていくうちに、いくつか気づいた点や理解する上で混乱した部分などをメモとして残しておきたいと思います。

移流項の内積は「発散ではない」

移流項は以下のように表されています。


-(\vec{u} \cdot \nabla)\vec{u}

 \nabla内積のように作用させたものを「発散」と呼びますが、上記式の \vec{u} \cdot \nablaは発散ではなく、単なる微分演算子との掛け算なので混同しないよう注意が必要です。

ちなみに「発散」は以下の式。


\nabla \cdot \vec{u}

掛ける順番が違いますね。

こちらのPDFで言及されていて気づきました。


余談

後編の実装編を見てもらうと分かりますが、これは「速度分少しだけ移動したところの値を取ってくる」くらいの意味だと思います。(微分演算子ってことは「微小量」を表すので、「速度方向に少しだけ進んだ地点」を意味しているってことだと理解しています)

なお、こちらの「ナビエ・ストークス方程式を導出する」ではひとつずつ丁寧に導出してくれていて、その過程でどうしてこういう表記になったのかがイメージできました。

そこでの解説から式を引用させてもらうと以下のようになります。


\frac{Df}{Dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}\frac{\partial x}{\partial t} + \frac{\partial f}{\partial y}\frac{\partial y}{\partial t} + \frac{\partial f}{\partial z}\frac{\partial z}{\partial t}

なおここで、


\frac{\partial x}{\partial t} = u,  \frac{\partial y}{\partial t} = v, \frac{\partial z}{\partial t} = w

であることに注意すると、


\frac{Df}{Dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}u + \frac{\partial f}{\partial y}v + \frac{\partial f}{\partial z}w \tag{2-1}

となります。この物理量fが速度uだったらどうなりますか? Du/Dtは加速度ですよね。


f = \vec{u} = u\vec{i} + v\vec{j} + w\vec{k}

そうすると(2-1)式はベクトルになって、それぞれの方向の式が出来上がります。


\begin{eqnarray}
\frac{D u}{D t}\vec{i} &=& a_x \vec{i} = \biggl(\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z}\biggr) \vec{i} \\\
\frac{D v}{D t}\vec{j} &=& a_y \vec{j} = \biggl(\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z}\biggr) \vec{j} \\\
\frac{D w}{D t}\vec{k} &=& a_z \vec{k} = \biggl(\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z}\biggr) \vec{k} \tag{2-2}
\end{eqnarray}

そして加速度ベクトルは


\vec{a} = a_x \vec{i} + a_y \vec{j} + a_z \vec{k} \tag{2-3}

(2-2)式を一つの式にまとめようとするなら


\vec{a} = \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla) \vec{u} \tag{2-4}

となります。

ナビエ・ストークスの式に出てくる形ですね。

そしてこれはこうも考えることができます。
上記式は「ベクトルの勾配を取ったものとの発散を取る」ことになり、ランクは変わりません。(ランクについてはこちらの記事が分かりやすいです→ ベクトル解析の基本 発散 勾配 回転

つまり、 (u \cdot \nabla)uを計算してもベクトルの性質は残ったままになっている、というわけですね。
(ナブラ( \nabla)についてはこの記事の後半に使い方などをまとめています)

ナビエ・ストークス方程式を、できるだけ平易に解説してくれている記事があるので紹介しておきます。

shimaphoto03.com


連続の式

上記解説で出てくる「連続の式」については以下の記事を参照。
要は、流体は連続しており、突然増えたり減ったりしないことを前提にしている、ということだと思います。

shimaphoto03.com

連続の式は流体力学における基本の式で「流体が突然発生したり、消滅したりしない」ということをコンセプトとしています。流体が発生したり、消滅したりしないのを示すには、流体を非常に非常に小さい要素に分割して、その要素内に出入りする流体の量に過不足がないという方程式を立てればよいのです。

微分

微分とは、ある関数の接線を1次関数で近似することを言います。
以下の記事に詳しく載っているので詳細はそちらをご覧ください。

mathwords.net

大事な点だけを引用させてもらうと以下のようになります。

接線の傾き(つまり、関数を1次関数で近似したときの傾き)が微分です。すなわち、  f(x_0+h)=f(x_0)+ah+o(h) を満たす  a f(x) の x=x_0 における微分係数です。ただし、 o(h) は微小量 (\lim_{h \to 0} \frac{o(h)}{h}=0 を満たす関数)です。

さらに、

二変数関数の場合に確認してみます。全微分可能な場合、
 f(x_0 + h, y_0 + k) = f(x_0, y_0) + ah + bk + o(\sqrt{h^2 + k^2})
という式が成立しますが、このとき  x についての偏微分
 \lim_{h \to 0}\frac{f(x_0 + h, y_0) − f(x_0, y_0)}{h} = a
となります(極限値が存在します)。同様に、 y についての偏微分 b となります。

つまり、全微分偏微分の間には、
 a = \frac{\partial f}{\partial x}、b = \frac{\partial f}{\partial y}
という関係式が成立します。言い換えると、全微分(一次近似式)に登場する一次の各係数は偏微分になります。よって、全微分可能な場合、
 df = \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y}dy
などと書くことができます。

ナブラ( \nablaの使い方。

ナビエ・ストークス方程式にはこのナブラ( \nabla)がたくさん登場します。
特に、ナブラを2乗したものを「ラプラシアン」と呼んで区別することもあります。

使い方・意味については以下の通り。

使い方1. 勾配

関数 f(x, y, z)に対して各変数での微分を並べたベクトルを \nabla fと書く。具体的には、


\nabla f = \biggl( \frac{\partial f}{\partial  x}, \frac{\partial f}{\partial  y}, \frac{\partial f}{\partial  z} \biggr)

使い方2. 発散

 \nablaをベクトルに作用させ、内積のように足し合わせたもの。


\nabla \cdot V = \frac{\partial V_x}{\partial  x} + \frac{\partial V_y}{\partial  y} + \frac{\partial V_z}{\partial  z}

使い方3. 回転

 \nablaをベクトルに作用させ、外積のように計算したもの。


\nabla \times V = \biggl( \frac{\partial V_z}{\partial  y} - \frac{\partial V_y}{\partial  z}, \frac{\partial V_x}{\partial  z} - \frac{\partial V_z}{\partial  x}, \frac{\partial V_y}{\partial  x} -  \frac{\partial V_x}{\partial  y} \biggr)

使い方4. ラプラシアン

関数 f(x, y, z)に対して、各変数での2階微分の和を \nabla^2 fと書く。


\nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2}

ちなみに、こちらの記事で書かれているように勾配の発散と考えたほうがしっくりくるかも。


\nabla^2 f = \nabla \cdot (\nabla f)

演算子の位置について

演算子の作用:演算子と関数の順番を変えると意味が変わる(参考

  • 関数が演算子の右側にあるときは作用を受ける
  • 関数が演算子の左側にあるときは作用を受けない

以下の式を見ると分かりやすい。


\frac{\partial}{\partial x} fg = \frac{\partial fg}{\partial x}

f \frac{\partial}{\partial x} g = f \frac{\partial g}{\partial x}

\biggl(\frac{\partial}{\partial x}f \biggr) g = \frac{\partial f}{\partial x} g

参考にした記事