引言

在前面几节，我们完成了若干项任务：定义着色器实例、搭建管线、绑定SBT、完成模型纹理读入、完成二次弹射的追踪shader、完成实时降噪。好像我们的光追渲染器已经成型了，然而事实上还差的远，因为我们现在的光追就是：光线第一次求交得到命中物体，然后第二次直接定向找光源，这样弹射两次就结束了。真正的光追绝对不是弹射两次光线就能计算出结果的，因为这样的算法与光栅化差异不大（光栅化也是弹射两次，区别在于光栅化是从光源到物体再到眼睛，这里的光追是从眼睛到物体再直接追踪到光源）。因此我们要将其拓展为多次弹射、对路径作渲染方程积分估计的真正的光追。

不过从本节开始，已经脱离了官方课程optix course的范畴了，这意味着从本节开始每一句代码都需要自己摸索，并适当的参考pbrt、mitsuba等开源渲染器中的一些架构和api，因此本节开始的代码难免会出现一些小bug，如果大家发现了可以在评论区指出，万分感谢。

Part1. 理论准备

回忆一下我们在C++光追渲染器中是怎么定义光追程序的。我们定义了一个Ray结构体，用来记录当前追踪的光线的起点、终点；又定义了一个Intersection结构体，用来记录与三角面求交得到的交点位置、法线、模型材质等等信息。同理，optix路径追踪也离不开这俩东西，一会儿我们也要在shader里定义它们。

然后我们在C++光追渲染器写了一个递归程序，从相机视口的某个像素出发，递归式地找到每一个交点、直到最后找到光源（或求交失败）作为最深一层，然后我们从最深一层开始赋予光照能量，再回溯去求每一层的能量，最后得到有多少能量被相机捕捉到。我们可以分析一下中间的能量计算过程：

假如我们从某像素射出一根光线，先打到物体A，再打到物体B，再打到物体C，最后追踪到光源D，那么该像素对应的这根光线最终是多少能量呢？（下文用e表示自发光，用f表示渲染方程的积分部分）

按照正常的光路顺序，光线从D发出，能量是eD，打到C表面上后调用C的渲染方程，可以得到出射光强是eC + fC * eD，然后打到B表面调用渲染方程，得到出射光是eB + fB(eC + fC * eD)，以此类推……便可以得到e = eA + fA(eB + fB(eC + fC * eD))。所以我们倾向于给光追过程编写一个递归程序，先递归地找到最深一层（D）的出射能量，得到D的能量后回溯到C层计算C的出射能量，再回溯到B计算B的出射能量，以此类推……所以大部分光追程序都是写成递归的形式。

然而放在optix渲染器里就不一样了。尽管如今cuda已经支持递归了，但是在gpu端写递归函数终归会拖累性能，而且容易出现栈溢出。因此我们不得不将光追迭代从递归函数 强行改成for循环迭代的方式。那么问题来了，计算光照我们得从最深一层D开始计算，然后再算C、B、A，而for循环只能正向从A开始一步一步求交，当求交到D时，ABC的信息早就被抹了，还怎么一步步往回代入计算呢？

这里我们就将刚才的公式拆解一下：

e = eA + fA(eB + fB(eC + fC * eD))

e = eA + fA * eB + fA * fB * eC + fA * fB * fC * eD

看到这里，算法就很明朗了。我们直接在for循环中做弹射（每一步弹一次并求交），维护一个e作为总能量（初始=0），维护一个accum值作bsdf函数的累乘值（初始=1）。每求交到一个物体，我们直接给e+=该材质自发光*accum的值，然后给accum乘一个该材质的bsdf函数值，进入下一次迭代就可以了。

明白这部分后，基本就没有什么要补充的知识了，我们直接动工！

Part.2 朴素路径追踪shader

首先我们创建一个MyInteraction.h文件，放一下我们要定义的Ray和Interaction：

MyInteraction.h：

#include "gdt/math/vec.h"
#include "LaunchParams.h"
#include "math.h"

