読者です 読者をやめる 読者になる 読者になる

alembic for realtime rendering


ここ数ヶ月仕事で Alembic を扱ってきたので、Alembic に関してここまで得られた知見を書き残しておこうと思います。
想定しているシチュエーションは主に、Alembic のライブラリを用いて alembic ファイル (.abc) を読み込み、そのデータを D3D11 などのグラフィックス API に流し込んで描画する、というようなものです。
Alembic はゲームに使うのは難しいと思われますが、例えば VR コンテンツには使い出がありますし、近年のゲームエンジンを使った映像制作でも重要な役割を担う可能性があります。また、Alembic は書き込みも簡単にできるため、プレイ中のゲームのシーンを Alembic で 3D 録画し、それを映像制作などに利用、といった応用も考えられます。

Alembic 概要

http://www.alembic.io/

過去にも軽く触れましたが、Alembic とは主に映像業界で使われているデータフォーマットです。ポリゴンメッシュなどを全フレームベイクして格納するのに用いられます。映像業界ではスキニングや物理シムは全てベイクして Alembic に格納し、レンダラやコンポジットソフトに渡す、といった使い方がなされます。
全フレームベイクして格納するわけですから、当然ながらゲーム屋からするとデータ量は桁外れにでかいものとなります。

Alembic はいろんなデータを時系列に沿って格納できるようになっています。データの各要素は property と呼ばれ、格納されているデータ群は sample と呼ばれます。ポリゴンモデルなどの具体的なオブジェクトは schema と呼ばれるオブジェクトに格納されています。この schema はいろんな property の集合体になっています。
例えばポリゴンモデルは PolyMesh という schema に格納されており、この PolyMesh は position、uv、index などの各種 property で構成されています。そして、時間 (またはインデックス) を指定すると、Alembic のライブラリはその時間/インデックスに対応する sample (PolyMesh の場合モデルデータ) を返す、といった具合に機能します。ちなみに Alembic は補間は行いません。

Alembic 導入

いきなりですが、自分的にこの導入が Alembic を使うにあたっての最大の難関だと思っています。
Alembic はバイナリの配布がされておらず、使うには自力でビルドする必要があります (関連する HDF や ILMBase も含めて)。
ビルドシステムが CMake なので、CMake で VisualStudio のプロジェクトを生成してビルドするだけ…、といきたいところですが、残念ながらそう簡単にはいきません。映像業界では Linux がデファクトスタンダードだそうで、映像関係のオープンソースソフトウェアは Windows はまともにサポートされてないことが多いです。Alembic も素直にはビルドできず、CMakeList とソース両方にいくつか手を加えてビルドを通るようにする必要がありました。
その詳細は面倒なので省略しますが、ここらへんにヘッダファイルと VisualStudio2015 用のビルド済みバイナリを置いています (容量削減のため .lib ファイル郡は .7z に圧縮済み)。とりあえず手を付けてみたいという方はこれの使用をおすすめします。

入門用の簡単なサンプルプログラムが欲しいのであれば、Alembic のソースに含まれるテストコードが参考になるでしょう。(これ とか同ディレクトリにある他のテストプログラム) 意外と簡単そうだというのが見て取れると思います。
より実践的な例が欲しいのであれば、拙作 AlembicImporter のソースが参考になるかもしれません。下の Schema の項目で触れる問題点に一通り対処してあります。Unity のプラグインとして使う前提で作ってはいますが、このプラグイン自身は Unity には非依存です。

Alembic の archive には HDF と Ogawa の 2 種類がありますが、HDF はレガシーなフォーマットです。自分で .abc を生成する場合は Ogawa を選んだほうがいいでしょう。こちらの方が速くて小さいことが多く、最近の DCC ツールでも大体 Ogawa がデフォルトになっています。
また、Alembic はエラーは例外を投げて通知してくることが多いです。ゲーム畑の人は面食らうかもしれません。私は面食らいました。

Schema

先に触れたように、ポリゴンモデルなどの具体的なデータは schema というオブジェクトが保持します。ここではいくつかの schema について使用の際の注意点などを挙げます。Xform, Camera, Points, PolyMesh の 4 種です。
schema は他にも何種類かありますが、Light はそれ自身にはライトに関連するパラメータは入っていないし、Material は抽象的すぎて対応が困難な上 DCC ツール側もろくに対応していない有様で、それ以外は Subdivision や NURBS といったもので今回は保留です。Material に関してはインポートの後独自対応が必要になりますが、現状映像業界でもそんな感じなようです。

  • Xform

Xform とは位置、回転、スケール などを表すもので、ゲームでは同等品は Transform と呼ばれることが多いと思います。Alembic では 1 つのオブジェクトが複数の schema を持つことはできないため、ほとんどのケースでは Xform の子として他の schema がぶら下がる形になっています。
例:

root - Xform - Camera
             |- Xform - PolyMesh

この Xform、Transform 相当品と書くと何ら難しいことはなさそうに見えますが、内部処理を理解していないとハマる可能性があります。
Xform は座標やスケール値などを直接保持するのではなく、各オペレーションおよびその順番を保持しています。オペレーションは 移動、任意軸回転、X/Y/Z 軸回転、スケールなどの種類があり、これらを任意の数、任意の順番で保持できます。そして、XformSample::getMatrix() を呼ぶとこのオペレーションを全て実行して結果を行列で返します。いろんな DCC ツールに対応するための苦慮が伺える実装です…。

問題は回転を取得したい場合です。XformSample には getAxis(), getAngle() というのが用意されているのでそれで簡単に取れるように見えますが、罠があります。getAxis() getAngle() は、getMatrix() で全オペレーションを実行して行列を算出した後、その行列からクオータニオンを抽出し、axis/angle に変換して返します。これは遅いだけでなく、回転とスケールが同時にかかっている場合に元と違う結果になってしまいます。正確な回転を取得したい場合、自力で全オペレーションを巡回して回転オペレーションだけを見て算出しなければなりません。
XformSample::getScale() も全く同じ問題を含んでおり、これも getMatrix() の結果からスケール値を抽出しようとするので、正確な値を取得したければ自力で全オペレーションを巡回して算出する必要があります。

  • Camera (ICamera, OCamera)

カメラの情報が入った schema です。含まれる情報は実にプリレンダ向けといった趣きで、焦点距離、撮影距離、口径、シャッター開始時間/終了時間 などなど、そういう感じです。
幸い、CameraSample には getNearClippingPlane(), getFarClippingPlane(), getFieldOfView() といったリアルタイムでも馴染み深いのもちゃんとあります。ただし、getFieldOfView() は横方向の FoV を返してきます。リアルタイムの場合必要なのは縦方向のそれです。残念ながらこれは自力で算出する必要があります。算出の処理はこうなります。

    verticalFoV: 2.0 * RadiansToDegrees( atan( verticalAperture * 10.0 / (2.0 * focalLength)));

ちなみにこれは getFieldOfView() の実装をコピペ改変しただけです。CameraSample は FoV そのものは保持しておらず、aperture と focal length から算出するようになっています。このため、逆に Alembic へカメラ情報をエクスポートしたい場合、FoV から aperture か focal length を算出する必要があります。AlembicImporter では aperture を固定して focal length を算出するようにしています。focal length 算出は以下のような処理になります。

    focalLength:  (aperture * 10.0) / tan(DegreesToRadians(fieldOfView / 2.0)) / 2.0;

本格的なレンズエフェクトを使う場合、その他の詳細なカメラの情報も必要になると思われますが、その際は単位に注意が必要です。focus distance や aperture は cm、focal length は mm といった具合に、メートルじゃない上に統一されていません。

  • PolyMesh (IPolyMesh, OPolyMesh)

ポリゴンメッシュの schema です。多くの場合これをレンダリングするのが Alembic をインポートする主目的になるでしょう。そして、やはりこやつもプリレンダ向けなデータ構造をしており、リアルタイム用途では扱いづらいものとなっています。障害になるのは主に次の 2 点だと思われます。

