使用预计算散射的方式实时基于物理的渲染云

作者：< GPU Pro 6 > -----Egor

​ 真实感强烈的云，在游戏中一直是一种需求度很高的东西。云微观上是由会散射光的水分子构成的。在云的渲染中，难点是一束光从一片云射入到射出，要经历多次的散射，而要去模拟，抽样计算散射是比较昂贵的计算。因此，我们需要建立一个相对简单的计算模型来进行模拟。

​ 使用面向摄像机的广告牌技术，可能是最普遍的做法了。但是广告牌是一个平面，这在很多视角下会丧失它的立体感。这种方法还有其他限制：光照效果是依靠美术去画的，这也会使它在很多角度看上去会失去真实感。除此之外，体积渲染也是一种渲染云的方法。为了避免锯齿，在同一个方向上通常需要采样多个采样点，这会造成性能瓶颈，尤其是在高分辨率的机器上（PS:体积渲染技术的问题）。当然还有更多，在物理上更精确的方法，使得表现效果更佳，但很难降低计算量，性能消耗较为巨大。

​ 我们发现了一种新的基于物理的方法，且能有效率的渲染云。假定一朵云是由单个被我们称为“参考粒子”的粒子复制后通过扩张或者旋转构成的。在渲染之前，我们可能到达的摄像机位置，对视线上的方向，预计算其深度，和散射积分，并把数据保存。当运行时，我们加载数据，避免消耗大的步进光线采样或者slicing(极线采样)。比较从前的做法，我们做了如下的改进：

• 一个更好的实时着色模型，基于预计算查表
• 一个使用3D网格计算云内光照衰减的改善方法
• 一个新的粒子生成算法
• 效率优化，包括GPU粒子的排序

我们简要回顾了方法的主要概念，在这里我们主要讨论实现细节和改进。其他信息可去查看原始论文

2 光线传输原理

现在我们将简单介绍下在传播介质中光传输的主要概念。在传播媒介中，可以发现有三种现象：散射，吸收，以及射出。散射仅仅只改变光线传播的方向，吸收则是吸收光的能量，而射出则与之相反。这三个现象的强度我们用三个参数来进行描述：

吸收和散射都会减少介质中的光的能量。他的衰减参数:

在云内部，吸收和射出是可以忽略不计的，即：

所以散射和衰减可以表述为同一个参数，用：

表示。

​ 当一束光在云内部，由A点传输到B点的衰减为，我们称为光学深度(Optical Depth)：

​ 为了得到单一光照的被散射后的强度，我们需要使用步进采样进行积分。每一个步进方向：

在这个等式中，C是当前摄像机的位置,\(\vec{v} $ 是视角方向。P0和P1是云对于视线的入口和出口两点，\)L_{Sun}\( 是云外部的阳光强度，Q是阳光到达当前积分点P的点。\)P(\theta)\( 是定义 散射出去多少能量的phase函数，\)\theta\( 是观察方向PC和出射方向QP的夹角。注意，阳光在到达相机前，会衰减两次：即\)e^{-\tau(Q,P)}\( 从阳光入口Q，到达散射点P，和\)e^{-\tau(P,P0)}$ 从散射点P，到视线的入口P0

​ 云的phase函数是非常复杂的。在实时渲染方法中，通常使用Cornette-Shanks方法来近似它：

​ 使用单一散射强度\(L_{In}^{(1)}\) ，我们可以计算二次散射\(L_{In}^{(n)}\) ，然后三次散射也是如此。然后n次散射强度积分如下

在公式2.4中\(J^n(C,\vec{v})\) 是n-1阶散射的净强度

这是对P点的整个球面积分，\(\vec{w}\)是入射方向，\(\theta $ 是\)\vec{w} 和 \vec{v}$ 的夹角

总的散射强度计算如下：

摄像机最终测量的辐射度是总散射强度和背景辐射强度之和：

3 预计算解决方案