using namespace osc;

struct Ray
{
    vec3f origin;
    vec3f direction;
    float tmax = FLT_MAX;
};

struct Interaction
{
    float bias = 0.001f;
    float distance;
    vec3f position;
    vec3f geomNormal;
    vec3f mat_color;

    __forceinline__ __device__ Ray spawn_ray(const vec3f& wi) const
    {
        vec3f N = geomNormal;
        if (dot(wi, geomNormal) < 0.0f)
        {
            N = -geomNormal;
        }

        Ray ray;
        ray.origin = position + N * bias;
        ray.direction = wi;
        ray.tmax = FLT_MAX;
        return ray;
    }
};

这里参考了一个KDRay的架构设计，非常感谢这位大佬的分享 : )

首先定义一个Ray结构体，origin表示射线起点，direction表示射线方向，tmax表示射线能追踪的最远距离；然后定义一个Interaction结构体来记录求交信息，distance表示两次求交点之间的距离，position表示交点位置，geomNormal表示交点的法线，mat_color表示交点的材质颜色（这一节我们的bsdf都是超级简单的纯色材质）。额外定义的spawn_ray函数讲的是，根据出射方向来创建一个Ray，bias就是为了防止暗疮现象而加的出射偏移。

接下来我们来修改着色器。因为目前我们不需要阴影射线，所以我们直接把shadow相关的shader直接置空。然后我们看一下Raygen Shader：

extern "C" __global__ void __raygen__renderFrame()
    {
        // compute a test pattern based on pixel ID
        const int ix = optixGetLaunchIndex().x;
        const int iy = optixGetLaunchIndex().y;
        const auto& camera = optixLaunchParams.camera;

        int numPixelSamples = optixLaunchParams.numPixelSamples;

        PRD prd;

        vec3f pixelColor = 0.f;
        // normalized screen plane position, in [0,1]^2
        vec2f screen(vec2f(ix, iy) / vec2f(optixLaunchParams.frame.size));

        // generate ray direction
        vec3f rayDir = normalize(camera.direction
            + (screen.x - 0.5f) * camera.horizontal
            + (screen.y - 0.5f) * camera.vertical);

        for (int sampleID = 0; sampleID < numPixelSamples; sampleID++) {
            Ray myRay;
            myRay.origin = camera.position;
            myRay.direction = rayDir;
            vec3f radiance = 0.0f;
            vec3f accum = 1.0f;
            for (int bounces = 0; ; ++bounces)
            {
                if (bounces >= optixLaunchParams.maxBounce) {
                    radiance = 0.0f;
                    break;
                }
                Interaction isect;
                isect.distance = 0;
                unsigned int isectPtr0, isectPtr1;
                packPointer(&isect, isectPtr0, isectPtr1);
                optixTrace(optixLaunchParams.traversable,
                    myRay.origin,
                    myRay.direction,
                    0,    // tmin
                    1e20f,  // tmax
                    0.0f,   // rayTime
                    OptixVisibilityMask(255),
                    OPTIX_RAY_FLAG_DISABLE_ANYHIT,//OPTIX_RAY_FLAG_NONE,
                    RADIANCE_RAY_TYPE,            // SBT offset
                    RAY_TYPE_COUNT,               // SBT stride
                    RADIANCE_RAY_TYPE,            // missSBTIndex 
                    isectPtr0, isectPtr1);
                if (isect.distance == FLT_MAX)
                {
                    radiance += vec3f(1.0f) * accum;
                    break;
                }
                radiance += 0.0f;
                accum.x *= isect.mat_color.x;
                accum.y *= isect.mat_color.y;
                accum.z *= isect.mat_color.z;
                // 下一次反弹方向
                vec3f wi;
                vec3f rnd;
                prd.random.init(optixLaunchParams.frame.frameID * 234834 % 32849 + ix * 385932 % 82921,
                    optixLaunchParams.frame.frameID * 348593 % 43832 + iy * 324123 % 23415);
                rnd.x = prd.random() * 2 - 1;
                prd.random.init(optixLaunchParams.frame.frameID * 972823 % 12971 + ix * 743782 % 82013, 
                    optixLaunchParams.frame.frameID * 893022 % 28191 + iy * 918212 % 51321);
                rnd.y = prd.random() * 2 - 1;
                prd.random.init(optixLaunchParams.frame.frameID * 383921 % 48839 + ix * 572131 % 47128,
                    optixLaunchParams.frame.frameID * 389291 % 29301 + iy * 716271 % 63291);
                rnd.z = prd.random() * 2 - 1;
                wi = isect.geomNormal + normalize(rnd);
                wi = normalize(wi);
                myRay = isect.spawn_ray(wi);
            }
            pixelColor += radiance;
        }

        vec4f rgba(pixelColor / numPixelSamples, 1.f);
        rgba.x = powf(rgba.x, 1 / 2.2f);
        rgba.y = powf(rgba.y, 1 / 2.2f);
        rgba.z = powf(rgba.z, 1 / 2.2f);
        if (rgba.x > 1)rgba.x = 1.0f;
        if (rgba.y > 1)rgba.y = 1.0f;
        if (rgba.z > 1)rgba.z = 1.0f;
        if (rgba.w > 1)rgba.w = 1.0f;

        // and write/accumulate to frame buffer ...
        const uint32_t fbIndex = ix + iy * optixLaunchParams.frame.size.x;
        if (optixLaunchParams.frame.frameID > 0) {
            rgba
                += float(optixLaunchParams.frame.frameID)
                * vec4f(optixLaunchParams.frame.colorBuffer[fbIndex]);
            rgba /= (optixLaunchParams.frame.frameID + 1.f);
        }
        optixLaunchParams.frame.colorBuffer[fbIndex] = (float4)rgba;
    }

