Shadertoy--光线追踪与Nayar反射模型

前言

Shadertoy是一个神奇的网站，有无数图形学大神分享的用shader写的各种让人匪夷所思的效果，而代码仅仅只有数十行至数百行。虽然很多实现无法直接应用到我们开发的实时渲染上，但分析它们的实现可以让我们理解shader以及带来很多灵感。

OrenNayar是一种微表面反射模型，本文主要分析shadertoy中实现的一个简易的OrenNayar模型的实现。要展示OrenNayar模型的效果，必须首先要在着色器里面实现一个简单的光线追踪器（其实我主要是想分析一下作者实现的光线追踪器....）

当然分析过程纯属个人猜测，不代表原作者思路，原文地址如下：
https://www.shadertoy.com/view/ldBGz3

正文

Lambert和OrenNayar模型的对比

Lambert和OrenNayar都是一种光照反射模型，即描述特定材质的物体表面如何反射光线。（具体原理可以看《Physically Based Rendering》第五章和第八章相关内容）。
在Lambert模型中，所有方向的反射光的辐射率都相等，辐射率大小与入射光与表面法线的夹角的余弦值成正比。OrenNayar是一种微表面反射模型，即表面是由无数个微小的表面构成，每个微表面被假设为是理想Lambert模型，OrenNayar模型中考虑了微表面之间的遮挡，阴影，法线方向等因素。对OrenNayar的原理感兴趣的同学可以看《Oren_SIGGRAPH94》文档（反正我是没怎么看懂）。该文档的开头引入了一张图片

这张图片主要是展示理想Lambert模型的不足，左边为真实的圆柱体花瓶，右边是用理想Lambert模型渲染的圆柱体花瓶。可以看到，Lambert模型中，当靠近花瓶两侧边缘时，亮度显著下降，而真实的花瓶的亮度变化整体上基本是趋于平坦的。而用OrenNayar渲染的模型可以更接近真实。

可以看下shadertoy中渲染出来的两个球

可以看出，Lambert模型随着角度变化亮度会快速下降，而OrenNayar模型亮度变化趋于平缓。

光线追踪原理

简单画了个图，简单概况下光线追踪原理。我们之所以能看见物体，是因为光源发射的光经过物体表面的反射，最终有一部分进入了我们的眼睛。光线追踪为了模拟现象，从眼睛（或相机）位置出发，向屏幕上的每个像素上发射一条或多条射线，屏幕像素的颜色和亮度，取决于射线进入场景后与物体交点的反射辐射率（如果交点）或者背景颜色（如果没有交点）。

光线追踪的具体实现由许多方式（比如path tracing），大部分是采用递归的方式，即光线碰到物体表面后，计算反射方向后继续前进，直到遇到光源或递归次数达到最大值。本文要分析的光线追踪非常简单，只检测一次碰撞，但用了特殊的技巧来模拟遮挡和阴影，后面会讲到。

下面开始分析代码。

第一步，创建相机空间，放置光源，球体等

首先要明确一点，因为shadertoy中是在像素着色器中编程，像素着色器的代码本身就是对每一个像素进行操作的过程，因此我们不需要像在CPU中实现光追那样手动搭建一个虚拟相机的屏幕（上图的film平面），也不需要循环遍历每个像素（像素着色器在GPU中并行运行了）。

