CN102819860A - 基于辐射度算法的次表面散射物体的实时全局光照方法 - Google Patents

Abstract

本发明提出了一种基于辐射度算法的次表面散射物体的实时全局光照方法，主要分为三个步骤，第一步：预计算，利用已有的近似dipole模型，和辐射度算法预计算出场景的能量传输矩阵，以稀疏矩阵的形式保存。第二步，利用GPU的快速稀疏矩阵库计算能量传输过程，获取每个面片所获得的能力。第三步，利用快速光子收集算法，获得最终的全局光照效果。本发明在一定的预计算基础上，可以获得实时的次表面散射物体的全局光照效果。

Description

基于辐射度算法的次表面散射物体的实时全局光照方法

技术领域

本发明涉及一种基于辐射度算法的次表面散射物体的实时全局光照方法。

背景技术

在真实感渲染中，辐射度算法在物理逼真渲染中占有很重要的地位，尤其是在全局光照中，他将光的传播看成了一种能量的传播，将光解释成为离散的粒子态能量从而具有了很好的物理基础，通过各种改进手段可以逼真实现除了光的衍射和干涉之外的任意效果。

次表面散射(subsurface scattering)是指光从物体表面某点进入物体，经内部散射，最终从物体表面其他顶点出射的光线传递过程。具有次表面散射性质的物体一般具有较为圆润的外观色泽，这种现象在半透明材质上表现的更为明显。由于次表面物体存在多次散射让光路变得尤为复杂，传统的path-tracing，以及monte-carlo的光线跟踪可以很好的模拟次表面散射物体的全局光照效果，但是代价太高，得到一个较好的效果图往往需要数天甚至数月的时间。

近年来，专家学者们在次表面散射物体的全局光照方面取得了一定的研究成果：

1）次表面散射

次表面散射（BSSRDF）其在数学表达上是一个八个维度的函数：

其入射和出射方向角可以用参数

和

来表示，其入射位置和出射位置可以用参数p_i(x_i,y_i)和p_o(x_o,y_o)来表示，则其数学表达形式为：

$S(p_i,{\mathrm{\ω}}_i;p_o,{\mathrm{\ω}}_o)=\frac{\mathrm{dL}(p_o,{\mathrm{\ω}}_o)}{\mathrm{dE}(p_i,{\mathrm{\ω}}_i)}$

其中：表示BSSRDF的函数值，表示出射光线的强度强，表示入射光线的辐照度。

最初利用path-tracing以及ray-tracing来模拟次表面散射效果，虽然可以得到较好的效果，但是代价太大。2001年，Jensen等在假定次表面散射表面为无限大半平面的情况下，简化了漫散射模型，提出了偶极光源(Dipole)近似的方法将BSSRDF中光照模型分解为单次次表面散射和多次次表面散射，并在多次次表面散射中首先使用了偶极光源的方法进行了近似计算，可以获得较快的绘制速度。2003年，对于部分次表面散射物体都无法满足无限大半平面的假设，Craig Donner等提出了dipole模型的修正模型multi-dipole模型，采用多个偶极子来逼近新的临界条件。Craig Donner等在2009年提出了一种经验公式的BSSRDF的多次散射计算模型，通过测量得到的大量的数据来模拟次表面散射的真实状况，实验证明可以获得很好的仿真效果。

2）预计算辐射度算法

代预计算辐射传输方法是近年来用于解决复杂光源下物体真实感绘制的实时计算框架，由Sloan等最先提出。算法的基本思想是将绘制方程中的函数分解到指定的一组基函数上，用这组基函数的线性组合来逼近原函数，在实时绘制时，利用指定基函数之间的积分的良好性质来加速绘制方程的积分运算。

Hao等人提出了种基于球面谐波的实时绘制方法。该方法是对Jensen方法的进一步改进，算法依然采用两遍分级的方法进行绘制。在预计算环节算法利用球面谐波的球面径向基函数间积分易于计算的性质对预计算积分进行压缩在存储。为了弥补球面谐波压缩对高频分量的损失，文章设计了一种相关点算法来保存损失的高频分量。为了提高算法的绘制效率，文章对BSSRDF模型进行了适当的化简，假设半透明介质是均匀介质(BSSDF)，基础光照模型和绘制方程采用经典的Jensen算法进行模拟实现。文章提出了一种量化光源的概念来降低算法的时间复杂度，预计算部分主要计算量化光源的向量积分并对其进行存储。对于一个新的量化光源，采用采样邻居光源插值方式来计算其向量值。

