GLSLとVulkanで偏光レイトレレンダラーを作って遊ぶ
リポジトリ
今回のプロジェクトのリポジトリはこちらにアップロードしてあります．Vulkanのインストールが必要です．
本記事の目的
偏光（Polarization）とは、電場および磁場の振動方向が規則的な光のことを指します．
偏光の理論の解説はこちらで非常に詳しく解説されているので，本稿では省略します．
残念ながら人間の目には入射した偏光であるかどうかを識別する能力はありません．しかし，偏光計測は従来の光計測より多くの情報量を持っており，近年偏光の情報を使うことで現実物体の形(法線)や反射率(BRDF)が従来より高精度に取れるということが分かってきました[1, 2]．
その流れを受けて，有名な物理ベースレンダラーMitsuba2では，偏光を使ってレンダリングを行う機能が追加されています．
Mitsuba2は素晴らしいレンダラーであることは結果を見ても明らかですが，普段シェーダーを書く人間からするとグラフィックAPIを使ってレンダリングをしたいと思うものです．しかし，探してみた限りではそのようなソースは見つからなかったので，今回はVulkanを使って効率の良い早い偏光レイトレーシングレンダラーを書いてみようと思います．
個人的にはRGB-Dカメラの次に偏光カメラでゴリゴリに法線やBRDFを計測する未来が数年内に来ると思っているので，その際に本記事が生きたら幸いです．
無偏光レイトレと偏光レイトレの違い
まず偏光レイトレを始める前に，無偏光レイトレで解く式の復習をします．
無偏光レイトレでは，あるピクセルについて，以下のレンダリング方程式を解くことを目的としています．
ここで，
と表せます．ここで，
ここで，
それでは次に偏光レンダリングの説明に移ります．偏光レンダリングでは，光の偏光状態を4次元のストークスベクトル
と表せます．偏光レイトレでは，ミュラー行列を偏光BRDF(pBRDF)として扱います．このようにストークスベクトルで偏光状態を表現すると，パス上での偏光状態の変化はミュラー行列を用いて繰り返し線形変換を施すことで簡単に追跡できます．すなわち，あるパス
となります．ここで，
ミュラー行列の意味や各成分の値は，こちらの記事が日本語でとても丁寧に解説されています．本稿では以降，こちらの記事を前提として話を進めていきます．
実装
レイに載せるパラメータ
無偏光レイトレの場合，レイの情報は以下のような構造体で表すことが一般的です．（実際はもっと多くの情報を載せるので構造体はもう少し複雑になります．）
struct Ray{
vec3 o; //レイの原点
vec3 d; //レイの方向
}
それに対し，偏光レイトレでは，レイの方向
そこで，本稿では上記の構造体に対し，
struct Ray{
vec3 o; //レイの原点
vec3 d; //レイの方向
vec3 t; //接線ベクトル
}
pBRDFの実装
ここがこの記事のメインとなります．メインとは言えども，先程紹介した文献[1]でpBRDFの形は与えられているので，ほぼハードコーディングすれば容易に実装できます．
先程も書いたように，以降で使う記号や行列はこちらの記事で定義済みとします．
pBRDFの本体に入る前に，pBRDFに必要なパーツのコードを紹介します．まずはフレネル屈折率
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}
}
フレネル反射率
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;
}
次は偏光座標回転行列
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反射)行列
// 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;
}
それでは，ある物体面上でのミュラー行列
ここで，
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シェーダーで最後に処理するようにしています．
void main(){
//一部抜粋
vec4 storks=vec4(0.0);
for(int i=0;i<3;i++){
storks = prd.albedo[i] * prd.storks;
hitStorks[i] += storks;
}
}
pBRDFの引数として与えるマテリアルですが，本稿では以下のような構造体にてパラメータの定義を行っています．
struct Material {
vec4 albedo; //アルベド
float metallic; //金属性
float roughness; //粗さ
vec3 f0; //0度でのフレネル反射率(今回は未使用)
float anisotropy; //異方性（今回は未使用）
float ior; //空気との相対屈折率
};
次に，偏光拡散反射行列
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;
}
最後に，偏光鏡面反射行列
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と同一の重点的サンプリング手法を使い回すことができます．よって，本稿では詳細なサンプリング手法については割愛します．
偏光状態の可視化
残念ながら人間の目には偏光状態を識別する能力はなく，ストークスベクトル
- 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}=[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はほぼ
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.
Discussion