e.blog

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

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)を取っています。これはまさに微分の計算ですね。
そして\(x\)と\(y\)をそれぞれ微分したものを足し合わせているので、まさに発散の計算となっているわけです。

圧力の計算

次は圧力の計算です。圧力の計算では「ヤコビ法」と呼ばれる、反復による計算で解(の近似)を求める方法を使います。
そのためこの部分だけは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 $$

参考にした記事

シャボン玉風シェーダを実装してみたのでメモ

概要

今回の実装はこちらの論文(汎用的な構造色レンダリング手法の開発)を参考に実装したものになります。
ただ理解が足りていない部分や勘違いなどあるかもしれないので、もし間違いなどあれば指摘いただけると嬉しいです。

また今回の実装にあたってこちらのブログも参考にさせていただきました。

light11.hatenadiary.com

実装して実行した結果がこちら↓

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

github.com

論文概要

まずはざっくり概要を。
最近の3DCGのレンダリングでは様々な物質の光の反射をモデル化して、それを表す近似した関数を用いてレンダリングを行います。
レンダリング方程式、なんかでググると色々出てきます。いわゆる「物理ベースレンダリング」ってやつですね。

この記事が参考になるかも?↓

qiita.com

そして物理ベースレンダリングを調べていると「BRDF」という言葉が出てきます。
(ちなみに自分も以前、調べて記事を書いたのでよかったらこちらもどうぞ↓)

qiita.com

また、こちらのシリーズ記事もとても良いです↓

qiita.com

なぜ物理ベースレンダリング

なぜこれだけ物理ベースレンダリングの記事を紹介したかと言うと、今回の論文の中で「BRDF」という言葉が出てきます。

BRDFは「Bidirectional Reflectance Distribution Function」の頭文字を取ったものです。

BRDFについては以下の記事が詳細に書かれています。(数式めっちゃ出てくるのでうってなるかもですが・・)

rayspace.xyz

このBRDF、前述した「物質の光の反射をモデル化」したものです。上記記事から引用させていただくと、

BRDFは物体表面の位置、入射方向、反射方向を変数とする6次元関数で、反射に関する物体表面の特性を表す。

というもの。つまりは「関数」ってことです。

構造色をテクスチャフェッチで解決する

さて、BRDFは光の反射をモデル化したものと書きましたが、現在のレンダリング方程式ではすべての物質の光の反射を扱えるほど精巧には出来ていません。ある一定の条件でのみ成り立つものだったり、特定の表現のための関数などもあります。それを匠に使いこなしながら、今日のフォトリアリスティックな映像は作られているというわけですね。

そして今回の論文ではそのBRDF的な関数を提案するものではありません

複雑な関数から導き出すのではなくテクスチャに情報を格納し、そこから構造色(詳細は後述)を決定しよう、というものです。

構造色とは

論文の説明から引用させてもらうと、

構造色は微細構造がもたらす光学現象による発色であるため、視点位置や照明環境に応じて色が変化するという特徴がある。

シャボン玉やCDの記録面を見ると分かると思いますが、見る角度や照明の状態に応じて色が変化しますよね。ああいう、視点や光源によって変化する色を「構造色」と言うようです。

構造色の仕組み

構造色は前述のように、視点と光源位置によって見え方が変わるものです。ではなぜ、視点や光源位置に応じて色が変化するのでしょうか。

論文から画像を引用させてもらうと以下のような構造になっているものが「構造色」を表します。

f:id:edo_m18:20200123135709p:plain

膜やスリットなどに光が入り、入射した光と反射した光の「光路」に差が生じると、別の反射光と干渉を起こして色が変化する、というのが理由のようです。(光は波の性質も持っているのでこうした干渉が起こる)

光の反射

反射光の話が出ました。これはまさに、前述のBRDFの話と同じですね。
3DCGでは光の反射をどう扱うかで色が決まります。

そして今、入射光と反射光について以下のように定義します。

f:id:edo_m18:20200123140033p:plain

上図は論文から引用したものです。
これは、注目している点の光の入射・反射を\(\theta\)と\(\phi\)によって表した図です。
入射光の\(\theta\)と\(\phi\)、そして反射光の\(\theta\)と\(\phi\)の4つのパラメータがあります。

この4パラメータが大事になります。

構造色を決定する

さて、大まかに概観したところで実際に構造色を決定する部分を見ていきましょう。
前述のように、入射・反射の角度の4パラメータが大事なポイントです。

また論文から画像を引用させていただくと構造色を決定するフローは以下のようになります。

f:id:edo_m18:20200123170156p:plain

視点位置\(E\)、光源位置\(L\)、それによって決定する\((\theta_i, \phi_i)\), \((\theta_o, \phi_o)\)の4パラメータ、そしてUV値がどのように使われるかを示した図です。

続いて以下の図を見てください。

f:id:edo_m18:20200123172222p:plain

パラメータをどう扱うかを図示したものです。

詳細を説明する前にさらに次の図を見てください。

f:id:edo_m18:20200123172347p:plain

こちらは実際にレンダリングに使用するテクスチャです。
左上がシャボン玉の厚みテクスチャで、UV値を利用して単純にアクセスされるテクスチャです。

次に右上が光路差を格納したテクスチャで、視点・光源の角度によってアクセスされるテクスチャです。

そして下の最後のものが構造色テクスチャです。
最終的に色としてマッピングされるのはこの構造色テクスチャです。

つまり、下部の構造色テクスチャのどの位置をフェッチしたらいいか、を光路差情報を元に決定するのが実装概要となります。

構造色テクスチャのフェッチ位置を計算

さてではどのようにしてフェッチ位置を決定するのか。それが、ふたつ前に載せた画像です。再掲すると以下の画像です。

f:id:edo_m18:20200123172222p:plain

左の図は通常のUV値を利用して「厚みテクスチャ」にアクセスする様子を描いています。
場所依存の「厚み」を表現しているわけですね。ある意味で法線マップと同様のことをやっています。なのでここは細かい説明は不要でしょう。

右側のテクスチャが今回の肝である「光路差テクスチャ」です。
前述した3枚のテクスチャのうち、「視点・光源依存テクスチャ」と書かれているテクスチャがそれです。

