再探体渲染(二)--Ray Marching 从理论到实现

前文中探讨了基于物理的体渲染的理论以及相关结论的推导。在本文中，我将展示如何按照前文推导的公式，在unity中实现一套基于ray marching的”丐版“体积雾(只支持方形，全局只能有一个volume)。最后，我们会简单讨论该方案的一些优化方案。

理论回顾

关于体渲染的理论以及公式在上一篇文章中有详细介绍。

这里只列出几个重要的公式。

体渲染公式如下所示:

$$L_i(p,\omega) = L_o(p'',-\omega) Tr(p''\to p) + \int_{p\to p'} Tr(p'\to p)\sigma_t(p',\omega)L_s(p',\omega) dt (1)$$

其中

$$Tr(p\to p') =\frac{L_o(p',\omega)}{L_o(p,\omega)}= e^{\int_0^d - \sigma_t(p+\omega t,\omega) dt} (2)$$ $$L_s(p,\omega) = \frac{\sigma_\alpha(p,\omega)}{\sigma_t (p,\omega)} L_e(p,\omega) + \frac{\sigma_{s}(p,\omega)}{\sigma_t(p,\omega)} \int \rho (p,\omega_i,\omega) L_i(p,\omega_i) d \omega_i (3)$$

且满足:

$$Tr(p\to p') = Tr(p\to p'') Tr(p''\to p') (4)$$

从理论到实现

Ray Marching本质上就是用黎曼积分求解该方程。说人话就是： 将积分区域分成等长的小区间，计算小区间内的矩形面积得到。而Ray marching的每一步对应一个小矩形面积。

先看方程的左半部分

由于Ray Marching在后处理步骤进行，可以通过直接采样 光线对应的屏幕空间位置得到 的值。

而对于 ，假设光线步进的步长为 ,最多步进 步，且 ，则利用性质(3)可以得到:

$$Tr(p''\to p) = \prod_{i=1}^{n}Tr(q_i\to q_{i-1})= \prod_{i=1}^{n}e^{-\int_{(i-1)s}^{is}\sigma_t(p,\omega)dt}=\prod_{i=1}^{n} e^{-s\sigma_t(q_{i-1},\omega)} (5)$$

对应的代码为:

float3 step = rayDir * rayStepDistance;
float4 transmittance  = 1.0f;
for (int i = 0;i < RAY_MARCH_STEP_COUNT;i++)
{
     rayPos += step;
     rayDistance -= rayStepDistance;

     if (intersect(rayPos, rayDistance))
     {
           float texDensity = tex3D(_VolumeTex, rayPos * _VolumeScale + _VolumeOffset).r;
           transmittance  *=  exp(-_Density * texDensity * _AttenuationColor * rayStepDistance);
     }
}

其中，transmittance 为 项。 分为三个参数控制：控制总体浓度的_Density ,控制浓度关于位置变化的_VolumeTex ,控制烟雾颜色的 _AttenuationColor。 而光线步进长度 为rayStepDistance。

下面考虑右半部分方程

$$\int_{p\to p'} Tr(p'\to p)\sigma_t(p',\omega)L_s(p',\omega) dt=\sum_{i=1}^{n} Tr(q_i\to p)\sigma_t(q_i,\omega) L_s(q_i,\omega)s$$

其中可以观察到 等于上述代码中第 次循环中计算的 transmittance 。而 可以直接计算得到，因此我们只需解决如何计算L_s(q_i,\omega) 即可。

观察(3)，由于目前我们不考虑粒子的自发光，自发光项可以去掉。则

$$\sigma_t(q_i,\omega)L_s(q_i,\omega)=\sigma_{s}(p,\omega)\int \rho (q_i,\omega_i,\omega) L_i(q_i,\omega_i) d \omega_i (6)$$

关于 ，其示意图如下图所示。

