Monte-Carlo Ray Tracing System (一)原理以及架构设计

实时全局光照

对于现在的CG相关行业来说 就如同一个待跨越的圣杯一样。而在GPU不断进步的过程中，我们却对实现全局光照越来越没有信心 。性能不够也是我们最常在嘴上提到的词语
但是可以说其核心原理早就实现了 在这些年其核心原理并没有取得很大的突破这也是现阶段无法实现实时全局光照的原因之一。

---

对待这样一个系统来说 我们要先从原理说起也就是Monte-Carlo 数值积分。因此第一篇文章着重于算法和设计 具体的代码将不会很多。

---

Chapter （1）

从

∫10x2dx=13

说起
在通常的过程中 对于这样一个简单的定积分我们人脑的求解过程是以固定的微分积分公式计算，公式虽然能方便计算出定积分的精确值，但是有一个局限就是要首先通过不定积分得到被积函数的原函数。有的时候，求原函数是非常困难的，而有的函数，如

f(x)=sinxx

，已经被证明不存在初等原函数，这样，就无法用Newton-Leibniz公式，只能另想办法。

我们以y<(x2)作为判断依据去划分区域来判断区域划分 再通过随机生成的点去重复的做判断

当随机点数量到达一定之后 我们以两个区域的点的数量作为依据 即可得出结果 且点数越多误差越小
根据伯努利大数法则：事件发生的频率依概率收敛于事件的概率p

如下图：

而一般来讲Monte Carlo方法虽然可以解决很多疑难杂症 但是对于复杂度要求极高的光线追踪领域 其结果十分惊艳 但并不能够称之为一个十分简洁的算法

下面我们进入下一章

---

光

我们在其他文章中讲到了光照模型这一概念 然而应用在传统的实时渲染领域 光照模型则是一个相对高效但是低质的概念 所在实时渲染领域 其多数是考虑怎么尽量的使用障眼法去得到一个更逼近数值积分方法的结果 在洪培技术预计算的技术发展方面

类似于UE4的静态效果 Paris demo 在美术和实时光照技术的双重作用下 其效果已经达到了现有计算条件下画质的巅峰

但是相比真正的全局光照，效果仍旧差了一些。

也许会有人说 如此大的代价 去提升那20-30%的画质 值不值得 ？

当然值得！！

在讲解下一章之前我们来了解一下传统的光栅化渲染技术

---

什么是光栅化？什么是光栅？

光栅是光学中一种常见的概念 意为大量等距平行狭缝
然而在实际中光栅化的意义接近于像素化，离散化；用已知概念去理解 类似于透过纱窗去看外面的世界
传统构成三维观察方法常用的为等轴测投影和透视投影。

首先要讨论的是如何把一个三维模型的数据集绘制成相应的模型

这个阶段其实在《Real time Rending》这本书和之前的文章里介绍的很清楚了

渲染管线的基础结构分为：

1.Application应用阶段

需要渲染的几何体从一个固定的数据结构被传递到几何阶段。这传递过去的东西被称为图元，例如点、线、三角形，这些图元有可能最终显示到屏幕上。这是应用阶段最重要的任务。
由软件实现的应用阶段的问题可能在于，它不能像几何和光栅化阶段那样划分为多个子阶段。但是，为了提高性能，这一阶段经常是由多核并行执行。在设计中，这称之为超标量体系，

2.Geometry几何阶段

在这个阶段一般是对模型进行视图变换 对操作进行响应 对顶点进行着色 投影 裁剪和屏幕映射

相应的例如变换过程如下

如裁剪过程

裁剪之后的图像 还不能用来显示 现代操作系统的图像大多数运行在窗口之下 所以针对已经计算好的显示方式来说还需要最后一步 那就是图像帧对屏幕相应区域的映射

过程如下图：

从另一个角度我们也可以考虑到 为什么全屏运行的游戏 效果会更好一些 这与省略了映射阶段也有着一定的关系

3.Rasterizer光栅器

在这个阶段的细分任务中我们算是可以看出这些效果的成因

渲染的过程大致相当于：

至于牵扯到具体的光效 ，阴影 等信息着色器就开始担当大任啦
具体可见：着色器部分

