Niagara SimulationStageでSPHを実装する

# 前置き
https://qiita.com/marv_kurushimay/items/1dc485762c06bdef1c75

で
> 現在、こんな感じのSPH(Smoothed Particle Hydrodynamics)を使った表現を、UE4/Niagaraで実装すべく検証中です。
この内容については追々まとめる予定ですが、

なんて書いておいて、半年以上経ってしまいました。

https://qiita.com/marv_kurushimay/items/3e206a190c0646234a96

で解説した手法をUE4のNiagara SimulationStageを使って実装していきます。

だいぶ前に検証した内容なので、UnrealEngineのバージョンは __4.26__ です。

# 概要

細かい所は置いておいて、大まかな流れを示します。

1. SPH向けの係数などを初期化
1. ParticleをNeighborGrid3Dに入れる
1. コリジョン判定
1. サブステップループ
    1. 密度計算
    1. 圧力計算
    1. 加速度計算
    1. 加速度を速度/位置に適用
1. Grid2D上にラスタライズ
1. バイラテラルフィルタを使って、深度値などをぼかす
1. Grid2Dを元にマテリアルを使ってEmitterSpriteを描画

# １：SPH向けの係数などを初期化
![InitCoeff.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/4871983c-58b8-7687-c826-2005bec3b883.png)
ここでは、

* 質量
 * 固定値でも良いのですが、想定体積/パーティクル数で求めることである程度調整しやすくしています。
 * また、パーティクル毎の質量差は考慮せず、全てのパーティクルは均一な質量を持つようにしています。
*  各重み関数の係数
 * 	密度計算に使う$W_{poly6}(r)$
 * 	圧力項の計算に使う$\nabla W_{spiky}(r)$
 * 	粘性項の計算に使う$\Delta W_{visc}(r)$
* 影響範囲の2乗
 * 計算中に結構な頻度で使用するので、事前に変数化しておきます
* サブステップの経過時間(delta time)
 * 経過時間をサブステップ数で割るだけ

以上の各値を求めています。なお重み関数の係数については、粒子ごとの質量が一定であれば先に質量を乗算しておけるので、ここで乗算しておきます。

# ２：ParticleをNeighborGrid3Dに入れる
![FillNeighborGrid3D_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/c25f5525-1808-f593-6228-e6b0df4eba8e.png)
SimulationStageの設定には主に __IterationSource__ と __Iterations__ の二つがあります。
__IterationSource__ は、各モジュール（スクラッチパッド）がどの単位で実行されるか。 __Iterations__ は何回繰り返すか。です。

SimulationStageが強力である理由の一つとして、このイテレーションの仕組みにあります。この項目では、 __IterationSource__ はParticle、 __Iterations__ は1と普通にパーティクル毎に1回実行する設定になっています。
しかし、 __IterationSource__ をRenderTargetとすれば入力したRenderTargetのピクセル毎に実行しますし、 __IterationSource__ はParticleのまま __Iterations__ を10にすれば、パーティクル毎の実行を10回繰り返す事になります。 __Iterations__ があるおかげで、サブステップの処理が非常に簡単に組めるようになっています。

パーティクルをNeighborGrid3Dに入れる処理は以下の通りです。
![FillNeighborGrid3D_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/dbceacf6-1cea-75da-2db4-6c97f845ddeb.png)
パーティクルの位置を入力にNeighborGrid3D上のセル番号を計算します。SetParticleNeighborCount()やSetParticleNeighbor()は、引数にLinearIndexを必要としているので、適宜変換を行います。

# ３：コリジョン判定
## 主な設定
![MarvCollide.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/ed3c6e79-345f-1d5c-8a39-4929a430d10f.png)
標準のCollisionモジュールが使えます。が、そのままではちょっと問題が出たので、少しだけ改造しました。基本的には使い方は一緒なので、詳しくは公式ドキュメントなどに譲りますが、設定で私が変更したポイントを先に触れます。
### GPU Collision Type
GPU Distance Fieldsにした理由は、簡単にメッシュ形状のコリジョンを扱える為です。プロジェクト設定でディスタンスフィールドを有効にする必要がありますが、テスト段階でコリジョンメッシュを用意するのも大変だったので、ディスタンスフィールドを採用しました。
### Collision Radius - Radius Calculation Type
コリジョンの半径をどのパラメータから取得するかを設定します。
今回は、Particles.Scaleにサイズを設定していたので、__Mesh__ を選択しました。
### Rest - Enable Rest State
私も正確な意図を掴みかねているのですが、コリジョンにヒットした時、貫通させず表面上で停止させるかどうかのオプションです。
次の項で説明する改造ポイントとの兼ね合いもあって、ここではOFF（表面で停止させない）にしました。

