GLSL
CG
光学
BRDF
More than 1 year has passed since last update.

はじめに

この記事は「基礎からはじめる物理ベースレンダリング」の続きになります.ここでもう一度同じ説明をするものもありますが,復習のつもりで読み流してください.

まず,BRDFとはある領域に入射してくる光(放射照度)が,ある方向に出射する光(放射輝度)を割合として表す関数のことです.双方向反射分布関数と言います.双方向というのは入射方向と出射方向が逆になっても,同じ値になることを表しています.この性質をヘルムホルツの相反性ともいいます.

光が表面にぶつかると鏡面反射したり,物体の内部を透過して吸収されたり,散乱して表面から出射されたりします.表面のある地点に入射した光は透過して散乱し,入射した表面からまた出射される場合があります.入射した位置と出射した位置は必ずしも同じではありませんが,これを同じ位置として考えたものが BRDF でした.

光が表面にぶつかると,光は反射光と内部を透過する光の2つに分かれます.ただし,表面の材質によっては反射光や透過光だけの場合があります.内部を透過する光は入射する前と後の媒質によって屈折することがあります.入射光に対して,ある方向に反射する光の割合はBRDFで表し,ある方向に透過する光の割合はBTDFで表します.この BTDF は双方向透過分布関数といい,BRDF と BTDF を合わせて BSDF(双方向散乱分布関数)といいます.

コンピュータグラフィックスはその用途や環境によって,光の物理現象をどの程度再現するか変化します.リアルタイムレンダリングする場合は BRDF だけで計算したり,オフラインでレイトレースレンダリングするときは BSDF で計算するといったことがあるでしょう.

少し前置きが長くなってしまいましたが,「基礎からはじめる物理ベースレンダリング」に書かれているように,ここでは透過光は吸収されるか,入射した表面の位置から再び出射するものと仮定します.この出射光を拡散反射光といい,入射光のうち,透過する光が拡散反射する割合を表す関数を拡散反射BRDFといいます.

BRDF の性質

物理ベースレンダリングにおいて,BRDF は「ヘルムホルツの相反性」と「エネルギー保存則」に従います.「ヘルムホルツの相反性」は前述したとおり,入射光と出射光を入れ替えても反射の割合は変化しません.「エネルギー保存則」は出射光の量が入射光の量を超えないことです.

実際,物理ベースレンダリングにおいて使われるBRDFがこの性質を満たしているとは限りません.同じBRDFでも,渡すパラメータによってこの性質を満たさない場合があります.最終的なレンダリング画像の見栄えを重視し,かつそれを使用する側が操作しやすくするために,ある程度許容することが必要になる場合もあります.重要なのはそのBRDFの特徴を理解することです.

拡散反射

拡散反射を計算する場合は,ランバート反射がよく使われています.ランバート反射とは,ランバートの余弦則にもとづいて,どの角度から見ても輝度が等しくなることです.つまり,ある方向から入射した光が表面上で,半球面のすべての方向に対して同じ光の量を反射すると考えることができます.これはとても直感的でわかりやすく,計算も簡単です.

近年では,物理ベースレンダリングが使われ,鏡面反射には特に表面の粗さをシミュレートするために微小面モデルがよく使われています.拡散反射についても,表面の粗さによる影響を与えるために微小面モデルが使われることもあり,代表的なものに Oren-Nayar モデルがあります.また,フレネルの影響を拡散反射にも与えることも行われています.この代表的なものに Disney モデルがあります(本来,フレネルは反射光と透過光の割合を表すものです.例えば入射光を 100 としたとき,反射光が 80 で透過光が 20 だとしたら,フレネル反射率は 0.8 となります.ただし,物理ベースレンダリングでは,このフレネル反射率をメタルネスというパラメータで操作することが多く,ランバート反射ではフレネルの影響を受けません.ちなみに,鏡面反射では Torrance-Sparrowモデルをベースにした微小面モデルがよく使われており,これはフレネルの影響が含まれています).

本記事では特に Disney モデルと Oren-Nayar モデルについて説明します.また,最後に少しだけ最近の Diffuse BRDF モデルを紹介します.

Disney モデル

Burley氏の論文 "Physically-Based Shading at Disney" には,フレネルの影響を加えた拡散反射BRDFが記載されています.

まずは,以下の式見てください.

(1-F(\theta_{l}))(1-F(\theta_{d}))

$F(\theta)$ はフレネル反射率です.フレネル反射率は入射光に対する反射光の割合で 0.0 から 1.0 の間の値です.これを 1.0 から引いた値は透過光の割合になります.透過光は拡散反射光と仮定していましたので,拡散反射光の割合と見なすことができます.$\theta_{l}$ は入射角,$\theta_{d}$ はハーフベクトル $\vec{h}$ と入射ベクトル $\vec{l}$ との間の角度です.ハーフベクトルとは2つのベクトルの中間を指すベクトルのことで,ここでは入射ベクトルと反射ベクトル $\vec{v}$ とのハーフベクトルのことです.入射方向に対して反対のベクトルのことをライトベクトル,反射ベクトルのことを視線ベクトルともいいます.

h = \frac{\vec{v}+\vec{l}}{|\vec{v}+\vec{l}|}

粗さを考慮するために反射ベクトル方向ではなくて,ハーフベクトル方向に対して反射すると見なしています.入射と反射を乗算することで,入射方向と反射方向を入れ替えても変わらないようにしています.これはヘルムホルツの相反性を満たすためです.

フレネル反射率の計算には媒質の屈折率が必要ですが,ここでは無視します.代わりに表面の粗さによって変化するようにします.特に表面が粗く,グレージング角に近づくほど反射率を高くするように調整します.これらを考慮し,Schlick の近似式を使ってフレネル反射率を求める式は次のようになります.

F = (1+(F_{D90}-1)(1-\cos\theta_{l})^5)
(1+(F_{D90}-1)(1-\cos\theta_{v})^5)