粘贴相关代码：

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
   //第一步，创建相机空间，放置光源，球体等    
   vec2 p = (-iResolution.xy + 2.0*fragCoord.xy) / iResolution.y;  //屏幕坐标

        //1.1 光源
    vec3 lig = normalize( vec3( 0.6, 0.5, 0.4) );   //主光源
    vec3 bac = normalize( vec3(-0.6, 0.0,-0.4) ); //辅助光源
    vec3 bou = normalize( vec3( 0.0,-1.0, 0.0) );  //辅助光源
    
     // 1.2 相机矩阵
    float an = 0.6 - 0.5*iTime;  
    vec3 ro = vec3( 3.5*cos(an), 1.0, 3.5*sin(an) );   //相机位置，随时间移动
    vec3 ta = vec3( 0.0, 1.0, 0.0 );        //相机方向

    // camera matrix
    vec3 ww = normalize( ta - ro );
    vec3 uu = normalize( cross(ww,vec3(0.0,1.0,0.0) ) );
    vec3 vv = normalize( cross(uu,ww));
    // create view ray
    vec3 rd = normalize( p.x*uu + p.y*vv + 1.5*ww );  //(p.x，p.y, 1.5)是屏幕上的像素在相机空间的位置（或者说是相机空间下的射线方向，因为射线的原点是(0,0,0)），乘以lookat矩阵后，rd转换为世界空间的射线

    //1.3 放置球体
    vec4 sph1 = vec4(-1.2,1.0,0.0,1.0);
    vec4 sph2 = vec4( 1.2,1.0,0.0,1.0);

   //…….
}

1.1 光源

在场景中放置的三个光源，lig，bac，bou，它们都是方向光，向量就代表光源的方向。lig是主光源，对效果贡献最大，是产生直接反射和阴影的主要来源。bac和bou主要是为了模拟真实场景的间接反射而加入的，对效果贡献较小。

1.2 相机矩阵

我们需要相机矩阵（或者叫lookat矩阵）来进行相机坐标系和世界坐标系之间的转换。原理不说了，重要的一点是在相机空间中，观察方向（从视点望向目标点）是+z轴。如下图

1.3 放置球体

两个球体，sph1和sph2，其中 sph1.xyz 为球体的在世界空间的坐标，sp1h.w 为球体的半径。

第二步，判断光线是否与场景中物体是否相交，并求出遮挡系数。

上一步计算的向量rd，就是从摄像头出发，并穿过当前像素的射线，我们现在要逐一判断它和场景中的物体是否相交，也就是光线追踪。

1.1 判断射线与平面是否相交

// 1.1 判断射线与平面是否相交   
float h = (0.0-ro.y)/rd.y;
    if( h>0.0 ) 
    { 
           //求出两个球的遮挡影响系数occ
           //.....
        obj = 0.0;   //obj = 0.0标记光线与平面相交了
    }
}

除了两个球体外，场景中还有一个y == 0的平面。用 h = (0.0-ro.y) / rd.y 来判断射线和平面是否相交。如下图，比如左边的，ro.h > 0，rd.y < 0那么h = (0.0-ro.y)/rd.y 的结果h就是大于0，即射线会与平面相交。相反，如右边的ro和rd，计算出来h小于0，不相交。

如果是相交的，那么接下来要干嘛的，我们先看下面这个图：

首先要注意到，rd是从相机发出的射线，它实际上是经物体表面反射进入到相机的光线wo的取反。wo实际上是场景中所有光入射光线wi的总体作用（这就是BRDF）。这些wi，有可能是直接光照（图中的wi1），也有可能是间接光照（图中的wi2）。因为场景中还存在别的物体（两个球），光源有可能被球挡住了（比如图中的light2的一些光线就无法达到pos）。我们可以确定的是，如果直接光被遮挡得越多，那么到达pos的光的总能量肯定越少，导致wo强度越小。那么如何衡量遮挡的影响呢？作者用了一个方法来估计球体遮挡的影响，那就是oSphere函数。

float oSphere( in vec3 pos, in vec3 nor, in vec4 sph )
{
    vec3 di = sph.xyz - pos;
    float l = length(di);
    return 1.0 - max(0.0, dot(nor, di / l)) * sph.w * sph.w / (l * l); 
}