​ 上述等式，计算特别复杂，是不可能实时进行计算的。我们的解决这个问题的办法是，对参考体积粒子中的光传输进行建模，在预处理的时候求解出该粒子的所有方程式。而后，我们将结果保存起来，在实时着色时读取。

3.1 光学深度(Opitical Depth)

​ 考虑一些知道密度分布的不均匀体积粒子。我们的目标是通过粒子为每一个摄像机位置和视图方向预计算公式（2.1）中的光学深度积分。为了描述每一条穿过粒子的射线，我们需要四维向量参数。前两个参数维度是方位角，即球坐标的两个角度\( \varphi_S \in [0,2\pi] ,\theta_S \in[0,\pi] $，表示入射点S，后面两个参数也是方位角，即\)\varphi_\nu \in[0,2\pi],\theta_\nu \in [0,\frac{\pi}{2}]\( 表示在入口点处构建的切线空间中，视线的角度方向。该切线空间的构建中，轴指向球心。因为我们仅仅需要考虑穿过球的视线，所以另外一个角的最大值是 表示在入口点S处构建的切线空间中，视线的角度方向。该切线空间的构建中，z轴指向球心。因为我们仅仅需要考虑穿过球的视线，所以另外一个角的最大值是\)\frac{\pi}{2} $。

​ 为了预计算光学深度积分，我们遍历所有的\(\varphi_S,\theta_S,\varphi_\nu,\theta_\nu\) ，并利用公式2.1计算其积分值，在3.5中提供细节描述。

3.2单次散射(Single Scattering)

​ 对比Opitical Depth，我们无法预计算不均匀粒子的散射，原因是，我们需要考虑光线方向，这就需要5个参数了，是不切实际的。所以我们预计算均匀分布球形颗粒中的散射。我们假定光的方向是正z轴，因为光场，相对光是对称的，所以我们可以舍弃 \(\varphi_S $ ,另一方面，为了计算公式2.5，我们需要知道整个体积的光场，而不仅仅是球表面的光场。因此，我们将从球心到起点的距离作为第四个参数。所以我们计算单次散射的参数为\)\theta_S \in [0,\pi] ,\varphi_\nu \in[0,2\pi],\theta_\nu \in[0,\pi],r \in[0,1]。请注意因为光场需要覆盖整个球面角，所以 。请注意因为光场需要覆盖整个球面角，所以\theta_\nu的最大值是 的最大值是\pi$

​ 然后，是遍历所有参数值，使用2.2公式来预计算单次散射。因为假定了粒子是均匀的，\(\beta(P) \equiv \beta\) 。阳光强度\(L_{Sun}，和phase函数P(\theta),被分解出来，单独计算\) 细节在3.5.1提及。

3.3多次散射(Multiple Scattering)

​ 我们使用与单次散射相似的参数对多次散射预处理。为了预计算多次散射，我们逐步执行以下的步骤：

• 根据2.5，对所有参数项使用先前的阶数\(L_{In}^{n-1}\) 去计算J(n)项
• 根据2.4，计算当前次散射\(L_{In}^{n}\)
• 累加到当前散射序列和之中

实现细节，在3.5.1中描述

​ 在散射序列计算结束后，我们只保留表面上的光场，丢弃其余的数据。必须注意的是，相对于Optical Depth，在散射中，密集度与粒子半径不是线性相关的。在原始论文中，我们计算了许多密度的散射，并且在4D查找表中，对结果进行了编码，第四个参数是粒子密度的scale。后来发现，仅仅使用一个密度，也工作得挺好，而且简化了算法。

​ 这一步将在3.5.3中进一步讨论

3.4 Volume-Aware Blending(粒子融合)

​ 我们的云，由一系列的粒子构成，我们需要把他们混合成一个完整的传播媒介。实现这一目标的方法，通常是用alpha 融合。这种方法主要用于薄的物体，比如玻璃或者树叶。我们的粒子是体积渲染实体，没办法使用通用的alpha混合方式进行计算。为了解决这个问题，我们发现了一种新的技术，我们称为“Volume-Aware Blending”。这个技术的关键点是跟踪离摄像机最近的体积元素，针对每一个像素去混和每一个新的粒子。

