解説: Drops

気が付いたら年が明けていました。あけましておめでとうございます。

さて、この記事では年末から制作を始めて先日公開した Drops の解説と制作背景を書いていこうと思います。どうしても SIMD が使いたかったので、ついに WebAssembly に手を出すことになりました。

そのおかげでシミュレーション部分が CPU で高速に動くようになったので、余った GPU パワーを使ってレンダリングパートにも力を入れてみました。

ソースコードはこちらから見られます。

You can find the source code from here!

全体の構成

全体の構成は以下のような感じになっています。

• シミュレーションパート
  • シンプルな SPH + 表面張力モデル
    • 4 回の substep で 1 フレーム
  • コア部分は WebAssembly で計算
    • JS からでは使えない SIMD で計算をベクトル化
    • JS 側でデータを詰めて WASM に渡す
    • 計算後にデータを JS 側に持ってくる
• メッシュ生成パート
  • Marching Cubes 法で得たメッシュを後処理で滑らかに
    • ここも WASM
• レンダリングパート
  • caustics のために Photon Mapping ぽいことをする
    • front face と back face をレンダリングして、ランダムな波長のフォトンをテクスチャに投影
    • 副産物として影ができる
  • 最後に caustics 用テクスチャを参照しながらシーンをレイトレ
    • ここでも front face と back face をレンダリングしておく

各パートに詳しい人なら大体どんな感じかこれで分かるかと思います。特にめちゃめちゃすごいことはやっておらず、一方で高速化はかなり頑張りました。

上から詳しく解説していきます。

シミュレーションパート

最近は Position Based Dynamics (PBD) を使った作品をよく作ってきたわけですが、今回は一周回ってシンプルな SPH 法を使いました。

位置ベース物理について考え直す

というのも、恐らく現在の state of the art な方法で PBD を使ってシミュレーションをしようとすると、XPBD によって物理的に正しい固さの拘束を細かい substep で計算していくことになると思います。

ただ、最近はここまで来ると位置ベースである必要って実はそんなにないのでは、みたいなことをよく考えていました。

Substep XPBD に至るまでの非常にざっくりとした流れ

元々の位置ベース物理は、物体の挙動を位置に対する拘束条件として記述し、速度を無視して位置に関する拘束条件だけを Gauss-Seidel 的に反復して解き、位置の変化量から速度を更新する陰的な手法でした。

しかし、本当に位置に対する条件しか見ていなかったため、反復の回数によって拘束の固さに違いが出るという問題が生じます。

CG 系の人々はその辺をあまり気にしない^[1]のでこの問題はしばらく放置されてきましたが、拘束の柔らかさをきちんと定式化して、速度とポテンシャルから次の時刻の物理的に正確な位置を反復計算で求めることができる XPBD という手法が登場します。

しかし、陰解法由来のエネルギー減衰問題については XPBD でも変わらず存在していました。その後、XPBD と小さな substep を組み合わせることで、高速かつ高精度な PBD の state of the art と呼べる Substep XPBD が開発されます。ここでは、もはや反復による位置更新処理は残されておらず、各 substep で位置更新のループを 1 回だけ回します。

拘束の計算も位置と速度と時間刻み幅を全て気にしながら行い、位置の更新自体も反復せずに陽解法的に 1 回だけ、ということになれば、位置ベースでないその他陽解法の手法との区別が曖昧になってきて、「元々の位置ベース物理の旨味ってなんだったっけ？」という疑問がわいてきます。

Substep XPBD でも残っている「位置ベース」由来の成分

そして、これに対する答えは恐らく「計算が安定することが保証されている」です。

「位置を正しい場所に戻す」ということを Gauss-Seidel 的に繰り返すので、普通に力を求めて速度を更新する方法と違い、どれほど大きな時間刻み幅で計算しても、シミュレーションが発散しないことが保証されています。

ということは、裏を返せば「爆発しないように調整された陽解法の手法を用いれば、PBD でなくても同じレベルで高速・高品質なシミュレーションが作れる」のでは？と思い至ります。

PBD の品質を普通の陽的手法で実現する