这部分程序还是比较好读的，首先照常根据屏幕像素坐标，得到第一根出射光线的起始坐标和方向，这回我们把它的信息装入自定义的Ray中。

然后第一层循环是采样循环，表示同一个像素要采样多次，这里我们先把总能量radiance和bsdf积累量accum这些变量初始化好；第二层循环就是我们的光路迭代循环，表示现在追踪到了第几层。首先定义一个Intersection准备装求交信息，我们将这个Intersection的指针装入gpu的0号、1号槽位，用来接受hit shader和miss shader传递的信息。然后我们便开始发射射线！在miss shader和hit shader会把求交点的位置、法线、材质颜色都写入Intersection里。这样我们就得到了这一轮的所有求交信息，并更新radiance、accum这些变量。同时，我们也要根据法线和bsdf信息来决定下一次的弹射方向，因为我们这一节都是纯粗糙diffuse材质，因此我们直接进行均匀半球采样，用prd来撒种子求随机数并决定采样方向（随机数方案我目前是乱写的，之后会介绍采样器）。有了这个方向，我们便可以用spawn_ray函数来得到下一次弹射的Ray信息，进入下一个迭代。

我们预设天空是光源。当弹射次数超过预置次数时，表示这根光线追踪不到光源，因此直接给能量置0；当追踪到无限远的天空光源时，就直接给天空自发光设为天光能量值（例如rgb(1,1,1)），然后用我们前面说的公式就可以积累出最终的光线颜色。

对多次采样的结果进行取平均，然后要对颜色做一次gamma矫正（不然整体会偏暗），然后再处理溢出情况，最后和之前的帧结果进行混合，完成实时降噪。

怎么样，思路很清晰吧？我们再来看一下miss shader，更简单：

extern "C" __global__ void __miss__radiance()
    {
        uint32_t isectPtr0 = optixGetPayload_0();
        uint32_t isectPtr1 = optixGetPayload_1();
        Interaction* interaction = reinterpret_cast<Interaction*>(unpackPointer(isectPtr0, isectPtr1));
        interaction->distance = FLT_MAX;
    }