​ 算法第一步，清掉Closet element Buffer 和 back buffer。然后从后往前渲染体积粒子。针对当前最近的粒子，对每一个粒子进行测试。如果这个粒子比目前的最近粒子更接近相机，那么对最近粒子使用alpha混合，并将颜色写入back buffer，然后使新粒子取代最近粒子。如果新粒子比最近粒子远，那么新粒子就会混合到back buffer里头，最近粒子保持不变。

​ 如果最近粒子与新粒子相交，那么事情就会变得复杂一些。首先把尾部融合进buffer，然后使用密度加权平均的方法融合相交部分的颜色

\(C_0,C_1是非alpha预乘颜色，\rho_1,\rho_2是密度，d_i是相交的长度\) 。然后与头部的颜色进行混合后进buffer。

3.5 实现

例子是用C++的DXD11 的API实现的，全部的源代码可以在Github仓库中找到

3.5.1 预计算光传输路径

​ 我们将预计算数据做成一个3D一个4D的查找表。Optical Depth的积分保存在一个\(32 \times16\times32\times16(N_{\varphi s}=32,N_{\theta s}=16,N_{\varphi \nu}=32,N_{\theta \nu}=16)\) 的一个8位LUT(lookup table)中。多次散射保存在一个\(32 \times 64 \times16(N_{\theta S}=32,N_{\varphi \nu}=64,N_{\theta \nu}=16)\) 的16位 float的LUT中。第一个LUT需要0.25MB，第二个则需要64KB。注意与基础方法不同，我们不依赖预计算，而是使用别的方法来近似计算单次散射。

​ 因为当前图形设备不支持4D纹理，我们只能使用3D纹理来进行取代。比如一张\(X \times Y \times Z \times W的纹理，我们可以储存成X \times Y \times Z \cdot W\) 的3D纹理。我们对第四个坐标执行手动过滤。

#define SAMPLE_4D(tex3DLUT,LUT_DIM, f4LUTCoords, fLOD, Result)\
{\
    float3 f3UVW0;                                             \
    f3UVW0.xy=f4LUTCoords.xy;\
    float fQSlice=f4LUTCoords.w * LUT_DIM.w -0.5;\
    float fQ0Slice =floor(fQSlice);\
    float fQWeight=fQSlice-fQ0Slice;\
    f3UVW0.z=(fQ0Slice+f4LUTCoords.z)/LUT_DIM.w;\
    /* frac() assures wraparound filtering of w coordinate */
    float3 f3UVW1=frac(f3UVW0+float3(0,0,1/LUT_DIM.w));\
    Result = lerp(\
        tex3DLUT.SampleLevel(samLinearWrap,f3UVW0,fLOD),\
        tex3DLUT.SampleLevel(samLinearWrap,f3UVW1,fLOD),\
        fQWeight\
     );\
}

​ 注意\(\varphi_S 和 \varphi_\nu\) 坐标需要环绕滤波来避免伪影 。我们使用frac()函数，来让第四个坐标实现这个功能。也要注意z轴是不能使用环绕滤波模式的。

​ 预计算进程可以被总结如下：

• 预计算Optical Depth 积分
• 预计算单次散射将结果保存在32-bit的浮点型LUT
• 计算n次散射 从第2次到第N次：
• 计算\(J^n\) 部分
• 计算\(L_{In}^n\) 部分
• 累加\(L_{In}^n\) 并存储到LUT中
• 复制多次散射中，位于球面的结果，保存到最终的16-bit的LUT中

接下来对每一步的细节进行描述

Optical depth.