というわけで、計算手法としては一般的な陽解法に分類される SPH 法に、PBD における品質向上のエッセンスとも言える substepping を同じ方法で取り入れ、同レベルに高速で安定した計算を行うことに挑戦しました。

特にリアルタイム計算における SPH の難点として、

1. 圧縮しやすい
2. （1. を改善しようとして）圧力の係数を大きく設定すると爆発する
3. 強い外力に対して不安定になりがち

というものがありますが、これらは一般に陽解法に共通する問題でもあり、substepping の導入によって、陰解法における反復数のように「適当な分割数を決定する」ことで解消できます。2. については残念ながら完全な解決は難しいですが、そもそも爆発するほど大きな係数を設定しなくても十分非圧縮に振舞うようになるのであまり問題はありません。

さて、単純に時間刻み幅を N 分割するだけでは計算時間も N 倍になり、あまり効率的ではありません。

そこで、PBD と同じように、衝突や近傍粒子に関する情報は一つの step で一度だけ更新し、全ての substep 間で使い回します。衝突検出や近傍粒子探索は計算時間の中のかなりの割合を占める（陽解法なら尚更）ので、これで計算効率が大幅に向上します。

substep 内で多少変なことが起きても、結局ユーザーには見えないのでフレーム間を跨いで変な挙動が発生しなければ特に問題はありません。

が、さすがに全く同じように近傍粒子を計算するとsubstep 内で発生する新たな近傍の漏れが無視できない影響を及ぼすため、影響半径の 1.2 倍の粒子まで事前にペアを計算しておきます。その代わり、影響半径を通常の SPH よりもかなり小さい初期粒子間距離の 2.1 倍^[2]としています。

普通の SPH でここまで影響半径を小さくすると挙動が怪しくなってきますが、今は substepping によって実質的に SPH が陰解法化しているので、陰解法でないと厳しいレベルに challenging なパラメータを設定してもうまく動いてくれます。

粒子の個数と計算速度を考えた結果、substep 数を 4 としました。小さめの数値ですが、以前書いた記事で説明した通り、substepping は陰解法の iteration よりもはるかに強力なので、陰解法で反復を 10 回以上回したときの精度に匹敵する挙動が実現できているのではないかと思います。

今回用いた SPH 法の詳細

一口に SPH と言っても色々な派生があります。今回用いたのは Particle-based viscoelastic fluid simulation (PDF) で採用されている「距離に対して2乗で減少する density カーネル」と「距離に対して3乗で減少する near-density カーネル」を組み合わせたモデルです。ここからさらに独自に改変を加えた手法を用いています。

具体的には、以下のような手法で計算しています。

• 粒子の影響半径を r_e とします。また、ある 2 粒子 i, j 間の距離が r_{ij} であるとき、その 2 粒子間の重みを w_{ij} = 1 - r_{ij} と定義します。ただし、r_{ij}>r_e の場合は w_{ij}=0 とします。
• 粒子 i の近傍粒子とは、重み w_{ij} が正になるような粒子 j のことを指します。ただし、粒子 i は粒子 i の近傍粒子に含まれません。
• 粒子 i の粒子数密度 n_i を、粒子 i の全ての近傍粒子の重みの 2 乗の和 n_i = \sum_j w_{ij}^2 とします。
• 粒子 i の圧力 p_i を、標準粒子数密度 n^0 と圧力係数 k_p を用いて p_i=k_p(n_i-n^0) とします。
• 2 粒子 i, j 間の近傍がそれぞれの粒子に及ぼす圧力項は、2 粒子を結ぶ方向の単位ベクトルの w_{ij}(p_i+p_j+k_{p2}w_{ij}^3) 倍で求められます。ここで、k_{p2} は近-圧力係数です。

粘性項については省略しますが、エネルギーの減衰を最小限に抑えるため、中心力のみになるように調整しています（2 粒子を結ぶバネ-ダンパモデルを想像してください）。

論文中でも触れられていますが、実は、この計算だけである程度表面張力が再現できます。しかし、より大きな、あるいはより正確な表面張力を計算しようとすると、また別の工夫が必要になってきます。

表面張力モデル

