PBR:双向反射分布函数(BRDF)介绍与Cook-Torrance模型的实现

PBR:双向反射分布函数(BRDF)介绍与Cook-Torrance模型的实现

 

BRDF简介

再介绍BRDF之前我们要引入渲染方程这个东西:

其中L表示辐射率,其公式为:

它表示了一个拥有辐射强度Φ的光源在单位面积A,单位立体角ω上的辐射出的总能量。

辐射率是辐射度量学上表示一个区域平面上光线总量的物理量,它受到入射光线与平面法线间的夹角θ的余弦值cosθ的影响:当直接辐射到平面上和法线夹角越大时,光线就越弱。

我们继续看渲染方程,L_{o}w_{0}接收的光量由两部分组成,其中L_{e}是自发光在w_{0}方向的辐射光量,而’+‘号右边的公式则是我们的反射率方程。如果物体不会自发光,则:

该函数表示了表面A在接收了半球上所有的入射光量L_{i}乘上f_{r},然后它们的累加在反射方向w_{0}的光量为L_{o}。用一个最简单的例子,假设场景只有一个光源,且忽视间接照明:

我们有:

L_{o}(p,w_{o})=f_{r}(p,w_{i},w_{o})L_{i}(p,w_{I})n\cdot w_{i}

这里的f_{r}就是我们的BRDF函数。许多地方会看到BTDF(双向透射分布函数),BSDF(双向散射分布函数)或者SVBRDF(空间变化的BRDF)等。BTDF表明若有些材料如玻璃等需要考虑镜面透射则可以使用。BSDF则是BRDF和BTDF的结合使用。基本BRDF或者SVBRDF已经能满足我们的需求了。

我们现在来看一下BRDF。其中ωo和ωi分别是表面法线(N)和视线方向(V)与光方向之间的夹角。因此就有了“双向反射分布函数”的美称。BRDF可以近似的求出每束光线对一个给定了材质属性的平面上最终反射出来的光线所作出的贡献程度。举例来说,如果一个平面拥有完全光滑的表面(比如镜面),那么对于所有的入射光线ωi而言BRDF函数都会返回0.0,只有一束与出射光线ωo拥有相同(被反射)角度的光线会得到1.0这个返回值。到目前为止已经提出了很多的BRDF模型包括Phong,Blinn-Phong,Lafortune,Torrance-Sparrow,Cook-Torrance,Ward各向异性,Oren-Nayar等等。其中Phong模型是我们最常见的也是最简单,但是由于其是经验型的,不符合能量守恒,我们本文就介绍Cook-Torrance模型。

 

Cook-Torrance模型

Cook-Torrance BRDF兼有漫反射和镜面反射两个部分:

这里的k_{d}是早先提到过的入射光线中被折射部分的能量所占的比率,而k_{s}被反射部分的比率。BRDF的左侧表示的是漫反射部分,这里用f_{Lambert}来表示。这是Lambertian漫反射。右侧是我们的镜面反射部分。

 

Lambertian反射

Lambertian反射是用来模拟完美的漫反射模型。它以兰伯特余弦定律的名称而闻名。表面接收的光量与表面法线N和光方向L之间的角度成正比。

表明了:

漫反射颜色和cos(\theta )成正比(∝为正比)。这还需要知道反照率albedo:

即反射光量和入射光量的比。我们用符号\large \rho _{d}表示albedo。根据lambert模型我们就可以计算出漫反射的光量为:

\large Diffuse Surface Color=\rho _{d}*Incident Light Energy*N\cdot L

但实际上并不是这样的,我们实际的lanbert公式为:

为什么要除以\large \pi呢?根据lambert漫反射定义:

我们可以用积分计算出反射的光量:

即我们的反射总光量不能大于入射光量。求解出积分得:

但是我们定义的反照率属于[0,1]。所以看上式会有问题,我们只要将\large \rho _{d}除以\large \pi,便可以使式子正确,即:

所以我们的漫反射为:

 

BRDF的镜面反射部分

BRDF的镜面部分有如下形式:

Cook-Torrance BRDF的镜面反射部分包含三个函数,此外分母部分还有一个标准化因子 。字母D,F与G分别代表着一种类型的函数,各个函数分别用来近似的计算出表面反射特性的一个特定部分。三个函数分别为法线分布函数(Normal Distribution Function),菲涅尔方程(Fresnel Rquation)和几何函数(Geometry Function):

  • 法线分布函数:估算在受到表面粗糙度的影响下,取向方向与中间向量一致的微平面的数量。这是用来估算微平面的主要函数。
  • 几何函数:描述了微平面自成阴影的属性。当一个平面相对比较粗糙的时候,平面表面上的微平面有可能挡住其他的微平面从而减少表面所反射的光线。
  • 菲涅尔方程:菲涅尔方程描述的是在不同的表面角下表面所反射的光线所占的比率。

以上的每一种函数都是用来估算相应的物理参数的,而且你会发现用来实现相应物理机制的每种函数都有不止一种形式。它们有的非常真实,有的则性能高效。这里我们其中D使用Trowbridge-Reitz GGX,F使用Fresnel-Schlick近似(Fresnel-Schlick Approximation),而G使用Smith’s Schlick-GGX。

 

法线分布函数

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

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

如果我们把h当成是不同粗糙度参数下,平面法向量和光线方向向量之间的中间向量的话(实际运算H = normalize(V(观察方向) + L(光线方向))),我们可以得到如下图示的效果:

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

使用GLSL代码编写的Trowbridge-Reitz GGX法线分布函数是下面这个样子的:

float D_GGX_TR(vec3 N, vec3 H, float a)
{
    float a2     = a*a;
    float NdotH  = max(dot(N, H), 0.0);
    float NdotH2 = NdotH*NdotH;

    float nom    = a2;
    float denom  = (NdotH2 * (a2 - 1.0) + 1.0);
    denom        = PI * denom * denom;

    return nom / denom;
}

 

几何函数

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

与NDF类似,几何函数采用一个材料的粗糙度参数作为输入参数,粗糙度较高的表面其微平面间相互遮蔽的概率就越高。我们将要使用的几何函数是GGX与Schlick-Beckmann近似的结合体,因此又称为Schlick-GGX:

这里的k根据使用直接光照还是IBL(基于图像的照明)会有不同:

值得注意的是根据引擎把粗糙度转化为α的方式不同,得到α的值也有可能不同。为了有效的估算几何部分,需要将观察方向(几何遮蔽(Geometry Obstruction))和光线方向向量(几何阴影(Geometry Shadowing))都考虑进去。我们可以使用史密斯法(Smith’s method)来把两者都纳入其中:

使用史密斯法与Schlick-GGX作为Gsub可以得到如下所示不同粗糙度的视觉效果:

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

使用GLSL编写的几何函数代码如下:

float GeometrySchlickGGX(float NdotV, float k)
{
    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;

    return nom / denom;
}

float GeometrySmith(vec3 N, vec3 V, vec3 L, float k)
{
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx1 = GeometrySchlickGGX(NdotV, k);
    float ggx2 = GeometrySchlickGGX(NdotL, k);

    return ggx1 * ggx2;
}

 

菲涅尔方程

菲涅尔方程描述的是被反射的光线对比光线被折射的部分所占的比率,这个比率会随着我们观察的角度不同而不同。当光线碰撞到一个表面的时候,菲涅尔方程会根据观察角度告诉我们被反射的光线所占的百分比。利用这个反射比率和能量守恒原则,我们可以直接得出光线被折射的部分以及光线剩余的能量。

当垂直观察的时候,任何物体或者材质表面都有一个基础反射率(Base Reflectivity),但是如果以一定的角度往平面上看的时候所有反光都会变得明显起来。你可以自己尝试一下,用垂直的视角观察你自己的木制/金属桌面,此时一定只有最基本的反射性。但是如果你从近乎90度(是指和法线的夹角)的角度观察的话反光就会变得明显的多。如果从理想的90度视角观察,所有的平面理论上来说都能完全的反射光线。这种现象因菲涅尔而闻名,并体现在了菲涅尔方程之中。如下图:

我们发现近处的水面多呈现折射特性(观察方向和水面法线夹角小),远处呈现反射特性(观察方向和水面法线夹角大)。

在菲涅耳方程中,光由两个垂直波组成,我们称其为平行偏振光和垂直偏振光。我们需要使用两个不同的方程式(对于每种类型的波而言)计算这两个波的反射光比率,并对结果求平均值以找到解。两个菲涅耳方程为:

通过对这两个值求平均值,便可以获得我们反射比率:

项η1,η2是两种介质的折射率。项cosθ1和cosθ2分别是入射角和折射角。如前所述,由于能量守恒,折射光的比率可以简单地计算为:

