下のムービーは、512x512=262144個をソートする様子です(2MB)。背景が赤くなると、ソートは完了しています。
もういい加減、順序を気にしなくてよい半透明表示をハードウェアで実装してもらいたいところですが、なかなかそんなものは実現していません。
半透明を正確に実現するためには、GPU によるソートが必要不可欠です。
ソートは、異なる位置のピクセルを入れ替えるという作業で、現在のGPUでは行いにくい方法の１つです。
そのような状況で、Graphics hardware 2003 において、Timothy J. Purcell, Craig Donner, Mike Cammarano, Henrik Wann Jensen, Pat HanrahanによるPhoton Mapping on Programmable Graphics Hardwareが公開されました。
この論文は、フォトンマップをGPUで実装するという野心的な論文です。
この論文で使われている手法の１つに Bitonic sort というソート法があります(ここの説明が詳しいです)。
Bitonic sort はO((logN)2)パス数でソートを実行するソート法です。
1024x1024=1048576 の要素に対して、210パス（ポリゴンを描画するだけなので、最近のGPUによる実行速度1000FPSを考えれば、0.2秒で終了します）。
ここでは、Bitonic sort をDirectX9で実装をしてみます（元の論文にCg版が乗っているので、ほとんど何もやっていないのですが…）。
実は、私もソート法を考えていたのですが、私の考えた方法は N/2-1 の描画回数で終了するもので、今回の Bitonic sort に N=32 を越えた時点で負けるので、おとなしく Bitonic sort を受け入れることにしました。
今回のプログラムは、次のものです。
まぁ、いつものように適当にファイルが入っています。
APP WIZARD から出力されるフレームワークのファイルは紹介を省かせていただきます。
| hlsl.fx | シェーダの入ったエフェクトファイル |
| main.h | アプリケーションのヘッダ |
| main.cpp | アプリケーションのソース |
あと、実行ファイル、プロジェクトファイルが入っています。
Bitonic sort は、最初に隣り合った２つのデータをソートします。
次に、ソートされた２つのブロック（４つのデータ）をソートします。
あとは、ソートする塊をゆきだるま式に大きくして、すべてのデータをソートするまで続けます。
ここで、矢印は、ソートする２つのデータを指し示します。
下を向いた赤い矢印は、小さな値のデータを上に、上を向いた青い矢印は、小さな値のデータを下に動かします。
また、赤いブロックは上のデータほど値が小さくなるようなソートで、青いブロックは下のデータの値が小さくなるような逆順のソートです。
Bitonic sort の特徴として、途中のソートでは、正順と逆順を交互に繰り返す方法であるということがあげられます。
どうして、１つおきに並び替えのデータをかえるのでしょうか？
上のデータを正順、下のデータを逆順に並べておくと次のステップでの並び替えがスムーズに行われます。
たとえば、４つの正順のデータと４つの逆順のデータを組み合わせてソートする時を考えます。
上の図で、小さな値1～4は白い色で、大きな値5～8は灰色にしました。
ここで、４個おきにデータを比較、ソートしてみると、小さな値のデータは全て上に、大きな値のデータは全て下に移動してしまいました。
また、ソートした後も、２個づつデータはソートされるので、（特に上の４つに注目すれば、２つづつ正順、逆順が並んだ全く同じ構造になっています。）
同じステップをあと２回続ければ、データは完全にソートされます。
Bitonic sort のそれぞれのブロックでのソートは、ブロックのデータを2n個として、
最初に2n-1個おきのデータを比較、ソートして、
次に 2n-2個おきの比較をします。
後は、2n-3、2n-4…個おきと、
比較する間隔を数を減らしていけば、
最後にとなりのデータの比較までたどり着いたところでソートが完了します。
それぞれ、右隣に移るのに１回のレンダリングパスで実現することができ、
１つのブロックをソートするのに、n回のレンダリングパスが必要になります。
全体でN個のデータを持つ場合のソートは、それぞれのブロックのソートのパス数を合計して
log2N + log2(N/2) + log2(N/4) + … + 1 = log2N + log2N-1 + log2N-2 + … + 1 = log2N (1+log2N)/2
のパスで完了することになります。
ということで、今回のプログラムを紹介してきましょう。
GPU でデータをソートするので、使えるバッファはテクスチャです。
テクスチャは２次元ですが、テクスチャを左から右に、また上から下の順に番号付けして、
そのインデックスに関して、ソートをします。
１次元から２次元への変換は、次のようになります。
hlsl.fx
0039: // ------------------------------------------------------------
0040: // 1次元のインデックスからテクスチャ座標を算出する
0041: // ------------------------------------------------------------
0042: float2 convert1dto2d(float index)
0043: {
0044: float2 dst;
0045:
0046: dst.x = fmod (index, BufInfo.x) / BufInfo.x;
0047: dst.y = floor(index / BufInfo.x) / BufInfo.y;
0048:
0049: return dst;
0050: }
ここで、BufInfo.x と BufInfo.y は、テクスチャの幅と高さです。
index は、0～BufInfo.x * BufInfo.y - 1の範囲で、２次元のテクスチャ座標の値dstは、0～1の範囲を返します。
ちなみに、この逆に２次元の座標からインデックスへの変換は、次のようになります。
hlsl.fx 0058: // 次のインデックス化のために小数部を切り取る 0059: float2 elem2d = BufInfo * In.WPOS; 0060: elem2d = floor(elem2d); 0061: 0062: // 1次元のIndexを求める 0063: float elem1d = elem2d.y * BufInfo.x + elem2d.x;
In.WPOS は、0～1の範囲のテクスチャの座標、elem1d が１次元のインデックスになります。
Timothy らによるシェーダのHLSL版は、次のようになります。
hlsl.fx
0051: // ------------------------------------------------------------
0052: // Bitonic sort
0053: // ------------------------------------------------------------
0054: float4 BitonicSortPS(VS_OUTPUT In) : COLOR
0055: {
0056: float4 dst;
0057:
0058: // 次のインデックス化のために小数部を切り取る
0059: float2 elem2d = BufInfo * In.WPOS;
0060: elem2d = floor(elem2d);
0061:
0062: // 1次元のIndexを求める
0063: float elem1d = elem2d.y * BufInfo.x + elem2d.x;
0064:
0065: // 上と下のどちらのテクセルと比較するのか？
0066: float csign = (fmod(elem1d, stage)< offset) ? 1 : -1;
0067: // ソートの向き
0068: float cdir = (fmod(floor(elem1d/stepno),2)<=0.5) ? 1 : -1;
0069:
0070: // レンダリング位置のテクセルを読み込む
0071: float4 val0 = tex2D( SrcSamp, In.WPOS );
0072:
0073: // ソート対象のテクセルを読む込む
0074: float adr1d = csign*offset + elem1d;// 比較対照のテクセル
0075: float2 adr2d = convert1dto2d( adr1d );// インデックスからテクスチャ座標に変換
0076: float4 val1 = tex2D( SrcSamp, adr2d );
0077:
0078: // y 成分をソートのキーとして使用
0079: float4 cmin = (val0.y<val1.y) ? val0: val1;// 小さいほう
0080: float4 cmax = (val0.y<val1.y) ? val1: val0;// 大きいほう
0081:
0082: // ソートの向きとデータサンプルの向きを比較して、どちらの値を採用するか決める
0083: dst = (csign==cdir) ? cmin : cmax;
0084:
0085: return dst;
0086: }
最初に１次元のインデックス elem1d を計算します。
次に、比較対照のデータの位置に必要な参照するデータの向きcsignを計算します。csignが1なら、自分より後ろのインデックスのデータと比較して、-1なら、自分よりも前のインデックスのデータと比較します。csign*offset の位置が比較するデータになります。
次に、ソートする向きcdirを計算します。今までに出てきた図の赤い部分は1、青い部分は-1になります。
Timothy らのシェーダでは、この部分は等号で比較していましたが、HLSLではきちんと動かなかったので、不等号を使って、cdirを計算しました。
あとは、それぞれのデータを読み込んで、小さなデータcminと大きいデータcmaxを計算し、ソートの向きとピクセルに格納するデータの位置からどちらのデータを出力するか決定します。
パラメータ stepno, offset, stage は、HLSL内でグローバル変数（定数レジスタ）として扱います。それぞれのパスでのパラメータの値は、次のようになります。
アプリケーション側からパラメータを計算する部分は、次のようになります。
最初に、ソートするブロックの大きさ rank と、ブロック単位でのソートパス数stepを
計算しておき、そこからパラメータの動きをよく考えて各パラメータを作りました。
main.cpp
0352: //-------------------------------------------------
0353: // パラメータの計算
0354: //-------------------------------------------------
0355: int step = m_cnt;
0356: int rank;
0357: for(rank = 0; rank<step ; rank++){
0358: step -= rank+1;
0359: }
0360:
0361: float stepno = (float)(1<<(rank+1));
0362: float offset = (float)(1<<(rank-step));
0363: float stage = 2*offset;
いや、これすげーよ。
速い、速いって！！
210 回レンダリングは、ほとんどの人は非現実的だと思うでしょうが、
1048576個のデータをソートするなら個人的には、これぐらいは全然許容範囲内だと思います。
HWによる正確なアルファレンダリングシステムは、今回の方法と、Per-Pixel Ordering Table(ピクセル単位のｚソート（俺用語）)を用いれば、ps3.0で実現できると考えているのですが、もう２世代先では実現できてほしいなぁ…
ぜんぜん関係ないんだけど、9月号のCG WORLDの倉地紀子さんの記事で bitonic sort が紹介されていたけど、「あれじゃ全然わかんねぇよ！！」と思ったのは、私だけですか…（この説明もわかりやすいかどうか疑問ですが）