float2  PreComputeOpticalDepthPS(SQuadVSOutput IN):SV_Target
{
    float3 f3StartPos, f3RayDir;
    //Convert lookup table 4D coordinate into the start position and view direction
    OpticalDepthLUTCoordsToWorldParams(float4(ProjToUV(In.m_f2PosPS),g_Attribs.f4Params.xy),
        f3StartPos,f3RayDir);
    //Intersect the view ray with the unit sphere
    float2 f2RayIsecs;
    //x:first intersect length,y:second intersect length
    //f3StartPos is located exactly on the surface ; slightly move it inside the sphere to avoid       //precision issues
    GetRaySphereIntersection(f3StartPos + f3RayDir*1e-4,f3RayDir,0,1.f,f2RayIsecs);
    float3 f3EndPos=f3StartPos+f3RayDir*f2RayIsecs.y
    float fNumSteps=NUM_INTEGRATION_STEPS;
    float3 f3Step = (f3EndPos - f3StartPos) / fNumSteps;
    float fTotalDensity=0;
    for(float fStepNum=0.5; fStepNum<fNumSteps;++fStepNum)
    {
        float3 f3CurrPos = f3StartPos + f3Step * fStepNum;
        float fDensity=ComputeDensity(f3CurrPos);
        fTotalDensity+=fDensity;
    }
    return fTotalDensity/fNumSteps;
}

​ 需要使用步进算法，一步一步计算。首先计算输入4D坐标，使用OpticalDepthLUTCoordsToWorldParams()函数得出的开始位置和射线方向。4D参数的前两维来自 像素位置，另外两维被存储在g_Attribs.f4Params.xy 变量中。然后继续得到射线与单位圆的交点以及射线离开单位圆的出口。GetRaySphereIntersection()函数中，参数是射线开始坐标，和方向，球的中心点以及半径。然后就是根据公式2.1计算积分。我们不存储积分值，而是存储归一化的平均值。平均值正好能保存在8-bit的的table中。Optical Depth在使用时 将保存的平均数乘上距离即可。ComputeDensity()函数作用是，用pos采样3Dnoise去计算当前点的密度。

Single scattering

预计算单次散射的片元着色器代码如下。注意单次散射的计算要包含整个体积，而不仅仅是在表面上。使用一张4DLUT储存结果。坐标的第4维编码保存的是到球心的距离，由g_Attribs.f4Params.y提供。

float  PrecomputeSingleSctrPS(SQuadVSOutput In) : SV_Target
{
    float3 f3EntryPoint , f3ViewRay , f3LightDir;
    ScatteringLUTToWorldParams(
    float4(ProjToUV(In.m_f2PosPS), g_Attribs.f4Param.xy),
    g_Attribs.f4Param.z, f3EntryPoint, f3ViewRay, f3LightDir
    );
    float2 f2RayIsecs;
    GetRaySphereIntersection(f3EntryPoint, f3ViewRay, 0 , 1.f ,f2RayIsecs);
    float3 f3EndPos = f3EntryPoint + f3ViewRay * f2RayIsecs.y;
    float fNumSteps = NUM_INTEGRATION_STEPS;
    float3 f3Step = (f3EndPos - f3EntryPoint) / fNumSteps;
    float fStepLen = length(f3Step);
    float fCloudMassToCamera = 0;
    float fParticleRadius = g_Attribs.RefParticleRadius;
    float fInscattering = 0;
    for(float fStepNum=0.5,fStepNum < fNumSteps; ++fStepNum)
    {
        float3 f3CurrPos = f3EntryPoint + f3Step * fStepNum;
        GetRaySphereIntersection(f3CurrPos, f3LightDir, 0 , 1.f ,f2RayIsecs);
        float fCloudMassToLight = f2RayIsecs.x * fParticleRadius;
        float fAttenuation = exp( 
        -g_Attribs.fAttenuationCoeff *
        (fCloudMassToLight + fCloudMassToCamera)
        );
        fInscattering += fAttenuation * g_Attribs.fScatteringCoeff;
        fCloudMassToCamera += fStepLen * fParticleRadius;
    }
    return fInscattering * fStepLen * fParticleRadius;
}

