Triangle Cullingを実装してみました。

アイデア

出典はOptimizing the Graphics Pipeline with Compute [Graham16]です。
Radeon GPU（PS4/XboxOne世代）のプリミティブ処理能力が低いから、余っているALU性能を活用してプリミティブをカリングし、高速化しようぜという話が書かれています。

今回は、スライドの中盤に書かれているTriangle Cullingのうち、Backface Culling, Frustum Culling, Small Primitive Cullingの3つを実装しました。

ソースコードはこちらあとでアップする

目的

今回はDirect3D 11を使い、ベンダー拡張を含む高速化のアイデア（ballot, MBCNT, Ordered Count, Multi draw indirect, Asynchronous Computeなど）はすべて無視し、カリングができること自体に重きを置きます。

3つのカリング自体は容易に実装できます。
また、プロファイラーにかけ、ラスタライザのカウンターを確認することで、正しくカリングができていることを証明できます。
高速化の効果は測定しません。

実装

面倒なのは、このカリングを仕込むほうです。
カリング結果を格納するIndex BufferとIndirect Draw Argsを新しく用意し、Compute Shaderから適切に出力する必要があります。
後述しますが他にも面倒な点がいくつかあり、アイデアのシンプルさに反して必要なステップ数が多いです。
折衷案として、Geometry Shaderでカリングする方法[Iwasaki18]が提案されています。
残念ながら資料は非公開ですが、このアイデアは容易に実装できます。
双方の利点・欠点を表にしました。

シェーダステージ	CS	GS
容易か？	No（Dispatch Callが増える、レンダリングパイプラインの変更が必要）	Yes（シェーダステージが増える、シグネチャ受け渡しが必要）
追加のメモリ	必要	不要
Draw call	増えない（Multi drawの使用が必須）	増えない
非同期コンピュートへのオフロード	Yes	No
頂点再利用の効率	減る（Radeonならば256単位の分割数であれば影響がない）	？（Geometry Shaderステージがあっても効くか？）
テセレーターとの併用	Yes（Factorに応じた事前仕分け・分割も可能）	Yes
拡張性	Yes	No

評価はしていませんが、前者はAsynchronous Computeが十分高速な環境を要求し、後者はGeometry Shaderのスループットが十分高速な環境を要求しているように思えます。

今回は、最初にGeometry Shaderを使って実装し、後でCompute Shaderを使って同じ結果が得られるものを実装します。

実装

Backface Culling

行列式を使うと三角形の面積を求めることができ、表向きなら正、裏向きなら負、縮退ポリゴンなら0になることを利用しているようです。
どういうわけか、判定条件の正負を逆にしないと正しくカリングされませんでした。
理由は分かっていません。

float3x3 m = { input[0].position.xyw, input[1].position.xyw, input[2].position.xyw };
float d = determinant(m);
if (-d ≦ 0.0) {
	return;
}

実装して気づきましたが、浮動小数点数の誤差があり、ゼロ面積を期待するところで+0.6e7等になってしまい、カリングされないケースがありました。
0.0以外の閾値も検討した方がよいかもしれません。

表か裏かを判定するのならば、外積を使って面法線を求め、射影空間でのZ方向との内積を利用すれば同様に判定できるのではないかと考えました。
試したところ、問題なく動作しました。
座標系は射影空間でもカメラ空間でも問題ないようです（表・裏の関係は変わらない）。

float3 n = cross(positions[2] - positions[0], positions[1] - positions[0]);
float d = dot(float3(0, 0, 1), n);
if (d ≦ 0.0) {
	return;
}

Frustum Culling

ここから先は正規化デバイス座標系(NDC)で計算するので、w除算が必要です。
上に書いたコードを流用します。

Direct3Dでは、NDCでの視錐台の範囲は(-1, -1, 0)から(+1, +1, +1)なので、全ての頂点がこの範囲外ならばカリングできます。

float3 pmin = min(positions[0], min(positions[1], positions[2]));
float3 pmax = max(positions[0], max(positions[1], positions[2]));
if (any(pmax.xy ＜ -1) || any(pmin.xy ＞ 1) || pmax.z ＜ 0 || pmin.z ＞ 1) {
	return;
}

資料では、xyのみチェックしていて、zはチェックしていません。
確かにzがFar Clipを超える状況は限定的でしょうが、Near Clipを超える状況は割とある気がする（部屋や洞窟の壁とか）ので、入れてあげても良いのではないでしょうか。

Small Primitive Culling

NDCでxy座標が[-0.5, +0.5)に収まっていれば、ラスタライザのサンプリングポイントにヒットせず、何も描かれないので、カリングすることができます。
Viewportの大きさに依存するため、定数バッファ等を使ってシェーダの外から情報を与えてあげる必要があります。

