物理ベースレンダリング入門 その③ - PBRシェーダをUnityで実装する

Unity Shader Rendering

物理ベースレンダリングに入門するための連載の第3回です。
最終回の今回はPBRのシェーダをUnityで実装します。

Unity2019.3.3

はじめに

この連載では物理ベースレンダリングに関する基礎的な知識をまとめています。
第1回では物理ベースレンダリングの概念とその背後にある考え方についてまとめました。

light11.hatenadiary.com

また第2回ではBRDFについてまとめました。

light11.hatenadiary.com

第3回の本記事ではこれらの内容を踏まえて実際にUnityでPBRのシェーダを実装します。
内容は上に挙げた二つの記事の知識を前提としますので、必要に応じて参照してください。

なおBRDFによるシェーディングはイチから実装しますが、間接光に関してはUnityの関数を使って適用します。
また、簡易化のため色空間はリニアのみ想定し、ライトはディレクショナルライト一つのみを処理します。
レンダリングパイプラインに関してはSRPではなく従来のパイプラインを前提として記述します。

鏡面反射BRDF

それではまず鏡面反射BRDFから実装します。
BRDFのモデルとしては第2回で解説したCook-Torranceモデルを使います。

{\displaystyle f_{r}(v,l)=\frac{D(h,\alpha )G(v,l, \alpha)F(v,h,f_{0})}{4(n\cdot v)(n\cdot l)} \tag{1} }

ただし今回は幾何減衰項にHeight-Correlated Smithモデルを採用するのでBRDFは以下のように変形します。

{\displaystyle f_{r}(v,l)=D(h,\alpha )V(v,l,\alpha)F(v,h,f_{0}) \tag{2} }

以下でD項、V項、F項をそれぞれ個別に関数化していきます。

法線分布関数（D項）を実装する

D項にはGGX（Trowbrdige-Retiz）を使います。

{\displaystyle D_{GGX}(h,\alpha) = \frac{\alpha ^{2}}{\pi ( (n \cdot h)^{2}(\alpha ^{2}-1)+1)^{2}} \tag{3} }

これをソースコードで表すと以下のようになります。

inline float D_GGX(float ndoth, float alpha)
{
    float a2 = alpha * alpha;
    float d = (ndoth * a2 - ndoth) * ndoth + 1.0f;
    return a2 / (d * d + 0.0000001) * UNITY_INV_PI;
}

UNITY_INV_PIはUnityCG.cgincに定義されている、 1/{\pi}を表す値です。

またUnityの実装は上記のような元の式そのままの感じですが、Filament*1ではもう少し最適化を行っているようなので、これを参考にもう少しこの関数を書き換えていきます。
まず上式を以下のように変形します。

inline float D_GGX(float ndoth, float perceptualRoughness) {
    float a = ndoth * perceptualRoughness;
    float k = perceptualRoughness / (1.0 - ndoth * ndoth + a * a);
    return k * k * UNITY_INV_PI;
}

ここで半精度の浮動小数点を使えばより効率化できますが、
このとき 1-(n \cdot h)^{2}の計算を半精度で行うことで以下の二つの問題が浮上します。

  • (n \cdot h)^{2}の計算結果が1に近いとき（すなわちハイライト部分）に桁落ちが発生する
  • n \cdot hが1付近で十分な精度を持たない

これを解決するためにラグランジュの恒等式を使います。

{\displaystyle \left | a \times b \right |^{2}=\left| a \right|^{2}\left| b \right|^{2}-(a \cdot b)^{2} \tag{4} }

nと hはそれぞれ単位ベクトルなので \left| n \times h \right|^{2}=1-(n \cdot h)^{2}となります。
つまりこれを使えば 1-(n \cdot h)^{2}が外積を用いて半精度小数点で計算できます。
以下はこれを適用するとD項は以下の関数として表せます。

inline half D_GGX(half perceptualRoughness, half ndoth, half3 normal, half3 halfDir) {
    half3 ncrossh = cross(normal, halfDir);
    half a = ndoth * perceptualRoughness;
    half k = perceptualRoughness / (dot(ncrossh, ncrossh) + a * a);
    half d = k * k * UNITY_INV_PI;
    return min(d, 65504.0);
}