ここで

F_{D90} = 0.5 + 2\ \cos^2(\theta_{d})\ roughness

$roughness$ は表面の粗さを表すパラメータです.

この式をランバート反射に加えた拡散反射BRDF,$f_{d}$ は次のようになります.

f_{d} = \frac{\sigma}{\pi}(1+(F_{D90}-1)(1-\cos\theta_{l})^5)
(1+(F_{D90}-1)(1-\cos\theta_{v})^5)

$\sigma$ はリフレクタンスです.アルベドともいいます.下記は GLSL コードです.

float SchlickFresnel(float u, float f0, float f90)
{
  return f0 + (f90-f0)*pow(1.0-u,5.0);
}

float DisneyDiffuse(float albedo, vec3 N, vec3 L, vec3 V, float roughness)
{
  vec3 H = normalize(L+V);
  float dotLH = saturate(dot(L,H));
  float dotNL = saturate(dot(N,L));
  float dotNV = saturate(dot(N,V));
  float Fd90 = 0.5 + 2.0 * dotLH * dotLH * roughness;
  float FL = SchlickFresnel(1.0, Fd90, dotNL);
  float FV = SchlickFresnel(1.0, Fd90, dotNV);
  return (albedo*FL*FV)/PI;
}

正規化

この Disney拡散反射BRDF はエネルギー保存則を満たさない場合があります.

diffuse_brdf_figXX_plot_disney.png

この図はある領域に入射した光を前述のBRDFを使って表面上の半球すべての方向に対して反射した光の総和をグラフにしたものです.ラフネスと入射角を軸にしています.表面が粗く,グレージング角に近づくと 1.0 を超えているところが見つかります.

そこで,EA のゲームエンジン Frostbite ではこれを正規化しました.実際のコードは次のようになります.

float SchlickFresnel(float u, float f0, float f90)
{
  return f0 + (f90-f0)*pow(1.0-u,5.0);
}

float NormalizedDisneyDiffuse(float albedo, vec3 N, vec3 L, vec3 V, float roughness)
{
  float energyBias = mix(0.0, 0.5, roughness);
  float energyFactor = mix(1.0, 1.0/1.51, roughness);
  float Fd90 = energyBias + 2.0 * dotLH * dotLH * roughness;
  float FL = SchlickFresnel(1.0, Fd90, dotNL);
  float FV = SchlickFresnel(1.0, Fd90, dotNV);
  return (albedo*FL*FV*energyFactor)/PI;
}

このBRDFのグラフは次のようになります.

diffuse_brdf_figXX_plot_renormalized_disney.png

最大値が 1.0 に収まるようになっています.

Oren-Nayar モデル

ランバート反射面では表面の粗さを考慮していませんでした.Disneyモデルではフレネルの効果を反映させるために粗さを使用しました.Oren-Nayar モデルは表面が微小面で構成されていると考えた反射モデルです.この微小面モデルは Torrance-Sparrow モデルがベースとなっており,鏡面反射でよく使われている Cook-Torrance モデルと似ています.

V字型の凹み

物体表面のある微小領域をパッチ(1区画)といいます.パッチは複数の微小面が集まったものになります.このとき,微小面は図のようにV字型の凹みが並んでいるものと仮定します.

diffuse_brdf_figXX_rough_surface.png

2枚の微小面を組み合わせて1つの凹みを表します.

diffuse_brdf_figXX_vcavity.png

このとき,凹みは表面の上から押し出して出来たものと思ってください.つまり,V字型の凹みは断面を見ると2等辺三角形で,そのうち2頂点が表面上にあることになります.また,もう1頂点は表面上の2点の中点を通る法線軸上にあります.

diffuse_brdf_figXX_vcavity_vertices.png

放射定義

diffuse_brdf_figXX_brdf_coordinate.png

光源からの方向 $\vec{s} = (\theta_{i},\phi_{i})$ から表面積 $dA$ に入射し,$\vec{v} = (\theta_{r},\phi_{r})$ に出射されるとき,放射照度は次のようになります.

E(\theta_{i},\phi_{i}) = \frac{d\Phi_{i}(\theta_{i},\phi_{i})}{dA}

$\vec{v}$ 方向に出射する放射輝度は次のようになります.

L_{r}(\theta_{r},\phi_{r};\theta_{i},\phi_{i}) = 
\frac{
  d^2\Phi(\theta_{r},\phi_{r};\theta_{i},\phi_{i})
}{
  dA\cos\theta_{r}d\omega_{r}
}

BRDFは次のようになります.

f_{r}(\theta_{r},\phi_{r};\theta_{i},\phi_{i}) = \frac{
  dL_{r}(\theta_{r},\phi_{r};\theta_{i},\phi_{i})
}{
  dE(\theta_{i},\phi_{i})
}

表面が等方性の場合,BRDF は次のように簡単になります.

f_{r}(\theta_{r},\theta_{i},\phi_{r}-\phi_{i}) = \frac{
  dL_{r}(\theta_{r},\theta_{i},\phi_{r}-\phi_{i})
}{
  dE(\theta_{i})
}

粗い表面

diffuse_brdf_figXX_rough_microfacet.png

微小面の傾きと向きは $(\theta_{a},\phi_{a})$ です.法線 $\vec{a} = (\theta_{a},\phi_{a})$ を持つ微小面の数を表す分布は $N(\theta_{a},\phi_{a})$ です.ここで,表面積のうち微小面が占める割合を表したものを勾配面積分布 $P(\theta_{a},\phi_{a})$ といいます.微小面の数を表す分布と勾配面積分布は次のような関係になります.

P(\theta_{a},\phi_{a}) = N(\theta_{a},\phi_{a})da\cos\theta_{a}

表面が等方性の場合,$N(\theta_{a},\phi_{a}) = N(\theta_{a})$ と $P(\theta_{a},\phi_{a}) = P(\theta_{a})$ と表せます.