任意の n 角形を保持できる。混在もできる
3 角形だったり 4 角形だったり、100 角形である可能性すらなきにしもあらずで、それらが 1 つのモデルの中に混在しています。
何角形かという情報は IPolyMeshSchema::Sample::getFaceCounts() で取れるコンテナに保持されており、これが例えば [4,3,5] だった場合、index[0-3] が最初の 4 角系、index[4-6] が次の 3 角系、index[7-11] が最後の 5 角系、といった具合です。レンダリングの際は 3 角形化などの中間処理を入れることになるでしょう。

位置、法線、UV、それぞれ個別にインデックスを持てる
位置、法線、UV の要素数がそれぞれ異なる可能性があり、それぞれが個別にインデックスを持っている可能性があります。(インデックスの要素数は全て一致)
モデリングツールにおけるポリゴンの扱いを考えると、なるほど、という感じですが、レンダリングする側としては厄介です。インデックスが独立している場合、インデックスを展開するか、D3D11 世代以降の機能を利用して独自に頂点&インデックスバッファを用意して複数インデックスに対応させる必要があります。(インデックス展開 = index の要素数分の各要素の配列を用意して、for(int i=0; i<num_indices; ++i) { new_pos[i]=pos[index[i]]; } みたいに展開する)
前者はメモリ使用量が膨れ上がる上にたぶん速度も後者より劣りますが、既存のシェーダをそのまま使えます。一方後者は実行効率はいいものの頂点シェーダで特殊なことをやる必要があり、専用のシェーダを整備する必要が生じます。AlembicImporter では Unity 側で一般的なシェーダでレンダリングできるようにする必要があったため、インデックスを展開する方法を採っています。

他には、トポロジが変化する可能性があることもリアルタイム用途としてはやや特殊です。
IPolyMeshSchema::getTopologyVariance() でトポロジがどう変化するかという情報が取れます。これが kConstantTopology であれば頂点の位置もインデックスも一切変化しない static な mesh であり、kHomogenousTopology であれば頂点の位置は変化するもののインデックスは変化しない mesh であり、kHeterogenousTopology であれば頂点の位置もインデックスも変化する mesh です。
kHeterogenousTopology が一番取り扱いが面倒です。インデックスの内容だけでなく、インデックスや頂点の要素数も変化する可能性があります。事前に全サンプルを巡回して最大数を算出し、その分のバッファを割り当てておく、などの対応が必要になるかもしれません。

  • Points (IPoints, OPoints)

パーティクルなどを格納する schema です。sample が持つデータは位置、速度、ID、バウンディングボリュームだけで、リアルタイム用途でもほぼ中間処理なしで使えるでしょう。
ID はパーティクルの数が増減するときの識別用に使います。例えば位置や色をランダムにずらしたいような場合、パーティクルのインデックス (0-n) を乱数のシードにすると、パーティクルの増減に伴ってランダム値も変わってしまいます。こういう場合 ID を乱数のシードにすることで、パーティクルが増減してもランダム値に影響は出ないように改善できます。

Export

インポートができたらならエクスポートも簡単です。schema オブジェクトを作り、sample を内容を埋めて set するだけです。ただ、いくつか注意が必要な点もあります。

読み込み用のオブジェクト群と書き込み用のオブジェクト群ではなぜか class の関係が若干違っています。
読み込みの場合、IObject と schema class 群 (IXform など) は継承関係がなく、IObject が shcema を保持する形になっています。一方、OObject と schema class 群 (OXform など) は継承関係になっています。なので、読み込み用 class 群だけを見た場合 1 つのオブジェクトが複数の schema を持てるように見えますが、実際にはそういうケースはないと見ていいようです。

Alembic はデフォルトでは 1 sample / seconds (= 1 FPS) で記録するようになっています。これを変えるには TimeSampling を作り、各 schema にその TimeSampling のインデックスを与える必要があります。
TimeSampling とは、サンプルがどういう周期で記録されているかという情報です。これには 3 種類のモードがあります。

Uniform: 時間間隔が均一なモードです。開始時間とインターバルだけが記録され、n 番目のサンプルの時間は 開始時間 + インターバル * n になります。
Acyclic: 時間間隔が不均一なモードです。この場合全てのサンプルの時間がそのまま配列として記録されます。サンプル間には delta time が負になってはいけない以外に規則性はありません。
Cyclic: 時間間隔が不均一ながら周期的なモードです。インターバルとサンプル (配列) が記録されます。 例えばインターバルが 1.0、サンプルが [0.0, 0.1, 0.3] であれば、0.0, 0.1, 0.3, 1.0, 1.1, 1.3, ... という時間分布になります。

多くの場合は Uniform で事足りると思われますが、例えばゲームを 3D 録画したい場合は Acyclic にするとゲームへの影響を最小限にできるでしょう。
Cyclic は主にプリレンダにおけるモーションブラーのためのモードだそうで、シャッターオープン開始、シャッターフルオープン、シャッタークローズ開始、シャッターフルクローズ、のサイクルを繰り返し記録したいような時使うそうです。リアルタイムからは最も縁遠いモードと言えそうです。
TimeSampling は混在もでき、schema 単位はもちろん、property 単位で個別に設定できます。

object space raymarching

ProceduralModeling

夏コミ版 exception reboot で用いた、オブジェクトスペースでレイマーチする手法について解説してみます。(ここで言うレイマーチは厳密には sphere tracing のことですが、面倒なのでレイマーチで統一します)
アイデア自体は特に新しくも難しくもなく、シーン内に単純なポリゴンモデルを配置し、そのモデルの中でレイマーチを行うというものです。レイマーチにより G-Buffer を生成し、あとは通常通りライティングを行います。

レイマーチについては過去にこの blog でも簡単な入門を書きました。最近では日本語でも結構情報が出てくるくらい知名度が上がってきているように見受けられます。
raymarching for games - primitive: blog
レイマーチで G-Buffer を生成する手法もしばらく前にこの blog で紹介しました。
rendering fractals in Unity5 - primitive: blog

実装の詳細の大部分は既にこれらの記事中に書かれており、本記事ではそれらとかぶってる部分は省略します。また、成果物はこちらになります。上の画像のシーンがそのまま含まれています。(ProceduralModeling.unity)
https://github.com/i-saint/Unity5Effects
ちなみに作例は Unity で作られていますが、レンダリング方式が deferred shading かつ自作シェーダを動かせる環境であればどこでもこの手法は使えるはずです。

レイマーチで distance function を描く場合、fullscreen quad で画面全体を描くことが多いですが、それを小分けにしてオブジェクトスペースでやろうというのが今回の主題です。これにより、レイマーチの欠点である柔軟性のなさと計算量の多さの緩和を狙います。
レイマーチでは、シェーダに式を書くことで図形を表現するため、図形の変更はつまりシェーダソースの変更になります。外部パラメータ化によりある程度変更可能にはできますが、例えば特定の位置に図形を追加、みたいな変更に耐える柔軟性を持たせるのは結構大変です。今回のオブジェクトスペースのやり方であれば、オブジェクトの配置や位置の変更はゲームエンジンのエディタ上でオブジェクトを追加したり TRS をいじることで簡単に実現できます。
画面全体をレイマーチする場合、レイの開始点は常にカメラ位置 (もしくは view plane) であり、図形に到達するまでに相応のステップ数が必要になります。一方オブジェクトスペースの場合、レイの開始点がそれなりに目的の図形の近くであることが期待でき、ステップ数削減が見込めます。

以下実装の手順。
まず、描画に使うポリゴンモデルは立方体か球だと都合がいいです。これらの場合ピクセルシェーダ内の単純な計算でレイがオブジェクトの内側にいるかが判定できます。もっと複雑なモデルでも適用は可能ですが、その場合内外判定に裏面の depth の情報が必要になります (最大マーチ距離などで適当に近似するのもアリではありますが) 。内外判定を怠ると オブジェクトの裏面=無限遠 と扱われることになり、シーンによっては顕著にヘンな結果になります。

