0.背景介紹

我依稀記得自己寫過一個“用unity寫光線追蹤”的系列，寫了有幾篇吧，最新一篇的大體內容早已寫完，但始終無法解決網格模型在unity中的讀取問題，故擱置了下去。點數組、面數組與序數總是對不上號。

后來我就去看DirectX12和算法相關了，順帶買了本線性代數的書學習。一晃眼竟然過去了四個月多，前些日子偶然看到了samllpt這個項目，99行實現路徑追蹤，成品相當逼真。我自然是超級感興趣，尋思著99行代碼不消一小時就能敲完。

依照那個代碼臨摹了一遍，學到了很多，更感慨自己當時做光線追蹤的細節仍有欠缺。比如，我當時在球的檢測射線碰撞方法里，為了得到碰撞結果，返回了好幾個數據；而大神的代碼只返回所碰撞球的序號，使得程序簡潔了不少，可讀性也更高了；還有射線與球求交，他是解二次方程的，應用了二次方程的性質，降低了計算量；求折射方向時也是，應用了向量的性質，把結果雜糅到一起，最后的式子比我寫的簡潔的多。不知從哪本書上看到過，程序員進步的最好途徑，就是學習前人的代碼。深以為然。

所以，這篇文章的目的，啊這篇文章就沒什么目的，純粹是記錄以下學習光線追蹤的歷程，會夾雜一些思索和心得。希望以后的我能用得著。

1.單根光線，Whitted-Style

只用一根光線，就能實現折射和反射，相當簡單。誰讓它的來路比較專一呢。

而對應的，漫射的實現是個令人頭疼的問題。它需要采樣半球空間里的所有方向來計算顏色；咱們不可能去采樣完所有方向的，只能用蒙特卡羅法近似；而一旦近似，就不可避免因采樣次數不足而出現的噪點。

簡便起見，光線碰到漫反射物體，做一次檢測，看碰撞點能否受到光照？受到光照即做一次漫射著色：否則就返回黑色。檢測能否受到光照是用的碰撞點與光源位置的向量，與光源的體積無關。所以說最后的結果是個點光源的成像。沒有軟陰影哦。

而當光線碰到有折射、反射等屬性的物體時，那就遞歸起來，誰讓光路可逆呢。閆大佬關于光路可逆的描述很易懂：你看到的物體某一點的顏色，就等于你在這一點看到的光與這一點的屬性作的著色。

遞歸不能無止境對吧，最妙的一點——俄羅斯輪盤賭就是來解決這事的。具體做法：為遞歸方法添加一個深度參數，用來記錄當前是該光線的第幾次作用；當深度大于某一值時，就以概率p停止遞歸，返回黑色；而沒有停止的那1-p概率的光線，將其強度除以1-p。這是非常巧妙的操作，在有限（遞歸次數依概率趨近于0）的次數里，可以保證能量的守恒（雖然是依概率守恒）。對比我之前用的限定深度法，更感到智力的壓制。

2.代碼實現

我仿著smallpt的格式，寫C++代碼。用過C#之后，感覺C++好來勁啊。

目的：使用程序渲染出一張ppm格式的圖片。這種格式比較簡單，易于編程。

第一步，實現一些小東西，為正片做鋪墊。

首先是寫引用的頭文件。與DirectX12相比，引用的不多，就仨：

#include <math.h> #include <stdio.h> #include <stdlib.h>

math.h包含了用到的開平方sqrt()函數、求冪的pow()函數，這些數據計算函數用的非常多；stdio.h包含了寫入文件的fopen_s()和fsprintf()函數，以及輸出調試信息的printf()函數；stdlib.h則提供了生成隨機數的函數，我們在俄羅斯輪盤賭時要用到。

然后，把要用的簡易數據、函數寫一下：

