Unity Shader 简单实现一个PBR模型（不含IBL）

2020/6/6，答辩完成，虽然一大堆材料还没写，但是我已经很久没写博客，因为毕设弄的是和PBR有关的东西，所以回过头来看这篇文章，把错误的地方改正。
PBR究竟是个什么？我也弄不太明白，一开始我以为只是一个公式一个着色器罢了，但在看了大佬的专栏（中国特色技术美术）之后才知道这其实是一整套工序，但其实光公式部分就够让人头大的，特别是对于我这种数学弱鸡来说，推算部分看来看去看不明白。其他的我就先不管那么多了，现在在这里我只关注于Shader代码部分并且且不包含IBL（无IBL即无间接光照部分，本文章主要实现直接光照部分——2020/6/6添加）（这一部分还是需要点时间去理解，以后肯定会写的），推理部分这里贴一个大佬的文章topameng，已经写的十分详细了，也可以去LearnOpenGl看。

这是实现的效果，其实不用IBL的PBR和Blinn-Phong差别不大，贴图在这里可以找到Free PBR

这边的代码全是自己写的，可以说是几乎没用到Unity自带的东西，但其实大部分的东西Unity本身就已经实现了，我可以说是纯粹把LeranOpengGL的代码Copy了一遍,以下所有引用部分皆来自LearnOpenGl。
PBR最重要的理论基础是微表面和能量守恒（辐射度我是真的理不清,就直接忽略了），我们先来看看最终的公式：

BRDF，或者说双向反射分布函数，它接受入射（光）方向${\omega }_i$，出射（观察）方向${\omega }_o$，平面法线$n$以及一个用来表示微平面粗糙程度的参数$a$作为函数的输入参数。BRDF可以近似的求出每束光线对一个给定了材质属性的平面上最终反射出来的光线所作出的贡献程度。举例来说，如果一个平面拥有完全光滑的表面（比如镜面），那么对于所有的入射光线${\omega }_i$（除了一束以外）而言BRDF函数都会返回0.0 ，只有一束与出射光线${\omega }_o$拥有相同（被反射）角度的光线会得到1.0这个返回值。
BRDF基于我们之前所探讨过的微平面理论来近似的求得材质的反射与折射属性。对于一个BRDF，为了实现物理学上的可信度，它必须遵守能量守恒定律，也就是说反射光线的总和永远不能超过入射光线的总量。严格上来说，同样采用${\omega }_i$和${\omega }_o$作为输入参数的 Blinn-Phong光照模型也被认为是一个BRDF。然而由于Blinn-Phong模型并没有遵循能量守恒定律，因此它不被认为是基于物理的渲染。现在已经有很好几种BRDF都能近似的得出物体表面对于光的反应，但是几乎所有实时渲染管线使用的都是一种被称为Cook-Torrance BRDF模型。

$$L_o(p,{\omega }_o)={\int }_{\Omega }^{}(k_d\frac{c}{\pi }+k_s\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$
这就是Cook-Torrance反射率方程，即Cook-Torrance BRDF。

这个方程完整的描述了一个基于物理的渲染模型，它现在可以认为就是我们一般意义上理解的基于物理的渲染,也就是PBR。

现在我们需要把一步步拆开一步步转化成代码：
Cook-Torrance BRDF兼有漫反射和镜面反射两个部分（我们暂时忽略积分${\int }_{\Omega }^{}$以及$L_i(p,{\omega }_i)d{\omega }_i$部分，因为这部分貌似已经涉及到IBL了，$n\cdot {\omega }_i$不用忽略) （前面括号部分描述错误，在直接光照部分可以直接忽略积分项，因为渲染方程是对半球领域$\Omega$的积分—— 然而，当我们为一个表面上的特定的点$p$着色时，在其半球领域$\Omega$的所有可能的入射方向上，只有一个入射方向向量${\omega }_i$直接来自于该点光源。 假设我们在场景中只有一个光源，位于空间中的某一个点，因而对于p点的其他可能的入射光线方向上的辐射率为0——这是LearnOpenGL的原话，大致意思是在直接光照的情况下只需要计算一个方向的入射光，因为其他方向都是为零的（多个光源可以分开计算再将结果相加），所以积分就可以直接划去，而$L_i(p,{\omega }_i)$部分在直接光照的计算中是光源的强度衰减因子，在平行光时该值为1）：
$$f_r=k_d{f}_{lambert}+k_s{f}_{cook-torrance}$$
首先是漫反射部分，$f_{lambert}=\frac{c}{\pi }$。即表面颜色$\pi$,这是对光照进行标准化。
然后是镜面反射部分,$f_{cook-torrance}=\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)}$。${\omega }_i$和${\omega }_o$分别是光照方向、观察方向，$n$为法线。
D 正态分布函数（Normal Distribution Function）:估算在受到表面粗糙度的影响下，取向方向与中间向量一致的微平面的数量。这是用来估算微平面的主要函数。(也有人翻译成法线分布函数，其实法线分布更为准确——大佬们说的)
F 菲涅尔方程（Fresnel Rquation）：菲涅尔方程描述的是在不同的表面角下表面所反射的光线所占的比率。菲涅尔效果就是当视线与法线夹角越大折射效果减弱而反射效果增强，拿水来说我们近处的水透明见底可理解为光线发生了全折射而零反射，远处的水波光粼粼可理解为光线发生了零折射而全反射。
G 几何函数（Geometry Function）：描述了微平面自成阴影的属性。当一个平面相对比较粗糙的时候，平面表面上的微平面有可能挡住其他的微平面从而减少表面所反射的光线。