## 改造ポイント
### デフォルトのコリジョンモジュールで都合が悪い点
![Collision_Vanish.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/d0c61ece-9a9a-3b66-7a46-e9cfced9e256.gif)
このように標準のコリジョンモジュールだと、コリジョン判定させるオブジェクトを動かすたびにパーティクルがどんどん消えていきます。これは、DistanceFiledのコリジョン判定で、メッシュ内にある（Inside a mesh）と判定されたパーティクルが消されているためです。
メッシュ内に食い込んだパーティクルはどっちの方向に流すべきか、またどれぐらいの力で跳ね返すべきか、計算が難しい場合があり、それならいっその事消してしまえ。という事の様です。
パーティクルが継続的に発生している状況では良いのかもしれませんが、今回の場合は最初に発生させたパーティクルを使い続けたいので、消えてしまっては困ります。
### 修正点
![CollisionInsideamesh.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/e1374378-6eff-5ee4-c03d-20c6cd8a936b.png)
これを修正するには、2つの標準モジュールをコピーして使う必要がありました。
1つは、Niagaraから直接参照する __Collision__ モジュール。もう１つは、Collisionモジュールから参照されている __CollisionQueryAndResponse__ モジュールです。
上記の画像はCollisionQueryAndResponseモジュールの修正する部分です。変更前は赤で書き込んだ接続があって、Inside a meshならAliveをfalseにして、パーティクルを消していました。

上記のCollisionQueryAndResponseモジュールの修正と、Collisionモジュールで参照されているCollisionQueryAndResponseモジュールを今回コピーした物に差し替えることで、この問題を回避できます。
個人的には、Collisionモジュールのオプションとして正式につけてくれないかなと思うぐらいの変更点ですが、残念ながらそうではないので、独自の変更が必要でした。

# ４：サブステップループ
ここからがSPH処理の肝になる部分です
![Substep_1.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/0e3b0dec-6ce9-e9ca-d22b-61e160291ba3.png)
前の項目で説明した通り、Iterationsを10にして1フレームに10回サブステップを処理する事にしています。
本当はIterationsに変数を指定したい所ですが、少なくとも4.26ではできませんでした。そのため、サブステップの経過時間を計算した時に使った値と手動で一致させる必要があります。
それを逆手にとって、サブステップの経過時間の計算はそのままにIterationsを減らすとスロー再生が可能になります。粒子がすぐにどこかへ飛んで行ってしまう場合など、デバッグする時には重宝しました。

## 密度計算
![CalcDensity.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/891a56c9-6557-0402-0716-798adcdb4df9.png)
近傍の粒子（パーティクル）に対して計算を行いたいので、$3\times3\times3=27個$のNeighborGrid3Dセルに対して処理を行い、総和を求めています。もし、2Dで処理を行う場合は、$3\times3＝9個$のセルに対して行う事になります。

基本的に理論編の通りの実装ですが、重み関数$W_{poly6}(r)=c\left(h^2−|r|^2\right)^3$の辺りだけちょっと改変しています。

```hlsl
float DiffLenSq = 1.0 - LengthSq/InSmoothLenSq;
OutDensity += InDensityCoeff * DiffLenSq * DiffLenSq * DiffLenSq;
```
と、$(h^2-|r|^2)$の部分を$\left(1-\frac{|r|^2}{h^2}\right)$に変えました。この理由は、元の式のままだと粒子の数や対象エリアのスケールを変えた時に、密度の変化が激しくて調整が難航したためです。
そこで、差分ではなく __比率__ にすることで調整しやすくなったので採用しました。