粗い表面をV字型の凹みが並んでいると仮定しました.これが表面反射にどのような影響を与えるのか考えてみます.

diffuse_brdf_figXX_vcavity_lighting.png

上図はV字型の凹みに光源から入射していることを表しています.微小面の法線と光源への向きが近いほど光を多く受け取ります.2つの微小面のうち左側が明るく,右側が暗くなります.

diffuse_brdf_figXX_lighter_darker.png

左側の図では,視点からだと暗い微小面の方が多く,右側の図では明るい微小面の方が見える割合が多いことがわかります.このことから,光源の方向に視点が近づくほど表面が明るくなります.なお,微小面はランバート反射するものとします.これは,鏡面反射の微小面モデルだと微小面は鏡面反射するものとしますので違いがあります.

あるパッチ(表面積は $dA$)はV字型の凹みが並んでいます.V字型の凹みは2つの微小面で構成されています.パッチの放射発散度(パッチから出射する放射輝度の総和)は,パッチに含まれている微小面の反射した放射束の総和になります.微小面の表面積は $da$ で,投影面積は $da\cos\theta_{a}$ です.パッチから出射する放射輝度は微小面の投影面積を使います.式で表すと次のようになります.

L_{rp}(\theta_{a}, \phi_{a}) = \frac{
  d^2\Phi(\theta_{a}, \phi_{a})
}{
  (da\cos\theta_{a})\cos\theta_{r}d\omega_{r}
}

あるパッチの放射発散度は,パッチに含まれる微小面からの反射 $L_{rp}(\theta_{a},\phi_{a})$ の平均となります.

L_{r}(\theta_{r},\phi_{r};\theta_{i},\phi_{i}) = 
\int_{\theta_{a} = 0}^{\pi/2}
  \int_{\phi_{a} = 0}^{2\pi}
    P(\theta_{a},\phi_{a})L_{rp}(\theta_{a},\phi_{a})\sin\theta_{a}d\phi_{a}d\theta_{a}

単方向勾配分布モデル

最初はすべての微小面が同じ方向 $\theta_{a}$ を向いているモデルについて考えていきます.V字型の凹みを構成する2枚の微小面は $\phi_{a}$ と $\phi_{a}+\pi$ を向いて隣接しています.これは視点によって見た目が変わるため,異方性モデルになりますので,異方性モデルから等方性モデルに対応し,ガウス分布モデルに発展していきます.

シャドウイング,マスキングともに行われないランバート面とすると,微小面の放射輝度は放射照度に比例し,$\frac{\rho}{\pi}E(\theta_{a},\phi_{a})$ に等しいです.微小面の放射照度は $E(\theta_{a},\phi_{a}) = E_{0}(\vec{s} \cdot \vec{a})$ です.ここで $E_{0}$ は微小面が光源に向いたとき($\vec{s} = \vec{a}$) の放射照度です.微小面から視点方向に反射される放射束は

d^2\Phi = \frac{\rho}{\pi}E_{0}(\vec{s}\cdot\vec{a})(\vec{v}\cdot\vec{a})dad\omega_{r}

これを式に代入すると

L_{rp}(\theta_{a},\phi_{a}) = \frac{\rho}{\pi}
  E_{0}\frac{
    (\vec{s}\cdot\vec{a})(\vec{v}\cdot\vec{a})
  }{
    (\vec{a}\cdot\vec{n})(\vec{v}\cdot\vec{n})
  }

式に $\vec{v}$ が含まれているので,視点によって変化します.

幾何減衰項

微小面によって発生するシャドウイング(自己陰影),マスキング(自己遮蔽)を表したものが幾何減衰項です.

diffuse_brdf_figXX_shadowing_masking.png

左図がシャドウイングで入射光が遮られています.右図がマスキングで反射光が遮られています.パターンとしてはシャドウイングとマスキングなし,シャドウイングのみ,マスキングのみ,シャドウイングとマスキングありの4種類が考えられます.

diffuse_brdf_figXX_local_occulusion.png

微小面の長さと幅をそれぞれ $l$ と $w$ とします.シャドウィング,マスキングされたときの幅はそれぞれ $m_{s}, m_{v}$ です.視点から照明された微小面の領域は $l\cdot\min(w-m_{s}, w-m_{v})$ です.幾何減衰項はこの式から微小面の面積 $wl$ で除算することで求められます.

G = \min(1-\frac{m_{s}}{w}, 1-\frac{m_{v}}{w})

diffuse_brdf_figXX_shadowing_masking2.png

入射時のシャドウイング幅 $m_{s}$ と入射角($\theta_{i}$)(または反射角 $(\theta_{r})$)を使って幾何減衰項を表現します.

diffuse_brdf_figXX_shadowing_01.png

x_{w} = w\cos\theta_{a} \qquad
z_{w} = w\sin\theta_{a} \qquad
x_{s} = m_{s}\cos\theta_{a} \qquad
z_{s} = m_{s}\sin\theta_{a}

diffuse_brdf_figXX_shadowing_02.png

x_{n} = x_{w} + x_{s} \qquad
z_{n} = z_{w} - z_{s}
\begin{align}
n\sin\theta_{i} &= x_{n} = m_{s}\cos\theta_{a} + w\cos\theta_{a} \\[2ex]
n\cos\theta_{i} &= z_{n} = -m_{s}\sin\theta_{a} + w\sin\theta_{a}
\end{align}

最初の式に $\cos\theta_{i}$ を,2番目の式に $-\sin\theta_{i}$ を掛けると

\begin{align}
n\sin\theta_{i}\cos\theta_{i} &= m_{s}\cos\theta_{a}\cos\theta_{i} + w\cos\theta_{a}\cos\theta_{i} \\[2ex]
-n\cos\theta_{i}\sin\theta_{i} &= m_{s}\sin\theta_{a}\sin\theta_{i} - w\sin\theta_{a}\sin\theta_{i}
\end{align}