正态分布函数

正态分布函数D，或者说镜面分布，从统计学上近似的表示了与某些（中间）向量h取向一致的微平面的比率。举例来说，假设给定向量$h$，如果我们的微平面中有35%与向量$h$取向一致，则正态分布函数或者说NDF将会返回0.35。目前有很多种NDF都可以从统计学上来估算微平面的总体取向度，只要给定一些粗糙度的参数以及一个我们马上将会要用到的参数Trowbridge-Reitz GGX：

$$ND{F}_{GGXTR}(n,h,\alpha )=\frac{{\alpha }^2}{\pi ((n\cdot h{)}^2({\alpha }^2-1)+1{)}^2}$$

在这里$h$表示用来与平面上微平面做比较用的中间向量，而a表示表面粗糙度。

如果我们把$h$当成是不同粗糙度参数下，平面法向量和光线方向向量之间的中间向量的话，我们可以得到如下图示的效果：

当粗糙度很低（也就是说表面很光滑）的时候，与中间向量取向一致的微平面会高度集中在一个很小的半径范围内。由于这种集中性，NDF最终会生成一个非常明亮的斑点。但是当表面比较粗糙的时候，微平面的取向方向会更加的随机。你将会发现与h向量取向一致的微平面分布在一个大得多的半径范围内，但是同时较低的集中性也会让我们的最终效果显得更加灰暗。

Shader代码

float DistributionGGX(float3 N,float3 H,float roughness)
{
     float a2=roughness*roughness;
     a2=a2*a2;
     float NdotH=saturate(dot(N,H));
     float NdotH2=NdotH*NdotH;

     float denom=(NdotH2*(a2-1.0)+1.0);
     denom=UNITY_PI*denom*denom;
     return a2/denom;
}

菲涅尔函数

学过Shader的基本都知道的菲涅尔简化方程$F_{Schlick}(h,v,F_0)=F_0+(1-F_0)(1-(h\cdot v){)}^5$,平面对于法向入射的响应或者说基础反射率，可以在这里找到更多数据，不同于我以前用的$F_0$都是单一的float,这里是一个Vector

这里可以观察到的一个有趣的现象，所有电介质材质表面的基础反射率都不会高于0.17，这其实是例外而非普遍情况。导体材质表面的基础反射率起点更高一些并且（大多）在0.5和1.0之间变化。此外，对于导体或者金属表面而言基础反射率一般是带有色彩的，这也是为什么F0要用RGB三原色来表示的原因（法向入射的反射率可随波长不同而不同）。这种现象我们只能在金属表面观察的到。

金属表面这些和电介质表面相比所独有的特性引出了所谓的金属工作流的概念。也就是我们需要额外使用一个被称为金属度(Metalness)的参数来参与编写表面材质。金属度用来描述一个材质表面是金属还是非金属的。

通过预先计算电介质与导体的F0值，我们可以对两种类型的表面使用相同的Fresnel-Schlick近似，但是如果是金属表面的话就需要对基础反射率添加色彩。我们一般是按下面这个样子来实现的：

float3 F0=float3(0.04,0.04,0.04);
float metallic=tex2D(_Metallic,i.uv).r;
F0=lerp(F0,albedo,metallic);

我们为大多数电介质表面定义了一个近似的基础反射率。F0取最常见的电解质表面的平均值，这又是一个近似值。不过对于大多数电介质表面而言使用0.04作为基础反射率已经足够好了，而且可以在不需要输入额外表面参数的情况下得到物理可信的结果。然后，基于金属表面特性，我们要么使用电介质的基础反射率要么就使用F0来作为表面颜色。因为金属表面会吸收所有折射光线而没有漫反射，所以我们可以直接使用表面颜色纹理来作为它们的基础反射率。