当miss shader被调用，说明该光线已经打到无限远的天空上了。因此我们取出0号、1号槽装的intersection指针，将距离信息（FLT_MAX无限远）写入，这样raygen shader的光路迭代部分就知道这一轮弹射打到天空光源了。

再来看一下closesthit shader：

extern "C" __global__ void __closesthit__radiance()
    {
        uint32_t isectPtr0 = optixGetPayload_0();
        uint32_t isectPtr1 = optixGetPayload_1();
        Interaction* interaction = reinterpret_cast<Interaction*>(unpackPointer(isectPtr0, isectPtr1));
        const TriangleMeshSBTData& sbtData
            = *(const TriangleMeshSBTData*)optixGetSbtDataPointer();
        PRD& prd = *getPRD<PRD>();

        // ------------------------------------------------------------------
        // gather some basic hit information
        // ------------------------------------------------------------------
        const int   primID = optixGetPrimitiveIndex();
        const vec3i index = sbtData.index[primID];
        const float u = optixGetTriangleBarycentrics().x;
        const float v = optixGetTriangleBarycentrics().y;

        // ------------------------------------------------------------------
        // compute normal, using either shading normal (if avail), or
        // geometry normal (fallback)
        // ------------------------------------------------------------------
        const vec3f& A = sbtData.vertex[index.x];
        const vec3f& B = sbtData.vertex[index.y];
        const vec3f& C = sbtData.vertex[index.z];
        vec3f Ng = cross(B - A, C - A);
        vec3f Ns = (sbtData.normal)
            ? ((1.f - u - v) * sbtData.normal[index.x]
                + u * sbtData.normal[index.y]
                + v * sbtData.normal[index.z])
            : Ng;

        const vec3f pos = (1.f - u - v) * A + u * B + v * C;
        interaction->position = pos;
        // ------------------------------------------------------------------
        // face-forward and normalize normals
        // ------------------------------------------------------------------
        const vec3f rayDir = optixGetWorldRayDirection();

        if (dot(rayDir, Ng) > 0.f) Ng = -Ng;
        Ng = normalize(Ng);

        if (dot(Ng, Ns) < 0.f)
            Ns -= 2.f * dot(Ng, Ns) * Ng;
        Ns = normalize(Ns);

        interaction->geomNormal = Ns;
        // ------------------------------------------------------------------
        // compute diffuse material color, including diffuse texture, if
        // available
        // ------------------------------------------------------------------
        vec3f diffuseColor = sbtData.color;
        if (sbtData.hasTexture && sbtData.texcoord) {
            const vec2f tc
                = (1.f - u - v) * sbtData.texcoord[index.x]
                + u * sbtData.texcoord[index.y]
                + v * sbtData.texcoord[index.z];

            vec4f fromTexture = tex2D<float4>(sbtData.texture, tc.x, tc.y);
            diffuseColor *= (vec3f)fromTexture;
        }

        // start with some ambient term
        vec3f pixelColor = (0.1f + 0.2f * fabsf(dot(Ns, rayDir))) * diffuseColor;
        interaction->mat_color = pixelColor;
    }

别看hit shader这么长，90%都是上一节写的计算插值法线、计算插值uv、获取纹理、代入lambert模型计算光照颜色。只是这里我们需要将求交点坐标、法线以及最终的材质颜色写入Intersection中即可。（这里我的设计有点不严谨，严格意义上讲，计算bsdf的代码不应该放在hit shader，因为bsdf依赖于入射、出射两根光线的信息，然而调用hit shader的时候还没开始采样出射光线。只是因为本例是diffuse模型，出射光对bsdf值无贡献，所以我才能把bsdf的计算放在hit shader中。日后但凡有了新材质，肯定都要把bsdf计算搬到光路迭代的主循环中去算）

最后我们把弹射深度设为5~10，便可以看到有深度的光追了。接下来我们导入一个车辆的渲染场景：