この記事はレイトレ(Raytracing) Advent Calendar 2024 - Qiita Advent Calendar 2024 - Qiitaの21日目の記事です。リアルタイム向けの面光源の計算方法であるLinearly Transformed Cosine(LTC)という手法について理論から実装までを解説します。

注意事項

今回の解説では明示的に書かない限り論文に従ってz-up座標系を前提にしています。

LTCの理論

LTCとは

Linearly Transformed Cosine(LTC)[Heitz, et al. 2016]とは面光源の証明計算を近似的に行うリアルタイム向けの手法です。代表点サンプリングなどの手法に比べると非常にGround Truthに近く、またTexture Lightにも対応しているのが特徴です。

そもそもの話として面光源の計算は難しいです。なぜなら積分を解く必要があるためです。面光源の立体角を[$ P]とした時、面光源による出射方向 $\omega_v$ への出射光の輝度 $L_v(\omega_v)$ は次のように積分として表現されます。

$\displaystyle{ L_v(\omega_v) = \int_P L(\omega) f(\omega_v, \omega) \cos{\theta} d\omega }$

$f(\omega_v, \omega_i)$ はBSDF、 $L(\omega)$ は $\omega$ 方向から来る入射光を表しています。BSDFは多くの場合複雑な関数であり、面光源がTexture Lightであった場合は $L(\omega)$ もやってみないとわからない関数となるため、一般にこれは解析的に解くことができません。

数値積分的に解ける方法としてはパストレーシングが有名ですが、一発で計算できるというわけではなく何回も計算を行いそれの平均値を求めるという手法です。10回とか20回、下手したら100回の計算を行わないときれいな画像が出ないということもあり、計算量的な面でリアルタイムでは使うことができません(最近はアルゴリズムの発展もあって限定的に使われるようになりましたが)。

そんな面光源ですが、BSDFがLambertかつ、面光源が一様(Texture Lightではない)でPolygonである場合、解析的に解く方法が見つかっています。Lambert反射における面光源の積分はLambert BSDFが

$\displaystyle{ f(\omega_v,\omega) = \frac{1}{\pi} }$

であることを考えれば、次のようになります。

$\displaystyle{ L_v(\omega_v) = L \int_P \frac{1}{\pi} \cos{\theta} d\omega }$

面光源が一様なので定数として $L$ としておけるので積分はかなりシンプルな形状になっていることが分かります。 $N$ 個の頂点を持つ面光源の各頂点を $v_1, v_2, v_3, ...$ とし、各頂点を単位球面上に投影した頂点を $p_1, p_2, p_3,...$ と書くことにすると、この積分は次のような式が成り立つことが知られています(z-up座標系を前提、証明は後述)。

$\displaystyle{ \int_P \cos{\theta} d\omega = \sum_i^N \frac{1}{2} \arccos{(p_i \cdot p_{i+1})} (\frac{p_i \times p_{i+1} }{|p_i \times p_{i+1}|} \cdot (0,0,1) ) }$

これは頂点が少なければ十分リアルタイムに計算することが可能です(Quadぐらいなら全然余裕)。つまるところDiffuse反射なら解析的な手法が既に発見されていたわけです。

しかしながら、一般のPBRマテリアルにおいてはSpecular反射も存在します。Specular反射に使われるMicrofacet BRDF $f_{mf}$ は非常に複雑な式をしており、一般的に使用されるGGX NDFを使ったMicrofacet BRDFは

$\displaystyle{ f_{mf}(\omega_v, \omega) =\frac{D_{ggx}(\omega_m) G(\omega_v,\omega) F(\omega_v, \omega_m, F0)}{4 \cos{\theta_v} \cos{\theta}} }$

$\displaystyle{ D_{ggx}(\omega_m) = \frac{1}{\pi \alpha^2 (\frac{x_m^2}{\alpha^2} + \frac{y_m^2}{\alpha^2} + z_m^2)^2} }$

$\displaystyle{ G(\omega_v, \omega) = \frac{1}{1 + \Lambda(\omega_v) + \Lambda(\omega)}, \Lambda(\omega) = \frac{-1 + \sqrt{1 + \frac{\alpha^2 x^2 + \alpha^2 y^2}{z^2} } }{2} }$

$\displaystyle{ F(\omega_v, \omega_m, F0) = F0 + (1 - F0) (\omega_v \cdot \omega_m)^5 }$

という形をしています。とてもじゃないですがLambertのように解析解が見つかるとは思えません。

そこでLTCの出番というわけです。LTCは線形変換(Linearly Transform)を用いることで、Lambertの積分を任意のBSDFにおける積分に帰着させるという形で積分の計算を近似的に行います。

LTCの概要

LTCではBSDFとコサインを合わせたものを分布(Distribution)と呼びます。ここでは任意の分布は $\overline{D}$ と表現することとします。BSDFは出射ベクトル $\omega_v$ に依存するので、出射ベクトルも明示的に書く場合は $\overline{D_{\omega_v}}(\omega)$ と書くことにしておきます。

$\displaystyle{ \overline{D_{\omega_v}}(\omega) = f(\omega_v,\omega) \cos{\theta} }$

対して、Lambert BSDFにおける分布は特別に $D_o(\omega)$ と表記し、cosine分布と呼びます。

$\displaystyle{ D_o(\omega) = \frac{1}{\pi} \cos{\theta} }$

問題を簡単にするため、一旦ここからの話は面光源は一様であるとして話を進めます(Texture Lightは後述)。この時私たちが求めたい面光源の積分というのは

$\displaystyle{ L \int_P \overline{D_{\omega_v}}(\omega) d\omega }$

という形で書けます。 $L$ は定数なので本質的には

$\displaystyle{ \int_P \overline{D_{\omega_v}}(\omega) d\omega }$

という問題をどう解くかというのが手法の目的となります。

先述したようにLTCはcosine分布の積分で目的の積分を近似することを考えます。具体的にどうやるかというと、線形変換 $M$ を用いてcosine分布 $D_o$ から次のような新しい分布 $D$ を考えます。ここで $\omega_o$ は $\omega_o = M^{-1}\omega$ です。

$\displaystyle{ D(\omega) = D_o(\frac{M^{-1} \omega}{||M^{-1} \omega||}) \frac{|M^{-1}|}{||M^{-1} \omega||^3} = D_o(\omega_o) \frac{d \omega_o}{d \omega} }$

このような関係を持つ新しい分布 $D$ は $D_o$ の積分と等価な積分を作ることが可能という性質を持ちます(解説は後述)。 $D_o$ における積分範囲を $P_o$ と置くと、 $D_o$ と $D$ の積分は次のような関係を持ちます。

$\displaystyle{ \int_{P_o} D_o(\omega_o) d\omega_o = \int_{M^{-1} P_o} D(\omega) d\omega }$

この式の嬉しいところは $M$ によってcosine分布 $D_o$ は別の分布 $D$ に積分の値を変えずに変換できるということです。線形 $M$ は(inverseがある以外)特に条件はないわけですから、作れる分布 $D$ のは数多くあります。それらの $D$ は全てcosine分布として計算することが可能であるということになります。cosine分布から作られる $D$ をLTC(linealy transfomed cosine)分布と呼ぶことにします。

目的の分布 $\overline{D}$ がLTC分布に含まれていればそのままcosine分布にすることが可能ですが、少なくとも線形変換は非線形には対応できないのでLTC分布の範囲には限界があります。そのため、LTCは目的の分布 $\overline{D}$ に対して、それに近いLTC分布 $D$ を用意して近似することを考えます。

$\displaystyle{ \overline{D(\omega)} \approx D(\omega) }$

このように考えれば、対応する線形変換 $M_{\overline{D}}$ さえわかっていれば $D$ を介す形でcosine分布の積分に近似することができます。ここでは変換後の積分領域を $P_o = M_{\overline{D}}^{-1}P$ と置いています。

$\displaystyle{ \int_P \overline{D(\omega)} d\omega \approx \int_P D(\omega) d\omega = \int_{P_o} D_o d\omega_o \tag{3} }$

つまり、まとめれば目的の積分は対応する線形変換 $M_{\overline{D}}$ さえわかっていればcosine分布の積分へと近似することができるということです。LTCはこれによって出てきた近似のcosine分布を計算することで、面光源の計算をできるようにしたという感じです。

ということで、LTCは対象となる分布に対して各パラメーターごとの(実際使うのは逆行列だけなので)線形変換 $M^{-1}$ を事前計算しておき、実際の計算時に(3)を計算します。アルゴリズムの流れとしては以下のような形です。

