实时渲染：使用OpenGL进行蒙特卡洛路径追踪

先上个成果图：

使用OpenGL和GLSL实现了一个基于物理的蒙特卡洛路径追踪，接下来简述一下实现的原理和思路。仓库地址：RealtimeRayTracing

OpenGL端

在OpenGL端，采用了简单的ping-pong texture和一个final shading来实现了逐帧增加spp数量的效果：

这一段的核心代码大致如下：

// 绘制pathtracer framebuffer
if (m_cam_move_flag) {
    glClear(GL_COLOR_BUFFER_BIT);
    glClearColor(GL_RGBA_BLACK);
    m_cam_move_flag = false;
}

m_rtrt_shader.use();
m_rtrt_shader.setTexture("uLastFrame", 0, m_screen_texture[pingpong]);
m_rtrt_shader.setUniform("uPingpong", pingpong ^ 1);

m_rtrt_shader.setUniform("uRand3", randomVec3());
m_spp += m_light_samples;
m_quad.draw();

// 绘制screen framebuffer
m_fbo.unbind();
m_final_shader.use();
m_final_shader.setTexture("uRenderTexture", 0,
                          m_screen_texture[pingpong ^ 1]);
m_final_shader.setUniform("uSPP", std::max(1, (int)m_spp));
m_quad.draw();
pingpong = pingpong ^ 1;

需要注意的是，这里应该是必须要使用pingpong texture交替渲染才能实现，因为同时读写一个texture是未经过定义的行为。

总的来说，由于我暂时还没有实现渲染CPU端读入的模型，在OpenGL端并没有什么实现难度。

Shader实现

顶点着色器

在顶点着色器中，需要做两件事：

1. 计算出四个顶点的光线的origin和dir，在fragment shader中插值使用。
2. 根据摄像机uv向量和framebuffer大小计算出每个像素的随机采样范围。

这里就需要搞清楚几个概念：Camera, Film, Framebuffer，首先，我们是通过摄像机发射光线，形成透视的光线到世界中进行采样，那么应该从什么方向去发射光线呢？类似于摄影的原理，我们使用一张“胶片”来保存所有的信息，也就是所谓的Film，我们规定好摄像机光线起点的位置作为origin，在摄像机正前方放上一张胶片，对于胶片上的每个点发射光线即可，而后将胶片上的各个光线着色点映射到Framebuffer的每个像素上，这样就可以很清楚明白地解释成像的原理了！

首先定义摄像机的几个重要属性：

struct Camera {
    vec3 position;
    float fov;
    float aspectRatio;
    vec3 front;
    vec3 up;
};

这里的up并不一定是垂直于front方向的，我们只需要利用这个向量来计算一个摄像机的uvw坐标系（这里的w和摄像机的front方向是相反的）：

vec3 w = normalize(-uCamera.front);
vec3 u = normalize(cross(uCamera.up, w));
vec3 v = cross(w, u);

接下来，我们假设摄像机到胶片的距离为1.0，则可以根据fov和aspect ratio来计算出胶片的大小：

vec2 calcViewportSize() {
    // distance == 1.0
    float h = tan(uCamera.fov / 2.0);
    float height = 2.0 * h;
    float width = uCamera.aspectRatio * height;
    return vec2(width, height);
}

根据uvw构成的坐标系，可以计算出胶片的lowerLeft坐标，再将lowerLeft坐标加上根据纹理坐标得到的胶片像素偏移（因为OpenGL里传入的模型为一个屏幕quad），即可得到每个胶片像素的位置：

vec3 lowerLeft =
        uCamera.position - 0.5 * u * filmSize.x - 0.5 * v * filmSize.y - w;
rayOrigin = uCamera.position;
rayDir = lowerLeft + u * aTexCoord.x * filmSize.x +
          v * aTexCoord.y * filmSize.y - rayOrigin;