但本文使用Fresnel-Schlick近似法求得菲涅尔方程近似解:

F0表示平面的基础反射率,它是利用所谓折射指数(Indices of Refraction)或者说IOR计算得出的。然后正如你可以从球体表面看到的那样,我们越是朝球面掠角的方向上看(此时视线和表面法线的夹角接近90度)菲涅尔现象就越明显,反光就越强:

菲涅尔方程还存在一些细微的问题。其中一个问题是Fresnel-Schlick近似仅仅对电介质(非金属)表面有定义。对于导体(Conductor)表面)金属),使用它们的折射指数计算基础折射率并不能得出正确的结果,这样我们就需要使用一种不同的菲涅尔方程来对导体表面进行计算。由于这样很不方便,所以我们预先计算出平面对于法向入射(F0)的反应(处于0度角,好像直接看向表面一样)然后基于相应观察角的Fresnel-Schlick近似对这个值进行插值,用这种方法来进行进一步的估算。这样我们就能对金属和非金属材质使用同一个公式了。

平面对于法向入射的响应或者说基础反射率可以在一些大型数据库中找到,比如:

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

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

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

vec3 F0 = vec3(0.04);
F0      = mix(F0, surfaceColor.rgb, metalness);

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

Fresnel Schlick近似可以用代码表示为:

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

 

Cook-Torrance BRDF中所有元素都已经介绍完毕,我们现在可以将基于物理的BRDF纳入到最终的反射率方程当中去了:

 

完善模型

知道如何计算反射率方程后,我们需要完善我们的BRDF模型,我们列出我们需要的属性代码(片段着色器):

#version 330 core
out vec4 FragColor;
in vec2 TexCoords;
in vec3 WorldPos;
in vec3 Normal;

uniform vec3 camPos;

uniform vec3  albedo;
uniform float metallic;
uniform float roughness;
uniform float ao;

我们这里看到了ao,这是我们的环境光遮蔽参数,因为我们没有精确的计算场景的间接照明,所以我们可以使用ao来粗略模拟环境光遮蔽。根据上面讲述我们给出片段着色器代码:

#version 410 core

out vec4 FragColor;
in vec2 TexCoords;
in vec3 WorldPos;
in vec3 Normal;
// 材质参赛
uniform sampler2D albedoMap;
uniform sampler2D normalMap;
uniform sampler2D metallicMap;
uniform sampler2D roughnessMap;
uniform sampler2D aoMap;

uniform float uRoughness;
uniform float uMetallic;

// 灯光
uniform vec3 lightPositions[4];
uniform vec3 lightColors[4];

uniform vec3 camPos;

const float PI = 3.14159265359;
// ----------------------------------------------------------------------------
float DistributionGGX(vec3 N, vec3 H, float roughness)
{
    float a = roughness*roughness;//采用平方更加自然
    float a2 = a*a;
    float NdotH = max(dot(N, H), 0.0);
    float NdotH2 = NdotH*NdotH;
    
    float nom   = a2;
    float denom = (NdotH2 * (a2 - 1.0) + 1.0);
    denom = PI * denom * denom;
    
    return nom / (denom); // 避免分母为0
}
// ----------------------------------------------------------------------------
float GeometrySchlickGGX(float NdotV, float roughness)
{
    float r = (roughness + 1.0);
    float k = (r*r) / 8.0;
    
    float nom   = NdotV;
    float denom = NdotV * (1.0 - k) + k;
    
    return nom / denom;
}
// ----------------------------------------------------------------------------
float GeometrySmith(vec3 N, vec3 V, vec3 L, float roughness)
{
    float NdotV = max(dot(N, V), 0.0);
    float NdotL = max(dot(N, L), 0.0);
    float ggx2 = GeometrySchlickGGX(NdotV, roughness);
    float ggx1 = GeometrySchlickGGX(NdotL, roughness);
    
    return ggx1 * ggx2;
}
// ----------------------------------------------------------------------------
vec3 fresnelSchlick(float cosTheta, vec3 F0)
{
    return F0 + (1.0 - F0) * pow(1.0 - cosTheta, 5.0);
}
// ----------------------------------------------------------------------------

vec3 getNormalFromMap(){
    vec3 tangentNormal=texture(normalMap,TexCoords).xyz*2.0-1.0;
    
    vec3 Q1=dFdx(WorldPos);
    vec3 Q2=dFdy(WorldPos);
    vec2 st1=dFdx(TexCoords);
    vec2 st2=dFdy(TexCoords);
    
    vec3 N=normalize(Normal);
    vec3 T=normalize(Q1*st2.t-Q2*st1.t);
    vec3 B=-normalize(cross(N,T));
    mat3 TBN=mat3(T,B,N);
    return normalize(TBN*tangentNormal);
}