そしてこのテクスチャから値を取得するわけですが、「どこの値を取り出すのか」を決定するのが、前段で大事だと話した4つのパラメータ、つまり\((\theta_i, \phi_i)\), \((\theta_o, \phi_o)\)です。

論文では以下のように説明されています。

視点・光源依存テクスチャは、視点および光源位置による光路差の変化を表現するテクスチャである。テクスチャは図3.3の座標軸をもち、横軸に法線成分、縦軸に接線成分が対応している。ただしTは接線ベクトルであり、テクスチャ座標のUの方向を持つ単位ベクトルである。左半分が光源位置による光路差であり、光源\(L\)により\(L \cdot N\)および\(L \cdot T\)が決まる。同様に右半分が視点位置による光路差であり、視点\(E\)により\(E \cdot N\)および\(E \cdot T\)が決まる。よってテクスチャの中心がオブジェクトの面に垂直な位置、つまり法線上に視点・光源がある場合に対応する。

特に、

テクスチャの中心がオブジェクトの面に垂直な位置、つまり法線上に視点・光源がある場合に対応する。

を見ると、テクスチャの中心が視点・光源がともに法線上にある場合に対応することになります。

計算してみましょう。法線上ということは(\(L \cdot N = 1\)、\(L \cdot T = 0\))同様に(\(E \cdot N = 1\)、\(E \cdot T = 0\))となります。
この値が、視点・光源依存テクスチャのどこにあるかを見てみると確かにテクスチャの中心にマッピングされますね。

さて、ではこの事実をどう利用するのでしょうか。考え方はこうです。
視点・光源依存テクスチャは0 ~ 1の正規化された値が格納されています。値はただのfloatですね。なのでスカラーです。しかしテクスチャから色をフェッチするためにはUV、つまりfloat2のベクトルとしての値が必要です。

つまり、入射光の\((\theta_i, \phi_i)\)のふたつのパラメータによってひとつのスカラーが得られます。同様に反射光の\((\theta_o, \phi_o)\)のふたつのパラメータによってひとつのスカラーが得られます。

このふたつのスカラーを利用すればUVとしてのベクトルを得ることができますね。

論文の別の画像を引用させてもらうと以下のように説明されています。

f:id:edo_m18:20200124091431p:plain

つまり入射光の角度から算出されたスカラーをU軸に、反射光の角度から算出されたスカラーをV軸に割り当てることで構造色が決定できる、というわけですね。

そしてそれを実際に実行した結果がこちら↓

f:id:edo_m18:20200124084234g:plain

値を計算しているところのコードを抜粋すると以下になります。

float d = /* Calculate thickness */

float u, v;

// Calculate U.
{
    float ln = dot(i.lightDir, i.normal);
    ln = (ln * 0.5) * d;

    float lt = dot(i.lightDir, i.tangent);
    lt = ((lt + 1.0) * 0.5) * d;

    u = tex2D(_LETex, float2(ln, lt)).x;
}

// Calculate V.
{
    float en = dot(i.viewDir, i.normal);
    en = ((1.0 - en) * 0.5 + 0.5) * d;

    float et = dot(i.viewDir, i.tangent);
    et = ((et + 1.0) * 0.5) * d;

    v = tex2D(_LETex, float2(en, et)).x;
}

Uの値は入射光(i.lightDir)によって計算され、Vの値は反射光(i.viewDir)によって計算されていますね。
ここで求まったUVの値を使って構造色テクスチャから色をフェッチしてあげれば晴れて、視点・光源に依存した色、つまり構造色を決定することができる、というわけです。

最後に

構造色の計算については以上です。
今回はそれ以外にもシンプルなPhong反射とリムライティングを追加しています。なので全部が正確なレンダリング、というわけではないですがそれなりに説得力のある絵になったかなと思います。

コード全体は長くないので全文を掲載しておきます。実際に動作するものを見たい場合はGitHubからダウンロードしてご覧ください。

Shader "Unlit/Bubble"
{
    Properties
    {
        [PowerSlider(0.1)] _F0("F0", Range(0, 1)) = 0.02
        _RimLightIntensity ("RimLight Intensity", Float) = 1.0
        _Color ("Color", Color) = (1, 1, 1, 1)
        _MainTex ("Texture", 2D) = "white" {}
        _DTex ("D Texture", 2D) = "gray" {}
        _LETex ("LE Texture", 2D) = "gray" {}
        _CubeMap ("Cube Map", Cube) = "white" {}
    }

    SubShader
    {
        Tags
        {
            "RenderType"="Transparent"
            "Queue"="Transparent"
        }
        LOD 100

        Pass
        {
            Blend SrcAlpha OneMinusSrcAlpha

            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            #include "UnityCG.cginc"
            #include "./PerlinNoise.cginc"

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

            struct v2f
            {
                float2 uv : TEXCOORD0;
                float4 vertex : SV_POSITION;
                float3 normal : TEXCOORD1;
                float3 tangent : TEXCOORD2;
                float3 viewDir : TEXCOORD3;
                float3 lightDir : TEXCOORD4;
                half fresnel : TEXCOORD5;
                half3 reflDir : TEXCOORD6;
            };

            sampler2D _MainTex;
            sampler2D _DTex;
            sampler2D _LETex;

            UNITY_DECLARE_TEXCUBE(_CubeMap);

            float _F0;
            float _RimLightIntensity;
            float4 _MainTex_ST;
            fixed4 _Color;

            v2f vert (appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.uv = TRANSFORM_TEX(v.uv, _MainTex);
                o.normal = v.normal;
                o.tangent = v.tangent;
                o.viewDir  = normalize(ObjSpaceViewDir(v.vertex));
                o.lightDir = normalize(ObjSpaceLightDir(v.vertex));
                o.fresnel = _F0 + (1.0h - _F0) * pow(1.0h - dot(o.viewDir, v.normal.xyz), 5.0);
                o.reflDir = mul(unity_ObjectToWorld, reflect(-o.viewDir, v.normal.xyz));
                return o;
            }

            fixed4 frag(v2f i) : SV_Target
            {
                i.uv = pow((i.uv * 2.0) - 1.0, 2.0);
                //float d = tex2D(_DTex, i.uv + _Time.xy * 0.1);
                float d = perlinNoise((i.uv + _Time.xy * 0.1) * 3.0);

                float u, v;

                // Calculate U.
                {
                    float ln = dot(i.lightDir, i.normal);
                    ln = (ln * 0.5) * d;

                    float lt = dot(i.lightDir, i.tangent);
                    lt = ((lt + 1.0) * 0.5) * d;

                    u = tex2D(_LETex, float2(ln, lt)).x;
                }

                // Calculate V.
                {
                    float en = dot(i.viewDir, i.normal);
                    en = ((1.0 - en) * 0.5 + 0.5) * d;

                    float et = dot(i.viewDir, i.tangent);
                    et = ((et + 1.0) * 0.5) * d;

                    v = tex2D(_LETex, float2(en, et)).x;
                }

                float2 uv = float2(u, v);
                float4 col = tex2D(_MainTex, uv);

                float NdotL = dot(i.normal, i.lightDir);
                float3 localRefDir = -i.lightDir + (2.0 * i.normal * NdotL);
                float spec = pow(max(0, dot(i.viewDir, localRefDir)), 10.0);

                float rimlight = 1.0 - dot(i.normal, i.viewDir);

                fixed4 cubemap = UNITY_SAMPLE_TEXCUBE(_CubeMap, i.reflDir);
                cubemap.a = i.fresnel;

                col *= cubemap;
                col += rimlight * _RimLightIntensity;
                col += spec;

                return col;
            }
            ENDCG
        }
    }
}