const double _pi = 1 / 3.1415926535898; double erand() {return rand()*1.0 / RAND_MAX; } double Clamp(double x) { return x < 0 ? 0 : x>1 ? 1 : x; } int ToInt(double x) { return (int)(pow(Clamp(x),1/2.2) * 255); }

數據_pi為

，我們在計算漫射著色時要用到。

erand()函數會隨機生成范圍在[0, 1)內的數，注意在分子那乘以1.0來使其化為double類型，不然每次的結果都是0。函數名是我抄samllpt的，能抄到它是我的榮幸。

Clamp()函數會把數x限定在[0, 1]內，ToInt()函數用于把[0, 1]內的數轉化為{0,1,...,255}內的數，誰讓ppm格式存儲的顏色就是這樣呢。其中，對數據進行1/2.2次乘方是gamma校正，把我們計算用的顏色數值映射為視覺正確的顏色。

第二步，數據類型。只用到三種數據類型：Vector3向量，Ray射線，Sphere球。

struct Vector3 {double x, y, z, length2, length;Vector3(double a = 0, double b = 0, double c = 0) :x(a), y(b), z(c) {length2 = a * a + b * b + c * c;length = sqrt(length2);}Vector3 Normalized(){return fabs(length) < 1e-4 ? 0 : *this/length;}Vector3 operator/(double k) { return Vector3(x/k,y/k,z/k);}Vector3 operator*(double k) { return Vector3(x * k, y * k, z * k); }Vector3 operator*=(double k) { return *this = *this*k; }Vector3 operator+=(Vector3 v) { return *this = *this + v; }Vector3 operator+(Vector3 v) { return Vector3(x + v.x, y + v.y, z + v.z); }Vector3 operator-(Vector3 v) { return Vector3(x - v.x, y - v.y, z - v.z); }Vector3 operator*(Vector3 v) { return Vector3(x * v.x, y * v.y, z * v.z); }Vector3 Cross(Vector3 v) { return Vector3(y*v.z-v.y*z,v.x*z-x*v.z,x*v.y-y*v.x); }double Dot(Vector3 v) { return v.x * x + v.y * y + v.z * z; } };

Vector3類型的參數和方法還是挺簡潔的，就單是解析幾何課上老師講的那幾種。只有兩點，一是Normalized()函數中判斷向量長度為0切不可直接==，要判斷其與0的距離小于某一數，這點很重要；二是定義了兩個Vector3類型的乘法，既非點乘也不是叉乘，而是相應位置上的數相乘，用于顏色的處理。

順便說一句，以后用到向量的話，我還是調用DX12庫吧，省的自己造輪子，還有SIMD硬件加速。

struct Ray {Vector3 origin, direct;double n;Ray(Vector3 O, Vector3 D, double n = 1) :origin(O), direct(D.Normalized()), n(n) {};Vector3 GetPosition(double t) { return origin + direct * t; } };

Ray類型由一個起點，一個方向和所在介質的折射率組成。起點和方向好說，GetPosition()函數也就是依距離在射線上取值罷了。所在介質的折射率n，在折射時用的到。

struct Sphere {Vector3 position, albedo;double radius, n, roughness;Sphere(Vector3 p, double r, Vector3 color, double n = INFINITY, double a = 1) :position(p), radius(r), albedo(color), n(n), roughness(a){}double Hit(Ray ray){Vector3 OP = ray.origin - position;double b = OP.Dot(ray.direct);double c = OP.length2 - radius * radius;double delta = b * b - c;if (delta <= 1e-4)return 0;double t = -b - sqrt(delta);if(t<=1e-4) t = -b + sqrt(delta);return t>1e-4?t:0;} };const int ballsNumber = 7; Sphere balls[ballsNumber] = {Sphere(Vector3(0,1,-1),0.2,Vector3(1,1,1)),Sphere(Vector3(0,0,1e3),999,Vector3(0.25,0.75,0.25)),Sphere(Vector3(0,0,0),0.28,Vector3(),1.33),Sphere(Vector3(-0.7,0,0),0.3,Vector3(0.25,0.25,0.75)),Sphere(Vector3(0.7,0,0),0.3,Vector3(0.75,0.25,0.25)),Sphere(Vector3(0,-0.7,0),0.3,Vector3(1,1,0),INFINITY,0),Sphere(Vector3(0,0.7,0),0.3,Vector3(0, 1, 0.4),INFINITY,0) };