ちなみにperceptualRoughnessとはユーザが入力した値としてのラフネスです。
BRDFで粗さを扱うときにはこれを二乗して \alphaとすることが多いですが、
これと明確に区別するためにこのような変数名にしています。

幾何減衰項(V項)を実装する

次に幾何減衰項を実装します。
幾何減衰項にはHeight-Correlated Smithモデルを使います。
このモデルは以下の式で表されます。

{\displaystyle V(v,l,\alpha)= \frac{0.5}{\Lambda _{l}+\Lambda_{v}} \tag{5} }

{ \displaystyle \Lambda _{l}=n \cdot v \sqrt{(n \cdot l)^{2}(1.0-\alpha ^{2})+\alpha ^{2}} \tag{6} }

これを基にしてシェーダを書くと以下のようになります。

inline float V_SmithGGXCorrelated(float ndotl, float ndotv, float alpha)
{
    half alpha2 = alpha * alpha;
    half lambdaV = ndotl * sqrt((-ndotv * alpha2 + ndotv) * ndotv + alpha2);
    half lambdaL = ndotv * sqrt((-ndotl * alpha2 + ndotl) * ndotl + alpha2);
    return 0.5f / (lambdaV + lambdaL + 0.0001);
}

この計算式にはsqrt関数が入っていて重いので、もう少し最適化します。
式（6）を見ると以下のことが見て取れます。

  • 平方根に含まれるすべての項が平方
  • それらすべてが0～1の値を取る

これらから、以下のざっくりとした近似式が得られます。

{\displaystyle \Lambda _{l}=n \cdot v(n \cdot l(1- \alpha)+\alpha) \tag{7} }

この近似は数学的には誤りですが、今回のケースでは費用対効果を考えると十分です。
これを元にシェーダを実装すると以下のようになります。

inline float V_SmithGGXCorrelated(float ndotl, float ndotv, float alpha)
{
    float lambdaV = ndotl * (ndotv * (1 - alpha) + alpha);
    float lambdaL = ndotv * (ndotl * (1 - alpha) + alpha);
    return 0.5f / (lambdaV + lambdaL + 0.0001);
}

UnityのStandardシェーダでもFilamentでもこのように実装されているようです。

フレネル項(F項)を実装する

フレネル反射はSchlickの近似式を使用します。

{\displaystyle F_{Schlick}(v,h,f_{0},f_{90})=f_{0}+(f_{90}-f_{0})(1-v \cdot h)^{5} \tag{8} }

シェーダは以下のように記述します。
なお今回のシェーダでは f_{90}=1としています。

inline half3 F_Schlick(half3 f0, half cos)
{
    return f0 + (1 - f0) * pow(1 - cos, 5);
}

拡散反射モデルを実装する

それではまず拡散反射の実装を行います。
拡散反射モデルはDisneyのものを使います。
式としては以下のようになります。

{\displaystyle f_{d}(v,l)=\frac{\sigma}{\pi}F_{Schlick}(n,l,1,f_{90})F_{Schlick}(n,v,1,f_{90}) \tag{9} }

{\displaystyle f_{90}=0.5+2 \alpha cos^{2}(\theta_{d}) \tag{10} }

これをシェーダで実装すると以下のようになります。

inline half Fd_Burley(half ndotv, half ndotl, half ldoth, half roughness)
{
    half fd90 = 0.5 + 2 * ldoth * ldoth * roughness;
    half lightScatter = (1 + (fd90 - 1) * pow(1 - ndotl, 5));
    half viewScatter = (1 + (fd90 - 1) * pow(1 - ndotv, 5));

    half diffuse = lightScatter * viewScatter / UNITY_PI;
    return diffuse;
}

UNITY_PIはUnityCG.cgincに定義されている、 {\pi}を表す値です。
実装の都合上、ディフューズアルベドはこの関数を使う側で乗算します。

また、UnityのStandardシェーダではレガシーなシェーダとの見え方の互換性を取るために {\pi}による除算を省いています。