シェーダ内の処理。頂点シェーダの出力に、頂点のワールド座標と法線を追加します。このワールド座標 (=ポリゴン表面位置) をレイマーチの開始位置とするわけです。
ピクセルシェーダでレイマーチを行いますが、この中でレイをオブジェクトのローカル座標系に変換する処理を挟みます。具体的には、レイをオブジェクトの TRS の逆行列で変換し、scale だけ掛けなおします (TR の逆行列でもいいと思われますが。Unity の場合 TRS の逆行列はエンジン側が提供しているのでこうしました)。レイの方向は normalize(頂点のワールド座標-カメラ位置) で算出できます。
注意すべき点として、レイの開始点が既に図形の内側だった場合の対処が必要です。これを怠ると法線が反転した面が出てきてライティングがおかしくなります。対処法としては、レイマーチの結果得られた距離が負 (=レイは図形の内側) である場合、元ポリゴンモデルの法線を出力する、というものになります。

あとは通常通り distance function を捏ねていい感じに見える形状を構築します。distance function 以外は定型的な処理なので、distance function だけ書いてパラメータを設定すればいい感じにレンダリングできるようなフレームワークを作るといいでしょう。私の場合この 2 つがフレームワークになっています:ProceduralModeling.cginc Framework.cginc

レイがオブジェクトの外に行ってしまったら discard、とすることで、元ポリゴンモデルと輪郭が異なる形状を表現できます。(前述の 立方体か球だと都合がいい のはこの内外判定が容易なため)
SV_Depth を用いてレイの到達点の depth を出力すれば、影や立体交差も正しく表現できるようになり、SSAO などもいい感じに載るようになります。ただし、early Z culling が効かなくなるため顕著に重くなります。また、元ポリゴンモデルとの形状がかけ離れるほど、正しく形状を表現するのに多くのステップ数が必要になり、重くなります。ここらへんはバランスを考える必要がありそうです。
https://raw.githubusercontent.com/i-saint/Blog/master/ObjectSpaceRaymarching/depth.jpg
例:元ポリゴンモデルとの形状の差が激しい状況。左は元ポリゴンモデルの depth をそのまま出力したもの。右は SV_Depth でレイの到達位置を出力したもの。見ての通り左は交差部分がおかしいとか色々エラーがあります。が、右よりだいぶん速いです。

作例:
https://raw.githubusercontent.com/i-saint/Blog/master/ObjectSpaceRaymarching/cube.jpg
レイマーチを試したことがある方なら一目で式が分かると思います。使い古された形状ですが、画面密度を増やすには割と効果的なんじゃないかと。経過時間で縦方向の座標を動かせばそれっぽくアニメーションになります。

https://raw.githubusercontent.com/i-saint/Blog/master/ObjectSpaceRaymarching/hex.gif
みんな大好き六角形。歯抜けにしたりランダム高低差をつけたりすると見た目の華やかさが上がります。六角形の中心が元図形の範囲外だったら描画しない、とすることで、元ポリゴンモデルをスケールするだけで六角模様が追随するようになります (=元モデルの境界で六角形がぶった切られない)。

https://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/Metaball.gif
メタボール。これはレイマーチだとすっごく単純な式で実現できます (実装のコア部分)。ただし、これだとボールの数に比例してすごい勢いで重くなっていきます。数百や数千を想定する場合もっと高次のアルゴリズムを考える必要がありそうです。


冒頭の画像のシーンは、これらにポストエフェクトやパーティクルを盛りまくることでできあがりました。

Unity ちゃん以外の背景オブジェクトは全て立方体をピクセルシェーダで加工したものです。
下図はシェーダによる差を示したもので、両方とも同じシーンで左はシェーダを通常の Standard Shader に差し替えたものになります。
https://raw.githubusercontent.com/i-saint/Blog/master/ObjectSpaceRaymarching/before.pnghttps://raw.githubusercontent.com/i-saint/Blog/master/ObjectSpaceRaymarching/after.png

ここまでやってみた実感として、レンダリング結果や編集のしやすさは良好なんですが、この方法でも図形が複雑になると重くて、もう一工夫入れないとちょっと実用に耐えられない印象です。上記画像のシーンは六角模様の床が明らかに危険域の重さになっています…。
いくつか最適化のネタはあるものの、まだ検証中の段階でそちらの詳細はまた別の機会に書く予定です。まあ adaptive sub-sampling、temporal、screen space normals、early-Z culling モドキ など、ネタとしてはごくありふれたものになりそうです。
(ちなみに夏コミ版 exception reboot には時間がなくて一つも入れられませんでした…。つまり最終的に今よりは描画負荷は軽くなる見込みです)


今後ハードウェアの進化に伴ってこういうポリゴン以外の 3D 形状の表現がだんだん市民権を得ていくんじゃないかと思います。
ただ、今回みたいなプロシージャル的なネタはなかなか一般化が難しく、タイトル固有の実装にならざるを得ないところが多いです。逆に言うと、ゲームエンジンが当たり前に使える時代にプログラマーが最も輝けるのがこういう領域なのかもしれません。

before C88

C88 で頒布する exception reboot 開発途中版の PV です。R23b でお待ちしております。
今回も体験版です…。前回同様 100 円でソースコード (Unity プロジェクト一式) 同梱です。カワノさんによる新曲も追加されています。
残念ながらゲーム内容は冬コミ版とほとんど変わっていません。しかしグラフィックは大きく変わりました。見ての通りだいぶん派手になっています。相応に重くもなっていますが、もう低スペックは切り捨てて全力で自分の目指す絵を追求する方向に舵を切りました。ただ、冬コミ版は非 D3D11 環境では色がおかしかったりパーティクルが出なかったりしたのが、今回は一応一通り機能するようになっています。


Unity 5 移行に伴い、レンダリングに絡む処理は全て書き直しています。このせいでえらい手間取ってしまいましたが、標準レンダリングパイプラインで描画するように変更したことで、他のプロジェクトへ簡単にアセットを流用できるようになりました。労力に見合う成果は出せたんじゃないかと思います。(冬コミ版の時はレンダリングパイプラインを全て独自に実装してオレオレ deferred shading で描いていたので、他プロジェクトへの流用は困難でした)
レンダリング関連のソースはほとんどこちらで公開しているので、ソースが目的であればこちらで事足りると思います。
https://github.com/i-saint/Unity5Effects

今回の絵作りの肝となっているのが、オブジェクトスペースでレイマーチするという小技。要するに立方体をピクセルシェーダで加工して複雑な図形に見せかけるというものですが、思ってた以上にいい結果が得られたので、後日詳細を書いてみようと思います。
exception reboot

playing with Unity 5's deferred shading pipeline

WaterSurface2
Unity4 の時に書いた自作 deferred shading 用の機能を Unity5 用に再実装しています。その経過の記録です。成果物はこちら https://github.com/i-saint/Unity5Effects

Screen Space Boolean

https://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/Boolean.gifhttps://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/BooleanDesc.gif
ステンシルや ZTest Greater の組み合わせでスクリーンスペースでブーリアン演算をやるというトリック。以前解説したものですが、Unity5 でこれをやるのは苦労しました。
実装の方針として、ブーリアン演算の段階では depth だけを出力し、実際に G-Buffer を書き出す段階では ZTest Equal を使用するように変更を加えた Standard Shader を用います。これは後述のステンシルの制限の回避や、必要なシェーダのバリエーションを最小限にするためです。代償として drawcall の数が増えています。

まずステンシル。deferred の G-Buffer 生成からライティングの段階では、カメラのデフォルトのレンダーターゲットに対するステンシルの書き込み指定は無視されるようです。(おそらく Unity が内部的にライトとオブジェクトの組み合わせをステンシルで実現しているため) これをなんとかするため、レンダーターゲットを別に用意してそこでブーリアン演算を行い、depth を元のレンダーターゲットにマージしています。

ZTest Equal を使う場合注意が必要で、Unity で SV_Depth で depth を出力する場合、デフォルトの depth と出力を一致させることはたぶんできないようです。そんなはずはないと思いたいのですが、方法を見つけられませんでした。(D3D11 に限れば SV_POSITION の z を使えば一致しましたが、他プラットフォームでは使えません) このせいでずいぶん悩むことになりました。
ブーリアンの減算は、貫通していたら depth を初期化する処理が必要になります。これをストレートに実装すると

if(減算する側のモデルの depth > 減算される側のモデルの裏面の depth) {
    output_depth = 1.0; // 貫通しているので depth 初期化
}
else {
    output_depth = 減算する側のモデルの depth;
}

