Unity实现GPU光追——Part2全局光照与自发光

maltadirk

本文的主要源于http://three-eyed-games.com/2018/05/12/gpu-path-tracing-in-unity-part-2/
旨在对于该文章的学习解读，但由于原文仅有三章，之后将会是本人自身的一些内容，下面正篇开始
一、渲染方程

首先我们先来看一下渲染方程的主体




我们可以拆解一下该方程，我们的最终目的是计算屏幕中一个像素的亮度。
渲染方程可以解析为以下几项
L(x,\vec{\omega}_0) ，从x点射入人眼的光亮
L_e(x,\vec{\omega}_0) ，物体在x点向人眼方向的自发光
后面是个半球积分，即法向半球内所有击中x点的光对L的影响（暂时不考虑透明物体的情况）
第一部分 f_r 是BRDF（双向反射分布函数）用以描述物体的材质，即入射方向到出射方向的反射比例，不同材质的wi到w0的光量不同。
第二部分 \vec{\omega}_i \cdot \vec{n} 是入射光衰减部分，此处与兰伯特相同
最后 L(x,\vec{\omega}_i) 是递归部分，wi的光来源于与 L(x,\vec{\omega}_0) 相同
所以该方程是个无限递归的方程
使用蒙特卡洛积分

通过蒙特卡洛积分可以使我们进行有限的采样来求解上面的方程，它可以保证结果的收敛正确，它的基本形式如下


采样n个样本，除以每个样本被采样的概率。这样经常选中的样本权重就会小于很少选择的样本，由于半球的球面角为2pi，所以概率密度为1/2pi。将蒙特卡洛积分替换原来的半球积分，我们可以得到如下公式


二、准备工作

样本累积

由于Unity无法将HDR纹理作为渲染目标，所以我们的格式是 R8G8B8A8_Typeless ，但这造成了精度问题，无法添加多个样本。所以我们创建一个RenderTexture _converged，这是输出图像前的高精度缓冲，在Render函数中如下修改。
            if (_target != null)
            {
                _target.Release();
                _converged.Release();
            }
            
            _target = new RenderTexture(Screen.width, Screen.height, 0,
                RenderTextureFormat.ARGBFloat, RenderTextureReadWrite.Linear);
            _converged = new RenderTexture(Screen.width, Screen.height, 0,
                RenderTextureFormat.ARGBFloat, RenderTextureReadWrite.Linear);
            
            _target.enableRandomWrite = true;
            _converged.enableRandomWrite = true;
            _target.Create();
            _converged.Create();
在渲染环节，先渲染到_converged中，再输出到dest
        Graphics.Blit(_target,_converged,_addMaterial);
        Graphics.Blit(_converged,dest);固定场景

我们之前的渲染中，场景是通过Onenable生成的，但有对比会更好，我们可以通过加入Seed来控制随机数，在创建场景前的Random中加入Seed
//构造函数
        public CreateSpheres(int Seed)
    {
        this.SphereSeed = Seed;
    }

//Create下
Random.InitState(SphereSeed);着色器随机值

在Shader中我们根据当前像素和传入的Seed作为参数生成随机值
float rand()
{
    float result = frac(sin(_Seed / 100.0f * dot(_Pixel, float2(12.9898f, 78.233f))) * 43758.5453f);
    _Seed += 1.0f;
    return result;
}三、半球采样

我们接下来需要根据法向量生成一个在半球上的随机向量，通过它来进行采样，我们先以z轴正方向为法向量生成采样，再对其进行 切线空间-世界空间转换 即可。如果了解过SSAO的原理想必会十分熟悉。
float3x3 GetTangentSpace(float3 normal)
{
    float3 helper = float3(1,0,0);
    if(abs(normal.x) > 0.99f)
        helper = float3(0,0,1);

    //生产切线和副切线
    float3 tangent = normalize(cross(normal,helper));
    float3 binormal = normalize(cross(normal,tangent));
    return float3x3 (tangent,binormal,normal);
}

float3 SampleHemisphere(float3 normal)
{
    //根据两个角度确定向量
    float cosTheta = rand();
    float sinTheta = sqrt(max(0.0f, 1.0f - cosTheta * cosTheta));
    float phi = 2 * PI * rand();
    float3 tangentSpaceDir = float3(cos(phi) * sinTheta, sin(phi) * sinTheta,cosTheta);

    //坐标转换
    return mul(tangentSpaceDir,GetTangentSpace(normal));
}四、兰伯特漫反射