inline half Fd_Burley(half ndotv, half ndotl, half ldoth, half roughness)
{
    half fd90 = 0.5 + 2 * ldoth * ldoth * roughness;
    half lightScatter = (1 + (fd90 - 1) * pow(1 - ndotl, 5));
    half viewScatter = (1 + (fd90 - 1) * pow(1 - ndotv, 5));

    half diffuse = lightScatter * viewScatter;
    // 本来はこのDiffuseをπで割るべきだけどUnityではレガシーなライティングとの互換性を保つため割らない
    //diffuse /= UNITY_PI;
    return diffuse;
}

BRDFとしては正しくないですが、本記事ではこのStandardShaderに合わせて実装を採用しています。

BRDFを使って色を求める関数を作る

さてそれではここまで作った関数を使ってBRDFを実装します。
シェーダは以下のように記述します。

half4 BRDF(half3 albedo, half metallic, half perceptualRoughness, float3 normal, float3 viewDir, float3 lightDir, float3 lightColor)
{
    float3 halfDir = normalize(lightDir + viewDir);
    half ndotv = abs(dot(normal, viewDir));
    float ndotl = max(0, dot(normal, lightDir));
    float ndoth = max(0, dot(normal, halfDir));
    half ldoth = max(0, dot(lightDir, halfDir));
    half reflectivity = lerp(_DielectricF0, 1, metallic);
    half3 f0 = lerp(_DielectricF0, albedo, metallic);
    
    // Diffuse
    half diffuseTerm = Fd_Burley(ndotv, ndotl, ldoth, perceptualRoughness) * ndotl;
    half3 diffuse = albedo * (1 - reflectivity) * lightColor * diffuseTerm;
    
    // Specular
    float alpha = perceptualRoughness * perceptualRoughness;
    float V = V_SmithGGXCorrelated(ndotl, ndotv, alpha);
    float D = D_GGX(perceptualRoughness, ndotv, normal, halfDir);
    float3 F = F_Schlick(f0, ldoth); // マイクロファセットベースのスペキュラBRDFではcosはldothが使われる
    float3 specular = V * D * F * ndotl * lightColor;
    // 本来はSpecularにπを掛けるべきではないが、Unityではレガシーなライティングとの互換性を保つため、Diffuseを割らずにSpecularにPIを掛ける
    specular *= UNITY_PI;
    specular = max(0, specular);

    half3 color = diffuse + specular;
    return half4(color, 1);
}

DiffuseについてはアルベドにDisneyDiffuseで求めた値を乗算することで反射率が求められます。
また、鏡面反射成分が大きいほど拡散反射は小さくなるので1 - reflectivityを乗算しています。

SupecularについてはD,V,F項の関数を使って求め、それぞれを乗算します。
フレネル項にはldothを渡している点に注意が必要です。

また、拡散反射の節でUnityのレガシーなシェーダとの整合性を取るために \piの除算を省略しましたが、
さらにこれと整合性を取るために鏡面反射成分に本来乗算しない \piを乗算しています。

PBRシェーダ全文

ここまでのシェーダをまとめると以下のようになります。

