e.blog

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

TensorFlow Lite Unity Pluginを利用してDeepLabを動かすまで

DeepLearning AI DeepLearning Python Unity C# 機械学習

概要

最近、ARでSemantic Segmentaionを試そうと色々やっているのでそのメモです。
最終的には拾ってきたDeepLabモデルを変換して実際に動かすまでをやってみようと思います。

今回の記事はこちらを参考にさせていただきました。

asus4.hatenablog.com

元々TensorFlow LiteにはUnity Pluginがあるようで、上記ブログではそれを利用してUnityのサンプルを作ってるみたいです。
サンプル自体もTensorFlowにあるものを移植したもののようです。

公開されているUnityのサンプルプロジェクトは以下です。
今回の記事もこれをベースに色々調べてみたものをまとめたものです。

github.com

今回試したものは、以下のような感じで動作します。

youtu.be

全体を概観する

まずはTensorFlowについて詳しくない人もいると思うのでそのあたりから書きたいと思います。

TensorFlow Liteとは

いきなりサンプルの話にはいる前に、ざっくりとTensorFlow Liteについて触れておきます。

TensorFlow Liteとは、TensorFlowをモバイルやIoTデバイスなど比較的非力なマシンで動作させることを目的に作られたものです。
つまり裏を返せば、通常のTensorFlowはモバイルで動作させるには重くて向いていないということでもあります。

詳細についてはTensorFlow Liteのサイトをご覧ください。

www.tensorflow.org

こちらの記事も参考になります。

note.com

TensorFlow Liteのモデル

TensorFlow Liteのモデルは.tfliteという拡張子で提供され、主に、TensorFlow本体で作られた（訓練された）データを変換することで得ることができます。サイトから引用すると以下のように説明されています。

To use a model with TensorFlow Lite, you must convert a full TensorFlow model into the TensorFlow Lite format—you cannot create or train a model using TensorFlow Lite. So you must start with a regular TensorFlow model, and then convert the model.

特にこの「you must convert a full TensorFlow model into the TensorFlow Lite format.」というところからも分かるように、TensorFlow LiteではTensorFlowのすべての機能を使えるわけではないということです。
逆に言えば、軽量化・最適化を施して機能を削ることでエッジデバイス（スマホやIoTデバイス）でも高速に動作するというわけなんですね。

DeepLabとは

今回試したのは『DeepLab』と呼ばれる、Googleが提案した高速に動作するセマンティックセグメンテーションのモデル（ネットワーク）です。
セマンティックセグメンテーションについては後日書く予定ですが、こちらのPDFがだいぶ詳しく解説してくれているのでそちらを見るとより理解が深まるかと思います。

DeepLabについてはGoogleのブログでも記載があります。

ai.googleblog.com

ものすごくざっくり言うと、セマンティックセグメンテーションとは、画像の中にある物体を認識しその物体ごとに色分けして分類する、という手法です。
手法自体はいくつか提案されており、今回使用したDeepLabもそうした手法のうちのひとつです。

特徴としてはモバイルでも動作するほど軽く、高速に動くという点です。
そのためARで利用することを考えた場合に、採用する最有力候補になります。（なので今回試した）

ちなみにGitHubにも情報が上がっています。

github.com

ひとまず、前提知識として知っておくことは以上となります。
次からは実際にUnityで動かす過程を通して、最終的には拾ってきたDeepLabモデルを変換して実際に動かすまでをやってみようと思います。

Unityサンプルを動かす

Unityのサンプルを動かしてみましょう。
まずはGitHubからダウンロードしてきたものをそのまま実行してみます。

youtu.be

動画を見てもらうと自転車と車に対して反応しているのが分かるかと思います。

サンプルプロジェクトに含まれるモデルはこちらで配布されているものです。

スターターモデルと書いてあるように、若干精度は低めかなという印象です。
本来的にはここから自分の目的にあったものにさらに訓練していくのだと思います。

サンプル以外のモデルを試す

サンプルに含まれているモデルが動くことが確認できました。
次は別のモデルを試してみることにします。

モデルを変換する

さて、別のモデルといってもイチから訓練するものではなく、すでに訓練され公開されているものを探してきてそれを変換したものを使ってみたいと思います。
これによってサンプルとの違いができ、より深く動作を確認できると思います。

ということで参考にさせてもらったのが以下のリポジトリです。

github.com

ここのREADMEの下の方にあるリンクからモデルをダウンロードしてきます。

TensorFlow Liteのモデルへの変換について

TensorFlowで訓練したモデルをLite版へと変換していきます。

ちなみにTensorFlowのモデルにはいくつかの形式があり、以下の形式が用いられます。

  • SavedModel
  • Frozen Model
  • Session Bundle
  • Tensorflow Hub モジュール

これらのモデルについては以下のブログを参照ください。

note.com

今回はこの中の『Frozen Model』で保存されたものを扱います。

ちなみにもし自分で訓練データを作成する場合は注意しなければならない点があります。