こんな感じになります。しかし、depth を自力で算出して出力すると精度がズレるようで、その後 ZTest Equal で同モデルを書くと細かい穴があいてしまいました。しょうがないので貫通処理は 2 パスに分けることで対処しました。最初のパスではごく普通に depth を出力。2 パス目では

if(減算する側のモデルの depth > 減算される側のモデルの裏面の depth) {
    output_depth = 1.0;
}
else {
    discard;
}

とすることで、貫通箇所以外は SV_Depth を経由しないようにします。

あとは影。これは対処不能という結論に至りました。
Unity の影は描画順がランダムのようで、こういう 不透明だけど描画順に依存する 代物はまともな方法では対処できないようです。CommandBuffer にも影バッファ生成の前後に差し込めるイベントはなく、たぶんライティング処理を全部自力で書いて影バッファ生成をフルコントロールしないと実現不可能だと思われます。

Screen Space Shadow

https://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/ScreenSpaceShadows.gif
スクリーンスペース全方位影。これを実現するにはライティング処理を自力で書く必要がありますが、幸いそこはいい公式サンプルがあるのでごっそりコピペしました。あとはごく普通にライトの中心からピクセル位置までレイマーチして途中で遮られたら光量を減らすだけです。

一つ問題になったのが、カメラが HDR ではないときへの対処。これは上記公式サンプルでも未対処でした。
HDR が有効なときとそうでないときでは内部処理が変わります。HDR が有効な場合、G-Buffer の emission buffer にはカメラのデフォルトのターゲットが使われており、HDR 無効の場合専用のバッファが用意されます。よって、レンダーターゲットを自力で G-Buffer に設定したいような場合、以下のような場合分けが必要になります。

Camera cam = GetComponent<Camera>();
CommandBuffer cb = new CommandBuffer();
cb.SetRenderTarget(new RenderTargetIdentifier[] {
    BuiltinRenderTextureType.GBuffer0,
    BuiltinRenderTextureType.GBuffer1,
    BuiltinRenderTextureType.GBuffer2,
    cam.hdr ? BuiltinRenderTextureType.CameraTarget : BuiltinRenderTextureType.GBuffer3
}, BuiltinRenderTextureType.CameraTarget);

また、HDR が無効な場合、emission は logarithmic encoding という方法でエンコードされた状態で保持されます。これはより精度を保つためのエンコード方法だそうで、色を exp2(-color) で加工した状態で保持し、後で -log2() で復元する、というものです。(詳細)
この影響で、HDR の有無でブレンド方法とピクセルシェーダの出力を変える必要があります。HDR 無効な場合のブレンド方法は Blend DstColor Zero、有効な場合 Blen One One です。Properties で外部から変えられるようにした方がいいでしょう。ピクセルシェーダには以下のような処理を加えます。

half4 frag(v2f I]N) : SV_Target
{
    half4 color;
    // ...
#ifndef UNITY_HDR_ON
    color = exp2(-color);
#endif
    return color;
}
#pragma multi_compile ___ UNITY_HDR_ON

まあしかし、古いモバイルデバイスを視野に入れない限りは "強制的に HDR をオンにする" が一番簡単な対処法であると思われます。

あと、今回はやっていませんが、こういう放射状のレイマーチには epipolar sampling という強力な最適化方法があるのを最近知りました。
アルゴリズムの詳細はこの記事の 3. Epipolar sampling 以降で詳しく解説されています。こちらの Unity 用 Light shaft はこのテクニックを用いてるそうです。いずれ試してみたいところです。

Screen Space Reflections

https://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/ScreenSpaceReflections.jpg
以前実装の詳細を書いたやつです。幸いほぼそのまま移植できました。
障害になったのが、depth からピクセルの位置を算出する方法。ライティングの処理の中に同等処理があるのですが、なんだかよくわからない難解な処理になってる上、ポストエフェクトの場合これは機能しないように見えます。
結局ストレートに view projection の逆行列をシェーダに渡して対処しました。 projection 行列は Unity が内部的に加工を施すので、逆行列を求める際にも同じことをやる必要があります。

// C# 側処理
var cam = GetComponent<Camera>();
Matrix4x4 view = cam.worldToCameraMatrix;
Matrix4x4 proj = cam.projectionMatrix;
// Unity が内部でやってる projection 行列の加工
proj[2, 0] = proj[2, 0] * 0.5f + proj[3, 0] * 0.5f;
proj[2, 1] = proj[2, 1] * 0.5f + proj[3, 1] * 0.5f;
proj[2, 2] = proj[2, 2] * 0.5f + proj[3, 2] * 0.5f;
proj[2, 3] = proj[2, 3] * 0.5f + proj[3, 3] * 0.5f;
Matrix4x4 inv_view_proj = (proj * view).inverse;
Shader.SetGlobalMatrix("_InvViewProj", inv_view_proj);
// shader 側処理
sampler2D_float _CameraDepthTexture;
float4x4 _InvViewProj;

float4 GetPosition(float2 uv)
{
    float2 screen_position = uv * 2.0 - 1.0;
    float depth = tex2D(_CameraDepthTexture, uv).x;
    float4 pos4 = mul(_InvViewProj, float4(screen_position, depth, 1.0));
    return pos4 / pos4.w;
}

また、D3D9 ではなぜか _CameraDepthTexture はバイリニアフィルタがかかっているようで、そのままだとなんかヘンな結果になります。ポイントフィルタに切り替えたいところですが、まともな方法ではできそうにないので、ピクセルの中心をサンプリングすることで対処しました。(_ScreenParams.zw-1.0)*0.5 がピクセルサイズの半分になるのでこれを利用しています。

ちなみに Screen Space Reflection は近い先に標準搭載される予定です。この分野でたぶん世界一詳しい人が実装を担当しているのでそれを待っていたのですが、当初の予定から伸びに伸びている上、近々に必要になったので結局自作のを移植しました。
Roadmap によると 5.2 (2015/09/08 リリース) にリストされていますが、delayed になっています…。

Rim Light

https://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/RimLight.jpg
法線と カメラ -> ピクセル位置 の角度が浅いところを明るくするアレです。全く根拠レスな処理なのになんかカッコよく見えるようになります。
前述のピクセル位置の復元さえできれば特に難しいところはないですが、今更ながらフレネル反射の式 (正確にはそれの簡易版) を今回導入してみました。

// ...
float3 N = tex2D(_CameraGBufferTexture2, uv) * 2.0 - 1.0;
float3 I = normalize(pixel_pos.xyz - _WorldSpaceCameraPos.xyz);
float fresnel = saturate(_FresnelBias + pow(dot(I, N) + 1.0, _FresnelPow) * _FresnelScale);

いい感じになった気がしますが、最大のメリットはパラメータの調整がアーティストフレンドリーになることでしょうか。

Water Surface & Caustics Field

https://raw.githubusercontent.com/i-saint/Unity5Effects/master/doc/WaterSurface.gif
Tokyo Demo Fest 2015 用に作った demo の水面のリファイン版。3D ノイズとレイマーチでいい感じに見せかけるストレートな実装です。
こちらもピクセル位置の復元以外 Unity5 固有の難しい問題はありませんでした。フラクタル図形や Rim Light と組み合わせると実にいい感じになります。


そんなわけで、予想より大分苦労しましたが、主要な機能は大体移植できました。既存プロジェクトへの組み込みが簡単になったのが大きな成果です。

あと、夏コミ参加します。場所は 日曜日 R 23b です。
今回も exception reboot の製作途中版でソース同梱の予定です。Unity 5 移行に伴い、レンダリング部分を全面的に書きなおしています。結果としてゲーム内容は前回からあまり変わらないものになりそうです…。

recent works

最近作ったもの紹介。全て Unity から使う前提で作ってはいますが、BatchRenderer 以外は Unity 依存度は低く、一応他のプログラムから使うこともできるはずです。また、全て MIT ライセンスか CC BY ライセンスなので、ソースレベルで切り貼りして使うのも問題ないようになっています。

PatchLibrary