​ 算法是依据公式2.2实现的。注意phase函数\(P(\theta ) 和 阳光强度 L_{Sun}\)被忽略了。因此着色器在每一步，需要去计算该式的积分:\(\beta(P) \cdot e^{-\tau(Q,P)} \cdot e^{-\tau(P,P_0)}\) ，散射衰减因数\(\beta(P)\)被看作是一个常数，由g_Attribs.fScatteringCoeff 变量提供。 我们使用的 \(\beta =0.07\) 作为散射衰减系数。一个参考粒子（reference particle）的半径为200米。衰减\(e^{-\tau(Q,P )}\) 是从当前点到光线与粒子球的交点入口的衰减。衰减\(e^{-\tau(P,P_0)}\) 朝向摄像机部分的衰减是由云的所有mass表示的，使用fCloudMassToCamera变量进行累加。

Multiple scattering

在单次散射之后，我们计算了N=18次的多次散射。在这一阶段，我们使用了3张-4D 32-bit 的浮点型LUT：一张存储\(J^n\)部分，一张存储当次散射\(L_{In}^n\) 。第三张，则存储所有次散射混合的结果。

关于\(J^n\) 的计算代码如下;

float GatherScatteringPS(SQuadVSOutput In) :SV_Target
{
    float3 f3StartPos, f3ViewRay , f3LightDir;
    ScatteringLUTToWorldParams
    (
        float4(ProjToUV(In.m_f2PosPS),g_Attribs.f4Param.xy),
        f3StartPos, f3ViewRay , f3LightDir
    );
    float3 f3LocalX, f3LocalY, f3LocalZ;
    ConstructLocalFrameXYZ(-normalize(f3StartPos), f3LightDir, 
    f3LocalX,f3LocalY,f3LocalZ);
    float fJ=0;
    float fTotalSolidAngle=0;
    const float fNumZenithAngles = SCTR_LUT_DIM.z;
    const float fNumAzimuthAngles = SCTR_LUT_DIM.y;
    const float fZenithSpan = PI;
    const float fAzimuthSpan = 2* PI;
    for(float Zen = 0.5;Zen < fNumZenithAnles; ++Zen)
    {
        for(float Az = 0.5; Az < fNumAzimuthAngles; ++Az)
        {
            float fZenith= Zen/fNumZenithAngles * fZenithSpan;
            float fAzimuth = (Az / fNumAzimuthAngles - 0.5) * fAzimuthSpan;
            float3 f3CurrDir = GetDirectionInLocalFrameXYZ(f3LocalX, f3LocalY, f3LocalZ, fZenith, fAzimuth);
            float4 f4CurrDirLUTCoords = WorldParmasToScatteringLUT(f3StartPos,f3CurrDir,f3LightDir);
            float fCurrDirSctr = 0;
            SAMPLE_4D(g_tex3DPrevSctrOrder , SCTR_LUT_DIM, f4CurrDirLUTCoords, 0 , fCurrDirSctr);
            if(g_Attribs.f4Param.w == 1)
                fCurrDirSctr *=HGPhaseFunc(dot(-f3CurrDir,f3LightDir),0.9);
            fCurrDirSctr *=HGPhaseFunc(dot(f3CurrDir,f3ViewDir),0.7);
            float fdZenithAngle = fZenithSpan / fNumZenithAngles;
            float fdAzimuthAngle = fAzimuthSpan / fNumAzimuthAngles *            sin(ZenithAngle);
            float fDiffSolidAngle = fdZenithAngle * fdAzimuthAngle;
            fTotalSolidAngle += fDiffSolidAngle;
            fJ += fCurrDirSctr * fDiffSolidAngle;
        }
    }
    fJ *= 4* PI / fTotalSolidAngle;
    return fJ;
}

这部分代码，第一步是从4Dtexture中取到世界空间下的参数(3-6行)，与之前的代码类似。之后的步骤是构建局部转换矩阵（8-10行）。这个函数需要两个向量当作输入来构建正交基。第一个向量使用z轴，注意z轴是指向球心的。