最初の式と2番目の式を足すと

0 = m_{s}\cos\theta_{a}\cos\theta_{i} + w\cos\theta_{a}\cos\theta_{i} + m_{s}\sin\theta_{a}\sin\theta_{i} - w\sin\theta_{a}\sin\theta_{i}

$m_{s}$ と $w$ で括ると

0 = m_{s}(\cos\theta_{a}\cos\theta_{i} + \sin\theta_{a}\sin\theta_{i}) + w(\cos\theta_{a}\cos\theta_{i} - \sin\theta_{a}\sin\theta_{i})

括弧の中は加法定理で次のように変形できます.

0 = m_{s}\cos(\theta_{a}+\theta_{i}) + w\cos(\theta_{a}-\theta_{i})

さらに $m_{s}/w$ の形に整理すると

\frac{m_{s}}{w} = -\frac{\cos(\theta_{a}+\theta_{i})}{\cos(\theta_{a}-\theta_{i})}

最終的に

\begin{align}
1-\frac{m_{s}}{w} &= 1-\left(-\frac{
\cos(\theta_{a}+\theta_{i})
}{
\cos(\theta_{a}-\theta_{i})
}\right) \\[2ex]
&= \frac{
\cos(\theta_{a}-\theta_{i})
}{
\cos(\theta_{a}-\theta_{i})
}+\frac{
\cos(\theta_{a}+\theta_{i})
}{
\cos(\theta_{a}-\theta_{i})
} \\[2ex]
&= \frac{
\cos(\theta_{a}-\theta_{i})+\cos(\theta_{a}+\theta_{i})
}{
\cos(\theta_{a}-\theta_{i})
} \\[2ex]
&= \frac{
2\cos\theta_{a}\cos\theta_{i}
}{
\cos(\theta_{a}-\theta_{i})
}
\end{align}

マスキングについては $\theta_{i}$ を $\theta_{r}$ に変えることで計算できます.シャドウイングとマスキングを計算し,値が大きい(照明されている領域が多い)側を選択します.また,幾何減衰項は0~1の値になりますので,その範囲内に収まるようにします.

G = \min\left(
  1, \max\left(
    0, \frac{
      2\cos\theta_{i}\cos\theta_{a}
    }{
      \cos(\theta_{i}-\theta_{a})
    }, \frac{
      2\cos\theta_{r}\cos\theta_{a}
    }{
      \cos(\theta_{r}-\theta_{a})
    }
  \right)
\right)

なお,この式では微小面の法線 $\vec{a}$ は任意の方向であり,光源と視点方向とのハーフベクトルである必要はありません.

一般的な幾何減衰項

パッチの法線に対して直交していて,パッチの法線方向を含む平面上に光源と視点があるという状況に限定していましたが,通常は光源方向と視線方向は任意の方向を向いています.そのため,光源方向と視線方向を垂直な平面に投影して計算します.

diffuse_brdf_figXX_gaf_proj_01.png

この図では,パッチの法線は $z$ 方向で,パッチの平面は xy 平面になっています.視線方向を $\vec{v}$ とし,このベクトルは正規化されているとします.これを xz 平面に投影することを考えます.

まず,視線方向と z 方向の平面において

diffuse_brdf_figXX_gaf_proj_02.png

$\vec{v}$ の長さは正規化されているので 1 です.底辺は $\cos\theta$,対辺は $\sin\theta$ になります.

diffuse_brdf_figXX_gaf_proj_03.png

次に $\phi$ について.$\phi$ を含む三角形の底辺は $\theta$ を含む三角形の対辺と同じなので $\sin\theta$ になります.対辺は $\sin\theta\cos\phi$ になります.

diffuse_brdf_figXX_gaf_proj_04.png

最後に xz 平面に投影します.$\theta_{p}$ を含む三角形の底辺 $a$ は $\cos\theta$,対辺 $b$ は $\sin\theta\cos\phi$ ですので,斜辺 $c$ はピタゴラスの定理から次のように求められます.

c = \sqrt{cos^2\theta + \sin^2\theta\cos^2\phi}

$\theta_{p}$ についてまとめると

\begin{align}
\sin\theta_{p} &= \frac{b}{c} = \frac{\sin\theta\cos\phi}{
  \sqrt{cos^2\theta + \sin^2\theta\cos^2\phi}
} \\[2ex]
\cos\theta_{p} &= \frac{a}{c} = \frac{\cos\theta}{
  \sqrt{cos^2\theta + \sin^2\theta\cos^2\phi}
} \\[2ex]
\tan\theta_{p} &= \frac{b}{a} = \frac{\sin\theta}{\cos\theta}\cos\phi = \tan\theta\cos\phi
\end{align}

投影された角度 $\theta_{p}$ を垂直V字型の式にある $\theta_{i}, \theta_{r}$ のところに代入します.ここで $c = \sqrt{cos^2\theta + \sin^2\theta\cos^2\phi}$ を使うと簡単です.

\begin{align}
\frac{
  2\cos\theta_{i}\cos\theta_{a}
}{
  \cos(\theta_{i}-\theta_{a})
} &= \frac{
  2\cos\theta_{i}\cos\theta_{a}
}{
  \cos\theta_{i}\cos\theta_{a}+\sin\theta_{i}\sin\theta_{a}
} \\[2ex]
&= \frac{
  \frac{2\cos\theta_{i}\cos\theta_{a}}{c}
}{
  \frac{
    \cos\theta_{i}\cos\theta_{a}+\sin\theta_{i}\sin\theta_{a}\cos\phi
  }{c}
} \\[2ex]
&= \frac{2\cos\theta_{i}\cos\theta_{a}}{c}
\frac{c}{
 \cos\theta_{i}\cos\theta_{a}+\sin\theta_{i}\sin\theta_{a}\cos\phi
}
 \\[2ex]
