Preface
往后看了幾章，對這本書有了新的理解
上一篇，我們第一次嘗試把MC積分運用到了Lambertian材質中，當然，第一次嘗試是失敗的，作者發現它的渲染效果和現實有些出入，所以結尾處聲明要通過實踐，改進當前的效果
於是乎，就有了后面的章節，幾乎整本書都在講，如何一步一步地改進上一篇的畫質，使其更加符合現實，上一篇其實是拋磚引玉
這本書的小標題名為the rest of your life
通過前面幾章，我們可以更好地理解這句話：我們通過MC積分優化效果，采用的是pdf函數，之前說過，這就是一場游戲：尋找一個pdf函數，使得使用它進行重要性采樣得到的渲染圖形更加貼合實際，其實它是沒有止境的，比如pdf是一次曲線、二次曲線、高次曲線、正態分布、高斯分布等等，對應的研究方法也是沒有止境的，比如：你可以通過對光源進行pdf采樣實現最終目的（比如在雙向追蹤中，光源也要發射光線），你也可以通過對不同材質表面的反射狀態進行pdf采樣，進而使得表面顏色變化更光滑更柔和更貼合實際。
上述為個人理解，可能有些出入，吾姑妄言之，汝姑妄聽之，便罷。
Ready
上述說到拋磚引玉，但是好像我們用的不是一張圖，思量再三，還是先把磚整一個，畢竟之后都是圍繞那塊磚評說效果的，另辟蹊徑可能不是明智之舉
所以，我們先把磚搞到手
造磚的代碼：
void Cornell(intersections** scene, camera** cam, rtvar aspect) { intersect ** list = new intersect*[8]; size_t cnt = 0; material * red = new lambertian(new constant_texture(rtvec(0.65, 0.05, 0.05))); material * white = new lambertian(new constant_texture(rtvec(0.73, 0.73, 0.73))); material * green = new lambertian(new constant_texture(rtvec(0.12, 0.45, 0.15))); material * light = new areaLight(new constant_texture(rtvec(15, 15, 15))); list[cnt++] = new flip_normal(new yz_rect(0, 555, 0, 555, 555, green)); list[cnt++] = new yz_rect(0, 555, 0, 555, 0, red); list[cnt++] = new xz_rect(213, 343, 227, 332, 554, light); list[cnt++] = new flip_normal(new xz_rect(0, 555, 0, 555, 555, white)); list[cnt++] = new xz_rect(0, 555, 0, 555, 0, white); list[cnt++] = new flip_normal(new xy_rect(0, 555, 0, 555, 555, white)); list[cnt++] = new translate(new rotate_y(new box(rtvec(), rtvec(165, 165, 165), white), -18), rtvec(130, 0, 65)); list[cnt++] = new translate(new rotate_y(new box(rtvec(), rtvec(165, 330, 165), white), 15), rtvec(265, 0, 295));
*scene = new intersections(list, cnt);
rtvec lookfrom(278, 278, -800); rtvec lookat(278, 278, 0); rtvar dist_to_focus = 10.0; rtvar aperture = 0.; rtvar vfov = 40.0; *cam = new camera(lookfrom, lookat, rtvec(0, 1, 0), vfov, 1, aperture, dist_to_focus, 0., 1.); }
為了清晰點，sample改為了250，圖片為200*200，為了方便重復做實驗，所以參數就這樣吧，能看清即可，"高清大圖"實在是熬不起，都是夜啊，各位見諒~
上一篇結束之后的代碼均不變，只是改一下Cornell 函數
我們得到
這就是上一篇最后得到的效果，沒錯
我們評說一下，這張圖不僅沒有減少噪點，而且高的長方體表面的顏色趨於均勻一致，與實際有偏差
評說好壞當然要有個參考，我們放上兩張第二本書中得到的圖形（未加入MC積分）
上面三張圖，無論你看哪張圖都能發現，正對我們的那個長方體表面是從上到下又黑到白漸變的，而且在正方體上表面高度處對應的長方體表面有一抹正方體上表面反射的白色光，而MC圖形基本上呈均勻色調，甚至可能上面部分還稍白一些
所以我們進入今天的這一篇
Content
今天我們講兩章，第一章是講三維隨機方向向量的生成，這有什么用呢，它給我們的三維空間內進行重要性采樣用，MC需要它！！
Chapter5 Generating Random Direction
在這一章和后續的兩章，我們需要加強我們的理解和我們手中的工具，搞明白什么才是正確的Cornell Box
讓我們首先從如何生成隨機方向開始說起。
為了方便，我們假定z軸為表面法線，之后再轉換坐標系，與此同時，我們規定θ為從法線張開的角度
我們將僅僅處理關於z軸旋轉對稱的分布，所以其他相關的量，均設為均勻分布
給定一個和方向相關的pdf，p(direction) = f(θ)，一維pdf中θ 和 φ如下：
g(φ) = 1/(2π) (均勻分布)
h(θ) = 2πf(θ)sinθ
對於兩個隨機生成的均勻變量r1和r2，我們在第三篇內容中推導出
r1 = ∫0->φ 1/(2π) dθ = φ/(2π)
得 φ = 2πr1
r2 = ∫0->θ 2πf(t)sin(t) dt
t只是一個虛擬的量，用以代替變化的θ，然后輔助實現從0~θ對我們之前的被積函數f(θ)積分（變限積分量不能相同，所以引入t）
我們從下面開始嘗試不同的f()
首先，我們考慮球體上采取均勻密度采樣，因為整個球面為4πsr，故單位球體表面均勻密度采樣的pdf函數為：p(direction) = 1/(4π)，所以得到：
r2 = ∫0->θ sin(t)/2 dt
= （-cos(t)/2）|0->θ
= （1-cosθ）/2
cosθ = 1 - 2r2
一般cosθ更常用一些，就不進一步求θ了，用的時候再acos()
為了在笛卡爾坐標系下生成一個指向（θ，φ）的單位向量，這就涉及到我們的球坐標代換方程了：
x = cosφ sinθ
y = sinφ sinθ
z = cosθ
它對應的球坐標如下
此處r為1
我們開始推導
把關於θ和φ的式子帶入x,y,z中
x = cos(2πr1)*sqrt(1-cos2θ)
y = sin(2πr1)*sqrt(1-cos2θ)
z = 1-2r2
解得：
x = cos(2πr1)*sqrt(r2*(1-r2))
y = sin(2πr1)*sqrt(r2*(1-r2))
z = 1-2r2
然后我們來做個圖，驗證一下它是否如我們所願，生成的點聚攏為單位球面
作者給出的畫圖方法是plot.ly，可以在線畫，但是它要求注冊賬號，流程還挺麻煩的，GitHub授權賬號老是沒響應（可能我網不好），如果有興趣的可以直接到這里
其實matlab最好使了，對不對啊
於是乎，我們先寫入txt,千萬不要寫入xls，還是比較麻煩的
stds ofstream outfile("random_direction.txt"); for (int i = 0; i < 200; ++i) { double r1 = lvgm::rand01(); double r2 = lvgm::rand01(); double x = cos(2 * π * r1) * 2 * sqrt(r2 * (1 - r2)); double y = sin(2 * π * r1) * 2 * sqrt(r2 * (1 - r2)); double z = 1 - 2 * r2; outfile << x << "\t" << y << "\t" << z << stds endl; } outfile.close();
創建一個random_direction.xls的文件，點菜單欄中的數據->自文本，選擇random_direction.txt，一路回車就加載完畢了，可見xls加載文本容易多了
假設你的路徑是
E:\OpenGL\光線追蹤\code\ray tracing 1-3\ray tracing 1-3\random_direction.xls
數據所在表的表名為sheet1
則matlab代碼和效果如下圖
效果還是很好的
我們可以看到，它是均勻隨機的，達到了我們的預期
吶，我們接下來試一下我們第二常用的pdf
即p(direction) = cosθ/π
r2 = ∫0->θ 2π（cos(t)/π）sint dt
= (-cos2t)|0->θ
= 1 - cos2θ
cosθ = sqrt(1-r2)
於是乎我們得到下述的x,y,z
x = cos(2πr1)*sqrt(r2)
y = sin(2πr1)*sqrt(r2)
z = sqrt(1-r2)
下面我們就生成隨機向量
inline rtvec random_cosine_direction() { double r1 = lvgm::rand01(); double r2 = lvgm::rand01(); double z = sqrt(1 - r2); double φ = 2 * π*r1; double x = cos(φ) * 2 * sqrt(r2); double y = sin(φ) * 2 * sqrt(r2); return rtvec(x, y, z); }
我們按照原來的方法把第二個pdf函數產生的隨機情況模擬一下：
那么我們用pdf做一個數值模擬
void estimate2() { int n = 1000000; double sum = 0.; for (int i = 0; i < n; ++i) { rtvec v = ::random_cosine_direction(); sum += pow(v.z(), 3) / (v.z() / π); } stds cout << "π/2 = " << π / 2 << stds endl; stds cout << "Estimate = " << sum / n << stds endl; }
還有一個模擬半球得到的近似，pdf = 1/（2π）
有興趣的可以自己推一下，這里直接給代碼，為了和上面做對比
void estimate() { int n = 1000000; double sum = 0.; for (int i = 0; i < n; ++i) { double r1 = lvgm::rand01(); double r2 = lvgm::rand01(); double x = cos(2 * π * r1) * 2 * sqrt(r2 * (1 - r2)); double y = sin(2 * π * r1) * 2 * sqrt(r2 * (1 - r2)); double z = 1 - r2; sum += z*z*z / (1. / (2.*π)); } stds cout << "π/2 = " << π / 2 << stds endl; stds cout << "Estimate = " << sum / n << stds endl; }
主函數
int main() { //build_1_1(); estimate(); estimate2(); }
結果
可以看到模擬半球的pdf函數效果沒有p() = cos/π 這個好
這一章中所有的都是基於z軸的，而下一章我們真正基於物體表面法線進行
Chapter6 Ortho-normal Bases
ONB是三個相互正交的單位向量集合，笛卡爾坐標系中的xyz軸也是一種ONB
我們這一章的目的就是把上一章基於z軸的隨機方向轉換到基於表面法線的
假設，我們有一個原點o和笛卡爾坐標系向量 x/y/z，當我們描述一個位置
比如：（3，-2,7）時，我們這樣計算：
location = o + 3x - 2y +7z
假如有另外一個坐標系，它的原點為o'，基向量為 u/v/w
我們將要用如下方式表示（u，v，w）這個位置：
location = o' + uu + vv + ww
如果學過計算機圖形學，那么你就可以用矩陣變換去完成坐標系，但是我們這里不需要這種操作。
我們需要的是生成具有相對於表面法線的集合分布的隨機方向。 我們不需要原點，因為向量無起點。
我們引用光線追蹤1-8中所述的相機坐標來計算ONB的三個基向量
我們現在擁有的是法向量n，我們可以把它看做是lookfrom->lookat向量
我們需要一個vup,假設法向量n本身幾乎平行於特定軸，那么vup就取和n垂直的基向量，反之，我們就使用特定軸作為vup
代碼是描述思維最好的方式之一:(假設特定軸為x軸)
if(fabs(n.x())>0.9) //單位法向量n的x分量大於0.9，認為n//x
vup = (0,1,0)
else
vup = (1,0,0)
我們設ONB基向量由 s,t,n 組成
那么 t = cross(vup,n).單位化
s = cross(t,n)
理解可參考上圖,t為圖中的u, s為圖中的v
所以我們的坐標系中的任一點（x,y,z）表示如下
隨機向量 = xs + yt + zn
至此，我們上代碼：
/// onb.hpp // // ----------------------------------------------------- // [author] lv // [begin ] 2019.3 // [brief ] ONB // ----------------------------------------------------- #pragma once namespace rt { class onb { public: onb() { } inline const rtvec& operator[](int index)const { return axis[index]; } inline const rtvec& u()const { return axis[0]; } inline const rtvec& v()const { return axis[1]; } inline const rtvec& w()const { return axis[2]; } public: inline rtvec local(double a, double b, double c)const; inline rtvec local(const rtvec& v)const; void build_from_w(const rtvec&); private: rtvec axis[3]; }; inline rtvec onb::local(double a, double b, double c)const { return a*u() + b*v() + c*w(); } inline rtvec onb::local(const rtvec& v)const { return local(v.x(), v.y(), v.z()); } inline void onb::build_from_w(const rtvec& V) { axis[2] = V.ret_unitization(); rtvec a; if (fabs(w().x()) > 0.9) a = rtvec(0, 1, 0); else a = rtvec(1, 0, 0); axis[1] = lvgm::cross(w(), a).ret_unitization(); axis[0] = lvgm::cross(w(), v()); } } // rt namespace
我們的Lambertian材質的scatter函數需要做相應的改動
原書沒有do-while循環，但是可能會出現pdf為零的除零錯誤，所以我們生成的隨機方向必須能夠算出一個非零的pdf，這並不妨礙我們的算法，因為這一切都是隨機的，隨機錯誤就再隨機一次
bool lambertian::scatter(const ray& rIn, const hitInfo& info, rtvec& alb, ray& scattered, rtvar& pdf)const { onb uvw; uvw.build_from_w(info._n); do { rtvec direction = uvw.local(random_cosine_direction()); scattered = ray(info._p, direction.ret_unitization(), rIn.time()); pdf = dot(uvw.w(), scattered.direction()) / π; } while (pdf == rtvar(0)); alb = _albedo->value(info._u, info._v, info._p); return true; }
然后我們跑一遍場景
感覺還不是很真實，好像和以前的差不多（正常，這不能算鴿~，因為第八章以后差不多才可以。。）
探索過程也是很值得借鑒的嘛，畢竟這是在一步一步引入新的技術，算是循循善誘~
下面是原書結尾，下一篇我們講牛逼的技術——直接光源采樣
感謝您的閱讀，生活愉快~