そのための表面張力モデルは、Versatile Surface Tension and Adhesion for SPH Fluids (PDF) を参考に、というかほぼそのままのモデルを計算に取り入れています。このモデルの素晴らしい点として、曲率を用いた外力計算を行わないので、全体の運動量が保存されるというものがあります。

計算方法もとてもシンプルで、近傍間で重み付けした差分の足し合わせで粒子ごとの法線ベクトルを求め、それを正規化せずに再び近傍間で差分を取って表面積最小化項として適用する、というものです。

法線ベクトルを正規化しないので特異点・不連続点が発生せず、流体内部では対称性から自然に法線ベクトルが消失するので表面だけに力がかかります。さらに表面粒子を検出する必要もなく、pairwise な力なので運動量が保存されるといういいことずくめなモデルです。

難点として、陰的な拘束条件ではなく外力なので大きな表面張力を適用することが難しいことが挙げられますが、今の SPH は陰解法化していて最強なので、普通に遠慮なくデカい値をぶち込めば大きな表面張力が得られます。これはとても大きな魅力です。

例えば、substep なしではどう頑張ってもこのくらいの表面張力の強さしか獲得できません（これ以上強くしようとすると爆発します）。

しかし、substep を 6 まで増やすとここまで強力な表面張力を実現できます^[3]。

あるいは、もう少し力を小さくしてエネルギーがよりよく保存されるようなパラメータを選ぶこともできます。

ただし、いくら強くなったとはいえ根は陽解法なので、値がデカすぎると普通に爆発します。

Drops のページで Surface Tension のスライダーを最大まで上げると、substep が 4 回のときの、かなり限界に近い強さの表面張力を見ることができます。

今まで substepping ではなく陰解法での iteration を主に触ってきた感覚からすると、たった 4 回の substep でここまで強力な力を扱えるというのはかなり驚きがありました。

WebAssembly による高速化

この作品が可能になった最も大きな要因の一つに、WebAssembly によるシミュレーションの高速化が挙げられます。

と言っても、ただ単に WebAssembly を使えば速くなるわけではなく、SIMD による計算のベクトル化が必須です。自分の感覚では、これがないならわざわざ処理を WebAssembly に書き換えるメリットはあまりないように感じます。

というのも、実は JavaScript の実行エンジンの性能は本当に凄まじく、高度に最適化された JavaScript のコードはネイティブの機械語を吐く言語と比べて遜色ない速度で動作するからです。

しかし、SIMD 命令だけはどうあがいても JavaScript サイドからは実行できないので、これをうまく活用できるのであれば、既に極限まで最適化された JS のコードと比べても、WASM でうまく書き直すことでプログラムが数倍速く動くようになる、ということがあり得ます。

実際、今回のプログラムでは「もうこれ以上は無理」と思った JavaScript のコードと比べて 3 倍以上の高速化を達成することができました。

note that the JS part is already this optimized and I doubt there’s any room for further optimization.

but the WASM part is also extreme. I used AssemblyScript to write this but this is as good as writing assembly directly
maybe I should use Rust next time pic.twitter.com/odF8ifyt6I

— saharan / さはら (@shr_id) December 25, 2024

ただし、今回 WebAssembly を出力するための言語として AssemblyScript を選択してしまったというのもあって、WASM 側のコードも相当大変なことになっています。感覚的には直にアセンブリを書いているに等しいです。

AssemblyScript を触ってみて

触ってみて感じたのは次のような点です。

• TypeScript を踏襲しているので、C や C++ の struct に相当する概念が存在しない
  • オブジェクトを「作って」しまうと、それはヒープに確保され、オブジェクトへの参照はポインタに化けることになります
• 生 array が使えない
  • どうやっても「ポインタの配列」しか作れないので、memory.alloc で領域を確保して、直にデータを store & load することで解決
• ランタイムが充実しており、GC が組み込まれている
  • ほとんどの人にとっては嬉しい話ですが、CUDA のカーネルを書く、みたいな気持ちだと微妙かも
  • サポートされているランタイムは何種類かあり、GC どころか free すらできないものもある