&= \frac{
  2\cos\theta_{i}\cos\theta_{a}
}{
  \cos\theta_{i}\cos\theta_{a}+\sin\theta_{i}\sin\theta_{a}\cos\phi
}
\end{align}

一般的な幾何減衰項の式は次のようになります.

G = \min\left(
  1, \max\left(
    0, \frac{
      2\cos\theta_{i}\cos\theta_{a}
    }{
      \cos\theta_{i}\cos\theta_{a}+\sin\theta_{i}\sin\theta_{a}\cos\phi
    }, \frac{
      2\cos\theta_{r}\cos\theta_{a}
    }{
      \cos\theta_{r}\cos\theta_{a}+\sin\theta_{r}\sin\theta_{a}\cos\phi
    }
  \right)
\right)

この式をベクトルで表すと次のようになります.

G = \min\left(
  1, \max\left(
    0, \frac{
      2(\vec{s}\cdot\vec{n})(\vec{a}\cdot\vec{n})
    }{
      \vec{s}\cdot\vec{a}
    }, \frac{
      2(\vec{v}\cdot\vec{n})(\vec{a}\cdot\vec{n})
    }{
      \vec{v}\cdot\vec{a}
    }
  \right)
\right)

相互反射

Oren-Nayar モデルの特徴にこの相互反射があります.相互反射とは微小面同士で反射して放射輝度に寄与するものです.

diffuse_brdf_figXX_interreflection.png

図を見ると本来マスキングで見えなくなる部分が反対側の微小面でさらに反射して視線方向に出射しています.Oren-Nayarモデルではこの現象を計算するために,最大2回の反射をモデル化しました.

Oren-Nayar モデルでは,この相互反射と直接照明を別々に考えることにし,区別するために前者の放射輝度を $L_{rp}^2$,後者の放射輝度を $L_{rp}^1$ とします.

相互反射の放射輝度 $L_{rp}^2$ は次のようになります.

L_{rp}^2 = \frac{\sigma}{\pi}\int\int K(\vec{x},\vec{y}) L_{rp}^1(\vec{y})d\vec{y}

ここで,$\vec{x}$ は隣接する微小面上のすべての点 $\vec{y}$ の放射輝度を受け取って反射する反対側の微小面上の位置です.$K(\vec{x},\vec{y})$ はカーネルと呼ばれていて,$\vec{x}$ と $\vec{y}$ の関係を表す関数です.V字型の凹みのような形状の場合は1次元の形状と考えることができるので,1次方程式で書き表せます.

L_{rp}^2 = \frac{\sigma}{\pi}\int K'(x,y) L_{rp}^1(y)dy

ここで,$x$ と $y$ は2つの微小面の交点からそれぞれの微小面上の最短位置となります.$K'$ は次のようになります.

K'(x,y) = \frac{\pi\sin^2(2\theta_{a})}{2}\frac{
  xy
}{
  (x^2+2xy\cos(2\theta_{a})+y^2)^{3/2}
}

相互反射を考慮した放射輝度は,光源から照らされる領域と,視線から見える領域に限定されます.つまり $w-m_{s}$, $w-m_{v}$ の範囲を積分します.

L_{rp}^2 = \frac{
  l(\vec{a}\cdot\vec{v})
}{
  da(\vec{a}\cdot\vec{n})(\vec{v}\cdot\vec{n})
} \int_{x=m_{v}}^{w}L_{rp}^2(x)dx = \left(
\frac{\rho}{\pi}\right)^2 E_{0}\frac{
  l(\vec{a}\cdot\vec{s})(\vec{a}\cdot\vec{v})
}{
  da(\vec{a}\cdot\vec{n})(\vec{v}\cdot\vec{n})
} \int_{x=m_{v}}^w \int_{y=m_{s}}^w K'(x,y)dxdy

ここで $t = \frac{x}{w}; r=\frac{y}{w}$ として置換積分すると

L_{rp}^2 = \left(
\frac{\rho}{\pi}\right)^2 E_{0}\frac{
  (\vec{a}\cdot\vec{s})(\vec{a}\cdot\vec{v})
}{
  (\vec{a}\cdot\vec{n})(\vec{v}\cdot\vec{n})
} \int_{t=\frac{m_{v}}{w}}^1 \int_{r=\frac{m_{s}}{w}}^1 K'(t,r)drdt

先ほど定義したカーネル関数を使うと,上記の積分は次のようになります.

\int_{t=\frac{m_{v}}{w}}^1 \int_{r=\frac{m_{s}}{w}}^1 K'(t,r)drdt = 
\frac{\pi}{2}\left[
  d(1,\frac{m_{v}}{w}) + d(1,\frac{m_{s}}{w}) - d(\frac{m_{s}}{w},\frac{m_{v}}{w}) - d(1,1)
\right] \\
d(x,y) = \sqrt{x^2+2xy\cos(2\theta_{a})+y^2}

この $\frac{\pi}{2}\left[d(1,\frac{m_{v}}{w}) + d(1,\frac{m_{s}}{w}) - d(\frac{m_{s}}{w},\frac{m_{v}}{w}) - d(1,1)\right]$ を相互反射係数と呼びます.

微小面から放射される量は直接反射による放射輝度 $L_{rp}^1$ と相互反射による放射輝度 $L_{rp}^2$ を合わせたものになります.

L_{rp}(\theta_{a},\phi_{a}) = L_{rp}^1(\theta_{a},\phi_{a}) + L_{rp}^2(\theta_{a},\phi_{a})

V字型の凹みは2つの微小面で構成されていますので,ここからの放射輝度は法線 $(\theta_{a},\phi_{a})​$ と $(\theta_{a},\phi_{a}+\pi)​$ を向いている微小面からの放射輝度の平均になります.

L_{r}(\theta_{r},\phi_{r};\theta_{i},\phi_{i};\theta_{a},\phi_{a}) = \frac{
  L_{rp}(\theta_{a},\phi_{a}) + L_{rp}(\theta_{a},\phi_{a}+\pi)
}{2}