## 圧力計算
![Substep_CalcPressure.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/a76a2601-0f8a-2b25-7ab2-97071ebb8e10.png)
これも、理論編の理想気体の圧力状態方程式 $p_i=k(\rho_i-\rho_0)$ の通りに、密度から圧力を求めています。係数$k$はPressureStiffness、静止密度$\rho_0$はRestDensityというパラメータで入力しています。

## 加速度計算
![CalcAccelaration.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/bd001053-e256-dcbf-1576-1e9ce7ace0f9.png)
また長いCustomHLSLが出てきましたが、近傍の粒子に対して行う処理など基本的には密度計算と一緒です。近傍の粒子に対して理論編で紹介した圧力項と粘性項を計算して総和を求めています。

なお、ここでも重み関数に使用する距離の差分は、 __比率__ に変えました。

```hlsl
float Length = length(DiffPos);
float diffLength = 1.0 - (Length / InSmoothLength);
```
この辺、理論を厳密にとらえると正しくは無いと思いますが、動作的に問題なければOKとしています。

## 加速度を速度/位置に適用
ここまでの計算で、各粒子の加速度が求まりました。
![Substep_CalcPosition.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/12e74878-7aaf-f8b9-edaf-1791fbe97534.png)
これも理論編の通り、加速度から速度、速度から位置を計算して更新します。

# ５：Grid2D上にラスタライズ
![Grid2D_DrawSpheres.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/ddfcb1e0-db82-3dde-1fc3-505a2cda8e1b.png)
Grid2D上にラスタライズ（レンダリング）する手法がいくつか考えられますが、ひとまず対Sphereのレイチェックで実装してみました。この結果はGrid2Dの要素ごとに

* ワールド位置
* 深度
* ワールド法線
* デバッグ用カラー
* OpacityMask

を出力しています。
これだと正確にラスタライズできるのですが、いくつか懸念点があります。

## 問題点
### 負荷の問題
この状態では、Grid2Dの各要素に対してパーティクル数分のループを回しているので、負荷がかなり高いです。その為、最終的にはここでもNeighborGrid3Dを使い、チェック対象パーティクルの絞り込みを行いました。
しかし、それでもまだ負荷が高いので何か対策が必要になっています。

なお、実装内容が分かりづらくなるので、前段の画像はNeighborGrid3Dの対応を入れる前の状態です。
### 法線がきつすぎる
正確な法線を計算して出力すると、法線がきつすぎてビー玉や吸水ポリマーが転がっている様にしか見えませんでした。その為、この例ではそれを緩和する処理を入れています。

```hlsl
OutVelocityMask = saturate((length(Velocity)-InVelocityRange.x)/(InVelocityRange.y-InVelocityRange.x));
OutVelocityMask = max(OutVelocityMask, saturate((OutParticlePosition.z-200.0)));
vNormal.xy *= OutVelocityMask;
OutNormal = (vNormal.z > 1e-6) ? normalize(vNormal) : float3(0, 0, 1);
```
![SPH_CorrectNormal.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/c284e140-2b7e-30b5-7e08-4ce8762e1d3f.png)
上図は静止状態の水面ですが、左半分は補正処理を入れる前。右半分が補正処理を入れた後です。分かりやすくするために、次で説明するバイラテラルフィルタはOFFにしています。
このように、元の法線では静止状態でも水面が波打ってしまい、コースティクスなどが見えません。その為、Velocityを参照して一定速度以下なら徐々に法線の強度を下げて平面化するようにしています。

:::note info
このセクションで出したスクショですが、同じパラメータを用意しても配布版の4.26では同じ見た目になりません。SingleLayerWater関連で屈折の計算とPixelDepthOffsetの処理にバグがあり、その結果正しい見た目を得られなくなっています。

具体的には、
<a href="https://issues.unrealengine.com/issue/UE-102259">UE-102259 SingleLayerWater with PixelDepthOffset gives wrong refraction</a>
<a href="https://issues.unrealengine.com/issue/UE-130479">UE-130479 AbsoluteWorldPosition and PDO not working correctly in SingleLayerWater calculations</a>
でバグ対応されましたが、残念ながらどちらも4.27ではなく<strong>5.0</strong>からの適用になっています。
:::