Unityで簡易TexturePackerを実装してみた

はじめに

こちらの記事を参考にUnityで簡易的なTexturePackerを実装してみたのでそのまとめです。

blackpawn.com

ちなみにランダムに生成したものを配置したのがこれ↓

実装しようと思った理由は今実装中の機能に必要になったためです。
今、下の動画のようなスタンプ機能を作っているのですが、プレビュー用には1枚のテクスチャしか渡せないため複数設定されたテクスチャを自前でパッキングして利用したいと思ったのが理由です。

※ ただ、実際のスタンプ実装はまったく別の方法に切り替えたのでテクスチャパックは必要なくなりましたが・・w

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

github.com

概要

今回参考にした記事はライトマップをパックするものを解説しているようです。
ライトマップは特に、ひとつのテクスチャに大量のライト情報をパックして利用するためこうした方法が利用される最適なケースでしょう。

上でも書いたように、今回の実装理由はライトマップと意図は同じで複数のテクスチャをひとつにまとめてシェーダに送りたい、というのが発端です。

実装のフロー

まずはざっくり全体像を概観してみます。
参考にした記事にはコンセプトの説明と図解、そして疑似コードが掲載されています。
今回の実装はその疑似コードを元に実装したものになります。


実装は二分木構造になっていて、最初はルートノードのみが存在します。

そしてテクスチャをパックするたびに、そのテクスチャがピッタリ収まるノードに分割していき、そこへ画像を設定していく、というイメージです。

もう少し具体的なフローにすると、

  1. とあるノードに対して、パックしようとしているテクスチャが収まるかチェックする(だいたいの場合はまったく収まらないか完全に収まるかの2択)
  2. 入る場合かつ完全フィットしていない場合はその領域を、対象画像が収まるサイズとそれ以外の領域に分ける(一度に縦横を分割するのではなく、そのテクスチャの長辺方向にだけ分割するのがポイント
  3. 上記で分割したノードに対して再びフィットするかをチェック(前段で長辺側を区切りに分割しているため、必ず短辺はノードにフィットする)
  4. (2)と同様のチェックを行うと、今度は必然的に短辺側に区切ることになる。(最終的にはそのテクスチャがぴったり収まるサイズに分割される)
  5. 対象画像が完全にフィットするようになったノードに対してImageIDを設定しリーフノード扱いにする
  6. 次のテクスチャを追加する
  7. ルートノードから子ノードを辿り空いているノードImageIDがないノード)を探す
  8. 以後、(1)〜(7)をテクスチャの数だけ繰り返す

図解すると以下のようになります。

f:id:edo_m18:20191222112136j:plain

見てもらうと分かるように、ひとつのノードに対して分割は「2回」発生します。
最初自分が勘違いしていた部分なのですが、テクスチャのサイズに同時に縦横を分割するのではなく、あくまで一度のチェックでは「長辺方向」にだけ分割します。

該当部分のコードを抜粋すると以下のようになります。

float dw = Rectangle.width - image.Width;
float dh = Rectangle.height - image.Height;

if (dw > dh)
{
    Child[0].Rectangle = new Rect(Rectangle.x, Rectangle.y, image.Width, Rectangle.height);
    Child[1].Rectangle = new Rect(Rectangle.x + image.Width + 1, Rectangle.y, Rectangle.width - image.Width - 1, Rectangle.height);
}
else
{
    Child[0].Rectangle = new Rect(Rectangle.x, Rectangle.y, Rectangle.width, image.Height);
    Child[1].Rectangle = new Rect(Rectangle.x, Rectangle.y + image.Height + 1, Rectangle.width, Rectangle.height - image.Height - 1);
}

ここでRectangleは確認しているノードの矩形を表していて、imageがパックしようとしている画像を表しています。

そしてそれぞれの幅と高さを引き、どちらが大きいかで分岐しています。

分岐内のコードはほぼ同じで、長辺方向に分割しています。

するとひとつのノードがふたつのノードに分割されます。
そして、プログラムのロジック的には再帰処理的に、分割された片方のノードに対してテクスチャの挿入処理を繰り返します。

すると、自動的に該当テクスチャがぴったりと収まるサイズに分割されたノードが出来上がります。

コードで言うと以下のところが、2回目のチェック時にはdw(あるいはdh)が0となり、前段で分割した方向とは異なる方向が分割方向として採用され、結果、対象テクスチャがぴったり収まるノードができあがる、というわけです。

float dw = Rectangle.width - image.Width;
float dh = Rectangle.height - image.Height;