https://github.com/i-saint/PatchLibrary
実行中のプログラムにロードされている dll を強引に更新するツールです。プラグイン開発などの際、VisualStudio ビルド後イベントでこれを呼び、実行中のホストプログラムを再起動することなく更新する、というような使い方を想定しています。
以前作った DynamicPatcher の簡易版的なもので、dllexport されている関数のテーブルを書き換えて新しい dll の関数にリダイレクトさせることで更新を実現しています。このため、dll で作成したコンテキストを保持しているような場合に途中でこれで更新すると死ぬなどの問題がありますが、その点を理解して使えばそれなりに有用です。以降のツール群はこれ使いながら開発しています。

FrameCapturer

https://github.com/unity3d-jp/FrameCapturer
https://raw.githubusercontent.com/unity3d-jp/FrameCapturer/master/Screenshots/gif_capturer.gif

ゲーム画面をキャプチャして gif アニメや exr などに出力する代物です。後述の TweetMedia を同梱しており、画面のキャプチャから Twitter への投稿までを全てゲーム内からできるようになっています。
開発ブランチでは mp4 キャプチャも実装中で、まだ若干動作が怪しいものの、キャプチャから Twitter への投稿まで機能するようになっています。(Twitter は web からは動画は投稿できないものの、API レベルでは mp4 を受け付けるようになっています) こちらも近い先ちゃんと公開する予定です。

TweetMedia

https://github.com/unity3d-jp/TweetMedia
前述の通り、Twitter にメッセージと画像/mp4 動画を投稿できるようにする代物です。
既にいい実装があればそれを使いたかったんですが、画像/動画まで対応したのは見つけられなかったので仕方なく自作しました。(と言っても大部分は twitcurl を流用していますが)


BatchRenderer

https://github.com/i-saint/BatchRenderer
https://raw.githubusercontent.com/i-saint/BatchRenderer/master/Screenshots/procedural_gbuffer.gif
こちらは単なる機能追加ですが。Unity 5.1 で CommandBuffer.DrawProcedural() が加わり、やっと Unity の描画パイプラインからインスタンシング描画をできるようになりました。これを使って DrawProcedural() で G-Buffer を生成する機能を追加した、というものです。
本物のインスタンシング描画なので当然ながら何十万インスタンスであろうと 1 drawcall で描けはするのですが、残念ながら速度面では擬似インスタンシングとあまり変わらないようです。DrawProcedural() はその不思議な仕様上、頂点データへのアクセスに間接参照が必要になるため、理想的なケースに比べると一歩パフォーマンスが落ちる感触です。それでも実用に足る速さは出てるとは思うんですが、surface shader が使えなくなるデメリットを覆せるほどのメリットが擬似インスタンシングに対してあるかと問われると、現状返答に詰まる感じです…。
スクリーンショットは SD Unity ちゃんを 2000 体 (約 16,000,000 triangles) 出してみたもので、このくらいがうちの GeForce GTX 970 で 60 FPS 出せるギリギリのラインになっています。

AlembicImporter

https://github.com/unity3d-jp/AlembicImporter
https://raw.githubusercontent.com/unity3d-jp/AlembicImporter/master/Screenshots/alembic_example.gif

Unity 上で Alembic ファイルを再生する代物です。
Alembic とは映像業界でよく使われているファイル、およびそのライブラリで、シーン上のオブジェクトをフレーム毎に全て頂点キャッシュへ bake して格納するのに使われます。映像業界ではこれを用いてシミュレーション結果を bake してレンダラなどに渡すような使い方をするそうです。



当分仕事が映像関係になる予定なので、今後しばらくこの blog の内容もそれっぽいものになるかもしれません。AlembicImporter はまさに映像関係者向けですし、FrameCapturer も G-Buffer を exr でエクスポートしたいという要望を受けたのが始まりでした。映像方面はゲーム業界とは近いようで全く違う文化圏を持っていて、おもしろい発見が多いです。

introdunction to SIMD programming

Screenshot
Unite 2015 Tokyo の講演で詳細を話せなかったのが心残りだったので、大量のオブジェクトの更新処理についてこの場で書いてみます。

主に C++ で、簡単なパーティクルエンジンを作り、それを SIMD を用いて高速化する手順を解説します。
話を簡単にするため、以下の前提を設けます。
・x86 環境のみ考慮
・パーティクルは位置と速度のみを保持
・パーティクル同士の相互衝突は総当たりで計算
総当たりなので超遅いですが、実装は容易で SIMD による恩恵を受けやすく、題材として手頃です。

この記事の中で引用されているソースの元は こちら、ビルド結果 (上のスクリーンショットのデモプログラム) は こちら になります。

相互衝突するパーティクルを実装する場合、お互いの距離を計算し、当たっていたらめり込み具合に応じて押し返す、というのがよくある実装だと思います。まずはそれをストレートに C++ で実装してみます。

const int m_particle_count = 4096;
Particle m_particles[m_particle_count];
float m_particle_size = 0.1f;
float m_pressure_stiffness = 500.0f; // 押し返し係数

float dt = 1.0f / 60.0f;

// 速度を更新。パーティクル同士の押し返し
for (int i = 0; i < m_particle_count; ++i) {
    float3 &pos1 = m_particles[i].position;
    float3 accel = { 0.0f0.0f0.0f };

    for (int j = 0; j < m_particle_count; ++j) {
        float3 &pos2 = m_particles[j].position;
        float3 diff = pos2 - pos1;
        float dist = length(diff);
        if (dist > 0.0f) { // dist==0.0 : 自分自身との衝突なので省略
            float3 dir = diff / dist;
            // めり込んでいたら反対方向に加速度を追加
            float overlap = std::min<float>(0.0f, dist - m_particle_size*2.0f);
            float3 a = dir * (overlap * m_pressure_stiffness);
            accel = accel + a;
        }
    }

    float3 &vel = m_particles[i].velocity;
    vel = vel + accel * dt;
}

// 位置を更新。マルチスレッド対応を考慮して速度更新とは分ける
for (int i = 0; i < m_particle_count; ++i) {
    float3 &pos = m_particles[i].position;
    float3 &vel = m_particles[i].velocity;
    pos = pos + (vel * dt);
}

最適化を有効にしてビルドし、とりあえずシングルスレッドで 4096 個のパーティクルを相互に衝突させてみたところ、所要時間は以下のようになりました。

Plain C++  (ST) : 140.0ms

これがスタートラインです。これをがんばって速くしていきます。
ついでに、適当に重力加速と床とのバウンド処理も入れて Unity で可視化してみるとこうなりました。単純実装ながらなかなかいい感じではないでしょうか。

SIMD

近年では大抵の CPU には 1 命令で複数のデータを同時に計算する命令群が備わっています。これを使うと、例えば 4 つの float の足し算や掛け算を一度に行うことができます。これらは SIMD (Single Instruction Multiple Data) と呼ばれ、CPU プログラミングでシビアに速度を求める場合避けては通れない道になります。 

x86 系 CPU には SSE という SIMD 命令群が備わっており、今回の例ではこれを用います。(x86 用なので ARM など他の CPU では動きません。ARM の CPU も大抵は NEON という SIMD 命令群が備わっており、マルチプラットフォーム対応するには個別に実装するなりライブラリを用いるなりする必要があります)
SSE の命令群は intrinsic と呼ばれる関数を用いることで C++ から使うことができます。最初のストレートな C++ 実装を、ベクトル演算を SSE の命令で置き換える形で書きなおすと以下のようになります。

__m128 particle_size2 = _mm_set_ps1(m_particle_size*2.0f);
__m128 pressure_stiffness = _mm_set_ps1(m_pressure_stiffness);
__m128 dt = _mm_set_ps1(dt_);
__m128 zero = _mm_setzero_ps();

// パーティクル同士の押し返し
for (int i = 0; i < m_particle_count; ++i) {
    __m128 pos1 = _mm_load_ps((float*)&m_particles[i].position);
    __m128 accel = zero;
    for (int j = 0; j < m_particle_count; ++j) {
        __m128 pos2 = _mm_load_ps((float*)&m_particles[j].position);
        __m128 diff = _mm_sub_ps(pos2, pos1);
        __m128 dist = length(diff);
        __m128 dir = _mm_div_ps(diff, dist);
        __m128 overlap = _mm_min_ps(zero, _mm_mul_ps(_mm_sub_ps(dist, particle_size2), pressure_stiffness));
        __m128 a = _mm_mul_ps(dir, overlap);
        accel = _mm_add_ps(accel, select(zero, a, _mm_cmpgt_ps(dist, zero)));
    }

    __m128 vel = _mm_load_ps((float*)&m_particles[i].velocity);
    vel = _mm_add_ps(vel, _mm_mul_ps(accel, dt));
    _mm_store_ps((float*)&m_particles[i].velocity, vel);
}