因为我们希望在之后的片段着色器中实现多次采样的时候能够实现在像素中随机地采样光线，我们还需要将u、v方向的变化范围传入fragment shader。已知的是，OpenGL对于每个像素着色时，实际着色的坐标为(x + 0.5, y + 0.5)，我们希望在生成[-1, 1]的随机数的时候，光线在(x, y) - (x + 1.0, y + 1.0)的范围中随机采样，于是我们的range范围便应该是\(frac{uv * 0.5 * filmeSize}{framebufferSize}\)：

// flat表示输出到fragment shader中不需要进行插值
flat out vec3 vHorizontalRange;
flat out vec3 vVerticalRange;
vHorizontalRange = u * 0.5 * filmSize.x / uFramebufferSize.x;
vVerticalRange = v * 0.5 * filmSize.y / uFramebufferSize.y;

片段着色器

首先这里要讲解一下蒙特卡洛路径追踪的几大特征：

• 直接空间中根据几何信息进行求交（和RayMarching的步进策略不同），得到的是准确的几何解。
• 使用SPP从摄像机发射大量光线，代替从反射点大量发射光线，来实现蒙特卡洛采样。
• 递归追踪光线，直到俄罗斯轮盘赌判定为结束/没有击中任何物体/达到最大反射次数时结束（这里需要提及的是，如果按照正确的渲染方程，在击中有\(L_e\)自发光的物体的情况下仍然需要进行采样和追踪，但是我们为了简化追踪流程，在光线击中了灯光后不会继续反射，而是直接结束弹射）

首先祭出祖师爷渲染方程：

\[L_o = L_e + \int_\Omega f_r\cdot L_i\cdot(\omega_i\cdot \vec{n})\text{d}\omega_i\]

其中\(f_r\)项就是大名鼎鼎的BRDF函数了，这个函数接受入射方向、反射方向、法线方向以及材质信息，计算光线经过这一次弹射的能量变换（颜色的出现本质上也是不同材料对于不同波长光的吸收比率不同而有所不同，因此颜色的不同本质上也是brdf的不同），但是我们这里为了偷懒简化模型，和ray tracing in one weekend中一样，认为\(brdf = color\)，即，认为采样的每一条反射/折射光线都只被吸收了vec3(1.0) - vec3(color)的能量，虽然并不是完全的物理正确，但是对于一个简单的路径追踪也是够用的。

剩下的就很简单了，这里引用一张GAMES101的图片来说明整个路径追踪流程：

即，我们需要递归地对光线打到的每个点进行着色，然后对结果进行累乘，得到的即是最终的结果。

但是，问题就来了：路径追踪的核心显而易见：递归，glsl中并不支持递归（甚至不是不支持，而是glsl的内存模型禁止递归，这很好理解，GPU的显存+并行化的模式不允许有过深的函数栈存在）。那么怎样才能实现这样的递归方法呢？我们观察路径追踪的流程，发现其实整个过程中并不需要有辅助栈的存在——即使实现为递归，由于每一次递归调用均在函数末尾，因此编译器也会进行尾递归优化。我们只需要使用两个辅助变量：color和brightness即可使用循环替代整个递归流程，伪代码如下：

eval

vec3 brightness = vec3(1.0), color = vec3(0.0);
float p_RR = 0.8;
Intersection inter;
for(int i = 0; i < MAX_BOUNCE; i++){
  if (random() > p_RR) break;
  if(hasIntersection(ori, dir, inter)){
    if(inter.isLight){
      color = inter.color * brightness;
    }else{
      vec3 new_dir = randomSample(pdf);
      float cosine = dot(inter.normal, new_dir);
      brightness = brightness * inter.color * cosine / pdf / p_RR;
      ori = inter.position;
      dir = new_dir;
    }
  }else{
    break;
  }
}
return color;

glsl的实现如下：