if (dw > dh)
{
    // Check witch edge is short.
}

コード全容

実装は上記の図の通りにノードを分割していき、パックする位置が決定したらそのノードに対してImageIDを付与しリーフノードとする、という処理をテクスチャ分だけ繰り返すのみです。

あとはコードを見たほうが理解が早いと思うのでコード全体を載せておきます。

public class Node
{
    public Node[] Child = new Node[2];
    public Rect Rectangle;
    private int _imageID = -1;

    private bool _isLeafNode = true;

    private bool CheckFitInRect(IPackImage image)
    {
        bool isInRect = (image.Width <= Rectangle.width) &&
                        (image.Height <= Rectangle.height);
        return isInRect;
    }

    private bool CheckFitPerfectly(IPackImage image)
    {
        bool isSameBoth = (image.Width == Rectangle.width) &&
                            (image.Height == Rectangle.height);
        return isSameBoth;
    }

    public Node Insert(IPackImage image)
    {
        if (!_isLeafNode)
        {
            Node newNode = Child[0].Insert(image);
            if (newNode != null)
            {
                return newNode;
            }

            return Child[1].Insert(image);
        }
        else
        {
            if (_imageID != -1)
            {
                return null;
            }

            if (!CheckFitInRect(image))
            {
                return null;
            }

            if (CheckFitPerfectly(image))
            {
                return this;
            }

            _isLeafNode = false;

            Child[0] = new Node();
            Child[1] = new Node();

            float dw = Rectangle.width - image.Width;
            float dh = Rectangle.height - image.Height;

            if (dw > dh)
            {
                Child[0].Rectangle = new Rect(Rectangle.x, Rectangle.y, image.Width, Rectangle.height);
                Child[1].Rectangle = new Rect(Rectangle.x + image.Width + 1, Rectangle.y, Rectangle.width - image.Width - 1, Rectangle.height);
            }
            else
            {
                Child[0].Rectangle = new Rect(Rectangle.x, Rectangle.y, Rectangle.width, image.Height);
                Child[1].Rectangle = new Rect(Rectangle.x, Rectangle.y + image.Height + 1, Rectangle.width, Rectangle.height - image.Height - 1);
            }

            return Child[0].Insert(image);
        }
    }

    public void SetImageID(int imageID)
    {
        _imageID = imageID;
        _isLeafNode = true;
    }

    public int GetImageID()
    {
        return _imageID;
    }

    public Node Find(int imageID)
    {
        if (imageID == _imageID)
        {
            return this;
        }

        if (_isLeafNode)
        {
            return null;
        }

        Node child = Child[0].Find(imageID);
        if (child != null)
        {
            return child;
        }

        child = Child[1].Find(imageID);

        return child;
    }
}

シェーダによる書き込みと読み込み

今回のPackerはシェーダによって書き込みを行っています。
また書き込んだだけでは使えないので、実際に利用する際の読み込み用のシェーダも必要となります。

次はそれらシェーダについて解説します。

シェーダによる書き込み

CPU側から適切にパラメータを設定したのち、シェーダで書き込みます。
C#側の処理は以下のようになっています。

private void Pack(IPackImage image, Rect rect)
{
    Vector4 scaleAndOffset = GetScaleAndOffset(rect);

    _material.SetVector("_ScaleAndOffset", scaleAndOffset);
    _material.SetTexture("_PackTex", image.Texture);

    Graphics.Blit(_current, _next, _material);

    SwapBuffer();
}

_ScaleAndOffsetがパック先テクスチャの「どの位置にどれくらいのサイズで」書き込むかのパラメータで_PackTexが書き込むテクスチャです。

さて、これを利用しているシェーダはどうなっているかというと、

fixed4 frag (v2f i) : SV_Target
{
    float2 puv = i.uv;
    puv -= _ScaleAndOffset.zw;
    puv *= _ScaleAndOffset.xy;

    fixed4 col = tex2D(_MainTex, i.uv);

    if (puv.x < 0.0 || puv.x > 1.0)
    {
        return col;
    }

    if (puv.y < 0.0 || puv.y > 1.0)
    {
        return col;
    }

    return tex2D(_PackTex, puv);
}

フラグメントシェーダによって渡されたUV値(つまりパック用テクスチャ全体のUV値)を、C#側から設定したスケールとオフセットを用いて変換し、それを利用して書き込みを行っています。

UV0より下、あるいは1より上の場合は書き込もうとしているテクスチャの範囲外なので現在のテクスチャの色(※)を返します。

※ ... シェーダによる書き込みは、一度にすべてのテクスチャをパックするのではなく、「ひとつずつ」パックを行っています。その際、ダブルバッファを利用して交互にRenderTexuterを交換しながら書き込んでいるため、上記UV値をはみ出した部分は以前に書き込まれたテクスチャの色がある可能性があるためこうしています。

シェーダによる読み込み

さて、上記まででパックが終わりました。あとはそれを利用して画面に表示する必要があります。
つまりパックされたテクスチャから、該当のテクスチャ部分を読み出す必要があります。

読み出し用のシェーダは以下のようになっています。

fixed4 frag (v2f i) : SV_Target
{
    float2 uv = i.uv;
    uv /= _ScaleAndOffset.xy;
    uv += _ScaleAndOffset.zw;

    uv = clamp(uv, 0.0, 1.0);

    fixed4 col = tex2D(_MainTex, uv);
    return col;
}

書き込み用のシェーダと見比べてもらうと分かりますが、オフセット/スケールの計算が逆になっただけですね。

書き込み時は「オフセットさせて」から「スケールさせる(掛ける)」という手順でしたが、読み込み時は「スケールさせて(割って)」から「オフセットさせる」という手順で復元しています。

最後に

今回の実装では「回転」は考慮していません。
回転すればより良い位置があったとしても、それは考慮していない、ということです。

とはいえ冒頭に載せたようにそこまで隙間が目立つわけではないのでちょっとした利用くらいなら問題ないかなと思っています。

ちなみに読み込み/書き込みした動画がこちら↓

ちょっとした容量削減だったり、レンダリング負荷軽減目的なら意外と使えそうです。