菲涅尔简化方程代码

float3 FresnelSchlick(float cosTheta,float3 F0)
{
	return F0+(1.0-F0)*pow(1.0-cosTheta,5.0);
}

几何函数

几何函数从统计学上近似的求得了微平面间相互遮蔽的比率，这种相互遮蔽会损耗光线的能量。

与NDF类似，几何函数采用一个材料的粗糙度参数作为输入参数，粗糙度较高的表面其微平面间相互遮蔽的概率就越高。我们将要使用的几何函数是GGX与Schlick-Beckmann近似的结合体，因此又称为Schlick-GGX：
$$G_{SchlickGGX}(n,v,k)=\frac{n\cdot v}{(n\cdot v)(1-k)+k}$$
这里的$k$是$\alpha$基于几何函数是针对直接光照还是针对IBL光照的重映射(Remapping) :
$$k_{direct}=\frac{(\alpha +1{)}^2}{8}$$
$$k_{IBL}=\frac{{\alpha }^2}{2}$$

这里我们只针对直接光照不考虑IBL，所以只用第一个就好

为了有效的估算几何部分，需要将观察方向（几何遮蔽(Geometry Obstruction)）和光线方向向量（几何阴影(Geometry Shadowing)）都考虑进去。我们可以使用史密斯法(Smith’s method)来把两者都纳入其中：
$$G(n,v,l,k)=G_{sub}(n,v,k)G_{sub}(n,l,k)$$
使用史密斯法与Schlick-GGX作为$Gsub$可以得到如下所示不同粗糙度的视觉效果：

几何函数是一个值域为[0.0, 1.0]的乘数，其中白色或者说1.0表示没有微平面阴影，而黑色或者说0.0则表示微平面彻底被遮蔽。

代码

float3 GeometrySchlickGGX(float NdotV,float roughness)
{
	float r=roughness+1.0;
	float k=r*r/8.0;
	
	float denom=NdotV*(1.0-k)+k;
	return NdotV/denom;
}

float GeometrySmith(float3 N,float3 V,float3 L,float roughness)
{
	float NdotV=saturate(dot(N,V));
	float NdotL=saturate(dot(N,L));
	float ggx1=GeometrySchlickGGX(NdotV,roughness);
	float ggx2=GeometrySchlickGGX(NdotL,roughness);
	
	return ggx1*ggx2;
}

以下是完整代码

