【RAY TRACING THE REST OF YOUR LIFE 超详解】 光线追踪 3-5 random direction & ONB
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; }
然后我们跑一遍场景
感觉还不是很真实，好像和以前的差不多（正常，这不能算鸽~，因为第八章以后差不多才可以。。）
探索过程也是很值得借鉴的嘛，毕竟这是在一步一步引入新的技术，算是循循善诱~
下面是原书结尾，下一篇我们讲牛逼的技术——直接光源采样
感谢您的阅读，生活愉快~