World SpaceのCanvasにWorld SpaceからRaycastする

概要

ARやVRなどの開発を行っているとGUIも3D空間に配置する必要があります。
しかしUnityのuGUIは2Dで扱うことを想定しており、通常のGraphicRaycastはスクリーンスペースの位置から判定を行うものになっています。

つまり、3D空間に置かれた(World Spaceな)Canvasに対して、3D空間上にあるポインタ(例えば指やコントローラなど)からuGUIに対してRaycastしようとすると簡単には行きません。

そこで、ARやVRでも利用できるRaycastの仕組みを作らないとなりません。
以前、VRTKというフレームワークを参考に自分でEventSystemを拡張したのですが、その際に利用したWorld SpaceなCanvasへの3DオブジェクトからのRaycastをするために必要な部分について、(改めて必要になったので)メモとして残しておきたいと思います。

実際に動かしたやつはこんな感じです↓

処理フロー

まずはざっくり概観してみましょう。

  1. Canvasに配置されているGraphicすべてを判定対象にする
  2. Rayの方向とGraphicforward方向が同じかチェックする
  3. Rayの位置とGraphicの位置の距離を計算する
  4. Ray.GetPointを利用してワールドの位置を求める
  5. 求めたワールド位置をカメラのスクリーンスペース位置に変換する
  6. 変換した位置を元に、RectTransformがそれを含んでいるかチェックする
  7. 含んでいたらRaycast成功

コード

さて、全体の流れを説明したところでコードを見てみます。
コード自体はそこまで長くありません。

private void Raycast(Canvas canvas, Camera eventCamera, Ray ray)
{
    if (!canvas.enabled)
    {
        return;
    }

    if (!canvas.gameObject.activeInHierarchy)
    {
        return;
    }

    IList<Graphic> graphics = GraphicRegistry.GetGraphicsForCanvas(canvas);
    for (int i = 0; i < graphics.Count; i++)
    {
        Graphic graphic = graphics[i];

        if (graphic.depth == -1 || !graphic.raycastTarget)
        {
            continue;
        }

        Transform graphicTransform = graphic.transform;
        Vector3 graphicForward = graphicTransform.forward;

        float dir = Vector3.Dot(graphicForward, ray.direction);

        // Return immediately if direction is negative.
        if (dir <= 0)
        {
            continue;
        }

        float distance = Vector3.Dot(graphicForward, graphicTransform.position - ray.origin) / dir;

        Vector3 position = ray.GetPoint(distance);
        Vector2 pointerPosition = eventCamera.WorldToScreenPoint(position);

        // To continue if the graphic doesn't include the point.
        if (!RectTransformUtility.RectangleContainsScreenPoint(graphic.rectTransform, pointerPosition, eventCamera))
        {
            continue;
        }

        // To continue if graphic raycast has failed.
        if (!graphic.Raycast(pointerPosition, eventCamera))
        {
            continue;
        }

        Debug.Log($"Raycast hit at {graphic.name}", graphic.gameObject);
    }
}

すべてのGraphicを評価する

IList<Graphic> graphics = GraphicRegistry.GetGraphicsForCanvas(canvas);によって対象CanvasにあるGraphicオブジェクトのリストを得ることが出来ます。

ワールド空間の位置からスクリーンスペースの位置を求める

ワールド空間からスクリーンスペースの位置を求める方法です。
今回は特定オブジェクトがポインタの役割を果たして、そのオブジェクトからのRayを使いたいわけなので、少しだけ計算が必要となります。

コード部分は以下。

Transform graphicTransform = graphic.transform;
Vector3 graphicForward = graphicTransform.forward;
float dir = Vector3.Dot(graphicForward, ray.direction);

// Return immediately if direction is negative.
if (dir <= 0) { continue; }

float distance = Vector3.Dot(graphicForward, graphicTransform.position - ray.origin) / dir;

まず、対象Graphicforwardray.direction内積を取りその値で分岐します。
これはuGUIのオブジェクトのforward方向はこちらを向いている面の反対側になります。
そのため、Rayの方向と同じ場合を向いている場合(つまり内積が0より上の場合)だけ処理すればいいことになります。

なので内積結果が0以下だった場合は無視しているわけです。

そして同じ方向を向いていた場合はレイの位置とレイがぶつかった位置の距離を計算します。
コードにすると以下の部分。(分かりやすいように改行を入れています)

float distance = Vector3.Dot(graphicForward, graphicTransform.position - ray.origin);
distance /= dir;

Graphicforward方向」と「Graphicオブジェクトの位置からレイの位置を引いたベクトル」の内積を計算し、dirで割っています。

本来はVector3.Distance(from, to);で求めてもいいのですが、距離計算は負荷が高めなのと、上記計算で内積と、すでに計算済みのdirとの除算のみで求まるためそちらを利用しています。(dirを求める必要があるため一石二鳥、というわけです)

なぜこれで距離が求まるのかと言うと、差分ベクトルとの内積によってGraphicforward方向、すなわちGraphicが存在する平面との最短距離が求まります。

そしてそれを実際にGraphicが存在する位置の長さとするためには三角関数を利用します。
つまり、半径 / cos(θ)が実際に求めたい距離です。

そして実はdir内積の結果はまさにこのcos(θ)の値となっているため、それで割ることで距離が求まっていた、というわけです。

図にすると以下です。

f:id:edo_m18:20191220104604j:plain

RectTransformUtilityを利用して当たり判定

そして最後に、取得した全Graphicオブジェクトに対して当たり判定をしてやればOKなわけです。
当たり判定自体は関数が用意されていてRectTransformUtility.RectangleContainsScreenPoint(graphic.rectTransform, pointerPosition, eventCamera)を使います。

引数には対象となるGraphicRectTransform、スクリーンスペースのposition、そして判定対象となるCameraです。

対象オブジェクトにイベントを送る

今回は主にuGUIのオブジェクトをワールド空間から把握するためのもののメモなので、ちゃんとしたイベント周りの構築については以下の過去の記事を参考にしてみてください。

edom18.hateblo.jp