Shader "MyShader/PBRTexture"
{
    Properties
    {
        _MainTex ("Albedo", 2D) = "white" {}
        [NoScaleOffset]_BumpTex("Normal Map",2D)="bump"{}
        [NoScaleOffset]_Metallic("Metallic",2D)="metallic"{}
        _Color("Color",Color)=(1,1,1,1)
        _BumpSacle("Bump Sacle",Range(-1,1))=1
        _Roughness("Roughness",Range(0,1))=0.1
        _AO("AO",Range(0,1))=0.1
    }
    SubShader
    {
        Tags { "RenderType"="Opaque" "LightMode"="ForwardBase"}
        LOD 100

        CGINCLUDE
        #include "UnityCG.cginc"
        #include "Lighting.cginc"

        float3 FresnelSchlick(float cosTheta,float3 F0)
        {
            return F0+(1.0-F0)*pow(1.0-cosTheta,5.0);
        }

        float DistributionGGX(float3 N,float3 H,float roughness)
        {
            float a2=roughness*roughness;
            a2=a2*a2;
            float NdotH=saturate(dot(N,H));
            float NdotH2=NdotH*NdotH;

            float denom=(NdotH2*(a2-1.0)+1.0);
            denom=UNITY_PI*denom*denom;
            return a2/denom;
        }

        float3 GeometrySchlickGGX(float NdotV,float roughness)
        {
            float r=roughness+1.0;
            float k=r*r/8.0;

            float denom=NdotV*(1.0-k)+k;
            return NdotV/denom;
        }

        float GeometrySmith(float3 N,float3 V,float3 L,float roughness)
        {
            float NdotV=saturate(dot(N,V));
            float NdotL=saturate(dot(N,L));
            float ggx1=GeometrySchlickGGX(NdotV,roughness);
            float ggx2=GeometrySchlickGGX(NdotL,roughness);

            return ggx1*ggx2;
        }
        ENDCG

        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            struct appdata
            {
                float4 vertex : POSITION;
                float2 uv : TEXCOORD0;
                float3 normal:NORMAL;
                float4 tangent:TANGENT;
            };

            struct v2f
            {
                float2 uv : TEXCOORD0;
                float4 vertex : SV_POSITION;
                float3 lightDir:TEXCOORD1;
                float3 viewDir:TEXCOORD2;
            };

            sampler2D _MainTex;
            float4 _MainTex_ST;
            sampler2D _BumpTex;
            float _BumpSacle;
            float _Roughness;
            sampler2D _Metallic;
            float _AO;
            fixed4 _Color;

            v2f vert (appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.uv = TRANSFORM_TEX(v.uv, _MainTex);
				//用法线要用的，基本操作了
                TANGENT_SPACE_ROTATION;
                o.lightDir=mul(rotation,ObjSpaceLightDir(v.vertex).xyz);
                o.viewDir=mul(rotation,ObjSpaceViewDir(v.vertex).xyz);
                return o;
            }

            fixed4 frag (v2f i) : SV_Target
            {

                float4 packNormal=tex2D(_BumpTex,i.uv);
                float3 normal=UnpackNormal(packNormal);
                normal.xy*=_BumpSacle;
                normal.z=sqrt(1.0- saturate(dot(normal.xy,normal.xy)));
                float3 lightDir=normalize(i.lightDir);
                float3 viewDir=normalize(i.viewDir);
                float3 halfDir=normalize(lightDir+viewDir);
				
				//把色彩转到线性空间，我不太确定Unity的光照颜色是线性的还是Gamma的......
                fixed4 lightColor=pow(_LightColor0,2.2);
                fixed4 color=pow(_Color,2.2);

                //fixed4 albedo = pow(tex2D(_MainTex, i.uv),2.2)*_LightColor0*_Color;
                fixed4 albedo = pow(tex2D(_MainTex, i.uv),2.2)*lightColor*color;
                float3 F0=float3(0.04,0.04,0.04);

                float metallic=tex2D(_Metallic,i.uv).r;
                F0=lerp(F0,albedo,metallic);
                float3 Lo=float3(0,0,0);

                float NDF=DistributionGGX(normal,halfDir,_Roughness);
                float G=GeometrySmith(normal,viewDir,lightDir,_Roughness);
                float3 F=FresnelSchlick(clamp(dot(halfDir,viewDir),0,1),F0);

                float3 nom=NDF*G*F;
                float3 denom=4*max(dot(normal,viewDir),0)*saturate(dot(normal,lightDir));
                float3 specular=nom/max(denom,0.001);//max避免denom为零

                float3 Ks=F;
                float3 Kd=1-Ks;
                Kd*=1.0-metallic;
                float NdotL=max(dot(normal,lightDir),0.0);
                Lo+=(Kd*albedo/UNITY_PI+specular)*NdotL;//
                float3 ambient=0.03*albedo*_AO;
                float4 finalColor=float4(ambient+Lo,1.0);
				
				//转回Gamma空间
                finalColor=finalColor/(finalColor+0.1);
                finalColor=pow(finalColor,1.0/2.2);
                finalColor.a=1;
                
                return finalColor;
            }

            ENDCG
        }
    }
}

最终的结果Lo，或者说是出射光线的辐射率，实际上是反射率方程的在半球领域Ω的积分的结果。但是我们实际上不需要去求积，因为对于所有可能的入射光线方向我们知道只有4个方向（我们这里只用了一个灯光）的入射光线会影响片段(像素)的着色。因为这样，我们可以直接循环N次计算这些入射光线的方向(N也就是场景中光源的数目)。
比较重要的是我们没有把kS乘进去我们的反射率方程中，这是因为我们已经在specualr BRDF中乘了菲涅尔系数F了，因为kS等于F，因此我们不需要再乘一次。

即
$$L_o(p,{\omega }_o)={\int }_{\Omega }^{}(k_d\frac{c}{\pi }+k_s\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$
并非完全数学意义上的正确，因为F（菲涅尔）已经等于$k_s$,意味着镜面反射部分已经包含$k_s$，说以最终可以改成
$$L_o(p,{\omega }_o)={\int }_{\Omega }^{}(k_d\frac{c}{\pi }+\frac{DFG}{4({\omega }_o\cdot n)({\omega }_i\cdot n)})L_i(p,{\omega }_i)n\cdot {\omega }_{i}d{\omega }_i$$

这个只是个用来学习的Shader，最后写成这个样子也是受我自己的能力所限，而且我自己也有许多的东西没用弄明白，如果有错误非常欢迎指正。感谢你的阅读。