普通に TypeScript を書くくらいの気持ちで使いたい場合にはとてもよさそうですが、アセンブリのレベルでゴリゴリに最適化されたコードを書きたい場合にはちょっと辛そうでした。

Rust でも WebAssembly のコードを出力できるらしいので、これを機に次回は Rust を使ってみようと思っています。

あるいは、WebAssembly の仕様自体がすごくシンプルなスタックマシンなので、何なら HGSL と同じように、Haxe のマクロ機能を使って言語の subset からコンパイル時に WASM のバイトコードを出力する、というのもある程度現実的かもしれません。最適化とか結構大変そうではありますが……

JS-WASM 間のデータの受け渡し

一部の処理を WASM に任せるにあたって割と気にしていたのがデータの受け渡しとその速度だったんですが、JavaScript と WebAssembly 双方からアクセス可能な線形メモリ (ArrayBuffer) が存在するので、そこから余計なオーバーヘッドなくデータを読み書きできます。

ただし、さすがにオブジェクトのプロパティを読み書きするようには書けないので、ポインタとオフセットから計算される位置に Float32Array などの typed array を用いてアクセスする必要があります。

メッシュ生成パート

メッシュは Marching Cubes 法を使って CPU 側で生成しています。これも結構時間のかかる処理なので、WASM 側でできる限りベクトル化を行いながら処理しています。

メッシュを生成する上では全てのセルを見る必要はなく、書き込まれた重みが閾値を超えたセルの周辺だけ見ることである程度の高速化になります。

生成されたメッシュの修正

MC 法で愚直に生成されたメッシュは、解像度が低かったりグリッドに入った重みの品質が低かったりすると結構なアーティファクトが出ます。

小さい液滴は仕方ないですが、大きな液滴にもグリッドの境界がハッキリ見えてしまっています。

そこで、生成したメッシュの頂点を平滑化して、法線ベクトルを再計算します。

だいぶマシになりましたが、屈折エフェクトを適用すると表面にグリッドに沿った凹凸が見えてしまいます。そのため、最後に法線ベクトルだけを追加で平滑化する処理を行います。

なぜ頂点ごと平滑化しないかというと、MC 法で生成されるメッシュは時系列方向に見たときにトポロジーが不連続に変化するので、移動したときに振動が生じてしまいます。

トポロジーに依存しにくい平滑化手法を使うことも考えらえますが、数万ポリゴンをリアルタイムで処理するには負荷が大きすぎるので、多少の振動については諦めることにしました。

最終的に法線を平滑化した後は以下のようなメッシュが得られます。

ここまで綺麗であれば問題なさそうですね！

余談: マーチングキューブ法で穴が開くのは解像度のせいではない

よくある誤解として「MC 法ではメッシュの張り方に曖昧さが生じるせいで穴が開くことがある。グリッドの解像度を上げれば解決する」というものがあります。

後半はまあ間違いではないのですが、この問題は解像度に本質があるのではなく、用いるテーブルの一貫性の欠如に本質があります。

一貫性のあるテーブルを用いてメッシュを作れば、解像度によらず絶対に穴が開かないようになります。興味がある方は以前書いた記事を読んでみてください。

レンダリングパート

シミュレーションとメッシュ構築が全部 CPU 側で動いているので、余った GPU パワーを惜しみなくレンダリングに使うことができます！

ということで以前から試してみたかった綺麗な caustics に挑戦してみます。

Caustics 計算

コースティクスというのは光の屈折・反射によって物体表面に届く光の量に濃淡が生じた結果生じる効果を指します。身近な例ではプールなどの水底に生じる光の網目模様などがあります。

綺麗な caustics で有名なデモに WebGL Water があります。

(Screenshot of "WebGL Water")

Caustics だけでなく波のシミュレーションや ambient occlusion も高速かつ高精度に計算されており、自分も見た当初大変な衝撃を受けたのを覚えています。

ここで用いられている caustics の計算手法は作者自身による解説記事が公開されており、グリッド状のメッシュを屈折させながら vertex shader で投影し、fragment shader で微分演算を用いて屈折による変形前後の面積比を計算することで光の強さを求めています^[4]。