ただ、ざっくりでいいからイベントを送りたい、というケースもあるかと思います。
その場合は以下のように、Rayがヒットしたオブジェクトに対してイベントを送ってやればOKです。

// graphic object is detected by ray casting.
ExecuteEvents.Execute(graphic.gameObject, new BaseEventData(eventSystem), ExecuteEvents.submitHandler);
````

Compute ShaderとGraphics.DrawMeshInstancedIndirectを使ったレンダリングを理解する

概要

f:id:edo_m18:20191124215851p:plain

今回はGraphics.DrawMeshInstancedIndirectメソッドを使ってGPUパーティクルをレンダリングする方法をまとめます。
レンダリングに利用するパーティクルの位置計算はコンピュートシェーダで行います。

コンピュートシェーダについては過去に2つ記事を書いているので以下を参照ください。

edom18.hateblo.jp

edom18.hateblo.jp

ドキュメントはこちら。今回はこのドキュメントに沿って実装しつつ、できるだけミニマムに実装しました。

docs.unity3d.com

実際に動作するサンプルはGitHubにアップしてあります。

github.com

実際に動作している動画はこちら↓

ComputeShaderを準備する

まずはComputeShaderを用意します。
ComputeShaderを利用してGPUインスタンシングされたパーティクルの位置計算を行います。
今回実装するのはシンプルに、ランダムな位置から徐々に元の位置に戻るアニメーションです。

コード

ComputeShaderのコードを掲載します。

#pragma kernel ParticleMain

struct Particle
{
    float3 basePosition;
    float3 position;
    float4 color;
    float scale;
};

RWStructuredBuffer<Particle> _ParticleBuffer;

float _DeltaTime;

[numthreads(8, 1, 1)]
void ParticleMain(uint3 id : SV_DispatchThreadID)
{
    const int index = id.x;

    Particle p = _ParticleBuffer[index];
    p.position += (p.basePosition - p.position) * _DeltaTime;

    _ParticleBuffer[index] = p;
}

コンピュートシェーダ自体の説明は割愛します。詳細については上で紹介した記事を参照ください。
ここではGraphics.DrawMeshInstancedIndirectを利用する上で重要な点に絞って解説します。

システム全体で利用する構造体を定義する

今回はGPUパーティクルなので、パーティクルシステムで利用する構造体を定義します。
これはComputeShader内だけでなく、C#側のコードと、そしてGPUパーティクルをレンダリングするシェーダ側でも定義する必要があります。

Particle構造体は以下。

struct Particle
{
    float3 basePosition;
    float3 position;
    float4 color;
    float scale;
};

注意点は構造体のレイアウトをすべて同じにすることです。
構造体は宣言された順番にメモリにレイアウトされます。

この順番が異なっていると各プログラム間で整合性が取れなくなってしまうので気をつけてください。

処理を少しだけ説明しておくとRWStructuredBuffer<Particle>型のバッファから、カーネルParticleMain)の引数に渡されるスレッドID(SV_DispatchThreadIDセマンティクス)を使って計算対象となるParticle構造体を取り出します。

そして現在位置(Particle.position)から元の位置(Particle.basePosition)へ徐々に近づけていきます。

C#コード側でデータを準備し計算を実行する

次にC#側のコードを見てみましょう。
C#側で行うのは各データ(バッファ)の準備と計算の実行開始(Dispatch)、そしてコンピュートシェーダで計算されるバッファをシェーダ(マテリアル)に対して紐付けることです。

順に見ていきましょう。

コンピュートシェーダで利用するものと同じレイアウトの構造体を定義する

C#側でもコンピュートシェーダで利用するのと同じレイアウトの構造体を定義します。
構造体の定義は以下。

private struct Particle
{
    public Vector3 basePosition;
    public Vector3 position;
    public Vector4 color;
    public float scale;
}

若干 型が違っていますが配置の順番が同じ点に注目してください。(floatN型はそれぞれVectorN型になります)

各データ用変数を定義

C#側のタスクとして、データの準備と毎フレームごとの更新処理リクエストがあります。

ということでまずはデータの宣言について見てみましょう。

[SerializeField] private int _count = 10000;
[SerializeField] private ComputeShader _computeShader = null;

private ComputeBuffer _particleBuffer = null;
private ComputeBuffer _argBuffer = null;
private uint[] _args = new uint[5] { 0, 0, 0, 0, 0, };

コンピュートシェーダで利用する変数についてのみ抜粋しました。

これをセットアップしていきます。

データのセットアップ

データのセットアップをしているところを抜粋します。

List<Vector3> vertices = new List<Vector3>();
_targetMeshFilter.mesh.GetVertices(vertices);

Particle[] particles = new Particle[_count];

for (int i = 0; i < _count; i++)
{
    particles[i] = new Particle
    {
        basePosition = vertices[i % vertices.Count],
        position = vertices[i % vertices.Count] + Random.insideUnitSphere * 10f,
        color = _color,
        scale = Random.Range(0.01f, 0.02f),
    };
}

_particleBuffer = new ComputeBuffer(_count, Marshal.SizeOf(typeof(Particle)));
_particleBuffer.SetData(particles);

_computeShader.SetBuffer(_kernelId, "_ParticleBuffer", _particleBuffer);

今回のサンプルは指定したオブジェクト(メッシュ)の頂点位置からランダムに離れた位置にパーティクルを初期配置し、あとは元の頂点位置に徐々に戻っていくというものです。
そのためセットアップはターゲットとなるメッシュの各頂点位置を取得し、それを元に初期状態をバッファに設定しています。

バッファにセットする用の配列にデータを詰め込みそれを元にComputeBufferを生成、それをコンピュートシェーダにセットします。

Indirect(間接実行)用のバッファを用意する

ちょっとまだしっかりとは把握しきれていないのですが、このバッファはインスタンスを管理するために用いられるもののようです。
具体的にはインスタンス数や頂点数を渡し、適切な数頂点シェーダなどが起動するようにするもののようです。

ちなみに、凹みさんの記事から引用させてもらうと以下のように説明されていました。(Draw Procedural Indirectの箇所を参照)