​ 然后构建两个循环，迭代zenith \( \theta\) 和 azimuth \(\varphi\) (18-19行) 。在每一步迭代中，着色器使用(\( \theta,\varphi $) 构建了采样方向。从转换这个方向为可采样的参数（26-28行）读取该方向上第n-1次散射(29-31行)。记得我们在计算单次散射的时候，没有计算phase 函数的。如果有必要，我们在这一步添加（32-34行）g_Attribs.f4Param.w的值是如果这是第二次散射则为1否则为0。之后我们用公式2.5计算phase 函数（35行），对于单次散射的各向异性系数我们使用0.9，而多次散射我们则使用0.7。最后计算\)dw=d\theta \cdot d \varphi \cdot \sin{\theta}$ 部分(37-40行)

​ 在\(J^n\) 部分计算完成后，我们接着计算第n次散射。这部分与之前的单次散射十分类似，只需要把\(J^n\) 读取出来，使之代替单次散射中的阳光衰减，我们还使用了梯度积分来提高准确性。

​ 第三阶段，便是，取出各次散射，并进行加法混合，并将累积结果存储起来。

3.5.2 粒子生成

​ 当我们渲染云的时候，我们想要能有效率的做LOD和对近处的云提供叫高的真实度。为了完成这个，我们使用了一个网格嵌套数据结构，这种结构的灵感来自几何图形贴图。每一个外环中的粒子是其内环中粒子的两倍大小，其粒子间距也是两倍。我们将这种结构称为一个单元网格。一个网格中的单元包含一个预定义数目的层数。每一个三维结构结果中的体素可能包含一个粒子。我们称这种结构是粒子格子。为了加速粒子生成和光照

我们维护了两个3D数据结构：云密度3D格子 和 光照衰减3D格子。这两个结构的分辨率每一维都是一个粒子格子的两倍，而且是使用3D贴图实现的。

粒子生成的过程如下：

• 处理一张2D 单元网格，去建立一个有效非空单元列表，并计算每一个单元的属性
• 计算位于非空单元中每一个云密度格子中的每一个voxel
• 处理位于非空单元中每一个光照衰减格子中的每一个可见的voxel
• 处理粒子格子，为可见的单元且密度超过了阈值的粒子格子生成粒子。
• 处理粒子并储存光照信息
• 对粒子排序

以上步骤的实施是通过一个compute shader去实现的。我们基于GPU去实现它，CPU不知道每一个Kernel，GPU需要多少的线程去运行。我们使用Dispatch Indirect() 函数来让GPU自己给自己分配。这个函数的参数和Dispatch()是一样的，但是这些参数是被储存在GPU Buffer中的，这就允许GPU控制它自己。我们接下来讨论，每一步的细节·

Processing cell grid

处理单元网格，是使用compute shader来完成的，为每一个单元分配一个线程。它负责基于摄像机的世界空间坐标计算每一个单元的中心点以及大小。使用计算出的单元中心位置，通过组合两个2Dnoise函数计算单元的基础密度。如果求出的结果超过阈值，那么就说明这个单元是合法的。着色器将所有合法的单元Append进Buffer中（g_ValidCellsAppendBuf）。最后能得到一个未排序的列表。如果单元在摄像机裁剪矩阵内，即是可见的，那么着色器将会把这些单元Append进另一个保存着所有可见的单元Buffer中(g_VisibleCellsAppendBuf)。

Processing cloud density lattice

下一个阶段，我们需要处理在合法单元网格内的那些lattice中的voxel。为了计算这些，需要一定数量的GPU线程，我们执行了一个简单的单线程compute shader：

RWBuffer<uint> g_DispatchArgsRW : register(u0);
[numthread(1,1,1)]
void ComputeDispatchArgsCS()
{
    uint s = g_GlobalCloudAttribs.uiDensityBufferScale;
    g_DispatchArgsRW[0] = (g_ValidCellsCounter.Load(0) * s *s *s *
                          g_GlobalCloudAttribs.uiMaxLayers + THREAD_GROUP_SIZE-1)
        /THREAD_GROUP_SIZE;
}