讲完了光栅化体系 我们可以看得到 其效果都是由理论模型模拟效果 效率在现在的光栅化芯片加持下 还算不错 在好的硬件中如GTX1080Ti 甚至能把一个4K 的普遍场景渲染到120Hz

但是

在大多数作品中 我们不可能达到非常好的效果 因为即使类似Unreal Paris demo这种 预渲染的效果 已经计算好了光照 美术和设计还有着色器方面的设计要求实在是太高了

那么 我们实现实时全局光照的意义在哪里？

举个例子 如果在一个场景之下 我们对一个像素需要30条光线来收敛 以最小的720p 也就是1280*720分辨率 再加上60hz的刷新率 可以看到每秒钟需要处理的光线数量达到了16亿5000万条

那我们以一个正常效果的收敛来计算
收敛采样数为200 1080p分辨率 60hz刷新率
每秒钟需要处理的光线数量达到了248亿

不得不承认

在传统的SIMD GPU架构上光线追踪的效率被大大折扣 因此移动GPU巨头 多年为苹果设计GPU的Imagination甚至收购了一些企业 制作了光线追踪加速卡 以至于我们可以在移动端的功耗前提下实现稳定的光线追踪技术
以下为相应的架构图

而我们再看看SIMD

因此实时光线追踪绝不是遥不可及的技术

---

但是仅仅这些我们不足以对光线追踪产生如此大的兴趣 因此实时光线追踪的好处甚至包含了准确反映光的衍射，色散，等等光学特性 这在传统体系之下几乎是无法做到的

最终结合计算物理中的流体模拟，动力学模拟，甚至量子物理模拟 在计算机中建立一个完全拟真的世界 永远是人类为之努力的目标

很庆幸的是 这其中很多技术都取得了关键突破

Chapter（2）

1979年，TurnerWhitted在光线投射的基础上，加入光与物体表面的交互，是光线在物体表面沿着反射，折射以及散射方式上继续传播，直到与光源相交这一方法后来也被称为经典光线跟踪方法、递归式光线追踪（Recursive Ray Tracing）方法，或 Whitted-style 光线跟踪方法。

其主要思想是从视点向成像平面上的像素发射光线，找到与该光线相交的最近物体的交点，如果该点处的表面是散射面，则计算光源直接照射该点产生的颜色；如果该点处表面是镜面或折射面，则继续向反射或折射方向跟踪另一条光线，如此递归下去，直到光线逃逸出场景或达到设定的最大递归深度。
by浅墨

原理：

我们如果能看到一个物体的某个点 这个点必然反射/折射了光线

而对于光线追踪这个过程本身来讲就像是自然世界的逆过程，朝着我们能看到的所有点发射光线，追踪次光线（shadow ,reflection ,refraction）。必然能够回到光源。

图中的几个部分分别说明了 对象的表示 光线的求解方式
相交的求解方式 等等 我们分公式来举例说明一下每个公式怎样用代码求解

基础定义类

三位向量运算类