ComputeBuffer.CopyCount() では追加した要素数を調べることが出来ます。この際、第 2 引数に渡している ComputeBufferType.IndirectArguments フォーマットのバッファは int 4 つ分のバッファです。GetData() をして一旦 CPU 側に値を持ってきてみると、何個の要素数が追加されたか確認できます。

tips.hecomi.com

そしてコード例として載っているのが以下。

var args = new int[4];
cbDrawArgs.GetData(args);
Debug.LogFormat("vert: {0}, {1}, {2}, {3}", args[0], args[1], args[2], args[3]);
// --> 939, 1, 0, 0

数値が格納されているのが第1要素と第2要素のみです。
さて、それを踏まえた上で今回のセットアップコードを見ると以下のようになっています。

int subMeshIndex = 0;

// Indirect args.
_args[0] = _targetMeshFilter.mesh.GetIndexCount(subMeshIndex);
_args[1] = (uint)_count;
_args[2] = _targetMeshFilter.mesh.GetIndexStart(subMeshIndex);
_args[3] = _targetMeshFilter.mesh.GetBaseVertex(subMeshIndex);

_argBuffer = new ComputeBuffer(1, sizeof(uint) * _args.Length, ComputeBufferType.IndirectArguments);
_argBuffer.SetData(_args);

第1要素がターゲットとなるメッシュのインデックス数、そして第2要素がパーティクルの数です。
ドキュメントには第3、第4要素も値が設定されているのですが試しにこれを0にしても正常に動きました。

なので3、4要素目はなにを意味しているのかちょっとまだ分かっていません。

分かっているのは、第1要素はインスタンスとしてレンダリングされる対象の頂点数、そして第2要素がインスタンスの数になることです。

このバッファを、DrawMeshInstancedIndirectを実行する際に引数に渡してやることで無事、レンダリングが行われます。

レンダリング用シェーダを作成する

データの準備、計算、更新処理が完成したら最後はパーティクルをレンダリングするシェーダを準備します。
それほど長くないのでまずはコード全体を掲載します。

Shader "Unlit/Particle"
{
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        LOD 100

        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            #include "UnityCG.cginc"

            struct Particle
            {
                float3 basePosition;
                float3 position;
                float4 color;
                float scale;
            };

            struct appdata
            {
                float4 vertex : POSITION;
            };

            struct v2f
            {
                float4 vertex : SV_POSITION;
                float4 color : COLOR;
            };

            StructuredBuffer<Particle> _ParticleBuffer;

            v2f vert (appdata v, uint instanceId : SV_InstanceID)
            {
                Particle p = _ParticleBuffer[instanceId];

                v2f o;

                float3 pos = (v.vertex.xyz * p.scale) + p.position;
                o.vertex = mul(UNITY_MATRIX_VP, float4(pos, 1.0));
                o.color = p.color;

                return o;
            }

            fixed4 frag (v2f i) : SV_Target
            {
                return i.color;
            }

            ENDCG
        }
    }
}

まず見てもらいたいのが、コンピュートシェーダとC#側で定義してきたParticle型の構造体があることです。
こちらも同様にレイアウトが同じになっています。

そして普段シェーダを書いているとあまり見慣れないStructuredBuffer<>という型が宣言されています。
これがまさにコンピュートシェーダで計算されたバッファが受け渡される変数です。

頂点位置を計算する

通常、シェーダはとあるオブジェクト(メッシュ)をレンダリングするために用いられます。
つまりレンダリングする対象はひとつのみです。

しかし今回は対象のメッシュはパーティクルひとつひとつを対象とし、それぞれがインスタンシングされるため複数オブジェクトを描くようにしないとなりません。

頂点シェーダのコードを抜粋してみましょう。

v2f vert (appdata v, uint instanceId : SV_InstanceID)
{
    Particle p = _ParticleBuffer[instanceId];

    v2f o;

    float3 pos = (v.vertex.xyz * p.scale) + p.position;
    o.vertex = mul(UNITY_MATRIX_VP, float4(pos, 1.0));
    o.color = p.color;

    return o;
}

普段見慣れないSV_InstanceIDというセマンティクスがついた引数があります。
これはGPUインスタンシングされた情報に対して何番目のオブジェクトかを表すIDとなります。
このIDを利用してコンピュートシェーダで計算された情報にアクセスします。
つまり以下の部分です。

Particle p = _ParticleBuffer[instanceId];

こうすることで起動された頂点シェーダが今、どのパーティクルを処理しているかが分かるわけです。

頂点の座標を変換する

次に対象パーティクルの頂点の座標変換処理を見てみましょう。

float3 pos = (v.vertex.xyz * p.scale) + p.position;
o.vertex = mul(UNITY_MATRIX_VP, float4(pos, 1.0));

v.vertex.xyzはパーティクルメッシュの頂点情報です。(今回はCubeを利用しています)
そしてその頂点をまずスケーリングし、そのあとで、コンピュートシェーダによって計算された位置を足し込んでパーティクルの位置としています。

Cubeの各頂点が一律同じように移動されるので、結果的にパーティクルとしてのCubeが移動するというわけですね。
そして実質、これがモデル行列を掛けた状態、つまりワールド座標空間での位置になります。

なのであとはビュー行列とプロジェクション行列を掛けてあげれば無事、クリッピング座標系に変換されるというわけです。

o.vertex = mul(UNITY_MATRIX_VP, float4(pos, 1.0));

として頂点シェーダの出力にしているわけですね。
これで無事、パーティクルが画面に表示されるようになります。

まとめ

シンプルな実装ならばコード量はさほど多くありません。
たくさんのパーティクルを出して絵を盛るもよし、複雑な処理をインスタンシングの力で打破するもよし。
使い方が分かれば色々と応用が効きそうです。

もしもっと複雑な処理をしたい場合は、以下のカールノイズの記事を参考に実装してみるといいかもしれません。

edom18.hateblo.jp

edom18.hateblo.jp

Unityで簡易的なドローイングツールを作ってみたので実装についてまとめ

概要

Unityで簡易的なドローイング機能がプロジェクトで必要になったので作ってみたもののまとめです。
実装はシェーダを利用して描いていて、ブラシと色とサイズを変更できるようになっています。