前節で書いたように、TensorFlow Liteはいくつかの機能を削減して軽量化を図るものです。
そのため、TensorFlowでは問題なく動いていたものが動かないことも少なくありません。

このあたりについてはまだ全然詳しくないのでここでは解説しません。（というかできません）

以下の記事が変換についてとても詳しく書かれているので興味がある方は参照してみてください。
ちなみに以下の記事は「量子化（Quantization）」を目的とした変換に焦点を絞っています。

qiita.com

以下、余談。

量子化（Quantization）とは

上記記事では以下のように説明されています。

上記のリポジトリから Freeze_Graph (.pb) ファイルをダウンロードします。 ココでの注意点は ASPP などの特殊処理が入ったモデルは軒並み量子化に失敗しますので、なるべくシンプルな構造のモデルに限定して取り寄せることぐらいです。

※ ... ちなみにASPP - Atrous Spatial Pyramid Poolingの略です。

ASPPがなにかを調べてみたら以下の記事を見つけました。

37ma5ras.blogspot.com

記事には以下のように記載されています。

  1. 画面いっぱいに写っていようと画面の片隅に写っていようと猫は猫であるように，image segmentationではscale invarianceを考慮しなければならない． 著者はAtrous Convolutionのdilationを様々に設定することでこれに対処している(fig.2). 著者はこの技法を”atrous spatial pyramid pooling”(ASPP)と呼んでいる.

さらに該当記事から画像を引用させてもらうと、

おそらくこの画像の下の様子がピラミッドに見えることからこう呼んでいるのだと思います。
が、量子化に際してこの技法が含まれていると変換できないようなのでここではあまり深堀りしません。

簡単に量子化について触れておくと、以下の記事にはこう説明されています。

情報理論における量子化とは、アナログな量を離散的な値で近似的に表現することを指しますが、本稿における量子化は厳密に言うとちょっと意味が違い、十分な（=32bitもしくは16bit）精度で表現されていた量を、ずっと少ないビット数で表現することを言います。

ニューラルネットワークでは、入力値とパラメータから出力を計算するわけですが、それらは通常、32bitもしくは16bit精度の浮動小数点（の配列）で表現されます。この値を4bitや5bit、もっと極端な例では1bitで表現するのが量子化です。1bitで表現する場合は二値化（binarization）という表現がよく使われますが、これも一種の量子化です。

量子化には、計算の高速化や省メモリ化などのメリットがあります。

developer.smartnews.com

要するに、本来はintやfloatなど「大きな容量（16bit ～ 32bit）」を使うものをより小さい容量でも同じような精度を達成する方法、ということでしょうか。
これにはメモリ的なメリットや、単純に演算数の削減が見込めそうです。

モバイルではとにかくこのあたりの最適化は必須になるので量子化は必須と言ってもいいと思います。
（そういう意味で、量子化モデルの変換記事はめちゃめちゃ濃いので一度目を通しておくといいと思います）

が、後半で説明するtfliteへの変換でQuantizationのパラメータを設定するとうまく動作しなかったのでさらなる調査が必要そうです・・・。

閑話休題。

モデルの入力・出力を調べる

ではモデルを変換していきましょう。
モデルを変換するためにはいくつかの情報を得なければなりません。

ここでは詳細は割愛しますが、ディープラーニングではニューラルネットワークへの入力と出力（Input / Output）が大事になってきます。

ものすごくざっくり言えば、ディープラーニングはひとつの関数です。
プログラムでも、関数を利用したい場合はその引数と戻り値がなにかを知る必要があることに似ています。

自分でネットワークを構築し訓練したものであればすでに知っている情報かもしれませんが、だいたいの場合はどこからか落としてきたモデルや、あるいは誰かの構築したネットワークを利用するというケースがほとんどでしょう。

そのため入力と出力を調べる必要があります。
調べ方については上のほうで紹介したこちらの記事が有用です。

そこで言及されていることを引用させていただくと、

INPUTは Input、 形状と型は Float32 [?, 256, 256, 3]、 OUTPUTは ArgMax、 形状と型は Float32 [?, 256, 256] のようです。 なお一見すると ExpandDims が最終OUTPUTとして適切ではないか、と思われるかもしれませんが、 実は Semantic Segmentation のモデルにほぼ共通することですが ArgMax を選定すれば問題ありません。

とあります。
ここで言及されている名前はネットワーク次第なので毎回これになるとは限りません。
しかし

実は Semantic Segmentation のモデルにほぼ共通することですが ArgMax を選定すれば問題ありません。

というのはとても大事な点なので覚えておきましょう。

ちなみにモデルの中身を可視化するツールがあります。
以下の『Netron』というものです。

Webサービスを利用して見てもいいですし、アプリもあるのでよく使う場合はインストールしておいてもいいでしょう。

lutzroeder.github.io

では先ほど落としてきたモデルを読み込ませて見てみましょう。

入力はInputという名前で、入力の型はfloat32、入力のサイズは256 x 256の3チャンネルというのが分かります。
とても簡単に可視化できるのでとてもオススメのツールです。