float2 vmin = pmin.xy * ViewportScale.xy + ViewportScale.zw;
float2 vmax = pmax.xy * ViewportScale.xy + ViewportScale.zw;
if (any(round(vmin) == round(vmax))) {
	return;
}

シンプルすぎて直感で理解しにくいコードですが、資料の図を見ると簡単に理解できます。
私は図を見たとき目からウロコでした。
これを見て実装しようと思い立ったほどです。
ただし、資料ではfalse negativeが発生するケースが存在することを説明しています。

MSAAを使う場合、特にcentroid指定をしない通常のMSAAでは注意が必要です。
サンプリングポイントが0.5固定ではなくなるため、false positiveが発生します。
資料では、保守的にピクセル中心からサブピクセル位置の差を判定条件として使うよう紹介していますが、これではfalse negativeが大量発生しそうに思えます。
2xと4xに限れば、ピクセル中心を軸としてサブピクセルを回転させれば、xy軸に沿って等間隔に配置させることができるので、厳密に判定できそうですが、明らかに面倒そうなのでパスします。

Geometry Shaderでの実装

Direct3D 11では、シェーダとそのConstant Bufferの設定の追加のみで対応できます。
（正確にはVertex ShaderのSV_Positionセマンティクスを別名に変える必要もあります。）
シグネチャを一致させて値を受け渡すコードを用意するのが面倒ですが、既存のコードを破壊せず、レンダリングパイプラインの変更も不要なので、手軽です。

Compute Shaderでの実装

Compute ShaderではIndex BufferとVertex Bufferを自分で取得し、カリングを通り抜けた三角形のみを新しいIndex Bufferへ書き出します。
プリミティブ速度の改善という当初の目標を達成するため、新しいIndex Bufferはカリングを通り抜けた三角形が隙間なく詰まっている状態にする「コンパクション」が必要があります。
縮退ポリゴン扱いにしてしまっては意味がありません。
今回は、1スレッドが1三角形（3インデックス）を担当するよう実装します。
資料では、コンパクションを256三角形（768インデックス）単位で行っていますが、今回はデバッグの都合で128単位にしました。

頂点の座標値が必要になるので、Vertex Bufferのフォーマットが定義されていなければなりません。
Compute ShaderにInput Layoutという便利なものはないので、自分で与えてやる必要があります。
普通は先頭12byte（offset 0、float3）が座標値なので、今回はその仮定をもとにstrideだけ与えました。
資料では、頂点レイアウトをAoSからSoAに変更しているので、これならoffsetもstrideも自明ですが、アイデアの本質ではないのでAoSのまま使います。

さらに、Draw Indirect Argsも用意します。
これは、コンパクションによりインデックス数が可変になるため、DrawIndexed()ではなくDrawIndexedInstancedIndirect()を使う必要があるからです。
インデックス数はコンパクションの途中で勝手に求まるので、Indirect ArgsのInstance Countを1に固定することを忘れなければ、atomic演算1回の追加で実装できます。
ただし、今回は資料に沿ってより最適なコードを実装します。

groupthread uint triangleCounter;
// ...
if (groupThreadID == 0)
{
	triangleCounter = 0;
}
sync();
// ...
if (triangle is survived) {
	InterlockedAdd(triangleCounter, 1, myThreadID); // 本質的にはこのコードだけで十分
}
sync();
// ...
if (groupThreadID == 0)
{
	IndirectArgs.Store(20 * groupID, triangleCounter * 3);
}

Geometry Shaderでは（既に完了しているため）不要だったWVP変換が必要になります。
2度WVP変換を行うことになり、とても微妙です。
もしスキニングを行っていれば、これも2回必要になり、非常に微妙です。

2byte Index Bufferの扱いは面倒です。
HLSLでは4byte単位でないByteAddressBufferへのオフセットは利用できないので、工夫が必要です。
読み込みの場合、8byteまとめて読んだ後に2byte×3へ分解します。
書き込みの場合、groupsharedメモリにキャッシュし、1スレッドが2スレッド分のデータを出力します（2×2byte×3=12byte）。
お隣のスレッドとのデータ交換なので、資料のDepth Tile Cullingの項で説明している最適化がここでも適用できそうですが、Direct3D 11では使えないので見送ります。

uint indices[3];
uint2 temp;
temp = IndexBuffer.Load2((dtid * 6) & ~3);
indices[0] = (groupThreadID & 1) ? (temp.x ≫ 16) : (temp.x & 0xFFFF);
indices[1] = (groupThreadID & 1) ? (temp.y & 0xFFFF) : (temp.x ≫ 16);
indices[2] = (groupThreadID & 1) ? (temp.y ≫ 16) : (temp.y & 0xFFFF);

float3 vertices[3];
vertices[0] = asfloat(VertexBuffer.Load3(indices[0] * VertexStride));
vertices[1] = asfloat(VertexBuffer.Load3(indices[1] * VertexStride));
vertices[2] = asfloat(VertexBuffer.Load3(indices[2] * VertexStride));

