下の画像をクリックすると、影が動くアニメーションが見られます(注意:3.2MB)。
ということで、PRTです。PRTとは、Peter-Pilke Sloan らがSIGGRAPH 2002 の論文Precomputed Radiance Transfer for Real-Time Rendering in Dynamic, Low-Frequency Lighting Environmentsで提案したグローバルイルミネーションを高速に反映させる手法です。物体表面から反射される光を任意の入射光に関してあらかじめ求めておいて、その結果を球面調和関数の係数として展開することによって、現在のGPUにとってはそれほど負荷を掛けることなくグローバルイルミネーションを実現します。
ただし、ここでは拡散光しか扱いませんし、相互反射もありません。さらに言えば、光源も平行光源なので、非常に手抜きをしたものになっています。
で、いつものようにプログラムです。大事なファイルは次のものです。
ソースには、いつものように適当にファイルが入っています。
| main.h | アプリケーションのヘッダ |
| main.cpp | アプリケーションのソース |
| hlsl.fx | シェーダプログラム |
カーソルキーで、カメラが回って、zやxでズームが変えられます。 また、aキーで「イラディアンスのみ」、「デカールのみ」、「デカール*イラディアンス」の切り替えができます。
なお、実行すると、Athron 2000+ と RADEON 9700 Proの環境で６分ぐらい待たされると思います。これが前計算の放射輝度伝達を計算してる部分です。
実際のアプリケーションを組む場合には、あらかじめ処理しておいてテクスチャとして持つのがよいでしょう。
さて、実際にＰＲＴを実現するにはどうするのかということですが、拡散光に関するレンダリング方程式を思い浮かべてみましょう。
この式は、点xで反射された光L(x)を求めるための式です。
L(ω)は入射光で、V(ω)は、ω方向がさえぎられていたら0、遮られていなかったら1となる遮蔽項です。Nは、xでの法線の向きで、積分はxを取り囲む全天球で行います。ρは拡散反射係数です。
ここで、入射光L(ω)が球面調和関数展開されていた
としましょう。
ちなみに、入射光L(ω)が分かっているときには、球面調和関数展開の係数llmは、
によって求めることができます。
この球面調和関数展開を利用すると、レンダリング方程式は、
のような線型結合で記述することができます。
ここでMlmは、伝達ベクトル（Transfer vector）と呼ばれる量で、
になります。
伝達ベクトルは、遮蔽項と法線ベクトルという物体の幾何学的な情報だけで決まる量なので、あらかじめ計算してしまうことができます。これがPRTの前計算です。
一度この量を計算しておけば、任意のライティングに対して、その球面調和関数展開の係数を求めて線型合成することによって、射出放射輝度を計算することができます。
ちなみに、拡散光でないBRDFの場合には、反射分布が向きに依存するようになりますが、このときには、反射方向に関してさらに球面調和関数展開することによって、任意の角度から見たときの放射輝度を簡単に計算することができます。ただし、伝達ベクトルは、入射方向と反射方向に関する２つの係数を持つようになるので行列（伝達行列：Transfer matrix）に拡張する必要があります。
今回、どのような方法を用いて伝達ベクトルを計算しているのか解説しましょう。
伝達ベクトルは、物体表面の位置に応じて決まるデータ列なので、テクスチャに保存するのが便利です。
テクスチャのR,G,B,A各成分に伝達ベクトルの１つの係数を記録します。より細かく説明すると、球面調和関数の２次元のインデックス(l,m)を「n=l*(l+1) + m」によって１次元のインデックスに変換することができます。この１次元のインデックスの並びを、RGBAの各４成分とテクスチャの番号によって対応させて、テクスチャに情報を記録します。ちなみに、０次から３次までの球面調和関数展開のインデックス数は１６個になるので、４枚のテクスチャで記録することができます。ちなみに、このレベルならば、ピクセルシェーダ1_Xの世代のＧＰＵでもリアルタイムに展開することができます（他のテクスチャは張れないけどね）。
なお、伝達ベクトルの各テクセルの計算は、法線方向をｚ軸にしたサーフェス空間で半球をレンダリングして、それをミップマップの要領で１テクセルまで縮小して張り付けます。
ここで、外側のがくがくしたのが自分自身をレンダリングした効果による遮蔽項で、その中のグラデーションが球面調和関数のそれぞれの(l,m)値に応じた関数値になります。 その向きは、全てのサーフェス表面に関して共通でなくてはならないので、ワールド座標系の値を用いらなくてはなりません（ワールドは動かないので、ここではローカル座標系を使います）。その計算は、ローカル座標系からサーフェス座標系への変換の逆変換をかましてやればよく、サーフェス座標系のテクセル座標から大きさが１になるように３次元ベクトルを作成し、それをローカル-サーフェス座標系の変換行列の逆行列で変換してやればローカル座標系での向きを求めることができます（実際には、法線ベクトルの変換なので、転置行列を使います）。
ちなみに、球面調和関数の値ですが、単位ベクトルを用意しておくと簡単に求めることができます。球面調和関数の３角関数とルジャンドル倍関数の部分は、単位ベクトルを使うと、それぞれの係数の積和で表現することができます。 これらは、PS_2_0レベルで計算することができるので、per-pixelで球面調和関数の値を求めることができます。 残りの規格化の部分は定数としてあらかじめ求めておけばよいです。なお、今回の場合は計算結果が0～1の範囲に収まるような規格化を行っています。
今回の伝達ベクトルのテクスチャを作成するのが、イラディアンスを計算したときと違うのは、遮蔽項を曲面座標系でレンダリングする前のマスクとなるテクスチャを描画する部分に球面調和関数の値を乗算することです（また、今回は球面調和関数の値が負になるので、色の強さが0.5の部分を光が届いていない遮蔽された部分とみなすよう変更しました）。
たとえば、l=2,m=-1の球面調和関数の値を乗算してマスクを作成するピクセルシェーダは次のようになります。
マスクのテクスチャを張る位置のテクスチャ座標からサーフェス座標系に置ける方向Nを求めておき、それをローカル座標系に変換してから、その値をピクセルの見ている向きとして球面調和関数の値を計算します。
なお、使っている球面調和関数の値は黄色の部分で、先ほど紹介した方向ベクトルを使った球面調和関数の値を利用しています。
hlsl.fx
0484: // ------------------------------------------------------------
0485: // SH(2,-1)
0486: // ------------------------------------------------------------
0487: float4 PS_SH21 (VS_OUTPUT_SH In) : COLOR0
0488: {
0489: float4 N = In.Position;
0490:
0491: float l = dot(N,N);
0492: if(1.0<l) N *= rsqrt(l);
0493: N.z = sqrt(1-dot(N,N));
0494: // サーフェス座標系からローカル座標系に変換
0495: N = mul( N, mST );
0496:
0497: float val = -N.y * N.z;
0498:
0499: return 0.5+0.5 * (val * tex2D( SrcSamp, In.Tex ));
0500: }
この計算では、計算した球面調和関数の値の絶対値が1を超えないように規格化定数を調整しています。具体的には、普段つけている規格化定数はつけずに、さらに(2l-1)!!で割った値を使います。ルジャンドル倍関数を評価すると、この値が1を超えないことが分かるので（さらにsinやcosを掛けるが、これらを掛けても絶対値が１以下の数字は１以下である）この規格化定数を用いれば値が飽和して計算がおかしくなることがないと保証できます。
この規格化定数を変えたしわ寄せは入射光を球面調和関数展開するときに現れて、このときには逆に今省いた規格化定数を余計に掛ける処理が必要になります。
なお、できたテクスチャを良く見ると、lの次数が上がるにつれて色の変化が少なくなっているように見えます。これは、「今回の計算では高次になるとアンダーフローが置きやすくなる」とみなす事ができるので、今回の規格化定数も完璧なものではないことが分かります（浮動小数点バッファを使えばいいですが、あれば双線型フィルタリングができないですしねぇ…）。
前計算したラディアンス伝達を反映させるのは、次のような方法で行います。
ベクトルの配列vSHをfloat型の配列とみなして平行光源を球面調和関数で展開した係数を代入します。平行光源は方向に関してはDiracのデルタ関数で処理できるので、球面調和関数のその向きでの値を使えばよいことになります（光の強さや反射係数に関する定数倍は必要ですが）。
ここで、l=2や3の時には球面調和関数の係数を定数倍します。これは、テクスチャに書き込んだ球面調和関数の非対称な規格化定数の部分(2l-1)!!を元に戻すための処理です。
あとはテクスチャを設定した後描画します。
main.cpp
0975: D3DXVECTOR4 vSH[TEX_MAX];
0976:
0977: float ref = 3.0f;
0978:
0979: float phi = atan2f(LightDir.z, LightDir.x);
0980: float theta = acosf(LightDir.y);
0981:
0982: for(int l=0; l<=L_MAX; l++){
0983: for(int m=-l; m<=l; m++){
0984: int n = l*(l+1) + m;
0985: ((float*)&vSH[0])[n] = ref * SphericalHarmonics(l, m, theta, phi);
0986: // テクスチャの展開の補正のため、(2l-1)!! の補正をつける
0987: if(2==l){
0988: ((float*)&vSH[0])[n] *= 3;
0989: }
0990: if(3==l){
0991: ((float*)&vSH[0])[n] *= 15;
0992: }
0993: }
0994: }
0995:
0996: m_pEffect->SetVectorArray("vSH", vSH, (L_MAX+1)*(L_MAX+1)/4);
0997:
0998: // PRT テクスチャの設定
0999: for(int tex = 0; tex < TEX_MAX; tex++){
1000: char str[256];
1001: sprintf(str, "PrtTex%d", tex);
1002: m_pEffect->SetTexture(str, this->m_pFinalTex[tex] );
1003: }
1004: // 背景の描画
1005: pMtrl = m_pMeshBg->m_pMaterials;
1006: for( i=0; i<m_pMeshBg->m_dwNumMaterials; i++ ) {
1007: m_pEffect->SetTexture(m_htSrcTex, m_pMeshBg->m_pTextures[i] );
1008: m_pMeshBg->m_pLocalMesh->DrawSubset( i ); // 描画
1009: pMtrl++;
1010: }
シェーダプログラムでは、0～1の範囲で記録されたPRTテクスチャを-1～1の範囲に復元し、先ほど求めた入射光の球面調和関数に関する係数を掛けて放射ラディアンスを計算します。
あとはデカールを張るなら、デカールの色にラディアンスの強さを掛けたりします。
hlsl.fx
0107: // ------------------------------------------------------------
0108: // PRT テクスチャから入射光を復元する
0109: // ------------------------------------------------------------
0110: float GetPRT(float2 Tex)
0111: {
0112: return dot(2.0*tex2D( PrtSamp0, Tex )-1.0, vSH[0])
0113: + dot(2.0*tex2D( PrtSamp1, Tex )-1.0, vSH[1])
0114: + dot(2.0*tex2D( PrtSamp2, Tex )-1.0, vSH[2])
0115: + dot(2.0*tex2D( PrtSamp3, Tex )-1.0, vSH[3])
0116: ;
0117: }
0118:
0119: // ------------------------------------------------------------
0120: // 全部込みピクセルシェーダプログラム
0121: // ------------------------------------------------------------
0122: float4 PS (VS_OUTPUT In) : COLOR
0123: {
0124: float4 Out;
0125:
0126: // 色
0127: Out = tex2D( SrcSamp, In.Tex0 ) * GetPRT(In.Tex0);
0128:
0129: return Out;
0130: }
やってみて実感したのが、球面調和関数展開では鋭い影を作るのが難しいということでしょうかねぇ。
考えてみれば３次程度では４５度ぐらいの広がりを持った光までしか表現できないので当たり前なのですが、地平線マップと比較してもあまりにもボケボケなので、ちょっと面食らってしまいました。
入射光をHDRにしたり、光沢反射とか、相互反射とかまだまだまだまだやらなきゃいけないことは多いですな。