続けて出力も見てみましょう。

参考にした記事に習ってArgMaxの部分を見てみます。
確認すると出力はArgMaxという名前でint32型、256 x 256の出力になるようです。

モデルの変換には「tflite_convert 」コマンドを使う

情報が揃ったのでダウンロードしてきたモデルを変換します。

変換にはtflite_convertコマンドを利用します。

今回の変換では以下のように引数を指定しました。

tflite_convert ^
  --output_file=converted_frozen_graph.tflite ^
  --graph_def_file=frozen_inference_graph.pb ^
  --input_arrays=Input^
  --output_arrays=ArgMax ^
  --input_shapes=1,256,256,3 ^
  --inference_type=FLOAT ^
  --mean_values=128 ^
  --std_dev_values=128

引数に指定している--input_arraysと--output_arraysが、先ほど調べた名前になっているのが分かります。
さらに--input_shapesには入力の形として先ほど調べた256 x 256を指定しています。
これを指定することで適切に変換することができます。

では変換されたモデルをNetronで可視化してみましょう。

内容が変化しているのが分かります。
これをUnityのサンプルプロジェクトに入れて利用してみます。

これであとは実行するだけ・・・にはいきません。

drive.google.com ※ 変換したモデルを念の為公開しておきます。

モデルに合わせてC#を編集する

実はサンプルで用意されているDeepLabスクリプトはサンプルに含まれているモデルに合わせて実装されているため、今回のケースの場合は少し修正をしなければなりません。

なのでNetronによって確認できる型や形状に定義を変更します。
ということでDeepLabクラスを修正します。

今回修正した内容のdiffは以下です。

さて、さっそく変換したモデルを読み込ませて使ってみましょう。
・・・と勢い込んでビルドしてみるものの動かず。Logcatで見てみると以下のようなエラーが表示されていました。