如下图，di的长度l表示球与交点pos的距离，shp.w表示球的半径，dot(nor, di / l)表示di与平面法线的夹角。返回的遮挡系数等于 1.0 - max(0.0, dot(nor, di / l)) * sph.w * sph.w / (l * l) ，如果仔细观察这个公式不难发现，当夹角越小(dot越大)，或者球半径sph.w越大，或者距离l越小时，max(0.0, dot(nor, di / l)) * sph.w * sph.w / (l * l) 越大，也就是 1.0 - max(0.0, dot(nor, di / l)) * sph.w * sph.w / (l * l) 这个遮挡系数越小。换句话说，当球半径越大，或者距离越近，或者球体位置接近平面的法线时，球体对于交点遮挡的光线就越多，返回的遮挡系数就越小。

tip：个人觉得这只是一种估计遮挡系数的技巧，一般的光线追踪都是通过生成shadow光线来判断可见性的（可以看《PBRT》的相关章节）

知道了oSphere函数的作用后，我们应该就能理解当下面的代码了

// 1.1 判断射线与平面是否相交   
float h = (0.0-ro.y)/rd.y;
    if( h>0.0 ) 
    {       
                tmin = h; 
        nor = vec3(0.0,1.0,0.0);   //平面法线
        pos = ro + h*rd;               //射线与平面交点
        occ = oSphere( pos, nor, sph1 ) *
              oSphere( pos, nor, sph2 );  //求出两个球的遮挡影响系数occ
        obj = 0.0;   //obj = 0.0标记光线与平面相交了
    }
}

因为场景中一共有两个球，遮挡系数occ是sph1和sph2遮挡系数的乘积。

1.2 判断射线与球体是否相交

    h = iSphere( ro, rd, sph1 );
    if( h>0.0 && h<tmin ) 
    { 
        tmin = h; 
        pos = ro + h*rd;
        nor = normalize(pos-sph1.xyz); 
        occ = (0.5 + 0.5*nor.y) *
              oSphere( pos, nor, sph2 );
        obj = 1.0;
    }
    h = iSphere( ro, rd, sph2 );
    if( h>0.0 && h<tmin ) 
    { 
        tmin = h; 
        pos = ro + h*rd;
        nor = normalize(pos-sph2.xyz); 
        occ = (0.5 + 0.5*nor.y) *
              oSphere( pos, nor, sph1 );
        obj = 2.0;
    }

iSphere函数判断射线与球体是否相交：

float iSphere( in vec3 ro, in vec3 rd, in vec4 sph )
{
    vec3 oc = ro - sph.xyz;
    float b = dot( oc, rd );
    float c = dot( oc, oc ) - sph.w*sph.w;
    float h = b*b - c;
    if( h<0.0 ) return -1.0;
    return -b - sqrt( h );
}

推导过程不分析了，有兴趣的话可以看下我前面写的《pbrt笔记--第三章 Shapes》3.2.2节的推导。当iSphere函数返回值大于0时，则表示射线与球体相交。

当相交时，同样要判断场景中其他物体的遮挡影响：occ = (0.5 + 0.5*nor.y) *
oSphere( pos, nor, sph2 );其中(0.5 + 0.5*nor.y) 是平面遮挡的影响（总感觉怪怪的），oSphere( pos, nor, sph2 )是sph2的遮挡影响。另一个球的相交测试也类似，不赘述了。

上面还有一个地方要注意的是，如果射线既与平面相交，又与球体相交，那么怎么处理呢？如下图：

假设射线和平面有一个交点pos2，对应的射线的系数为h2；射线与球体也有一个交点pos1，对应的射线的系数为h1。这个时候，如果h1比h2小，说明射线先和球体相交，那么射线前进的方向就被挡住了，也就不存在pos2这个交点。代码中是通过tmin来记录的。

    h = iSphere( ro, rd, sph1 );
    if( h>0.0 && h<tmin )   //tmin是由上一步计算交点得出的，只有当h小于tmin时，球和射线才有交点
    { 
///...
    }

1.3 计算阴影