先前被写入buffer的元素可以用CopyStructureCount()函数拷贝进一个合适被读取的资源中。然后将先前作为UAV绑定到g_DispatchArgsRW的缓冲区传递给DispathIndirect()函数，以生成所需数量的线程。每一个线程读取g_ValidCellsUnorderedList Buffer中的合法单元，在上阶段被填满的cell内容，找出它在这个单元的位置。然后这个shader结合两个3Dnoise函数，创建一个单元密度的体积噪声。该噪声的幅度，随着高度降低，以产生典型的体积云形状，具有更宽的底部，和更窄的顶部。

Light attenuation

光照衰减是在每一个可见单元中的每一个voxel中被计算的。为了完成计算，我们需要分配一定数量的线程。使用与前阶段类似的方式，但是这次我们提供了具有合法性和可见性的单元buffer g_ValidCellsCounter变量。光照衰减是投一从voxel的中心朝光照方向的射线，通过密度lattice进行步进采样而计算出的。我们执行了16次步进。我们存储云衰减的mass系数，而不是直接存储光照衰减，是因为云衰减的mass能更好的被插值。

Particle generation

下一阶段，是处理cloud lattice中的有效可见的voxel，并为其中一些生成粒子。执行这一步，需要分配一定数目的线程。我们再次使用了简单的单线程compute shader。粒子生成器shader加载了云的密度，从密度lattice上，如果它不是0，那么就要创建粒子。shader随机的替换voxel中心的粒子，并加上一个随机的旋转和scale，以消除样式的重复性。着色器会将属性，比如粒子的坐标位置，密度，大小，写入粒子信息buffer中，并将粒子索引写入另一个append buffer中(g_VisibleParticlesAppendBuf)。

Processing visible particles

这一步，需要计算光照信息。特别是，我们在计算射到粒子中心的阳光颜色时，忽略了云的遮挡和周围天光的强度。我们也会采样光照衰减mass的贴图去计算光照遮挡。我们使用粒子表面的值去计算多次散射的衰减，用粒子中心的值去计算单次散射的衰减。此外，当计算多次散射的衰减时，我们让光照-衰减mass乘上一个0.25的系数去计算前向的强散射。

Sorting

在粒子能被正确渲染前，将粒子从后往前排序是最后的阶段了。在我们最开始的论文时，我们将粒子lattice的voxel排序是用CPU完成的，然后只有可见有效的voxel会从GPU中流出（用GPU做voxel的可见性和有效性判定）。这导致了一定数目的Drawbacks，首先他需要激活CPU-GPU通话。第二，由于随机的offset，粒子顺序和体素顺序可能有些微不同。但是主要的问题是，所有的voxel都会被排序，即使它们大部分是空的，这导致了CPU的overhead

我们现在将排序粒子完全用GPU去做，使用Satish et al 的合并排序算法（一个简化的合并算法程序）。我们开始时，将可见列表的粒子细分为128个粒子的子序列，并使用bitonic排序对每个子序列进行排序。然后我们执行合并操作，将已排好序的子序列合并成一个序列。当执行二进制搜索寻找其序号时，我们直接访问全局内存。因为需要排序的粒子相对较少通常不会超过50000，整个列表可以存入缓存，解释我们不适用共享内存，合并仍然非常有效。

重要的一点是，我们不知道GPU上生成了多少粒子以及我们需要执行多少次合并。因此我们执行了足够多数目的Pass，来对最大可能的粒子数目进行排序。当不需要再完成其他任务时，compute shader提前退出工作，性能成本时非常低的。

3.5.3 渲染

当前期的粒子生成，预处理和排序完成后，就可以进行渲染了。因为只有GPU知道我们生成了多少的粒子，所以我们使用DrawInstancedIndirect。它的使用与DrawInstance类似，只是它的参数，来自GPU Buffer。我们读取每一个应该可见的粒子，并从它的属性中读取数据，生成它的BoundingBox，最后将其送入光栅化流程。

