fopen的路径怎么写_用C++写光线追踪：单根光线的渲染

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;
		//折射率修正。当光线的折射率与球的折射率相同时，表明是球中的射线，要射出球。
		//所以射入介质的折射率是1
		if (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;
		else
			return (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个月。