// 位置を更新
for (int i = 0; i < m_particle_count; ++i) {
    __m128 pos = _mm_load_ps(m_particles[i].position.v);
    __m128 vel = _mm_load_ps(m_particles[i].velocity.v);
    pos = _mm_add_ps(pos, _mm_mul_ps(vel, dt));
    _mm_store_ps(m_particles[i].position.v, pos);
}

__m128 というのが SSE 用のデータ型で、float が 4 つ格納できる 16 byte のデータになっています。
_mm で始まる関数群が SSE の intrinsic で、大体はこの intrinsic が SSE の命令に直接対応するようになっています。アンダースコアまみれで見にくいですが、関数名に mul とか add とかついているので初見の方でも大体何やってるのかわかるんじゃないかと思います。

注意が必要な点がいくつかあります。
_mm_load_ps() はメモリから SSE のレジスタにデータを移す命令ですが、ロード先のアドレスは 16 byte にアラインされている必要があります。されていない場合クラッシュします。アライン不要の _mm_loadu_ps() というのもありますが、若干速度にペナルティがあるので今回は使っていません。これらは SSE のレジスタからメモリにデータを移す _mm_store_ps() にも同じことが言えます。このため、今回の Particle は position も velocity も 4 次元ベクトル (というか x,y,z + パディング) としています。
SIMD プログラミングでは、例えば x+y+z みたいなベクトルの水平方向の計算をする場合ややこしいことをやる必要があります。こういう処理の典型的な例が内積です。3 次元の内積のストレートな C++ 実装はこうですが、

inline float dot(float3 v1, float3 v2)
{
    return v1.x*v2.x + v1.y*v2.y + v1.z*v2.z;
}

SSE だとこうなります。

inline __m128 dot(__m128 v1, __m128 v2)
{
    __m128 d = _mm_mul_ps(v1, v2);                            // d: x,y,z,w
    __m128 t = _mm_shuffle_ps(d, d, _MM_SHUFFLE(2, 1, 2, 1)); // t: z,y,z,y
    d = _mm_add_ss(d, t);                                     // d: xz,...
    t = _mm_shuffle_ps(t, t, _MM_SHUFFLE(1, 1, 1, 1));        // t: y,y,y,y
    d = _mm_add_ss(d, t);                                     // d: xyz,...
    return _mm_shuffle_ps(d, d, _MM_SHUFFLE(0, 0, 0, 0));     // xyz,xyz,xyz,xyz
}

// 上の SSE 化したパーティクルコードの length()
inline __m128 length(__m128 v)
{
    return _mm_sqrt_ps(dot(v,v));
}

_mm_shuffle_ps() が要素を交換する命令で、これを使って要素を入れ替えつつ足しあわせています。見ての通り実装がめんどくさい上に並列性も少なく、コストの高い処理になります。これは SSE に限らず SIMD プログラミングではつきまとう問題で、できるだけ避けた方がいいですが必要な状況も多いです。(SSE 3 で水平方向計算の命令が追加され、SSE 4 では内積を計算するそのものズバリな命令までありますが、今回は一応一般的な SIMD プログラミングという体裁なので使いません)
あとは分岐。SIMD の場合、分岐は両パス計算して結果を選択した方がずっと速いケースが多いです。比較命令 (_mm_cmpgt_ps() = compare greater than など) は結果が真の要素は全ビットが立っている (0xffffffff) ので、これを使ってビット演算で結果を選択します。

inline __m128 select(__m128 v1, __m128 v2, __m128 mask)
{
    __m128 t1 = _mm_andnot_ps(mask, v1);
    __m128 t2 = _mm_and_ps(v2, mask);
    return _mm_or_ps(t1, t2);
}

例えば v1 が {0,1,2,3}、v2 が {4,5,6,7} で mask が {0, 0xffffffff, 0xffffffff, 0} の場合、この select() は {0,5,6,3} を返します。上の SSE 化したパーティクルのコードもこれで結果を選択しています。

SSE 化したソースをビルドして最初のケースと同条件 4096 パーティクル シングルスレッドで動かしてみました。

Plain C++  (ST) : 140.0ms
SIMD (SSE) (ST) : 46.6ms (new!)

実に 3 倍も速くなりました。

SoA

SIMD 化により大きく速くはなりましたが、上記の SSE 化したコードには色々無駄があります。position も velocity も実際には x,y,z の 3 要素しか使われていません。なので SIMD 演算の際 w 要素部分の計算は完全に無駄になっており、計算リソースを 25% も無駄にしていることになります。また、先述の dot() / length() の中の並べ替えも大きなロスになります。
これらを解決する手法として、SoA (Structure of Arrays) と呼ばれるデータ構造があります。

C++ で多数のオブジェクトを扱う場合、各オブジェクトに付随するデータを class や struct にまとめて配列にする、というのが一般的だと思われます。

struct Particle
{ 
     float pos_x, pos_y, pos_z;
     float vel_x, vel_y, vel_z;
};
Particle particles[particle_count];

このようなデータの持ち方は AoS (Array of Structures) と呼ばれることもあります。SoA はこれとは逆で、データの各要素を配列にする、というデータの持ち方になります。

float pos_x[particle_count];
float pos_y[particle_count];
float pos_z[particle_count];
float vel_x[particle_count];
float vel_y[particle_count];
float vel_z[particle_count];

直感的ではないデータの持ち方ですが、SIMD プログラミングではこれは大きなメリットがあります。4 つのオブジェクトを同時に計算できるようになるのです。先ほどの続きで、内積を例に考えてみましょう。SoA なデータに対して SSE で内積と距離を計算する場合、こうなります。

inline __m128 soa_dot(__m128 x1, __m128 y1, __m128 z1,  __m128 x2, __m128 y2, __m128 z2)
{
    __m128 tx = _mm_mul_ps(x1, x2);
    __m128 ty = _mm_mul_ps(y1, y2);
    __m128 tz = _mm_mul_ps(z1, z2);
    return _mm_add_ps(_mm_add_ps(tx, ty), tz);
}

inline __m128 soa_length(__m128 x, __m128 y, __m128 z)
{
    return _mm_sqrt_ps(soa_dot(x,y,z, x,y,z));
}

やってることはストレートな C++ 実装と同じですが、4 つのベクトルを同時に計算するようになっています。前述のベタな SSE 化の dot() と違い、shuffle は必要ないですし w 要素の計算が無駄になることもありません。ベタな SSE 化はベクトル演算それぞれを SSE 命令に置き換えていましたが、SoA の場合 4 つのパーティクルを同時に行うことで並列化するイメージです。
SoA にはスケーラビリティの面でも大きなメリットがあります。近年の x86 CPU には AVX という SIMD 命令群が追加されており、これを使うと 8 つの float を同時に計算できます。しかし 8 次元ベクトルなんてゲーム開発ではまず出てこないので、単純にベクトル演算を SIMD 命令に置き換えるアプローチではこれを活かすのは難しいです。一方、SoA であれば 8 オブジェクト並列に計算することで最大限に恩恵を受けることができますし、将来 16 並列や 32 並列 SIMD が出てきたとしても簡単に対応できます。

実際にパーティクルを SoA 化して計算してみます。まずパーティクルデータを SoA 化する必要があります。AoS <-> SoA の転置を行うには、SSE では以下のような処理を行います。

union soafloat4
{
    struct { __m128 x, y, z, w; };
    __m128 v[4];
};