Sphere類型肯定要有位置和半徑；再加上Vector3的albedo，即表面反射率（可以理解為顏色），折射率n（這個應該也是Vector3的以模擬色散，因為不同波長的光折射率是不同的），粗糙度roughness。

主要是這個Hit()函數，我照抄的smallpt，非常巧妙。其返回射線與球最近碰撞點的距離，依照這個距離可獲取碰撞點的位置；如果此距離為負，說明球在射線后。思路是解二次方程，方程為：

為射線起點， 為射線方向， 即為射線上距起點 距離的點。 是球心， 是球半徑。方程的意義很明確：射線上距球心的距離是球半徑的點。根據球面的定義——所有距球心為球半徑的點的集合，所以算出來的點是在球面上的。

然后就是解方程的步驟了：

因為

是單位向量，其長度為1，所以，二次項 的系數也就是1。方程化為標準型：

注意到形如

的二次方程，其解可表示為

所以射線與球交點的距離為：

取其最小值即可。

題圖中用到了七個球，除了背景的大球和其前的五個小球，還有一個作光源的球在視野外。這些都是在坐標那設置了。而不透明球的折射率，我設其為無窮大（實際上不可能是無窮大，后續我會加個消光系數，目前先這么搞）。

第三步，光線追蹤。這是本文的核心，也是我們最最了解的部分（真的，大多數人對這部分的了解遠勝于之上的C++代碼，盡管那非常簡單）。

我們的球有四種類型：

光源球：balls[0]；

漫射（材質）球：balls[1]，balls[3]，balls[4]；

透射（材質）球：balls[2]；

反射球：balls[5]，balls[6]；

在程序中，這四種球對光的作用各不相同，光源只發光，漫射球只作光照著色，這兩者都會終止遞歸；反射球把所有光都反射走，而透射球既會反射一部分光，也會折射一部分光，這兩者會遞歸光線。