其中红色线代表光线从光源到中间点 的路径，黑色线代表光线从 到 的路径。黑色方框代表体积雾，虚线部分代表光线经过路径上穿过体积雾的部分。 为从光源到 的光线与体积雾区域的交点。

对于给定方向 根据公式(2)

$$L_i(q_i,\omega_i) = Tr(q'\to q_i) L_i(q',\omega_i)$$

带入(6)中得

$$= \sigma_s(p,\omega) \int \rho(p,\omega_i,\omega) Tr(q'\to q_i) L_i(q',\omega_i)d\omega_i$$

考虑方向为 ,强度为 的方向光源，则 其中 为在 方向上非0的delta分布。带入上式得

$$=\sigma_s(p,\omega) \int \rho(p,\omega_i,\omega) Tr(q'\to q_i) \delta(\omega_l;\omega_i) L_l d\omega_i$$ $$=\sigma_s(p,\omega) \rho(p,\omega_l,\omega) Tr(q'\to q_i) L_l (7)$$

(7)中唯一无法直接计算的只有 项。其计算方法分为以下几步

1. 发射光线 计算其与体积雾区域的交点
2. 利用公式(5)以及ray marching计算transmittance。

最终(7)的计算代码如下

float4 sourceFunction(float3 rayPos, float3 w)
{
    // 获取主光源的方向
    float3 dirToLight = GetMainLight().direction;
    // 获取主光源的强度
    float4 intensity = float4(GetMainLight().color, 1.0f);

    float3 invRayDir = 1.0f / (dirToLight + 1e-6f);

    float _, dstInsideBox;
    // 计算交点位置，其中dstInsideBox为光线在体积雾内部传播距离。
    rayBoxDst(_BoundMin.xyz, _BoundMax.xyz, rayPos, invRayDir, _, dstInsideBox);

    // 根据距离计算步长
    float rayStepDistance = dstInsideBox / SOURCE_RAY_MARCH_STEP_COUNT;
    float3 step = rayStepDistance * dirToLight;

    float4 transmittance = 1.0f;
    float4 sigma = 1.0f;
    for (int i = 0; i < SOURCE_RAY_MARCH_STEP_COUNT; i++)
    {
        rayPos += step;
        

        float texDensity = tex3D(_VolumeTex, rayPos * _VolumeScale + _VolumeOffset).r;
        transmittance *= exp(-_Density * texDensity * rayStepDistance * _AttenuationColor);

        // 计算 sigma_t (p,w)的值
        if (i == 0)
        {
            sigma = _Density * texDensity * _AttenuationColor;
        }
    }

    return sigma * phase(rayPos, dirToLight, w) * transmittance * intensity;
}

带入方程右半部分结果，主循环中ray marching的代码可以改为

float3 step = rayDir * rayStepDistance;
float4 transmittance  = 1.0f;
float4 source = 0.0f;
for (int i = 0;i < RAY_MARCH_STEP_COUNT;i++)
{
     rayPos += step;
     rayDistance -= rayStepDistance;

     if (intersect(rayPos, rayDistance))
     {
           float texDensity = tex3D(_VolumeTex, rayPos * _VolumeScale + _VolumeOffset).r;
           transmittance  *=  exp(-_Density * texDensity * _AttenuationColor * rayStepDistance);
           source += sourceFunction(rayPos, rayDir) * rayStepDistance * transmittance;
     }
}

最后完整的shader代码如下所示，unity 版本使用2021.3.39f1c1, 项目模板使用URP。该shader用于后处理pass中。