inline soafloat4 soa_transpose(__m128 v0, __m128 v1, __m128 v2, __m128 v3)
{
    __m128 r1 = _mm_unpacklo_ps(v0, v1);
    __m128 r2 = _mm_unpacklo_ps(v2, v3);
    __m128 r3 = _mm_unpackhi_ps(v0, v1);
    __m128 r4 = _mm_unpackhi_ps(v2, v3);
    soafloat4 r = {
        _mm_shuffle_ps(r1, r2, _MM_SHUFFLE(1, 0, 1, 0)),  // {v0.x, v1.x, v2.x, v3.x}
        _mm_shuffle_ps(r1, r2, _MM_SHUFFLE(3, 2, 3, 2)),  // {v0.y, v1.y, v2.y, v3.y}
        _mm_shuffle_ps(r3, r4, _MM_SHUFFLE(1, 0, 1, 0)),  // {v0.z, v1.z, v2.z, v3.z}
        _mm_shuffle_ps(r3, r4, _MM_SHUFFLE(3, 2, 3, 2)) };// {v0.w, v1.w, v2.w, v3.w}
    return r;
}

これを用いてパーティクルを SoA 化し、衝突計算を行います。

__m128 particle_size2 = _mm_set_ps1(m_particle_size*2.0f);
__m128 pressure_stiffness = _mm_set_ps1(m_pressure_stiffness);
__m128 dt = _mm_set_ps1(dt_);
__m128 zero = _mm_setzero_ps();

// パーティクル同士の押し返し
for (int i = 0; i < m_particle_count; ++i) {
    __m128 pos1x = _mm_set_ps1(m_soa.pos_x[i]);
    __m128 pos1y = _mm_set_ps1(m_soa.pos_y[i]);
    __m128 pos1z = _mm_set_ps1(m_soa.pos_z[i]);
    __m128 accelx = zero;
    __m128 accely = zero;
    __m128 accelz = zero;

    for (int j = 0; j < m_particle_count; j += 4) {
        __m128 pos2x = _mm_load_ps(&m_soa.pos_x[j]);
        __m128 pos2y = _mm_load_ps(&m_soa.pos_y[j]);
        __m128 pos2z = _mm_load_ps(&m_soa.pos_z[j]);

        __m128 diffx = _mm_sub_ps(pos2x, pos1x);
        __m128 diffy = _mm_sub_ps(pos2y, pos1y);
        __m128 diffz = _mm_sub_ps(pos2z, pos1z);
        __m128 dist = soa_length(diffx, diffy, diffz);

        __m128 dirx = _mm_div_ps(diffx, dist);
        __m128 diry = _mm_div_ps(diffy, dist);
        __m128 dirz = _mm_div_ps(diffz, dist);

        __m128 overlap = _mm_min_ps(zero, _mm_mul_ps(_mm_sub_ps(dist, particle_size2), pressure_stiffness));
        __m128 ax = _mm_mul_ps(dirx, overlap);
        __m128 ay = _mm_mul_ps(diry, overlap);
        __m128 az = _mm_mul_ps(dirz, overlap);
        __m128 gt = _mm_cmpgt_ps(dist, zero);
        accelx = _mm_add_ps(accelx, select(zero, ax, gt));
        accely = _mm_add_ps(accely, select(zero, ay, gt));
        accelz = _mm_add_ps(accelz, select(zero, az, gt));
    }

    __m128 velx = _mm_set_ss(m_soa.vel_x[i]);
    __m128 vely = _mm_set_ss(m_soa.vel_y[i]);
    __m128 velz = _mm_set_ss(m_soa.vel_z[i]);
    velx = _mm_add_ss(velx, _mm_mul_ss(reduce_add(accelx), dt));
    vely = _mm_add_ss(vely, _mm_mul_ss(reduce_add(accely), dt));
    velz = _mm_add_ss(velz, _mm_mul_ss(reduce_add(accelz), dt));
    _mm_store_ss(&m_soa.vel_x[i], velx);
    _mm_store_ss(&m_soa.vel_y[i], vely);
    _mm_store_ss(&m_soa.vel_z[i], velz);
}

// 位置を更新
for (int i = 0; i < m_particle_count; i += 4) {
    __m128 posx = _mm_load_ps(&m_soa.pos_x[i]);
    __m128 posy = _mm_load_ps(&m_soa.pos_y[i]);
    __m128 posz = _mm_load_ps(&m_soa.pos_z[i]);

    __m128 velx = _mm_load_ps(&m_soa.vel_x[i]);
    __m128 vely = _mm_load_ps(&m_soa.vel_y[i]);
    __m128 velz = _mm_load_ps(&m_soa.vel_z[i]);

    posx = _mm_add_ps(posx, _mm_mul_ps(velx, dt));
    posy = _mm_add_ps(posy, _mm_mul_ps(vely, dt));
    posz = _mm_add_ps(posz, _mm_mul_ps(velz, dt));

    _mm_store_ps(&m_soa.pos_x[i], posx);
    _mm_store_ps(&m_soa.pos_y[i], posy);
    _mm_store_ps(&m_soa.pos_z[i], posz);
}

すっごく冗長になってしまいましたが、__m128  を float に置換してみると意味は明瞭になるんじゃないかと思います。
速度更新フェーズの内側のループで 4 パーティクルとの距離計算を同時に行い、最後に 4 つの結果を reduce add して統合しています。recuce_add() の中身はそのまんま __m128 の中の 4 要素を足しあわせているだけです。

inline __m128 reduce_add(__m128 v1)
{
    __m128 v2 = _mm_movehl_ps(v1, v1);                      // z,w,z,w
    v1 = _mm_add_ps(v1, v2);                                // xz,yw,zz,ww
    v2 = _mm_shuffle_ps(v1, v1, _MM_SHUFFLE(2, 2, 2, 2));   // yw, yw, yw, yw
    v1 = _mm_add_ss(v1, v2);                                // xzyw, ...
    return v1;
}

実行結果。

Plain C++  (ST) : 140.0ms
SIMD (SSE) (ST) : 46.6ms
SIMD SoA   (ST) : 20.0ms (new!)

ベタな SSE 化からさらに倍以上速くなりました。Update のたびに毎回 SoA に並べ替えて計算して AoS に戻しているのでトータル処理量は増えているはずですが、それでもこれだけ効果があります。

ISPC

先に少し触れたように、近年の x86 CPU には AVX と呼ばれる命令群が備わっており、これを使うと 8 つの float を同時に計算できます。理想的なケースでは本当に SSE のほぼ倍の速度が出せます。ぜひ使いたいところです。
しかし、AVX が使える CPU は Intel Sandy Bridge (2011) 以降であり、まだ十分に普及しているとは言えません。なので AVX を使うにしても SSE 実装と AVX 実装両方を用意して実行時に使える方を選ぶのが実用的なやり方になります。
とはいえ、両方実装するのは面倒です。今回の例くらい単純なケースでは SSE, AVX 両方実装するのも大したコストではありませんが、ここでは紹介もかねて ISPC を使ってラクをしてみます。

http://ispc.github.io/
ISPC は SIMD に特化したプログラミング言語です。
SIMD の各レーンを実行単位と見立てた記述をするようになっています。ISPC で記述してコンパイルすると、SSE, AVX, ARM (NEON) それぞれに最適化されたコードを出力することができます。SSE と AVX に関しては、両方の実装を用意して実行時に使える方を選ぶのを自動的にやってくれます。コンパイル結果は .obj と .h として出力され、簡単に C/C++ のプログラムにリンクできます。Windows はもちろん、Mac、Linux、PS4 向けにもクロスコンパイルできます。
詳しい解説はこちらに委ねますが、ここまで読んでくれた方であれば、最初のストレートな C++ 実装に近い記述で、上の SIMD & SoA と同じ内容にコンパイルできる言語、という説明がわかりやすいんじゃないかと思います。

SoA <-> AoS 化は C++ 側で行い、それを ISPC で書いたカーネルに渡して処理します。ISPC のコードは以下のようになります。

// パーティクル同士の押し返し
for(uniform int i=0; i<particle_count; ++i) {
    uniform float3 pos1 = {pos_x[i], pos_y[i], pos_z[i] };
    float3 accel = {0.0f, 0.0f, 0.0f};
    foreach(j=0 ... particle_count) {
        float3 pos2 = {pos_x[j], pos_y[j], pos_z[j] };
        float3 diff = pos2 - pos1;
        float dist = length(diff);
        float3 dir = diff / dist;
        float3 a = dir * (min(0.0f, dist-(particle_size*2.0f)) * pressure_stiffness);
        accel = accel + (a * select(dist > 0.0f, 1.0f, 0.0f));
    }

    uniform float3 vel = {vel_x[i], vel_y[i], vel_z[i] };
    vel = vel + reduce_add(accel) * dt;
    vel_x[i] = vel.x;
    vel_y[i] = vel.y;
    vel_z[i] = vel.z;
}