相关代码如下：

    vec3 col = vec3(0.93);
    if( tmin<100.0 )  //如果有交点
{
         pos = ro + tmin*rd;   //计算出交点位置

        // shadows
        float sha = 1.0;
        sha *= ssSphere( pos, lig, sph1 );
        sha *= ssSphere( pos, lig, sph2 );
}

经过上面的步骤，如果射线与场景中的物体存在交点（if( tmin<100.0 )），我们求出最近的交点pos后，就开始制作阴影效果。

阴影效果主要由ssSphere函数返回的sha系数来决定的。我们先不看ssSphere的实现，简单说下如何得出阴影效果，如下图：

假设经过前面的计算，我们得出平面上一个交点pos（这个交点是相机发出的射线和平面的交点，而不是光源和平面的交点）。pos是否处于阴影中，实际上就是看主光源lig的光是否可以直接到达pos上。图中左边的pos被球挡住了，那么pos就应该处于阴影之中，反之，右边的pos被lig直接照射到，应该具有较强的亮度。现在我们看等式sha *= ssSphere( pos, lig, sph1 )；pos为交点，lig为光源，sph1为球体1，我们大概已经猜到，ssSphere函数的作用就是lig和pos的光线是否和shp1有交点！（有交点说明被遮挡了）

float ssSphere( in vec3 ro, in vec3 rd, in vec4 sph )
{
    vec3 oc = sph.xyz - ro;
    float b = dot( oc, rd );  //rd和oc的夹角的余弦
    
    float res = 1.0;
    if( b>0.0 )   //
    {
        float h = dot(oc,oc) - b*b - sph.w*sph.w;  //注意h小于0的时候是有交点的
        res = clamp( 16.0 * h / b, 0.0, 1.0 );
    }
    return res;
}

首先会判断oc和rd的夹角的余弦值，当b小于0.0时，说明夹角超过了90度或者小于0度，由下图可以看出，球无论如何都无法遮挡住光线，返回值res就是1.0.

如果夹角小于90度时，球就有可能挡住要到达ro的光线。那么就判断从ro出发，方向为rd的射线是否和球有交点。如果有交点（h<0时），那么公式clamp( 16.0 * h / b, 0.0, 1.0 );就会返回0。如果没有交点（h>0时，），res的值就为16.0 * h / b（为啥是这个还不是很清楚）。

将各个光源的效果进行叠加

代码如下：

        // integrate irradiance with brdf times visibility
        vec3 diffColor = vec3(0.18);
        if( obj>1.5 )   //sph1
        {
            lin += vec3(0.5,0.7,1.0)*diffColor*occ;
            lin += vec3(5.0,4.5,4.0)*diffColor*OrenNayar( lig, nor, rd, 1.0 )*sha;
            lin += vec3(1.5,1.5,1.5)*diffColor*OrenNayar( bac, nor, rd, 1.0 )*occ;
            lin += vec3(1.0,1.0,1.0)*diffColor*OrenNayar( bou, nor, rd, 1.0 )*occ;
        }
        else   //sph2
        {
            lin += vec3(0.5,0.7,1.0)*diffColor*occ;
            lin += vec3(5.0,4.5,4.0)*diffColor*Lambert( lig, nor )*sha;
            lin += vec3(1.5,1.5,1.5)*diffColor*Lambert( bac, nor )*occ;
            lin += vec3(1.0,1.0,1.0)*diffColor*Lambert( bou, nor )*occ;
        }

有代码可知，对sph1使用的是OrenNayar模型的BRDF，对sph2使用的是Lambert模型的BRDF。

gamma校正

因为我们现在计算的亮度是在线性空间，而显示器会对该亮度调整到gamma2.2空间，所以我们在把亮度传给显示器前要把它做一次gamma校正。（gamma校正的问题可以参看https://zhuanlan.zhihu.com/p/66558476
）

    // gamma    
    col = pow( col, vec3(0.45) );