在片元着色器中，我们对视线与粒子的bounding box求交，如果不相交，就Discard产生这条射线的这个像素。此外我们的着色模型是基于预计算的LUT的，正如下述代码。

//Compute lookup coordinates
float4 f4LUTCoords;
WorldParamsToOpticalDepthLUTCoords(f3EntryPointUSSpace, f3ViewRayUSSpace, f4LUTCoords);
//Randomly rotate the sphere
f4LUTCoords.y += ParticleAttrs.fRndAzimuthBias;
//Get the normalized density along the view ray
float fNormalizedDensity = 1.f;
SAMPLE_4D_LUT(g_tex3DParticleDensityLUT, OPTICAL_DEPTH_LUT_DIM, f4LUTCoords, 0, fNormalizedDensity);
//Compute actual cloud mass by multiplying the normalized
//density with ray length
fCloudMass = fNormalizedDensity * fRayLength;
fCloudMass *= ParticleAttrs.fDensity; 

//Compute transparency
fTransparency = exp( -fCloudMass * g_Attribs.fAttenuationCoff);

//Evaluate phase function for single scattering
float fCosTheta =dot(-f3ViewRayUSSpace, f3LightDirUSSpace);
float PhaseFunc = HGPhaseFunc(fCosTheta, 0.8);

float2 f2Attenuation = ParticleLighting.f2SunLightAttenuation;
//Compute intensity of single scattering
float3 f3SingleScattering = fTransparency * ParticleLighting.f4SunLight.rgb *
f2Attenuation.x * PhaseFunc;

//Compute lookup coordinates for multiple scattering
float4 f4MultSctrLUTCoords = WorldParamsToScatteringLUT(f3EntryPointUSSpace ,
    f3ViewRayUSSpace, f3LightDirUSSpace);
//Load multiple scattering from the lookup table
float fMultipleScattering = g_tex3DScatteringLUT.SampleLevel( samLinearWrap,f4MultSctrLUTCoords.xyz, 0);

float3 f3MultipleScattering = 
    (1-fTransparency) * fMultipleScattering * 
    f2Attenuation.y * ParticleLighting.f4SunLight.rgb;

//Compute ambient light
float3 f3EarthCentre = float3(0, -g_Attribs.fEarthRadius, 0);
float fEnttryPointAltitude length(f3EntryPointWS - f3EarthCentre);
float fCloudBottomBoundary =
    g_Attribs.fEarthRadius + g_Attribs.fCloudAltitude - g_Attribs.fCloudThickness / 2.f;
float fAmbientStrength =
    (fEnttryPointAltitude - fCloudBottomBoundary) / g_Attribs.fCloudThickness;
fAmbientStrength = clamp(fAmbientStrength,0.3,1);
float3 f3Ambient = (1-fTransparency) * fAmbientStrength * ParticleLighting.f4AmbientLight.rgb;

我们的第一步是计算密度,使用OpticalDepthLUT,我们随机的围绕垂直轴旋转一个角度，来减少重复度。f3EntryPointUSSpace和f3ViewRayUSSpace 是粒子空间下的起点和射线方向。（因为是单位圆空间，所以后缀是US）。接下来我们计算透明度

我们的实时渲染由三部分组成：单次散射，多次散射，和环境光。我们在第20-27行计算了单次散射。他是日照强度和日照衰减以及相位函数的乘积。因为单次散射在云密度低的地方最为明显，所以我们乘上一个透明度。

接下来我们计算多次散射。我们乘上了光照衰减强度。因为多次散射发生在云的密集区域，所以我们还要乘上不透明度（1-fTransprency）。

最后我们对环境光效果做一个近似。环境光的强度，在顶部边界是最强的，越往底部，越小。下图显示的是每一过程的结果

从左到右，依次是单次散射，多次散射，环境光，以及所有结果之和