Shader "Example"
{
    Properties
    {
        _Color("Color", Color) = (1,1,1,1)
        _MainTex("Albedo", 2D) = "white" {}
        _Roughness("Roughness", Range(0.0, 1.0)) = 0.5
        _Metallic("Metallic", Range(0.0, 1.0)) = 0.0
    }

    SubShader
    {
        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #pragma multi_compile_fwdbase

            #include "UnityCG.cginc"

            struct appdata
            {
                float4 vertex : POSITION;
                half2 texcoord : TEXCOORD0;
                half3 normal : NORMAL;
            };

            struct v2f
            {
                float4 pos : SV_POSITION;
                half2 uv : TEXCOORD0;
                float3 worldPos : TEXCOORD1;
                half3 worldNormal : TEXCOORD2;
                half3 viewDir : TEXCOORD3;
            };

            float3 _Color;
            sampler2D _MainTex;
            float4 _MainTex_ST;
            float _Metallic;
            float _Roughness;
            float3 _LightColor0;

            // 誘電体の反射率（F0）は4%とする
            #define _DielectricF0 0.04

            inline half Fd_Burley(half ndotv, half ndotl, half ldoth, half roughness)
            {
                half fd90 = 0.5 + 2 * ldoth * ldoth * roughness;
                half lightScatter = (1 + (fd90 - 1) * pow(1 - ndotl, 5));
                half viewScatter = (1 + (fd90 - 1) * pow(1 - ndotv, 5));

                half diffuse = lightScatter * viewScatter;
                // 本来はこのDiffuseをπで割るべきだけどUnityではレガシーなライティングとの互換性を保つため割らない
                //diffuse /= UNITY_PI;
                return diffuse;
            }
            
            inline float V_SmithGGXCorrelated(float ndotl, float ndotv, float alpha)
            {
                float lambdaV = ndotl * (ndotv * (1 - alpha) + alpha);
                float lambdaL = ndotv * (ndotl * (1 - alpha) + alpha);

                return 0.5f / (lambdaV + lambdaL + 0.0001);
            }

            inline half D_GGX(half perceptualRoughness, half ndoth, half3 normal, half3 halfDir) {
                half3 ncrossh = cross(normal, halfDir);
                half a = ndoth * perceptualRoughness;
                half k = perceptualRoughness / (dot(ncrossh, ncrossh) + a * a);
                half d = k * k * UNITY_INV_PI;
                return min(d, 65504.0h);
            }

            inline half3 F_Schlick(half3 f0, half cos)
            {
                return f0 + (1 - f0) * pow(1 - cos, 5);
            }

            half4 BRDF(half3 albedo, half metallic, half perceptualRoughness, float3 normal, float3 viewDir, float3 lightDir, float3 lightColor)
            {
                float3 halfDir = normalize(lightDir + viewDir);
                half ndotv = abs(dot(normal, viewDir));
                float ndotl = max(0, dot(normal, lightDir));
                float ndoth = max(0, dot(normal, halfDir));
                half ldoth = max(0, dot(lightDir, halfDir));
                half reflectivity = lerp(_DielectricF0, 1, metallic);
                half3 f0 = lerp(_DielectricF0, albedo, metallic);
    
                // Diffuse
                half diffuseTerm = Fd_Burley(ndotv, ndotl, ldoth, perceptualRoughness) * ndotl;
                half3 diffuse = albedo * (1 - reflectivity) * lightColor * diffuseTerm;
    
                // Specular
                float alpha = perceptualRoughness * perceptualRoughness;
                float V = V_SmithGGXCorrelated(ndotl, ndotv, alpha);
                float D = D_GGX(perceptualRoughness, ndotv, normal, halfDir);
                float3 F = F_Schlick(f0, ldoth); // マイクロファセットベースのスペキュラBRDFではcosはldothが使われる
                float3 specular = V * D * F * ndotl * lightColor;
                // 本来はSpecularにπを掛けるべきではないが、Unityではレガシーなライティングとの互換性を保つため、Diffuseを割らずにSpecularにPIを掛ける
                specular *= UNITY_PI;
                specular = max(0, specular);

                half3 color = diffuse + specular;
                return half4(color, 1);
            }
            
            v2f vert(appdata v)
            {
                v2f o = (v2f)0;
                o.pos = UnityObjectToClipPos(v.vertex);
                o.uv = TRANSFORM_TEX(v.texcoord, _MainTex);
                o.worldPos = mul(unity_ObjectToWorld, v.vertex).xyz;
                o.worldNormal = UnityObjectToWorldNormal(v.normal);
                o.viewDir = UnityWorldSpaceViewDir(o.worldPos);

                return o;
            }

            half4 frag(v2f i) : SV_Target
            {
                half3 albedo = tex2D(_MainTex, i.uv).rgb * _Color.rgb;
                half metallic = _Metallic;
                half perceptualRoughness = _Roughness;

                i.worldNormal = normalize(i.worldNormal);
                i.viewDir = normalize(i.viewDir);

                half4 c = BRDF(albedo, metallic, perceptualRoughness, i.worldNormal, i.viewDir, _WorldSpaceLightPos0.xyz, _LightColor0.rgb);
                return c;
            }

            ENDCG
        }
    }
}