对预计算光源向量积分使用球面谐波进行压缩存储，作为第二遍的直接输入部分。算法采用Cohn等提出的三维网格化简方法来得到相关点网格，得到最终全频预计算压缩数据。

Wang于2005也提出了一种基于预计算辐射传输全频光线算法。文章通过对单次次表面散射进行预计算和对相位角函数进行相似估计，很大程度的提高了最终的绘制效率。算法首先对一个的立方图光源进行参数化。算法把照度和光传输渲染到OPENGL的pbuffer中，对于多次次表面散射，需要把实体模型的顶点均衡的光栅化到立方图的六个面上，并对应的设置模板缓冲区，最后通过绘制方程得到多次次表面散射的出射光照强度。对于单次次表面散射，使用光线跟踪算法对每一个采样顶点预计算出其积分上限值，使用的超采样把化简模型渲染到pbuffer中，对于每一个片段输出顶点法线和深度纹理。最后对每一个采样顶点按照绘制方程进行叠加。对于不出的相位函数采用不同的预计算方式，为了达到更高的绘制效率，上述过程均在GPU中完成。

在国内，Wang也提出了一种基于预计算的半透明物体实时绘制方法。算法通过主元分析将偶极光源近似模型中的漫散射(多次次表面散射)材质分解为与形状相关函数和半透明材质相关函数的乘积形式。利用此分解表示，在预计算辐射传输的实时真实感绘制框架下，通过对散射传输的预计算来实现在多种光源环境下对半透明物体材质的实时编辑绘制。同时文章还提出了一种对预计算辐射传输的数据在空域上进行二次小波压缩的方法，利用表面点在空间分布位置的相关性，在保证绘制质量的前提下，大大压缩了数据，提升了绘制效率。

发明内容

本发明要解决的技术问题是：目前所提的方法中，path-tracing以及ray-tracing虽然可以解决次表面物体的全局光照问题，但是速度难以接受，预计算的辐射度算法，只能处理无限远的光以及环境光下的次表面物体的全局光照问题。而本发明可以解决任意光源光照下的次表面散射物体的全局光照效果，并速度上可以达到实时。

本发明解决上述技术问题的技术方案为：一种基于辐射度算法的次表面散射物体的实时全局光照方法，包括如下步骤：

步骤（1）、利用面片均匀的场景和辐射度算法以及dipole近似多次散射模型，建立一个全局的能量传输矩阵；所述的全局的能量传输矩阵为每对面片间能量传输的衰减变化情况，该矩阵的具体计算方法为：次表面物体到次表面物体面片间的内部能量传输矩阵计算利用dipole近似多次散射模型来，面片之间的外部能量传输矩阵则使用辐射度算法来计算；

步骤（2）、从光源开始对场景做到场景的初始能量计算，再利用步骤（1）中获取的能量矩阵，通过矩阵运算获得最终的能量分布计算结果；其中光源到场景的初始能量计算是利用GPU的光栅化的过程，首先将面片投影到光源坐标系下的屏幕空间中，通过面片所占的像素点的个数来得到面片从光源可以获取的初始光照值；统计每个面片的像素点总数，我们利用CUDA来并行完成；而计算最终能量分布的矩阵运算我们则利用cuSparse库来完成；

步骤（3）、将步骤（2）中获得的每个面片包含的能量，虚拟成在面片中间的能量光子，利用已有的屏幕空间的光子收集算法进行光子收集过程；从而计算最终能量分布。

进一步的，所述的次表面物体到次表面物体面片的计算则利用dipole近似多次散射模型具体为：计算次表面散射面片间的能量传输矩阵计算记为S，采用diffusion approximation的近似模型dipole模型来计算，dipole模型将一个多维的BSSRDF函数简化为了一个只与距离有关的dipole模型，计算公式为：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{\mathrm{tr}}=\sqrt{3{\mathrm{\σ}}_a{\mathrm{\σ}}_t^{\′}}\\ d_r=\left|\right|x-x_r\left|\right|\\ d_V=\left|\right|x-x_V\left|\right|\end{array}\right.$

