UnityのCompute ShaderでCurl Noiseを実装（流体編）

概要

今回は「流体編」と題しました。
というのも、発表されたカールノイズ自体には剛体との当たり判定を行う部分の記載もあるのですが、そこがうまく動いてないのでまずは流体部分の説明にフォーカスして書いていこうと思います。

Curl Noice。なんとなく流体っぽく見えるノイズ。
Curlとは「カール」、つまり「渦巻く」という意味からつけられた名前だそう。

（自分の浅い知識ながら）もう少し具体的に言うと、ベクトル場を回転（Curl）させることで作るノイズ、ということで「Curl Noise」というと思っています。

いつもお世話になっている物理のかぎしっぽさんで、この「ベクトル場の回転」についての記載があります。

ちなみに今回の主な実装は以下の記事を参考にさせていただきました。

qiita.com

また、発表された論文はこちら（PDF）

ソースコード

今回のサンプルはGithub上にUnityプロジェクトとしてアップしてあります。興味がある方はDLして見てみてください。

github.com

（自分の顔アイコンになっちゃうのどうにかならんか・・）

今回の記事で書くこと

今回の記事は、主にCurl Noiseについて触れます。
サンプルはComputeShaderを使ってUnityで実装しました。ComputeShaderによる実装についても少し触れます。

ただ、ComputeShaderなどの使い方などは知ってる前提であまり詳しくは解説しません。
ComputeShader自体については前の記事を読むなどしてください。

大まかな解説

自分はあまり物理に明るくありませんが、今回の実装を通して色々と学びがあったので、それをまとめる意味でも少しだけ解説を書きたいと思います。
(注) あくまで個人的な理解なので、間違っている可能性があります。

非圧縮性を確保する

流体が流体っぽく見える所以は、この非圧縮性が必要だという認識です。
非圧縮性の説明をWikipediaから引用させてもらうと、

非圧縮性流れ（ひあっしゅくせいながれ）とは流体力学において、流体粒子の内部で密度が一定の流体である。縮まない流体とも呼ばれる[1][2]。連続体力学における非圧縮性の概念を流体に適用したものである。
言い換えると、非圧縮性とは流体の速度の発散が 0 になることである（この表現が等価である理由は後述）。

とのこと。
ここで、流体の速度の発散が0になることであるとあります。

物理のかぎしっぽさんの記事のdivから、「発散」についての説明を引用させてもらうと、

スカラー場の勾配を考えたとき，ベクトル微分演算子 $\nabla = \left(\frac{\partial}{\partial x_{1}}, \frac{\partial}{\partial x_{2}}, \frac{\partial}{\partial x_{3}}\right)$ というものを導入しました．そして，この $\nabla$ をスカラー関数に作用させたものを勾配（ ${\rm grad}$ ）と呼びました． $\nabla$ をベクトル関数と内積を取る形で作用させたものを発散と呼びます．英語で発散を divergence と言うので，記号 ${\rm div}$ を使う場合もあります．

さらに別の箇所から引用すると、

ベクトル場を，水（非圧縮流体）の流れだと考えると状況がイメージしやすいでしょう．この直方体領域は流れの中に置かれていますから，絶えず水が流れ込んだり出て行ったりしています．しかし，水は非圧縮流体だと仮定していますので，普通なら，入ってくる水量と出て行く水量は同じはずです．

と書かれています。
つまり、（自分の理解では）この「発散が0になる」とは、流れの中で流入と流出が等しく起こる＝差分がゼロ、ということだと思います。

そしてカールノイズはこの「発散0を実現してくれるための手法」ということのようです。
理由として自分が理解しているのは、以下のように考えています。
まず言葉で説明すると、

あるベクトル場から得たベクトルを「回転」させ、「内積」を取った結果がゼロである。

ということ。式にすると、

$\nabla \cdot \nabla \times \vec{v} ≡ 0$

カール（回転）は、どうやら90°回転させることのようです。
つまり、「とあるベクトルを90°回転させたベクトルとの内積＝0」ということです。
そもそも内積は、垂直なベクトルとの計算は0になりますよね。

なので、ある意味0になるのは自明だった、というわけです。（という認識です）

発散の定義が、 $\nabla$ との内積を取ったものが $0$ になる、ということなので、「だったら最初からゼロになるように回転しておこうぜ」というのが発想の基点なのかなと理解しています。

実際にUnity上で実行した動画はこんな感じ↓

ひとまずカールノイズ（衝突判定なし）のキャプチャ pic.twitter.com/b4boyMaBFC
— edom18@AR / MESON (@edo_m18) 2018年1月17日