ここまでのレンダリング結果

さてここで一度レンダリング結果を見てみます。

f:id:halya_11:20200305213232p:plain:w600
ここまでのレンダリング結果

上段は非金属、下段は金属です。
また左から右に向けてラフネスの値を上げています。

直接光に関してはうまく実装できていそうです。

間接光を反映する

最後に間接光を反映します。
少々長くなりますが、まずシェーダ全文を載せます。

Shader "Example3"
{
    Properties
    {
        _Color("Color", Color) = (1,1,1,1)
        _MainTex("Albedo", 2D) = "white" {}
        _Roughness("Roughness", Range(0.0, 1.0)) = 0.5
        _Metallic("Metallic", Range(0.0, 1.0)) = 0.0
    }

    SubShader
    {
        Pass
        {
            // SH求めるのに必要
            Tags{ "LightMode" = "ForwardBase" }

            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #pragma multi_compile_fwdbase

            #include "UnityCG.cginc"
            #include "UnityStandardUtils.cginc"

            struct appdata
            {
                float4 vertex : POSITION;
                half2 texcoord : TEXCOORD0;
                half3 normal : NORMAL;
            };

            struct v2f
            {
                float4 pos : SV_POSITION;
                half2 uv : TEXCOORD0;
                float3 worldPos : TEXCOORD1;
                half3 worldNormal : TEXCOORD2;
                half3 viewDir : TEXCOORD3;
                half4 ambient : TEXCOORD4;
            };

            float3 _Color;
            sampler2D _MainTex;
            float4 _MainTex_ST;
            float _Metallic;
            float _Roughness;
            float3 _LightColor0;

            #define _DielectricF0 0.04

            inline half Fd_Burley(half ndotv, half ndotl, half ldoth, half roughness)
            {
                half fd90 = 0.5 + 2 * ldoth * ldoth * roughness;
                half lightScatter = (1 + (fd90 - 1) * pow(1 - ndotl, 5));
                half viewScatter = (1 + (fd90 - 1) * pow(1 - ndotv, 5));

                half diffuse = lightScatter * viewScatter;
                return diffuse;
            }
            
            inline float V_SmithGGXCorrelated(float ndotl, float ndotv, float alpha)
            {
                float lambdaV = ndotl * (ndotv * (1 - alpha) + alpha);
                float lambdaL = ndotv * (ndotl * (1 - alpha) + alpha);

                return 0.5f / (lambdaV + lambdaL + 0.0001);
            }

            inline half D_GGX(half perceptualRoughness, half ndoth, half3 normal, half3 halfDir) {
                half3 ncrossh = cross(normal, halfDir);
                half a = ndoth * perceptualRoughness;
                half k = perceptualRoughness / (dot(ncrossh, ncrossh) + a * a);
                half d = k * k * UNITY_INV_PI;
                return min(d, 65504.0h);
            }

            inline half3 F_Schlick(half3 f0, half cos)
            {
                return f0 + (1 - f0) * pow(1 - cos, 5);
            }
            
            half4 BRDF(half3 albedo, half metallic, half perceptualRoughness, float3 normal, float3 viewDir, float3 lightDir, float3 lightColor, half3 indirectDiffuse, half3 indirectSpecular)
            {
                float3 halfDir = normalize(lightDir + viewDir);
                half ndotv = abs(dot(normal, viewDir));
                float ndotl = max(0, dot(normal, lightDir));
                float ndoth = max(0, dot(normal, halfDir));
                half ldoth = max(0, dot(lightDir, halfDir));
                half reflectivity = lerp(_DielectricF0, 1, metallic);
                half3 f0 = lerp(_DielectricF0, albedo, metallic);
                
                half diffuseTerm = Fd_Burley(ndotv, ndotl, ldoth, perceptualRoughness) * ndotl;
                half3 diffuse = albedo * (1 - reflectivity) * lightColor * diffuseTerm;
                // Indirect Diffuse
                diffuse += albedo * (1 - reflectivity) * indirectDiffuse;
                
                float alpha = perceptualRoughness * perceptualRoughness;
                float V = V_SmithGGXCorrelated(ndotl, ndotv, alpha);
                float D = D_GGX(perceptualRoughness, ndotv, normal, halfDir);
                float3 F = F_Schlick(f0, ldoth); 
                float3 specular = V * D * F * ndotl * lightColor;
                specular *= UNITY_PI;
                specular = max(0, specular);
                
                // Indirect Specular
                half surfaceReduction = 1.0 / (alpha * alpha + 1.0);
                half f90 = saturate((1 - perceptualRoughness) + reflectivity);
                specular += surfaceReduction * indirectSpecular * lerp(f0, f90, pow(1 - ndotv, 5));
                
                half3 color = diffuse + specular;
                return half4(color, 1);
            }
            
            v2f vert(appdata v)
            {
                v2f o = (v2f)0;
                o.pos = UnityObjectToClipPos(v.vertex);
                o.uv = TRANSFORM_TEX(v.texcoord, _MainTex);
                o.worldPos = mul(unity_ObjectToWorld, v.vertex).xyz;
                o.worldNormal = UnityObjectToWorldNormal(v.normal);
                o.viewDir = UnityWorldSpaceViewDir(o.worldPos);

                // SH
                o.ambient.rgb = ShadeSHPerVertex(o.worldNormal, o.ambient.rgb);

                return o;
            }

            half4 frag(v2f i) : SV_Target
            {
                half3 albedo = tex2D(_MainTex, i.uv).rgb * _Color.rgb;
                half metallic = _Metallic;
                half perceptualRoughness = _Roughness;

                i.worldNormal = normalize(i.worldNormal);
                i.viewDir = normalize(i.viewDir);

                // Indirect Diffuse
                half3 indirectDiffuse = ShadeSHPerPixel(i.worldNormal, i.ambient, i.worldPos);
                
                // roughnessに対応する鏡面反射のミップマップレベルを求める
                half3 reflDir = reflect(-i.viewDir, i.worldNormal);
                half mip = perceptualRoughness * (1.7 - 0.7 * perceptualRoughness);
                // 間接光の鏡面反射（リフレクションプローブのブレンドとかは考慮しない）
                mip *= UNITY_SPECCUBE_LOD_STEPS;
                half4 rgbm = UNITY_SAMPLE_TEXCUBE_LOD(unity_SpecCube0, reflDir, mip);
                half3 indirectSpecular = DecodeHDR(rgbm, unity_SpecCube0_HDR);

                half4 c = BRDF(albedo, metallic, perceptualRoughness, i.worldNormal, i.viewDir, _WorldSpaceLightPos0.xyz, _LightColor0.rgb, indirectDiffuse, indirectSpecular);
                return c;
            }

            ENDCG
        }
    }
}

