GPU光线跟踪算法加速结构研究
摘要：基于GPU的光线跟踪算法是当前图形学研究的一个热点，也是将来用于广告、电影、游戏等娱乐产业的关键技术。本文论述了如何对基于GPU的光线跟踪算法进行实现，以及利用各种加速结构，加速算法实现，提高算法执行效率，并对各种加速结构的效果进行了比较研究。
关键词：GPGPU 光线跟踪 BVH KD-Tree
1．引言
近年来，CPU无论在运算能力，还是在可编程性上都得到了大幅的提高，GPU已经在需要大量运算的密集运算领域发挥了举足轻重的作用。各种基于CPU的密集运算被移植到GPU上，以利用GPU巨大的运算能力，加速整个算法的运算过程。光线跟踪算法是生成真实感图形的一种非常重要的方法，在电影、游戏、广告等产业，获得广泛的应用，而光线跟踪算法也是典型的密集运算算法，利用原始的基于CPU的光线跟踪渲染一幅图片是非常耗时的操作。因此，如果能够将CPU上的光线跟踪算法，映射到CPU上，加速光线跟踪算法的执行时间，将会带来巨大的经济效益。因此，基于CPU的光线跟踪算法已成为国内外科研人员的研究热点。
2．基于GPU的光线跟踪
2.1 相关工作
当前，主要由两种方法利用CPU来加速光线跟踪算法。第一种是Carr等人提出来的，将CPU转换为一个蛮力的执行光线一三角形求交的计算器，而将任何的光线生成以及着色过程在CPU上完成。这就需要CPU依然执行绝大部分的渲染工作。C arr等人指出，在ATI Radeon 8500上，每秒最快能够执行1亿2千万次的光线一三角形求交。同时，作者也指出，由于GPU的单精度浮点的限制，图片上依然存在一些不太真实的地方。
第二种方法由Purcell等人提出的，改种方法将整个光线跟踪器都移植到CPU上进行实现。从光线的产生，加速结构的遍历，到最后的着色过程都在GPU上执行。此后，有很多相同的项目都是基于Purcell的模型上进行的。
2.2 GPU上的光线跟踪算法的映射方式
将传统的CPU上执行的光线跟踪算法，映射成为一个GPU协助的，或者基于GPU的光线跟踪器有众多方法。下面重点介绍Purcell提出的映射模型，以及在本文的实现中提出的一个基于CPU的Whitted模型的光线跟踪器。该光线跟踪器的布局如图2.1所示：
在Purcell的论文中，它将光线一三角形求交，以及遍历过程分离成两个独立的遍历内核和求交内核。本文的实现中，也按照上述模型图，将光线跟踪算法分解成光线生成，光线一三角形求交，着色这三个步骤。
在对光线进行跟踪之前，需要生成从视点指向屏幕的原始光线( primary ray)。在一个GPU上，能够使用光栅器的插值的能力，在一个单一的内核调用中，产生所有的原始光线。
给定观察矩形（被采样用于产生图片的投影平面的一部分）的四个角，以及视点，首先计算出这个视锥体的四条边线。如果让光栅器在这4条光线之间，按照512×512规格，在这四条光线之间按照方向进行插值，最终就可以获得能够产生一幅512×512图片（一个像素一个采样点）的所有原始光线的方向。同时能够将这些方向存储在一个纹理里，并把它作为求交内核的输入。所有的原始光线具有相同的起始点，但是仍然把它存储在一个同方向纹理具有相同维度的纹理内。因为当生成阴影光线或者反射光线的时候，光线的原点会发生改变。
求交内核把光线的原点，方向，以及场景的描述作为输入数据。在内核被调用数次之后，我们对于每一个像素输出一个击中记录。如果一条光线击中了场景中的某个三角形，返回击中点的3个重心坐标，以及相关的被击中的三角形。此外，还将输出被发现的交点沿光线的距离，以及被击中三角形的材质。这就需要使用5个浮点数值组成一个击中记录。纹理只能够支持4个颜色通道( RCBA)，所以，如果能把击中记录裁减到4个值，那么将是非常有益的。
观察发现，只需要3个重心坐标的两个，因为在三角形内部，它们相加的和总是1。这就使得在一个单独的RGBA纹理中存储交点记录是可行的，并且它的维度同其它两个光线纹理的维度相同。
Moller和Trumbore提出了一个高效的光线一三角形求交算法，使用这个算法，并利用CPU在向量计算上的优势来进行求交计算。下面列出了求交的代码，这个代码也展示了如何利用向量指令来提高效率。
当所有的原始光线都已经计算出了相交的状态的时候，就能够查询着色过程所需要的表面法线和材质的信息。每一个击中记录都存储了一个指向材质纹理的索引，这个材质纹理包含了三角形的法线，材质颜色以及类型。三个顶点的法线根据击中记录的中心坐标进行了插值。最终的颜色能够按(N-L)C进行计算，此处Ⅳ是法线，L是光源的方向，G是三角形的颜色。
现在根据击中的三角形所具有的材质的类型（漫反射材质，或者镜面反射材质），需要产生二次光线，以此来计算阴影和反射。
1)如果一条光线射出场景之外，像素就被赋予全局的背景颜色。
2)如果一条光线击中了一个漫反射材质表面，就发射一条阴影射线( shdow ray)。这些光线的起始点在击中点，方向为从击中点指向光源。
3)如果一条光线击中了一个镜面反射材质表面。就发射一条镜面反射光线。镜面发射光线的起始点也在击中点，但是它的方向是在击中点处关于入射光线和插值后的法线对称的方向。一个真正的Whitted类型的光线跟踪器也支持透明材质，从而能够产生折射光线。但由于主要是研究加速结构，所以在本文的实现中，没有考虑折射光线。
4)如果阴影光线击中了某个几何体，这就说明在光源和击中点之间，存在某个几何体，所以这个像素就应该是黑色（处于阴影中）。当跟踪阴影光线的时候，不关心最近的那个击中点，更加关心的是是否存在这样的击中点。因此，当有一个交点被发现，就可以停止整个求交过程，从而加速算法的处理过程。在本文的实现中，以相同的方式跟踪阴影光线和反射光线，因此，就没有使用到这个优化策略。
已经对每一个像素产生了正确二次光线，如果需要，就能够执行另外一趟遍历/求交过程，对上述的二次光线进行跟踪。每一次调用着色程序就能够对每一个像素返回一个颜色值和一条新的光线。着色内核也可以将前一次着色程序的输出当作本次着色程序的输入。这就使得能够在跟踪连续的光线的时候合并这些连续的镜面反射的颜色。
同Carr等人的程序不同，本文所采用的程序不存在浮点精度太低的问题，因为Ceforce 7300在整个管线中支持真正的32位浮点操作。
3．加速结构的实现和比较
3.1均匀栅格
均匀栅格是第一个在GPU上实现的加速结构。Purcell给出了很多选择均匀栅格作为加速结构的理由，但是Purcell没有详细的说明为什么均匀网格对于硬件实现而言比其它的加速结构要更加的简单。当在探讨了均匀栅格的一些主要特性的时候，更加清晰的知道了均匀栅格为什么会成为一个好的GPU机速结构。
首先，只用使用简单的算术运算，就能够对于每个体素的遍历在常量时间能被定位和存取。这就消除了对树的遍历的需要，以及重复的纹理查找工作，而纹理查找是相当耗时的。
其次，体素的遍历是通过递增算术运算来完成的。这就消除了对堆栈的需要，使得我们能够从光线的起始点开始，以距离递增的顺序访问体素成为可能。
再其次，由于对于体素的访问是沿着光线，以距离递增的方式遍历的，所以，一旦在一个被访问的体素中报道发现有一个交点，就可以停止这条光线对体素的遍历过程，从而提高整个遍历过程的速度。
最后，用于遍历的代码非常适合用向量编写，而向量形式的编码风格又非常适合GPU的指令集。
然而，均匀栅格的缺点就是由于它是空间细分结构的一种特殊情况，多个体素可能包含相同三角形的多个引用。由于无法使用mailbox技术，这就意味着需要对于相同的光线和三角形之间进行不止一次的相交测试。
3.2 KD-tree
最近，Havran等人对基于CPU的光线跟踪算法的加速结构进行了比较，得出的结论是对于众多不同类型的测试场景，平均而言，KD-tree是最快的。所以，有必要考察一下对于基于KD-tree的GPU光线跟踪算法，是否也会有相似的结论。
就像均匀栅格一样，KD-tree也是一种空间细分结构。同均匀网格不同的是，KD-tree利用一个二叉树将场景表示成一个层次结构。
在二叉树中，我们将内部节点和叶子节点区分开。叶子节点用来表示体素和与之相关的保存在该体素内的三角形的引用。一个内部节点用来表示空间区域的某个部分。所以，内部节点包含一个分裂面的两个子树的引用，而叶子节点只包含一个三角形列表。
KD-tree的创建过程从上而下，根据一个评价函数，通过放置一个分离平面，递归的将场景分离成两个体素。我们能够以递归的方式遍历KD-tree，但是由于GPU没有堆栈结构，所以无法应用递归的策略。取而代之的是，我们能够通过记住我们沿着光线前进了多远来向上或者向下遍历树。这种策略消除了需要堆栈的限制，使得用CPU来完成对KD-tree结构的遍历成为可能。
当使用GPU对KD-tree进行遍历的时候，KD-tree像均匀栅格那样被表示成一个纹理的集合。这就意味着有一个保存树数据的纹理，一个保存三角形列表的纹理，和一个保存实际的三角形数据的纹理。GPU的遍历首先调用一个初始化内核，然后按照需要，多次调用合并后的遍历和求交内核。
3.3 包围体层次(BVH)
给定一些随机的光线，通过计算遍历包围体层次的平均花费，就可以测量出该包围体层次的质量。迄今为止，还没有构建最优的包围体层次的算法，也就是说，如何准确的测量一个包围体层次的平均遍历时间还不是很明显。
Goldsmith和Salmon提出了一个评价函数，通常被称为表面积启发式函数。他们通过父节点和孩子节点的表面积之比来形式化的表述这个关系，此评价函数如下所示：
此处，hit(n)是光线击中节点n的情况，Sn是节点n的表面积，c和p分别表示父节点和孩子节点。
这个评价函数给出了，当用一条随机的光线同层次结构求交的时候，成本上的估计。由于没有最优的方法去有效的构造一个最优的BVH，提出了不同的构造技巧。下面，将列出比较通用的方法。
在实践中，对于包围体应用的最广泛的就是轴对齐包围盒(AABB)。AABB易于实现，并且同光线的求交测试非常快。大多数有关BVH的论文在描述BVH的创建的时候，通常分别以Kay和Kajiya，或者Goldsmith和Salmon这两种基本的想法为基础。Kay和Kajiaya建议以自上而下递归的方式进行BVH的创建。
Goldsmith和Salmon提出了一个更加复杂的自底向上的构造方式。Goldsmith和Salmon指出，BVH的质量同作为输入传人的三角形的顺序有关。因此，他们建议在构造BVH之前，随机打乱三角形的顺序。下述算法就是利用Kay/Kajiya的思想创建某个场景的包围体层次的方法：
4．结束语
本文成功的在GPU上实现了用于光线跟踪算法中的各种加速结构，并对这些加速结构在GPU上的加速效果进行了比较。均匀栅格作为第一个在CPU上实现的光线跟踪器的加速结构，也被证明是最慢的，除非是只包含一个单独的物体的场景的情况。均匀栅格不适合几何体的密度非常高的场景。另外，对于均匀栅格的CPU上的遍历表示，也需要大量的数据。Foley和Sugerman认为，对于大多数场景，KD-tree的效率要比均匀栅格高。但是，在KD-tree的遍历过程中，无论是重置阶段还是回退阶段，片元程序都非常的复杂，但这种复杂性也使得其能够在场景的几何体的密度改变的时候做出适当的调整。本文实现的BVH被证明在加速效果上要超过均匀栅格和KD-tree，在现阶段，BVH是在GPU上实现的最快的加速结构。并且在GPU上实现BVH加速结构要比实现其他加速结构更加的简单。
参考文献：
[1]Randima Femado编，姚勇，王小琴译．GPU精粹一实时图形编程的技术，技巧和技艺[M]．北京：人民邮电出版社，2006.
[2] Matt Pharr编著，龚敏敏译.GPU精粹2-高性能图形芯片和通用计算编程技巧[M].北京：清华大学出版社．
[3]昊恩华，柳有权．基于图形处理器(GPU)的通用计算叨．计算机辅助设计与图形学学报，2004，16(5): 601-[4] Philip J.Schneider，David H.Eberly著，周长发译，计算机图形学几何工具算法详解[M].北京：电子工业出版社，2005.
[5] Martin Christen. Implementing ray tracing on GPU. Master´sthesis, University of Applied Sciences Basel,
地层面的拟合是多源地质建模中最为重要的步骤。无论是通过Delaunay细分方法增加节点还是通过网格分级增加的节点，都需要进一步求取其高程值。因此，必须借助插值方法来对有限的数据点信息所形成的初始地层面进行精确、光滑处理。插值是指在根据已知的数据计算未知值的过程，结果是形成一个连续分布的数据场(Spragueetal.，2005)。目前存在的插值方法众多，但适合三维数据场表达的插值方法还是比较有限的(胡小红等，2007)，常用的有距离加权反比法(Inverse-Distance Weighting，IDW)和普通Kriging法，其中距离加权反比法属于一种确定性差值方法，Lu et al.(2008)对该方法进行扩展，根据样本的数量和分布密度特征，使其能够根据样本的特征来确定其参数的取值；而Kriging法属于一种不确定性差值方法，Jessell(2001)对其进行了深入研究，基于该方法提出一种势场(potential-field)的插值方法，能够处理存在断层的不连续数据场。这里，仅对这两种方法进行讨论。
5.3.3.1 距离加权反比法
距离反比加权法是最常用的地质数据插值方法之一。它首先由气象学家及地质学工作者提出，后由D.Shepard进行改进，故该方法被称为Shepard方法。
距离反比加权法的基本思想是将插值函数f(P)定义为各数据点函数fk的加权平均，它认为与待插值点距离最近的若干个已知采样点对待插值点的值贡献最大，其贡献与距离的某次幂成反比。
距离反比加权法的基本原理可用下式表示：
数字地下空间与工程三维地质建模及应用研究
式中：f(P)是待插值点P的估计值；fi是第i(i=1，…，n)个已知采样点Pi的样本值；di是第i个样本点Pi与待插值点P的广义距离；v是距离的幂，它显著影响内插的结果，它的选择标准是最小平均绝对误差。相关研究结果表明，幂越高，内插结果越具有平滑的效果(Lu et al.，2008)。在幂指数为2时，不仅能得出较满意的内插结果，而且具有容易计算的优点。在实际应用中，一般用距离平方反比法来求待估算值。
传统的距离反比加权法是非常简单和自然的，但应用于实际的地质特征插值时，却存在以下明显的缺陷：
(1)当数据点的数目非常庞大时，f(P)的计算量将变得十分巨大，计算量的庞大甚至可能导致该方法变得无法实现。
(2)该方法只考虑了从Pi到P的距离，而没有考虑其方向。事实上，只考虑距离的大小是不充分的，有的已知离散点虽然距离待插值点较近，但它对待插值点的影响可能会被其他点屏蔽掉。
(3)在已知采样点Pi的邻域内，由于di≈0计算的误差将变得非常敏感，尤其是当两个项形式占优而又符号相反时，计算的误差更是如此。
因此，在实际应用时，必须对传统的距离反比加权法加以改进。考虑到地质特征数据的空间相关性，对距离反比加权法可附加地质体结构和影响距离两个限制条件，以提高空间几何或属性数据插值的合理性和精度。具体改进如下：
(1)地质体几何结构限制条件。在空间特征插值过程中，在选择影响待估计值的原始样本数据时只选取同一地层岩性地质内的样本；不是同一地层岩性地质体内的样本，即便距离很近，也不采用。即按照地质构造划分空间单元，并只在同一地层岩性地质体内选取样本点。
(2)邻近样本点的选择条件。选择待插值点的邻近点时，可考虑三个原则：一是距离原则，即根据地质属性数据的特征给出一距离r，在该距离之外的样本点对待插值点的估算无影响；二是点数原则，即给定一数据m，以距离待插值点最近m个样本点进行估算；三是利用Voronoi图求取待插值点的邻近点。
(3)建立数据点索引表，提高待插值点周围样本点的搜索效率，从而大大减少大数据量队的计算量。
5.3.3.2 普通Kriging插值方法
设研究区域为A，区域化变量(即欲研究的物理属性变量)为 {Z(x)∈A}，x表示空间位置。Z(x)在采样点xi(i=1，2，…，n)处的属性值(或称为区域化变量的一次实现)为Z(xi)(i=1，2，…，n)，则根据普通Kriging插值原理，未采样点x0处的属性值Z(x0)估计值是n个已知采样点属性值的加权和，即：
数字地下空间与工程三维地质建模及应用研究
其中，λi(i=1，2，…，n)为待求权系数。
假设区域化变量Z(x)在整个研究区域内满足二阶平稳假设：
(1)Z(x)的数学期望存在且等于常数：E[Z(x)]=m(常数)。
(2)Z(x)的协方差Cov(xi，xj)存在且只与两点之间的相对位置有关。或满足本征假设：
(3)E[Z(xi)-Z(xj)]=0。
(4)增量的方差存在且平稳：Var[Z(xi)-Z(xj)]=E[Z(xi)-Z(xj)]2。
依据无偏性要求：E[Z*(x0)]=E[Z(x0)]。
推导可得： 。
在无偏条件下使估计方差达到最小，即：
min{ Var[Z*(x0)-2μ( (λi-1))}，其中μ为拉格朗日乘子。
可求得求解权系数λii(=1，2，…，n)的方程组：
数字地下空间与工程三维地质建模及应用研究
求出诸权系数λii(=1，2，…，n)后，就可以求出采样点x0处的属性值Z*(x0)。
上述求解权系数λii(=1，2，…，n)的方程组中协方差Cov(xi，xj)若用变异函数γ(xi，xj)表示时，形式为：
数字地下空间与工程三维地质建模及应用研究
变异函数的定义为：
数字地下空间与工程三维地质建模及应用研究
由Kringing插值所得到的方差为：
数字地下空间与工程三维地质建模及应用研究
或
数字地下空间与工程三维地质建模及应用研究
A.Kringing插值法中的变异函数
变异函数是Kringing插值法插值的基础。插值中需要首先确定所研究的区域化变量的变异函数。假设研究的区域为A，区域A中有一区域化变量Z(x)，它在位置xi(i=1，2，…，N)上的一次采样为Z(xi)(i=1，2，…，N)，则Z(x)的变异函数的定义为：
数字地下空间与工程三维地质建模及应用研究
一个空间变量的空间变异性是指这个变量在空间中如何随着位置的不同而变化的性质。变异函数通过其自身的结构及其各项参数从不同的角度反映空间变异性，确定变异函数的过程就是一个对空间变异性进行结构分析的过程。
设h是一个模为r=|h|，方向为a的向量，如果存在着被向量h所隔开的Nh对观测数据点，则在a方向上相应于向量h的实验变异函数γ*(h)可表示为如下形式：
数字地下空间与工程三维地质建模及应用研究
其中，z(xi+h)和z(xi)分别位于点xi+h和xi(i=1，2，…，Nh)上的观测数据。
B.变异函数理论模型
当获取实验变异函数值后，需要先选择变异函数理论模型，然后对所选择的变异函数理论模型进行参数拟合，这一过程被称为“结构分析”。
变异函数理论模型参数一般包括：变程(range，一般用a表示)、基台(sill，一般用C(0)表示)、拱高(一般用C表示)、块金常数(Nugget，一般用C0表示)，如图5.9所示。
图5.9 理论变差函数
变程a表示了从空间相关性状态(|h|＜a)向不存在相关性状态(|h|＞a)转变的分界线；变异函数在原点处的间断性称为块金效应，相应的常数C0= (h)称为块金常数；基台C(0)具有协方差函数：C(h)=C(0)-γ(h)的二阶平稳区域化变量Z(x)的先验方差：Var{Z(x)}=C(0)=γ(h)；拱高C为变异函数中基台C(0)与块金常数C0之差：C=C(0)-C0。
C.变异函数理论模型分类
变异函数理论模型一般分为有基台值和无基台值两大类。有基台值的变异函数理论模型包括球状模型、指数模型、高斯模型等(图5.10)。最常用的是球状模型。球状模型公式为：
数字地下空间与工程三维地质建模及应用研究
图5.10 变异函数中的球状模型、指数模型、高斯模型
D.变异函数对空间变异性结构的反映
变异函数作为定量描述空间变异性的一种统计学工具，通过其自身的结构及其各项参数，从不同角度反映了空间变异性结构。利用变异函数可以对空间变量的连续性、相关性、变量的影响范围、尺度效应、原点处的间断性、各向异性等要素进行描述。
E.变异函数理论模型参数拟合
变异函数理论模型参数拟合就是利用原始采样点数据或实验变异函数取值对所选定的理论模型参数以特定的方法进行估计。拟合方法一般采用手工拟合法。
手工拟合就是依据实验变异函数的取值，一方面通过观察实验变异函数图；另一方面对所研究的区域化变量进行必要的分析，采用肉眼观察来确定变异函数模型参数，并对参数反复进行交叉验证，最终确定模型参数。其拟合的大致过程如下：
(1)首先对所研究的区域化变量进行必要的结构、背景等方面的分析，结合专家经验，确定变异函数理论模型。
(2)利用实验变异函数散点图确定变异函数参数中的块金常数、基台值、变程、各向异性角度以及各向异性比值。
(3)交叉验证。
数学领域中的一些著名悖论及其产生背景
如果那个问题采用我的的话 这个分是不是也给我呢