パーティクルの位置の計算についてはComputeShaderを使って計算しています。
ComputeShaderについて前回2回に分けて記事を書いているので、よかったら見てみてください。

edom18.hateblo.jp

それから、描画周りの実装については、凹みTipsで過去に凹みさんが書かれていた以下の記事を参考に実装しています。
~~（というかほぼそのままです）~~
いつもありがとうございます。

tips.hecomi.com

用語解説

いくつか数学・物理的なワードが出てくるので、簡単にまとめておきます。

ベクトル場

ものすごくざっくりいうと、「ベクトル」がちらばった空間です。
分かりやすい例で言うと、風速ベクトルを空間にマッピングしたもの、と考えるといいかもしれません。

風速はその場所その場所で当然向きや強さが違います。
そしてその瞬間の風速ベクトルを記録したのがベクトル場、という感じです。

プログラム的に言うと、ベクトルを引数に取ってベクトルを返す関数、と見ることが出来ます。（厳密に言うと風速は時間によっても変化するので、時間と位置による関数となります）

float3 windVec = GetWind(pos, time);

みたいな感じですね。
今回のCurlNoiseについては時間的変化は考えず、以下のように速度ベクトルを求めています。

float3 x = _Particles[id].position;
float3 velocity = CurlNoise(x);

ベクトル場については、物理のかぎしっぽさんのところにとても分かりやすい例と図があるので、そちらを見てみてください。

カール（Curl・rot）

ベクトル場の回転を表す。
前述の、物理のかぎしっぽさんの記事から引用させてもらうと、

三次元ベクトルの場合，ナブラを外積を取る形で作用させることも可能です．これを回転と呼び， $\nabla \times$ ， rot もしくは curl という記号で表現します．この記事では，回転について基本的な意味や性質を考えます．

ramp関数

ramp関数 | Wikipediaとは、グラフが斜路（ramp）に見えることから名付けられたそう。

論文で書かれている式をグラフにしてみると以下の形になります。
なんとなく「斜路」という意味が分かりますね。

f:id:edo_m18:20180117234619p:plain

本題

さて、ここからが、今回実装した内容の解説です。
まず、Curl Noiseについてざっと解説してみます。（自分もまだ理解が浅いですが）

前述のように、ベクトル場を（ベクトル解析でいう）回転させることによって得られるノイズを「Curl Noise」と呼んでいるようです。
元々は、アーティスティックな流体表現を実現するために考案されたもののようです。

物理のかぎしっぽさんの記事を参考に式を表すと以下のようになります。

f:id:edo_m18:20171014165055p:plain

$\nabla$ は「ナブラ」と読み、「微分演算子をベクトルに組み合わせたナブラというベクトル演算子」ということのようです。

記号を見ると偏微分の記号なので、要は、ベクトル場の各ベクトルに対してちょっとずつ位置をずらして計算したもの、と見ることができます。

CurlNoiseを生成している関数を見てみると、以下のようになります。

float3 CurlNoise(Particle p)
{
    const float e = 0.0009765625;
    const float e2 = 2.0 * e;
    const float invE2 = 1.0 / e2;

    const float3 dx = float3(e, 0.0, 0.0);
    const float3 dy = float3(0.0, e, 0.0);
    const float3 dz = float3(0.0, 0.0, e);

    float3 pos = p.position;

    float3 p_x0 = SamplePotential(pos - dx, p.time);
    float3 p_x1 = SamplePotential(pos + dx, p.time);
    float3 p_y0 = SamplePotential(pos - dy, p.time);
    float3 p_y1 = SamplePotential(pos + dy, p.time);
    float3 p_z0 = SamplePotential(pos - dz, p.time);
    float3 p_z1 = SamplePotential(pos + dz, p.time);

    float x = (p_y1.z - p_y0.z) - (p_z1.y - p_z0.y);
    float y = (p_z1.x - p_z0.x) - (p_x1.z - p_x0.z);
    float z = (p_x1.y - p_x0.y) - (p_y1.x - p_y0.x);

    return float3(x, y, z) * invE2;
}

eが、「ちょびっと動かす」部分の値で、それを各軸に対して少しずつ変化させる量としてdx、dy、dzを定義しています。
そしてそれぞれ変化した値をSamplePotential関数に渡し、さらにそれの差分を取ることで計算を行っています。

ベクトル場の回転について、論文に書かれているものは以下です。

$\vec{v}(x, y, z) = \left( \frac{∂\psi_3}{∂ y} − \frac{∂\psi_2}{∂ z}, \frac{∂\psi_1}{∂ z} − \frac{∂\psi_3}{∂ x}, \frac{∂\psi_2}{∂ x} − \frac{∂\psi_1}{∂ y} \right)$