$\overline{D_{\omega_v}}$ に対する線形変換 $M^{-1}_{\overline{D}}$ を取得する
面光源 $P$ の各頂点を $M^{-1}_{\overline{D}}$ で線形変換して $P_o$ を得る
面光源 $P_o$ でCosine分布の積分を行う

積分と線形変換

アルゴリズムの概要を説明しましたがもう少し理論的な部分もちょっと追っていきましょう。

LTCのモチベーションは何かというと「今考えている積分を同じ計算結果をもつCosine分布の積分にしたい」というものです。しかし、そんな積分が必ず存在するのでしょうか？、あってもどうやって作ればいいかわかるものでしょうか?

高校数学をやった方なら実はこれとほぼ同じ問題をやったことがあります。例えば以下のような積分の問題を考えましょう。

$\displaystyle{ \int_{0}^{3} \sqrt{x + 1} (x+2) dx }$

ぱっと見ただけだと手計算できんのかという気持ちになりますが、ここで置換積分というテクニックを使えば簡単に解くことができます。

置換積分というのは何だったかというと、被積分関数 $f(x)$ に対して、 $x = g(t)$ となるような新しい変数 $t$ を導入した時、成立する積分の関係式

$\displaystyle{ \int_a^b f(x) dx = \int_{\alpha}^{\beta} f(g(t)) \frac{dg(t)}{dt} dt }$

を利用して解く方法です。証明とかはこの辺を見て頂ければと思います。

manabitimes.jp

さっきの例でやれば、 $t$ という変数を次のように定義すれば

$\displaystyle{ x = t^2 - 1 }$

置換積分の公式によって、形が違うものの等価で簡単な積分を導くことができます。

$\displaystyle{ \int_{0}^{3} \sqrt{x + 1} (x+2) dx = \int_{1}^{2} t^3 + t dt }$

左辺は見てわかるように手計算も可能なレベルの簡単な積分です。置換積分を使えばこのように難しい積分も変換次第で解析的に求められる積分へと変形することができます。

LTCで行うことはこれと全く同じことです。変数である入射方向ベクトル $\omega_o$ に対して、線形変換の逆行列 $M$ を使った $\omega$ による置換積分を行っていた、というのがLTCの積分変換の流れです。

$\displaystyle{ \omega = \frac{M \omega_o}{|M \omega_o|} }$

この変換におけるヤコビアン $d\omega / d\omega_o$ は以下のようになります(Appendix参照)。

$\displaystyle{ \frac{d \omega}{d\omega_o} = \frac{|M|}{|M \omega_o|^3} }$

また、置換積分は積分範囲も変換に応じて変えなければいけないのを思い出すと、面光源の立体角(積分領域) $P$ も適切に変えてあげる必要があります。元々 $P$ の範囲にいた方向ベクトルは $M$ をかけられているため、返還後の立体角 $P_o$ は $P$ に対して全体的に $M$ で動かしたものであることが言えるわけです。

例としてPolygon $v_1,v_2,v_3,...$ の面光源を考えますと、変換後の立体角 $P$ は各頂点に $M$ を作用させた $Mv_1, Mv_2, Mv_3, ....$ というPolygonの面光源の立体角になります。ここでは便宜上 $P = M P_o$ と書くことにします。

ここで、 $\omega_o$ から $\omega$ への置換積分を考えてみると先ほどの公式を使えば以下のように書くことができます

$\displaystyle{ \int_{P_o} D_o(\omega_o) d \omega_o = \int_P D_o(\omega_o) \frac{d\omega_o}{d\omega} d\omega }$

この時、変換後の被積分関数Dは以下のような関係があるわけです。これがLTC分布の定義であったわけです。

$\displaystyle{ D(\omega) = D_o(\omega_o) \frac{d\omega_o}{d\omega} }$

もう少し一般的にして、cosine分布から変換できる被積分関数の範囲について考えてみましょう。ある変換 $\omega_o = f_j(\omega_j)$ があった場合、置換積分の公式より

$\displaystyle{ \int D_o(\omega_o) d\omega_o = \int D_o(f_j(\omega_j)) \frac{df_j(\omega_j)}{d\omega_j} d\omega_j }$

と書くことができます。ここで変換後の被積分関数をTransfomed Cosine分布 $\overline{D_j}$ と書くことにします。

$\displaystyle{ \overline{D_j}(\omega_j) = D_o(f_j(\omega_j)) \frac{df_j(\omega_j)}{d\omega_j} }$

変換には様々あると思いますが逆関数 $f_j^{-1}$ があるようなものである場合、 $\omega_j = f^{-1}(\omega)$ と置くことができます(全単射であれば)。そのため、 $D_j(\omega_o)$ は正確にcosine分布に戻すことが可能です。

$\displaystyle{ \int D_j(\omega_j) d\omega_j = \int D_o(\omega_o) d\omega_o }$

逆に言えば、これ以外の関数ではcosine分布に戻すことは(少なくとも)保証はされないはずです。つまりはcosine分布の積分にすることができるのは性質のいい逆関数 $f_j^{-1}$ を持つ変換 $f_j$ による分布 $D_j$ であると考えることができます(多分？)。

何の変換を使うかという話なりますが、少なくとも実用性を考えると以下のような性質を持っていてほしいです。

逆関数が一意的に存在すること
自由度が高い
計算が簡単

線形変換はこれらの性質を満たしています。線形変換 $M$ に対する逆関数は $M^{-1}$ で表せますし、自由度はそこそこ高く(後述)、なおかつ既存のフレームワークで計算することが可能です。なのでLTCでは線形変換を採用しています。

まとめると、LTCではcosine分布を線形変換Mで置換積分したものであり、その変換後の分布はいわば線形変換(Linearly Transformed) cosineと言えるでしょう。これがLTCの名前の由来であり、 $M$ がcosine分布に対する変換であることが分かったと思います。

Fitting

LTCに必要なことは $\overline{D}$ 分布に対する線形変換 $M_{\overline{D}}^{-1}$ の事前計算です。実装上は既にLUTが配布されているので使うだけなんですが、どうやってこの $M_{\overline{D}}^{-1}$ を求めているのかという話も見ていこうと思います。

線形変換 $M$ と行ってきましたが、具体的にはこれは3x3の行列として表現されます。この時、 $M$ の成分とそれによるLTC分布 $D$ の関係性どういう感じか見てみましょう。

LTC分布 $D$ と $M$ の関係性は次のような関係性がありました。

$\displaystyle{ D(\omega) = D_o(\frac{M^{-1} \omega}{||M^{-1} \omega||}) \frac{|M^{-1}|}{||M^{-1} \omega||^3} }$

まずは $M$ が対角行列の場合でx,yにかかる部分をそれぞれ $\alpha_x,\alpha_y$ として、 $D$ がどんな形状になるか上の式で計算してみましょう。

$\displaystyle{ M = \begin{pmatrix} \alpha_x&0&0\\ 0&\alpha_y&0\\ 0&0&1 \end{pmatrix} }$

実際にかかる逆行列 $M^{-1}$ は次のようになります

$\displaystyle{ M^{-1} = \begin{pmatrix} 1 / \alpha_x&0&0\\ 0&1 / \alpha_y&0\\ 0&0&1 \end{pmatrix} }$

これに対するLTC分布は0に近づくほど分布が一つの場所に集まるような形になります。一方で値が大きくなればなるほど分布は外に広がる感じになります。また、 $\alpha_x,\alpha_y$ はそれぞれ動かすことが可能なので縦横で分布を引き延ばすことをこの線形変換で表現できます(Roughness)。

また、 $z$ 成分の方でも同じように変換してみると $\alpha_x = \alpha_y$ と同じように分布を縮める変換をすることができます。

$\displaystyle{ M^{-1} = \begin{pmatrix} 1&0&0\\ 0&1&0\\ 0&0&1/\alpha_z \end{pmatrix} }$

今度は右上と左下に値を入れてみたらどうなるか見てみましょう。これは直感的に理解するのは難しいですが、これらは下の図のように分布を移動するような作用を持ちます(Skewnes)。

$\displaystyle{ M = \begin{pmatrix} 1&0&s_1\\ 0&1&0\\ s_2&0&1 \end{pmatrix} }$

LTCはこれらの操作を用いて分布のFittingを行います。実際に使われる行列はデータ量的な都合や発散防止(Appendix参照)のため、以下のような4つのパラメーターでFittingを行います。これをLUTに保存し、実際に使う際はLUTからこの行列を復元して使用します。