前節までのものからの主な変更点はコメントに記載した部分になります。

フラグメントシェーダを見ると、まずShadeSHPerPixel()を使って間接光の拡散反射成分を取得しています。

またその下ではunity_SpecCube0という名前のキューブマップをサンプリングしています。
これは間接光の鏡面反射成分に当たります。
前回の記事で説明したように鏡面反射成分はラフネスが大きいほどぼやけた感じになるので、
ラフネスの値をミップマップのレベルに影響させることでこれを実現しています。

レンダリング結果

それでは間接光も反映させた最終的なレンダリング結果を見ていきます。
色空間の設定はリニア、パイプラインは非SRPであることを確認してください。

f:id:halya_11:20200305214201p:plain:w600
最終的なレンダリング結果

間接光の影響が加わったことで質感が大分加わったことが確認できます。
もう少ししっかり現実的な材質設定を行うと以下のような見た目になります。

f:id:halya_11:20200305220249p:plain:w600
最終的なレンダリング結果（調整）

物理ベースの質感がちゃんと出ていることが確認できました。

連載一覧

light11.hatenadiary.com

light11.hatenadiary.com

参考

qiita.com

google.github.io

*1:Romain Guy and Mathias Agopian. 2018. Physically Based Rendering in Filament . https://google.github.io/filament/Filament.md.html#materialsystem/specularbrdf/normaldistributionfunction(speculard)