08-07 10:50:49.160: E/Unity(6127): Unable to find libc
08-07 10:50:49.163: E/Unity(6127): Following operations are not supported by GPU delegate:
08-07 10:50:49.163: E/Unity(6127): ARG_MAX: Operation is not supported.
08-07 10:50:49.163: E/Unity(6127): BATCH_TO_SPACE_ND: Operation is not supported.
08-07 10:50:49.163: E/Unity(6127): SPACE_TO_BATCH_ND: Operation is not supported.
08-07 10:50:49.163: E/Unity(6127): 53 operations will run on the GPU, and the remaining 31 operations will run on the CPU.
08-07 10:50:49.163: E/Unity(6127): TensorFlowLite.Interpreter:TfLiteInterpreterCreate(IntPtr, IntPtr)
08-07 10:50:49.163: E/Unity(6127): TensorFlowLite.Interpreter:.ctor(Byte[], InterpreterOptions)
08-07 10:50:49.163: E/Unity(6127): TensorFlowLite.BaseImagePredictor`1:.ctor(String, Boolean)
08-07 10:50:49.163: E/Unity(6127): TensorFlowLite.DeepLab:.ctor(String, ComputeShader)
08-07 10:50:49.163: E/Unity(6127): DeepLabSample:Start()

いくつかのオペレーションがGPUに対応していないためのエラーのようです。
ということで、以下の部分をfalseにして（GPU未使用にして）ビルドし直してみます。
（ちなみにiOSではGPUオンの状態でも問題なく動いたのでAndroid版の問題のようです）

public DeepLab(string modelPath, ComputeShader compute) : base(modelPath, true)
{
    // ... 後略

これを、

public DeepLab(string modelPath, ComputeShader compute) : base(modelPath, false)
{
    // ... 後略

こうします。
この第二引数がGPUを使うかどうかのフラグの指定になっています。

あとはこれをビルドし直して動かすだけです。

youtu.be

（すでに冒頭でも載せていますが）これで無事に動きました！

最後に

これがゴールではなくむしろスタート地点です。
ここから、独自のモデルの訓練をして目的に適合する結果を得られるように調整していかなければなりません。

が、ひとまずはTensorFlowで作られたモデルをtflite形式に変換して動作させるところまで確認できたので、あとはこの上に追加して作業をしていく形になります。

ここまで来るのは長かった・・・。
やっと本題に入れそうです。

ランキング参加中
Unity

Oculus Questでハンドトラッキングを使ってみる

C# VR Oculus Quest ハンドトラッキング

概要

Oculus Questのハンドトラッキングが利用できるようになったので使ってみたいと思います。
そこで、実際に使用するにあたってセットアップ方法とどういう情報が取れるのか、どう使えるかなどをまとめておきます。

ドキュメントは以下です。

developer.oculus.com

ちなみにドキュメントにも注意書きが書かれていますが、Oculus Linkでのハンドトラッキングは開発用でのみ動作するようです。

注:Oculus QuestとOculus Linkを使用する場合、PC上でのハンドトラッキングの使用はUnityエディターでサポートされています。この機能は、Oculus Quest開発者の反復時間短縮のため、Unityエディターでのみサポートされています。

Table of Contents

セットアップ

まずはプロジェクトをセットアップしていきます。

ここはHand Gesture用ではなく普通のOculusのセットアップです。
なのですでに知っている方は読み飛ばしてもらって大丈夫です。

XR Managementのインストール

Oculusを利用する場合は、Package ManagerからXR Managementをインストールします。

インストールすると、Project SettingsにXR Managementの項目が追加されるので、そこからOculus用のPluginをインストールします。

こちらもインストールが終わると画面が以下のように変化するので、Plugin ProvidersにOculus Loaderを追加します。

※ ちなみに、Unityエディタ上でOculus Linkを使って開発を行う場合はPC向けにもOculus Loaderを設定する必要があります。

Oculus Integrationをインポート

次に、Unity Asset StoreからOculus Integrationをインポートします。

インポートが終わったらOVRCameraRigをシーン内に配置します。
その際、シーン内のカメラと重複するので元々あったほうを削除します。

配置したらOVRCameraRigにアタッチされているOVR ManagerのHand Tracking SupportのリストからControllers and Handsを選択します。

※ ちなみにPlatformの設定がAndroidなっていないと設定できないようなので注意してください。

ハンドトラッキングを有効にする

ハンドトラッキングの機能を利用するためにはOculus Quest本体側の設定も必要になります。
ドキュメントから引用すると以下のように設定します。

ユーザーが仮想環境で手を使用するには、Oculus Quest上でハンドトラッキング機能を有効にする必要があります。

Oculus Questで、[Settings(設定)] > [Device(デバイス)]に移動します。 トグルボタンをスライドすることによってハンドトラッキング機能を有効にします。 手とコントローラーの使用を自動で切り替えられるようにするには、トグルボタンをスライドすることによって手またはコントローラーの自動有効化機能を有効にします。

シーンへの手の追加

ドキュメントによると、

手を入力デバイスとして使用するには、手動でシーンに追加する必要があります。手はOVRHandPrefab prefaにより実装されています。

とのことなので、対象のPrefabを配置します。
対象のPrefabはOculus/VR/Prefabsにあります。

それを、シーンに配置したOVRCameraRig以下のLeftHandAnchorとRightHandAnchorの下に配置します。

手のタイプを設定する

配置したOVRHandPrefabは両手用になっているので、以下の3つのコンポーネントの設定を適切な手のタイプ（左手 or 右手）に変更します。

  • OVRHand
  • OVRSkeleton
  • OVRMesh

以上でセットアップは終了です。
あとはOculus Linkでつないで再生ボタンを押すと以下のように手が表示されるようになります。

※ バージョンによってはこれで正常に動作しない場合があるかもしれません。その場合はOculus Integrationを最新にしてみてください。

ハンドトラッキングによるデータの取得

セットアップが終わったので、あとはハンドトラッキングシステムから得られるデータを用いて様々なコンテンツを作っていくことができます。
ここではいくつかのデータの取得方法をまとめておこうと思います。

各指のピンチ強度を測る

OVRHandには各指のピンチ強度(*)を測るAPIがあるので簡単に測ることができます。

  • ... ピンチ強度は各指が『親指とどれくらい近づいているかを測る値』です。曲がり具合ではないので注意です。
float thumbStr = _rightHand.GetFingerPinchStrength(OVRHand.HandFinger.Thumb);
float indexStr = _rightHand.GetFingerPinchStrength(OVRHand.HandFinger.Index);
float middleStr = _rightHand.GetFingerPinchStrength(OVRHand.HandFinger.Middle);
float ringStr = _rightHand.GetFingerPinchStrength(OVRHand.HandFinger.Ring);
float pinkyStr = _rightHand.GetFingerPinchStrength(OVRHand.HandFinger.Pinky);

このあたりのデータを組み合わせれば、簡単なジェスチャーなどは検知できそうですね。

ボーン情報を取得する

ボーンの各情報を得るためにはOVRSkeletonクラスを利用します。
OVRSkeletonにはOVRBone構造体を保持するリストがあり、そこから情報を取り出します。

リストのどこにどのボーン情報が入っているかはOVRSkeleton.BoneId enumによって定義されており、それをintに変換して利用します。

なお、どこにどの情報が入っているかはドキュメントに記載されています。引用すると以下のように定義されています。

Invalid          = -1
Hand_Start       = 0
Hand_WristRoot   = Hand_Start + 0 // root frame of the hand, where the wrist is located
Hand_ForearmStub = Hand_Start + 1 // frame for user's forearm
Hand_Thumb0      = Hand_Start + 2 // thumb trapezium bone
Hand_Thumb1      = Hand_Start + 3 // thumb metacarpal bone
Hand_Thumb2      = Hand_Start + 4 // thumb proximal phalange bone
Hand_Thumb3      = Hand_Start + 5 // thumb distal phalange bone
Hand_Index1      = Hand_Start + 6 // index proximal phalange bone
Hand_Index2      = Hand_Start + 7 // index intermediate phalange bone
Hand_Index3      = Hand_Start + 8 // index distal phalange bone
Hand_Middle1     = Hand_Start + 9 // middle proximal phalange bone
Hand_Middle2     = Hand_Start + 10 // middle intermediate phalange bone
Hand_Middle3     = Hand_Start + 11 // middle distal phalange bone
Hand_Ring1       = Hand_Start + 12 // ring proximal phalange bone
Hand_Ring2       = Hand_Start + 13 // ring intermediate phalange bone
Hand_Ring3       = Hand_Start + 14 // ring distal phalange bone
Hand_Pinky0      = Hand_Start + 15 // pinky metacarpal bone
Hand_Pinky1      = Hand_Start + 16 // pinky proximal phalange bone
Hand_Pinky2      = Hand_Start + 17 // pinky intermediate phalange bone
Hand_Pinky3      = Hand_Start + 18 // pinky distal phalange bone
Hand_MaxSkinnable= Hand_Start + 19
// Bone tips are position only. They are not used for skinning but are useful for hit-testing.
// NOTE: Hand_ThumbTip == Hand_MaxSkinnable since the extended tips need to be contiguous
Hand_ThumbTip    = Hand_Start + Hand_MaxSkinnable + 0 // tip of the thumb
Hand_IndexTip    = Hand_Start + Hand_MaxSkinnable + 1 // tip of the index finger
Hand_MiddleTip   = Hand_Start + Hand_MaxSkinnable + 2 // tip of the middle finger
Hand_RingTip     = Hand_Start + Hand_MaxSkinnable + 3 // tip of the ring finger
Hand_PinkyTip    = Hand_Start + Hand_MaxSkinnable + 4 // tip of the pinky
Hand_End         = Hand_Start + Hand_MaxSkinnable + 5
Max              = Hand_End + 0

なお、以下の記事から画像を引用させていただくと、各IDの割り振りはこんな感じになるようです。

Finger's ID

qiita.com

簡単なハンドコントロール

OVRHandは、ホーム画面などで利用される手と同じ機能を簡単に利用するためのAPIを提供してくれています。
その情報にアクセスするにはOVRHand.PointerPoseプロパティを利用します。

これはTransform型で、手をポインタとして見た場合の位置と回転を提供してくれます。

視覚化してみたのが以下の動画です。

youtu.be

個人的にはやや直感に反する挙動だなと思っています。
手の指の向きは参考にされていないようで、基本的にポインタ方向は『手の高さ』によって算出されているような印象を受けます。

手のひらの法線を計算する

今回、個人プロジェクトで『手のひらの法線』が必要になり、それを求めるプログラムを書いたので参考までに載せておきます。

using System.Collections;
using System.Collections.Generic;

using UnityEngine;

namespace Conekton.ARUtility.Input.Application
{
    public class HandPoseController : MonoBehaviour
    {
        [SerializeField] private OVRHand _targetHand = null;
        [SerializeField] private OVRSkeleton _rightSkeleton = null;
        [SerializeField] private Transform _palmNormalTrans = null;
        [SerializeField] private float _detectLimit = 0.5f;

        private OVRSkeleton.BoneId[] _forPalmCalcTargetList = new[]
        {
            // First two of them are used for calculating palm normal.
            OVRSkeleton.BoneId.Hand_Index1,
            OVRSkeleton.BoneId.Hand_Pinky0,

            OVRSkeleton.BoneId.Hand_Middle1,
            OVRSkeleton.BoneId.Hand_Ring1,
            OVRSkeleton.BoneId.Hand_Pinky1,
            OVRSkeleton.BoneId.Hand_Thumb0,
        };

        private OVRSkeleton.BoneId BoneIDForNormalCalculation1 => OVRSkeleton.BoneId.Hand_Index1;
        private OVRSkeleton.BoneId BoneIDForNormalCalculation2 => OVRSkeleton.BoneId.Hand_Pinky0;

        public bool TryGetPositionAndNormal(out Vector3 position, out Vector3 normal)
        {
            if (!_targetHand.IsTracked)
            {
                position = Vector3.zero;
                normal = Vector3.up;
                return false;
            }

            Vector3 center = Vector3.zero;

            foreach (var id in _forPalmCalcTargetList)
            {
                OVRBone bone = GetBoneById(id);
                center += bone.Transform.position;
            }

            center /= _forPalmCalcTargetList.Length;

            position = center;

            OVRBone bone0 = GetBoneById(BoneIDForNormalCalculation1);
            OVRBone bone1 = GetBoneById(BoneIDForNormalCalculation2);

            Vector3 edge0 = bone0.Transform.position - center;
            Vector3 edge1 = bone1.Transform.position - center;

            normal = Vector3.Cross(edge0, edge1).normalized;

            return true;
        }

        private OVRBone GetBoneById(OVRSkeleton.BoneId id)
        {
            return _rightSkeleton.Bones[(int)id];
        }
    }
}

考え方はシンプルで、ハンドトラッキングから得られるボーンの位置をいくつか選び、それらの平均の位置を手のひらの位置としています。
また法線については、求めた手のひらの位置とふたつのボーンの位置との差分ベクトルを取り、それの外積を取ることで求めています。

まとめ

Oculus Questのハンドトラッキングの精度は驚異と言っていいと思います。

過去に、Leapmotionを使ったコンテンツを開発したことがありますが、Leapmotionよりも精度が高い印象です。（それ用デバイスより精度が高いって・・）
実際に体験すると、本当にVR内に自分の手があるかのように思えるくらいなめらかにトラッキングしてくれます。

ハンドトラッキングを用いたUI/UXはさらに発展していくと思うので、これからとても楽しみです。

Unityの推論エンジン『Barracuda』を試してみたのでそのメモ

C# Unity 機械学習 DeepLearning

概要

以下の記事を参考に、最近リリースされたUnity製推論エンジンを試してみたのでそのメモです。

qiita.com

note.com

Barracudaとは？

Barracudaとは、ドキュメントにはこう記載されています。

Barracuda is lightweight and cross-platform Neural Net inference library. Barracuda supports inference both on GPU and CPU.

軽量でクロスプラットフォームなニューラルネットワークの推論ライブラリということですね。
そしてこちらのブログによるとUnity製のオリジナルだそうです。

セットアップ

BarracudaはPackage Managerから簡単にインストールできます。
インストールするにはWindow > Package ManagerからBarracudaを選択してインストールするだけです。

モデルの準備

これでC#からBarracudaを利用する準備は終わりです。
しかし、Deep Learningを利用して処理を行うためには訓練済みのモデルデータが必要です。
これがなければDeep Learningはなにも仕事をしてくれません。

Barracudaで扱えるモデル

Deep Learningを利用するためのフレームワークは多数出ています。有名どころで言えばTensorFlowなどですね。
そしてこうしたフレームワークごとにフォーマットがあり、そのフレームワークで訓練したデータは独自の形式で保存されます。
つまり言い換えればフレームワークごとにデータフォーマットが異なるということです。

そしてBarracudaでもそれ用のデータ・フォーマットで保存されたモデルデータが必要になります。
しかしBarracudaではONNX形式のモデルデータも扱うことができるようになっています。

ONNXとは？

ONNXはOpen Neural Network eXchangeの略です。（ちなみに『オニキス』と読むらしいです）
Openの名がつく通り、Deep Learningのモデルを様々なフレームワーク間で交換するためのフォーマット、ということのようです。

前述の通りBarracudaでも利用できます。

ONNXモデルの配布サイト

以下のリポジトリからいくつかの学習済みモデルがDownloadできます。

github.com

モデルの変換

Barracudaでは主要なフレームワークのモデルからBarracuda形式およびONNX形式に変換するためのツールを提供してくれています。　

note.com

詳細は上記の記事を見てもらいたいと思いますが、どういう感じで変換を行うのかのコマンド例を載せておきます。

$ python tensorflow_to_barracuda.py ../mobilenet_v2_1.4_224_frozen.pb ../mobilenet_v2.nn

Converting ../mobilenet_v2_1.4_224_frozen.pb to ../mobilenet_v2.nn
Sorting model, may take a while...... Done!
IN: 'input': [-1, -1, -1, -1] => 'MobilenetV2/Conv/BatchNorm/FusedBatchNorm'
OUT: 'MobilenetV2/Predictions/Reshape_1'
DONE: wrote ../mobilenet_v2.nn file.

モデルを利用して推論する（Style Chnage）

冒頭で紹介したこちらの記事からStyle Changeの方法を見ていきます。
（これがおそらく一番短くて分かりやすい例だと思います）

using UnityEngine;
using Unity.Barracuda;

public class StyleChange : MonoBehaviour
{
    [SerializeField] private NNModel _modelAsset = null;
    [SerializeField] private RenderTexture _inputTexture = null;
    [SerializeField] private RenderTexture _outputTeture = null;

    private Model _runtimeModel = null;
    private IWorker _worker = null;

    private void Start()
    {
        // Load an ONNX model.
        _runtimeModel = ModelLoader.Load(_modelAsset);

        // Create a worker.
        // WorkerFactory.Type means which CPU or GPU prefer to use.
        _worker = WorkerFactory.CreateWorker(WorkerFactory.Type.Compute, _runtimeModel);
    }

    private void Update()
    {
        Tensor input = new Tensor(_inputTexture);
        Inference(input);
        input.Dispose();
    }

    private void Inference(Tensor input)
    {
        _worker.Execute(input);
        Tensor output = _worker.PeekOutput();
        output.ToRenderTexture(_outputTeture, 0, 0, 1/255f, 0, null);
        output.Dispose();
    }

    private void OnDestroy()
    {
        _worker?.Dispose();
    }
}

だいぶ短いコードですね。これで映像を変換できるのだから驚きです。
ただし、複雑なネットワークを使っている場合はその分処理が重くなるので、リアルタイムなポストエフェクトとしては利用できないと思います。

このコードを実行すると以下のような結果が得られます。
（もちろん、設定するモデルによって出力の絵は変わります）

InputにRenderTextureを与え、OutputもRenderTextureで受け取っていますね。
入出力ともに画像なので利用イメージがしやすいと思います。

しかし（書いておいてなんですが）こうしたスタイルの変更というのはゲームではあまり使用されないかもしれません。
それよりも、物体検知やセグメンテーションなどでその真価を発揮するのではないかなと思っています。

ということで、次は物体検知についても書いておきます。

モデルを利用して推論する（物体検知）

以下のコードはこちらのGitHubのものを参考にさせていただいています。
（ファイルへの直リンク）

using System;
using Barracuda;
using System.Linq;
using UnityEngine;
using System.Collections;
using System.Collections.Generic;
using System.Text.RegularExpressions;

public class Classifier : MonoBehaviour
{
    public NNModel modelFile;
    public TextAsset labelsFile;

    public const int IMAGE_SIZE = 224;
    private const int IMAGE_MEAN = 127;
    private const float IMAGE_STD = 127.5f;
    private const string INPUT_NAME = "input";
    private const string OUTPUT_NAME = "MobilenetV2/Predictions/Reshape_1";

    private IWorker worker;
    private string[] labels;


    public void Start()
    {
        this.labels = Regex.Split(this.labelsFile.text, "\n|\r|\r\n")
            .Where(s => !String.IsNullOrEmpty(s)).ToArray();
        var model = ModelLoader.Load(this.modelFile);
        this.worker = WorkerFactory.CreateWorker(WorkerFactory.Type.ComputePrecompiled, model);
    }


    private int i = 0;
    public IEnumerator Classify(Color32[] picture, System.Action<List<KeyValuePair<string, float>>> callback)
    {
        var map = new List<KeyValuePair<string, float>>();

        using (var tensor = TransformInput(picture, IMAGE_SIZE, IMAGE_SIZE))
        {
            var inputs = new Dictionary<string, Tensor>();
            inputs.Add(INPUT_NAME, tensor);  
            var enumerator = this.worker.ExecuteAsync(inputs);

            while (enumerator.MoveNext())
            {
                i++;
                if (i >= 20)
                {
                    i = 0;
                    yield return null;
                }
            };

            // this.worker.Execute(inputs);
            // Execute() scheduled async job on GPU, waiting till completion
            // yield return new WaitForSeconds(0.5f);

            var output = worker.PeekOutput(OUTPUT_NAME);

            for (int i = 0; i < labels.Length; i++)
            {
                map.Add(new KeyValuePair<string, float>(labels[i], output[i] * 100));
            }
        }

        callback(map.OrderByDescending(x => x.Value).ToList());
    }


    public static Tensor TransformInput(Color32[] pic, int width, int height)
    {
        float[] floatValues = new float[width * height * 3];

        for (int i = 0; i < pic.Length; ++i)
        {
            var color = pic[i];

            floatValues[i * 3 + 0] = (color.r - IMAGE_MEAN) / IMAGE_STD;
            floatValues[i * 3 + 1] = (color.g - IMAGE_MEAN) / IMAGE_STD;
            floatValues[i * 3 + 2] = (color.b - IMAGE_MEAN) / IMAGE_STD;
        }

        return new Tensor(1, height, width, 3, floatValues);
    }
}

ちなみにStartメソッドで読み込んでいるテキストは分類の名称が改行区切りで入っているただのテキストファイルです。
以下はその一部。（ファイルへの直リンク）

background
tench, Tinca tinca
goldfish, Carassius auratus
great white shark, white shark, man-eater, man-eating shark, Carcharodon carcharias
tiger shark, Galeocerdo cuvieri
hammerhead, hammerhead shark
electric ray, crampfish, numbfish, torpedo
stingray
cock
hen
...

入力画像から識別したラベルに対する確率を取得する

これらの動作の基本的な流れは以下です。

  1. 画像をテンソル化して入力データとする
  2. ニューラルネットワークを通して出力を得る
  3. 出力はラベル数と同じ数の1階テンソルで、値はそれぞれの分類の確度（%）
  4. あとは確度に応じて望む処理をする

という具合です。
そしてその出力部分を抜粋すると以下のようになっています。

var output = worker.PeekOutput(OUTPUT_NAME);

for (int i = 0; i < labels.Length; i++)
{
    map.Add(new KeyValuePair<string, float>(labels[i], output[i] * 100));
}

output[i] * 100としている箇所が確率を%に変換している部分ですね。
つまり、該当番号（i）のラベルである確率をDictionary型の変数に入れて返しているというわけです。

取得するテンソルはvar output = worker.PeekOutput(OUTPUT_NAME);でアクセスしています。
そしてこのOUTPUT_NAMEはprivate const string OUTPUT_NAME = "MobilenetV2/Predictions/Reshape_1";と定義されています。

この文字列は訓練されたモデルのネットワークの変数名でしょう。
つまり推論の結果をこれで取得している、というわけですね。

モデルのInput / Outputを確認する

前述したように、モデルへのInputとOutputを明示的に指定する必要があり、そのための文字列を知る必要があります。

そのための情報はインスペクタから確認することができます。
インポートしたモデルファイルを選択すると以下のような情報が表示されます。

この図のInputs (1)とOutputs (1)がそれに当たります。
その下のLayersはニューラルネットワークのレイヤーの情報です。
つまりどんなネットワークなのか、ということがここで見れるわけですね。

画像として出力されるモデルかどうか確認する

ちなみにStyle Changeのところで書いたように、画像自体を変換して出力するネットワークなのかどうかは以下のInputs / Outputsを確認すると分かると思います。

入力が画像サイズで、出力も画像サイズになっている場合はそれが画像として出力されていると見ることが出来ます。
出力が[1, 224, 224, 3]となっているのは224 x 224サイズの3チャンネル（RGB）を1つ出力することを意味しているわけですね。

まとめ

Deep Learning（ニューラルネットワーク）は、極論で言えば巨大な関数であると言えると思います。
入力となる引数も巨大であり、そこから目的の出力を得ること、というわけですね。

y = f(x)

の、 xが入力（つまり今回の場合は画像）で、その結果（どのラベルの確率が高いか）が y、というわけです。

これを上記の物体検知に当てはめてC#で書き直すと、

// inferenceは「推論」を意味する英単語
Dictionaly<string, float> map = Inference(_inputTexture);

という感じですね。

最近はDeep Learningについてずっと調査をしています。

そこでの自分の大まかな理解を書いておくと、この巨大な関数を解くためのパラメータを『機械学習』で調整させる。
その調整する方法が『パーセプトロン』をベースとする『ニューラルネットワーク』を利用して行っている。

そしてこの調整されたパラメータとネットワーク構成がつまりは訓練済みモデル（データ）というわけです。
なのでそのパラメータを利用した関数を通すとなにかしらの意味がある出力が得られる、というわけなんですね。

もちろん、そのネットワークをどう組むか。それがどう実現しているのか。基礎を学ぼうとすると膨大な知識が必要になります。

しかし、こと利用する視点だけで見ればなんのことはない、ただの関数実行だ、と見ると利用がしやすいのではないかと思います。

これを機に色々とUnity上でディープラーニングを使って色々とやっていきたいと思います。

ランキング参加中
Unity

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

ComputeShader C# Unity シェーダ レンダリング パーティクル アニメーション

概要

今回は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

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

このグループクラスでは、子どもに持つ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のコードではグループ化の最適化を行っています。具体的には、サブグループという概念を取り入れてグループデータはひとつにまとめ、サブグループ単位でパーティクル位置を指定できるようにしました。

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

ランキング参加中
Unity

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

C# Unity UniRx

概要

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

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

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

qiita.com

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

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

ドラッグ処理

最近実装したものでドラッグ処理です。
例で使われている_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で実装する

C# レンダリング 数学 シェーダ Unity

はじめに

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

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 API（Quaternion.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.CalculateFrustumPlanesとGeometryUtility.TestPlanesAABBメソッドです。

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

ズームも行える

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

www.youtube.com

最後に

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

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

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

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

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

C# Unity レンダリング 数学 パーティクル シェーダ ComputeShader アニメーション

はじめに

以前に、格子法を用いた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}

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

余談

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

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

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

そして全微分の係数（a）は偏微分で求まります。（同様に、多数変数の場合はb、cと変数分、係数が増え、それぞれの変数の偏微分で求まります）
（例えば 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っていうことですね。

ちなみに記事内で言及されている以下の点、

例えば、３変数関数 f(x, y, z) = x^2 + y^3 + yz を (x_0,y_0,z_0) の近くで一次近似すると、
f(x, y, z) ≒ f(x_0, y_0, z_0) + 2x_0(x − x_0) + (3y^2_0 + z_0)(y − y_0) + y_0(z − z_0) となります。
（なぜ上式のようになるのか、計算は後ほど「全微分と偏微分の関係」で述べます）

x − x_0 = \Delta x などと置いて上式を移項すると、 f(x_0 + \Delta x, y_0 + \Delta y, z_0 + \Delta z) − f(x_0, y_0, z_0) ≒ 2x_0 \Delta x + (3y^2_0 + z_0)\Delta y + y_0 \Delta z となります。

これはこう解釈できます。

f(x_0, y_0, z_0)の近くで一次近似ということは、その面に対しての接線ということです。
つまり、その地点においては線とみなせる＝微小量変化した、というわけですね。

そして「xをどれくらい変化させたら関数にどういう変化をもたらすか」というのは偏微分で求まる微分係数です。

偏微分を用いて変化率を求めると \frac{\partial f}{\partial x} = 2xです。
そして増加量は \Delta xなので 2x_0 \Delta xになりますね。

そしてそれぞれの変数分同じことをやれば上記式になります。

閑話休題。

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

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

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

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

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

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

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


\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

https://www.epii.jp/articles/note/physics/continuum/derivativewww.epii.jp