等方性モデル

これまですべての微小面が同じ傾き($\theta_{a}​$)として考えてきましたが,ここで $\theta_{a}​$ 方向に一様に分布する微小面を考えます.法線 $\vec{n} = (\theta_{a},\phi_{a})​$ の微小面の直接照明による放射輝度は $L_{rp}^1(\theta_{a},\phi_{a})​$ です.したがって等方性表面において, 直接照明による放射輝度は全方位の放射輝度の積分となります.

L_{rp}^1(\theta_{a}) = \frac{1}{2\pi}\int_{\phi_{a}=0}^{2\pi}L_{rp}^1(\theta_{a},\phi_{a})d\Phi

この積分は複雑なのでまともに計算するのは難しいです.そこで,$\phi_{a}$ を3つの範囲に考えて計算し,加重和をとる方法が提示されました.範囲と計算式については論文を参照してください.

結果として次のような式になります.

L_{rp}^1(\theta_{a}) = \frac{\rho}{\pi}\cos\theta_{i}\cos\theta_{a}\left[
  1+
  \cos(\phi_{r}-\phi_{i})(A_{1}(\alpha;\theta_{a})\tan\beta+A_{2}(\beta,\phi_{r}-\phi_{i};\theta_{a}))+
  (1-|\cos(\phi_{r}-\phi_{i})|)A_{3}(\theta_{r},\theta_{i};\theta_{a})
\right]

$\alpha = \max(\theta_{i},\theta_{r})$, $\beta = \min(\theta_{i},\theta_{r})$ で,$A_{1},A_{2},A_{3}$ は次のようになります.

\begin{align}
A_{1}(\alpha;\beta) &= \tan\theta_{a}\frac{2\sin\phi_{c}^\alpha}{\pi}+
\frac{1}{2}\tan^2\theta_{a}\tan\alpha(1-\frac{
  2\phi_{c}^\alpha + \sin(2\phi_{c}^\alpha)
}{\pi}) \\[2ex]
A_{2}(\beta,\phi_{r}-\phi_{i};\theta_{a}) &= \begin{cases}
  \frac{2\phi_{c}^\beta}{\pi}-\tan\theta_{a}\tan\theta_{\beta}\frac{2\sin\phi_{c}^\beta}{\pi} & \text{if $\cos(\phi_{r}-\phi_{i})<0$} \\[2ex]
0 & \text{if $\cos(\phi_{r}-\phi_{i})\geq 0$}
\end{cases} \\[2ex]
A_{3}(\theta_{r},\theta_{i};\theta_{a}) &= \begin{cases}
0 & \text{if $\phi_{c}^i + \phi_{c}^r \leq \frac{\pi}{2}$} \\[2ex]
\frac{1}{2} -
\frac{\phi_{c}^i + \phi_{c}^r}{\pi} +
\frac{
  \sqrt{\tan^2\theta_{r}\tan^2\theta_{a}-1}
}{\pi} +
\frac{
  \sqrt{\tan^2\theta_{i}\tan^2\theta_{a}-1}
}{\pi} -
\frac{
  \tan\theta_{a}\sqrt{
    \tan^2\theta_{i}+\tan^2\theta_{r}
  }
}{\pi} & \text{if $\phi_{c}^i + \phi_{c}^r > \frac{\pi}{2}$}
\end{cases}
\end{align}

次に相互反射の方です.相互反射による放射輝度も直接照明と同じです.

L_{rp}^2(\theta_{a}) = \frac{1}{2\pi}\int_{\phi_{a}=0}^{2\pi}L_{rp}^2(\theta_{a},\phi_{a})d\Phi

相互反射係数については次の近似式が使われます.

\frac{\pi}{2}\left[
  d(1,\frac{m_{v}}{w}) + d(1,\frac{m_{s}}{w}) - d(\frac{m_{s}}{w},\frac{m_{v}}{w}) - d(1,1)
\right] \approx
\pi(1-\cos\theta_{a})(1-m_{s})(1-m_{v})

相互反射の計算式は次のようになります.

L_{rp}^2 \approx \frac{\rho^2}{\pi}E_{0}\cos\theta_{i}\cos\theta_{a}(1-\cos\theta_{a})\left[
  1-
  \cos(\phi_{r}-\phi_{i})\left(
    \frac{2\phi_{c}^\beta}{\pi}+
    2\tan\theta_{a}\tan\beta\frac{
      \sin\phi_{c}^\alpha-\sin\phi_{c}^\beta
    }{\pi}+
    \frac{1}{2}\tan^2\theta_{a}\tan\alpha\tan\beta(1-\frac{
      2\phi_{c}^\alpha+\sin(2\phi_{c}^\alpha)}{\pi}
    )
  \right)
\right]

ガウス分布モデル

等方性モデルでは微小面の傾きは一様分布でしたが,現実の表面に近づけるために,Oren-Nayar モデルではガウス分布を使用します.勾配分布関数 $P(\theta_{a},\phi_{a})$ を使ってモデル化します.等方性モデルでは微小面が $\phi_{a}$ に一様分布しているので,$\theta_{a}$ で記述することができます.

L_{r}(\theta_{r},\theta_{i},\phi_{r}-\phi_{i}) = \int_{0}^{\frac{\pi}{2}}P(\theta_{a})L_{rp}(\theta_{a})\sin\theta_{a}d\theta_{a}

勾配分布関数はガウス分布なので次のようになります.

P(\theta_{a}) = ce^{-\frac{\theta_{a}^2}{2\sigma^2}}

正規化定数 $c$ は次のようになっています.

\frac{1}{c} = \int_{\theta_{a}=0}^{\frac{pi}{2}}
  \int_{\phi_{a}=0}^{2\pi}
    e^{-\frac{\theta_{a}^2}{2\sigma^2}}
    \sin\theta_{a}d\phi_{a}d\theta_{a}

