精华内容
下载资源
问答
  • BRDF
    2019-10-21 20:43:56

    BRDF分为漫反射和高光反射部分

    漫反射部分

    unity采用的是Disney漫反射模型

    half DisneyDiffuse(half NdotV, half NdotL, half LdotH, half roughness)

    {

        half nlPow5 = Pow5 (1-NdotL);

        half nvPow5 = Pow5 (1-NdotV);

        half Fd90 = 0.5 + 2 * LdotH * LdotH * roughness;

        half disneyDiffuse = (1 + (Fd90-1) * nlPow5) * (1 + (Fd90-1) * nvPow5);

     

        half diffuseTerm = disneyDiffuse * NdotL;

        return diffuseTerm;

    }

    高光反射部分

    spec term比较多,分GDF三大部分

    //G

    //------------------------------------------------------------------------------------------------------------------------

    //Geometric Shadowing Functions

     

    float ImplicitGeometricShadowingFunction (float NdotL, float NdotV){

        float Gs = (NdotL*NdotV);

        return Gs;

    }

     

    float AshikhminShirleyGeometricShadowingFunction (float NdotL, float NdotV, float LdotH){

        float Gs = NdotL*NdotV/(LdotH*max(NdotL,NdotV));

        return (Gs);

    }

     

    float AshikhminPremozeGeometricShadowingFunction (float NdotL, float NdotV){

        float Gs = NdotL*NdotV/(NdotL+NdotV - NdotL*NdotV);

        return (Gs);

    }

     

    float DuerGeometricShadowingFunction (float3 lightDirection,float3 viewDirection, float3 normalDirection,float NdotL, float NdotV){

        float3 LpV = lightDirection + viewDirection;

        float Gs = dot(LpV,LpV) * pow(dot(LpV,normalDirection),-4);

        return (Gs);

    }

     

    float NeumannGeometricShadowingFunction (float NdotL, float NdotV){

        float Gs = (NdotL*NdotV)/max(NdotL, NdotV);

        return (Gs);

    }

     

    float KelemenGeometricShadowingFunction (float NdotL, float NdotV, float LdotH, float VdotH){

        // float Gs = (NdotL*NdotV)/ (LdotH * LdotH); //this

        float Gs = (NdotL*NdotV)/(VdotH * VdotH); //or this?

        return (Gs);

    }

     

    float ModifiedKelemenGeometricShadowingFunction (float NdotV, float NdotL, float roughness)

    {

        float c = 0.797884560802865; // c = sqrt(2 / Pi)

        float k = roughness * roughness * c;

        float gH = NdotV * k +(1-k);

        return (gH * gH * NdotL);

    }

     

    float CookTorrenceGeometricShadowingFunction (float NdotL, float NdotV, float VdotH, float NdotH){

        float Gs = min(1.0, min(2*NdotH*NdotV / VdotH, 2*NdotH*NdotL / VdotH));

        return (Gs);

    }

     

    float WardGeometricShadowingFunction (float NdotL, float NdotV, float VdotH, float NdotH){

        float Gs = pow( NdotL * NdotV, 0.5);

        return (Gs);

    }

     

    float KurtGeometricShadowingFunction (float NdotL, float NdotV, float VdotH, float alpha){

        float Gs = (VdotH*pow(NdotL*NdotV, alpha))/ NdotL * NdotV;

        return (Gs);

    }

     

    //F

    //------------------------------------------------------------------------------------------------------------------------

    //schlick functions

    float SchlickFresnel(float i){

        float x = clamp(1.0-i, 0.0, 1.0);

        float x2 = x*x;

        return x2*x2*x;

    }

    float3 FresnelLerp (float3 x, float3 y, float d)

    {

        float t = SchlickFresnel(d);

        return lerp (x, y, t);

    }

     

    float3 SchlickFresnelFunction(float3 SpecularColor,float LdotH){

        return SpecularColor + (1 - SpecularColor)* SchlickFresnel(LdotH);

    }

     

    float SchlickIORFresnelFunction(float ior,float LdotH){

        float f0 = pow((ior-1)/(ior+1),2);

        return f0 + (1 - f0) * SchlickFresnel(LdotH);

    }

    float SphericalGaussianFresnelFunction(float LdotH,float SpecularColor)

    {

        float power = ((-5.55473 * LdotH) - 6.98316) * LdotH;

        return SpecularColor + (1 - SpecularColor) * pow(2,power);

    }

     

    //D

    //------------------------------------------------------------------------------------------------------------------------

    //Normal Distribution Functions

     

    float BlinnPhongNormalDistribution(float NdotH, float specularpower, float speculargloss){

        float Distribution = pow(NdotH,speculargloss) * specularpower;

        Distribution *= (2+specularpower) / (2*3.1415926535);

        return Distribution;

    }

    float PhongNormalDistribution(float RdotV, float specularpower, float speculargloss){

        float Distribution = pow(RdotV,speculargloss) * specularpower;

        Distribution *= (2+specularpower) / (2*3.1415926535);

        return Distribution;

    }

     

    float BeckmannNormalDistribution(float roughness, float NdotH)

    {

        float roughnessSqr = roughness*roughness;

        float NdotHSqr = NdotH*NdotH;

        return max(0.000001,(1.0 / (3.1415926535*roughnessSqr*NdotHSqr*NdotHSqr))* exp((NdotHSqr-1)/(roughnessSqr*NdotHSqr)));

    }

     

    float GaussianNormalDistribution(float roughness, float NdotH)

    {

        float roughnessSqr = roughness*roughness;

        float thetaH = acos(NdotH);

        return exp(-thetaH*thetaH/roughnessSqr);

    }

     

    float GGXNormalDistribution(float roughness, float NdotH)

    {

        float roughnessSqr = roughness*roughness;

        float NdotHSqr = NdotH*NdotH;

        float TanNdotHSqr = (1-NdotHSqr)/NdotHSqr;

        return (1.0/3.1415926535) * sqr(roughness/(NdotHSqr * (roughnessSqr + TanNdotHSqr)));

        // float denom = NdotHSqr * (roughnessSqr-1)

     

    }

     

    float TrowbridgeReitzNormalDistribution(float NdotH, float roughness){

        float roughnessSqr = roughness*roughness;

        float Distribution = NdotH*NdotH * (roughnessSqr-1.0) + 1.0;

        return roughnessSqr / (3.1415926535 * Distribution*Distribution);

    }

     

    float TrowbridgeReitzAnisotropicNormalDistribution(float anisotropic, float NdotH, float HdotX, float HdotY){

        float aspect = sqrt(1.0h-anisotropic * 0.9h);

        float X = max(.001, sqr(1.0-_Glossiness)/aspect) * 5;

        float Y = max(.001, sqr(1.0-_Glossiness)*aspect) * 5;

        return 1.0 / (3.1415926535 * X*Y * sqr(sqr(HdotX/X) + sqr(HdotY/Y) + NdotH*NdotH));

    }

     

    float WardAnisotropicNormalDistribution(float anisotropic, float NdotL, float NdotV, float NdotH, float HdotX, float HdotY){

        float aspect = sqrt(1.0h-anisotropic * 0.9h);

        float X = max(.001, sqr(1.0-_Glossiness)/aspect) * 5;

        float Y = max(.001, sqr(1.0-_Glossiness)*aspect) * 5;

        float exponent = -(sqr(HdotX/X) + sqr(HdotY/Y)) / sqr(NdotH);

        float Distribution = 1.0 / ( 3.14159265 * X * Y * sqrt(NdotL * NdotV));

        Distribution *= exp(exponent);

        return Distribution;

    }

     

    更多相关内容
  • 软件介绍: brdf explorer是一迪士尼开发的brdf的浏览工具,功能非常强大,这个是V1.0版本,绿色软件不用安装。能够以各种视图模式下查看。
  • 参见 Example.m 文件
  • pbrt-重要性-采样 PBRT 扩展对 brdf 数据使用重要性采样算法用于数据压缩(brdf 功能):Douglas Peucker
  • 研究了沙漠地区地表起伏程度及太阳天顶角对BRDF 的影响,结果表明,在太阳几何一定的情况下,不同起伏地形的BRDF 随观测方位角的变化呈不对称分布,且近垂直主平面方向的BRDF 比近主平面方向整体偏大;对于同一地形...
  • BRDF模型介绍

    2017-11-22 22:09:21
    BRDF模型介绍得非常详细,photemetric stereo多光技术的非朗伯反射方法大多基于BRDF模型,有了一个清晰的认识
  • matlab集成c代码 /系统各模块说明*****/ 一、系统初始化 1、材质台就位 2、相机、光源、材质台控制测试 二、预处理 1、工业相机标定及校正预处理 ...8)逐角度计算图像平均亮度/逐像素亮度,保存BRDF材质数据
  • 提出了一种快速测量双向反射分布函数(BRDF)的子孔径扫描傅里叶变换系统。子孔径扫描傅里叶变换方法,使用单个傅里叶变换透镜进行测量,具有扩展单次测量视场、工作距离长、成本低以及系统简单等优点,可以在大空间频率...
  • 在此将褶皱看作一种材质,提出基于宏观光学散射截面测量的双向反射分布函数(BRDF)生成方法,求得褶皱材质的BRDF数据,进一步利用误差逆传播(BP)神经网络建立了褶皱材质的BRDF模型,代替了复杂的褶皱建模过程,大大简化了...
  • 由于双向反射分布函数( BRDF)经验模型与半经验模型对材质散射特性描述时存在局限性,导致其拟合结果与实测数据的误差较大。针对此问题,基于深度神经网络(DNN)构建了一种适用于具有不同散射特性空间目标材质的BRDF...
  • 将均匀度指数引入林木榨间分布格对双向反射分布函数(BRDF)影响的分析中,利川光学遥感模型定最分析林木空间分布格局对BRDF的影响、结果表明:对于红光波段,在大于树冠尺度上均匀度口以以有效地描述林木空间分布局...
  • 开发这些程序是为了利用 Eldim 开发的 EZ-Contrast 仪器的亮度或 BRDF 数据。 它们适用于使用 ccd 传感器进行的任何测量。 开发的库包括: - 最大值研究- ccd 截面的极坐标天顶图- ccd 截面的极坐标图- 根据天顶角...
  • 对目标表面双向反射分布函数(BRDF)进行研究,在各向异性Gaussian模型的基础上引入波长因子,提出了改进的BRDF各向异性Gaussian计算模型。运用该模型对卫星太阳能板及包裹多层隔热材料进行了BRDF计算,结果表明该模型...
  • brdf浏览工具

    2018-06-16 11:13:58
    非常好的brdf的浏览工具,迪士尼开发,超强大,非常好用
  • disney brdf viwer

    2017-06-26 14:33:59
    disney brdf viwer 可视化
  • Anisotropic Phong BRDF Model

    2019-02-20 09:10:20
    Physically Based Shading技术中的各向异性BRDF模型技术原理
  • 文章目录光照、照明(Illumination)预备知识球面坐标(Spherical Coordinate)立体角(Solid ...BRDF的定义与性质BRDF 的定义可逆性(Reciprocity)能量守恒性质渲染方程(Rendering Equation)BRDF模型经验模型(Em


    双向反射分布函数(BRDF:Bidirectional Reflectance Distribution Function),描述的是物体表面对光的反射性质。

    光照、照明(Illumination)

    光照可分为局部和全局两种
    局部光照(local illumination) 和物体直接被光源照射的情况相关
    全局光照(global illumination) 和物体被光源以外的其它地方来的光线照射的情况相关,包括其它物体反射和投射过来的光线。

    今天的主题:一个关于光线如何被表面反射的物理描述,也就是 BRDF。

    预备知识

    球面坐标(Spherical Coordinate)
    立体角(Solid Angle)
    投影面积(Foreshortened Area)
    光能(Radiant Energy)
    光通量(Radiant Flux)
    辉度(Irradiance)
    发光强度(Intensity)
    光亮度(Radiance)

    球面坐标(Spherical Coordinate)

    由于光线主要是通过方向来表达,通常用球面坐标表达它们比用笛卡尔坐标更方便。
    如图所示,球面坐标中的向量用三个元素来指定:
    r r r 表示向量的长度;
    θ \theta θ 表示向量和 z z z 轴的夹角;
    ϕ \phi ϕ 表示向量在 x − y x-y xy 平面上的投影和 x x x 轴的逆时针夹角。
    在这里插入图片描述
    球面坐标与三维笛卡尔坐标(Cartesian Coordinates)之间的对应关系:

    C a r ( x , y , z ) = S p h ( r , θ , ϕ ) Car(x, y, z) = Sph(r, \theta, \phi) Car(x,y,z)=Sph(r,θ,ϕ)

    它们的转换公式为:
    r = x 2 + y 2 + z 2 θ = arccos ⁡ ( z / r ) ϕ = arctan ⁡ ( y / x ) \begin{aligned} r &= \sqrt{x^2 + y^2 + z^2} \\ \theta &= \arccos(z/r) \\ \phi &= \arctan(y/x) \end{aligned} rθϕ=x2+y2+z2 =arccos(z/r)=arctan(y/x)

    z = r cos ⁡ ( θ ) y = r sin ⁡ ( θ ) s i n ( ϕ ) x = r sin ⁡ ( θ ) c o s ( ϕ ) \begin{aligned} z &= r \cos (\theta) \\ y &= r \sin (\theta) sin(\phi) \\ x &= r \sin (\theta) cos(\phi) \end{aligned} zyx=rcos(θ)=rsin(θ)sin(ϕ)=rsin(θ)cos(ϕ)

    立体角(Solid Angle)

    立体角描述了从原点向一个球面区域张成的视野大小,是平面角在三维的自然推广。立体角的最大值为全角: 4 π 4\pi 4π,该最大值可以在区域为整个球面时取到。
    在这里插入图片描述
    立体角 ω \omega ω 具有如下微分形式:
    d ω = d A r 2 d \omega = \frac{d A}{r^2} dω=r2dA
    由于面积微元在球面坐标系下可以写成:
    d A = ( r d θ ) ( r sin ⁡ θ d ϕ ) = r 2 sin ⁡ θ d θ d ϕ dA = (rd\theta)(r \sin \theta d\phi)=r^2\sin\theta d\theta d\phi dA=(rdθ)(rsinθdϕ)=r2sinθdθdϕ
    因此:
    d ω = d A r 2 = sin ⁡ θ d θ d ϕ d\omega = \frac{dA}{r^2}=\sin \theta d\theta d\phi dω=r2dA=sinθdθdϕ

    投影面积(Foreshortened Area)

    投影面积描述了一个物体表面的微小区域在某个视线方向上的可见面积。对于面积微元 A A A,则沿着与法相夹角为 θ \theta θ 方向的 A A A 的可见面积为:
    a r e a = A cos ⁡ θ area = A \cos \theta area=Acosθ在这里插入图片描述

    光能(Radiant Energy)

    光能表示的是(一个区域中)光子能量的总和。光能通常使用符号 Q Q Q 来表示,其单位是焦耳(J)在这里插入图片描述

    光通量(Radiant Flux)

    光能并不会停留和储存在某个位置,而是在始终不断运动着的。光通量描述的是在单位时间穿过截面的光能(通常用符号 Φ \Phi Φ 来表示,单位是瓦特 W W W
    Φ = d Q d t \Phi = \frac{dQ}{dt} Φ=dtdQ
    在这里插入图片描述

    辉度(Irradiance)

    辉度表示的是物体表面受光能的影响程度,它等于单位面积上的光通量(通常用符号 E E E来表示,单位是 W / m 2 W/m^2 W/m2
    E = d Φ d A E=\frac{d\Phi}{dA} E=dAdΦ
    在这里插入图片描述

    发光强度(Intensity)

    对一个点(例如点光源)来说,发光强度表示单位立体角上的光通量(发光强度简称光强,使用符号 I I I来表示,单位是 W W W
    I = d Φ d ω I=\frac{d\Phi}{d\omega} I=dωdΦ在这里插入图片描述

    光亮度(Radiance)

    光亮度表示物体表面沿某一方向的明亮程度,它等于单位投影面积在单位立体角上的光通量,使用符号 L L L来表示,单位是 W / m 2 W/m^2 W/m2在这里插入图片描述
    一种直观的光亮度的理解方法是:将光亮度理解为物体表面的微面元所接收的来自于某方向光源的单位面积的光通量,因此截面选用垂直于该方向的截面,其面积按投影面积计算。
    在这里插入图片描述
    光亮度的微分形式:
    L = d Φ d A cos ⁡ θ d ω L=\frac{d\Phi}{dA\cos\theta d \omega} L=dAcosθdωdΦ
    光亮度使用物体表面沿目标方向上的投影面积,而不是面积

    辉度和光亮度的关系

    辉度可以写成光亮度在入射光所形成的半球上的积分:
    d Φ d A = E = ∫ Ω L ( ω ) cos ⁡ θ ( ω ) d ω \frac{d\Phi}{dA} = E = \int_{\Omega}L(\omega) \cos \theta(\omega) d\omega dAdΦ=E=ΩL(ω)cosθ(ω)dω

    • Ω \Omega Ω 是入射光所形成的半球
    • L ( ω ) L(\omega) L(ω)是沿 ω \omega ω 方向的光亮度

    BRDF的定义与性质

    BRDF 的定义
    BRDF 的性质:
     - 可逆性(Reciprocity)
     - 能量守恒性质
     - 基于 BRDF 的渲染方程(Rendering Equation)

    BRDF 的定义

    BRDF 描述的是物体表面将光能从任何一个入射方向反射到任何一个视点方向的反射特性。BRDF 模型是绝大多数近代图形学算法中用于描述光反射现象的基本模型。
    BRDF 是关于入射光方向和反射光方向的四维实值函数,它等于反射方向的光亮度和沿入射方向的入射光的辉度之比:
    f ( ω i → ω r ) = d L r ( ω r ) d E i f(\omega_i \rightarrow \omega_r) = \frac{dL_r(\omega_r)}{dE_i} f(ωiωr)=dEidLr(ωr)
    可以写成关于入射光的光亮度的形式:
    f ( ω i → ω r ) = d L r ( w r ) L i ( ω i ) cos ⁡ ω i d ω i f(\omega_i \rightarrow \omega_r)=\frac{dL_r(w_r)}{L_i(\omega_i)\cos \omega_i d\omega_i} f(ωiωr)=Li(ωi)cosωidωidLr(wr)
    在这里插入图片描述

    可逆性(Reciprocity)

    BRDF 的可逆性源自于 Helmholtz 光路可逆性(Helmholtz Reciprocity Rule)
    BRDF 的可逆性是说:交换入射光与反射光的角色,并不会改变 BRDF 的值:
    f ( ω i → ω r ) = f ( ω r → ω i ) f(\omega_i \rightarrow \omega_r)=f(\omega_r \rightarrow \omega_i) f(ωiωr)=f(ωrωi)

    能量守恒性质

    BRDF 需要遵循的另一个物理定律是能量守恒定律。能量守恒定律指出:入射光的能量与出射光的总能量应该相等。
    能量守恒方程如下:
    Q i n c o m i n g = Q r e f l e c t e d + Q a b s o r b + Q t r a n s m i t t e d Q_{incoming}=Q_{reflected}+Q_{absorb}+Q_{transmitted} Qincoming=Qreflected+Qabsorb+Qtransmitted
    由此我们知道:
    Q r e f l e c t e d ≤ Q i n c o m i n g Q_{reflected} \leq Q_{incoming} QreflectedQincoming
    因此 BRDF 必须满足如下的积分不等式,也就是所谓的能量守恒性质:
    ∫ Ω f ( ω i → ω r ) cos ⁡ ω r d ω r ≤ 1 \int_{\Omega}f(\omega_i \rightarrow \omega_r) \cos \omega_r d \omega_r \leq 1 Ωf(ωiωr)cosωrdωr1

    渲染方程(Rendering Equation)

    渲染方程(Rendering Equation)用于计算环境光照明下的反射光的光亮度,它可以写成不同角度下入射光的光亮度乘以 BRDF 的积分:
    L r = ∫ Ω f ( ω i → ω r ) L ( ω i ) cos ⁡ ω i d ω i L_r = \int_{\Omega}f(\omega_i \rightarrow \omega_r)L(\omega_i)\cos \omega_i d\omega_i Lr=Ωf(ωiωr)L(ωi)cosωidωi

    BRDF模型

    为了方便和高效地使用 BRDF 数据,它们往往被组织成为参数化的数值模型。
    BRDF 的数值模型具有如下三类:
      - 经验模型(Empirical Models)
      - 基于物理的模型(Physical-based Models)
      - 数据表达的模型(Data-driven Models)

    • 经验模型,使用基于实验提出的公式对 BRDF 做快速的估计
    • 基于物理的模型,根据物体表面材料的几何以及光学属性建立反射方程,从而计算 BRDF
    • 数据表达的模型,将 BRDF 按照实测数据建立查找表,以便于快速的查找和计算。
    经验模型(Empiriccal Models)

    经验模型提供简洁的公式以便于反射光线的快速计算;
    经验模型不考虑材质特性,仅仅提供一个反射光的粗糙近似。
    经验模型不一定满足物理定律,如 Helmholtz 可逆性或能量守恒定律等。
    经验模型因为其简洁和高效性被广泛运用。

    经验模型1:Lambertian 模型
    • Lambertian模型是最基本的反射模型:
        - 入射光线被均匀地反射到各个方向。
        - 沿不同方向地 BRDF 是一个常数。
    • 反射率(Albedo): ρ \rho ρ
        - 反射率是反射光辉度与入射光辉度之比。
      L r ( ω r ) = ∫ Ω f L i ( ω i ) cos ⁡ ω i d ω i = f ∫ ω L i ( ω i ) cos ⁡ ω i d ω i = f E i \begin{aligned} L_r(\omega_r) &= \int _{\Omega}f L_i(\omega_i)\cos \omega_i d\omega_i\\ &= f\int_{\omega}L_i(\omega_i) \cos \omega_i d\omega_i \\ &= fE_i \end{aligned} Lr(ωr)=ΩfLi(ωi)cosωidωi=fωLi(ωi)cosωidωi=fEi
      ρ = E r E i = ∫ Ω L r ( ω r ) cos ⁡ ω r d ω r E i = f E i π E i = π f \rho = \frac{E_r}{E_i}=\frac{\int_\Omega L_r(\omega_r)\cos \omega_r d\omega_r}{E_i}=\frac{fE_i \pi}{E_i} = \pi f ρ=EiEr=EiΩLr(ωr)cosωrdωr=EifEiπ=πf
      在这里插入图片描述

    该模型能够很好地用于描述那些包含存粹漫反射的物体(例如:纸张)。
    Lambertian 漫反射模型不能表现出材质的镜面反射效果,而镜面反射对于金属材质非常重要。
    由于 Lambertian 模型的简洁以及对漫反射良好的描述特性,它常常在其它的经验模型(例如:Phong 模型)中作为分量的形式被包含。

    经验模型2:Phong模型

    Phong 模型是在 Lambertian 漫反射模型的基础上,添加了一项镜面反射项,以表达在反射角上的镜面反射效果:

    f ( l → v ) = ρ d + ρ s ( r ⋅ v ) s ( n ⋅ l ) f(l\rightarrow v) = \rho_d + \rho_s \frac{(r\cdot v)^s}{(n\cdot l)} f(lv)=ρd+ρs(nl)(rv)s
    其中, ρ d \rho_d ρd ρ s \rho_s ρs 分别表示漫反射光和镜面反射光线的反射率, s s s表示发光指数,用于描述镜面反射的锐利度。

    Phong 模型不满足可逆性:
    f ( l → v ) ≠ f ( v → l ) f(l\rightarrow v) \neq f(v\rightarrow l) f(lv)=f(vl)
    在这里插入图片描述
    尽管 Phong 模型缺乏物理解释,并且对于某些金属材质,它并不十分准确;但 Phong 模型仍是目前计算机图形学中被最为广泛使用的基于经验的反射模型。

    Phong 模型的优势在于它的简洁和高效性,以及能够同时表现漫反射和镜面反射的特征。

    Phong 模型的扩展

    Phong 模型的大部分扩展是为了实现进一步的加速。

    Blinn-Phong 模型
      - 通过引入入射方向 l l l 和视点方向 v v v 的角平分线 h h h,使用 h h h 和法向 n n n 的点积替代原先 Phong 模型中 r r r v v v 的点积,可以简化运算:
    f ( l → v ) = ρ d + ρ s ( n ⋅ h ) s ( n ⋅ l )    w h e r e   h = v + l 2 f(l\rightarrow v) = \rho_d + \rho_s \frac{(n\cdot h)^s}{(n\cdot l)}\ \ where\ h=\frac{v+l}{2} f(lv)=ρd+ρs(nl)(nh)s  where h=2v+l

    快速 Phong 绘制(Fast Phong Shading)
    由 Bishop 和 Weimer 发表在 1986年 的 ACM SIGGRAPH 年会上。
    思想是使用制表和插值的方法对指数项 ( r ⋅ v ) s (r\cdot v)^s (rv)s进行快速计算,从而实现绘制的加速。

    可逆的 Phong 模型(Modified Phong Model)
    通过去掉 Phong 模型中镜面反射分量的分母项,从而使得修改后的 Phong 模型能够满足可逆性:
    f ( l → v ) = ρ d + ρ s ( r ⋅ v ) s f ( l → v ) = f ( v → l ) \begin{aligned} f(l \rightarrow v) &= \rho_d+\rho_s(r\cdot v)^s\\ f(l \rightarrow v) &= f(v \rightarrow l) \end{aligned} f(lv)f(lv)=ρd+ρs(rv)s=f(vl)
    Phong模型的示例效果
    在这里插入图片描述

    物理模型(Physical-Based Models)

      - 经验模型(Empirical Models)源于设计者的直觉和实践经验,而物理模型则建立在有关光的相互作用的科学知识上。通过包含材料的各种几何及光学性质来尽可能精确的近似现实世界中的材料。
      - 物理模型(Physical-Based Models)通常建立在被成为表面粗糙度(surface roughness)的细节几何结构上。
       粗糙度:从微观角度来看,几乎没有完全光滑的表面;微观尺度的表面几何是通过一组微平面(microfacets)集合来建模的;粗糙度通过微平面法向的统计分布来表达的
    在这里插入图片描述
    菲涅尔项(Fresnel Term)
    在实际应用中,我们发现,单向反射性在擦地角(grazing angles)附近增大
    在这里插入图片描述
    入射光的反射量是由从麦克斯韦电磁波方程组中得到的菲涅尔公式得到的。

    菲涅尔项(Fresnel Term)
    在这里插入图片描述
    定义如下:
    F s = a 2 + b 2 − 2 a cos ⁡ θ + cos ⁡ 2 θ a 2 + b 2 + 2 a cos ⁡ θ + cos ⁡ 2 θ F p = F s a 2 + b 2 − 2 a sin ⁡ θ tan ⁡ θ + sin ⁡ 2 θ tan ⁡ 2 θ a 2 + b 2 + 2 a sin ⁡ θ tan ⁡ θ + sin ⁡ 2 θ tan ⁡ 2 θ 2 a 2 = ( η 2 − κ 2 − sin ⁡ 2 θ ) 2 + 4 η 2 κ 2 + ( η 2 − κ 2 − sin ⁡ 2 θ ) 2 b 2 = ( η 2 − κ 2 − sin ⁡ 2 θ ) 2 + 4 η 2 κ 2 − ( η 2 − κ 2 − sin ⁡ 2 θ ) η = η m η a + κ m κ a η a 2 + κ a 2 κ = η a κ m − η m κ a η a 2 + κ a 2 \begin{aligned} F_s &= \frac{a^2+b^2-2a\cos\theta +\cos^2\theta}{a^2+b^2+2a\cos\theta +\cos^2\theta}\\ F_p &= F_s\frac{a^2+b^2-2a\sin\theta\tan\theta+\sin^2\theta\tan^2\theta}{a^2+b^2+2a\sin\theta\tan\theta+\sin^2\theta\tan^2\theta} \\ 2a^2 &= \sqrt{(\eta^2-\kappa^2-\sin^2\theta)^2+4\eta^2\kappa^2} +(\eta^2 - \kappa^2 - \sin^2 \theta) \\ 2b^2 &= \sqrt{(\eta^2-\kappa^2-\sin^2\theta)^2+4\eta^2\kappa^2} -(\eta^2 - \kappa^2 - \sin^2 \theta)\\ \eta &= \frac{\eta_m\eta_a+\kappa_m\kappa_a}{\eta^2_a+\kappa^2_a}\\ \kappa &= \frac{\eta_a\kappa_m-\eta_m \kappa_a}{\eta^2_a+\kappa^2_a} \end{aligned} FsFp2a22b2ηκ=a2+b2+2acosθ+cos2θa2+b22acosθ+cos2θ=Fsa2+b2+2asinθtanθ+sin2θtan2θa2+b22asinθtanθ+sin2θtan2θ=(η2κ2sin2θ)2+4η2κ2 +(η2κ2sin2θ)=(η2κ2sin2θ)2+4η2κ2 (η2κ2sin2θ)=ηa2+κa2ηmηa+κmκa=ηa2+κa2ηaκmηmκa
    菲涅耳反射率: F = F p + F s 2 F=\frac{F_p+F_s}{2} F=2Fp+Fs

    Cook-Torrance 模型

    由 Cook 和 Torrance 提出[ A Reflectance Model for Computer Graphics, ACM SIGGRAPH, 1981]
    在图形学中使用最早的 BRDF 物理模型
    是应用物理学家开发的 Torrance-Sparrow 模型的一个应用版本。[Theory for Off-Specular Reflection from Roughened Surfaces, J.Optical Society of America, 1975]

    假设微平面是镜面反射的
    微平面被假定按如图所示 V V V形沟槽排列。
    在这里插入图片描述
    结合 Lambertian 漫反射项与微平面反射的镜面反射项:
    f λ = d ρ d π + s F λ D G 4 π ( n ⋅ l ) ( n ⋅ v ) f_\lambda = d\frac{\rho_d}{\pi}+s F_\lambda\frac{DG}{4\pi(n\cdot l)(n\cdot v)} fλ=dπρd+sFλ4π(nl)(nv)DG
    这里的 F λ F_\lambda Fλ 是菲涅尔因子, D D D是微平面法向的分布函数, G G G 是几何衰减因子, s s s, d d d 是镜面反射和漫反射系数。

    微平面法向的分布函数 D D D
    由于微平面只有镜面反射,只有那些法向沿着平分线方向 h h h的才对镜面反射起作用。
    采用 Bechmann 分布来描述:
    D ( h ) = 1 m 2 cos ⁡ 4 β e − tan ⁡ 2 β m 2 D(h) = \frac{1}{m^2\cos^4\beta}e^{-\frac{\tan^2\beta}{m^2}} D(h)=m2cos4β1em2tan2β
    其中, m m m 是表面的粗糙度值, β \beta β n n n h h h 的夹角
    在这里插入图片描述
    几何衰减系数 G G G(Geometric Attenuation Factor)
    考虑来自视角方向的遮挡效果和来自光线方向的阴影效果
    用于排除被遮挡的(masked)和阴影覆盖的(shadowing)微平面,G定义如下:
    G ( n , l , v ) = min ⁡ { 1 , 2 ( n ⋅ h ) ( n ⋅ v ) ( v ⋅ h ) , 2 ( n ⋅ h ) ( n ⋅ l ) ( v ⋅ h ) } G(n,l,v)=\min \left \{1, \frac{2(n\cdot h)(n\cdot v)}{(v\cdot h)}, \frac{2(n\cdot h)(n\cdot l)}{(v\cdot h)} \right\} Gn,l,v)=min{1,(vh)2(nh)(nv),(vh)2(nh)(nl)}
    在这里插入图片描述
    这种对微平面遮挡和阴影的考虑使得可以产生某些反射特性
    一个是反射方向不在镜面方向的反射(off-specular reflection)
    另一个是逆反射(retroreflection)
    在这里插入图片描述
    反射方向不在镜面方向的反射(Off-specular reflection)
      - 反射方向不在镜面反射方向上
      - 这是典型的粗糙表面的特征
    逆反射(Retroreflection)
      - 很大一部分光线沿着反方向返回
      - 例如下图中的满月图片,边缘看起来和中心一样明亮
    在这里插入图片描述
    Cook-Torrance模型与Phong模型的结果对比
    在这里插入图片描述
    更多的 Cook-Torrance 模型示例结果
    在这里插入图片描述

    其它的物理模型

    BRDF可分为两类

    • 各向同性(Isotropic)
        - 反射不受与给定表面法向夹角的约束
        - 随机表面微结构
    • 各项异性(Anisotropic)
        - 反射比随着与某个给定的表面法向之间的夹角而变化
        - 图案的表面微结构
        - 金属丝,绸缎,毛发
    Ward 模型

    然而,Phong 和 Cook-Torrance BRDF 模型都不能处理各项异性的效果
    现在我们介绍另一种 BRDF 模型:Ward 模型
    在这里插入图片描述
     - Ward 模型由 Ward 于 1992 年提出[Measuring and Modeling Anisotropic Reflection, ACM SIGGRAPH, 1992]
     - 介绍了一种更一般的表面法向表达方式:通过椭圆体(ellispsoids)这种允许各向异性反射的形式来表达
     - 然而,由于没有考虑菲涅尔因子(Fresnel factor)和几何衰减因子(Geometric attenuation factor),该模型更像是一种经验模型

    各向同性的 Ward 模型定义为:
    ρ b d   i s o ( θ i , ϕ i , θ i , ϕ r ) = ρ d π + ρ s ⋅ 1 cos ⁡ θ i cos ⁡ θ r ⋅ e − tan ⁡ 2 σ α 2 4 π α 2 \rho _{bd\ iso}(\theta_i,\phi_i, \theta_i,\phi_r)=\frac{\rho_d}{\pi}+\rho_s \cdot \frac{1}{\sqrt{\cos \theta_i \cos \theta_r}}\cdot \frac{e^{-\tan^2{\frac{\sigma}{\alpha^2}}}}{4\pi \alpha^2} ρbd iso(θi,ϕi,θi,ϕr)=πρd+ρscosθicosθr 14πα2etan2α2σ
    其中:
    ρ d \rho_d ρd是漫反射系数
    ρ s \rho_s ρs是镜面反射系数
    σ \sigma σ n ^ \hat{n} n^ h ^ \hat{h} h^的夹角
    α \alpha α 是表面坡度的标准差

    将菲涅尔因子和几何衰减因子替换为一个用于保证分布在整个半球内积分的简单归一化项。

    各向异性的 Ward 模型定义为:
    ρ b d ( θ i , ϕ i , θ r , ϕ r ) = ρ d π + ρ s ⋅ 1 cos ⁡ θ i cos ⁡ θ r ⋅ exp ⁡ [ − tan ⁡ 2 σ ( cos ⁡ 2 ϕ / α x 2 + sin ⁡ 2 ϕ / α y 2 ) ] 4 π α x α y \rho_{bd}(\theta_i,\phi_i,\theta_r,\phi_r)=\frac{\rho_d}{\pi}+\rho_s\cdot \frac{1}{\sqrt{\cos \theta_i\cos\theta_r}}\cdot \frac{\exp[-\tan^2\sigma(\cos^2\phi/\alpha_x^2 + \sin^2\phi / \alpha_y^2)]}{4\pi\alpha_x\alpha_y} ρbd(θi,ϕi,θr,ϕr)=πρd+ρscosθicosθr 14παxαyexp[tan2σ(cos2ϕ/αx2+sin2ϕ/αy2)]
    其中:

    ρ d \rho_d ρd是漫反射系数
    ρ s \rho_s ρs是镜面反射系数
    σ \sigma σ n ^ \hat{n} n^ h ^ \hat{h} h^的夹角
    α x \alpha_x αx x ^ \hat{x} x^方向表面坡度的标准差
    α y \alpha_y αy y ^ \hat{y} y^方向表面坡度的标准差
    ϕ \phi ϕ h ^ \hat{h} h^投影到表面的方位角

    Word 模型示例结果
    在这里插入图片描述

    更多的BRDF类模型

    Oren-Nayar模型[Generalization of Lambert’s Reflectance Model ACM SIGGRAPH, 1994]
      - 将微平面当做 Lambertian 反射体而不是镜子
    Poulin-Fournier 模型[ A Model for Anisotropic Reflection, ACM SIGGRAPH, 1990]
      - 将微平面的法向表示为一组平行的圆柱体,同样处理各项异性的情况
    波动光学(Wave optic)相关模型
      - 由波动光学原理发展而来
      - 材料的微平面大小和光的波长相当
      - 基于不同的衍射(diffraction)理论,有两个工作[He, etal., A Comprehensive Physical Model for Light Reflection, ACM SIGGRAPH, 1991] 和 [Stam, Diffraction Shaders, ACM SIGGRAPH, 1999]
      -虽然有很强的描述能力,但由于模型本身过于复杂而限制了它们的应用
    由波动光学相关的模型得到的效果:
    在这里插入图片描述

    数据驱动的模型

    度量一个大的材料集合的 BRDF,并将其记录为高维向量
    利用降维方法从这些数据中计算一个低维流型
    代表性工作:Matusik et al. [A Data-Driven Reflectance Model, ACM SIGGRAPH, 2003]

    BRDF 模型表达形式的对比

    经验模型(Empirical Models):
      - 计算简单,可以产生视觉上可接受的结果
    物理模型(Physical-Based Models):
      - 参数具有物理意义
      - 基于科学知识构建
    数据驱动模型(Data-driven Models):
      - 灵活,对材料属性没有假设限定
      - 很大的数据集,通常需要数据降维方法来压缩数据
      - Needs interpolating data, missing high light parts

    BRDF度量、评价

    动机
      - 为了对具有未知反射属性的材料进行建模,并生成具有高度真实感的结果
      - 恢复 BRDF 和其它场景属性的过程有时也被称作逆渲染
      - BRDF 度量概述
        - 度量设备
        - 实际获取中的问题
    度量设备
      - 由于 BRDF 是一个关于光线和视角的函数,其度量可通过对 2D 光照空间和2D视角空间进行采样获取。
    在这里插入图片描述
    测量设备:Light Stage
      - Debevec et al. [SIGGRAPH 2000]
      - 从固定物体(或人)上采集 BRDF
      - 从少数的几个视角拍摄 2000 多个不同光源方向条件下的图像
      - 从所有图像点上同时获取清楚的BRDF
    在这里插入图片描述
    测量设备:Gonioreflectometer
    在这里插入图片描述
    在这里插入图片描述

    展开全文
  • 光照模型和BRDF

    2021-11-02 23:49:51
    本文目录一、预备知识二、BRDF三、Phong光照模型 本文参考了本篇博客,同时加入了Phong光照模型 一、预备知识 在介绍光照模型和BRDF(双向反射分布函数)时,我们要先理解一些基本概念,这些概念与辐射度学有关 1. ...

    本文参考了本篇博客,同时加入了Phong光照模型

    一、预备知识

    在介绍光照模型和BRDF(双向反射分布函数)时,我们要先理解一些基本概念,这些概念与辐射度学有关

    1. 能量

    每个光子都具有一定的能量(Energy),能量大小与频率有关 E = h v E = hv E=hv。用符号 Q Q Q表示,单位是焦耳 J J J

    2. 功率

    功率、辐射通量、通量

    功率(Power),也被称为辐射通量(Radiant Flux)或者通量(Flux),指的是单位时间内通过表面或者空间区域的能量的总量。用符号 Φ \Phi Φ表示, Φ = d Q d t \Phi = \frac{dQ}{dt} Φ=dtdQ,单位为瓦特 W W W,或者焦耳/秒 J / S J/S J/S

    3. 辐照度和辐出度

    辐照度(辐射通量密度)

    辐出度、辐射出射度、辐射度(辐射通量密度)

    • 辐照度(Irradiance)指单位时间内到达单位面积的辐射能量,或单位面积的辐射通量,即通量对于面积的密度。用符号 E E E表示, E = d Φ d A = d Q d t d A E = \frac{d\Phi}{dA} = \frac{dQ}{dtdA} E=dAdΦ=dtdAdQ,单位 W / m 2 W/m^{2} W/m2
    • 辐出度(Radiant Existance),也称为辐射出射度、辐射度(Radiosity)。与辐照度不同的是,辐出度是衡量离开表面的通量密度,二者都可以被称为辐射通量密度(Radiant Flux Density)。用符号 M M M表示

    注意辐射通量和辐射通量密度之间的区别

    如下图所示

    一束平行光之间的间距为 d d d,一个垂直入射,一个与表面法线成 θ i \theta_{i} θi的角度入射,观察表面的间距,其中左图间距为 d d d,右图为 d c o s θ i \frac{d}{cos\theta_{i}} cosθid,相对于左图,右图的间距增大了。由于辐射通量是单位时间内通过表面或空间区域内的能量总和,所以辐射通量相同;而辐射通量密度是单位时间内通过单位面积的能量,右边面积大,所以右边辐射通量密度变小

    让我们来看一下二者辐射通量密度的差异:假设右图中不垂直于光线方向的表面面积为 A A A,将他投影到垂直于直线方向得到一个虚拟表面,面积为 A ⊥ = A c o s θ i A^{\perp } = Acos\theta_{i} A=Acosθi,通过这两个面积的通量相同,设为 Φ \Phi Φ,则表面接收到的辐射通量密度 E = Φ A E = \frac{\Phi}{A} E=AΦ,虚拟表面的辐射通量密度 E ⊥ = Φ A ⊥ = Φ A c o s θ i E_{\perp} = \frac{\Phi}{A^{\perp}} = \frac{\Phi}{Acos\theta_{i}} E=AΦ=AcosθiΦ,得到 E = E ⊥ c o s θ i E = E_{\perp}cos\theta_{i} E=Ecosθi

    同理对于点光源,假想以点光源为中心不同半径的求包围着点光源,穿过这些球的辐射通量是相通的,均为 Φ \Phi Φ,球的表面积为 4 π r 2 4\pi r^{2} 4πr2,得到辐射通量密度 E = Φ 4 π r 2 E = \frac{\Phi}{4\pi r^{2}} E=4πr2Φ,即通量密度与距离的平方成反比,即光的衰减与距离的平方成反比

    4. 辐射强度

    辐射强度

    立体角:度量三维角度的量,用符号 ω \omega ω表示,立体角计算为 ω = s r 2 \omega = \frac{s}{r^{2}} ω=r2s,单位球的表面积为 4 π r 2 = 4 π 4\pi r^{2} = 4\pi 4πr2=4π,所以球面的立体角也为 4 π 4\pi 4π

    辐射强度(Radiant Intensity),通过单位立体角的辐射能量,用符号 I I I表示, I = d Φ d w I=\frac{d\Phi}{dw} I=dwdΦ,单位 W / s r W/sr W/sr

    引入辐射强度的原因是:有时候需要度量通过一个点的通量的密度,因为点的面积是0,无法直接使用辐射通量密度,所以引入辐射强度。辐射强度不会随距离的增大而衰减,这是因为立体角不会随距离的变化而变化

    5. 辐射率

    辐射率

    辐射率(Radiance),每单位面积每单位立体角的辐射通量密度,用符号 L L L表示, L = d Φ d ω d A ⊥ L = \frac{d\Phi}{d\omega dA^{\perp}} L=dωdAdΦ,单位为 W / m 2 s r W/m^{2}sr W/m2sr,其中 d A ⊥ dA^{\perp} dA是微分面积 d A dA dA在垂直于光线方向的投影,如图

    辐射率可以看成我们眼睛看到(相机拍到)的物体上一点的颜色,基于物理着色时,计算表面一点的颜色就是计算它的辐射率,辐射率不会随距离变化而衰减,这和我们日常感受一致,在没有雾霾的干扰时,我们看到的物体表面上一点的颜色并不会随距离变化而变化。为什么辐照度会随距离增大而衰减,但是我们看到的颜色却不会衰减呢?这是因为随着距离变大,我们看到的物体上的一块区域到达视网膜的通量密度会变小,同时这块区域在视网膜表面上的立体角也会变小,正好抵消了通量密度的变化

    二、BRDF

    我们看到一个表面,实际上是周围环境的光照射到表面上,然后表面将一部分光反射到我们眼睛里。双向反射分布函数BRDF(Bidirectional Reflectance Distribution Function)就是描述表面入射光和反射光关系的

    对于一个方向的入射光,表面会将光反射到表面上半球的各个方向,不同方向反射的比例是不同的,我们用BRDF来表示指定方向的反射光和入射光的比例关系,BRDF定义为:

    f ( l , v ) = d L o ( v ) d E ( l ) f(l, v) = \frac{dL_{o}(v)}{dE(l)} f(l,v)=dE(l)dLo(v)

    • f f f:BRDF, l l l是入射光方向, v v v是观察方向,即反射光方向
    • d L o ( v ) dL_{o}(v) dLo(v):表面反射到 v v v方向的反射光的微分辐射率。表面反射到 v v v方向的反射光的辐射率为 L o ( v ) L_{o}(v) Lo(v),来自于表面上半球所有方向的入射光线的贡献,而微分辐射率 d L o ( v ) dL_{o}(v) dLo(v)特指来自方向 l l l的入射光贡献的反射辐射率
    • d E ( l ) dE(l) dE(l):表面上来自入射光方向 l l l的微分辐照度。表面接收到的辐照度为 E E E,来自上半球所有方向的入射光线的贡献,而微分辐照度 d E ( l ) dE(l) dE(l)特指来自于方向 l l l的入射光

    表面对不同频率的光反射率可能不一样,因此BRDF和光的频率有关。在图形学中,将BRDF表示为RGB向量,三个分量各有自己的 f f f函数

    BRDF需要处理表面上半球的各个方向,如下图使用球坐标系定义方向更加方便。球坐标系使用两个角度来确定一个方向:

    • 方向相对法线的角度 θ \theta θ,称为极角(Polar Angle)或天顶角(Zenith Angle)
    • 方向在平面上的投影相对于平面上一个坐标轴的角度 ϕ \phi ϕ,称为方位角(Azimuthal Angle)

    所以BRDF也可以表示成 f ( θ i , ϕ i , θ o , ϕ o ) f(\theta_i, \phi_i, \theta_o, \phi_o) f(θi,ϕi,θo,ϕo)。对于各向同性材质,当 l l l v v v同时绕法线 n n n旋转时, f f f值保持不变,此时可以用 l l l v v v在平面投影的夹角 ϕ \phi ϕ来代替 ϕ i \phi _i ϕi ϕ o : f ( θ i , θ o , ϕ ) \phi _o:f(\theta_i, \theta_o, \phi) ϕof(θi,θo,ϕ)

    为什么BRDF要定义成辐射率和辐照度的比值,而不是直接定义为辐射率和辐射率比值,有两种解释

    第一种解释参考BRDF为什么要定义为一个单位是sr-1的量?

    结合下面辐射通量密度 ( A ) (A) (A)和辐射率 ( B ) (B) (B)测量仪的示意图来看:辐照度(辐射通量密度)测量仪 ( A ) (A) (A)接受平面上半球的所有光线,可以测量一个较小面积来自于四面八方的所有光通量,光通量 Φ \Phi Φ除以传感器面积 A A A就可以得到辐照度(辐射通量密度) E E E。辐射度测量仪 ( B ) (B) (B)则有一个长筒控制光线只能从一个很小的立体角进入测量仪,光通量 Φ \Phi Φ除以传感器面积 A A A和立体角 ω \omega ω就可以得到辐射率 L L L

    测平面上一点在某一个方向的出射辐射率很简单,只需要用仪器 ( B ) (B) (B)从该方向对准该点就可以了。而测平面一点入射的辐射率则没有那么简单,必须保证光源正好覆盖测量仪开口立体角,大了该点会接受到比测量值更多的光照,导致测量值比实际值小,小了则与仪器的设计立体角不一致,可在实际中是基本做不到光源大小正好覆盖测量仪开口立体角的。而测表面的辐照度则简单得多,只要保证光源很小,而且没有来自其他方向的光干扰,这时候测到的辐照度就是平面上来自光源方向的微分辐照度 d E dE dE

    第二种解释从数学的角度出发

    对于现实世界中的非光学平面,一束光线射到表面上后,被表面反射到各个方向,其中一个出射方向的光通量只是整个反射光通量极小的一部分,当出射方向立体角趋于0时, lim ⁡ ω o → 0 d L o L i = 0 \lim_{\omega_{o} \rightarrow 0}{\frac{dL_{o} }{Li} } = 0 limωo0LidLo=0,所以在实际计算中使用辐射率和辐射率比值是没有意义的。而如果分母改成表面上接收到的来自光源方向的微分辐照度,我们知道 d E = L i ( l ) d ω i c o s θ i dE = L_i(l) d\omega _{i} cos \theta _{i} dE=Li(l)dωicosθi,由于给入射辐射率乘了一个趋于零的微分立体角, d E dE dE的值会小很多,比值 d L o d E \frac{dL_o}{dE} dEdLo是有意义的,而不是0

    下面我们来看看怎么用BRDF来计算表面辐射率

    我们考虑来自方向 l l l的入射光辐射率 L i ( l ) L_i(l) Li(l),由辐射率和辐照度的定义:

    L i ( l ) = d Φ d ω i d A ⊥ = d Φ d ω i d A c o s θ i = d E ( l ) d ω i c o s θ i L_i(l) = \frac{d \Phi }{d\omega_i d A^{\bot } } = \frac{d \Phi }{d\omega_i dA cos \theta_i } = \frac{dE(l) }{d\omega_i cos \theta_i } Li(l)=dωidAdΦ=dωidAcosθidΦ=dωicosθidE(l)

    则照射到表面来自于方向l的入射光贡献的微分辐照度:

    d E ( l ) = L i ( l ) d ω i c o s θ i dE(l) = L_i(l) d\omega_i cos \theta_i dE(l)=Li(l)dωicosθi

    表面反射到 v v v方向的由来自于方向l的入射光贡献的微分辐射率:

    d L o ( v ) = f ( l , v ) ⊗ d E ( l ) = f ( l , v ) ⊗ L i ( l ) d ω i c o s θ i dL_o(v) = f(l, v) \otimes dE(l) = f(l, v) \otimes L_i(l) d\omega_i cos \theta_i dLo(v)=f(l,v)dE(l)=f(l,v)Li(l)dωicosθi
    符号 ⊗ \otimes 表示按向量的分量相乘,因为 f f f L i L_i Li都包含RGB三个分量

    要计算表面反射到v方向的来自上半球所有方向入射光线贡献的辐射率,可以将上式对半球所有方向的光线积分:

    L o ( v ) = ∫ Ω f ( l , v ) ⊗ L i ( l ) c o s θ i d ω i L_o(v) = \int_{\Omega }^{} f(l, v) \otimes L_i(l) cos \theta_i d\omega_i Lo(v)=Ωf(l,v)Li(l)cosθidωi
    上式称为反射方程(Reflectance Equation),用来计算表面反射辐射率

    对于点光源、方向光等理想化的精准光源(Punctual Light),计算过程可以大大简化。我们考察单个精准光源照射表面,此时表面上的一点只会被来自一个方向的一条光线照射到(而面积光源照射表面时,表面上一点会被来自多个方向的多条光线照射到),则辐射率:

    L o ( v ) = f ( l , v ) ⊗ E L c o s θ i L_o(v) = f(l, v) \otimes E_L cos \theta_i Lo(v)=f(l,v)ELcosθi
    对于多个精准光源,只需简单累加就可以了:

    L o ( v ) = ∑ k = 1 n f ( l k , v ) ⊗ E L k c o s θ i k L_o(v) = \sum_{k = 1}^{n}{f(l_k, v) \otimes E_{L_k} cos \theta_{i_k}} Lo(v)=k=1nf(lk,v)ELkcosθik
    这里使用光源的辐照度,对于阳光等全局方向光,可以认为整个场景的辐照度是一个常数,对于点光源,辐照度随距离的平方衰减,用公式 E L = Φ 4 π r 2 E_{L} = \frac{\Phi }{4\pi r^2} EL=4πr2Φ就可以求出到达表面的辐照度, Φ \Phi Φ是光源的功率,比如100瓦的灯泡,r是表面离光源的距离

    回头看看反射方程,是对表面上半球所有方向的入射光线积分,这里面包含了来自精准光源的光线,也包括周围环境反射的光线。处理来自周围环境的光线可以大幅提高光照的真实程度,在实时图形学中,这部分光照可以用基于图像的光照(Image Based Lighting)来模拟。我们将在下篇文章讨论基于图像的光照。

    上面给出了BRDF的定义和使用BRDF计算表面反射辐射率的公式。但这个定义实际上是无法直接用于计算表面反射辐射率的,我们还要建立一个能模拟真实光照的模型,使得输入入射方向和出射方向, f ( l , v ) f(l, v) f(l,v)能输出表面反射微分辐射率和入射微分辐照度的比率。

    1967年Torrance-Sparrow在Theory for Off-Specular Reflection From Roughened Surfaces中使用辐射度学和微表面理论建立了模拟真实光照的BRDF模型,1981年Cook-Torrance在A Reflectance Model for Computer Graphics中把这个模型引入到计算机图形学领域,现在这个模型已经成为基于物理着色的标准,被称为Cook-Torrance模型。下面我们来看看微表面理论和Cook-Torrance模型的推导过程。

    三、Phong光照模型

    光照到物体表面时,物体对光会发生反射、透射和吸收。简单光照模型只考虑物体对直接光照的反射作用,而物体间的光反射作用,只用环境光来表示。此时光源被假定为点光源,反射作用被细分为镜面反射和漫反射

    Phong光照模型

    根据Phong光照模型,物体被感知的亮度由环境光、漫反射光及镜面反射光组成,物体表面的反射光强度为
    I = K a I a + K d I l c o s ( θ ) + K s I l c o s n ( α ) I = K_{a}I_{a} + K_{d}I_{l}cos(\theta) + K_{s}I_{l}cos^{n}(\alpha) I=KaIa+KdIlcos(θ)+KsIlcosn(α)

    • K a   K d   K s K_{a}~K_{d}~K_{s} Ka Kd Ks:分别为环境光、漫反射光及镜面反射光系数
    • I a I_{a} Ia:入射的环境光强度
    • I l I_{l} Il:光源的入射光强度
    • θ \theta θ:入射光与物体表面法线之间的夹角
    • α \alpha α:反射光与视线之间的夹角
    • n n n:镜面反射参数

    1. 漫反射

    理想漫反射当光源来自一个方向时,漫反射光均匀向各方向传播,与视点无关,他是由表面的粗糙不平引起的,因而漫反射光的空间分布是均匀的。记入射光强为 I p I_{p} Ip,物体表面上点 P P P的法向量为 N N N,从点 P P P指向光源的向量为 L L L,当 L 、 N L、N LN为单位向量时,理想漫反射光强表示为:
    I d = I p K d ( L ⋅ N ) I_{d} = I_{p}K_{d}(L\cdot N) Id=IpKd(LN)

    K d K_{d} Kd:与物体有关的漫反射系数,取值 ( 0 , 1 ) (0, 1) (0,1)。在RGB模型下,漫反射系数 K d K_{d} Kd有三个分量 K d r 、 K d g 、 K d b K_{dr}、K_{dg}、K_{db} KdrKdgKdb分别表示RGB三原色的漫反射系数,反应物体颜色,通过调整三个分量可以改变物体的颜色,也可以吧入射光强 I I I设为三个分量 I r 、 I g 、 I b I_{r}、I_{g}、I_{b} IrIgIb通过调整这些分量来调整光源的颜色

    2. 镜面反射

    镜面反射光对于理想镜面,反射光集中在一个方向,并遵守反射定律。对一般光滑表面,反射光集中在一个范围内,且由反射定律决定的反射方向光强最大。因此对于同一点来说,从不同位置所观察到的镜面反射光强是不同的,镜面反射光强可以表示为:
    I s = I p K s c o s n α     α ∈ ( 0 , π 2 ) I_{s} = I_{p}K_{s}cos^{n}\alpha~~~\alpha \in (0, \frac{\pi}{2}) Is=IpKscosnα   α(0,2π)

    • K s K_{s} Ks:与物体有关的镜面反射系数
    • α \alpha α:视线方向 V V V与反射方向 R R R的夹角
    • n n n:反射指数,反映了物体表面的光泽程度,一般为1~2000,数目越大表示物体表面越光滑

    镜面反射光将会在反射方向附近形成很亮的光斑,称为高光现象,将 V 和 R V和R VR看做单位向量,镜面反射光强可以重写为:
    I s = I p K s ( R ⋅ V ) n I_{s} = I_{p}K_{s}(R\cdot V)^{n} Is=IpKs(RV)n

    其中 R = 2 N ( N ⋅ L ) − L R = 2N(N\cdot L)-L R=2N(NL)L,镜面反射光产生的高光区域只反映光源的颜色,例如在红光照射下,一个物体的高光域是红光,镜面反射系数 K s K_{s} Ks是一个与物体颜色无关的参数,在简单光照明模型中,只能通过改变物体的漫反射系数来控制物体的颜色

    3. 环境光

    环境光指光源间接对物体的影响,是在物体和环境之间多次反射,最终达到平衡时的一种光,可以近似地认为同一环境下的环境光,其光强分布是均匀的,他在任何一个方向上的分布都相同,在简单光照模型中,用一个常数来模拟环境光,用式子表示为:
    I e = I a K a I_{e} = I_{a}K_{a} Ie=IaKa

    展开全文
  • BRDF理论及shader实现(下)

    千次阅读 2020-10-06 15:14:40
    BRDF理论及shader实现(上) Specular BRDF 对于specular分量来说,fmf_mfm​是一个遵循菲涅尔反射定律的镜面BRDF项,此时的fmf_mfm​满足([3]和[21]有详细的推导): fm(l,v,m)=F(v,m)δωm(h,m)4(l⋅h)2 f_m({\bf...

    接上篇:
    BRDF理论及shader实现(上)

    Specular BRDF

    对于specular分量来说, f m f_m fm是一个遵循菲涅尔反射定律的镜面BRDF项,此时的 f m f_m fm满足([3]和[21]有详细的推导):

    f m ( l , v , m ) = F ( v , m ) δ ω m ( h , m ) 4 ( l ⋅ h ) 2 f_m({\bf{l}},{\bf{v}},{\bf{m}}) = F({\bf{v}},{\bf{m}})\frac{\delta_{\omega_m}({\bf{h}}, {\bf{m}})}{4({\bf{l}}\cdot{\bf{h}})^2} fm(l,v,m)=F(v,m)4(lh)2δωm(h,m)

    h {\bf{h}} h表示half vector,是 v {\bf{v}} v l {\bf{l}} l的平均;这里分母上第一次出现了 4 4 4,这也是后面specular BRDF公式的分母上的 4 4 4的来源。 δ ω m ( s , o ) \delta_{\omega_m}({\bf{s}}, {\bf{o}}) δωm(s,o)记为:

    ∫ Ω g ( o ) δ ω o ( s , o ) d ω o = { g ( s ) , if  s ∈ Ω 0 , otherwise \int_\Omega g({\bf{o}})\delta_{\omega_o}({\bf{s}}, {\bf{o}})d\omega_o = \begin{cases} g({\bf{s}}), & \text {if $s\in\Omega$} \\ 0, & \text{otherwise} \end{cases} Ωg(o)δωo(s,o)dωo={g(s),0,if sΩotherwise

    此时, f r f_r fr可以化简为为:

    f r ( l , v ) = D ( h , α ) G ( v , l , α ) F ( v , h , F 0 ) 4 ( n ⋅ v ) ( n ⋅ l ) f_r({\bf{l}},{\bf{v}}) = \frac{D({\bf{h}},\alpha)G({\bf{v}},{\bf{l}},\alpha)F({\bf{v}},{\bf{h}},F_0)}{4({\bf{n}}\cdot{\bf{v}})({\bf{n}}\cdot{\bf{l}})} fr(l,v)=4(nv)(nl)D(h,α)G(v,l,α)F(v,h,F0)

    其中, α \alpha α取决于我们前面提到的粗糙度roughness,具体为

    α = r o u g h n e s s 2 \alpha = roughness^2 α=roughness2

    可以理解为 α \alpha α是对粗糙度roughness的一个映射, α \alpha α将多次被用到。

    菲涅尔项及 F 0 F_0 F0会在后面详细介绍,这里暂时略过。

    shader实现:

    // #define saturate(x) clamp(x, 0, 1)
    // N     = normal;
    // V     = normalize(camPos - WorldPos);
    // L     = normalize(LightPos - WorldPos));
    // H     = normalize(V + L);
    // NdotV = saturate(dot(N, V));
    // NdotL = saturate(dot(N, L));
    // NdotH = saturate(dot(N, H));
    // LdotH = saturate(dot(L, H));
    // VdotH = saturate(dot(V, H));
    vec3 SpecularBRDF(float NdotV, float NdotL, float NdotH, float LdotH, float VdotH, vec3 F0, float roughness) {
        float r   = roughness * roughness;
    
        float D = Distribution(NdotH, r);
        float V = Geometry(NdotV, NdotL, r);
        vec3  F = Fresnel(VdotH, F0);
    
        return D * V * F / (4 * NdotV * NdotL);
    }
    

    注意,这里的NdotLNdotV等都是clamp到0到1的。

    接下来具体看 f r f_r fr中每个分量的可能形式都有哪些。

    法向分布函数 D

    这一部分介绍3个法向分布函数的公式,以及一个推广。

    Beckmann

    来源[5], D B e c k m a n n D_{Beckmann} DBeckmann假设微表面的法向分布是以 n {\bf{n}} n为均值的高斯分布,也即 h {\bf{h}} h n {\bf{n}} n越接近,反射的光线越多, D B e c k m a n n D_{Beckmann} DBeckmann越大。再结合 D B e c k m a n n D_{Beckmann} DBeckmann的积分约束,求得:

    D B e c k m a n n ( h , α ) = χ + ( n , h ) π α 2 ( n ⋅ h ) 4 e ( n ⋅ h ) 2 − 1 α 2 ( n ⋅ h ) 2 D_{Beckmann}({\bf{h}}, \alpha) = \frac{\chi^+({\bf{n}},{\bf{h}})}{\pi\alpha^2({\bf{n}}\cdot{\bf{h}})^4}e^{\frac{({\bf{n}}\cdot{\bf{h}})^2-1}{\alpha^2({\bf{n}}\cdot{\bf{h}})^2}} DBeckmann(h,α)=πα2(nh)4χ+(n,h)eα2(nh)2(nh)21

    这个公式有一个很致命的缺陷,那就是当roughness接近于1的时候, D B e c k m a n n D_{Beckmann} DBeckmann n ⋅ h ∈ [ 0 , 1 ] {\bf{n}}\cdot{\bf{h}}\in[0,1] nh[0,1]区间内,不是单调递增的。表现在渲染上,就是在高光最强的中心点会产生一个暗斑。

    DistributionBeckmann_graph

    shader实现:

    float DistributionBeckmann(float NdotH, float r) {
        float NdotH2   = NdotH * NdotH;
        float r2       = r * r;
        float r2NdotH2 = r2 * NdotH2;
        return exp((NdotH2 - 1) / (r2NdotH2)) / (PI * r2NdotH2 * NdotH2);
    }
    

    BlinnPhong

    来源[6],BlinnPhong公式纯粹是基于经验的,在恰当选取参数的情况下,它的函数曲线非常接近于Beckmann。

    BlinnPhong原始的模型是:

    D B l i n n ( h , α ) = χ + ( n , h ) α p + 2 2 π ( n ⋅ h ) α p D_{Blinn}({\bf{h}}, \alpha) = \chi^+({\bf{n}},{\bf{h}})\frac{\alpha_p + 2}{2\pi}({\bf{n}}\cdot{\bf{h}})^{\alpha_p} DBlinn(h,α)=χ+(n,h)2παp+2(nh)αp

    其中, α p \alpha_p αp表示粗糙系数,或者准确的说,是光滑系数—— α p \alpha_p αp越大,表示物体表面越光滑。

    • α p = ∞ \alpha_p=\infty αp=的时候,表示绝对光滑的物体,此时 D B l i n n ( h , α ) D_{Blinn}({\bf{h}}, \alpha) DBlinn(h,α)只有在 h = n {\bf{h}} = {\bf{n}} h=n,即入射角等于出射角的时候为 ∞ \infty ,否则为0。
    • α p = 0 \alpha_p=0 αp=0的时候,表示绝对粗糙的物体, D B l i n n ( h , α ) = 1 π D_{Blinn}({\bf{h}}, \alpha) = \frac{1}{\pi} DBlinn(h,α)=π1,这个式子也是后面会提到的diffuse的式子。

    α p = ( 2 α 2 − 2 ) \alpha_p = (\frac{2}{\alpha^2} - 2) αp=(α222),则有:

    D B l i n n ( h , α ) = χ + ( n , h ) π α 2 ( n ⋅ h ) ( 2 α 2 − 2 ) D_{Blinn}({\bf{h}}, \alpha) = \frac{\chi^+({\bf{n}},{\bf{h}})}{\pi\alpha^2}({\bf{n}}\cdot{\bf{h}})^{(\frac{2}{\alpha^2} - 2)} DBlinn(h,α)=πα2χ+(n,h)(nh)(α222)

    这个公式即接近于Beckmann的法向分布公式,也是常用的BlinnPhong形式。

    shader实现:

    float DistributionBlinnPhong(float NdotH, float r) {
        float a = r * r;
        return pow(NdotH, 2.0 / a - 2.0) / (PI * a);
    }
    

    GGX

    来源[3],GGX是根据实测数据拟合出来的一个公式:

    D G G X ( h , α ) = χ + ( n , h ) ⋅ α 2 π ( ( n ⋅ h ) 2 ( α 2 − 1 ) + 1 ) 2 D_{GGX}({\bf{h}}, \alpha) = \frac{\chi^+({\bf{n}},{\bf{h}})\cdot\alpha^2}{\pi(({\bf{n}}\cdot{\bf{h}})^2(\alpha^2-1)+1)^2} DGGX(h,α)=π((nh)2(α21)+1)2χ+(n,h)α2

    shader实现:

    float DistributionGGX(float NdotH, float r) {
        float a2     = r * r;
        float NdotH2 = NdotH * NdotH;
        float nom    = a2;
        float denom  = (NdotH2 * (a2 - 1.0) + 1.0);
        denom        = PI * denom * denom;
        return nom / max(denom, 0.001);
    }
    

    除了这三种公式,还有更多更复杂的法向分布函数D,具体可以参考[17]。但是其实最常用的还是GGX(及其各向异性模式),无论是游戏还是影视行业都比较喜欢用GGX。

    GTR

    Burley通过对Berry(与GGX公式类似,分母上的指数为1)和GGX公式的观察,提出了广义的Trowbridge-Reitz(Generalized-Trowbridge-Reitz,GTR)法线分布函数:

    D G T R ( h , α ) = c ⋅ χ + ( n , h ) ( ( n ⋅ h ) 2 ( α 2 − 1 ) + 1 ) γ D_{GTR}({\bf{h}}, \alpha) = \frac{c\cdot\chi^+({\bf{n}},{\bf{h}})}{(({\bf{n}}\cdot{\bf{h}})^2(\alpha^2-1)+1)^\gamma} DGTR(h,α)=((nh)2(α21)+1)γcχ+(n,h)

    其中, c c c表示缩放系数,是一个常数; γ \gamma γ用于控制尾部的形状,当 γ = 1 \gamma=1 γ=1的时候, D G T R D_{GTR} DGTR就是Berry公式,当 γ = 2 \gamma=2 γ=2的时候, D G T R D_{GTR} DGTR就是 D G G X D_{GGX} DGGX

    γ \gamma γ的取值对 D G T R D_{GTR} DGTR的影响如下图所示。

    GTR

    以下是 γ = 1 \gamma=1 γ=1 γ = 2 \gamma=2 γ=2时的shader实现:

    float DistributionGTR1(float NdotH, float r)
    {
        if (r >= 1) return 1/PI;
        float a2 = r*r;
        float t = 1 + (a2-1)*NdotH*NdotH;
        return (a2-1) / (PI*log(a2)*t);
    }
    
    float DistributionGTR2(float NdotH, float r)
    {
        float a2 = r*r;
        float t = 1 + (a2-1)*NdotH*NdotH;
        return a2 / (PI * t * t);
    }
    

    效果对比

    Distribution

    可以看出,BlinnPhong和Beckmann的差异不大。而GGX有着更平滑的边缘和更小的峰值。除此之外,GGX运算压力更小,因为它没有指数操作。

    遮挡项 G

    和法向分布函数 D D D一样,遮挡项 G G G也是入射角、出射角和表面粗糙度的函数。

    有些文章会把遮挡项G和BRDF的分母 ( n ⋅ l ) ( n ⋅ v ) ({\bf{n}}\cdot{\bf{l}})({\bf{n}}\cdot{\bf{v}}) (nl)(nv)放在一起组成一项约分掉,这也是一种优化思路,因为G通常包含这两个cosine因子。这里我们约定本文的遮挡项 G G G是不约分 ( n ⋅ l ) ( n ⋅ v ) ({\bf{n}}\cdot{\bf{l}})({\bf{n}}\cdot{\bf{v}}) (nl)(nv) G G G

    Implicit

    来源[7],有些BRDF公式会忽略遮挡项G,将其跟分母上的 ( n ⋅ l ) ( n ⋅ v ) ({\bf{n}}\cdot{\bf{l}})({\bf{n}}\cdot{\bf{v}}) (nl)(nv)一起忽略掉,这就有了第一个隐式 G G G

    G I m p l i c i t ( l , v , h ) = ( n ⋅ l ) ( n ⋅ v ) G_{Implicit}({\bf{l}},{\bf{v}},{\bf{h}})=({\bf{n}}\cdot{\bf{l}})({\bf{n}}\cdot{\bf{v}}) GImplicit(l,v,h)=(nl)(nv)

    它的形态大概是,当且仅当视角和光源都垂直于物体表面的时候, G I m p l i c i t = 1 G_{Implicit}=1 GImplicit=1,光源、视角与物体表面法线的夹角越大, G I m p l i c i t G_{Implicit} GImplicit越小,直到衰减为0,这也是很符合常识的。

    shader实现:

    float GeometryImplicit(float NdotV, float NdotL) {
        return NdotL * NdotV;
    }
    

    但是隐式遮挡项 G I m p l i c i t G_{Implicit} GImplicit最大的问题在于,它随着视角的衰减速度太快,这会使得高光区域太窄。为了解决这个问题,我们继续看显式的遮挡项 G G G

    Cook-Torrance

    来源[9], G C o o k − T o r r a n c e G_{Cook-Torrance} GCookTorrance解决了 G I m p l i c i t G_{Implicit} GImplicit衰减速度太快的问题:

    G C o o k − T o r r a n c e ( l , v , h ) = min ⁡ ( 1 , 2 ( n ⋅ h ) ( n ⋅ v ) v ⋅ h , 2 ( n ⋅ h ) ( n ⋅ l ) v ⋅ h ) G_{Cook-Torrance}({\bf{l}},{\bf{v}},{\bf{h}})=\min{\left(1, \frac{2({\bf{n}}\cdot{\bf{h}})({\bf{n}}\cdot{\bf{v}})}{{\bf{v}}\cdot{\bf{h}}}, \frac{2({\bf{n}}\cdot{\bf{h}})({\bf{n}}\cdot{\bf{l}})}{{\bf{v}}\cdot{\bf{h}}}\right)} GCookTorrance(l,v,h)=min(1,vh2(nh)(nv),vh2(nh)(nl))

    shader实现:

    float GeometryCookTorrance(float NdotV, float NdotL, float VdotH, float NdotH) {
        float ct1 = 2 * NdotH * NdotV / VdotH;
        float ct2 = 2 * NdotH * NdotL / VdotH;
        return min(1, min(ct1, ct2));
    }
    

    Kelemen

    来源[10],也是解决 G I m p l i c i t G_{Implicit} GImplicit衰减速度太快的问题,同时 G K e l e m e n G_{Kelemen} GKelemen G C o o k − T o r r a n c e G_{Cook-Torrance} GCookTorrance的效率更高:

    G K e l e m e n ( l , v , h ) = ( n ⋅ l ) ( n ⋅ v ) ( v ⋅ h ) 2 G_{Kelemen}({\bf{l}},{\bf{v}},{\bf{h}})=\frac{({\bf{n}}\cdot{\bf{l}})({\bf{n}}\cdot{\bf{v}})}{({\bf{v}}\cdot{\bf{h}})^2} GKelemen(l,v,h)=(vh)2(nl)(nv)

    shader实现:

    float GeometryKelemen(float NdotV, float NdotL, float VdotH) {
        return NdotV * NdotL / (VdotH * VdotH);
    }
    

    Neumann

    来源[8], G N e u m a n n G_{Neumann} GNeumann用另一种方式解决了 G I m p l i c i t G_{Implicit} GImplicit衰减速度太快的问题:

    G N e u m a n n ( l , v , h ) = ( n ⋅ l ) ( n ⋅ v ) max ⁡ ( n ⋅ l , n ⋅ v ) G_{Neumann}({\bf{l}},{\bf{v}},{\bf{h}})=\frac{({\bf{n}}\cdot{\bf{l}})({\bf{n}}\cdot{\bf{v}})}{\max{({\bf{n}}\cdot{\bf{l}}, {\bf{n}}\cdot{\bf{v}}})} GNeumann(l,v,h)=max(nl,nv)(nl)(nv)

    shader实现:

    float GeometryNeumann(float NdotV, float NdotL) {
        return (NdotL * NdotV) / max(NdotL, NdotV);
    }
    

    但是,以上三个解决方案也不够完美。前面提到过,遮挡项G应该是入射角、出射角和表面粗糙度的函数,而以上四个G,包括隐式遮挡项都与粗糙度无关。

    Smith

    Smith家族[13]都是采用了前面介绍的 G 1 G_1 G1相乘的形式:

    G 2 ( l , v , h ) = G 1 ( l ) G 1 ( v ) G_2({\bf{l}},{\bf{v}},{\bf{h}})=G_1({\bf{l}})G_1({\bf{v}}) G2(l,v,h)=G1(l)G1(v)

    他们之间的区别就是 G 1 G_1 G1的选取不同。

    Beckmann

    Beckmann的 G G G是跟 D D D一起提出的,前面介绍过 G G G是可以从 D D D推导出来的,因此Beckmann的 Λ \Lambda Λ为:

    c = n ⋅ v α 1 − ( n ⋅ v ) 2 Λ ( v ) = erf ( c ) − 1 2 + 1 2 c π exp ⁡ ( − c 2 ) \begin{aligned} c & = \frac{{\bf{n}}\cdot{\bf{v}}}{\alpha\sqrt{1-({\bf{n}}\cdot{\bf{v}})^2}} \\ \Lambda({\bf{v}}) & = \frac{\text{erf}(c)-1}{2}+\frac{1}{2c\sqrt{\pi}}\exp(-c^2) \end{aligned} cΛ(v)=α1(nv)2 nv=2erf(c)1+2cπ 1exp(c2)

    但是由于有 erf \text{erf} erf函数的存在,计算起来过于复杂,因此通常用如下的近似形式:

    Λ ( v ) ≈ { 1 − 1.259 x + 0.396 c 2 3.535 c + 2.181 c 2 , if  c < 1.6 0 , if  c ≥ 1.6 \Lambda({\bf{v}}) \approx \begin{cases} \frac{1-1.259x+0.396c^2}{3.535c+2.181c^2}, & \text{if }c<1.6 \\ 0, & \text{if }c\geq1.6 \end{cases} Λ(v){3.535c+2.181c211.259x+0.396c2,0,if c<1.6if c1.6

    因此,Beckmann的 G 1 G_1 G1

    G B e c k m a n n ( v ) ≈ { 3.535 c + 2.181 c 2 1 + 2.276 c + 2.577 c 2 , if  c < 1.6 1 , if  c ≥ 1.6 G_{Beckmann}({\bf{v}}) \approx \begin{cases} \frac{3.535c+2.181c^2}{1+2.276c+2.577c^2}, & \text{if }c<1.6 \\ 1, & \text{if }c\geq1.6 \end{cases} GBeckmann(v){1+2.276c+2.577c23.535c+2.181c2,1,if c<1.6if c1.6

    shader实现:

    float GeometryBeckmann(float NdotV, float r) {
        float c  = NdotV / (r * sqrt(1 - NdotV * NdotV));
        float c2 = c * c;
        if (c < 1.6)
            return (3.535 * c + 2.181 * c2) / (1 + 2.276 * c + 2.577 * c2);
        else
            return 1.0;
    }
    float GeometrySmithBeckmann(float NdotV, float NdotL, float r) {
        float ggx2 = GeometryBeckmann(NdotV, r);
        float ggx1 = GeometryBeckmann(NdotL, r);
        return ggx1 * ggx2;
    }
    
    GGX

    GGX[3]跟Beckmann类似,都是从法向分布函数推导出来的:

    c = n ⋅ v α 1 − ( n ⋅ v ) 2 Λ ( v ) = − 1 + 1 + 1 c 2 2 \begin{aligned} c & = \frac{{\bf{n}}\cdot{\bf{v}}}{\alpha\sqrt{1-({\bf{n}}\cdot{\bf{v}})^2}} \\ \Lambda({\bf{v}}) & = \frac{-1+\sqrt{1+\frac{1}{c^2}}}{2} \end{aligned} cΛ(v)=α1(nv)2 nv=21+1+c21

    对应的 G 1 G_1 G1定义为

    G G G X ( v ) = 2 ( n ⋅ v ) ( n ⋅ v ) + α 2 + ( 1 − α 2 ) ( n ⋅ v ) 2 G_{GGX}({\bf{v}}) = \frac{2({\bf{n}}\cdot{\bf{v}})}{({\bf{n}}\cdot{\bf{v}})+\sqrt{\alpha^2+(1-\alpha^2)({\bf{n}}\cdot{\bf{v}})^2}} GGGX(v)=(nv)+α2+(1α2)(nv)2 2(nv)

    shader实现:

    float GeometryGGX(float NdotV, float r) {
        float r2 = r * r;
        return (2 * NdotV) / (NdotV + sqrt(r2 + (1 - r2) * NdotV * NdotV));
    }
    float GeometrySmithGGX(float NdotV, float NdotL, float r) {
        float ggx2 = GeometryGGX(NdotV, r);
        float ggx1 = GeometryGGX(NdotL, r);
        return ggx1 * ggx2;
    }
    
    GGX Joint

    前面提到的GGX用的是 G 2 = G 1 ∗ G 1 G_2=G_1*G_1 G2=G1G1的separable G,如果用height-correlated G,那么 G 2 G_2 G2变为:

    G 2 − G G X J o i n t ( l , v , m ) = 1 1 + Λ ( l ) + Λ ( v ) = 2 ( n ⋅ v ) ( n ⋅ l ) ( n ⋅ l ) ⋅ α 2 + ( 1 − α 2 ) ( n ⋅ v ) 2 + ( n ⋅ v ) ⋅ α 2 + ( 1 − α 2 ) ( n ⋅ l ) 2 \begin{aligned} G_{2-GGXJoint}({\bf{l}},{\bf{v}},{\bf{m}}) & =\frac{1}{1+\Lambda({\bf{l}})+\Lambda({\bf{v}})}\\ & =\frac{2({\bf{n}}\cdot{\bf{v}})({\bf{n}}\cdot{\bf{l}})}{({\bf{n}}\cdot{\bf{l}})\cdot\sqrt{\alpha^2+(1-\alpha^2)({\bf{n}}\cdot{\bf{v}})^2} + ({\bf{n}}\cdot{\bf{v}})\cdot\sqrt{\alpha^2+(1-\alpha^2)({\bf{n}}\cdot{\bf{l}})^2}} \end{aligned} G2GGXJoint(l,v,m)=1+Λ(l)+Λ(v)1=(nl)α2+(1α2)(nv)2 +(nv)α2+(1α2)(nl)2 2(nv)(nl)

    shader实现:

    float GeometrySmithGGXJoint(float NdotV, float NdotL, float r) {
        float r2 = r * r;
        float Vis_SmithV = NdotL * sqrt(NdotV * (NdotV - NdotV * r2) + r2);
    	float Vis_SmithL = NdotV * sqrt(NdotL * (NdotL - NdotL * r2) + r2);
    	return 2 * NdotV * NdotL / (Vis_SmithV + Vis_SmithL);
    }
    

    为了提高计算效率,UE4对GGX Joint方法做了一个近似,公式为:

    G 2 − G G X J o i n t ( l , v , m ) = 1 1 + Λ ( l ) + Λ ( v ) ≈ 2 ( n ⋅ v ) ( n ⋅ l ) ( n ⋅ l ) ⋅ ( α + ( 1 − α ) ( n ⋅ v ) ) + ( n ⋅ v ) ⋅ ( α + ( 1 − α ) ( n ⋅ l ) ) \begin{aligned} G_{2-GGXJoint}({\bf{l}},{\bf{v}},{\bf{m}}) & =\frac{1}{1+\Lambda({\bf{l}})+\Lambda({\bf{v}})}\\ & \approx\frac{2({\bf{n}}\cdot{\bf{v}})({\bf{n}}\cdot{\bf{l}})}{({\bf{n}}\cdot{\bf{l}})\cdot(\alpha+(1-\alpha)({\bf{n}}\cdot{\bf{v}})) + ({\bf{n}}\cdot{\bf{v}})\cdot(\alpha+(1-\alpha)({\bf{n}}\cdot{\bf{l}}))} \end{aligned} G2GGXJoint(l,v,m)=1+Λ(l)+Λ(v)1(nl)(α+(1α)(nv))+(nv)(α+(1α)(nl))2(nv)(nl)

    shader实现:

    float GeometryGGXJointApprox(float NdotV, float NdotL, float r) {
        return (NdotV) / (NdotL * (r + (1 - r) * NdotV));
    }
    float GeometrySmithGGXJointApprox(float NdotV, float NdotL, float r) {
    	float Vis_SmithV = NdotL * ( NdotV * ( 1 - r ) + r );
    	float Vis_SmithL = NdotV * ( NdotL * ( 1 - r ) + r );
    	return 2 * NdotV * NdotL / ( Vis_SmithV + Vis_SmithL );
    }
    
    Schlick-Beckmann

    Schlick[11]的 G 1 G_1 G1定义为
    k = α 2 π G S c h l i c k ( v ) = n ⋅ v ( n ⋅ v ) ( 1 − k ) + k k=\alpha\sqrt{\frac{2}{\pi}} \\ G_{Schlick}({\bf{v}})=\frac{{\bf{n}}\cdot{\bf{v}}}{({\bf{n}}\cdot{\bf{v}})(1-k)+k} k=απ2 GSchlick(v)=(nv)(1k)+knv
    shader实现:

    float GeometrySchlickBeckmann(float NdotV, float r) {
        float k     = (r)*sqrt(2.0 / PI);
        float nom   = NdotV;
        float denom = NdotV * (1.0 - k) + k;
        return nom / denom;
    }
    float GeometrySmithSchlickBeckmann(float NdotV, float NdotL, float r) {
        float ggx2 = GeometrySchlickBeckmann(NdotV, r);
        float ggx1 = GeometrySchlickBeckmann(NdotL, r);
        return ggx1 * ggx2;
    }
    
    Schlick-GGX

    Schlick-GGX[12]曾经是UE4所采用的的一个模型,跟Schlick有些类似, G 1 G_1 G1定义为
    k = α 2 G S c h l i c k ( v ) = n ⋅ v ( n ⋅ v ) ( 1 − k ) + k k=\frac{\alpha}{2} \\ G_{Schlick}({\bf{v}})=\frac{{\bf{n}}\cdot{\bf{v}}}{({\bf{n}}\cdot{\bf{v}})(1-k)+k} k=2αGSchlick(v)=(nv)(1k)+knv
    shader实现:

    float GeometrySchlickGGX(float NdotV, float r) {
        float k     = r * 0.5;
        float nom   = NdotV;
        float denom = NdotV * (1.0 - k) + k;
        return nom / denom;
    }
    float GeometrySmithSchlickGGX(float NdotV, float NdotL, float r) {
        float ggx2 = GeometrySchlickGGX(NdotV, r);
        float ggx1 = GeometrySchlickGGX(NdotL, r);
        return ggx1 * ggx2;
    }
    

    这里面还有一个细节,那就是迪士尼后来提出了对粗糙粗roughness做一个remapping,使得它更接近于真实:
    α ′ = ( r o u g h n e s s + 1 2 ) 2 \alpha' = (\frac{roughness + 1}{2})^2 \\ α=(2roughness+1)2
    其他的部分不变。这样shader实现为:

    float GeometrySmithSchlickGGX(float NdotV, float NdotL, float roughness) {
        float r = (roughness + 1.0) * 0.5; // remapping roughness
        r = r * r
        float ggx2 = GeometrySchlickGGX(NdotV, r);
        float ggx1 = GeometrySchlickGGX(NdotL, r);
        return ggx1 * ggx2;
    }
    

    注意,此时GeometrySmithSchlickGGX的输入参数不是r,而改为了roughness

    优化

    考虑到几乎所有 G G G都带有 ( n ⋅ v ) ( n ⋅ l ) ({\bf{n}}\cdot{\bf{v}})({\bf{n}}\cdot{\bf{l}}) (nv)(nl)项,可以跟 f r ( l , v ) f_r({\bf{l}},{\bf{v}}) fr(l,v)的分母约分,因此在实现时,可以考虑定义
    G ′ = G ( n ⋅ v ) ( n ⋅ l ) G'=\frac{G}{({\bf{n}}\cdot{\bf{v}})({\bf{n}}\cdot{\bf{l}})} G=(nv)(nl)G
    节省一部分计算。

    这样做不只是出于性能的考虑,也是出于精度的考虑。如果 n ⋅ v {\bf{n}}\cdot{\bf{v}} nv n ⋅ l {\bf{n}}\cdot{\bf{l}} nl的乘积接近于0,那么specular项的分母会非常小,严重影响其精度,极端的情况下会在过渡区域产生一道割裂的分界线。下图展示了 ( n ⋅ v ) ∗ ( n ⋅ l ) ∗ 10 ({\bf{n}}\cdot{\bf{v}})*({\bf{n}}\cdot{\bf{l}})*10 (nv)(nl)10(左)、未优化时(中)、优化后(右)的效果,可以看出左侧两张图的分界线非常吻合。优化后则没有颜色割裂的问题。

    Geometry-bugs

    效果对比

    roughness = 0.9,计算球体的遮挡项效果为:

    Geometry

    最后一排是常用的几种方法,差异并不大,边缘的过渡也比较好。

    菲涅尔项 F

    菲涅尔项描述的是物体表面的反射、折射现象。一般我们会采用常量 F 0 F_0 F0来计算菲涅尔项 F ( v , h , F 0 ) F({\bf{v}},{\bf{h}},F_0) F(v,h,F0)

    要说明白菲涅尔项,得从光学在介质表面的折射反射现象说起。我们知道光线会在介质表面产生不连续性,具体表现为一部分光线反射——遵循光线反射定律,入射角等于反射角;另一部分光线会折射进入介质——遵循光线折射定律,折射角取决于入射角的大小以及构成交界面的两种材质,即斯涅耳定律(Snell’s law):

    n 1 sin ⁡ ( θ i ) = n 2 sin ⁡ ( θ t ) n_1\sin(\theta_i)=n_2\sin(\theta_t) n1sin(θi)=n2sin(θt)

    斯涅耳定律描述的仅仅是光线的角度,但是图形学研究的其实是光线的radiance/irradiance,所以我们要更进一步。定义Fresnel reflectance R F R_F RF为反射光线的radiance占入射光线radiance的比例, R F R_F RF是入射角 θ i \theta_i θi的函数。那么对于入射光线 L i L_i Li,在角度 θ i \theta_i θi时反射光线的radiance为 R F ( θ i ) L i R_F(\theta_i)L_i RF(θi)Li。再考虑折射部分,根据能量守恒,没有反射的能量都会被折射(不考虑被吸收的能量),因此折射的flux占入射flux的比例是 1 − R F 1-R_F 1RF。这里需要强调的是,radiance定义的是“irradiance每立体角”,它的大小跟角度有关系,因此折射光线的radiance L t L_t Lt不能简单用 1 − R F 1-R_F 1RF乘上 L i L_i Li,而要转换角度:
    L t = ( 1 − R F ( θ i ) ) sin ⁡ 2 θ i sin ⁡ 2 θ t L i L_t = (1-R_F(\theta_i))\frac{\sin^2\theta_i}{\sin^2\theta_t}L_i Lt=(1RF(θi))sin2θtsin2θiLi
    将斯涅耳定律带入上式,得到:
    L t = ( 1 − R F ( θ i ) ) n 2 2 n 1 2 L i L_t = (1-R_F(\theta_i))\frac{n_2^2}{n_1^2}L_i Lt=(1RF(θi))n12n22Li
    介绍了这么多 R F R_F RF的相关知识,其实关键点还是前面说的, R F R_F RF是入射角 θ i \theta_i θi的函数。我们再回头考虑这个 R F R_F RF与入射角 θ i \theta_i θi的关系。当 θ i = 90 ° \theta_i=90\degree θi=90°的时候,即 R F ( 90 ° ) R_F(90\degree) RF(90°),此时入射光平行于平面,垂直于法向,不存在折射光线, R F ( 90 ° ) = 1 R_F(90\degree)=1 RF(90°)=1;当 θ i = 0 ° \theta_i=0\degree θi=0°的时候,即 R F ( 0 ° ) R_F(0\degree) RF(0°),此时反射光线占比最低,根据不同的材质这个 R F ( 0 ° ) R_F(0\degree) RF(0°)有不同的值,Real-time Rendering[14]给出了常见的材质的 R F R_F RF θ i \theta_i θi的关系曲线:

    Fresnel

    为了近似这个曲线,采取的策略是利用 R F ( 0 ° ) R_F(0\degree) RF(0°),也就是前面说的 F 0 F_0 F0
    R F ( θ i ) ≈ R F ( 0 ° ) + ( 1 − R F ( 0 ° ) ) ( 1 − cos ⁡ θ i ) 5 R_F(\theta_i)\approx R_F(0\degree) + (1-R_F(0\degree))(1-\cos\theta_i)^5 RF(θi)RF(0°)+(1RF(0°))(1cosθi)5
    这里有一个默认的假设是, R F ( 90 ° ) = 1 R_F(90\degree)=1 RF(90°)=1,如果, R F ( 90 ° ) R_F(90\degree) RF(90°)未知, R F ( θ i ) R_F(\theta_i) RF(θi)应该写为:
    R F ( θ i ) ≈ R F ( 0 ° ) + ( R F ( 90 ° ) − R F ( 0 ° ) ) ( 1 − cos ⁡ θ i ) 5 R_F(\theta_i)\approx R_F(0\degree) + (R_F(90\degree)-R_F(0\degree))(1-\cos\theta_i)^5 RF(θi)RF(0°)+(RF(90°)RF(0°))(1cosθi)5
    这个 R F ( 90 ° ) R_F(90\degree) RF(90°)也就是 F 90 F_{90} F90

    最后,我们看一下 F 0 F_0 F0怎么计算。对于dielectrics来说, F 0 F_0 F0的值取决于折射率,公式为:
    F 0 = 0.16 ⋅ r e f l e c t a n c e 2 F_0=0.16\cdot reflectance^2 F0=0.16reflectance2
    其中, r e f l e c t a n c e reflectance reflectance由物体表面的材质定义。

    对于conductors, F 0 F_0 F0通过金属度metallic和basecolor来计算:
    F 0 = b a s e C o l o r ⋅ m e t a l l i c F_0=baseColor\cdot metallic F0=baseColormetallic
    综合dielectrics和dielectric,得到:

        vec3 F0 = 0.16 * reflectance * reflectance * (1.0 - metallic) + baseColor.xyz * metallic;
    

    说明白了 F 0 F_0 F0,我们接下来看看菲涅尔函数 F F F有哪些形式。

    简单形式

    最简单的情况,直接令菲涅尔函数等于 F 0 F_0 F0
    F N o n e ( v , h ) = F 0 F_{None}({\bf{v}},{\bf{h}})=F_0 FNone(v,h)=F0
    shader实现:

    vec3 Fresnel(vec3 F0) {
        return F0;
    }
    

    Schlick

    来源[11],公式:
    F S c h l i c k ( v , h ) = F 0 + ( 1 − F 0 ) ( 1 − ( v ⋅ h ) ) 5 F_{Schlick}({\bf{v}},{\bf{h}})=F_0+(1-F_0)(1-({\bf{v}}\cdot{\bf{h}}))^5 FSchlick(v,h)=F0+(1F0)(1(vh))5
    也就是我们前面说到的对 R F R_F RF的拟合。shader实现:

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

    如果引入 F 90 F_{90} F90,则变成:
    F S c h l i c k ( v , h ) = F 0 + ( F 90 − F 0 ) ( 1 − ( v ⋅ h ) ) 5 F_{Schlick}({\bf{v}},{\bf{h}})=F_0+(F_{90}-F_0)(1-({\bf{v}}\cdot{\bf{h}}))^5 FSchlick(v,h)=F0+(F90F0)(1(vh))5
    shader实现:

    vec3 FresnelSchlick(float VdotH, vec3 F0, vec F90) {
        return F0 + (F90 - F0) * pow(1.0 - VdotH, 5.0);
    }
    

    对specular来说, F 90 F_{90} F90可以从 F 0 F_0 F0计算得来[1]:

        float F90 = saturate(dot(F0, vec3(50.0 * 0.33)));
    

    Cook-Torrance

    来源[9],公式:
    η = 1 + F 0 1 − F 0 c = v ⋅ h g = η 2 + c 2 − 1 F C o o k − T o r r a n c e ( v , h ) = 1 2 ( g − c g + c ) 2 ( 1 + ( ( g + c ) c − 1 ( g − c ) c + 1 ) 2 ) \begin{aligned} \eta & =\frac{1+\sqrt{F_0}}{1-\sqrt{F_0}} \\ c & = {\bf{v}}\cdot{\bf{h}} \\ g & = \sqrt{\eta^2+c^2-1} \\ F_{Cook-Torrance}({\bf{v}},{\bf{h}}) & =\frac{1}{2}\left(\frac{g-c}{g+c}\right)^2\left(1+\left(\frac{(g+c)c-1}{(g-c)c+1}\right)^2\right) \end{aligned} ηcgFCookTorrance(v,h)=1F0 1+F0 =vh=η2+c21 =21(g+cgc)2(1+((gc)c+1(g+c)c1)2)
    shader实现:

    float FresnelCookTorrance(float VdotH, float F0) {
        float sqrtF = sqrt(F0);
        float Eta   = (1.0 + sqrtF) / (1.0 - sqrtF);
        float g     = sqrt(Eta * Eta + VdotH * VdotH - 1.0);
        return 0.5 * pow((g - VdotH) / (g + VdotH), 2) *
               (1 + pow(((g + VdotH) * VdotH - 1.0) / ((g - VdotH) * VdotH + 1.0), 2));
    }
    

    Diffuse BRDF

    相比于繁琐的specular部分,diffuse部分就简单的多。diffuse部分由baseColor和diffuse系数相乘得到,即:
    L d ( v ) = c d i f f ⋅ f d L_d({\bf{v}})={\bf{c}}_{diff}\cdot f_d Ld(v)=cdifffd
    shader实现:

        vec3 colorDiffuse = baseColor * DiffuseBRDF(NdotV, NdotL, LdotH, roughness);
    

    接下来看一下 f d f_d fd的可能取值。

    Lambert

    Lambert模型认为既然diffuse是漫反射,不如简单地认为各个方向都是一样的值,即出射光线的radiance与入射光线的角度无关。

    f d = 1 π f_d = \frac{1}{\pi} fd=π1

    这个实现相当于,采用BlinnPhong的法向分布 D B l i n n ( h , α ) = 1 π D_{Blinn}({\bf{h}}, \alpha) = \frac{1}{\pi} DBlinn(h,α)=π1,同时令遮挡项为隐式形式,并且菲涅尔项为1。虽然简单,但是已经足够近似现实了,效果还不错。

    shader实现:

    float DiffuseLambert() {
        return 1.0 / PI;
    }
    

    虽然Lambert模型已经足够接近真实情况,但是它还是不够理想。我们前面提到过,diffuse分量本质上是光线折射进入物体表面,经过多次反射再折射出来的现象,也就是它不是物理上真实存在的一个光学现象。而在讨论specular菲涅尔项的时候又提到过,反射部分会随着入射光线的角度变化,那么折射部分相应的也会随着入射角度变化,既然如此,来自于折射部分的diffuse分量肯定也是会随着入射光线的角度而改变的!也就是说, f d f_d fd是入射角 l {\bf{l}} l的函数: f d ( l ) f_d({\bf{l}}) fd(l)

    同时, f d f_d fd也应该是出射角 v {\bf{v}} v的函数[14]: f d ( l , v ) f_d({\bf{l}},{\bf{v}}) fd(l,v)。因为菲涅尔项考虑的是镜面反射,入射角等于出射角,而diffuse项的入射角不一定等于出射角,因此两个角度都会影响 f d f_d fd

    再者,前面影响specular分量的参数当中, r o u g h n e s s roughness roughness也会影响 f d f_d fd。根据常识,不同粗糙程度的物体的diffuse是有明显的不同的。即