Shader "Hidden/VolumeRenderer"
{
    Properties
    {
        // 后处理传入的全屏幕RT
        _MainTex("Texture", 2D) = "white" {}
        // 体积雾 Volume 的 Bounding Box
        _BoundMin("Bound Min", Vector) = (-5, -5, -5, 0)
        _BoundMax("Bound Max", Vector) = (5, 10, 5, 0)
        // 体积雾的总体密度
        _Density("Density", float) = 0.01
        // 体积雾的颜色
        _AttenuationColor("Attenuation Color", float) = (1, 1, 1, 1)

        // 体积雾的噪声图
        _VolumeTex("Volume Texture", 3D) = "white"{}
        _VolumeOffset("Volume Texture Offset", Vector) = (1, 1, 1, 0)
        _VolumeScale("Volume Texture Scale", Vector) = (1, 1, 1, 0)
    }
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        LOD 100

        Pass
        {
            ZTest Always Cull Off ZWrite Off

            // HLSL 代码块。Unity SRP 使用 HLSL 语言。
             HLSLPROGRAM
            // 此行定义顶点着色器的名称。
            #pragma vertex vert
            // 此行定义片元着色器的名称。
            #pragma fragment frag

            // Core.hlsl 文件包含常用的 HLSL 宏和
            // 函数的定义，还包含对其他 HLSL 文件（例如
            // Common.hlsl、SpaceTransforms.hlsl 等）的 #include 引用。
            #include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Core.hlsl"

            // DeclareDepthTexture.hlsl 文件包含用于对摄像机深度纹理进行采样的
            // 实用程序。
            #include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/DeclareDepthTexture.hlsl"
            #include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Lighting.hlsl"
            // #include "TerrainFog.hlsl"

            #pragma enable_d3d11_debug_symbols

            // 结构定义将定义它包含哪些变量。
            // 此示例使用 Attributes 结构作为顶点着色器中的
            // 输入结构。
            struct Attributes
            {
                // positionOS 变量包含对象空间中的顶点
                // 位置。
                float2 uv : TEXCOORD0;
                float4 positionOS   : POSITION;
            };

            struct Varyings
            {
                // 此结构中的位置必须具有 SV_POSITION 语义。
                float4 positionHCS  : SV_POSITION;
                float2 uv : TEXCOORD0;
            };

            sampler _MainTex;
            sampler _VolumeTex;

            float4 _BoundMin;
            float4 _BoundMax;
            float _Density;

            float4 _VolumeOffset;
            float4 _VolumeScale;
            float4 _AttenuationColor;

            // 顶点着色器定义具有在 Varyings 结构中定义的
            // 属性。vert 函数的类型必须与它返回的类型（结构）
            // 匹配。
            Varyings vert(Attributes IN)
            {
                // 使用 Varyings 结构声明输出对象 (OUT)。
                Varyings OUT;
                // TransformObjectToHClip 函数将顶点位置
                // 从对象空间变换到齐次裁剪空间。
                OUT.positionHCS = TransformObjectToHClip(IN.positionOS.xyz);
                OUT.uv = IN.uv;
                // 返回输出。
                return OUT;
            }

            int intersect(float3 rayPos, float rayDistance)
            {
                return rayPos.x < _BoundMax.x&& rayPos.x > _BoundMin.x && rayPos.y < _BoundMax.y&& rayPos.y > _BoundMin.y
                    && rayPos.z < _BoundMax.z&& rayPos.z > _BoundMin.z && rayDistance > 0.0f;
            }

            float phase(float3 rayPos, float3 wi, float3 wo)
            {
                return 1.0 / (4 * 3.1415926f);
            }

            // 光线与方形bounding box求交函数
            void rayBoxDst(float3 boundsMin, float3 boundsMax,
                //世界相机位置      光线方向倒数
                float3 rayOrigin, float3 invRaydir, out float dstToBox, out float dstInsideBox)
            {
                float3 t0 = (boundsMin - rayOrigin) * invRaydir;
                float3 t1 = (boundsMax - rayOrigin) * invRaydir;
                float3 tmin = min(t0, t1);
                float3 tmax = max(t0, t1);

                float dstA = max(max(tmin.x, tmin.y), tmin.z); //进入点
                float dstB = min(tmax.x, min(tmax.y, tmax.z)); //出去点

                dstToBox = max(0, dstA);
                dstInsideBox = max(0, dstB - dstToBox);
            }

#define SOURCE_RAY_MARCH_STEP_COUNT (8)
            // 计算 Ls项
            float4 sourceFunction(float3 rayPos, float3 w)
            {
                // 获取主光源的方向
                float3 dirToLight = GetMainLight().direction;
                // 获取主光源的强度
                float4 intensity = float4(GetMainLight().color, 1.0f);

                float3 invRayDir = 1.0f / (dirToLight + 1e-6f);

                float _, dstInsideBox;
                // 计算交点位置，其中dstInsideBox为光线在体积雾内部传播距离。
                rayBoxDst(_BoundMin.xyz, _BoundMax.xyz, rayPos, invRayDir, _, dstInsideBox);

                // 根据距离计算步长
                float rayStepDistance = dstInsideBox / SOURCE_RAY_MARCH_STEP_COUNT;
                float3 step = rayStepDistance * dirToLight;

                float4 transmittance = 1.0f;
                float4 sigma = 1.0f;

                float texDensity = tex3D(_VolumeTex, rayPos * _VolumeScale + _VolumeOffset).r;
                sigma = _Density * texDensity * _AttenuationColor;

                // ray marching 计算 transmittance
                for (int i = 0; i < SOURCE_RAY_MARCH_STEP_COUNT; i++)
                {
                    rayPos += step;
                    
                    float texDensity = tex3D(_VolumeTex, rayPos * _VolumeScale + _VolumeOffset).r;
                    transmittance *= exp(-_Density * texDensity * rayStepDistance * _AttenuationColor);
                }

                return sigma * phase(rayPos, dirToLight, w) * transmittance * intensity;
            }

#define RAY_MARCH_STEP_COUNT (32)

            // 计算体渲染方程，其中返回值对应Tr项，inScatter对应右半部分Ls的积分
            float4 VolumeRayMarching(float3 worldPos, out float4 inScatter)
            {
                float3 rayPos = _WorldSpaceCameraPos;
                float3 rayDir = worldPos - rayPos;
                float  rayDistance = length(rayDir);
                rayDir /= (rayDistance + 1e-6f);

                // float rayStepDistance = 0.25f;

                float3 invRayDir = 1.0f / (rayDir + 1e-6f);
                
                float dstToBox, dstInsideBox;
                rayBoxDst(_BoundMin.xyz, _BoundMax.xyz,rayPos, invRayDir, dstToBox, dstInsideBox);

                float maxStepDistance = max(_BoundMax.x - _BoundMin.x, max(_BoundMax.y - _BoundMin.y, _BoundMax.z - _BoundMin.z));
                float rayStepDistance = maxStepDistance / RAY_MARCH_STEP_COUNT;


                rayPos += dstToBox * rayDir;
                rayDistance -= dstToBox;

                float3 step = rayDir * rayStepDistance;
                float4 transmittance = 1.0f;
                inScatter = 0.0f;

                for (int i = 0;i < RAY_MARCH_STEP_COUNT;i++)
                {
                    rayPos += step;
                    rayDistance -= rayStepDistance;

                    if (intersect(rayPos, rayDistance))
                    {
                        float texDensity = tex3D(_VolumeTex, rayPos * _VolumeScale + _VolumeOffset).r;
                        transmittance *=  exp(-_Density * texDensity * rayStepDistance * _AttenuationColor);
                        inScatter += sourceFunction(rayPos, rayDir) * rayStepDistance * transmittance;
                    }
                }
                
                return transmittance;
            }

            // 片元着色器定义。
            half4 frag(Varyings i) : SV_Target
            {
                
                #if UNITY_REVERSED_Z
                    real depth = SampleSceneDepth(i.uv);
                #else
                    // 调整 z 以匹配 OpenGL 的 NDC
                    real depth = lerp(UNITY_NEAR_CLIP_VALUE, 1, SampleSceneDepth(i.uv));
                #endif

                // float3 worldPos = ComputeWorldSpacePosition(i.uv, depth, UNITY_MATRIX_I_VP);
                float4 ndc = float4(i.uv * 2 - 1, depth, 1);

                #if UNITY_UV_STARTS_AT_TOP
                    ndc.y *= -1;
                #endif

                float4 worldPos = mul(UNITY_MATRIX_I_VP, ndc);
                worldPos /= worldPos.w;

                
                float4 volumeInScattering = 0.0f;
                float4 attenduation = VolumeRayMarching(worldPos, volumeInScattering);
                // 将体渲染方程的两部分相加，其中_MainTex采样对应Lo项。
                half4 customColor = tex2D(_MainTex, i.uv) * attenduation + volumeInScattering;//lerp( tex2D(_MainTex, i.uv) , half4(1, 1, 1, 1), density);
                return customColor;
            }
            ENDHLSL
        }
    }
}