现在我们有了半球随机采样，这下我们就可以实施第一个BRDF了。我们先来尝试最简单的Lambert BRDF。它的表述为 f_r(x,\vec{\omega}_i,\vec{\omega}_i)= \cfrac{k_d}{\pi} ,kd为表面反照率albedo，由此渲染方程变为了


们在shader的Shade函数中，将hit.distance < 1. #INF中的代码替换如下。反射出的方向为随机反射，能量为2*albedo，由于我们是在时间上累加，所以我们不需要考虑N的值
        ray.origin = hit.position + hit.normal * 0.001f;
        ray.direction = SampleHemisphere(hit.normal);
        ray.energy *=  2 * hit.albedo * sdot ( hit.normal, ray.direction ) ;
        return 0;sdot为一个工具函数
float sdot(float3 x, float3 y, float f = 1.0f)
{
    return saturate(dot(x, y) * f);
}回顾我们做了什么。CMain发射光线，调用Shade，如果表面被击中就生成新的射线（半球采样），同时将BRDF和余弦值分解为能量，如果击中天空就对天空球采样，最后我们对样本进行混合。
记得在创建场景时将Metal设置为Random.value < 0.0f，禁止金属类材质生成



效果图

五、Phong高光

在光栅化中我门也是用过Phong模型作为高光，不给过原始通过次方进行高光的方法不是能量守恒的，但仅需简单的修改即可重新满足条件，公式如下


下面是它的函数图像


https://www.desmos.com/calculator/h7m3esuf28?lang=zh-CN
我们将其带入渲染方程可以得到


我们还可以进一步加上兰伯特光照


将其写入我们的Shade中
        //考虑浮点数精确性进行偏移
        ray.origin = hit.position + hit.normal * 0.001f;
        //反射方向
        float3 reflected = reflect(ray.direction,hit.normal);
        ray.direction = SampleHemisphere(hit.normal);
        //漫反射部分
        float3 diffuse = 2 * min(1.0f - hit.specular,hit.albedo);
        //粗糙度
        float a = 300.0f;
        //Phong着色
        float3 specular = hit.specular * (a + 2) * pow(sdot(ray.direction, reflected), a);
        ray.energy *=  (diffuse + specular) * sdot(hit.normal,ray.direction);
        return 0;然后重启金属类材质，最终在a=15和a=300的效果如下



a=15



a=300

可以看出，a=300的情况下有许多的噪点。
这是因为我们所使用的的渲染方程是个方差极大的方程，以a=300为例，当采样的光线防线与反射光线方向完美重合时，该值可以达到数百，因此需要极长的时间才能将其拉回正常值，而在等在中又会产生新的噪点。这需要当N极大时才能正常收敛
六、更好的采样

这里我们暂时放弃Phong光照，考虑一种能将完全镜面反射和漫反射结合的方法，因为完全镜面反射是绝对正确切不会变动的。
我们根据ks和kd来计算两者的能量占比，通过随机数来随机采样。但这样总能量会大幅下降，因为我们只渲染了其中一部分，因此我们还需要补充另一部分的能量。
但这样，我们也遇到了新的问题，我们失去了模拟粗糙金属的方法。
我们用energy来计算能量
float energy(float3 color)
{
    //平均颜色
    return dot(color,1.0f / 3.0f);
}
        //考虑浮点数精确性进行偏移
        ray.origin = hit.position + hit.normal * 0.001f;

        //计算漫反射和镜面反射的几率
        hit.albedo = min(1.0f - hit.specular,hit.albedo);
        float specChance = energy(hit.specular);
        float diffChance = energy(hit.albedo);
        float sum = specChance + diffChance;
        specChance /= sum;
        diffChance /= sum;

        //随机选择渲染类型
        float roulette = rand();
        if(roulette < specChance)
        {
            //发生镜面反射
            ray.direction = reflect(ray.direction,hit.normal);
            ray.energy *= (1.0f / specChance) * hit.specular * sdot(hit.normal , ray.direction);
        }
        else
        {
            //发生漫反射
            ray.direction = SampleHemisphere(hit.normal);
            ray.energy *= (1.0f / diffChance) * 2 * hit.albedo * sdot(hit.normal,ray.direction);
        }
        return 0;