リアルタイムの caustics 計算でまず最初に思い浮かんだのがこの手法だったので、とりあえずこの方法を真似てみることにしましたが、残念ながら水滴の境界付近での屈折が激しすぎたのかあまりいい見た目にはなりませんでした。

色々調節すれば動いた可能性はありますが、大変形した多量の三角形の描画には時間もかかるので、 Photon Mapping をベースにしたもう一つの手法を試してみることにしました。

フォトンマッピングはオフラインレンダリングの世界でも使われている手法で、きちんと計算すればとても物理的に正確な絵を出すことができます。しかし、リアルタイムで（それも WebGL という制約がかかった上で）フォトンを全て3次元空間上に格納して衝突判定を……というのは現実的ではないのでここでも色々な近似や制約を入れることになります。

この辺をリアルタイムでやろうとする研究はいくつか行われており、今回自分が最も参考にしたのは Caustics Mapping: An Image-space Technique for Real-time Caustics (PDF) という論文で、似たようなアプローチを取っています。

相違点や特筆すべき点は大体以下のような感じです。雰囲気で改造したので、既に同じ手法が出ているかもしれませんがあまり調べてません。

• モデルの頂点ではなく、出力先のテクセルに対して固定数のフォトン (vertex) を放った
  • 3 vertices per texel
  • 解像度 Medium (256×256) で 200k 近いフォトンが発射されるが、処理が軽いためか意外と耐えている
• 投影には vertex shader を使い、結果は加算合成を有効にすることで fragment shader 側で集計
  • Triangles ではなく Points で描画するとよい
  • MPM でのグリッド上の密度計算（これの下の方）でも使っていて、WebGL で和を計算したいときの常套手段として使える
    • ただし 2025 年 1 月現在、32-bit float ではアルファブレンドが iOS で使えないので注意！
• 影は個別に計算しておらず、フォトンが届かなかったところは暗くなる運用にしている
• 投影先が平面なので、明るさの計算はかなり適当
  • 一方で、綺麗な虹を出したかったので色味の計算は頑張った（後述）

虹を出す

特にリアルタイムだと caustics で虹が出ている（分光）ものは少なかったりします。Web だとざっと見た感じほぼ無かったので作ることにしました。

虹が出る原理

いろんなところでよく解説されているので詳細は省きますが、我々が普段見ている光というのはパッと見一色に見えてもいろいろな波長（≒色）の光が混ざり合っており、様々な媒質を屈折するときに異なる波長の光が異なる屈折率で屈折して別々の場所に届くことで虹が発生します。

これを再現するのは（方法としては）割と簡単で、異なる屈折率で異なる色の caustics を計算して、最後に加算合成をすればよいです。

この方法は先ほど触れた論文でも利用されており、図を引用すると次のような見た目のものが得られます。

(Image from Fig. 9 right of "Caustics Mapping: An Image-space Technique for Real-time Caustics")

いい感じに見えますね！　ただし、RGB それぞれで caustics を計算する必要があるので計算量は 3 倍になります。

もっと綺麗に虹を出す

一度これで実装をしてみて、それなりに納得がいく見た目にはなったのですが、画像を拡大してみると、意外と RGB がくっきり分かれています。

Caustics map にぼかし処理を入れているので、多少分光した程度では黄色や水色など中間色が滑らかに出現しますが、屈折が大きい部分では RGB それぞれの色がくっきりと分かれてしまいます。

あともう一つ、虹の七色目（日本では）である藍色が出現しません。

これはなぜかというと、人間が持つ「長波長（≒赤）の成分」を感じる細胞は、実は短波長（≒青）側にも多少の感度を持っており、このせいで RGB で各波長の色を表現しようとしたとき、長波長側と短波長側に 2 つの R 成分の山ができてしまいます。そのため単純に RGB で屈折率を分けたとすると、ぼかしを入れたとしても藍色の実現に必要な R の短波長側の山が無視されてしまうからなんですね。

この辺の話は、他にも「光の原色が三つなのは人間の都合」とか「ピンク色の単色光は存在しない」とか面白い話題があるので、興味がある方はぜひ深掘りしてみてください。

「虹色テクスチャ」を作る