// 位置を更新
for(i=0 ... particle_count) {
    float3 pos = {pos_x[i], pos_y[i], pos_z[i] };
    float3 vel = {vel_x[i], vel_y[i], vel_z[i] };

    pos = pos + (vel * dt);
    pos_x[i] = pos.x;
    pos_y[i] = pos.y;
    pos_z[i] = pos.z;
}

実にすっきりしました。これを SSE をターゲットにしてコンパイルすると上の SoA & SSE のコードとほぼ同じ内容のバイナリになりますし、AVX をターゲットにすれば AVX 実装のできあがりです。
AVX をターゲットにする際の注意点として、SSE が 16 byte align が必要だったように、AVX だと 32 byte align が必要になります。なので今回は最初からパーティクルは 32 byte align でメモリを確保しています。
これをビルドして実行してみます。

Plain C++  (ST) : 140.0ms
SIMD (SSE) (ST) : 46.6ms
SIMD SoA   (ST) : 20.0ms
ISPC       (ST) : 12.3ms (new!)

さらに 1.6 倍も速くなりました。処理内容は SSE & SoA と同じなので、このスピードアップは AVX 化によるものと言えます。
ISPC の便利さは SIMD プログラミングを長いことやってる人ほど身にしみるんじゃないかと思います。一頃に比べるとコンパイラのバグリーな動作は少なくなり、対応プラットフォームも増えたので、実戦投入も検討に値するレベルになってきたんじゃないでしょうか。
余談ですが、ISPC の "SIMD の各レーンを実行単位に見立てる" という考え方はシェーダプログラミングにも当てはまり、近年の NVIDIA や AMD の GPU ではシェーダプログラムは 32 の SIMD レーン (NVIDIA が "Warp"、AMD が "Wavefront" と呼んでるもの) それぞれで並列に実行されるようになっています。なので、ISPC はシェーダプログラミングのパラダイムを CPU の SIMD で実装したもの、とも言えると思います。


そんなわけで、色々頑張った結果最初のストレートな C++ 実装から 10 倍くらい速くなりました。これをマルチスレッド化までやると最近の Core i7 の比較的いいやつだと 8192 パーティクルが 60 FPS で動くようです。(= 67108864 回の距離計算が 16ms 以内で終わる!)
Unite 2015 Tokyo でデモした MassParticle はハッシュグリッド法で近隣パーティクルのみと衝突計算する最適化を入れていますが、衝突計算自体は今回の例と全く同じになっています。

解説はここまでです。記事中のソースは簡略化しているので、より詳細に実装を知りたい場合 github に上がってるソースを見ることをおすすめします。SIMD はアルゴリズムが単純でメモリアクセスが少ないケースで特に効果的です。弾幕やパーティクルなど、単純なオブジェクトの大量処理は SIMD に向いた問題なので、数出したい方は足を踏み入れてみるといいでしょう。
あと、もっと速くする方法がありましたらご一報ください。

C#

おまけ。最初のストレートな C++ 実装を C# で再実装してみました。
ソース
実行結果:

Plain C#   (ST) : 1196.1ms (new!)
Plain C++  (ST) : 140.0ms
SIMD (SSE) (ST) : 46.6ms
SIMD SoA   (ST) : 20.0ms
ISPC       (ST) : 12.3ms

Unity でスクリプトに速度を求めるならプラグイン書いた方がいい、というのがどういうことなのかよくわかると思います。
一応 C# を弁護しておくと、こういう数値計算の類は C# が苦手な処理の最たるものであり、ここまで差が出るケースはレアなはずです。また、C# でもバックエンドが最新の .Net だったりするともうちょっとマシな結果になると思います。
あと、元の C++ 実装を最適化 無効 でビルドすると C# の方が速くなりました。この結果を見ると IL2CPP のアプローチは意外と合理的なんじゃないかと思えてきます。

Compute Shader

おまけ 2。Compute Shader でも実装してみました。
ソース
実行結果:

CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
GPU: NVIDIA GeForce GTX 970

Plain C#   (ST) : 1196.1ms
Plain C++  (ST) : 140.0ms
SIMD (SSE) (ST) : 46.6ms
SIMD SoA   (ST) : 20.0ms
ISPC       (ST) : 12.3ms
Plain C#   (MT) : 224.7ms
Plain C++  (MT) : 33.4ms
SIMD (SSE) (MT) : 12.1ms
SIMD SoA   (MT) : 5.1ms
ISPC       (MT) : 3.2ms
Compute Shader  : 1.45ms (new!)

Compute Shader と比較する場合、CPU 側もマルチスレッド化していないとフェアじゃないのでそうしています。Compute Shader の処理時間は、Dispatch() を呼ぶ直前から計算結果を CPU 側に書き戻すまで、で測定しています。(Dispatch() から制御が戻った時点では実際の処理はまだ終わってないので) なので CPU 側にデータを移す時間分 Compute Shader に不利に働く条件になっているはずですが、これだけの性能が出ます。NVIDIA NSight で純粋な処理時間も見てみたところ、大体 1.1ms くらいという感じでした。
nsight
ゲームでは GPU は大抵レンダリングで忙しいので GPGPU は使いどころが難しいですが、グラフィックはそこそこに留めてこの計算パワーを活かしたゲームデザインを模索してみるのも面白いかもしれません。

Unite 2015 Tokyo

Unity ユーザー向け開発者イベント、Unite 2015 Tokyo で公演しました。上はそのスライドになります。タイトルの通り、基本的にこの blog の過去の検証記事をまとめた内容になっています。

Unite に来る人でこういう内容を求めている人はそんなにいない (シェーダー書く Unity ユーザーの割合はかなり少ない) のがわかっていたので、資料作りにはかなり悩みました。最終的に、あまり深くなり過ぎない程度にレンダリングとアップデート両方を解説した…つもりです。なので、実装について詳しく知りたい場合は元となったこの blog の検証記事とソースを見た方がいいでしょう。(これ これ これ)
少なくともデモはそれなりにインパクトがあったのが見て取れたので、Unity でプログラミングやゲーム開発を始めたような層がこういう領域に興味を持つきっかけになればいいなあ…という感じです。
unity_particle_flood
そんな私の悩みをよそに、本国から来日してきた Unity エンジニアによるセッションでは遠慮なくレンダリング方程式が出てきたりしていたそうで、Unite 終了後脱力していました。まあ、次話す機会があったとしてもそこまでマニアックな話はしないつもりです。

自分的に今回の目玉は MEVIUS FINAL FANTASY の話。奇しくもスクエニから Unity に移った直後にちょっとだけこのプロジェクトに関わる機会に恵まれたのですが、内側と外側両方から見て、改めてスクエニの開発力の高さを実感させられました。公開資料からもその片鱗が伺えるのではないかと思います。
十分な開発力を持ってるところでも、誰が実装しても同じ結果になるような部分にかける手間を省くため、苦手な領域の最低限のクオリティを担保するため、などを目的にゲームエンジンを導入する意義はあると思います。やりたくない部分は既にあるものに任せて、得意な部分により注力する、という考え方です。その得意な部分もゲームエンジンに表現を縛られる必要はなくて、こだわりがある領域や独自性を出したい部分は遠慮無く独自実装していいし、むしろ積極的にやるべきだと思っています。
(私がプライベートでわざわざ Unity 使って Unity っぽくないゲームとかを作ってるのも、それを示したいというのが 1 つのモチベーションとしてあります)

あとは、今回の Unite で 10 年ぶりに専門学校時代の先生に会って「名前と顔で記憶に引っかかって、パーティクルを見て確信した」的なことを言われたのが嬉しいやら感慨深いやらでした。CPU に火を注ぎ続けてもう 10 年経つんだなあと。10 年後の自分も相変わらず CPU に火を注ぎ続けていることを望みます。
IMG_20150412_170553-PANOIMG_20150415_175646IMG_20150415_180708