ベクトル微分演算子を作用させるので、微分を取る目的で以下のように計算しています。

float3 p_x0 = SamplePotential(pos - dx, p.time);
float3 p_x1 = SamplePotential(pos + dx, p.time);
float3 p_y0 = SamplePotential(pos - dy, p.time);
float3 p_y1 = SamplePotential(pos + dy, p.time);
float3 p_z0 = SamplePotential(pos - dz, p.time);
float3 p_z1 = SamplePotential(pos + dz, p.time);

それぞれの方向（x, y, z）に対してちょっとずつずらした値を取得しているわけですね。
そして続く計算で、前述の「ベクトル場の回転」を計算します。

float x = (p_y1.z - p_y0.z) - (p_z1.y - p_z0.y);
float y = (p_z1.x - p_z0.x) - (p_x1.z - p_x0.z);
float z = (p_x1.y - p_x0.y) - (p_y1.x - p_y0.x);

[2018.06.29追記]
すみません、計算ミスをしていました。以前は以下のように書いていましたが、最後のパラメータの微分の計算が「マイナス」になるのが正解です。

// 間違っていた記述
float x = (p_y1.z - p_y0.z) - (p_z1.y + p_z0.y);
float y = (p_z1.x - p_z0.x) - (p_x1.z + p_x0.z);
float z = (p_x1.y - p_x0.y) - (p_y1.x + p_y0.x);

↑カッコでくくった右側の計算が「プラス」になっていました。

前述の数式と見比べてもらうと同じことを行っているのが分かるかと思います。

誤解を恐れずに言えば、カールノイズのコア部分はここだけです。
つまり、

ベクトル場からベクトルを得る
得たベクトルを回転させる
得たベクトルを速度ベクトルとして計算する
速度ベクトルから次の位置を決定する

だけですね。

あとはこの、算出された回転済みのベクトル（速度）をパーティクルに与えるだけで、なんとなく流体の流れっぽい動きをしてくれるようになります。

ノイズを生成する

さて、カールノイズをもう少し具体的に見ていきます。
カール「ノイズ」というくらいなので、当然ノイズを生成する必要があります。

色々なカールノイズのサンプルを見ると、シンプレクスノイズかパーリンノイズを利用しているケースがほとんどです。
論文でもパーリンノイズに言及しているので、今回はパーリンノイズを利用してカールノイズを作成しました。

ノイズを生成する上で、論文には以下のように書かれています。

To construct a randomly varying velocity field we use Perlin noise $N(\vec{x})$ in our potential. In 2D, ψ = N. In 3D we need three components for the potential: three apparently uncorrelated noise functions (a vector $vec{N}(\vec{x})$ do the job, which in practice can be the same noise function evaluated at large offsets.

Note that if the noise function is based on the integer lattice and smoothly varies in the range [−1,1], then the partial derivatives of the scaled N(x/L) will vary over a length-scale L with values approximately in the range O([−1/L,1/L]). This means we can expect vortices of diameter approximately L and speeds up to approximately O(1/L): the user may use this to scale the magnitude of ψ to get a desired speed.

ざっくり訳してみると、

ランダムな変化するベクトル場を構築するのに、パーリンノイズ（ $\vec{x})$ をポテンシャルのために利用する。2Dの場合は「 $\psi = N$ 」。

3Dの場合はポテンシャルのために3つの要素が必要となる。
3つの、明らかに無関係なノイズ関数（ $\vec{N}(\vec{x})$ を使用する。それらは同じノイズ関数を、大きなオフセットを持たせたもの。

Note: もし、ノイズ関数が整数格子にもとづいており、[-1, 1]の範囲で滑らかに変化するのであれば、 $N(x/L)$ にスケーリングされた偏微分は、おおよそO([−1/L,1/L])の範囲の値でスケールLに渡って変化する。

これは、ほぼ直径Lの渦とO(1/L)までの速度を期待することができることを意味する。ユーザは $\psi$ の値にスケールされた（望みの）スピードを使えるかもしれない。

ということで、ノイズとスケールを考慮した処理は以下のようになります。

float3 SamplePotential(float3 pos)
{
    float3 s = pos / _NoiseScales[0];
    return Pnoise(s);
}

_NoiseSclaesは、いくつかのスケール値が入った配列です。
が、今回はひとつのみ利用しているので_NoiseScales[0]だけを使っています。

Pnoiseはパーリンノイズを利用したベクトル場の計算です。
実装は論文に書かれているように、大きなオフセットを持たせた同一のパーリンノイズ関数を利用して3要素を計算しています。

