作者：i_dovelemon
日期：2024-06-16
主題：Lightmap, PathTracer, Compute Shader
引言
一直以來，我都對離線 bake lightmap 操作很著迷。一方面，這個方案歷久彌新，雖然很古老，但是一直在實際專案中都有使用；另一方面，它能夠產生非常高質量的 GI 效果。
很久之前，我寫過一個基於 Radiosity 演算法的 baker 。那個演算法很簡單，但是卻不夠吸引人。當時我就想以後一定要寫一個基於 PathTracer 的 baker，所以，就有了這個 demo 的產生。
做這個 demo，我按照一定的步驟，迭代式的進行開發：先實現一個 CPU 版本的 PathTracer，然後基於此實現一個 GPU 版本的 PathTracer。PathTracer 相關操作完成之後，實現一個 CPU 版本的 Lightmap Baker，然後基於 CPU 版本，移植到 GPU 上。
首先申明，我這裡只是做概念上的驗證，很多功能都沒有支援，效能也未作任何最佳化，只能算是一個可工作的版本。主要目的，是為了掌握一般性的流程和一些核心的操作。
這裡的 Lightmap 是最基礎的，只支援 diffuse 材質，和兩種解析光源：平行光和點光源。因為這幾個屬性是實時遊戲最常使用到的，也是最容易實現的，所以就只支援了這些。
構造 BVH
以前我寫 ray tracing 相關演算法的時候，都是使用簡單的 sphere 進行場景的構建就能夠滿足要求。但是這個 demo，我們最終是想實現一個 lightmap 的功能，所以需要支援任意的模型。
而 ray tracing 中最核心的一個操作，就是射線和場景的相交性檢測。雖然對於任意的模型，我們可以透過遍歷三角形的方式，依次進行相交性檢測。但是我這裡想嘗試一些最佳化的方案來實現。
這裡主要是透過對模型，建立 BVH 來加速射線與場景的相交性檢測，主要參考 這系列[2]文章來實現。這個系列文章作者講的十分詳細，這裡就不再贅述了。
路徑追蹤 - PathTracer
路徑追蹤（PathTracer）是光線追蹤的一種演算法，它能夠渲染出非常逼真的影像，而且演算法原理十分簡單，所以我使用了它來作為 lightmap baker 的核心演算法。
網上有很多關於 PathTracer 的演算法教程，這裡推薦 Global Illumination and Path Tracing - ScratchAPixel 這篇來做基本的瞭解。我最初的版本就是基於它這個來實現的。
Brute Force 版本
我首先實現了一個最簡單粗暴的版本，流程上大概如下：
- 在一個畫素 quad 中，隨機選擇一個取樣點，作為 camera primary ray 的起始點 S，並且與 eye position 連線，構建 primary ray
- 使用構建出來的 ray 與場景進行相交性檢測
- 獲取到相交點的 position 和 normal 資訊，計算當前畫素的 direct lighting (使用 direct light sampling 的方式)
- 然後根據當前點的 normal 方向，隨機產生一些新的射線 second ray（參考 [3] ），然後對所有這些新的射線，依次繼續以上的 2-4 的操作，直到達到預定的 depth 停止
幾何上，整個追蹤流程類似如下圖所示：
很多關於 PathTracer 的教程裡面，都是使用一個有面積的幾何體作為光源來照亮場景。在這種配置下，射線有可能碰撞到光源幾何體，從而照亮整個場景。但是我這裡因為只支援解析光源 (Point, Directional)，射線永遠無法碰撞到一個點或者平行線，依靠粗暴的碰撞光源幾何體的方式來照亮場景，顯然行不通。所以我這裡使用了 direct light sampling 的方式來實現，即：對於每一個碰撞到的點，都主動的發射一條 shadow ray 到解析光源上（點光源就是向點光源中心點發射，平行光就是以平行光方向進行發射），從而判斷是否被遮擋。如果被遮擋，這個點就沒有 direct light 的貢獻；反之就存在 direct light 的貢獻。透過這種方式就能夠支援解析光源。
如下是這部分的程式碼：
lmb_color li_brute_force(lmb_bake_data& data, lmb_ray ray, int32_t depth, int max_depth, bool ray_as_hit = false) { if (depth >= max_depth) { return lmb_color(0.0f); } lmb_intersection intersection; if (!intersect_scene(data, ray, intersection)) { return data.ambient; } lmb_vec3 hit_pos = ray.o + ray.d * intersection.t; lmb_vec3 hit_normal = intersection.n; // direct lighting lmb_color direct_light = li_direct(data, hit_pos, hit_normal); // indirect lighting lmb_color indirect_light(0.0f); { float pdf = 1.0f / (2.0f * 3.1415926f); for (uint32_t i = 0; i < data.sample_count; i++) { lmb_vec3 sample; sample.x = rand_float(); sample.y = rand_float(); lmb_vec3 sample_dir = map_sample_to_direction(hit_normal, sample); lmb_ray indirect_ray; indirect_ray.o = hit_pos + sample_dir * 0.001f; indirect_ray.d = sample_dir; float ndotl = sqrtf(sample.x); indirect_light = indirect_light + li_brute_force(data, indirect_ray, depth + 1, max_depth) * (ndotl / pdf); } indirect_light = indirect_light * (1.0f / data.sample_count); } lmb_color albedo = data.albedo; return (direct_light + indirect_light) * albedo * (1.0f / 3.1415926f); }
進行直接光照的部分程式碼，如下所示：
lmb_color li_direct(lmb_bake_data& data, lmb_vec3& hit_pos, lmb_vec3& hit_normal) { lmb_color result(0.0f); if (data.light.type == LMB_LIGHT_TYPE_SUN) { float ndotl = max(0.0f, lmb_vec3::dot(hit_normal, data.light.dir_or_pos)); float shadow_atten = 0.0f; if (ndotl > 0.0f) { lmb_ray shadow_ray; shadow_ray.o = hit_pos + hit_normal * 0.001f; shadow_ray.d = data.light.dir_or_pos; lmb_intersection shadow_intersection; shadow_atten = intersect_scene(data, shadow_ray, shadow_intersection) ? 0.0f : 1.0f; } result = data.light.color * (shadow_atten * ndotl); } else if (data.light.type == LMB_LIGHT_TYPE_POINT) { lmb_vec3 dir = data.light.dir_or_pos - hit_pos; float l = dir.length(); dir.normalize(); float dist_atten = 1.0f / (l * l); dist_atten = min(1.0f, dist_atten); float ndotl = max(0.0f, lmb_vec3::dot(hit_normal, dir)); float shadow_atten = 1.0f; if (ndotl > 0.0f) { lmb_ray shadow_ray; shadow_ray.o = hit_pos + hit_normal * 0.001f; shadow_ray.d = dir; lmb_intersection shadow_intersection; bool b = intersect_scene(data, shadow_ray, shadow_intersection); if (b && shadow_intersection.t < l) { shadow_atten = 0.0f; } } result = data.light.color * (dist_atten * shadow_atten * ndotl); } return result; }
這個版本的演算法，簡單粗暴，非常容易理解（除了 蒙特卡洛積分 部分，這部分需要比較複雜的數學知識，可以透過我之前的文章 [4], [5], [6], [7] 來了解，也可以參考圖形學聖經 PBRT[8] 來了解，[3] 中也給出了詳細的解釋）下圖，是此方案下，渲染一張 1024x1024，sample count = 4 的圖的結果：
Single Sum 最佳化
上個章節中的 Brute Force 版本，有個很明顯的問題，隨著遞迴的深度越來越深，需要計算的光線呈指數級爆炸往上漲。雖然它的計算邏輯非常符合直覺，但是計算量也特別的大。而且，很多我們網上看到的教程，都是單條路徑進行追蹤，並不是這樣的一個指數爆炸式光線追蹤的版本。我當初一直很困惑，單條路徑追蹤的版本，怎麼可能和 brute force 等價了？很多教程也沒有從 brute force 版本到單條路徑追蹤版本的理論推導，只是告訴你這樣就是正確的。我找了很多資料來解答我的疑惑，最終找到了這個影片 [9]，它詳細的推導了從 brute force 版本到單條路徑追蹤版本的理論公式，稱之為 single sum。我這裡簡單的列舉下推導的過程。
首先我們給出最基礎的渲染方程公式：
我們將中間的 $L(x, \omega)$ 展開，得到如下：
我們只展開了一次，但是從中可以看到一些規律出來：
公式左右兩邊，隨著中間 $L(x, \omega)$ 展開，相同的模式一個接一個的出現。我們使用蒙特卡洛積分重寫上述公式，得到如下公式：
這個公式就是上面 Brute Force 版本所使用的公式。可以看到每個求和裡面，都還巢狀著另外的求和公式，如此不斷遞迴迭代，從而導致計算非常的慢。所以，透過對公式進行一定的變形，我們能夠找到另外最佳化的方法。我們重寫渲染方程為如下的形式：
這裡就是簡單的將連乘展開，變成連續的加法運算。仔細看下每個加法部分的內容，我們可以得到如下的結論：
從上面公式可以看到，自發光我們不用取樣，直接計算得到；直接光照需要對整個場景進行 N 次取樣；第一次反彈的間接光照需要對場景取樣 NxN 次；第二次反彈的間接光照需要對場景取樣 NxNxN 次；如此遞迴下去。但是很明顯這樣很不合理。因為對於最終的畫面來說，隨著反彈次數越來越高，間接光的貢獻就越來越小。但是上面的公式反應出，隨著反彈的次數變高，我們需要貢獻越來越大的計算量去計算間接光。所以，這裡就是我們該最佳化的地方。
我們可以看到，第一次間接光照的項是一個雙重積分，而第二次間接光照的項是一個三重積分。也就是說，所有間接光照的積分都是一個多重積分。有沒有什麼辦法可以簡化這個多重積分的計算？答案是有的，這裡就展現出了，為什麼我們要選擇蒙特卡羅爾積分求解演算法。蒙特卡羅爾積分演算法它是一個維度無關的演算法，也就是說它不會隨著維度的增加而需要增加更多的取樣數（詳細的講解參考 [8] 中關於蒙特卡羅爾積分演算法的說明）。使用蒙特卡羅爾，我們可以將多重積分的間接光項，轉化為一個高維（就是有很多輸入引數的函式）單重積分函式，從而得到如下新的公式：
然後帶入蒙特卡羅爾積分求和演算法，得到如下：
從上面的公式，我們可以看到，我們可以使用同樣的取樣次數 N 來計算所有的間接光和直接光照項，從而避免了指數爆炸式的計算間接光，卻對最終畫面貢獻不大的情況。因為所有的計算都是使用同樣次數的取樣，我們將相同的操作提取出來，得到如下：
也就是這樣一個計算公式：
這樣，我們就將渲染方程，透過蒙特卡羅爾轉變成了一個只需要單次求和公式計算的公式，也就是 single sum。這個演算法幾何上看起來如下圖所示：
每一次只會對單條路徑進行追蹤，也就是路徑追蹤的由來。
數學上的解釋上面已經完備了，這裡概念性的說一下，如果 brute force 和 single sum 都是產生無窮多條射線，那麼實際所有的光線路徑這兩個演算法都能夠覆蓋掉，也就是說這兩個演算法在最終的結果上是完全等價的，就看收斂到最終結果的速度了。
以下是這個方案所對應的程式碼：
lmb_color
li_single_sum(lmb_bake_data& data, lmb_ray ray, int32_t depth, int max_depth, bool ray_as_hit = false) {
if (depth >= max_depth) {
return lmb_color(0.0f);
}
lmb_intersection intersection;
if (!intersect_scene(data, ray, intersection)) {
return data.ambient;
}
lmb_vec3 hit_pos = ray.o + ray.d * intersection.t;
lmb_vec3 hit_normal = intersection.n;
// direct lighting
lmb_color direct_light = li_direct(data, hit_pos, hit_normal);
// indirect lighting
lmb_color indirect_light(0.0f);
{
float pdf = 1.0f / (2.0f * 3.1415926f);
lmb_vec3 sample;
sample.x = rand_float();
sample.y = rand_float();
lmb_vec3 sample_dir = map_sample_to_direction(hit_normal, sample);
lmb_ray indirect_ray;
indirect_ray.o = hit_pos + sample_dir * 0.001f;
indirect_ray.d = sample_dir;
float ndotl = sqrtf(sample.x);
indirect_light = li_single_sum(data, indirect_ray, depth + 1, max_depth) * (ndotl / pdf);
}
lmb_color albedo = data.albedo;
return (direct_light + indirect_light) * albedo * (1.0f / 3.1415926f);
}這個版本已經能夠比較快速的渲染出一張圖了，欣賞下 1024x1024，sample count = 128 的圖吧：
Russian Roulette
上面兩個演算法，我們都透過限定一個最大追蹤深度來控制遞迴不會無限進行下去。但是較小的深度值，會導致一些高貢獻度的路徑提前終止了；較大的深度值，又浪費了很多計算量。所以最好有個機制，能夠自動根據當前射線所能產生的貢獻度判斷是否需要終止。這裡就引入了 russian roulette 隨機演算法。
我們在追蹤一條路徑的時候，計算一個當前路徑所能夠產生的貢獻度，然後根據這個貢獻度來隨機判斷是否需要終止。隨著路徑不停的往前推進，反彈的光線所產生的貢獻度也越來越小，就越有可能終止它。
因為我們提前終止了路徑的計算，為了保證結果的正確性（即保證演算法是無偏的 unbais），我們需要對結果進行補償，補償的具體操作原理可以看 [9] 中的描述，這裡不再贅述。
如下是支援了 russian roulette 方案的程式碼：
lmb_color
li_russian_roulette(lmb_bake_data& data, lmb_ray ray, int32_t depth, lmb_color throughput, bool ray_as_hit) {
lmb_intersection intersection;
// avoid recursion too deep
if (depth > 64) {
return data.ambient;
}
if (!ray_as_hit && !intersect_scene(data, ray, intersection)) {
return data.ambient;
}
float rr_prob = throughput.max_coefficient();
rr_prob = max(rr_prob, 0.99f); // avoid bounce forever
if (depth < 2) {
rr_prob = 1.0f; // direct and first indirect bounce, always track
}
// russian roulette terminate
if (rand_float() >= rr_prob) {
return data.ambient;
}
lmb_color albedo = data.albedo;
lmb_color brdf = albedo * (1.0f / 3.1415926f);
lmb_vec3 hit_pos = ray.o + ray.d * intersection.t;
lmb_vec3 hit_normal = intersection.n;
if (ray_as_hit) {
hit_pos = ray.o;
hit_normal = ray.d;
}
// direct lighting
lmb_color direct_light = li_direct(data, hit_pos, hit_normal);
// indirect lighting
lmb_color indirect_light(0.0f);
{
float pdf = 1.0f / (2.0f * 3.1415926f);
lmb_vec3 sample;
sample.x = rand_float();
sample.y = rand_float();
lmb_vec3 sample_dir = map_sample_to_direction(hit_normal, sample);
lmb_ray indirect_ray;
indirect_ray.o = hit_pos + sample_dir * 0.001f;
indirect_ray.d = sample_dir;
float ndotl = sqrtf(sample.x);
throughput = throughput * brdf * (ndotl / (pdf * rr_prob));
indirect_light = li_russian_roulette(data, indirect_ray, depth + 1, throughput, false) * (ndotl / (pdf * rr_prob));
}
return (direct_light + indirect_light) * brdf;
}這裡解釋下，因為上面的實現是基於遞迴實現的，而且是隨機終止，理論上可能存在遞迴太深導致棧記憶體用完的情況，所以我加了一個保護避免崩潰。等到後面使用迴圈代替遞迴演算法之後，就可以避免這個操作了。
這個演算法在 single sum 的基礎上又提升了一些速度，這次咱換個燈光效果，老樣子 1024x1024, sample count = 128:
GPU Path Tracer
前面的章節我們已經實現了一個基本的 path tracer，現在我們需要基於這個實現，實現一個 GPU 版本。
從上面的實現可以看出，每個畫素的計算都是獨立的，所以我們能夠基於這個特性，使用 GPU 大量執行緒並行執行的特性來加速整個 path tracer 的計算。在這裡我是基於 vulkan，使用 compute shader 來實現整個 GPU 版本。基本來說，我們就是直接將 CPU 版本的程式碼直接用 GLSL 移植一遍。但是由於 GPU 本身的限制，有一些特殊的地方需要額外注意。
遞迴轉迴圈
基於棧結構（stack）轉迴圈
GPU 上是不支援類似 CPU 的遞迴呼叫的，所以我們需要將遞迴轉化為非遞迴的形式。
一般來說，遞迴我們都可以透過額外新增一個棧結構，將它轉變為一個迴圈。這是因為遞迴函式的呼叫，本質上就是將相關引數壓棧，遞迴函式返回就是出棧，所以完全可以手動使用一個棧結構來模擬，如下是在 CPU 端實現的基於棧的非遞迴版本：
lmb_color
li_russian_roulette_stack_based(lmb_bake_data& data, lmb_ray ray, int32_t depth, lmb_color throughput, bool lightmap) {
static const uint32_t k_stack_size = 64;
lmb_color stack_0[k_stack_size];
lmb_color stack_1[k_stack_size];
uint32_t stack_pos = 0;
lmb_color result = data.ambient;
while (true) {
if (stack_pos >= k_stack_size) {
// too deep, just ignore
break;
}
lmb_intersection intersection;
if (!lightmap && !intersect_scene(data, ray, intersection)) {
break;
}
float rr_prob = throughput.max_coefficient();
rr_prob = max(rr_prob, 0.99f); // avoid bounce forever
if (depth < 2) {
rr_prob = 1.0f; // direct and first indirect bounce, always track
}
// russian roulette terminate
if (rand_float() >= rr_prob) {
break;
}
lmb_color albedo = data.albedo;
lmb_color brdf = albedo * (1.0f / 3.1415926f);
lmb_vec3 hit_pos = ray.o + ray.d * intersection.t;
lmb_vec3 hit_normal = intersection.n;
if (lightmap) {
hit_pos = ray.o;
hit_normal = ray.d;
}
// direct lighting
lmb_color direct_light = li_direct(data, hit_pos, hit_normal);
// indirect lighting
lmb_color indirect_light(0.0f);
{
float pdf = 1.0f / (2.0f * 3.1415926f);
lmb_vec3 sample;
sample.x = rand_float();
sample.y = rand_float();
lmb_vec3 sample_dir = map_sample_to_direction(hit_normal, sample);
lmb_ray indirect_ray;
indirect_ray.o = hit_pos + sample_dir * 0.001f;
indirect_ray.d = sample_dir;
float ndotl = sqrtf(sample.x);
indirect_light = (ndotl / (pdf * rr_prob));
// setup next iteration parameters
ray = indirect_ray;
throughput = throughput * brdf * (ndotl / (pdf * rr_prob));
}
stack_0[stack_pos] = direct_light * brdf;
stack_1[stack_pos] = indirect_light * brdf;
stack_pos++;
depth++;
}
for (int32_t i = stack_pos - 1; i >= 0; i--) {
result = stack_0[i] + stack_1[i] * result;
}
return result;
}Iterative 方法
基於棧結構的方案的確能夠解決問題，但只是從 c++ 的壓棧實現，換成了我們自己的棧實現，理論上還是會存在棧過深的問題。最好能夠將棧結構給完全去除掉。
這裡就介紹一種被稱之為 iterative 的方法，這個方法我是從 IQ 大神的 simple path tracer [10] 中學到的。
首先，我們仔細看下 li_russian_roulette 遞迴版本的程式碼，可以發現整個函式可以定義為如下的形式：
其中，$f$ 表示當前追蹤到的點所對應材質的 brdf；$d$ 表示當前點進行 direct light sampling 計算得到的直接光照；$p\left(\omega\right)$ 表示當前材質的機率密度函式；$t$ 表示使用 russian roulette 下的補償機率。而 $L_i$ 表示入射光線的 irrandiance，它同樣也能夠被這個公式所定義。
我們簡單的展開這個公式得到如下形式：
繼續將公式帶入到 $L_i$ 中，不斷帶入，得到如下一系列的公式：
從這個公式裡面，能看到一些規律，我們可以一次計算一個加法項，而後一個項中很多部分是前一個項中已經計算過的：
根據這個規律，我們可以使用一個變數來儲存中間結果：
初始時，$a = 1.0$。同時需要注意雖然公式裡面寫的是 $r$，但是實際上表示的是上一個加法項中的 $r$，我們可以認為初始的 $r_0 = 1.0$。有了 $a$ 儲存中間的累計操作，我們就可以得到最後的結果為：
$o$ 的初始值為 0.0，透過上面兩個操作就能夠不需要棧結構儲存中間結果的方式來不斷迭代式的計算路徑追蹤。如下是透過這個變換之後得到的程式碼：
lmb_color
li_russian_roulette_iterative(lmb_bake_data& data, lmb_ray ray, int32_t depth, lmb_color throughput, bool lightmap) {
lmb_color accum = lmb_color(1.0f);
lmb_color result = lmb_color(0.0f);
float r = 1.0f;
bool first = true;
while (true) {
lmb_intersection intersection;
if (!(lightmap && first)) {
if (!intersect_scene(data, ray, intersection)) {
break;
}
}
float rr_prob = throughput.max_coefficient();
rr_prob = max(rr_prob, 0.99f); // avoid bounce forever
if (depth < 2) {
rr_prob = 1.0f; // direct and first indirect bounce, always track
}
// russian roulette terminate
if (rand_float() >= rr_prob) {
break;
}
lmb_color albedo = data.albedo;
lmb_color brdf = albedo * (1.0f / 3.1415926f);
lmb_vec3 hit_pos = ray.o + ray.d * intersection.t;
lmb_vec3 hit_normal = intersection.n;
if (lightmap && first) {
hit_pos = ray.o;
hit_normal = ray.d;
}
// direct lighting
lmb_color direct_light = li_direct(data, hit_pos, hit_normal);
accum = accum * brdf * r;
result = result + accum * direct_light;
// indirect lighting for next iteration parameters
lmb_color indirect_light(0.0f);
{
float pdf = 1.0f / (2.0f * 3.1415926f);
lmb_vec3 sample;
sample.x = rand_float();
sample.y = rand_float();
lmb_vec3 sample_dir = map_sample_to_direction(hit_normal, sample);
lmb_ray indirect_ray;
indirect_ray.o = hit_pos + sample_dir * 0.001f;
indirect_ray.d = sample_dir;
float ndotl = sqrtf(sample.x);
indirect_light = (ndotl / (pdf * rr_prob));
ray = indirect_ray;
r = ndotl / (pdf * rr_prob);
throughput = throughput * brdf * r;
}
depth++;
first = false;
}
result = result + accum * r * data.ambient;
return result;
}這裡需要注意，不管是因為路徑追蹤射向了背景，還是因為 russian rolette 提前中止，最後一次迭代的直接光照，我們都是使用 ambient 環境光來表示的，所以程式碼後面有最後一次的補充計算，千萬不能忘記。
漸進式渲染 - Progressive Rendering
即使我們將計算搬到了 GPU 上，需要的計算量依然十分龐大。特別是當 sample count 特別大的時候，如果指望一幀內渲染完畢，那勢必需要等待很久才能夠完成。所以為了能夠更加快速的看到渲染的結果，特別開發了漸進式渲染的方案（Progressive Rendering）。
這裡漸進式的方式，主要是將一個畫素的多個 sample ，分攤到到多幀裡面進行計算。在我這裡就是 1 幀計算一個 sample。也就是說，對於 1024x1024，sample count = 4096 的配置來說，需要分攤到 4096 幀去計算完成。
這裡就存在一個問題，之前我們的計算方法是計算每一個 sample 取樣出來的路徑，並且將結果儲存到一個列表中去。等到全部結果計算完畢，我們得到了一個長度為 4096 的列表，我們將整個列表中的所有結果相加，再除以 sample 數量求平均值，從而得到最終的顏色值。很顯然，這種計算方法在 GPU 上不太友好。
解決這個問題也不是很複雜，我們可以儲存從開始到現在累積的顏色值到另外一張圖上，然後新增一個額外的 pass 從這張圖中獲取值，再除以當前所累計的 sample 數量，得到當前 sample 數量下的平均值。這個方法能夠解決問題，但是需要額外的資料儲存和操作。這裡提供過一個更加簡便的方法來解決。
Welford 方法
上面的問題在其他領域也經常遇到。比如計算一組取樣的方差。這不僅需要當前取樣的平均值，而且需要儲存每個取樣資料才行。但是現實情況下，往往很難確定需要取樣多少次，也就是說我們需要一種流式的處理方式，不管有多少樣本，我都能夠簡單的不用儲存額外資料的方式計算得到，而這個方法就是 Welford 方法。
因為我這裡只是為了計算平均數，所以它的形式非常簡單，我們可以直接手動推匯出來，主要核心操作是使用前一次計算出來的結果。
上圖是前後兩次計算平均數的一般性公式，我們可以將 $a_n$ 用 $a_{n-1}$ 來進行表示：
這樣，新的 $a_n$ 就只依賴與之前一次計算的平均數和當前 sample 計算的顏色 $A_n$ 以及當前累計取樣總數 $n$ 了。透過這個公式，我們就可以不需要額外儲存的情況下，得到當前幀累計以來的平均值了，就可以漸進式的計算所有的 sample。
隨機函式
在 CPU 中，我是簡單的透過 c/c++ 庫自帶的 rand 函式來實現隨機數獲取的。但是 GPU 上，沒有自帶的隨機函式，需要自行定義一個隨機函式。這裡我使用了 [11] 中給出的 GPU 隨機函式，程式碼如下：
uint rng_state;
uint pcg_seed(uvec2 uv) {
return 19u * uv.x + 47u * uv.y + 101u;
}
uint rand_pcg() {
uint state = rng_state;
rng_state = rng_state * 747796405u + 2891336453u;
uint word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
return (word >> 22u) ^ word;
}
float rand_pcg_f() {
return float(rand_pcg()) * (1.0 / 4294967296.0);
}這裡我將 rng_state 設計為一個全域性變數，也就是說，只要在一個 compute shader thread 裡面不停的呼叫隨機函式，就能夠產生一序列分佈均勻的隨機數。同時為了支援漸進式渲染，噹噹前幀的計算結束之後，我們將當前的 rng_state 儲存到貼圖中去，從而可以在下一幀裡面繼續隨機過程，讓幀與幀之間也能夠產生均勻的隨機數序列。
這裡有個很重要的點，每一個隨機數都需要一個起始的種子，因為我的設計是每一個執行緒處理一個畫素，所以我就將當前畫素所在的位置作為起始隨機數種子設定給 rng_state，這樣可以保證畫素與畫素之間也能夠是均勻分佈的隨機數。隨機數設定操作如下所示：
void main() {
uvec2 cur_pixel_loc = uvec2(pushc.tile_pos + gl_GlobalInvocationID.xy);
if (cur_pixel_loc.x >= pushc.width || cur_pixel_loc.y >= pushc.height) return;
uint seed = imageLoad(seedmap, ivec2(cur_pixel_loc.xy)).x;
if (seed == 0) {
rng_state = pcg_seed(cur_pixel_loc.xy);
} else {
rng_state = seed;
}
......
imageStore(seedmap, ivec2(cur_pixel_loc.xy), uvec4(rng_state, 0u, 0u, 0u));
}結果
其他就沒有什麼特別的地方了，無非就是將 bvh 資料打包上傳到 GPU，指定一些渲染的配置引數（fov, light setup, albedo texture setup），然後一比一翻譯 CPU 端的 PathTracer。
雖然已經透過 sample 分攤的方式實現了漸進式渲染，但是對於大尺寸的渲染圖來說，一幀的計算量還是太多了，所以我額外增加了按 tile 進行分幀計算的方式，確保可以保持 60 fps 情況下進行渲染。
透過這樣一系列操作就能夠得到如下幾張圖，都是 1024x1024， sample count > 1024:
Lightmap 烘焙
Lightmap 是渲染中經常使用到的一個技術，隨著時代的演進，lightmap 也具有很多其他的含義。我們今天就只考慮最古老的那種情況，即：lightmap 中以畫素的形式儲存了物體表面最終的 GI 顏色值。我們使用 PathTracer 計算物體表面的 GI 顏色，並且將結果儲存在一張貼圖中便於後續訪問。
老樣子，先在 CPU 端進行基本流程的實現，然後搬到 GPU 上去。
Lightmap UV
第一步，我們需要為場景中的物體產生一個 lightmap uv，我們將使用這個 uv 來對 lightmap 進行取樣，從而得到 GI 光照顏色。由於 GI 光照資訊在物體的表面上都不一樣，一般並不能夠直接使用物體的第一套 uv。因為第一套 uv 可能為了充分利用貼圖，存在 uv 重疊的部分，而這在 lightmap uv 中是不允許的。所以，我們需要將場景物體進行全展開，得到一個合適的 uv 佈局。我這裡是透過 blender 建立了第二套 uv 作為 lightmap uv。使用 blender 建立 lightmap uv 的方式，網上有很多教程，這裡不再贅述。如下就是得到的 lightmap uv 佈局：
實際上，大部分的 lightmap baker 系統都是自己自行產生 lightmap uv，並且也不是像我這樣只處理單個物體，需要考慮很多其他的情況（lightmap 拆分，lightmap uv 合併，縫隙修正等等）。我這裡出於驗證核心烘焙流程的目的，簡化了這些操作，而且想要得到高質量的 lightmap，需要做很多額外的工作，比如文章 [12] 中介紹的各種產生瑕疵的情況就需要額外的操作進行修正。
Lightmap 畫素世界位置計算
有了 uv 之後，接下來我們就需要計算 lightmap 中每一個畫素的顏色。這裡我們就是遍歷每一個畫素，判斷畫素在物體表面的位置，然後從這個位置開始進行 PathTracer，計算當前點的 irrandiance 資訊，並將它儲存在當前貼影像素上。
判斷當前畫素在物體表面哪個位置，我們需要先判斷畫素在物體表面哪個三角形中。由於 lightmap uv 是全展開的，畫素在哪個三角形中是唯一的。因此，我們只要簡單的依次遍歷模型的所有三角形，判斷當前畫素是否在某個三角形中。而判斷一個點是否在某個三角形中，我們可以透過 edge function 來判斷實現，這個原理可以參考 [13] 中關於 edge function 的描述，如下是一個簡單程式碼：
float edge_func(vec2 p, vec2 v0, vec2 v1) {
return (p.x - v0.x) * (v1.y - v0.y) - (p.y - v0.y) * (v1.x - v0.x);
}
bool is_inside_tri(vec2 uv, tri_ex_t tri) {
float e0 = edge_func(uv, vec2(tri.uv1_0.x, tri.uv1_0.y), vec2(tri.uv1_1.x, tri.uv1_1.y));
float e1 = edge_func(uv, vec2(tri.uv1_1.x, tri.uv1_1.y), vec2(tri.uv1_2.x, tri.uv1_2.y));
float e2 = edge_func(uv, vec2(tri.uv1_2.x, tri.uv1_2.y), vec2(tri.uv1_0.x, tri.uv1_0.y));
return ((e0 >= 0.0f) && (e1 >= 0.0f) && (e2 >= 0.0f)) || ((e0 <= 0.0f) && (e1 <= 0.0f) && (e2 <= 0.0f));
}
int find_tri(float x, float y) {
float cur_ux = x / pushc.width;
float cur_uy = y / pushc.height;
for (uint i = 0; i < pushc.tri_count; i++) {
tri_ex_t tri = tri_ex_pool[tri_index_pool[i]];
if (is_inside_tri(vec2(cur_ux, cur_uy), tri)) {
return int(tri_index_pool[i]);
}
}
return -1;
}在知道了當前 lightmap 畫素取樣點所坐落的三角形之後，我們需要計算這個取樣點對應的 position,normal 和 albedo uv 等等資訊。在一個 2D UV 平面上，三角形三個點的 position，uv 和normal 資訊都是已知的情況下，求其中某個點的資訊，我們可以透過三角形重心座標（barycentric coordinate）來計算，相關原理在我最早的 lightmap baker 文章中有講述，參考文章 [1] 中的程式碼進行實現，最終我們能得到如下程式碼：
ray_t build_ray(uint tri_index, float x, float y, inout vec2 uv) {
tri_t tri = tri_pool[tri_index];
tri_ex_t tri_ex = tri_ex_pool[tri_index];
float cx = x / pushc.width;
float cy = y / pushc.height;
float x1 = tri_ex.uv1_0.x;
float y1 = tri_ex.uv1_0.y;
float x2 = tri_ex.uv1_1.x;
float y2 = tri_ex.uv1_1.y;
float x3 = tri_ex.uv1_2.x;
float y3 = tri_ex.uv1_2.y;
float lambda0 = ((y2 - y3) * (cx - x3) + (x3 - x2) * (cy - y3)) / ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
float lambda1 = ((y3 - y1) * (cx - x3) + (x1 - x3) * (cy - y3)) / ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
float lambda2 = 1.0f - lambda0 - lambda1;
// interpolote patch attributes
ray_t r;
r.o = tri.v0 * lambda0 + tri.v1 * lambda1 + tri.v2 * lambda2;
r.d = tri_ex.n0 * lambda0 + tri_ex.n1 * lambda1 + tri_ex.n2 * lambda2;
uv = tri_ex.uv0_0 * lambda0 + tri_ex.uv0_1 * lambda1 + tri_ex.uv0_2 * lambda2;
return r;
}透過上面的操作，我們最終就得到了當前需要追蹤畫素的 position，normal 和 albedo uv 資訊，有了這幾個資訊，我們就可以直接使用前面開發的 PathTracer 進行 GI 顏色的計算，結果計算完畢之後，就將顏色儲存在對應的畫素即可。
搬到 GPU
我這裡也是簡單粗暴的將實現直接搬到 compute shader 中，然後每一個 thread 計算 lightmap 中一個畫素的一個取樣，簡單直接的就能夠得到一個 GPU lightmap baker 了。文章首頁的圖就是使用當前系統烘焙出來的一張 lightmap。這裡再給一張 2048x2048 samplecount = 1024，燈光為 平行光的 lightmap 效果：
總結
至此一個簡單的 lightmap baker 功能就出來了。雖然它效能不好，功能不強，但是對於學習流程來說已經足夠。後面有機會，我會繼續完善這個系統，讓它變的更加實用，高效。
如果發現文中有錯誤之處，歡迎指出，感謝閱讀！
參考
[1] GraphicsLab Project 之光照貼圖
[2] How to build a BVH
[3] Global Illumination and Path Tracing
[4] 圖形學數學基礎之基本蒙特卡洛爾積分
[5] 圖形學數學基礎之重要性取樣
[6] 圖形學資料基礎之1D取樣分佈計算方法 - Inverse Method
[7] 圖形學資料基礎之取樣分佈對映
[8] PBRT Online
[9] Rendering Lecture 04 - Path Tracing Basics
[10] Simple Path Tracer
[11] Hash Functions for GPU Rendering
[12] Baking artifact-free lightmaps on the GPU
[13] Rasterization Practical Implementation - The Rasterization Stage