$\displaystyle{ M^{-1} = \begin{pmatrix} a&0&b\\ 0&1&0\\ c&0&d \end{pmatrix} }$

さて、論文のモチベーションはなんだったかというとSpecular反射における面光源の計算をしたいのでした。ということは、今回目的とする分布はSpecular反射の式として一般的に使われている等方性Microfacet BRDF $f_{ms}$ のものとなります。これはRoughnessがxとyで同じMicrofacet BRDFのことを指しており、次のように表されます。

等方性(isotopic)という意味はView DirectionとIncident Directionを方位角方向に回転させてもBRDFが変わらないという意味です。なのでくるくる回転する分には問題がないという性質を持ち、分布としてはView Directionの方位角には依存しません。

フレネルに関してはややこしい問題があるため、Fittingでは取り除いて話を進めていきます(後々ちゃんと扱います)。なので、まとめますとFittingに使用する等方性Microfacet分布はRoughnessとview directionの入射角の２つパラメーターにしか依存しない次のような関数として定義されています。

$\displaystyle{ \overline{D}_{\alpha,\theta_v}(\omega) = \frac{D(\omega_m) G(\omega_v,\omega)}{4 \cos{\theta_v}} }$

これが今回まとめたい目的の分布となります。対応するLTC分布 $D_{\alpha,\theta_v}$ が近似するようにパラメーターの状態ごとにそれぞれ対応した線形変換 $M_{\alpha,\theta_v}^{-1}$ を事前計算する形で、Fittingの処理が行われます。

$\displaystyle{ \overline{D}_{\alpha,\theta_v}(\omega) \approx D_{\alpha,\theta_v}(\omega) = D_o(\omega_o) \frac{d\omega_o}{d\omega}, \omega_o = M_{\alpha,\theta_v}^{-1} \omega }$

各パラメーターの線形変換 $M_{\alpha,\theta_v}^{-1}$ のLUTは以下の画像みたいな感じでHDR画像として保存されています。実装で使うLUTだとパラメーターとの対応として横が $roughness(=\sqrt{\alpha})$ , 縦が $\sqrt{1 - \cos{\theta_v}}$ となっています(ルート取ってるのは精度的な事情だった覚えがあります)。

Limitation

LTCは線形変換で十分表せるような分布という条件下なら基本的にBRDFに問わず使うことができます。なのでやろうと思えばDisney BRDFに対する面光源の計算も理論的には可能なはずです。しかしながら、LTCはLUTの次元が大きくなりやすいという実用上の問題を抱えています。

LTCで使用する線形変換 $M$ は分布 $\overline{D}$ ごとに用意しなくてはなりません。パラメーターが0.1でも違えば当然ながら対応する $M$ も変えなくてはいけなくなります。従って、LUTはパラメーターそれぞれの全部の状態で $M$ を記録したものでないといけなく、その次元数はBRDFのパラメーター数に等しくなります。

等方性Microfacetのパラメーター数はRoughnessと(Viewにも依存するので)Viewの入射角のたった2つなので、LUTは2次元Textureにすることができました。しかしながら、異方性MicrofacetではRoughnessがx,yで別れますし、View Directionについても方位角 $\phi$ も必要になり、計4つのパラメーターを要求します。愚直に考えればこれに対してLTCをするには4次元のLUTを要します。とてもじゃないですがこれは取り扱えません。

なので、LTCは何でも使えるわけではなくパラメーターが少ない、実用的には2つぐらいで済むようなBRDFでなくてはいけないというLimitationを持ちます。等方性 Microfacetが採用されているのはそういうわけです。

ただ、LUTを関数的に近似するorパラメーターを少なくするといった工夫によってLUTのデータ量問題が解決するのであれば異方性 Microfacetとかも使えるようにはなります。実際、異方性Microfacetに対するパラメーター問題を解決(?)してLTCに使えるようにした論文も出ていました(が、そんなに減っている感じでない・・・？)

https://aakashkt.github.io/ltc_anisotropic.html

フレネル

等方性Microfacet分布にはFresnelだけ取り除いた式であったわけですが、なんでこんなことをしたのかというと、これは分布のパラメーターを増やさないようにするためです。

一般的にFresnelの式にはF0値(Specularのカラー)と反射面の法線 $n$ 、入射ベクトル $\omega_i$ によって表現されるShlick近似式が用いられています。

$\displaystyle{ F(F0,\omega_i, n) = F0 + (1.0 - F0) (1.0 - \omega_i \cdot n)^5 }$

この式を見てわかるようにF0というカラーの値が入ってくるため、真面目に分布として考えてしまうとRoughnessなどのパラメーターにF0という新たなパラメーターを追加しなくてはいけません。これはLUTの次元的にも防ぎたいわけです(分布も複雑になりますし)。なのでFittingにおける分布ではフレネルを除いて考えていたわけです。

フレネルを計算しないというのも微妙な話なので考慮したいところですが、愚直な面光源のフレネルの計算はLTCが使えない以上難しそうに見えます。ですが、実は面光源のフレネル問題はEnviroment Lightingにおいて既に取り扱われていました。LTCのフレネル計算はこれをもとに行うことになっています。

Enviroment LightingというのはいわゆるSky Boxとかからの環境マップを光源としたライティングで、理論的にはクソデカTexture Lightからの証明計算という問題になります。

$\displaystyle{ L_v = \int_\Omega L(\omega_i) f(\omega_i,\omega_v) \cos{\theta_i} d\omega_i }$

今までの話で分かるように面光源の計算は非常に難しく、ましてやTexture Lightはまじめにやるのはリアルタイムでは無理です。なので、Enviroment Lightingでは積分をライトとBRDFで分離してしまうというとんでもないSplit Sumという近似を行って何とかします(これは理論とかはなくマジでアドホックな近似だと思います)。

$\displaystyle{ \int_\Omega L(\omega_i) f(\omega_i,\omega_v) \cos{\theta_i} d\omega_i \approx (\int_\Omega L(\omega_i) d\omega_i)(\int_\Omega f(\omega_i,\omega_v) \cos{\theta_i} d\omega_i) }$

意味的にはライトの合計値をBSDFを平均化したファクター的なので明るさを抑えるみたいな感じです。実際計算するときはライトの部分はMipmapなどで近似し、BRDF部分は事前計算から取得して、それを乗算するという形で行われます。正直これでうまく行くんかという気持ちではありますが、見た目的には意外とそこまで変わらないということが知られています。なのでリアルタイムではこの方法が使われているらしいです。

それでBRDFの事前計算はどうするかという話ですが、これがMicrofacet BRDFであればLTCと同様にF0が積分内にあるという問題が生じます。この手法ではフレネルをShlick近似式で展開して、F0を積分の外に出せるように分解してあげるということで解決しています。実際はこの二つの積分を事前計算することで以下の式をリアルタイムに評価できるようにしています。

$\displaystyle{ f(\omega_i,\omega_v) = \frac{F(\omega_i,\omega_m) D_{ggx}(\omega_m) G(\omega_v, \omega_i)}{4\cos{\theta_i}\cos{\theta_v}} }$

$\displaystyle{ \int_\Omega f(\omega_i,\omega_v) \cos{\theta_i} d\omega_i = F0 \int_\Omega \frac{f(\omega_i,\omega_v)}{F(\omega_i,\omega_m)} \cos{\theta_i} d\omega_i + (1.0 - F0) \int_\Omega \frac{f(\omega_i,\omega_v)}{F(\omega_i,\omega_m)} \cos{\theta_i} (1.0 - \omega_v \cdot \omega_m)^5 d\omega_i }$

LTCのフレネルもこれとほぼ同様に処理します。論文とかではちゃんと式がないので怪しいところではありますが、Enviroment Lightingと同様にライトの計算部分とBRDFによる減少量を分離して考えます。(若干両方にBRDFが入ってきているのでええんかという気持ちはありますが、式を見る限りこう解釈する感じなので・・・)

$\displaystyle{ \int_P L f_{mf}(\omega_i,\omega_v) \cos{\theta_i} d\omega_i \approx (\int_{P} L D_{\alpha,\theta_v}(\omega) d\omega_i )(\int_{\Omega} f_{mf}(\omega_i, \omega_v) \cos{\theta_i} d\omega_i ) }$

それでBRDFの減少量については、先ほどと全く同様にF0について展開して積分を二つに分け、それぞれ事前計算できるようにしています。