// パーリンノイズによるベクトル場
// 3Dとして3要素を計算。
// それぞれのノイズは明らかに違う（極端に大きなオフセット）を持たせた値とする
float3 Pnoise(float3 vec)
{
    float x = PerlinNoise(vec);

    float y = PerlinNoise(float3(
        vec.y + 31.416,
        vec.z - 47.853,
        vec.x + 12.793
        ));

    float z = PerlinNoise(float3(
        vec.z - 233.145,
        vec.x - 113.408,
        vec.y - 185.31
        ));

    return float3(x, y, z);
}

ちなみにPerlinNoiseは以前、JavaScriptで実装したときのものをCompute Shaderに移植したもので、以下のように実装しています。

float PerlinNoise(float3 vec)
{
    float result = 0;
    float amp = 1.0;

    for (int i = 0; i < _Octaves; i++)
    {
        result += Noise(vec) * amp;
        vec *= 2.0;
        amp *= 0.5;
    }

    return result;
}

float Noise(float3 vec)
{
    int X = (int)floor(vec.x) & 255;
    int Y = (int)floor(vec.y) & 255;
    int Z = (int)floor(vec.z) & 255;

    vec.x -= floor(vec.x);
    vec.y -= floor(vec.y);
    vec.z -= floor(vec.z);

    float u = Fade(vec.x);
    float v = Fade(vec.y);
    float w = Fade(vec.z);

    int A, AA, AB, B, BA, BB;

    A = _P[X + 0] + Y; AA = _P[A] + Z; AB = _P[A + 1] + Z;
    B = _P[X + 1] + Y; BA = _P[B] + Z; BB = _P[B + 1] + Z;

    return Lerp(w, Lerp(v, Lerp(u, Grad(_P[AA + 0], vec.x + 0, vec.y + 0, vec.z + 0),
                                    Grad(_P[BA + 0], vec.x - 1, vec.y + 0, vec.z + 0)),
                            Lerp(u, Grad(_P[AB + 0], vec.x + 0, vec.y - 1, vec.z + 0),
                                    Grad(_P[BB + 0], vec.x - 1, vec.y - 1, vec.z + 0))),
                    Lerp(v, Lerp(u, Grad(_P[AA + 1], vec.x + 0, vec.y + 0, vec.z - 1),
                                    Grad(_P[BA + 1], vec.x - 1, vec.y + 0, vec.z - 1)),
                            Lerp(u, Grad(_P[AB + 1], vec.x + 0, vec.y - 1, vec.z - 1),
                                    Grad(_P[BB + 1], vec.x - 1, vec.y - 1, vec.z - 1))));
}

パーリンノイズ自体の詳しい内容は、調べるとたくさん記事が出てくるのでそちらをご覧ください。

カールノイズのカーネル

さて最後に、カールノイズを計算しているカーネル（Compute Shaderのmain関数のようなもの）を見てみます。

[numthreads(8, 1, 1)]
void CurlNoiseMain(uint id : SV_DispatchThreadID)
{
    float3 pos = _Particles[id].position;

    float3 velocity = CurlNoise(_Particles[id]);

    _Particles[id].velocity = velocity * _SpeedFactor * _CurlNoiseIntencity;
    _Particles[id].position += _Particles[id].velocity * _DeltaTime;

    _Particles[id].time += _DeltaTime;
    float scale = 1.0 - (_Particles[id].time / _Particles[id].lifeTime);
    if (scale < 0)
    {
        _Particles[id].active = false;
        return;
    }

    _Particles[id].scale = scale;
}

_Particlesがパーティクルの配列になっていて、Particle構造体がそれぞれのパーティクルの位置や速度などを保持しています。
それを、毎フレームごとに更新し、それをレンダリングすることで流体のような動きを実現している、というわけです。

現在リポジトリにアップしているものは、これに加えて、剛体との衝突判定の実装を試みているのでもう少し複雑な状態になっていますが、衝突判定なしで流体のような動きだけがほしい場合は以上の実装で実現することができます。

ガチの流体計算に比べるとだいぶ簡単な処理で実現できているのが分かってもらえたかと思います。

次回予告 - 衝突判定

さて、冒頭でも書きましたが、このカールノイズ。
簡単な衝突判定を利用して、球体に沿わしたり、といったことができます。

やったー。カールノイズの衝突判定もちゃんと動いた。 pic.twitter.com/lMC8pgHjrc
— edom18@AR / MESON (@edo_m18) 2018年6月29日

球体が干渉しているのが分かるかと思います。
衝突判定が取れるようになれば、さらに流体感を増した表現になるので、ぜひマスターしたいところです。

e.blog

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