template<typename T>class Vec3{public:    T x, y, z;    Vec3() : x(T(0)), y(T(0)), z(T(0)) {}    Vec3(T xx) : x(xx), y(xx), z(xx) {}    Vec3(T xx, T yy, T zz) : x(xx), y(yy), z(zz) {}    //正规化    Vec3& normalize()    {        T nor2 = length2();        if (nor2 > 0) {            T invNor = 1 / sqrt(nor2);            x *= invNor, y *= invNor, z *= invNor;        }        return *this;    }    //运算定义    Vec3<T> operator * (const T &f) const { return Vec3<T>(x * f, y * f, z * f); }    //标量*向量    Vec3<T> operator * (const Vec3<T> &v) const { return Vec3<T>(x * v.x, y * v.y, z * v.z); }    //向量1*向量2    T dot(const Vec3<T> &v) const { return x * v.x + y * v.y + z * v.z; }    //点积    Vec3<T> operator - (const Vec3<T> &v) const { return Vec3<T>(x - v.x, y - v.y, z - v.z); }    //向量1-向量2    Vec3<T> operator + (const Vec3<T> &v) const { return Vec3<T>(x + v.x, y + v.y, z + v.z); }    //向量1+向量2    Vec3<T>& operator += (const Vec3<T> &v) { x += v.x, y += v.y, z += v.z; return *this; }    //向量自增    Vec3<T>& operator *= (const Vec3<T> &v) { x *= v.x, y *= v.y, z *= v.z; return *this; }    //向量自乘    Vec3<T> operator - () const { return Vec3<T>(-x, -y, -z); }    //求负    T length2() const { return x * x + y * y + z * z; }    //模^2    T length() const { return sqrt(length2()); }    //模    friend std::ostream & operator << (std::ostream &os, const Vec3<T> &v)    {        os << "[" << v.x << " " << v.y << " " << v.z << "]";        return os;    }};

当然 这是较为精简的写法和部分运算 下篇文章中 我们将着重来尝试优化渲染效率和效果 但是其中涉及的复杂运算不利于基础概念的理解 所以我们用比较基础的运算来说明

同时按照图中顺序 我们将对这整个过程作一个说明

---

Sphere / Ray intersection (给出光线和球的表达式 求相交)

Sphere equation/三维向量的球面表达式 可以想象球面点到球心的差的平方为半径的平方

(p⃗ −c⃗ )⋅(p⃗ −c⃗ )=r2

    Vec3f center;                           /// position of the sphere    float radius, radius2;                  /// sphere radius and radius^2    Vec3f surfaceColor, emissionColor;      /// surface color and emission (light)    float transparency, reflection;         /// surface transparency and reflectivity    Sphere(        const Vec3f &c,        const float &r,        const Vec3f &sc,        const float &refl = 0,        const float &transp = 0,        const Vec3f &ec = 0) :        center(c), radius(r), radius2(r * r), surfaceColor(sc), emissionColor(ec),        transparency(transp), reflection(refl)    { /* empty */ }

Ray equation

r⃗ (t)=o⃗ +td⃗ 

这个应该就不用说了 发射点向量和发射的方向向量

仅需要两个参数

    const Vec3f &rayorig;     const Vec3f &raydir;

Intersection

(1)

(o⃗ +td⃗ −c⃗ )⋅(o⃗ +td⃗ −c⃗ )=r2

or

t2(d⃗ ⋅d⃗ )+2(o⃗ −c⃗ )td⃗ +(o⃗ −c⃗ )⋅(o⃗ −c⃗ )−r2=0

值得注意的是 intersect是定义在结构体当中的 并非独立

而最终判断是否相交需要bool类型做判断 即判断给定光线和给定球体是否相交

    bool intersect(const Vec3f &rayorig, const Vec3f &raydir, float &t0, float &t1) const    {        Vec3f l = center - rayorig;        float tca = l.dot(raydir);        if (tca < 0) return false;        float d2 = l.dot(l) - tca * tca;        if (d2 > radius2) return false;        float thc = sqrt(radius2 - d2);        t0 = tca - thc;        t1 = tca + thc;        return true;    }

当然 这仅仅是与球体的相交检验 更多的 和各种几何体和各种网格的检验代码 我们将在以后的文章中提及

---

Illumination Equation(光照方程)

在上图中 Blin-Phone 光照方程如下

I=kaIa+Ii(kd(L⃗ ⋅N⃗ )+ks(V⃗ ⋅R⃗ )n)+ktIt+krIr

IaKa为递归元素

但是我们作为基本传参的方程样式应该是渲染方程中的

猛地一看这么长确实很懵逼 现在市面上很多的书籍教材都不会对参数做详解
所以就需要我们把这个渲染方程分开来看 看看每一部分到底代表什么

事实上 我们也可以发现 渲染方程都是分开求解的 最后的结果是所有光照类型部分结果的总和

1. （环境光）Iamb=kaIa

Ia是环境光的强度
ka代表表面环境光反射率在0-1之间
但是在基本的渲染方程中环境光被包含在漫反射中

2. （漫反射）Idiff=KdIpcos(θ)=KdIp(N⃗ ⋅L⃗ )+kaIa

Kd为表面漫反射率

Ip为点光源强度

N⃗ 为表面法向量

L⃗ 为入射光方向

#define MAX_RAY_DEPTH 5    Vec3f phit = rayorig + raydir * tnear; // point of intersection    Vec3f nhit = phit - sphere->center; // normal at the intersection point    float bias = 1e-4; // add some bias to the point from which we will be tracing 在追踪中加入偏移float mix(const float &a, const float &b, const float &mix){    return b * mix + a * (1 - mix);}   if ((sphere->transparency > 0 || sphere->reflection > 0) && depth < MAX_RAY_DEPTH) {        //如果 （球体不透明度或者反射率大于0）且深度小于最大递归深度        float facingratio = -raydir.dot(nhit);        // change the mix value to tweak the effect        float fresneleffect = mix(pow(1 - facingratio, 3), 1, 0.1);        // 计算反射方向        // compute reflection direction (not need to normalize because all vectors        // are already normalized)        Vec3f refldir = raydir - nhit * 2 * raydir.dot(nhit);        refldir.normalize();        Vec3f reflection = trace(phit + nhit * bias, refldir, spheres, depth + 1);//递归        Vec3f refraction = 0;        //此处应该包含前面的折射代码          // the result is a mix of reflection and refraction (if the sphere is transparent)          //结果混合折射和反射        surfaceColor = (reflection * fresneleffect + refraction * (1 - fresneleffect) * sphere->transparency)         * sphere->surfaceColor;        }

3. （高光Phong Model）Ispec=KsIpcosn(ϕ)=KsIp(R⃗ ⋅V⃗ )n

Ks代表表面高光反射率

Ip代表之前的点光源强度

n代表高光反射参数 设绝对高光镜面（即反射所有光线）为无穷

以此我们可以看得出 当在交点以法向量做判断求 反射光线如果能返回到光源 则这个区域都是高光区域如下图所示

---

当高光说完之后 我们对整个光照方程应该有了一个较为清晰的认识 在我之前的博客中也讲过Phone光照模型

I=Iamb+Idiff+Ispec=kaIa+Ii(kd(L⃗ ⋅N⃗ )+ks(V⃗ ⋅R⃗ )n)

---

Snell’s Law(折射定律)

sinΘ1sinΘ2=v1v2=n1n2

nairsini=nglasssinΘt

我们看到不同介质的折射率是不同的

        //如果球是透明  计算折射光线        if (sphere->transparency) {            float ior = 1.1, eta = (inside) ? ior : 1 / ior; // are we inside or outside the surface?            float cosi = -nhit.dot(raydir);            float k = 1 - eta * eta * (1 - cosi * cosi);            //1-(eta*sin(i))^2            Vec3f refrdir = raydir * eta + nhit * (eta *  cosi - sqrt(k));            refrdir.normalize();            refraction = trace(phit - nhit * bias, refrdir, spheres, depth + 1);        }        // the result is a mix of reflection and refraction (if the sphere is transparent)        surfaceColor = (            reflection * fresneleffect +            refraction * (1 - fresneleffect) * sphere->transparency) * sphere->surfaceColor;    }        // the result is a mix of reflection and refraction (if the sphere is transparent)        surfaceColor = (reflection * fresneleffect + refraction * (1 - fresneleffect) * sphere->transparency)         * sphere->surfaceColor;        }

当然还有很多可以细化的地方 比如高光的边缘可以有更多的虚化效果 还有不同材质的追踪效果 这些进阶的光线追踪处理方法 我们在下篇文章介绍

---

用以下一个伪代码来说明一下光线追踪的过程

for each pixel of the screen//循环 视平面中的每个像素{    Final color = 0;    Ray = { starting point, direction };// 常见的三维矢量 描述起点的位置 和光线的方向        Repeat      {        for each object in the scene //循环 场景中 物体            {                    determine closest ray object/intersection;                     //判定最近的和光线相交的物体              }        if intersection exists//如果相交存在        {                for each light in the scene//对场景中的每一束光                {                        if the light is not in shadow of another object//如果光不在其他物体的影子内                        {                            add this light contribution to computed color;//添加这束光的采样到已经计算的颜色中                        }                }        }        Final color = Final color + computed color * previous         //  最终的颜色+=已经计算的颜色*之前的颜色        reflection factor;//反射系数        reflection factor = reflection factor * surface reflectionproperty;        //反射系数*=表面反射属性        increment depth;        //深度增量      } until reflection factor is 0 or maximumdepth is reached      //直到反射系数为零或者到达最大值}

有了这段伪代码我们算是能够在一定程度上了解了光线追踪的过程 但是并不直观 也并不易懂 光线追踪的内核

我们用一段代码来帮助我们理解一下

void render(const std::vector<Sphere> &spheres){    unsigned width = 640, height = 320;    Vec3f *image = new Vec3f[width * height], *pixel = image;    float invWidth = 1 / float(width), invHeight = 1 / float(height);    float fov = 30, aspectratio = width / float(height);    float angle = tan(M_PI * 0.5 * fov / 180.);    // Trace rays    for (unsigned y = 0; y < height; ++y) {        for (unsigned x = 0; x < width; ++x, ++pixel) {             //像素遍历            float xx = (2 * ((x + 0.5) * invWidth) - 1) * angle * aspectratio;            float yy = (1 - 2 * ((y + 0.5) * invHeight)) * angle;            //坐标转化            Vec3f raydir(xx, yy, -1);            //方向            raydir.normalize();          //  cout<<" x: "<<x<<" y: "<<y<<endl;            *pixel = trace(Vec3f(0), raydir, spheres, 0);            //递归追踪  递增深度也就是上面伪代码中的        }    }     //完成追踪 缓存写入文件    // Save result to a PPM image (keep these flags if you compile under Windows)    std::ofstream ofs("./untitledHD.ppm", std::ios::out | std::ios::binary);    ofs << "P6\n" << width << " " << height << "\n255\n";    for (unsigned i = 0; i < width * height; ++i) {        ofs << (unsigned char)(std::min(float(1), image[i].x) * 255) <<               (unsigned char)(std::min(float(1), image[i].y) * 255) <<               (unsigned char)(std::min(float(1), image[i].z) * 255);            //    cout<<"Pixel ID"<<i <<"Red : "<< (unsigned char)(std::min(float(1), image[i].x) * 255) <<            //    "Green : "<<(unsigned char)(std::min(float(1), image[i].y) * 255) <<            //    "Blue : "<<(unsigned char)(std::min(float(1), image[i].z) * 255)<<endl;            //输出像素ID和颜色值    }    ofs.close();    delete [] image;}

还有一个可能会引起疑惑的问题
光源定义：

    std::vector<Sphere> spheres;    // position, radius, surface color, reflectivity, transparency, emission color    spheres.push_back(Sphere(Vec3f( 0.0, -10004, -20), 10000, Vec3f(0.20, 0.20, 0.20), 0, 0.0));    spheres.push_back(Sphere(Vec3f( 0.0,      0, -20),     4, Vec3f(1.00, 0.32, 0.36), 1, 0.5));    spheres.push_back(Sphere(Vec3f( 5.0,     -1, -15),     2, Vec3f(0.90, 0.76, 0.46), 1, 0.0));    spheres.push_back(Sphere(Vec3f( 5.0,      0, -25),     3, Vec3f(0.65, 0.77, 0.97), 1, 0.0));    spheres.push_back(Sphere(Vec3f(-5.5,      0, -15),     3, Vec3f(0.90, 0.90, 0.90), 1, 0.0));    // light 最后这项就是光源      spheres.push_back(Sphere(Vec3f( 0.0,     20, -30),     3, Vec3f(0.00, 0.00, 0.00), 0, 0.0, Vec3f(3)));

所以最后判断光线回到光源需要递归的次数来确定相交点的颜色 最后：

在此基础上 基础代码里没有完成的例如

• 对复杂几何网格的追踪，
• 对含有不同材质，纹理的追踪，
• BRDF双向反射分布函数
• 焦散
• 和路径追踪的区别

等等
我们将在下个文章中详细说明。
不仅如此我们后续将继续探讨如何利用深度学习来优化剪枝光线追踪

本人水平有限 又问题请及时指出 邮箱: [e-mail]