# ６：バイラテラルフィルタを使って、深度値などをぼかす
![BilateralFilter.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/c1d106b9-7fed-c0ae-8288-c5b4daa3ce5f.png)
Sphereを描画したそのままだと、あまりに球体然としすぎるので、バイラテラルフィルタでぼかします。
バイラテラルフィルタの概要は [Wikipedia](https://en.wikipedia.org/wiki/Bilateral_filter) を見てもらうとして、基本的にガウシアンフィルタです。ただ、ぼかし処理を行う二次元座標(Spatial Domain)と、その深度値の差(Range Domain)の2要素からフィルタを行うという点が画期的で、深度値に差がある部分はぼかさない事で輪郭を維持したままぼかし処理を行う事が出来ます。

![BilateralFilter.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/b3dc2d33-314f-d53e-b2a9-4bc5dd398681.png)
結果はこの通り、一目瞭然！
ちなみに、この結果画像では前述の法線の補正はOFFになっています。バイラテラルフィルタがONの状態でも、静止している水面が歪んでコースティクスが奇麗に出ていません。これにリフレクション要素が入ってくるとより顕著に歪みが見えてきますので、やはり前述の法線の補正も有った方が良さそうです。


# ７：Grid2Dを元にマテリアルを使ってEmitterSpriteを描画
ようやく最後の項目です。ここまでの処理でGrid2Dには法線とOpacityMaskなどがありますので、「[Grid2Dの中身をマテリアルで参照する方法](https://qiita.com/marv_kurushimay/items/1dc485762c06bdef1c75#grid2d%E3%81%AE%E4%B8%AD%E8%BA%AB%E3%82%92%E3%83%9E%E3%83%86%E3%83%AA%E3%82%A2%E3%83%AB%E3%81%A7%E5%8F%82%E7%85%A7%E3%81%99%E3%82%8B%E6%96%B9%E6%B3%95)」で紹介した方法で描画を行います。
ただ、他の項目もそうですが、あくまで私はこう実装した。という感じなので、他に実装の仕方はあると思います。各自の実装の踏み台程度に考えてください。

![Material_Water.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/431039/b9b3c315-9cf3-36ae-e9a8-72eada88cddd.png)
マテリアルの全景はこんな感じです。社内で構築した最終的なマテリアルから当り障りのない部分を抜き出したので、若干いびつな所があると思いますが、ご了承ください。
水という事で、ShadingModelは最近追加されたSingleLayerWaterを使いました。Spriteとして描画しているので、不要な部分を除去する為BlendModeはMaskedに。あと、Grid2Dで作成した法線はワールド法線なので、Tangent Space NormalのチェックボックスをOFFにします。
なお、SingleLayerWaterの説明は公式動画の[4.26 Water And Volumetrics Preview | Inside Unreal](https://www.youtube.com/watch?v=XHugDXx6c5c)が入門編としてはちょうどいいです。マテリアルについては21分辺り、コースティクスについては、2時間11分辺りから解説があります。

__WaterAttributes__ と __WaterCuaustics_SubUVAnimation_Custom__ のオリジナルは、Waterプラグインにあります。今回は、Grid2Dで生成した情報から計算したPixelDepthを使うために、マテリアル関数の入力としてPixelDepthを指定できるようカスタマイズしました

# 今後の改善に向けたネタ
## 残像を付与

https://youtu.be/-cKgZrrBJ2w?t=1624

上記は2020年に中国で行われたUnrealOpenDayの動画です。27:04付近から解説されている残像機能を実際に実装してみましたが、これがあると粒子が水滴っぽく見えるのでお勧めです。
実装の際は、以前の記事で書いた[ParameterMapSetの実行順](https://qiita.com/marv_kurushimay/items/1dc485762c06bdef1c75#parametermapset%E3%81%AE%E5%AE%9F%E8%A1%8C%E9%A0%86)の問題が有るので注意してください。

## ラスタライズ時の半径を調整
今の仕組みだと、跳ねて水滴状になったパーティクルも同じサイズで描き続けている為、丸く固まった水が飛んでいく事になっています。実際の水は、ちぎれて徐々に小さくなることが多いので、密度か圧力辺りでラスタライズ時の半径を調整したらいい見た目にならないかなと考えています。

## シミュレーションを安定させるためには…（CFL条件について）
SPHの調整をやっていると、よくパーティクルがはじけて飛んでいく事があります。これの原因の一つとして、CFL条件( [Wikipedia](https://ja.wikipedia.org/wiki/CFL%E6%9D%A1%E4%BB%B6) )を満たしていない事が考えられます。
UE4のWaterプラグインによるシミュレーションでも同様ですが、数値解析系のシミュレーションが安定するには

```math
\Delta t < \frac{\Delta x}{v}
```
を満たす必要があるそうです。SPHの場合、$\Delta t$は差分時間、$\Delta x$は影響範囲の距離、$v$は粒子の速さを表しています。
その為、

1. サブステップを刻んで（Iterationsを増やして）$\Delta t$を小さくする
1. 影響範囲を広げて$\Delta x$を大きくする
1. 粒子の速さを下げる。（もしくは、上限を設ける）

などを検討する必要があります。ただ、1.と2.については負荷に直結しますし、3.についてはシミュレーション結果の質に影響しますので、一筋縄ではいかない問題ですが。

# 参考情報
### Niagara Simulation Framework Overview
https://www.youtube.com/watch?v=-cKgZrrBJ2w

途中でも紹介した、UnrealOpenDayの動画です。中国で開催されたので音声は中国語ですが、公式の英語字幕が付いていますので、日本語への自動翻訳でも問題ない精度で翻訳してくれています。
今回のSPHをNiagaraを使って実装する研究は、この動画を見つけたことから始まりました。基本的な部分はこの動画を参考に実装しています。

### NiagaraSimulationStagesプラグイン
https://github.com/EpicGames/~~~/Engine/Plugins/FX/NiagaraSimulationStage

githubのUE4-masterブランチにあります。手続きをすれば誰でも見られますが、一応privateのリポジトリなのでURLは部分省略してます。
リリースされている物には入っていないExperimental以前の代物ですが、色々と参考になる実装があります。特にモジュール群は使えるものが多いので、SimulationStageを使うなら一度確認する事をお勧めします。

https://twitter.com/shiba_zushi/status/1355890887323131908

このTweetでも紹介されている通り、色々とサンプルもありますが、かなり粗くNiagaraSystemを開くだけで処理負荷が高すぎてエディタが落ちたりしますので、とりあえずNiagaraのプレビューは閉じておく事をお勧めします。


### Niagaraチュートリアルシリーズ
https://www.zhihu.com/column/c_1323697556422369280

中国語のサイトですが、Google翻訳で概ね問題ないレベルの翻訳をしてくれます。
NiagaraのSimulationStageやそれ以外の事についてもかなり詳しく書かれています。

### UE4 Niagara Visual Debug
https://zhuanlan.zhihu.com/p/349472913

これも中国語のサイトです。結局私は使わなかったんですが、Niagaraのデバッグ手法として覚えておいても損は無さそうです。ただ、UE4.27でNiagaraのデバッグ機能が拡充されたので、そちらで十分かもしれません。

### バイラテラルフィルタの解説
https://imagingsolution.net/imaging/bilateralfilter/

このページでは、Range Domainが深度差ではなく輝度差を使っていますが、考え方は一緒です。数式や画像を使って分かりやすく説明してくれています。

# あとがき
引き続き研究中の内容ですが、ここまでに実装した内容をまとめてみました。
この辺を研究していると、テクニカルアーティストの領分だと痛感してます。理論や技術が分かってもそれを元に見た目やテクスチャを作る能力が無いと、良い見た目を作り出せずやきもきしています。

NiagaraのSimulationStageがまだExperimentalという事もあり、ゲームで使うにはまだ時期尚早な感じもありますが強力なツールであることは間違いありません。UE5.0で正式リリース（せめてベータ）になると嬉しいんですが、EpicGamesの頑張りに期待！

# 追記

https://www.unrealengine.com/ja/blog/unreal-engine-5-is-now-available-in-preview

この記事を準備している間に、UE5.0 Preview1がリリースされました。

https://twitter.com/WataruVFX/status/1496325646494142466

そして、この中には __NiagaraFluid__ というプラグインが入っていて、今回実装したような処理が簡単に作れるようになってます。
私も中身を見るのはこれからで、答え合わせをしていこうと思います！