実際の動きはこんな感じ↓

ここで解説している機能のアセットはアセットストアで公開されているのでよかったら購入してください。

assetstore.unity.com

実装方針

実装方針は以下です。

  1. 描画対象となるImage要素上の位置を算出する
  2. 上記位置をImage要素内のUV値に変換する
  3. 該当位置に対して、設定されたブラシを描画する
  4. ブラシを描画する際は、さらにブラシ用のUV値を求める
  5. 描画ごとにバッファを差し替えて連続して描画する(ダブルバッファ)

という流れで実装しています。
実際のアセットではもう少し細かい調整を行っていますが、メインとなる処理は上記の通りです。

Image要素上の位置を算出する

Image(厳密にはRectTransform)上の位置を計算するには

public static bool ScreenPointToLocalPointInRectangle(RectTransform rect, Vector2 screenPoint, Camera cam, out Vector2 localPoint);

を利用します。

ドキュメント↓

docs.unity3d.com

実装は以下のような感じです。

private bool TryConvertToLocalPosition(Vector3 screenPos, out Vector2 pos)
{
    Camera camera = null;
    if (_drawImage.canvas.renderMode != RenderMode.ScreenSpaceOverlay)
    {
        camera = _camera;
    }

    return RectTransformUtility.ScreenPointToLocalPointInRectangle(_rectTrans, screenPos, camera, out pos);
}

RenderModeの判定を行っていますが、通常のモード(ScreenSpaceOverlay)の場合はCameraは関係ないのでnullを渡す必要がある点に注意してください。

この関数は、対象となる位置ベクトルがRectTrasnform内にある場合はtrueを返します。そして引数の最後で、ローカルの位置を格納して戻してくれます。
この位置ベクトルを次のUV値計算に利用します。

ちなみに、位置は関係なく、RectTrasnform内に入っているか、だけを判定する関数もあるので、入っているか否かだけを判断する場合はこちらの関数を利用するといいでしょう。

public static bool RectangleContainsScreenPoint(RectTransform rect, Vector2 screenPoint);
    // or
public static bool ScreenPointToLocalPointInRectangle(RectTransform rect, Vector2 screenPoint, Camera cam, out Vector2 localPoint);

要素上の位置をUV値に変換する

要素上の位置が求まったら、それをUV値に変換します。
といっても大した計算はしません。対象の要素(Image)のサイズでそれぞれ正規化してやればOKです。

// WidthとHeightはRectTransformのサイズ
// private float Width => _rectTrans.rect.width;
// private float Height => _rectTrans.rect.height;

private Vector4 NormalizePosition(Vector2 pos)
{
    float x = (pos.x + Width * 0.5f) / Width;
    float y = (pos.y + Height * 0.5f) / Height;
    return new Vector4(x, y, 0, 0);
}

注意点としては、RectTransform内の位置は中心を0として、プラスマイナスの値で返ってくるため、0 ~ 1の間になるように調整する必要があります。

ブラシのUVを算出する

位置が求められたので、次はブラシを描く位置とサイズを調整します。

と言ってもやっていることはむずかしくありません。
ひとつ問題なのは、描こうとしている対象のRenderTextureのUV値をそのまま利用することができません。
そのまま利用してしまうと、ただ単に、ブラシテクスチャをRenderTextureの範囲いっぱいに描いてしまうことになるからです。

なのでブラシ用に新しくUV値を計算する必要があるわけです。
といっても計算自体はそこまでむずかしくありません。

図にすると以下のようなイメージです。

f:id:edo_m18:20191110144113p:plain

今回の実装では以下のようにしました。

// 描きたい位置にまずオフセットする
float2 buv = i.uv - _Pos.xy;

// 次に、ブラシサイズになるようにUV空間をスケールする
buv *= _BrashSize;

// RenderTextureのアスペクト比を掛けて比率を調整する
buv.x *= _Aspect;

// 最後に、ブラシが位置の中央に描かれるようにオフセットする
buv += float2(0.5, 0.5);

オフセットのプラスマイナスは注意が必要です。図にすると以下のイメージ。

f:id:edo_m18:20191110144816p:plain

最後の計算が0.5になっているのは、この計算時点でブラシ用UV座標(0 ~ 1の範囲)に変換されているため、中心位置に変換するためには1.0の半分である0.5を足すだけで大丈夫です。

なお、ブラシUV座標に変換するためにブラシサイズを掛けていますが、ブラシサイズはC#側から計算して送っています。
具体的には以下のようにして計算しています。

private float NormalizeBrushSize(float brushSize)
{
    return Height / brushSize;
}

HeightRectTransformの高さですね。前の計算でアスペクト比があっているのでHeightを使っています。
そしてそれをブラシサイズで割った値を利用しています。これは、ブラシサイズになるように空間を「縮小」することで実現しています。

コードを見てみると計算自体は複雑ではありませんね。
ただ、UV座標を色々いじっていると頭の中が混乱することが多いです。
例えば、ブラシサイズを「小さく」しようとする場合は、スケールを「掛ける」必要があります。

uv *= 2.0; とすると「半分」のサイズになる。

よくよく考えれば当たり前なのですが、計算としては逆なのでたまに勘違いして計算を逆にしてしまうことがよくあります。

UVの計算については以前のこの記事を書いたときに、色々試していてだいぶ慣れてきました。

edom18.hateblo.jp

色を混ぜる

最後に、描こうとしている対象の色とブラシの色を混ぜ合わせて最終結果を作ります。
最初はたんに掛けたり足したりしていたのですが、それだけだと望んだ結果になりません。

最終的には、いわゆるアルファブレンディングと同じ方法でブレンドすることで解決しました。
コードは以下のような感じ。

fixed4 brash = tex2D(_BrashTex, buv);

// ブラシの黒いところに色をつけたいので反転させる
brash.rgb = 1.0 - brash.rgb;

// Blend SrcAlpha OneMinusSrcAlphaな感じでブレンドした色を返す
return lerp(col, (brash * _Color), brash.a);

このアセットに興味が出たらぜひ購入をお願いします!

assetstore.unity.com