效果图

七、重要性采样

让我们回到蒙特卡洛公式


这里每个样本要除上样本被选择的概率。由于我们目前使用的是半球均匀采样，因此p为常数1/2pi。但很明显，均匀的采样并非最佳的，例如在Phong BRDF的情况下，他在趋近反射方向的狭窄区间上有很大的值。
如果我们找到一个和被积函数分布完全相同的概率分布： p(x) = k \cdot f(x)
那么最终的结果将会是一个常数，这将使结果以很快的速度收敛。当然找到这样的分布是及其困难的，但我们仍可以从BRFD x cosine中找到一些分布。
余弦采样

我们用余弦的幂替换我们的均匀分布样本，我们的目的是按比例生成更少的样本。
Thomas Poulet的这篇文章描述了原理。我们想采样函数中加入alpha参数来确定余弦采样的强度，0为均匀采样，1为余弦采样，或更高的Phong指数
float3 SampleHemisphere(float3 normal,float alpha = 0)
{
    //根据两个角度确定向量,其中alpha是对采样的变换
    float cosTheta = pow( rand(),1.0f / (alpha + 1.0f));
    float sinTheta = sqrt(max(0.0f, 1.0f - cosTheta * cosTheta));
    float phi = 2 * PI * rand();
    float3 tangentSpaceDir = float3(cos(phi) * sinTheta, sin(phi) * sinTheta,cosTheta);

    //坐标转换
    return mul(tangentSpaceDir,GetTangentSpace(normal));
}重要性采样兰伯特

我们可以通过包含余弦因子来使采样有更好的效果，令a = 1，则我们的分布函数变为了 \frac{(\vec{\omega}_i,\vec{n})}{\pi} 带入渲染方程我们可以得到


            //发生漫反射
            ray.direction = SampleHemisphere(hit.normal,1.0f);
            ray.energy *= (1.0f / diffChance) * hit.albedo;重要性采样Phong

对于Phong BRDF也是类似的过程，我们有两个余弦的乘积，来自BRDF的幂余弦和渲染方程的正余弦


我们仅处理幂余弦，我们通过对 \vec{\omega}_0 方向将粗糙度作为值传入，可以化简得到以下公式，具体推导详见


            //发生镜面反射
            float alpha = 15.0f;
            ray.direction = SampleHemisphere(reflect(ray.direction,hit.normal),alpha);
            float f = (alpha + 2)/(alpha + 1);
            ray.energy *= (1.0f / specChance) * hit.specular * sdot(hit.normal,ray.direction,f);同样是300的粗糙度，可见效果好了许多


八、后续

接下来我们需要将粗糙度写入我们的球体，同时我们还可以加上一个自发光，记得在C#中也要更新数据结构以及步长
struct Sphere
{
    float3 position;
    float radius;
    float3 albedo;
    float3 specular;
    float smoothness;
    float3 emission;
};
同时不要忘记更新材质赋值
void IntersectSphere(Ray ray,inout RayHit bestHit,Sphere sphere)
{
    //计算射线与球体相交的距离
    //由线方程与求方程连列求解得
    //https://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
    float3 d = ray.origin - sphere.position;
    float p1 = -dot(ray.direction,d);
    float p2sqr = p1 * p1 - dot(d,d) + sphere.radius * sphere.radius;
    if(p2sqr < 0)
        return;
    float p2 = sqrt(p2sqr);
    float t = p1 - p2 > 0 ?
              p1 - p2 :
              p1 + p2 ;
    if(t > 0 && t< bestHit.distance)
    {
        bestHit.distance = t;
        bestHit.position = ray.origin + t * ray.direction;
        bestHit.normal = normalize(bestHit.position - sphere.position);
        bestHit.albedo = sphere.albedo;
        bestHit.specular = sphere.specular;
        bestHit.smoothness = sphere.smoothness;
        bestHit.emission = sphere.emission;
    }
}通过一个函数将粗糙度映射为我们用的alpha
float SmoothnessToPhongAlpha(float s)
{
    return pow(1000.0f,s * s);
}关于自发光，当物体有自发光时直接返回自发光的值即可