void main()
{
//    vec3 albedo=pow(texture(albedoMap,TexCoords).rgb,vec3(2.2));
//    float roughness=pow(texture(roughnessMap,TexCoords).r,2.2);
//    float metallic=pow(texture(metallicMap,TexCoords).r,2.2);
    vec3 albedo=vec3(1.0,0.71,0.29);
    float metallic=uMetallic;
    float roughness=uRoughness;
    float ao=texture(aoMap,TexCoords).r;
    
    vec3 N = getNormalFromMap();
    vec3 V = normalize(camPos - WorldPos);
    
    vec3 F0 = vec3(0.04);
    F0 = mix(F0, albedo, metallic);
    
    // 反射率方程
    vec3 Lo = vec3(0.0);
    for(int i = 0; i < 4; ++i)
    {
        // 计算每个光源的辐射
        vec3 L = normalize(lightPositions[i] - WorldPos);
        vec3 H = normalize(V + L);
        float distance = length(lightPositions[i] - WorldPos);
        float attenuation = 1.0 / (distance * distance);
        vec3 radiance = lightColors[i] * attenuation;
        
        // Cook-Torrance BRDF
        float NDF = DistributionGGX(N, H, roughness);
        float G   = GeometrySmith(N, V, L, roughness);
        vec3 F    = fresnelSchlick(clamp(dot(H, V), 0.0, 1.0), F0);
        
        vec3 nominator    = NDF * G * F;
        float denominator = 4 * max(dot(N, V), 0.0) * max(dot(N, L), 0.0);
        vec3 specular = nominator / max(denominator,0.001); // prevent divide by zero for NdotV=0.0 or NdotL=0.0
        
        // kS 是菲涅尔的值
        vec3 kS = F;

        vec3 kD = vec3(1.0) - kS;

        kD *= (1.0 - metallic);
        
        // 计算NdotL
        float NdotL = max(dot(N, L), 0.0);
        
        // 多光源累加
        Lo += (kD * albedo / PI + specular) * radiance * NdotL;
    }
    

    vec3 ambient = vec3(0.03,0.03,0.03) * albedo;
    
    vec3 color = ambient + Lo;
    
    // HDR Reinhard色调映射
    color = color / (color + vec3(1.0));
    // gamma矫正
    color = pow(color, vec3(1.0/2.2));
    
    FragColor = vec4(color, 1.0);
}

我们看到,我们还可以使用贴图形式给出我们的参数。该代码和之前讲述有点不同的是,我们执行了HDR和gamma矫正。由于直到现在,我们假设的所有计算都在线性的颜色空间中进行的,因此我们需要在着色器最后做gamma矫正。 在线性空间中计算光照是非常重要的,因为PBR要求所有输入都是线性的,如果不是这样,我们就会得到不正常的光照。另外,我们希望所有光照的输入都尽可能的接近他们在物理上的取值,这样他们的反射率或者说颜色值就会在色谱上有比较大的变化空间。Lo作为结果可能会变大得很快(超过1),但是因为默认的LDR输入而取值被截断。所以在伽马矫正之前我们采用色调映射使Lo从HDR的值映射为LDR的值。

完成结果如下图所示:

其中我们可以设定金的F0,图中从左到右粗糙度增加,从下到上金属性增加。

理论上来说,一个表面的金属度应该是二元的:要么是金属要么不是金属,不能两者皆是。但是,大多数的渲染管线都允许在0.0至1.0之间线性的调配金属度。这主要是由于材质纹理精度不足以描述一个拥有诸如细沙/沙状粒子/刮痕的金属表面。通过对这些小的类非金属粒子/刮痕调整金属度值,我们可以获得非常好看的视觉效果。

值得注意的是,Fresnel-Schlick近似只适用于非金属,我们虽然可以设定金属性我们可以将其定义为(0,1),去达到一些复合表面效果。但实际上电介质F0不会高于0.17,而金属的F0大多在0.5~1.0之间变化。而且如果一个材质是金属的话,是没有透射值的,即kD应赋值0。像上图的金F0参数,只有第一行才是正确的金材质。

  • 15
    点赞
  • 44
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值