// 遞歸，返回光線的顏色，以一條射線和一個深度值作參數 Vector3 Radiance(Ray ray, int depth) {// -->1.獲得光碰到的第一個球//定義距離。通過此距離的比較來獲取光線碰到的第一個球double t = INFINITY;//定義碰到球的序號int id = -1;//遍歷整個球數組，來獲取光碰到的球for (int i = 0; i < ballsNumber; i++){//還記得嗎，Sphere.Hit()函數獲得的就是射線到球最近的距離double distance = balls[i].Hit(ray);// 1e-4是誤差范圍，因為計算機嘛，精度有限。//當有球碰到光線，且碰撞點比當前的碰撞點要近時if (distance > 1e-4 && distance < t){//更新值t = distance;id = i;}}//表示沒碰到任何球，返回黑色if (id == -1) return Vector3();//碰到的球本球Sphere ball = balls[id];//return ball.albedo;// -->2.獲取相關數據//碰撞點坐標Vector3 point = ray.GetPosition(t); //碰撞點的法線，先設為向球外的方向，下面會校正。Vector3 pointNormal = (point - ball.position).Normalized(); //碰撞點法線與射線方向的內積。二者應該是相向的，也就是說此值應為負double DdotN = pointNormal.Dot(ray.direct);//如果光線與法線同向，則法線翻轉。//因為要考慮光從球內射出的現象。if (DdotN > 0){pointNormal *= -1;DdotN *= -1;}//反射方向。簡單的幾何學。Vector3 R = ray.direct - pointNormal * 2 * ray.direct.Dot(pointNormal);// -->3.折射//深度值更新，藉此記錄遞歸次數，方便俄羅斯輪盤賭。depth++;//用折射率判斷球是否該折射。if (ball.n<INFINITY) {//輪盤賭。不然遞歸就停不下來了。//獲取一隨機數，其范圍在 [0, 1) 內。q的值可以給定，也可以為隨機數。double p = erand(),q = 0.618;//當遞歸了一定次數后，進行一個概率判斷來終止遞歸。//有q的概率繼續遞歸下去，有1-q的概率直接返回。if (depth>5 && p > q) return Vector3();//這是射入介質的折射率。必須考慮光線從球射出的情形，不能用球的折射率。double bn = ball.n;//折射率修正。當光線的折射率與球的折射率相同時，表明是球中的射線，要射出球。//所以射入介質的折射率是1if (fabs(ray.n - ball.n) < 1e-4) bn = 1.0;//全反射現象，別說你沒學過。if (ray.n * sqrt(1 - DdotN * DdotN) / bn >= 1){//正如其名，全部反射出去，進行下一輪遞歸。return Radiance(Ray(point, R, ray.n), depth) / (depth>5?q:1);}//折射方向，直接用 入射光線方向、碰撞點法線、光線的折射率、射入介質的折射率表示了Vector3 T = (ray.direct*ray.n/bn - pointNormal * (DdotN * ray.n / bn+ sqrt(1 - ray.n*ray.n*(1- DdotN* DdotN)/(bn*bn)))).Normalized();//菲涅爾現象中的直視反射率。時下流行的菲涅爾公式近似式有用到double F0 = pow((bn - 1), 2) / pow((bn + 1), 2);//菲涅爾現象的反射率。double Fr = F0 + (1 - F0) * pow(1 + DdotN, 5);//當深度滿足一定關系時，返回值要除以繼續遞歸的概率q，以作為輪盤賭的補償。if(depth<5)return Radiance(Ray(point, T, bn), depth)*(1-Fr) + Radiance(Ray(point, R, ray.n), depth)* Fr;elsereturn (Radiance(Ray(point, T, bn), depth) * (1 - Fr) + Radiance(Ray(point, R, ray.n), depth) * Fr)/q;}// -->4.反射。這部分相當簡單。//這就有了反射，判斷粗糙度。這里的條件有些嚴苛，是為了反射全部光線。if (ball.roughness < 1e-4) {//還是輪盤賭，終止遞歸用。double p = erand(), q = 0.618;if (depth > 5 && p > q)return Vector3();//跟上面一樣，要補償輪盤賭的幸存者。if(depth<5) return Radiance(Ray(point, R, ray.n),depth);else return Radiance(Ray(point, R, ray.n), depth)/q;}// -->5.漫射。//獲取個光源方向先。Ray lightRay(point, balls[0].position - point);//看看從光源方向過去，能不能直達光源。t = INFINITY;id = -1;for (int i = 0; i < ballsNumber; i++){double distance = balls[i].Hit(lightRay);if (distance > 1e-4 && distance < t){t = distance;id = i;}}//因為光源的序號是0嘛，所以如果碰到的最近的球不是光源，那就返回黑絲。if (id != 0) return Vector3();//這樣得到有受到光源照射的部分，作著色計算即可，不再遞歸。return ball.albedo * balls[0].albedo * pointNormal.Dot(lightRay.direct) * _pi; }

這里的重中之重是折射方向的計算。不必費勁行數去寫什么正弦余弦，只需一個稍長的式子，折射方向就算出來了。推導過程如下：

從折射定律開始，

其中

表示折射率（可不是法線，法線是 ）， 是與法線的夾角（以上圖為準），下標 和 的作用是區分當前介質和目標介質。

我們知道光線方向

和法線方向 ，可知 ，也就是

根據定律馬上就有

我們可以把折射向量

表示為法線向量和另一個向量 的線性組合，這一個方向和光線方向、法線方向在同一平面內，且與法線方向垂直，與光線方向的內積為正。