vec3 recursivePathTracing(vec3 origin, vec3 dir) {
    vec3 rayColor = vec3(0.0);
    vec3 rayBrightness = vec3(1.0);
    float p_RR = 0.8;

    for (int i = 0; i < uLightBounceCount; i++) {
        if (random(randSeed.xy) > p_RR) break;
        Intersection inter;
        if (!hasIntersect(origin, dir, inter)) {
            break;
        }
        Material material = inter.material;
        if (material.isLight) {
            rayColor = material.color * rayBrightness;
            break;
        } else {
            float cosine = 1.0;
            float pdf = 1.0;
            // add a slight offset from hit position to prevent artifacts
            // https://computergraphics.stackexchange.com/questions/7789/weird-artifacts-in-my-ray-tracer
            origin = inter.position + inter.normal * EPS;
            vec3 rs = normalize(random3(randSeed));
            if (material.type == DIFFUSE) {
                dir = normalize(randomHemispherePoint(rs, inter.normal));
                pdf = 2.0 * PI;
                cosine = dot(dir, inter.normal);
            } else if (material.type == METALIC) {
                dir = reflect(dir, inter.normal);
                cosine = dot(dir, inter.normal);
            } else if (material.type == GLASS) {
                float si = dot(dir, inter.normal);
                float eta = 0.66;
                if (si < 0.0) {
                    float rd = random(randSeed.xy) + 1.0;
                    rd /= 2.0;
                    // incident
                    if (rd <= fresnelApprox(normalize(dot(-dir, inter.normal)),
                                            1.0, 1 / eta)) {
                        dir = reflect(dir, inter.normal);
                        cosine = dot(dir, inter.normal);
                    } else {
                        dir = refract(dir, inter.normal, eta);
                        origin = inter.position - inter.normal * EPS;
                        cosine = dot(dir, -inter.normal);
                    }
                } else {
                    // exitent
                    dir = refract(dir, -inter.normal, 1 / eta);
                    cosine = abs(dot(dir, inter.normal));
                }
            }
            cosine = max(cosine, 0.0);
            rayBrightness *= material.color * cosine / p_RR / pdf;
            randSeed = rs;
        }
    }
    return rayColor;
}

这里需要注意的是为了防止光线在打到物体表面后，由于float精度问题，弹射光线有可能被物体自己拦截而产生artifact。使用沿法线方向位移一小段距离的方式来防止自遮挡的问题。

小结

目前的路径追踪器虽然已经能运行出较为不错的效果了，但是离一个真正的路径追踪渲染器还差不少，接下来会继续添加如下的功能：

• 基于三角形的模型
• 低spp下对光源进行采样以优化效果
• 物理正确的brdf模型

PS：当初我入门图形学的第一个项目就是跟着raytracing in one weekend做的，现在回忆起来，当时做的是一知半解，连做个并行化都很费劲，但是现在自己可以独立地手写一个路径追踪器了，颇为感慨。

Update 2022-4-4

今天在为追踪器实现cosine-weight采样的时候发现了实现中的几个问题，这里需要提及一下：

BRDF & PDF & cosine

已知渲染方程是上面的形式，cosine项肯定是去不掉的，但是在之前的实现里，pdf和brdf的实现其实是有误的：pdf作为采样的概率密度函数，代表的应该是每个采样射线出现的概率，对于一个均匀半球采样，其值应该是\(1.0 / 2\pi\)，之前误写为了\(2\pi\)，但是如果直接更改会发现出现了颜色的过亮，这正是因为BRDF也不正确，diffuse材质下，Lambertian brdf模型规定brdf应该为\(\frac{\sigma}{\pi}\)——正好抵消了这个多出来一倍的\(\pi\)（虽然这么理解肯定是不正确的）。

随机采样函数

随机采样函数接收的种子是\((x, y, z) \in[-1, 1]\)，而且xyz三个分量应该都是独立均匀分布的，之前粗心直接用了一个normalize导致得到的范围是有问题的！

To be continued...