渲染效果如图所示

效果还是很不错的，就是很费。

优化，优化，优化

下面提几个我在网上看到以及自己想到的一些可以优化的地方。

泰勒展开

由于 时， 因此

$$Tr(q_i \to q_{i - 1}) = e^{-s\sigma_t(q_{i-1},\omega)}=1-s\sigma_t(q_{i-1},\omega)$$

因此可以将计算transmittance的代码改为:

...
//transmittance *=  exp(-_Density * texDensity * rayStepDistance * _AttenuationColor);
transmittance *=  1 - _Density * texDensity * rayStepDistance * _AttenuationColor;
...

理论上来说应该能减少n次exp 指令的计算(n=ray marching步数)。优化后和之前比确实无太大差别。

基于Volume

不再在后处理阶段渲染，而是直接渲染透明的Volume。可以避免全屏的ray marching计算，且可以支持多个volume的渲染。

TODO

大步长+随机抖动

一个简单粗暴的优化方案便是增大ray marching的步长。然而，加大步长往往会让渲染结果出现明显分层现象。

...
//#define RAY_MARCH_STEP_COUNT (32)
#define RAY_MARCH_STEP_COUNT (16)
...
//#define SOURCE_RAY_MARCH_STEP_COUNT (8)
#define SOURCE_RAY_MARCH_STEP_COUNT (4)
...