其中，||x′-x||为入射点到出射点间的距离，x′为入射点的位置，x为出射点的位置，σ′_a为散射系数，σ_a为吸收系数，σ′_t为湮灭系数是σ′_a与σ_a的和，σ_tr为有效湮灭系数，d_r为出射点到虚拟正光源的距离，d_v为出射点到虚拟负光源的距离，x_r为虚拟正光源的位置，x_v为虚拟正光源的位置，z_v为入射点到虚拟负光源的位置；利用dipole模型为只跟两点距离有关的特性，先将dipole曲线离散成有限个点，若需要计算的值在两点之间，则进行线性插值，通过这样的方式，预计算当前场景所有次表面散射面片两两间的距离是处于曲线哪两点间，在后面计算能量传输时通过索引索引矩阵就可以得到真正的次表面散射传输矩阵，从而达到实时材质的变化。

进一步的，所述的面片之间的外部能量传输矩阵的计算方法具体为：使用辐射度算法中的形状因子形状因子的推导公式为：

$F_{\mathrm{ij}}=\frac{1}{A_i}\∫\∫V(x,y)\frac{\cos \mathrm{\θ}\cos {\mathrm{\θ}}^{\′}}{\mathrm{\π}{r}^2}d{A}_{y}d{A}_x$

其中V（x,y）为Ax，Ay平面间任意两点间的可见性函数，假定平面的大小远小于平面的距离，则简化和离散推导公式后得到计算所用的公式为：

$F_{\mathrm{ij}}=A_j\frac{\cos \mathrm{\θ}\cos {\mathrm{\θ}}^{\′}}{\mathrm{\π}{r}^2}V(i,j)$

其中V（i，j）为面片i中心到面片j中心的可见性函数。

进一步的，所述的计算最终能量分布具体为通过第一部分通过索引矩阵构建次表面散射矩阵S，第二部分，用矩阵运算获取最终能量分布。

进一步的，首先在第一个部分，通过读入的材质，计算每个离散点的dipole值，通过建立的索引矩阵，利用CUDA并行计算出S矩阵中每项的值；

第二部分，获得的每个面片的初始能量建立一个初始能量向量L，利用cusparse库，以及下面的公式，用矩阵运算模拟真正的能量传输过程，

L_flnal＝(S×F)^N-1×S×L

其中，L_flnal则为光线传输N次后得到的最终能量值，N的大小由光源衰弱的情况来定，设定一个阈值，当低于这个阈值时，则停止迭代。

本发明的原理在于：

本发明首先预计算，利用已有的近似dipole模型，和辐射度算法预计算出场景的能量传输矩阵，以稀疏矩阵的形式保存。接着，从光源出发获得，利用GPU获得场景的初始能量分布，再利用GPU的快速稀疏矩阵库计算能量传输过程，获取每个面片所获得的最终能量。最后，利用快速光子收集算法，获得最终的全局光照效果。

本发明与现有技术相比的优点在于：

与原有算法相比，本算法的主要贡献体现在以下两点：一，利用预计算的结果，以及cusparse的GPU快速矩阵乘法库，可以达到实时的次表面散射物体的全局光照效果。二，由于实时渲染部分与材质的无关性，可以使得次表面散射物体的材质更换达到实时的效果。

附图说明

图1为次表面散射全局光照算法流程图；

图2为辐射度算法中的形状因子计算公式参数含义示意图；

图3为点光源示意图；

图4为次表面散射全局光照示意图；

图5为cornell-box场景效果图；

图6为Bowls and bunny场景效果图；

图7为Rose场景效果图。

具体实施方式

下面结合附图和具体实施方式进一步说明本发明。

本发明的基于辐射度算法的次表面散射物体的全局光照过程如图1，其实现方法具体如下：

1）预计算能量传输矩阵

预计算阶段主要分为两个步骤，步骤一是面片间的能量传输矩阵的计算记为F，该步骤使用辐射度算法中的形状因子计算公式可得。形状因子的推导公式为：

$F_{\mathrm{ij}}=\frac{1}{A_i}\∫\∫V(x,y)\frac{\cos \mathrm{\θ}\cos {\mathrm{\θ}}^{\′}}{\mathrm{\π}{r}^2}d{A}_{y}d{A}_x$

参数含义如图2所示。

其中V（x,y）为AxAy平面间任意两点间的可见性函数，我们假定平面的大小远小于平面的距离，则简化和离散推导公式后可得我们计算所用的公式为：