不管是幾何直覺，還是施密特正交化，總之可得此向量為

歸一化后其為

所以在

方向上 的分量為

正好把煩人的根號消去，很妙吧？

而在

方向上 的分量為

把這倆分量加到一起就組成了

：

好，就是這樣。寫成代碼就是那個式子了。

Vector3 T = (ray.direct*ray.n/bn - pointNormal * (DdotN * ray.n / bn+ sqrt(1 - ray.n*ray.n*(1- DdotN* DdotN)/(bn*bn)))).Normalized();

代碼里只是把第一項里的

拆開放到了第二項里罷了。最后調用歸一化方法，真的只是為了保險，萬一它就不是單位向量了呢。

其他就沒啥了，代碼里的注釋我可認真寫了的。

第四步，文件寫入。上面的東西可以叫做理科，現在就要開始工科的內容了。

int w = 1024, h = 1024; double stepAngle = 3.1415926535898 / (3*h); Vector3 View(double i, double j) {i = i - w / 2;j = j - h / 2;return Vector3(cos(stepAngle * j) * sin(stepAngle * i),sin(stepAngle* j),cos(stepAngle * j)*cos(stepAngle * i)); }

這些文件的參數我寫成了全局變量，因為有兩個函數都要用到它，其中第一個就是View()函數。View()函數接收像素的坐標，返回對應這個像素的射線方向。為什么參數是double類型的呢？要用到子像素來抗鋸齒，所以會有不是整數的像素坐標。

視野角度是縱向的角度，取60°也就是

，比較常見的視野角度了。每個像素的視角就是視野角度除以縱向像素數。

在View()函數中，取(w/2, h/2)的像素點為視野中心，用球面坐標系取得相應方向，僅此而已。

int main() {FILE *f;fopen_s(&f, "D:/ij.ppm", "w");fprintf(f, "P3n%d %dn255n", w, h);for (int j = h; j > 0; j--){for (int i = 0; i < w; i++){Vector3 color(0,0,0);color += Radiance(Ray(Vector3(0, 0, -1), View(i + 0.25, j)), 0);color += Radiance(Ray(Vector3(0, 0, -1), View(i, j + 0.25)), 0);color += Radiance(Ray(Vector3(0, 0, -1), View(i - 0.25, j)), 0);color += Radiance(Ray(Vector3(0, 0, -1), View(i, j - 0.25)), 0);color = color / 4;fprintf(f, "%d %d %d ", ToInt(color.x), ToInt(color.y), ToInt(color.z));}printf("已渲染%f...n", j * 1.0 / h);}printf("--渲染完畢--n"); }

打開文件、寫入文件倒沒啥說頭，任何一篇C++教程都比我說的好。

文件開頭寫的那個

fprintf(f, "P3n%d %dn255n", w, h);

是ppm的格式，好像也有其他的參數？不知道了。

值得注意的是這里對每個像素進行了四次采樣，取的是每個像素范圍內，上下左右四點，這樣的采樣不是最佳的；smallpt里的帳篷濾波器采樣是最佳的，以后再實現吧（但只是要用的話，把它的函數搬過來就行了）。

3.一些成品圖

改變球的位置、光源位置以及其他的一些參數，就能獲得不一樣的圖片，挺好玩的。

為隨機選取的球添加了反射效果。占據視野大半的大球充分放大了其正后方的灰色小球，而灰色小球周圍的球的圖像被折射壓縮到球的邊緣。

為隨機選取的球添加了不同折射率的折射效果。可以看到有的球折射率很接近1，看著像透明；而有的球折射率頗高，扭曲了其后的圖像。

隨機的反射+隨機的折射。若是有全局光照的話，效果會更上一層樓。

看著，還行。漫射，以后再說。希望不會隔4個月。

總結

以上是生活随笔為你收集整理的fopen的路径怎么写_用C++写光线追踪：单根光线的渲染的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。