$\displaystyle{ \begin{align} \int_{\Omega} f_{mf}(\omega_i, \omega_v) \cos{\theta_i} d\omega_i &= F0 \int_{\Omega} \overline{D}_{\alpha,\theta_v} d\omega_i + (1.0 - F0) \int_{\Omega} \overline{D}_{\alpha,\theta_v} (1 - \omega_v \cdot \omega_m)^5 d\omega_i \\ &= F0 n_D + (1.0 - F0) f_D \\ \end{align} }$

実際にLUTとして使われているのは以下のような画像です。Rが積分 $n_D$ ,Bが積分 $f_D$ を表しています(アルファになんか数値入ってるんですが別に使うわけではないみたいです、謎)。

cosine分布の解析解

LTC自体の理論は以上です。ここからはちょっとLTC本体ではないですが根幹となるcosine分布の解析解がいかに出たか簡単に解説します。

そもそもcosine分布における積分というのは思い出すとこんな感じの式でした。

$\displaystyle{ \int_P \cos{\theta} d\omega }$

この式の意味をもう少し考えてみましょう。 $d\omega_i$ というのは方向 $\omega_i$ の小さな立体角、言い換えれば球面上の小さな面積です。これに対して $\cos{\theta_i} d\omega_i$ というのは法線方向に面積 $\omega_i$ を投影した時の面積ということになります。

これを各微小立体角で行って集めるわけですから、その積分値はすなわち立体角 $P$ を法線方向 $n$ への投影面積に等しいことになります。つまるところcosine分布の積分は幾何的に見れば $P$ の投影面積 $P_{proj}$ を求める問題であることが分かります。こうして言われるとなんだかできそうな感じが出てきました。

$\displaystyle{ \int_P \cos{\theta} d\omega = P_{proj} }$

ちゃんとした証明はストークスの定理を使うと求めることができます。気になる人は下のノートを見ていただけると幸いです。

大体このノートに証明が書いてありますが、ある程度補足としてここでも解説をします。

結局のところ問題としては立体角 $P$ の投影面積 $P_{proj}$ の面積を求めるというものでした。一旦シンプルにxy平面だけを考えて、 $P_{proj}$ の面積を求める積分は次のように表現することができます。

$\displaystyle{ \int_{P_proj} dx dy }$

あとの話につなげるため、積分をいったん直交座標ではなく極座標系で考えることにします。

$\displaystyle{ \int_{P_proj} dx dy = \int_{P_proj} r dr d\theta }$

これに対して $P_{proj}$ の境界線を $C_{proj}$ とした時、この積分はグリーンの定理によって $C_{proj}$ の周回積分に変換することが可能です。グリーンの定理は関数 $P,Q$ があった時、領域 $D$ の積分は $D$ の境界線 $C$ の周回積分と次のような関係があるというものです。

$\displaystyle{ \int_D (\frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y}) dx dy = \int_C (P dx + Qdy) }$

$Q = 1/2 r^2$ , $P = 0$ と置くことで、右辺が先ほどの式に当てはまりますので $D_{proj}$ の式は $C_{proj}$ の周回積分へと変えることができます。

$\displaystyle{ \int_{P_{proj}} r dr d\theta = \int_{C_{proj}} 1/2 r^2 d\theta }$

この被積分関数を見てみると丁度半径 $r$ 、弧の長さ $r d\theta$ の扇形の面積に等しい値になっています。すなわちこの積分は $C$ の沿って座標を移動させた時、動径が作る符号付き面積を表すことになります。

とりあえずここで覚えるべきなのは動径の作る面積さえ求められれば良いということです。Polygon Lightの解析解もこれに従って設計されています。では実際にそうなってるのか、再び式を見て考えてみます。

$\displaystyle{ \sum_i^N \frac{1}{2} \arccos{(p_i \cdot p_{i+1})} (\frac{p_i \times p_{i+1} }{|p_i \times p_{i+1}|} \cdot (0,0,1) ) }$

総和の中身が今の頂点 $i$ と次の頂点 $i+1$ についての式であることから、この式は辺に関する値であることが推察できます。

まずは右側の $\arccos$ ですが、 $\arccos$ の中身にある内積は $p_i$ と $p_{i+1}$ のなす角 $\theta_i$ を用いれば $\arccos{(\cos{\phi_i})}$ つまりは $\phi_i$ そのものを求める式であることが分かります。それで $1/2 \phi_i$ という値が出てきますが、これは何かといえば $p_1,p_2$ が作る扇の面積 $S_i$ に相当します。

クロス積の関係より $p_i \times p_{i+1} / |p_i \times p_{i+1}|$ は扇 $S_i$ の法線を表していることから、左側の式は $S_i$ の法線とup-vector(0,0,1)の内積であることが分かります。これを $\cos{\theta_i}$ と表現するとしたら、

$\displaystyle{ \frac{1}{2} \arccos{(p_i \cdot p_{i+1})} (\frac{p_i \times p_{i+1} }{|p_i \times p_{i+1}|} \cdot (0,0,1) ) = S_i \cos{\theta_i} }$

というように扇 $S_i$ をxy平面に投影した面積であることが分かります。図でわかるようにこれは動径の作る面積に相当します。難しそうな式でしたが、こうして目的の面積を求めていたわけです。

周回積分なので向かう方向によって面積の符号が変わらなければなりませんが、これはクロス積の部分で表現されています。 $p_i$ と $p_{i+1}$ が逆転すればクロス積の性質によって符号が反転するので、ありがたいことに結果の面積にはちゃんと適切な符号が付けられることになります。

各辺ごとで総和を取ることで、周回積分の値は計算することができます。これがcosine分布の解析解の意味でした。

ここで一番覚えておいて欲しいのが、cosine分布の積分は本質的に立体角領域の周回積分であるということです。なので、渡す頂点はちゃんとループするようにしなくてはいけなく、外側をぐるっと回るような順番で渡すのが理想です。

一応、トライアングル単位で頂点を計算しても良いですが内部の重複するエッジの分は無駄ですし、回転の向きと合わせる必要があります。その辺に注意して計算していきます。

LTCの実装

それではここからはLTCの実装に入っていこうと思います。実装は[Heitz, et al. 2016]でも公開されていますが、後年にselfshadow氏(恐らく著者関係の人)によって発表された実装が一般的に用いられています。今回はselfshadow氏の実装をベースに解説していきます。

また、実装に関してはselfshadow氏の実装の他、gam0022氏の実装、shivaduke氏のkanikama shaderを参考にさせて頂きました。この場を借りて感謝を申し上げます。

github.com

今回の実装ではUnityを使用して制作しました。以下のリポジトリにUnityのプロジェクトを公開しています。

github.com

Area Lightに関して、理論的には如何なるPolygon Lightでも可能ですが、計算量や実装の制約のため(同一平面上にある)4つの頂点を持つQuad Meshを仮定しています。

Polygon Lightの入力

まずはShaderにAre Lightを結びつけれるようにC#でArea Lightの情報を送れるスクリプトを書くます。最低限Shaderに渡すべきPolygon Lightの情報は以下の3つです。

頂点情報
Emission
(Textureがあるなら) Emission Texture

このうちEmissionとTextureに関してはそのまま渡せばいいですが、頂点情報の転送だけは注意が必要です。先述したようにLTCの積分は境界線の周回積分であるため、外側の境界を作るVertexをループになるよう渡さなくてはいけません。

Quad Meshであればシンプルに各頂点をぐるっと一周するように渡してあげれば問題ありません。C#の実装はこんな感じでworldVerticesで頂点データをもらってshaderに受け渡しています。

        Mesh mesh = GetComponent<MeshFilter>().mesh;
        Vector3[] vertices = mesh.vertices;
        Vector3[] worldVertices = new Vector3[vertices.Length];

        for (int i = 0; i < vertices.Length; i++)
        {
            worldVertices[i] = transform.TransformPoint(vertices[i]);
        }

接空間変換

ここからはShader側の処理に入っていきます。今までの理論は接空間で考えてきましたので、まず最初にやることは送られてきた頂点に対して接空間へと変換する基底変換行列を作ることです。

ちょっとだけ注意なのが接空間の基底がView Vectorに沿うように設定されているという点です。今回対象とするBSDFは等方性(isotropic)という性質があり、方位角方向に全体を回転させても積分の値は変わらないため、x,y基底の取り方は自由です。View Vectorの方位角は考えなくて良いわけなので、統一的にView Vectorに沿うようなx基底を作り、その後に直行するy基底を作るというような形になっています。