$F_{\mathrm{ij}}=A_j\frac{\cos \mathrm{\θ}\cos {\mathrm{\θ}}^{\′}}{\mathrm{\π}{r}^2}V(i,j)$

其中V（i，j）为面片i中心到面片j中心的可见性函数。

步骤二，计算次表面散射面片间的能量传输矩阵计算记为S，这里我们采用diffusionapproximation的近似模型dipole模型来计算，dipole模型将一个多维的BSSRDF函数简化为了一个只与距离有关的dipole模型，计算公式为：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{\mathrm{tr}}=\sqrt{3{\mathrm{\σ}}_a{\mathrm{\σ}}_t^{\′}}\\ d_r=\left|\right|x-x_r\left|\right|\\ d_V=\left|\right|x-x_V\left|\right|\end{array}\right.$

其中，||x′-x||为入射点到出射点间的距离，x′为入射点的位置，x为出射点的位置，σ′_a为散射系数，σ_a为吸收系数，σ′_t为湮灭系数是σ′_a与σ_a的和，σ_tr为有效湮灭系数，d_r为出射点到虚拟正光源的距离，d_v为出射点到虚拟负光源的距离，x_r为虚拟正光源的位置，x_v为虚拟正光源的位置，z_v为入射点到虚拟负光源的位置。

利用dipole模型为只跟两点距离有关的特性，我们可以先将dipole曲线离散成有限个点，若需要计算的值在两点之间，则进行线性插值，通过这样的方式，可以预计算计算当前场景所有次表面散射面片两两间的距离是处于曲线哪两点间（称为计算索引矩阵）替代计算真正的矩阵，在后面计算能量传输时只需要通过索引这个索引矩阵就可以得到真正的次表面散射传输矩阵，这样就可以达到实时材质的变化。

2）计算从光源开始到场景的初始能量分布

这里分两类光源讨论，点光源，无限远平行光，聚光灯可以分为类别一，面光源单独作为类别二。

类型一的光源我们以点光源为例，如图3所示：

由点光源的特性可知，在一定的立体角上能量辐射是一定的，利用这个特性，将场景通过透视投影，投影到屏幕坐标系下，然后计算每个像素点的能量，然后根据面片所占像素点，计算出他从光源处应得的能量。具体流程如下：

（1）将场景投影到光源坐标系下；

（2）根据近裁剪面距离r，点光源的辐射公式以及近裁剪面实际大小A，在CUDA下计算并累加能量值到各个面片上。

对于无限远平行光源则只需要将透视投影改为平行投影即可。

对于类型二，面光源，我们可以把它当做一个已有初始能量的面片，而场景中其他面片的能量为0，则可以直接跳入计算最终能量分布阶段。

3）计算最终能量分布

主要分为两个部分，第一部分通过索引矩阵构建次表面散射矩阵S，第二部分，用矩阵运算获取最终能量分布。

首先在第一个部分，我们通过读入的材质，计算每个离散点的dipole值，通过1）中建立的索引矩阵，利用CUDA（计算统一设备架构）并行计算出S矩阵中每项的值。

第二部分，我们利用步骤2）中获得的每个面片的初始能量建立一个初始能量向量L，利用cuSparse库，以及下面的公式，用矩阵运算模拟真正的能量传输过程。

L_flnal-(S×F)^N-1×S×L

L_flnal则为光线传输N次后得到的最终能量值，N的大小由光源衰弱的情况来定，我们设定一个阈值，当低于这个阈值时，则停止迭代。示意图如图4所示：

4）最终效果生成

利用步骤3）生成的每个面片的能量，在每个面片中心生成一个光子，为了方便收集过程，我们以每个光子为中心，以收集半径r为一个外接圆，每个光子绘制一个正十六面体，这样一来，可以在屏幕空间收集光子而不需要建立kd-tree之类的数据结构。

我们最后对三个场景做了速度测试，结果如下表所示：

Claims (5)

1.一种基于辐射度算法的次表面散射物体的实时全局光照方法，其特征在于：包括如下步骤：

步骤（1）、利用面片均匀的场景和辐射度算法以及dipole近似多次散射模型，建立一个全局的能量传输矩阵；所述的全局的能量传输矩阵为每对面片间能量传输的衰减变化情况，该矩阵的具体计算方法为：次表面物体到次表面物体面片间的内部能量传输矩阵的计算利用dipole近似多次散射模型，而面片之间的外部能量传输矩阵则使用辐射度算法来计算；