これに対する対策は、RGB ごとに屈折率を複数割り当ててサンプリングするなどいくつか考えてみたんですが、最終的に波長に対応する RGB を記録したテクスチャを作り、波長を可視光の領域からランダムにサンプリングして対応する屈折率でフォトンを投影するという方法に落ち着きました。

波長に対する屈折率の変化は適当にやればそれっぽくなるんですが^[5]、波長から RGB を計算するのはそこそこ複雑な式が必要でテクスチャに焼き込んでしまった方が速かったです。

虹色テクスチャはこんな感じの 1 次元テクスチャです。

某所で「どうすっかな～」と思考を垂れ流していたところ、Mykhailo (@Michael_Moroz_) が色々助言をくれたりして最終的にこうなりました。　hey misha, thanks! :‌D

恐らく現実の液滴はこんなに分光を起こしませんが、めっちゃキレイですね。

液滴のレンダリング

実はメッシュができてしまえば難しくなくて、単純に表面の法線に従って視線のレイを屈折させ、skybox の色を取ってくるだけです。床にあたった場合は床の色を取ってきます。

ただし、液滴の裏面の法線も考慮したかったため、Efficient rendering of multiple refractions and reflections in natural objects (PDF) を参考に、裏面の法線も描画しておいて適当な割合（3～4 割程度）で混ぜることで対応しました。

フレネル反射も入れて、ついでに境界が見やすくなるようにリムライティングも入れています。この辺は他の処理に比べるとほぼノーコスト^[6]なので入れるだけ得です。

しかし、この方法では連続する屈折・反射を最大でも 2 回^[7]しか計算できないので、奥に隠れた液滴は存在しないことになってしまいます。多分これはきちんと計算しようとすると相当なコストを払わなければならない類のものなので、今回は諦めることにしました。

他にあったレンダリング方法の候補と却下された理由は以下の通りです。

密度マップを Ray Marching する

CPU 側でメッシュ生成に用いた密度マップがあるので、これを GPU 側に送ってレイトレ（レイマーチ）します。しかし、どう考えても全ピクセルを一定距離でレイマーチするのは処理が重すぎて caustics と両立しません。

SDF を作って Ray Marching する

次点でこれです。密度マップから renormalization の要領で SDF を作っておけば、後のレイマーチの負荷がだいぶ下がることが想定されます。

しかし、メッシュの解像度が 64x64x64 (=512×512) とそれなりに高いため、数回の renormalization を繰り返して SDF にするだけでも GPU にまあまあな負荷がかかることが予想されます。結局複数回の屈折と組み合わせると、割と厳しいか、ギリギリなのでは……？という試算になり、実装には踏み切れませんでした。

Stencil buffer を使って複数の表裏面の法線と深さを記録し、3回目以降の屈折は screen space で行う

やるとしたらこれです。で、これは多分それなりにうまく動いてくれます。

ただ、残念ながらここまで実装する体力と気力がありませんでした。。

それぞれの液滴に強めに屈折が働いている上に周囲を動き回るので、「見えるべきものが見えていない」ことが認知し辛く、正直今のままでもまあそこまで不自然には見えていないんじゃないかと思います、というのが言い訳です。多分みんな床の方に注目しているでしょうしね！

今回学んだこと

今回の作品ではいくつか新しい挑戦があり、それぞれの結果から学んだことは以下です。

Substepping は普通の陽解法でも強力

やっぱり強かったです。

これからは位置ベース物理に限らず適用していこうと思いました。

JS は速い、SIMD はもっと速い

JavaScript って最適化するとめちゃくちゃ速く動きます。

というのはまあ薄々気付いてはいた、というよりほぼ確信があったんですが、今回 WASM を書いてみて改めて JS（というより JS の実行エンジン）の底力を感じました。JS の限界と WASM の限界には思ったよりも差がありません^[8]。

そして、SIMD は上手くハマるとその限界を数倍レベルで突破できます。せっかく WebAssembly を使うならここを利用しない手はありません。

今後はボトルネックになりそうな処理で積極的に活用していきたいと思います。

Caustics はみんな好き

Caustics いいよね……