これもまた計算が複雑なので,次のような近似式が提示されました.

L_{r}^1(\theta_{r},\theta_{i},\phi_{r}-\phi_{i};\sigma) =
\frac{\rho}{\pi}E_{0}\cos\theta_{i}
\left[
  C_1{\rho}+
  \cos(\phi_{r}-\phi_{i})C_{2}(\alpha;\beta;\phi_{r}-\phi_{i};\sigma)\tan\beta+
  (1-|\cos(\phi_{r}-\phi_{i})|)C_{3}(\alpha;\beta;\sigma)
\right]

係数は次のようになっています.

\begin{align}
\alpha &= \max(\theta_{i},\theta_{r}) \\[2ex]
\beta &= \min(\theta_{i},\theta_{r}) \\[2ex]
C_{1} &= 1-0.5\frac{\sigma^2}{\sigma^2+0.33} \\[2ex]
C_{2} &= \begin{cases}
  0.45\frac{\sigma^2}{\sigma^2+0.09}\sin\alpha & \text{if $\cos(\phi_{r}-\phi_{i}) \geq 0$} \\[2ex]
  0.45\frac{\sigma^2}{\sigma^2+0.09}\left(
    \sin\alpha-(\frac{2\beta}{\pi})^3
  \right) & \text{otherwise} \\[2ex]
\end{cases} \\[2ex]
C_{3} &= 0.125\left(
  \frac{\sigma^2}{\sigma^2+0.09}
\right)\left(
  \frac{4\alpha\beta}{\pi^2}
\right)^2
\end{align}

同じように相互反射の放射輝度も近似式が提示されています.

L_{r}^2(\theta_{r},\theta_{i},\phi_{r}-\phi_{i};\sigma) = 
0.17\frac{\rho^2}{\pi}E_{0}\cos\theta_{i}\frac{\sigma^2}{\sigma^2+0.13}\left[
  1-\cos(\phi_{r}-\phi_{i})\left(
    \frac{2\beta}{\pi}
  \right)
\right]

この2つの成分を合わせたものが最終的な放射輝度になります.

L_{r}(\theta_{r},\theta_{i},\phi_{r}-\phi_{i};\sigma) =
L_{r}^1(\theta_{r},\theta_{i},\phi_{r}-\phi_{i};\sigma) +
L_{r}^2(\theta_{r},\theta_{i},\phi_{r}-\phi_{i};\sigma)

定性モデル

導出した反射モデルは処理が複雑です.そこで,もっと単純な式に近似したモデルも提示されました.

L_{r}(\theta_{r},\theta_{i},\phi_{r}-\phi_{i};\sigma) = 
\frac{\rho}{\pi}E_{0}\cos\theta_{i}(
  A + B \max[0, \cos(\theta_{r}-\theta_{i})] \sin\alpha\tan\beta
)
\begin{align}
A &= 1.0-0.5\frac{\sigma^2}{\sigma^2+0.33} \\[2ex]
B &= 0.45\frac{\sigma^2}{\sigma^2+0.09}
\end{align}

このモデルは精度を犠牲にして計算を少なくしています.また,相互反射成分も無くしています.Wikipedia に記載されている Oren-Nayar モデルはこのモデルになっていますね.ちなみに定性モデルで相互反射成分を補償したい場合は係数 $A$ を次のようにできると論文に書いてありました.

A = 1.0-0.5\frac{\sigma^2}{\sigma^2+0.57}

スーパーラフ

Oren-Nayar モデルは粗さの値 ($\sigma$) が 1.0 以上を扱えます.1.0 以上の粗さ(ラフネス)のことをスーパーラフというそうです.スーパーラフについては以下を参照してください.

マイクロファセットの分布関数 D(m) について

様々な Oren-Nayar モデル

Oren-Nayar モデルにはいくつかの問題を抱えています.それらに対して改良された様々なモデルが発表されています.最後にいくつかのモデルを紹介します.

Tiny Improved Oren-Nayar Model

Fujii Yasuhiro氏が発表しているモデルです.

A tiny improvement of Oren-Nayar reflectance model

定性モデルでは精度を犠牲にしているため,完全な Oren-Nayar モデルとはかなり異なっています.しかし,完全な Oren-Nayar モデルは計算コストが高いです.そこで,完全な Oren-Nayar モデルと同じような挙動をする近似モデルを導出しました.

\begin{align}
s &= (\vec{s}\cdot\vec{v})-(\vec{n}\cdot\vec{s})(\vec{n}\cdot\vec{v}) \\[2ex]
t &= \begin{cases}
1 & \text{if $s \leq 0$} \\[2ex]
\max(\vec{n}\cdot\vec{s}, \vec{n}\cdot\vec{v}) & \text{if $s > 0$}
\end{cases} \\[2ex]
L_{r}(\vec{n},\vec{s},\vec{v}) &= \rho(\vec{n}\cdot\vec{s})(A+B\frac{s}{t})
\end{align}

係数は次のとおりです.

\begin{align}
A &= \frac{1}{\pi}\left(1-0.5\frac{\sigma^2}{\sigma^2+0.33}+0.17\rho\frac{\sigma^2}{\sigma^2+0.13}\right) \\[2ex]
B &= \frac{1}{\pi}\left(0.45\frac{\sigma^2}{\sigma^2+0.09}\right)
\end{align}

さらに,粗さ($\sigma​$) を $0 \leq \sigma' \leq 1​$ に限定することで,計算を減らす係数も提示されています.