步骤（2）、从光源开始对场景做到场景的初始能量计算，再利用步骤（1）中获取的能量矩阵，通过矩阵运算获得最终的能量分布计算结果；其中光源到场景的初始能量计算是利用GPU的光栅化的过程，首先将面片投影到光源坐标系下的屏幕空间中，通过面片所占的像素点的个数来得到面片从光源可以获取的初始光照值；统计每个面片的像素点总数，我们利用CUDA来并行完成；而计算最终能量分布的矩阵运算我们则利用cuSparse库来完成；

步骤（3）、将步骤（2）中获得的每个面片包含的能量，虚拟成在面片中间的能量光子，利用已有的屏幕空间的光子收集算法进行光子收集过程；从而计算最终能量分布。

2.根据权利要求1所述的一种基于辐射度算法的次表面散射物体的实时全局光照方法，其特征在于：所述的次表面物体到次表面物体面片的计算则利用dipole近似多次散射模型具体为：计算次表面散射面片间的能量传输矩阵计算记为S，采用diffusion approximation的近似模型dipole模型来计算，dipole模型将一个多维的BSSRDF函数简化为了一个只与距离有关的dipole模型，计算公式为：

$\left\{\begin{array}{c}{\mathrm{\σ}}_{\mathrm{tr}}=\sqrt{3{\mathrm{\σ}}_a{\mathrm{\σ}}_t^{\′}}\\ d_r=\left|\right|x-x_r\left|\right|\\ d_Y=\left|\right|x-x_Y\left|\right|\end{array}\right.$

其中，||x′-x||为入射点到出射点间的距离，x′为入射点的位置，x为出射点的位置，σ′_a为散射系数，σ_a为吸收系数，σ′_t为湮灭系数是σ＇_a与σ_a的和，σ_tr为有效湮灭系数，d_r为出射点到虚拟正光源的距离，d_v为出射点到虚拟负光源的距离，x_r为虚拟正光源的位置，x_v为虚拟正光源的位置，z_v为入射点到虚拟负光源的位置；利用dipole模型为只跟两点距离有关的特性，先将dipole曲线离散成有限个点，若需要计算的值在两点之间，则进行线性插值，通过这样的方式，预计算当前场景所有次表面散射面片两两间的距离是处于曲线哪两点间，在后面计算能量传输时通过索引索引矩阵就可以得到真正的次表面散射传输矩阵，从而达到实时材质的变化。

3.根据权利要求1所述的一种基于辐射度算法的次表面散射物体的实时全局光照方法，其特征在于：所述的面片之间的外部能量传输矩阵的计算方法具体为：使用辐射度算法中的形状因子形状因子的推导公式为：

$F_{\mathrm{ij}}=\frac{1}{A_i}\∫\∫V(x,y)\frac{\cos \mathrm{\θ}\cos {\mathrm{\θ}}^{\′}}{\mathrm{\π}{r}^2}d{A}_{y}d{A}_x$

其中V（x,y）为Ax，Ay平面间任意两点间的可见性函数，假定平面的大小远小于平面的距离，则简化和离散推导公式后得到计算所用的公式为：

$F_{\mathrm{ij}}=A_j\frac{\cos \mathrm{\θ}\cos {\mathrm{\θ}}^{\′}}{\mathrm{\π}{r}^2}V(i,j)$

其中V（i，j）为面片i中心到面片j中心的可见性函数。

4.根据权利要求1所述的一种基于辐射度算法的次表面散射物体的实时全局光照方法，其特征在于：所述的计算最终能量分布具体为通过第一部分通过索引矩阵构建次表面散射矩阵S，第二部分，用矩阵运算获取最终能量分布。

5.根据权利要求4所述的一种基于辐射度算法的次表面散射物体的实时全局光照方法，其特征在于：首先在第一个部分，通过读入的材质，计算每个离散点的dipole值，通过建立的索引矩阵，利用CUDA并行计算出S矩阵中每项的值；

第二部分，获得的每个面片的初始能量建立一个初始能量向量L，利用cusparse库，以及下面的公式，用矩阵运算模拟真正的能量传输过程，

L_flnal＝(S×F)^N-1×S×L

其中，L_flnal则为光线传输N次后得到的最终能量值，N的大小由光源衰弱的情况来定，设定一个阈值，当低于这个阈值时，则停止迭代。

Images