実際に作っているのはこの部分です。tangentをx基底、binormalがy基底を表しています。tangentの式はグラムシュミットの正規直行化法というものを使っており、Vに対して同じ方向を向くような基底ベクトルを作ってくれます。binormalはシンプルにクロス積で作っています。

                float3 tangent, binormal;
                tangent = normalize(V - N * dot(V, N));
                binormal = cross(N, tangent);

これで接空間の基底を作ることができました。あとはこれで変換行列を作ればOKです。ちょっとややこしい点としてUnityはy-upなのでz-upの接空間に変換するように列ベクトルの位置が交換されています(なので以降はz-up座標系で考えてもらって大丈夫です)。

Matrixの取得

次はView Vectorの $\cos{theta_v}$ とroughnessの値から対応する線形行列 $M_{\alpha,\theta_v}^{-1}$ をLUTから取得します。selfshadow氏が提供しているLUTは64×64の4channel float Textureで、横が $roughness$ , 横が $\sqrt{1-\cos{theta_v}}$ に対応しています。

補間を前提としているため、texelの中心位置とパラメーターが合うようにちょっとだけuvに調整が必要です。この処理を踏まえてLUTのuvを取得する関数は次のように書くことができます。

static const float LUT_SIZE = 64;
static const float LUT_SCALE = (LUT_SIZE - 1.0) / LUT_SIZE;
static const float LUT_BIAS = 0.5 / LUT_SIZE;

float2 GetLUT_Texcoord(float cosine, float roughness)
{
    float2 uv;
    uv.x = roughness;
    uv.y = sqrt(1.0 - cosine);

    uv = uv * LUT_SCALE + LUT_BIAS;
    return uv;
}

このUVを使ってLUTから値を取得して、Matrixを生成します。

float3x3 GetLTC_InverseMatrix(float2 uv, sampler2D lutSampler)
{
    float4 lut = tex2D(lutSampler, uv);
    return float3x3(
        float3(lut.x, 0, lut.y),
        float3(0, 1, 0),
        float3(lut.z, 0, lut.w)
    );
}

Unity側の補足としてFloat Textureを扱う際はフォーマットにRGBA Floatを指定しておく必要があります。もしバグっぽい挙動をしていたらこの辺を見るといいかもしれません。

これで座標変換が整いました。実装では接空間の変換行列とLTCのMatrixを合わせて、Area Lightの座標変換をしています。

                Minv = mul(transpose(float3x3(tangent, binormal, N)), Minv);

                float3 v[5];
                v[0] = mul(points[0] - P, Minv);
                v[1] = mul(points[1] - P, Minv);
                v[2] = mul(points[2] - P, Minv);
                v[3] = mul(points[3] - P, Minv);
                v[4] = 0.0;

Clipping

これで座標変換の方に関する準備は終わりました。ここから実際にPolygon Lightに対する積分値を求める処理に入ります。以上の行列をかけてLTCにおけるPolygon Light $P_o$ を作ってもいいのですが、単にやるだけでは実は問題があります。

それは $P_o$ がshading点より下 $z < 0$ に行ってしまった場合に起きます。そもそもの前提として $z > 0$ の領域が積分対象であるため、それ以外の領域は計算の対象ではありません(=照らされることはない)。なのではみ出た部分のエッジを計算しないようにClippingという下処理を行います。

Clipping自体はシンプルで $z = 0$ で $P_o$ をバッサリ切り落として、新しい頂点を作るというものです。これによってQuadのArea LightはClipping後,最大5個の頂点を持つ可能性が存在します。

