ちょっと怪しいところがあるので、計算間違いがある可能性を捨てきれませんが、とりあえずシェーダーを実装してみて、ライトマップに焼き込んだ結果を比べることで、計算が正しそうかどうか判断することにします。
これが実装してみたシェーダーコードです。とりあえずテストコードなので最適化とかは考えてません。引数の p はライトローカル座標系での位置、n が法線ベクトル、halfAreaSize がエリアライトの縦横のサイズを２で割ったもので、戻り値が前のページでの積分結果 \(\frac{I}{L_0}\) を返します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 |
float AreaLightDiffuse(float3 p, float3 n, float2 halfAreaSize) { if (p.z < 0) { return 0; } // n.x, n.y が 0 以上になるように座標系を反転させる. float2 s = 2.0f * step(0, n.xy) - 1; p.xy *= s; n.xy *= s; float4 x = -p.xyxy; x.xy -= halfAreaSize; x.zw += halfAreaSize; // 各ケースごとに場合分け. float Dz = p.z * n.z; if (dot(x.zw, n.xy) - Dz <= 0) { // 前のページの最初のケース. return 0; } else { // 点 p を通り法線が n の平面がエリアライトと交差するときに、エリアを左右ではなく上下に分割するように // (前のページのx'_0 = x_0, x'_1 = x_1 のケースになるように) x 軸と y 軸を入れ替える. float2 d = n.x * x.xz + n.y * x.wy; if (d.x < d.y) { x = x.yxwz; n.xy = n.yx; d.xy = d.yx; } float D2 = p.z * p.z; float D11 = sqrt(D2 + dot(x.zw, x.zw)); float2 ln_B = D11 + x.zw; float arctan = -atan2(x.z * x.w, p.z * D11); if (d.x - Dz <= 0) { // 前のページの５番目(最後から２番目)のケース. x.xy = (Dz - n.yx * x.wz)/n.xy; float2 D01 = sqrt(D2 + x.xy*x.xy + x.wz*x.wz); float2 ln_A = D01 + x.xy; float sqnxy = sqrt(dot(n.xy, n.xy)); float lnxy = (sqnxy * D01.y + n.y * x.z - n.x * x.y)/(sqnxy * D01.x + n.y * x.x - n.x * x.w); arctan += atan2(x.x * x.w, p.z * D01.x); float Dx = p.z * n.x; arctan += atan2(Dx + n.z * x.z, n.y * D01.y); arctan -= atan2(Dx + n.z * x.x, n.y * D01.x); return n.x * log(ln_A.y/ln_B.y) + n.y * log(ln_A.x/ln_B.x) + sqnxy * log(lnxy) + n.z * arctan; } else { float D01 = sqrt(D2 + dot(x.xw, x.xw)); float2 ln_A = D01 + x.xw; arctan += atan2(x.x * x.w, p.z * D01); if (d.y - Dz <= 0) { // 前のページの３番目のケース. float invNy = 1.0f/n.y; x.yw = (Dz - n.x * x.xz) * invNy; float2 D00 = sqrt(D2 + x.xz*x.xz + x.yw*x.yw); ln_A.y *= D00.y + x.w; ln_B.y *= D00.x + x.y; float sqnxy = sqrt(dot(n.xy, n.xy)); float lnxy = (sqnxy * D00.y + n.y * x.z - n.x * x.w)/(sqnxy * D00.x + n.y * x.x - n.x * x.y); float Dx = p.z * n.x; arctan += atan2(Dx + n.z * x.z, n.y * D00.y); arctan -= atan2(Dx + n.z * x.x, n.y * D00.x); return n.x * log(ln_A.y/ln_B.y) + n.y * log(ln_A.x/ln_B.x) + sqnxy * log(lnxy) + n.z * arctan; } else { float D10 = sqrt(D2 + dot(x.yz, x.yz)); ln_A *= D10 + x.zy; arctan += atan2(x.z * x.y, p.z * D10); if (dot(x.xy, n.xy) - Dz < 0) { // 前のページの６番目(最後)のケース. x.w = (Dz - n.x * x.x)/n.y; x.z = (Dz - n.y * x.y)/n.x; float2 D00 = sqrt(D2 + x.xy*x.xy + x.wz*x.wz); ln_B *= D00.yx + x.zw; float sqnxy = sqrt(dot(n.xy, n.xy)); float lnxy = (sqnxy * D00.y + n.y * x.z - n.x * x.y)/(sqnxy * D00.x + n.y * x.x - n.x * x.w); float Dx = p.z * n.x; arctan -= atan2(x.z * x.y, p.z * D00.y); arctan += atan2(Dx + n.z * x.z, n.y * D00.y); arctan -= atan2(Dx + n.z * x.x, n.y * D00.x); return n.x * log(ln_A.y/ln_B.y) + n.y * log(ln_A.x/ln_B.x) + sqnxy * log(lnxy) + n.z * arctan; } else { // 前のページの２番目(エリア全体)のケース. float D00 = sqrt(D2 + dot(x.xy, x.xy)); ln_B *= D00 + x.xy; arctan -= atan2(x.x * x.y, p.z * D00); return n.x * log(ln_A.y/ln_B.y) + n.y * log(ln_A.x/ln_B.x) + n.z * arctan; } } } } } |
さて、上記の関数は前のページの \(\frac{I}{L_0}\) を返します。\(L_0\) はエリアライト上の点光源が単位立体角あたりに放出するライトの強さを表していたわけですが、これを Unity のエリアライトが持つ Intensity プロパティから計算しなければなりません。今、点光源が半球上の立体角に一様にライトを放出するとしているので、Intensity \(= 2\pi L_0\) と考えるのが妥当です。なので、上記関数の出力に Intensity\(/2\pi\) を掛けたものをシェーダーの出力とします。
このシェーダーの出力と、Unityでベイクしたエリアライトの結果を比べてみます。
上のシェーダーの全ての分岐がテストされるように、ライトの位置を動かしながら、いくつかスクリーンショットを撮りました。また、log の項と arctan の項がちゃんと入るように、ライトを少し傾けています。
各画像の左側がシェーダーでレンダリングしたもの、右側がライトマップにベイクしたものです。また、Player Settings で Color Space を Linear に設定しています。
なんかそれっぽい結果は得られたものの、全然違う！
これは計算間違いというよりも、そもそものエリアライトのモデルが違うんじゃなかろうか？顕著な違いが見られるのはライトの真横の部分。シェーダーで計算した方は真横にもくっきりライトが当たっちゃってる。
試しに上のシェーダーの出力に p.z*rsqrt(dot(p, p)) をかけてみた結果がこちら。
ちょっと近づきました。やっぱりモデルが違って、エリアライトが点光源の集まりだという前提が間違っていたみたいです。というわけで、これは失敗作。モデルを考え直して次回に続く。