而给每步增加一个随机偏移可以让这些分层的硬边缘模糊掉。

...// 随机函数
float rand(float x)
{
    return frac(sin(x) * 43758.5453);
}


...


 for (int i = 0;i < RAY_MARCH_STEP_COUNT;i++)
 {
    // 每步步进的距离不再是定值rayStepDistance, 而是在rayStepDistance的0.6倍到1.4倍之间浮动
    float ditterDistance = rayStepDistance * (1.4f -  rand(rayDistance) * 0.8f);
    // 随机偏移后的光线步进位置以及距离
    float  ditteredRayDistance = rayDistance - ditterDistance;
    float3 ditteredRayPos = rayPos + ditterDistance * rayDir;

    rayPos += step;
    rayDistance -= rayStepDistance;

    // 用随机偏移的位置和距离计算
    if (intersect(ditteredRayPos, ditteredRayDistance))
    {
        float texDensity = tex3D(_VolumeTex, ditteredRayPos * _VolumeScale + _VolumeOffset).r;
        transmittance *= 1 - _Density * texDensity * ditterDistance * _AttenuationColor;//exp(-_Density * texDensity * rayStepDistance * _AttenuationColor);
        inScatter += sourceFunction(ditteredRayPos, rayDir) * ditterDistance * transmittance;
    }
}

...

效果如下图所示，在远处没有明显的分层现象。取而代之的是近处的噪声。

降低渲染分辨率

降低渲染体积雾的RT的分辨率，可以显著减少ps调用的次数从而优化性能。

TODO