GLSLとVulkanで偏光レイトレレンダラーを作って遊ぶ

13 min read読了の目安(約12200字

リポジトリ

今回のプロジェクトのリポジトリはこちらにアップロードしてあります.Vulkanのインストールが必要です.

本記事の目的

偏光(Polarization)とは、電場および磁場の振動方向が規則的な光のことを指します.
偏光の理論の解説はこちらで非常に詳しく解説されているので,本稿では省略します.

残念ながら人間の目には入射した偏光であるかどうかを識別する能力はありません.しかし,偏光計測は従来の光計測より多くの情報量を持っており,近年偏光の情報を使うことで現実物体の形(法線)や反射率(BRDF)が従来より高精度に取れるということが分かってきました[1, 2].
その流れを受けて,有名な物理ベースレンダラーMitsuba2では,偏光を使ってレンダリングを行う機能が追加されています.

Mitsuba2は素晴らしいレンダラーであることは結果を見ても明らかですが,普段シェーダーを書く人間からするとグラフィックAPIを使ってレンダリングをしたいと思うものです.しかし,探してみた限りではそのようなソースは見つからなかったので,今回はVulkanを使って効率の良い早い偏光レイトレーシングレンダラーを書いてみようと思います.

個人的にはRGB-Dカメラの次に偏光カメラでゴリゴリに法線やBRDFを計測する未来が数年内に来ると思っているので,その際に本記事が生きたら幸いです.

無偏光レイトレと偏光レイトレの違い

まず偏光レイトレを始める前に,無偏光レイトレで解く式の復習をします.
無偏光レイトレでは,あるピクセルについて,以下のレンダリング方程式を解くことを目的としています.

\begin{aligned} L(\bm{p}) = \int_{\mathcal{P}} f(\bar{\bm{x}}) d\bar{\bm{x}} \tag{1} \end{aligned}

ここで,\bar{\bm{x}}=(\bm{p},\bm{x}_1,...,\bm{x}_d)\in \mathcal{P}は視点\bm{p}から始まって光源上の点\bm{x}_dで終わる,3D空間上で追跡可能な深さdのすべての光の経路(パス)表し,\bm{x}_iはパス上の各点の位置を表します.f(\bm{x})はパス\bar{\bm{x}}によって得られる光(放射輝度)の寄与で,

\begin{aligned} f(\bar{\bm{x}}) = \prod_{i=2}^{d-1}f_r(\bm{x}_{i-1}\leftarrow\bm{x}_{i}\leftarrow\bm{x}_{i+1}) L_e(\bm{x}_{d-1}\leftarrow\bm{x}_d) \end{aligned}

と表せます.ここで,f_r(\cdot)は点\bm{x}_{i}における双方向反射率分布関数(BRDF: Bidirectional Reflectance Distribution Function)を,L_e(\bm{x}_{d-1}\leftarrow\bm{x}_d)は光源の強さを表します.パス空間\mathcal{P}は無限にあるため,レイトレーシング(正確にはパストレーシング)では(1)式をモンテカルロ法によって解くのでした.

\begin{aligned} \braket{L(\bm{p})}_{N_s} = \cfrac{1}{N_s}\sum_{j=1}^{N_s} \cfrac{f(\bar{\bm{x}}_j)}{p(\bar{\bm{x}}_j)} \end{aligned}

ここで,p(\bar{\bm{x}}_j)はパス\bar{\bm{x}}_jをサンプルする確率密度関数を表し,\braket{L(\bm{p})}_{N_s}L(\bm{p})をサンプル数N_sで近似したモンテカルロ推定器であることを表します.この数式には,いくつかのパス\bar{\bm{x}}_jp(\bar{\bm{x}}_j)に沿ってサンプリングし,f(\bar{\bm{x}}_j)の期待値を取ればいつかは真値L(\bm{p})に収束するよねという意味が込められています.

それでは次に偏光レンダリングの説明に移ります.偏光レンダリングでは,光の偏光状態を4次元のストークスベクトル\bm{s}=[s_0, s_1, s_2, s_3]によって表現します.ストークスベクトルの各要素の意味についてはこちらを参照してください.偏光状態にある光が物体面上で反射すると,物体面の材質に依存してストークスベクトルが変化します.反射によるストークスベクトルの変化は,ミュラー行列\bm{M}\in\mathbb{R}^{4\times 4}を用いた線形変換として表されます.つまり,入射したストークスベクトルを\bm{s}_i,反射したストークスベクトルを\bm{s}_oとすると,

\bm{s}_o = \bm{M}\bm{s}_i

と表せます.偏光レイトレでは,ミュラー行列を偏光BRDF(pBRDF)として扱います.このようにストークスベクトルで偏光状態を表現すると,パス上での偏光状態の変化はミュラー行列を用いて繰り返し線形変換を施すことで簡単に追跡できます.すなわち,あるパス\bar{\bm{x}}上におけるストークスベクトルの寄与は,

\begin{aligned} \bm{f}(\bar{\bm{x}}) = \left[\prod_{i=2}^{d-1}\bm{M}(\bm{x}_{i-1}\leftarrow\bm{x}_{i}\leftarrow\bm{x}_{i+1})\right] \bm{s}(\bm{x}_{d-1}\leftarrow\bm{x}_d) \end{aligned}

となります.ここで,\bm{s}(\bm{x}_{d-1}\leftarrow\bm{x}_d)は光源の偏光状態を示します.光源にレイが当たるまで光源の情報が未知であるレイトレの場合でも,ミュラー行列を右から掛けて蓄積しておけばよいため,この偏光の表現方法はレイトレと非常に相性が良いですね.

ミュラー行列の意味や各成分の値は,こちらの記事が日本語でとても丁寧に解説されています.本稿では以降,こちらの記事を前提として話を進めていきます.

実装

レイに載せるパラメータ

無偏光レイトレの場合,レイの情報は以下のような構造体で表すことが一般的です.(実際はもっと多くの情報を載せるので構造体はもう少し複雑になります.)

struct Ray{
	vec3 o; //レイの原点
	vec3 d; //レイの方向
}

それに対し,偏光レイトレでは,レイの方向\bm{d}をz軸としたヨー方向の回転(方位角)によって偏光の状態を管理する必要があります.つまり,レイを単なる半直線として追うのでなく,向きを持った接平面と一緒に追っていかなければなりません.

そこで,本稿では上記の構造体に対し,\bm{d}に対する接線ベクトル\bm{t}を同時に保持してレイの状態を管理することにします.

struct Ray{
	vec3 o; //レイの原点
	vec3 d; //レイの方向
	vec3 t; //接線ベクトル
}

pBRDFの実装

ここがこの記事のメインとなります.メインとは言えども,先程紹介した文献[1]でpBRDFの形は与えられているので,ほぼハードコーディングすれば容易に実装できます.

先程も書いたように,以降で使う記号や行列はこちらの記事で定義済みとします.
pBRDFの本体に入る前に,pBRDFに必要なパーツのコードを紹介します.まずはフレネル屈折率T^{\perp}, T^{\parallel}です.

pBRDF.glsl
void fresnelTransparent(float cos_theta, float sin_theta, float ior, out float T_v, out float T_h) {
	float s_sin = sqrt(1.0f - (sin_theta / ior) * (sin_theta / ior));

	float nume = 4.0 * ior * cos_theta * s_sin;
	float v_denom = (cos_theta + ior * s_sin) * (cos_theta + ior * s_sin);

	T_v = nume / v_denom; \\T^{\perp}

	float h_denom = (s_sin + ior * cos_theta) * (s_sin + ior * cos_theta);
	T_h = nume / h_denom; \\T^{\parallel}
}

フレネル反射率R^{\perp}, R^{\parallel}は,以下のようにエネルギーの保存を考えれば屈折率を再利用して簡単に求められます.

pBRDF.glsl
void fresnelRefrect(float T_v, float T_h, out float R_v, out float R_h) {
	R_v = 1.0f - T_v;
	R_h = 1.0f - T_h;
}

次は偏光座標回転行列\bm{C}(\phi)です.これはz軸に対する回転行列で,簡単なハードコードで実装できます.

pBRDF.glsl
mat4 rot_z(float cos_theta, float sin_theta) {
	mat4 a;
	a[0][0] = 1.0f, a[0][1] = 0.0f, a[0][2] = 0.0f, a[0][3] = 0.0f;
	a[1][0] = 0.0f, a[1][1] = cos_theta, a[1][2] = sin_theta, a[1][3] = 0.0f;
	a[2][0] = 0.0f, a[2][1] = -sin_theta, a[2][2] = cos_theta, a[2][3] = 0.0f;
	a[3][0] = 0.0f, a[3][1] = 0.0f, a[3][2] = 0.0f, a[3][3] = 1.0f;
	return a;
}

最後にフレネル屈折(or反射)行列\bm{F}(\theta, \delta)です.これも定義通りに代入していきます.

pBRDF.glsl
// vert, horizは垂直,平行方向のフレネル屈折(or反射)率
mat4 fresnel_mat(float vert, float horiz, float delta) {
	float alpha = (vert + horiz) / 2.0f;
	float beta = (vert - horiz) / 2.0f;
	float chi = sqrt(vert * horiz);

	float c = cos(delta);
	float s = sin(delta);
	mat4 a;
	a[0][0] = alpha, a[0][1] = beta,  a[0][2] = 0.0f,   a[0][3] = 0.0f;
	a[1][0] = beta,  a[1][1] = alpha, a[1][2] = 0.0f,   a[1][3] = 0.0f;
	a[2][0] = 0.0f,  a[2][1] = 0.0f,  a[2][2] = chi*c,  a[2][3] = chi*s;
	a[3][0] = 0.0f,  a[3][1] = 0.0f,  a[3][2] = -chi*s, a[3][3] = chi*c;
	return a;
}

それでは,ある物体面上でのミュラー行列\bm{M}の中身を見ていきます.本稿では,ミュラー行列のモデルとして以下のものを使います.

\bm{M} = \rho (\bm{n}\cdot\vec{\omega}_i) (\bm{P}^d+\bm{P}^s)

ここで,\rhoはアルベド,(\bm{n}\cdot\vec{\omega}_i)は法線と入射方向による余弦項,\bm{P}^d, \bm{P}^sはそれぞれ偏光拡散(鏡面)反射行列を表します.実際のコードではこんな感じです.

pBRDF.glsl
void pBsdfEval(out mat4 muller,out vec3 albedo, in Ray ray, in vec3 N, in vec3 T, in vec3 B, in Material mat, in vec3 bsdfDir, vec3 bsdfTangent) {
	vec3 V = -ray.direction;
	vec3 L = bsdfDir;

	mat4 diffuse;
	pBRDF_diffuse(diffuse, bsdfTangent, ray.tangent, T, L, V, N, mat);

	mat4 specular;
	pBRDF_specular(specular, bsdfTangent, ray.tangent, T, L, V, N, mat);

	muller = diffuse + specular;
	albedo = mat.albedo.xyz;
}

アルベドは扱いが少々厄介なので,パス上でのミュラー行列の処理が終わった後にrgenシェーダーで最後に処理するようにしています.

pathtrace.rgen
void main(){
//一部抜粋
 vec4 storks=vec4(0.0);
 for(int i=0;i<3;i++){
     storks = prd.albedo[i] * prd.storks;
     hitStorks[i] += storks;
 }
}      

pBRDFの引数として与えるマテリアルですが,本稿では以下のような構造体にてパラメータの定義を行っています.

brdf.glsl
struct Material { 
	vec4 albedo;      //アルベド
	float metallic;   //金属性
	float roughness;  //粗さ
	vec3 f0;          //0度でのフレネル反射率(今回は未使用)
	float anisotropy; //異方性(今回は未使用)
	float ior;        //空気との相対屈折率
};

次に,偏光拡散反射行列\bm{P}^dの中身を見ていきます.

pBRDF.glsl
vec3 projectToPlane(vec3 a, vec3 n) {
    return normalize(a - dot(a, n) * n);
}
void pBRDF_diffuse(out mat4 result, vec3 T_i, vec3 T_o, vec3 T, vec3 L, vec3 V, vec3 N,
	in Material mat)
{
	float T_vi, T_hi, T_vo, T_ho;

	float cos_theta_i = clamp(dot(L, N),0.001,1.0); //天頂角
	float sin_theta_i = sqrt(1.0f - cos_theta_i * cos_theta_i);
	fresnelTransparent(cos_theta_i, sin_theta_i, mat.ior, T_vi, T_hi);

	float cos_theta_o = clamp(dot(V, N), 0.001, 1.0); //天頂角
	float sin_theta_o = sqrt(1.0f - cos_theta_o * cos_theta_o);
	fresnelTransparent(cos_theta_o, sin_theta_o, mat.ior, T_vo, T_ho);

	float phi_o = acos(dot(T_o, projectToPlane(N, V))); //方位角
	float phi_i = acos(dot(T_i, projectToPlane(N, L))); //方位角

	mat4 rot_o = rot_z(cos(2.0f * phi_o), sin(2.0f * phi_o));
	mat4 rot_i = rot_z(cos(-2.0f * phi_i), sin(-2.0f * phi_i));

	mat4 F_i = fresnel_mat(T_vi, T_hi, 0.0f);
	mat4 F_o = fresnel_mat(T_vo, T_ho, 0.0f);

	mat4 one;
	one_mat(one); //無偏光化行列

	mat4 i = F_i * rot_i;
	mat4 o = rot_o * F_o;

	result = o * one * i;
}

\theta_i\theta_oは入反射するレイと物体面上での座標系の間の天頂角のずれを表し,\phi_i, \phi_oは方位角のずれを表します.方位角のずれについては,物体面上での法線を,レイの方向を法線とする平面に射影してから,レイの接線ベクトルとの内積を取ることで求められます.

最後に,偏光鏡面反射行列\bm{P}^sについて見ていきます.

pBRDF.glsl
void pBRDF_specular(out mat4 result, vec3 T_i, vec3 T_o, vec3 T, vec3 L, vec3 V, vec3 N, in Material mat)
{
	vec3  H = normalize(L + V);
	float T_v, T_h, R_v, R_h;
	float cos_theta_d = clamp(dot(L, H), 0.001, 1.0); //天頂角
	float sin_theta_d = sqrt(1.0f - cos_theta_d * cos_theta_d);
	fresnelTransparent(cos_theta_d, sin_theta_d, mat.ior, T_v, T_h);
	fresnelRefrect(T_v, T_h, R_v, R_h);

	float phy_o = acos(dot(T_o, projectToPlane(H, V))); //方位角
	float phy_i = acos(dot(T_i, projectToPlane(H, L))); //方位角

	float theta_b = atan(mat.ior); //ブリュースター角
	float delta = (mat.metallic > 0.5f) && (acos(cos_theta_d) < theta_b) ? M_PI : 0.0f;

	mat4 rot_o = rot_z(cos(2.0f * phy_o), sin(2.0f * phy_o));
	mat4 rot_i = rot_z(cos(-2.0f * phy_i), sin(-2.0f * phy_i));

	mat4 F = fresnel_mat(T_v, T_h, delta);
	F = rot_o * F * rot_i;

	float alphaRoughness = max(0.001, mat.roughness);
	float V_g = V_GGX(abs(dot(N, V)), abs(dot(N, L)), alphaRoughness);
	float D_g = D_GGX(abs(dot(N, H)), max(0.001, alphaRoughness));

	result = V_g * D_g * F;
}

偏光鏡面反射行列はマイクロファセット理論をもとに出発しているため,天頂角や方位角はハーフベクトルから求めていることに注意してください.V_GGXはMasking-Shadowing項,D_GGXは法線分布関数を表し,無偏光BRDFと同一のものを使い回すことができます.

レイ方向のサンプリング

レイトレをする上で避けられない,パスを構築するためのレイ方向のサンプリングですが,本稿のpBRDFでは無偏光BRDFと同一の法線分布関数を採用しているため,無偏光BRDFと同一の重点的サンプリング手法を使い回すことができます.よって,本稿では詳細なサンプリング手法については割愛します.

偏光状態の可視化

残念ながら人間の目には偏光状態を識別する能力はなく,ストークスベクトル\bm{s}=[s_0, s_1, s_2, s_3]の成分もそのままではどのような物理的意味を持つのか分かりにくいため,直感的に偏光状態を可視化する手法が提案されています[3, 4].本稿ではそれらの手法に倣った可視化を実装しています.

  • DoP (Degree of Polarization): \mathrm{DoP} := \frac{\sqrt{s_1^2+s_2^2+s_3^2}}{s_0}により定義される偏光状態の強さ.
  • AoLP (Angle of Linear Polarization): \mathrm{AoLP} := \frac{1}{2}\tan^{-1}\left(\frac{s_2}{s_1}\right)により定義される,直線偏光における偏光面と基準面の間の角度.
  • ToP (Type of Polarization): \mathrm{ToP} := \frac{1}{2}\tan^{-1}\left(\frac{s_3}{\sqrt{s_1^2+s_2^2}}\right)により定義される偏光の種類. 0が直線偏光寄り,\piが円偏光寄りであることを示す.

結果

実験環境

  • CPU: Intel Core-i9 9900K
  • GPU: NVIDIA GeForce RTX3080
  • メモリ: 32GB
  • OS: Windows10

Cornell Box

まずはおなじみのコーネルボックスです.各結果は\bm{s}=[...]のようなタイトルをつけており,これは入力として用いた光源のストークスベクトルを表しています.

\bm{s}=[1.0, 0.0, 0.0, 0.0]

\bm{s}=[1.0, 1.0, 0.0, 0.0]

\bm{s}=[1.0, 0.0, -1.0, 0.0]

\bm{s}=[1.0, 0.0, 0.0, 1.0]

案の定ではありますが,拡散成分が大きいと偏光成分はほぼ打ち消されてDoPはほぼ0になります,それに対し,粗さ成分が低い金属球においては照明の一次反射をうまく拾って偏光状態が保存されていることが顕著に表れました.

Sponza

流石に拡散成分が多すぎるコーネルボックスではネタとして何も面白くなくてマズいので,これまたおなじみのSponzaシーンを偏光成分が強く出るように,粗さを下げたり金属製の高いオブジェクトを配置して改造したシーンを用意しました.

\bm{s}=[1.0, 0.0, 0.0, 0.0]

\bm{s}=[1.0, -1.0, 0.0, 0.0]

\bm{s}=[1.0, 0.0, 1.0, 0.0]

\bm{s}=[1.0, 0.0, 0.0, -1.0]

意図的に鏡面反射が強く出るように改造しただけあり,鏡面反射成分が強い領域ではそこそこ強い偏光状態が観察され,本稿の目的に対してそこそこ面白い結果が得られました.

おまけ

リアルタイムで動かしているデモです.YouTubeにリンクします.

まとめ

本稿では,偏光の情報を使って物体の情報を取るという流れが来ると勝手に踏んで偏光状態のビューワーをglslとvulkanを用いて実装しました.(ほとんどglslのコードしか載せていませんが...)

参考文献

[1] Baek et al. 2018. Simultaneous Acquisition of Polarimetric SVBRDF and Normals. ACM SIGGRAPH Asia 2018. ↩
[2] Baek et al. 2020. Image-based acquisition and modeling of polarimetric reflectance. ACM SIGGRAPH 2020. ↩
[3] Wilkie et al. 2010. A Standardised Polarisation Visualisation for Images. Proceedings of the 26th Spring Conference on Computer Graphics.
[4] Li et al. 2014. A method to calculate Stokes parameters and angle of polarization of skylight from polarized CIMEL sun/sky radiometers. Journal of Quantitative Spectroscopy and Radiative Transfer. Vol.149.