// Z-clip
void ClipQuadToHorizonNormalize(inout float3 L[5], out int n)
{
    // detect clipping config
    int config = 0;
    if (L[0].z > 0.0) config += 1;
    if (L[1].z > 0.0) config += 2;
    if (L[2].z > 0.0) config += 4;
    if (L[3].z > 0.0) config += 8;

    // clip
    n = 0;

    if (config == 0)
    {
        // clip all

    }
    else if (config == 1) // V1 clip V2 V3 V4

    {
        n = 3;
        L[1] = GetClipPoint(L[0], L[1]);
        L[2] = GetClipPoint(L[0], L[3]);
    }
    else if (config == 2) // V2 clip V1 V3 V4

    {
        n = 3;
        L[0] = GetClipPoint(L[1], L[0]);
        L[2] = GetClipPoint(L[1], L[2]);
    }
    else if (config == 3) // V1 V2 clip V3 V4

    {
        n = 4;
        L[2] = GetClipPoint(L[1], L[2]);
        L[3] = GetClipPoint(L[0], L[3]);
    }

  .....

Clippingを実装している部分は非常に長いですが、やっていることは単純で $z < 0$ の頂点を検出して、シンプルに全通り場合分けした実装をしています。新しく作る頂点を作る処理はエッジに対して $z = 0$ となるような内分点を求め、ループになるように頂点を格納するという内容になっています(後々正規化するので内分点の割り算は省略されています)。

これを正規化すればようやく解析解に使用する頂点 $p_0, p_1, ...$ が求められます。

                float3 p[5];
                p[0] = v[0];
                p[1] = v[1];
                p[2] = v[2];
                p[3] = v[3];
                p[4] = 0.0;


                int numPoint;
                ClipQuadToHorizon(p, numPoint);

                if (numPoint == 0)
                    return float3(0, 0, 0);

                p[0] = normalize(p[0]);
                p[1] = normalize(p[1]);
                p[2] = normalize(p[2]);
                p[3] = normalize(p[3]);
                p[4] = normalize(p[4]);

面光源の計算

面光源の計算は各エッジの回数分式(1.1)を計算すれば良いため、エッジの評価をする関数EvaluateEdgeを追加して計算すればOKです。

                float integration = 0.0;
                integration += IntegrateEdge(p[0], p[1]);
                integration += IntegrateEdge(p[1], p[2]);
                integration += IntegrateEdge(p[2], p[3]);
                if (numPoint >= 4)
                    integration += IntegrateEdge(p[3], p[4]);
                if (numPoint == 5)
                    integration += IntegrateEdge(p[4], p[0]);

                integration = abs(integration);

それでEvaluateEdgeの中身ですが、愚直に式(1.1)に基づいて書けば以下のように書くことができます。ちょっとわかりづらいですが、 $|v_i \times v_{i+1} | = \sin{\theta_i}$ ということを思い出すと、ちゃんとこの式になります。これをやると結構Ground Truthとあまり合わないという問題があります(あと私の環境ではわからなかったですが、バンディングが生じることも)。

            float IntegrateEdge(float3 v1, float3 v2)
            {
                float cosTheta = dot(v1, v2);
                float theta = acos(cosTheta);
                float res = cross(v1, v2).z * theta / sin(theta);
                
                return res;
            }

なぜ、これが生じるのかというとGPUで使用されるacos関数の精度が足りないためです。selfshadow氏の解説によるとGPUにおけるacos関数は内部的に3次関数のFittingによって実装されているそうです。

            float IntegrateEdge(float3 v1, float3 v2)
            {
                float x = dot(v1, v2);
                float y = abs(x);

                float a = 0.8543985 + (0.4965155 + 0.0145206 * y) * y;
                float b = 3.4175940 + (4.1616724 + y) * y;
                float v = a / b;

                float theta_sintheta = (x > 0.0) ? v : 0.5 * rsqrt(max(1.0 - x * x, 1e-7)) - v;

                return (cross(v1, v2) * theta_sintheta).z;

多くの場合では問題になりませんが端の部分では精度が落ちることが知られており、それが原因で今回の問題が起きています。なので対策としてselfshadow氏は $\theta / \sin{\theta}$ に対するフィッティング関数を提案しています。

            float IntegrateEdge(float3 v1, float3 v2)
            {
                float x = dot(v1, v2);
                float y = abs(x);

                float a = 0.8543985 + (0.4965155 + 0.0145206 * y) * y;
                float b = 3.4175940 + (4.1616724 + y) * y;
                float v = a / b;

                float theta_sintheta = (x > 0.0) ? v : 0.5 * rsqrt(max(1.0 - x * x, 1e-7)) - v;

                return (cross(v1, v2) * theta_sintheta).z;
            }

これによる結果はかなりGround Truthに近い感じになります。実装ではこちらを使います。

フレネル値の計算

フレネルに関しては先述したようにSplitSumによって計算します。LUTから同様に値を取得し、 $n_D$ , $f_D$ を取得します。あとはフレネルの節で説明したように

$\displaystyle{ F0 n_D + (1 - F0) f_D }$

を使って計算して積分値に対してかけてあげれば終わりです。

                float3 specular = EvaluateAreaLight(N, V, worldPos, Minv, L);

                float4 fresnelLut = tex2D(_Fresnel_LUT, lutUV);
                specular *= F0 * fresnelLut.r + (1.0 - F0) * fresnelLut.g;

まとめ

以上でLTCの実装ができました。LTCによる結果はパストレによるGroundTruthと比べてもほとんどわからないほど正確に出ています(若干明るいぐらい)。ゲーム用途なら十分な精度であることが分かります。

また、LTCのMatrixをIdentityにすればLambertとしてそのまま使えるのでDiffuse反射も同じフレームワークで計算することができます。

float3 diffuse = EvaluateAreaLight(N, V, worldPos, float3x3(1, 0, 0, 0, 1, 0, 0, 0, 1), L);

BaseColor, Metallic, Roughness Textureに対応すればかなり良い見た目の結果が得られます。これがリアルタイムに動かして見れるのでなかなか面白いです。

Texture Light

今までは面光源はどの位置でも輝度が変わらない一様光源を仮定してきましたが、次は輝度がTextureによって表現されるTexture Lightの場合について考えていきます。

今までは一様光源を改定していたおかげで輝度の関数を積分の外へと分離することができましたが、Texture Lightでは輝度が方向依存性を持つのでそれをすることはできません。なので証明計算の式は以下のようになります。

$\displaystyle{ L_v(\omega_v) = \int_P L(\omega) \overline{D}_{\omega_v}(\omega) d\omega }$

$L(\omega)$ は実際どんな値になるかはやってみないとわからないので、これを解析的に解くのは不可能です。なのでLTCではアドホックにライトの輝度を分離してガウシアンブラーでぼかした色 $L_{blur}$ をライトの光源として扱う形でなんとかそれっぽく見せるという手法を取っています。

これの理論的な背景を追ってみましょう。まず、Texture Lightの積分について通常通りLTCによる近似を考えます。

$\displaystyle{ \int_P L(\omega) \overline{D}_{\omega_v}(\omega) d\omega \approx \int_P L(\omega) D_{\omega_v}(\omega) d\omega }$

ちょっと急ですが次のような分離を考えます。

$\displaystyle{ \int_P L(\omega) D_{\omega_v}(\omega) d\omega = \frac{\int_P L(\omega) D_{\omega_v}(\omega) d\omega}{\int_P D_{\omega_v}(\omega) d\omega} \int_P D_{\omega_v}(\omega) d\omega }$

唐突なのでびっくりするかもしれませんが、これは単に分布だけの積分値を追加しただけで約分すれば何も変わっていません。右側の積分値については先ほどの一様光源でやったようにLTCによって計算することが可能です。

逆に右側の積分についてもう少し見てみましょう。分母の積分値は単なる定数なので分子の値の方に入れてみるとこんな形にすることができます

$\displaystyle{ \frac{\int_P L(\omega) D_{\omega_v}(\omega) d\omega}{\int_P D_{\omega_v}(\omega) d\omega} = \int_P L(\omega) \frac{D_{\omega_v}(\omega)}{\int_P D_{\omega_v}(\omega) d\omega} d\omega }$

この時、分布に関わる項はそれだけ積分すると(当たり前な気がしますが)1になるという性質を持ちます。つまり、正規化された関数であるわけです。

$\displaystyle{ \int_P \frac{D_{\omega_v}(\omega)}{\int_P D_{\omega_v}(\omega) d\omega} d\omega = \frac{\int_P D_{\omega_v}(\omega) d\omega}{\int_P D_{\omega_v}(\omega) d\omega} = 1 }$

正規化しているという性質から、これはある種のフィルターと見なすことも可能です。そのため、これを $Filter$ として考えると、

$\displaystyle{ Filter(\omega) = \frac{D_{\omega_v}(\omega)}{\int_P D_{\omega_v}(\omega) d\omega} }$

輝度 $L$ の積分はTextureのフィルタリング(畳み込み)をしているものと見なすことができます。

$\displaystyle{ \int_P L(\omega) \frac{D_{\omega_v}(\omega)}{\int_P D_{\omega_v}(\omega) d\omega} d\omega = \int_P L(\omega) Filter(\omega) d\omega =: L_{Filter} }$

なのでこの結果を $L_{Filter}$ と書けばPolygon LightはTextureのフィルター結果と分布の積分の掛け算という形で近似することができます。これがLTCにおけるTexture Lightの計算式です。

$\displaystyle{ L_v(\omega_v) \approx L_{Filter} \int_P D_{\omega_v}(\omega) \omega }$

つまるところTextureのFilteringさえ追加すればよいわけですが、結構複雑な式ですしこのFilterをまじめにやるのは厳しいところです(あと地味にTexture Spaceではないのでめんどい)。なので、LTCではもうFilterはTexture Spaceのガウシアンフィルターで何とかするという形になっています。

Filteringの仕方について

Filteringされた輝度 $L_{Filter}$ についてですが、元の式を考えるとこれはcosine分布に変換しても値は変化しません。

$\displaystyle{ L_{Filter} = \frac{\int_P L(\omega) D_{\omega_v}(\omega) d\omega}{\int_P D_{\omega_v}(\omega) d\omega} = \frac{\int_{P_o} L_o(\omega_o) D_o(\omega_o) d\omega_o}{\int_{P_o} D_o(\omega_o) d\omega_o} }$

ここでの $L_o(\omega_o)$ は $P_o$ の輝度分布を示しています。この式は $L_{Filter}$ は線形変換に対して不変であり、cosine分布でのTexture Lightの問題へと帰着できるということを示しています。なので以降はcosine分布での状況を考えていきます。

ガウシアンブラーをかけるでしたが、これにはブラーの中心位置と標準偏差(=ブラーの強さ)を指定しなくてはいけません。中心位置についてはぱっと考える限りは分布の強い方向の先のtexelを選択すればいいのではと思いつきます。ですが、確実にそちらの方向にArea Lightがある保証がないですし、Area Light面の傾きによってはフィルタリングの異方性も考えなくてはいけません。

なのでLTCでは確実性を重視して、Area Lightに対して平行な無限大のArea Light面を考えて、ライティング点をArea Light面に垂直に投影したtexelをブラーの中心とします。無限遠の平面を導入しましたが、この平面はArea Lightの領域だけTextureが張られていて、外側は黒というのを考えています。これならどのような状況下でも中心となるtexelを用意することができます。

こうしてとりあえず中心位置は決定できました。次はガウシアンの標準偏差 $\sigma$ の設定です。これは小さいほどブラーが弱く、大きいほどブラーが強くなるという関係があるのでボケるような状況下で値が大きくなるように設計してあげる必要があります。

LTCでは $\sigma$ についてArea Lightの面積 $A$ と距離 $r$ を用いて次のように定義しています。

$\displaystyle{ \sigma = \sqrt{\frac{r^2}{2 A}} }$

これはつまり、「Area Lightが大きければ大きいほど明瞭に、小さければ小さいほどボケる」「距離が離れるほどボケて、近づけば近づくほど明瞭になる」ということを表しています。なんとなく直感と会うような感じがあると思います(立体角的な感じ)。これは経験則的な面もあると思うのであまり気にしない方がいいと思います。

また、もう一つ重要な点としてArea Light外の領域のブラーは内部のTextureを含むように拡張するようにします。なのでFilteringした後のArea Light面はどこでも一応値が入っているような感じになります(ある種padding的なやつですね)。下の図では黒丸がブラーの半径を示しており、外部に行くほど半径を広くして外側を埋めていることがわかると思います。

[Heitz 2016]より引用、Textureの外の部分はTexture内部に触れるようにBlurを拡大する

なので実際のブラー半径はTexture Space上のArea Lightまでの距離 $r_{tex}$ を使うと次のようになります

$\displaystyle{ \sigma = \sqrt{\frac{r^2}{2 A}} + \max(0,r_{tex}) }$

実際、このpaddingの有無で結構輪郭が変わってきます。あった方がなんとなく自然な感じになります(根拠はあるんでしょうか？多分ヒューリスティックな気がする)。

これを用いてガウシアンブラーを行えばTextureのフィルタリングは終わりです。

実装

ではこれを実装していきましょう。Textureの実装は元論文にはないため結構アプリケーションによってまちまちで、今回はshivaduke氏のkanikama shaderの実装を元に解説させて頂きます(掲載許可を出していただき誠にありがとうございます)。いくつかわかりやすさのため、実装が異なる部分がありますが根本的にはkanikama shaderの実装に従っています。

github.com

実装的にはTextureのFilteringする関数を追加するだけです。LTC本体にはこの関数を呼び出し、最終出力に乗算するだけで良いです。

                Lo_i *= FilterTexture(_AreaLightTexture, v[0], v[1], v[2]);

Filteringする関数FilterTextureの中身はこんな感じです。

            float GaussianKernel(in float x, in float sigma)
            {
                float s = 1 / sigma;
                // 1/sqrt(2 * PI) = 0.39894
                return 0.39894 * exp(-0.5 * x * x * s * s) * s;
            }

            float GaussianInv(float y, float sigma)
            {
                // sqrt(2 * PI) = 2.50662
                return sigma * sqrt(-2 * log(2.50662 * sigma * y));
            }

            float SquareSDF(float2 uv)
            {
                uv -= 0.5;
                float2 st = abs(uv) - 0.5;
                return max(st.x, st.y);
            }

            float3 FilterTexture(sampler2D tex, float3 v0, float3 v1, float3 v2)
            {
                // Reference : Kanikama shader by shivaduke28
                // https://github.com/shivaduke28/kanikama

                // Orthogonal Orojection
                float3 V1 = v0 - v1;
                float3 V2 = v2 - v1;
                float Area = length(cross(V1, V2));
                float3 N = normalize(cross(V1, V2));
                
                float r = dot(v1, N);
                float3 P = r * N - v1;

                // Skew coordinates
                float dotP1V1 = dot(P, V1);
                float dotP1V2 = dot(P, V2);
                float dotV1V1 = dot(V1, V1);
                float dotV2V2 = dot(V2, V2);
                float dotV1V2 = dot(V1, V2);
                float delta = dotV1V1 * dotV2V2 - dotV1V2 * dotV1V2;

                float2 uv;
                uv.x = (dotV2V2 * dotP1V1 - dotV1V2 * dotP1V2) / delta;
                uv.y = (-dotV1V2 * dotP1V1 + dotV1V1 * dotP1V2) / delta;

                // Blur sigma
                float sigma = abs(r) / sqrt(Area);
                float add = max(0, SquareSDF(uv));
                sigma += add;

                // Approximate Gaussian function by step functions.
                // Texture's Filter Mode should be Trilinear.
                float y0 = GaussianKernel(0, sigma);
                float y1 = y0 * 0.75;
                float x1 = GaussianInv(y1, sigma);
                float y2 = y0 * 0.5;
                float x2 = GaussianInv(y2, sigma);
                float y3 = y0 * 0.25;
                float x3 = GaussianInv(y3, sigma);

                half4 col = 0;

                float2 dx = float2(0.5, 0);
                float2 dy = float2(0, 0.5);

                col += tex2Dgrad(tex, uv, dx * x3, dy * x3) * 0.333;
                col += tex2Dgrad(tex, uv, dx * x2, dy * x2) * 0.333;
                col += tex2Dgrad(tex, uv, dx * x1, dy * x1) * 0.333;
                return col;
            }

一つ一つ見ていきましょう。まずArea Light側の入力として線形変換後の頂点 $v_0, v_1, v_2$ をもらいます。これらがArea Lightを張っているとして、Area Lightの辺としてベクトル $V_1, V_2$ を以下のように生成します。位置関係としては図のような感じです

ここからArea Lightの情報を作ります。まず最初では法線 $N$ と面積 $A$ を求めています。

                float3 V1 = v0 - v1;
                float3 V2 = v2 - v1;
                float Area = length(cross(V1, V2));
                float3 N = normalize(cross(V1, V2));

その次はshading点 $o$ をArea Light面にOrthogonal Projectionし、Textureの中心となる点 $p$ を計算します。今回、 $o$ は接空間上では原点になるようにしているため(0,0,0)として、ポイント $p_1$ を通る法線 $N$ のArea Light面にそれを投影するみたいなことを考えます。

やり方は単純で $v_1 - o$ と $N$ の内積を取れば、その値というのは面から $o$ までの距離が得られるので、その分だけ $o$ から $N$ の方向に下すという感じです(符号とかは内積の符号がそのまま使える)。

$\displaystyle{ p = ((v1-o) \cdot N) N }$

実装上では後々 $v_1$ を中心に考えたいので $v_1$ の相対ベクトルとして作っておいています。ついでにここで $r$ も計算できます。

                float r = dot(v1, N);
                float3 P = r * N - v1;

これでtextureの中心点 $p$ が求まりました。次は $p$ のTexture UVを求めるわけですが、ちょっと面倒な話があります。 $V_1$ と $V_2$ に対する内積取ればええやんと思うんですが、線形変換後のQuadは一般的に歪ませられており、 $V_1, V_2$ は直交しているとは限りません。なのでゆがみに対して対応したUVを考えなくてはいけません。

これは基底が直行していない $V_1,V_2$ が張る斜交座標を求める問題として見なすことができます。UVを $u,v$ それぞれとした時、位置 $p$ とは

$\displaystyle{ p = u V_1 + v V_2 }$

という関係にあります。このUVを求めるという問題を考えます。

V1,V2で内積を取れば連立一次方程式が得られるので、

$\displaystyle{ p \cdot V_1 = u |V_1|^2 + v V_1 \cdot V_2 }$

$\displaystyle{ p \cdot V_2 = u V_1 \cdot V_2 + v |V_2|^2 }$

逆行列を使えば、 $u,v$ については次のように求められます。

$\displaystyle{ \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} |V_2|^2 & - V_1 \cdot V_2 \\ - V_1 \cdot V_2 & |V_1|^2 \end{pmatrix} \begin{pmatrix} p \cdot V_1 \\ p \cdot V_2 \end{pmatrix} }$

これを愚直に書くと以下のようになります。これでTexture UVを求めることができました(kanikama shaderではもうちょい違う方法でやっていらっしゃるみたいです)。ちなみに頂点の送り方によってUVは変わるのでその点はuvを反転させたりして合わせたりしてください。

                float dotP1V1 = dot(P, V1);
                float dotP1V2 = dot(P, V2);
                float dotV1V1 = dot(V1, V1);
                float dotV2V2 = dot(V2, V2);
                float dotV1V2 = dot(V1, V2);
                float delta = dotV1V1 * dotV2V2 - dotV1V2 * dotV1V2;

                float2 uv;
                uv.x = (dotV2V2 * dotP1V1 - dotV1V2 * dotP1V2) / delta;
                uv.y = (-dotV1V2 * dotP1V1 + dotV1V1 * dotP1V2) / delta;

次はBlurの $\sigma$ を求めていきます。 $A, r$ については既に求めらているので問題なく計算できます。ただし、Texture UVが外(uv < 0 or 1 < uv)にある場合はその分だけ $sigma$ を伸ばさなくてはいけないという話があったので、それを実現するためBoxのSDFを用いてTexture Space上でのArea Lightまでの距離を求めています。この分を加算して $\sigma$ は完成です。

            float SquareSDF(float2 uv)
            {
                uv -= 0.5;
                float2 st = abs(uv) - 0.5;
                return max(st.x, st.y);
            }

                // Blur sigma
                float sigma = abs(r) / sqrt(Area);
                float add = max(0, SquareSDF(uv));
                sigma += add;

これでようやくガウシアンブラーの準備が整いました。後はブラーの処理を書けばいいのですが、少し問題があります。それはガウシアンブラー自体が結構重いというものです。ShaderにおいてTexture Fetchは重い処理に分類されるもので、1Passガウシアンブラーは一般的にはやるべきではありません(今のGPUならそこまで神経質にならなくてもいいと思いますが)。

この辺はかなり実装が分かれるところで、現論文だとArea LightのTextureは $sigma$ ごと事前にいくつかフィルタリングしたものを使用するという感じでランタイムで行わないことを述べています。ですが、やはり自由度的にも事前フィルターは避けたいところです。kanikama shaderではmipmapを使ってガウシアンフィルターをstep functionとして分解して近似する手法を使っています。

この手法ではガウシアンについて、図のように $y_1, y_2, y_3$ のように $y$ の値で階段状に区分して考えます。これに対応する幅 $x_1, x_2, x_3$ を求め、それに対応するようなmipmapを取得し総和するという形でガウシアンフィルターを近似させます(詳しい理論はちょっと私にはわからないです)。平たく言えば、ガウシアンブラーっぽくなるmipmapの取り方をしているというような感じです。

実装は以下のようになっています。今回は $y_1, y_2, y_3$ は4分割を考え、それぞれ対応する $x_1,x_2,x_3$ を逆関数から求め、最終的に対応するmipmapをtex2Dgradでもらってきて3つのmipmapを合成しています。

            float GaussianKernel(in float x, in float sigma)
            {
                float s = 1 / sigma;
                // 1/sqrt(2 * PI) = 0.39894
                return 0.39894 * exp(-0.5 * x * x * s * s) * s;
            }

            float GaussianInv(float y, float sigma)
            {
                // sqrt(2 * PI) = 2.50662
                return sigma * sqrt(-2 * log(2.50662 * sigma * y));
            }

                float y0 = GaussianKernel(0, sigma);
                float y1 = y0 * 0.75;
                float x1 = GaussianInv(y1, sigma);
                float y2 = y0 * 0.5;
                float x2 = GaussianInv(y2, sigma);
                float y3 = y0 * 0.25;
                float x3 = GaussianInv(y3, sigma);

                half4 col = 0;

                float2 dx = float2(0.5, 0);
                float2 dy = float2(0, 0.5);

                col += tex2Dgrad(tex, uv, dx * x3, dy * x3) * 0.333;
                col += tex2Dgrad(tex, uv, dx * x2, dy * x2) * 0.333;
                col += tex2Dgrad(tex, uv, dx * x1, dy * x1) * 0.333;
                return col;

まとめ

こうしてTextureの実装が終わりました。ようやくこれで実装は終わりとなります。お疲れさまでした。

GroundTruthと比較してみるとやはりガウシアンブラーの近似などの影響でライト近辺では若干色合いなどが違いますが(ガウシアンを使うこと自体近似だし)、それでもリアルタイムでは十分なクオリティを持っていることが分かります。

Diffuse光を含めてみると正直これだけ見せられても正解と思うぐらいにはかなりきれいに出てくれます。これがぐりぐり動かせるのがなかなか面白いですね。

終わりに

ここまで見てくださってありがとうございます。

私がLTCを知ったのはVRChatで、クラブワールドなんかでよく画面の反射に使われていると聞いてからでした。見てみるとキレイなもんで仕組みをちゃんと知りたいなと思っていたのですが、一見するとわかりづらい理論(&謎の実装)で長年挫折していました。

最近になってgam0022さんがLT会でLTCについての発表をしており、いい機会だったのでがっつり調べてみるか！となり、やったのがこれです。当初はそこまで深い話をするつもりではなかったのですが、一旦書いてみると結構あれこれ気になる点がありつらつらと調べていたら結構でかい記事になってしまいました。

LTCは一見するとSpecularだけにしか使えないように誤解されがちですが、実は汎用性が高く、線形変換で近似可能であるならばどのようなBRDFでもArea Lightの計算を可能にしてくれます。実際、Sheenの計算にLTCを使用している論文もありました。思っていたよりLTCは強力な手法だなと感じます。

github.com

この記事をみてオレオレArea Lightシェーダーを書いてくれると幸いです。

~~レイトレアドベントカレンダーの記事なのにレイを飛ばしていない~~

Appendix

線形変換のヤコビアン

LTCの証明で線形変換による立体角のヤコビアンを用いていました。

これの導出をしていきましょう。立体角のヤコビアンはちっちゃい立体角の平面 $d\omega$ と変換によって変化した立体角 $d\omega_o$ の関係を求めて、その比を求めるのが一般的なやり方です。なので線形変換によって微小立体角 $d\omega$ がどう変化するかを見ていきます。

まず、微小立体角 $d\omega$ に対して立体角を構成するベクトル $\omega_1, \omega_2$ を考えます(接ベクトル)。この時、 $d\omega = |\omega_1 \times \omega_2|$ という関係が成り立つとしています。ここで何らかの線形変換 $M$ によって $d\omega$ が変換されるとして、変換後の面積を $A$ 、面の法線と球の方向の内積を $\cos{\theta}$ 、面積までの距離を $r$ とすると変換後の微小立体角 $d\omega_o$ は立体角の定義から

$\displaystyle{ d\omega_o = \frac{A \cos{\theta}}{r^2} }$

となります。すなわち、 $A$ , $\cos{\theta}$ と $r$ さえ求めてしまえばいいわけです。

ではまず $A$ についてですが、 $\omega_1, \omega_2$ とは変換後も接ベクトルであるという関係は変わらないので

$\displaystyle{ A = |M \omega_1 \times M \omega_2| }$

となります。この時、外積と行列の関係性の公式から

$\displaystyle{ M \omega_1 \times M \omega_2 = |M|M^{-T} (\omega_1 \times \omega_2) = |M| M^{-T} \omega d\omega }$

と求めることができます。次に $\cos$ についてですが、 $A$ の法線は

$\displaystyle{ \frac{M \omega_1 \times M \omega_2}{M \omega_1 \times M \omega_2} }$

と書くことができるので、

$\displaystyle{ \cos{\theta} =(\frac{M\omega}{|M\omega|}) \cdot (\frac{M \omega_1 \times M \omega_2}{|M \omega_1 \times M \omega_2|}) }$

となります。ここで内積と行列の公式 $(M v) \cdot (w) = (v) \cdot (M^{-T} w)$ を用いると

$\displaystyle{ \begin{align} (\frac{M\omega}{|M\omega|}) \cdot (\frac{M \omega_1 \times M \omega_2}{|M \omega_1 \times M \omega_2|}) &= (\frac{M\omega}{|M\omega|}) \cdot (\frac{M \omega_1 \times M \omega_2}{|M \omega_1 \times M \omega_2|}) \\ &= (\frac{M\omega}{|M\omega|}) \cdot (\frac{|M| M^{-T} \omega d\omega}{||M| M^{-T} \omega d\omega|}) \\ &= (\frac{M\omega}{|M\omega|}) \cdot (\frac{M^{-T} \omega}{|M^{-T} \omega|}) \\ &= \frac{1}{|M\omega|| M^{-T} \omega|} M \omega \cdot M^{-T} \omega \\ &= \frac{1}{|M\omega|| M^{-T} \omega|} \omega \cdot \omega \\ &= \frac{1}{|M\omega|| M^{-T} \omega|} \\ \end{align} }$

と書くことができます。 $r$ に関しては簡単な話 $r = |M \omega|$ なので、微小立体角 $d\omega_o$ は

$\displaystyle{ \begin{align} d\omega_o &= \frac{A \cos{\theta}}{r^2} \\ &= \frac{|M| |M^{-T} \omega|}{|M\omega|^2} \frac{1}{|M\omega| |M^{-T} \omega|} d\omega \\ &= \frac{|M|}{|M\omega|^3} d\omega \\ \end{align} }$

以上の関係より、ヤコビアンは

$\displaystyle{ \frac{d\omega_o}{d \omega} = \frac{|M|}{|M\omega|^3} }$

となります。 $M$ に関しては任意なので $M \rightarrow M^{-1}$ と置き換えれば理論で使用したヤコビアンが出てくるわけです。

[Heitz, et al. 2016]のLUT

[Heitz, et al. 2016]で公開されているLUTは今回使用したselfshadow氏のものとは結構形が違います。[Heitz, et al. 2016]のLUTでは行列について次のように設定し、Fittingを行っています。

$\displaystyle{ M^{-1} = \begin{pmatrix} 1 &0&b\\ 0&c &0\\ d&0&a \end{pmatrix} }$

マッピングに関してもちょっと異なり、 $1 - 2 \theta_v / \pi$ , $roughness$ という形になっています。

これを用いても大方はselfshadow氏の結果と変わらないのですが、roughnessが小さくて、なおかつ $\theta_v$ が高いときに極端に高い数値を持っていることが原因で、一部描画がおかしくなることが知られています。

selfshadow氏のLUTはFittingのパラメーターを変更することでこの問題を回避しています。なので基本的にselfshadow氏のLUTを使用するのが良いです。

参考、引用文献

[Heitz, et al. 2016] Real-Time Polygonal-Light Shading with Linearly Transformed Cosines
[Heitz 2017] Geometric Derivation of the Irradiance of Polygonal Lights
selfshadow
kanikamaGI
Unity ShaderPlaygroundURP