\begin{align}
A &= \frac{1}{\pi+\left(\frac{\pi}{2}-\frac{2}{3}\right)\sigma'} \\[2ex]
B &= \frac{\sigma'}{\pi+\left(\frac{\pi}{2}-\frac{2}{3}\right)\sigma'}
\end{align}

Normalized Oren-Nayar with Fresnel Model

tri-ace の五反田氏は Oren-Nayar モデルがエネルギー保存則を満たしていないため,フレネルも考慮して改良したモデルを導出しました [GOT14].

\begin{align}
F_{diff}(\vec{s},\vec{v}) &= \frac{21}{20}(1-f_{0})(1-(1-\vec{n}\cdot\vec{s})^5)(1-(1-\vec{n}\cdot\vec{v})^5) \\[2ex]
L_{r} &= \frac{\rho}{\pi}E_{0}(1-f_{0})\left[
  F_{diff}(\vec{s},\vec{v})(\vec{n}\cdot\vec{s})(1.0-0.5\frac{\sigma^2}{\sigma^2+0.65})+
  \left(
    (0.45\frac{\sigma^2}{\sigma^2+0.09})
    (\vec{v}\cdot\vec{s} - (\vec{n}\cdot\vec{v})(\vec{n}\cdot\vec{s}))
    \min(1, \frac{\vec{n}\cdot\vec{s}}{\vec{n}\cdot\vec{v}})
  \right)
\right]
\end{align}

Oren-Nayar with GGX+Smith Model

tri-ace の五反田氏は Oren-Nayar モデルで使われている法線分布関数と幾何減衰項の計算方法がガウス分布およびV字型の凹みになっているため,鏡面反射モデルでよく使われている GGX と Smith を適用したモデルを発表しました[GOT14].

L_{r}(e,l,\sigma,f_{0}) = \frac{21}{20\pi}(1-f_{0})(
  F_{r}(e,l,\sigma)L_{m}(l,\sigma) + V_{d}(e,l,\sigma)B_{p}(e,l)
) \tag{8}
F_{r}(e,l,\sigma) = \left(
  1-\frac{0.542026\sigma^2+0.303573\sigma}{\sigma^2+1.36053}
\right) \left(
  1-\frac{
    (1-n\cdot e)^{5-4\sigma^2}
  }{
    \sigma^2+1.36053
  }
\right) \left(
  (-0.733996\sigma^3+1.50912\sigma^2-1.16402\sigma)
  (1-n\cdot e)^{1+\frac{1}{39\sigma^4+1}}
  +1
\right)
L_{m}(l,\sigma) = (\max(1-2\sigma,0)(1-(1-n\cdot l)^5)+\min(2\sigma,1))
((1-0.5\sigma)n\cdot l+0.5\sigma(n\cdot l)^2)
V_{d}(e,l,\sigma) = \left(
  \frac{\sigma^2}{(\sigma^2+0.09)(1.31072+0.995584(n\cdot e))}
\right) \left(
  1-(1-n\cdot l)^{\frac{
    1-0.3726732(n\cdot e^2)
  }{
    0.188566+0.388410n\cdot e
  }}
\right)
B_{p}(e,l) = \begin{cases}
1.4(n\cdot e)(n\cdot l)(e\cdot l-(n\cdot e)(n\cdot l)) & \text{if} (e\cdot l - (n\cdot e)(n\cdot l)) < 0 \\[2ex]
e\cdot l - (n\cdot e)(n\cdot l) & \text{otherwise}
\end{cases}

GGX Approximation Model

Earl Hammon 氏によって発表されたモデルです[HAM17].微小面モデルをベースにしているところは Oren-Nayar と同じです.分布関数には GGX, 幾何減衰項は Smith を使用しています.また,高速に処理するように計算方法を工夫しているようです.ただし,Oren-Nayar モデルではスーパーラフを表現できますが,このモデルでは粗さ(ラフネス)は0~1までの間です.

\begin{align}
facing &= 0.5+0.5(\vec{s}\cdot\vec{v}) \\[2ex]
rough &= facing(0.9-0.4facing)\left(\frac{0.5+\vec{n}\cdot\vec{h}}{\vec{n}\cdot\vec{h}}\right) \\[2ex]
smooth &= 1.05(1-(1-\vec{n}\cdot\vec{s})^5)(1-(1-\vec{n}\cdot\vec
{v})^5) \\[2ex]
single &= \frac{1}{\pi}lerp(smooth, rough, \alpha) \\[2ex]
multi &= 0.1159\alpha \\[2ex]
diffuse &= albedo * (single + albedo * multi)
\end{align}

最後に

勉強不足でまだまだ不明がことが多いですが,わかる範囲でまとめてみました.拡散反射BRDFについてはランバート反射を使わないケースも多くなってきているかと思います.しかし,ランバート反射もよく使われているんじゃないでしょうか.それぞれの特徴を理解して使いこなすことが大事です.今回の内容はほとんどが英語論文から読み取ったものですので,間違いはあるかもしれません.気づいた方がいましたらご連絡いただけると幸いです.また,数式が複雑なため,おかしいと思ったら原文を参照することを強く推奨します.
次回は鏡面反射BRDFですかね.ではでは.

参考文献

  • TORRANCE, K. E., SPARROW, E. M. Theory for Off-Specular Reflection From Roughness Surfaces, 1966
  • Blinn, J. F. Models of light reflection for computer synthesized pictures, 1978
  • Schlick, C. An Inexpensive BRDF Model for Physically-based Rendering, 1994
  • Oren, M., Nayar, S. K. Generalization of Lambert's Reflection Model, 1994
  • Rusinkiewicz, S. M. A New Change of Variables for Efficient BRDF Representation, 1998
  • Burley, B. Physically-Based Shading at Disney, 2012
  • [GOT14] Gotanda, Y. Designing Reflectance Models for New Consoles, 2014
  • Lagarde, S., Rousiers, C. Moving Frostbite to Physically Based Rendering 3.0, 2014
  • [HAM17] Earl Hammon, Jr. PBR Diffuse Lighting for GGX+Smith Microsurfaces, 2017
  • Fujii, Y. A tiny improvement of Oren-Nayar reflectance model, link
  • hanecci, 「マイクロファセットの分布関数 D(m) について」, link
  • Pocol, 「超雑訳 Physically Based Shading at Disney」, link