groupshread uint3 indexBufferData[NUM_THREADS];
if (half threads) {
	uint3 index0 = indexBufferData[groupThreadID * 2];
	uint3 index1 = indexBufferData[groupThreadID * 2 + 1];

	uint3 indexPacked;
	indexPacked.x = index0.x | (index0.y ≪ 16);
	indexPacked.y = index0.z | (index1.x ≪ 16);
	indexPacked.z = index1.y | (index1.z ≪ 16);
	IndexBufferCulled.Store3(12 * groupThreadID + 6 * NUM_THREADS * groupID, indexPacked);
}

ちなみに、生き残った三角形が奇数個となった場合、グループの最後のスレッドがカリングされた三角形の最初のインデックス（index1.x）を参照してしまいますが、Index Countの範囲外になり描画には利用されないので、groupsharedメモリは未初期化のまま参照しても問題ありません。
ただ、ちょっと少し気持ち悪いので今回はgroupsharedメモリを初期化しています。

結果

絵が正しいかどうか確認したうえで、プロファイラーを使ってカリング直後の三角形数と実際にラスタライズされた三角形数を比較します。

Backface Culling

球をモデルとして使うと、深度バッファを使わずかつRasterizer StateのCullModeをNoneにしてレンダリングしたとき、もしBackface Cullingの実装が誤っていると表示がダブるので、テストにはちょうど良いです。
凸方でないモデルを使うとこの検証方法は破綻します。

auto rd = CD3D11_RASTERIZER_DESC(CD3D11_DEFAULT());
rd.CullMode = D3D11_CULL_NONE;
hr = mDevice-＞CreateRasterizerState(＆rd, ＆mRasterizerState));

カリング適用前

カリング適用後

実装は問題なさそうです。
ちなみに、Renderdocのフレームキャプチャからバッチを確認すると、カリングされて半分穴抜けになったメッシュが確認できます。
特に、Compute Shader版ではモデルがちょっとずつ描かれていく様子が分かるので面白いです。

では、ラスタライズされた三角形数を確認します。
Queryを使って実装してもよいのですが、面倒なのでRenderdocのプロファイラに数えてもらいます。
GPAなどの詳細な計測ツールを使ってもよいのですが、慣れていないのでRenderdocに頼ります。
ちなみに、Intel HD Graphics 620で動かしています。

計測のためにRasterizer StateのBackface CullingをCULL_BACKに戻します。
通常、1つのメッシュのうち約半数のポリゴンは裏を向いていますから、半分くらいの三角形が消滅するはずです。
なお、上の画像の通りDraw callが3つに分かれているので、プロファイラの結果も3つ出てきます。
脳内で加算してください。

カリング適用前
カリング適用後

予想通り、大体半分になりました。
PS Invokationsが全く同じなので、ズルはしていません。
しかしRasterized Primitivesの数字がBackface Cullingの結果を考慮してくれないのは予想外でした。
本当はちゃんとしたプロファイラを使った方が良いです。

Frustum Culling

さすがに上の方法では全て視錐台に入ってしまいテストにならないので、モデルを3倍スケールしてテストしました。
面倒なので比較画像はなしです。

カリング適用前
カリング適用後

こちらは分かりやすい数字です。
Rasterizer Invocationsの値がVertex Shader（あるいはGeometry Shader）からの出力、Rasterized Primitivesの値がカリングを通過してPixel Shaderに渡った三角形数です。
カリング無効のときは、ハードウェアでカリングされて数が減少しているのがわかります。
一方、カリング有効のときは、ハードウェアでカリングされたものが一つもなく、しかもカリング無効のときと同じ数字なので、期待できる最大の効率が得られています。

3番目のDraw callはtriangle countが0になってしまいました。
このような場合はDraw call自体を止めてしまいたいですが、果たして実装可能なのでしょうか。
NVIDIAの実装では、CPU側で数を指定するので、無理そうです。
AMDの実装では、CPU側で指定するものとBufferから指定するものがあるので、後者を使えば実現できそうです。
Direct3D 12の実装では、Comand Signatureに1つのdraw callを入れておき、ExecuteIndirect()でCounter bufferに数を指定すれば実現できそうです（ただし、Bundleは使えない）。

Small Primitive Culling

画面に極小ポリゴンが生まれる状況にするため、解像度を縦横半分、モデル分割数を約2.5倍にしました。
3倍スケールはしません。

カリング適用前
カリング適用後

カリングはできてはいますが、28しか減りませんでした。
かなり高密度なポリゴンを描画するケースでしか大きな効果はなさそうですね。
そしてDirect3D 11のQueryはここでも役立たないことを知りました。

〆

Trianlge Cullingを実装し、資料に書かれたアルゴリズが問題ないことを確認し、Geometry ShaderとCompute Shaderでのシンプルな実装方法を説明しました。