【Unity/大气渲染】单次散射的原理和简单实现

11.6

当初写的效果比较一般，只有一个大概的样子，完成度比较低，回过头来有点看不下去

而且查找表，拟合曲线之类的优化也没有讲

重读论文后决定将文章重写。

——————————————————————————————————————

这篇随笔将会简单实现一个基于物理的相对真实的大气渲染效果

如下图，太空中的星球、相对真实的天空盒

 

 

 

如果没有大气，太阳光没有直接照射到的地方将会是一片黑暗

而我们能从太空中看到星球表面泛起的蓝光，日出时的美丽景色，都得于太阳光在大气中的散射

地球的大气中充斥着诸如空气分子，小水滴，尘粒等微小分子

当太阳光遇到这些分子时，会改变辐射方向，我们将这种现象叫做散射

（通常辐射能不会转化为热能，而只是改变方向）

 

 

 

在现实中，散射可能会发生无数次，介于性能，我们只实现一次散射就可以得到一个相对不错的效果

下面是单次散射的简单原理

我们首先思考一下太阳光如何通过大气进入我们的眼睛

假设我们所在的位置为 A点

现在给定一个观察方向，我们看向天空，视线与大气顶部相较于B点，而迎着我们视线方向的太阳光都会进入我们的眼睛

（实际上我们知道，我们的视野不是一条线，而是可以假定为一个视锥体，所以一副画面将会有许多条视线方向充斥着视锥体）

由于我们只计算单次散射，所以只需要计算刚好散射到我们视线方向的太阳光就可以了

但其实！有一件事情

这看似只有一次散射，但散射实际上时时刻刻都在发生，太阳光每向前传播一点就会遇到大气粒子然后发生散射

也就是说光线在这条路径上不断发生着散射，而散射就会造成太阳光向其他方向分散，就可能不会进入我们的眼睛

而单次散射的实际意义是 只计算光的一次转向，两次之后转向的光我们统一按照衰减计算

所以固定方向光线的衰减就是类似的过程，光线在传播过程中不断地被散射，越来越弱

 

因此太阳光进入我们眼睛的整个流程是这样的

 对于给定的视线方向，太阳光沿着固定方向入射，在 CP段 不断被散射衰减，然后与 AB段 形成无数个交点，我们将其成为Pi（i = 1，2....），在 P处 散射一部分光刚好迎着我们的视线方向进入我们的眼睛A点，PA段也不断被散射衰减

 

为什么要首先给定视线方向，其实这是适应shader的编写，因为我们写shader都是对顶点和片元操作

而模型通常是空心的，我们不会直接取到球中的一个点，而是球最外层的一个点，也就是我们的B点

而我们也知道A点（摄像机坐标），所以直接计算视线方向上散射进摄像机的光

 

在这里我会简单说一下shader计算的具体思路，方便理解公式的推导

我们知道太阳光和 AB段 形成了无数个交点 Pi，在每个交点 Pi 上都有一部分光被散射进摄像机

介于性能，我们不会全部计算，而是将 AB段 均分成 N 份，也就是采样 N 次，计算每个采样点散射进摄像机的光

具体思路为，对于 球模型 的一个顶点，我们将其定为 B点，计算 AB段 的长度，除以采样数量 N，得到 步长

我们还可以得到 向量AB，进而得到 方向向量AB，乘以 步长 就可以得到一步前进的向量

从 摄像机A点 开始，循环 N 次，每次计算一个 采样点P，然后乘以 步长 代表这一段的结果

再加上 步长 得到下一个采样点，累加得到一条视线方向最终的散射结果

然后 顶点 / 片元 着色器会遍历顶点输出最终的画面效果

下面是公式的简单推导（不涉及具体证明）

我们首先计算到达 P点 的光强，由于只计算在 P点 的一次转向所以 CP 段只计算光的衰减

设定 C点 的入射光强为 Ic，到达 P点 的光强为 Ip

$Ip=Ic\cdot T\left(CP\right)$

T(CP) 为光线从 C点 到 P点 的衰减方程，我们又将其叫做 透射函数

到达 P点 后，设定到达 A点 的光强为 Ia

$Ia=Ip\cdot S(CP,PA)\cdot T\left(PA\right)$

T(PA) 和上面一样，是 PA段 的透射函数

S(CP, PA) 是沿 CP方向 传播的光向 PA 方向散射的散射函数

就是从 C点 到达 P点的光 又有多少刚好散射到 PA 方向

简单地代换一下我们就可以知道从 C点 到 A点 的光强

$Ia=Ic\cdot S(CP,PA)\cdot T\left(CP\right)\cdot T\left(PA\right)$

而我们知道 AB段 上有 N 个 P点

所以对于一个顶点，最终到达 摄像机的光强

$Ia=Ic\cdot \sum _{i=1}^{N}S(CP,PA)\cdot T\left(CP\right)\cdot T\left(PA\right)\cdot ds$

也就是说实际上我们只需要知道透射和散射函数就可以了

于是就轮到解释一些非常简单的物理公式

因为我们知道透射衰减实际上也是散射，我们只需要将所有角度的散射加起来（求球面积分）就可以得到衰减

所以来了解这次涉及到的两种主要散射，瑞利散射（Rayleight）和米氏散射（Mie）

 

瑞利散射，指光线在遇到比自身波长小的粒子时发生的散射，而且会更多地散射波长短的光，比如蓝光，所以我们能够看到地球表面有一层薄薄的蓝光

米氏散射，指光线在遇到和自身波长差不多或大一些的粒子时发生的散射，对所有所有波长的光都会同等散射，并且大部分光都会沿着原来原来的方向传播

所以光路的形成得于米氏散射，比如我们熟悉的丁达尔效应

 

瑞利散射函数公式

$S(\lambda ,\theta ,h)=\frac{{\pi }^2{(n^2-1)}^2}{2}\cdot \frac{\rho \left(h\right)}{D_0}\cdot \frac{1}{{\lambda }^4}(1+{\cos }^2\theta )$

其中 n 为空气折射系数 1.00029

ρ(h) 为密度比率，在海平面为 1，随高度增加指数下降

D0 为标准大气的分子数密度 2.504*10^25

θ 为偏转角度，如下图

而我们对散射函数进行球面积分，就得到了一点的透射函数

$\beta (\lambda ,h)={\int }_0^{2\pi }{\int }_0^{\pi }S(\lambda ,\theta ,h)\cdot \sin \theta \cdot d\theta \cdot d\phi =\frac{8{\pi }^3{(n^2-1)}^2}{3}\cdot \frac{\rho \left(h\right)}{D_0}\cdot \frac{1}{{\lambda }^4}$

 

下面我们对透射函数和散射函数进行拆解，这是计划的一部分（为了方便计算和编写）

 我们可以看到透射函数只有两个变量，波长和高度，而高度仅在 ρ(h) 中，我们将 ρ(h) 单独提出来，剩下的就是

$\beta \left(\lambda \right)=\frac{8{\pi }^3{(n^2-1)}^2}{3}\cdot \frac{1}{H_0}\cdot \frac{1}{{\lambda }^4}$

而对于给定波长的光，β(λ) 是一个常量，也就是说透射函数实际上只与高度相关

ρ(h) 是高度 h 的大气密度 与 海平面大气密度 的比值

而大气密度与高度的曲线如下图

 

 

 （图片来自知乎@爱吃兔子的草）

呈现指数衰减

所以我们可以用一个指数函数来拟合这条曲线

$\rho \left(h\right)=e^{-\frac{h}{h_0}}$

h0 为均质大气高度，对于瑞利散射，h0 = 5800m，对于米氏散射，h0 = 1200m

 

注意，上面我们得到的是一点的透射函数，并不是T(CP)

因为在 CP段 发生了无数次散射，我们也仿照在 AB段 的做法，将 CP段 分成 N份，进行 N次 采样

设定每一个采样点为 Li (i = 1, 2 ...）

对于第一个采样点 L1

$I_{l1}=Ic-Ic\cdot {\beta }_1$

而对于之后每个采样点

$I_{li}=I_{li-1}-I_{li-1}\cdot {\beta }_i$

βi 为该点的透射函数

在数学上我们知道

${\left(e^x\right)}^{{}^{\prime }}=e^x$

对于 |x| 很小，有

$f\left(x\right)=f\left(0\right)+f^{{}^{\prime }}\left(0\right)\cdot x$

所以有

$e^x=e^0+e^{0}x=1+x$

所以有下面推导

$I_{l1}=Ic\cdot e^{-{\beta }_1}$

$I_{l2}=I_1\cdot e^{-{\beta }_2}$

$......$

$I_{li}=I_{li-1}\cdot e^{-{\beta }_i}$

最终代换可得

$I_{pN}=Ic\cdot e^{-\sum _{i=1}^N{\beta }_i}$

所以 T(CP) 已经出来了

$T\left(CP\right)=e^{-\sum _{i=1}^N\beta (\lambda ,h)\cdot ds}$

而我们知道 β(λ,h) 可以拆分成 β(λ) * β(h)

所以公式可以进一步推导为

$T\left(CP\right)=e^{-\beta \left(\lambda \right)\sum _{i=1}^N\rho \left(h\right)\cdot ds}=e^{-\beta \left(\lambda \right)\sum _{i=1}^N{e}^{-\frac{h}{h_0}}\cdot ds}$

而下面这一项，通常被称为 光学深度

$D\left(h\right)=\sum _{i=1}^N{e}^{-\frac{h}{h_0}}\cdot ds$

 

到此为止，我们已经完全知道 T(CP) 怎么求了，T(PA) 同理

 

还剩一项 S(λ, θ, h)

这一项我们同样进行拆解，

$S(\lambda ,\theta ,h)=\beta (\lambda ,h)\cdot \gamma \left(\theta \right)$

其中 β(λ,h) 我们在上面已经推导过了

而 γ(θ) 只需要除一下

$\gamma \left(\theta \right)=\frac{S(\lambda ,\theta ,h)}{\beta (\lambda ,h)}=\frac{3}{16\pi }\cdot (1+{\cos }^2\theta )$

γ(θ) 这一项通常叫做 相位函数，它代表了光线的散射方向

同理 β(λ,h) 代表了光线的散射强度

我们可以发现，咦！？我不是已经知道 视线方向（AB）了，光照方向也知道（CP），对于每个顶点的计算这两个方向是固定的

所以 θ 也是一个定值，γ(θ) 也是一个常数

 

到此为止，我们终于可以把 Ia 等号右边所有函数的外衣都扒下来了

$Ia=\sum _{i=1}^N{I}_{pi}\cdot ds$

$Ia=Ic\cdot \sum _{i=1}^{N}S(CP,PA)\cdot T\left(CP\right)\cdot T\left(PA\right)\cdot ds$

$Ia=Ic\cdot \sum _{i=1}^N\beta (\lambda ,h)\cdot \gamma \left(\theta \right)\cdot e^{-\sum _{i=c}^p\beta (\lambda ,h)\cdot ds}\cdot e^{-\sum _{i=p}^a\beta (\lambda ,h)\cdot ds}$

$Ia=Ic\cdot \beta \left(\lambda \right)\cdot \gamma \left(\theta \right)\sum _{i=1}^N\rho \left(h\right)\cdot ds\cdot e^{-\beta \left(\lambda \right)\cdot (\sum _{i=c}^p\rho \left(h\right)\cdot ds+\sum _{i=p}^a\rho \left(h\right)\cdot ds)}$

其中 ρ(h) 我们已经知道了

（为了方便理解，我将 T(CP) 和 T(PA) 中的求和上下边换成了 端点）

 

 

现在我们理解为什么要将 S(λ, θ, h) β(λ,h) 拆解开了

因为对于 给定光线，给定摄像机位置，也就是给定视线和光线方向的情况下，我们实际上需要计算的变量只有一个，光学深度

也就是只与 高度h 有关，我们计算出的一次光学深度不仅可以用于 S(λ, θ, h)，也可以用于 β(λ,h)

 

到此为止，瑞利散射的原理以及公式简单地解释完毕了，而米氏散射其实是一样的原理，唯一不同的是 相位函数 和 均质大气高度（h0）

我们只需要简单地更改函数的不同参数就可以了

 

这里将给出 米氏散射 相位函数的表达式

$\gamma \left(\theta \right)=\frac{\frac{3}{8\pi }(1-g^2)\cdot (1+{\cos }^2\theta )}{{(1+g^2-2\cdot g\cdot \cos \theta )}^{1.5}\cdot (2+g^2)}$

g 的取值通常为 -0.75 到 -0.99

而当 g = 0 时，我们神奇地发现，这不就是瑞利散射的相位函数么！？

（实际上相位函数有多个版本，我也没学太明白，只知道 g表示散射的对称性，大于零，表示大部分向后散射，小于零，表示大部分向前散射）

 

好了好了，现在我们已经基本明白太阳光进入眼睛这件事是怎么回事了

下面是代码部分，在注释中进行了一定的解释

Shader "Custom/ProceduralCustom"
{
    Properties
    {
        //不含大气的星球半径
        _EarthGroundRadius("Earth Ground Radius", Float) = 6.0
        //大气厚度
        _AtmosphereThickness("Atmosphere Thickness", Float) = 3.0
        //颜色
        _RayleightColor("Rayleight Color", Color) = (1, 1, 1, 1)
        _MieColor("Mie Color", Color) = (1, 1, 1, 1)
        //均质大气高度
        _RayleightHomogeneousAtmosphere("Rayleight Homogeneous Atmosphere", Float) = 0.2
        _MieHomogeneousAtmosphere("Mie Homogeneous Atmosphere", Float) = 0.1
        //亮度
        _Brightness("Brightness", Float) = 10
        //采样数量
        _SampleNum("Sample Number", Int) = 500
        _SampleLightNum("Sample Light Number", Int) = 20
        //Mie相位函数系数
        _g("Mie Phase g", Float) = -0.78
        //Mie微调系数
        _MieAdjust("Mie Adjust", Float) = 5.1
    }
    SubShader
    {
        Tags{"LightMode" = "ForwardBase"}
        Pass
        {
            Cull Off
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #include "UnityCG.cginc"
            #include "Lighting.cginc"
            #define PI 3.1415926536
            //变量声明
            float _EarthGroundRadius;
            float _AtmosphereThickness;
            float _Brightness;
            float4 _RayleightColor;
            float4 _MieColor;
            float _RayleightHomogeneousAtmosphere;
            float _MieHomogeneousAtmosphere;
            float _SampleNum;
            float _SampleLightNum;
            float _g;
            float _MieAdjust;
            struct a2v
            {
                float4 vertex : POSITION;
                float2 texcoord : TEXCOORD0;
            };
            struct v2f 
            {
                float4 pos : SV_POSITION;
                //保存了顶点的世界坐标
                float4 worldPos : TEXCOORD0;
                float2 uv : TEXCOORD1;
            };
            //得到与大气的交点
            float2 GetPoint(float3 Pos, float3 ViewDir, float3 Center, float Radius)
            {
                Pos -= Center;
                float a = 1;
                float b = 2.0 * dot(Pos, ViewDir);
                float c = dot(Pos, Pos) - (Radius * Radius);
                float delta = b * b - 4 * a * c;
                if(delta < 0) return -1;
                else 
                {
                    float FinDelta = sqrt(delta);
                    return float2(-b - FinDelta, -b + FinDelta) / 2 * a;
                }
            }
            //两种散射的大气密度比率
            float RayleightGetAtmoRho(float Height)
            {
                return exp(-max(Height, 0) / _RayleightHomogeneousAtmosphere);
            }
            float MieGetAtmoRho(float Height)
            {
                return exp(-max(Height, 0) / _MieHomogeneousAtmosphere);
            }
            //两种散射的相位函数
            float RayleightPhase(float AngleCos)
            {
                return (0.1875 / PI) * (1 + AngleCos * AngleCos);
            }
            float MiePhase(float AngleCos)
            {
                float _g2 = _g * _g;
                float upTmp = (1.0 - _g2) * (1 + AngleCos * AngleCos);
                float downTmp1 = 2.0 + _g2;
                float downTmp2 = 1.0 + _g2 - 2.0 * _g * AngleCos;
                float downTmp = pow(downTmp2, 1.5) * downTmp1;
                return (3.0  / 8.0 / PI ) * (upTmp / downTmp);
            }
            //顶点着色器
            v2f vert(a2v v)
            {
                v2f o;
                o.pos = UnityObjectToClipPos(v.vertex);
                o.uv = v.texcoord;
                //世界坐标
                o.worldPos = mul(unity_ObjectToWorld, v.vertex);
                return o;
            }
            float4 frag(v2f i) : SV_TARGET0
            {
                //包含大气的星球半径
                float AllEarthRadius = _EarthGroundRadius + _AtmosphereThickness;
                //星球球心坐标
                //float3 EarthCenter = mul(unity_ObjectToWorld, float4(0, 0, 0, 1)).xyz;
                float3 EarthCenter = float4(0, 0, 0, 1);
                //摄像机位置和视角方向
                float3 CameraPos = _WorldSpaceCameraPos;
                float3 CameraViewDir = normalize(i.worldPos - CameraPos);
                //得到与大气的两个远近交点
                //x分量为近点，y分量为远点
                float2 PointLength = GetPoint(CameraPos, CameraViewDir, EarthCenter, AllEarthRadius);
                //没有交点
                if (PointLength.y <= 0 && PointLength.x <= 0) return float4(0.0, 0.0, 0.0, 0.0);
                //我们限制正值，这样同时可以用于摄像机在大气内的情况
                PointLength.x = max(PointLength.x, 0);
                PointLength.y = max(PointLength.y, 0);
                //采样点
                float3 SamplePoint = CameraPos + CameraViewDir * PointLength.x;
                //计算采样距离和采样步长
                float SampleLength = PointLength.y - PointLength.x;
                float SampleStep = SampleLength / _SampleNum;
                float3 SampleStepDir = CameraViewDir * SampleStep;
                
                float RayleightInScatterSum = 0;
                float MieInScatterSum = 0;

                float3 Rayleight = float3(0, 0, 0);
                float3 Mie = float3(0, 0, 0);

                for(int i = 0; i < _SampleNum; i++)
                {
                    float SamplePointHeight = length(SamplePoint) - _EarthGroundRadius;
                    if(SamplePointHeight < 0) break;

                    float OneRayleightInScatter = RayleightGetAtmoRho(SamplePointHeight) * SampleStep;
                    float OneMieInScatter = MieGetAtmoRho(SamplePointHeight) * SampleStep;

                    RayleightInScatterSum += OneRayleightInScatter;
                    MieInScatterSum += OneMieInScatter;

                    //_WorldSpaceLightPos0为光源（仅平行光）在世界空间的位置
                    float3 LightDir = normalize(_WorldSpaceLightPos0);
                    float2 LightDirPointLength = GetPoint(SamplePoint, LightDir, EarthCenter, AllEarthRadius);

                    float SampleLightLength = LightDirPointLength.y;
                    float SampleLightStep = SampleLightLength / _SampleLightNum;
                    float3 SampleLightStepDir = LightDir * SampleLightStep;
                    //用采样点做起点
                    float3 SampleLightPoint = SamplePoint;
                    //外散射衰减
                    //是光线从 大气入射点 到 顶点 期间的衰减结果
                    float RayleightOutScatterSum = 0;
                    float MieOutScatterSum = 0;

                    for(int j = 0; j < _SampleLightNum; j++)
                    {
                        float SampleLightPointHeight = length(SampleLightPoint) - _EarthGroundRadius;
                        RayleightOutScatterSum += RayleightGetAtmoRho(SampleLightPointHeight);
                        MieOutScatterSum += MieGetAtmoRho(SampleLightPointHeight);
                        SampleLightPoint += SampleLightStepDir;
                    }
                    RayleightOutScatterSum *= SampleLightStep;
                    MieOutScatterSum *= SampleLightStep;
                    float3 AllScatter = exp(-((RayleightInScatterSum + RayleightOutScatterSum) * _RayleightColor 
                                            + (MieInScatterSum + MieOutScatterSum) * _MieColor * _MieAdjust)); 
                    Rayleight += OneRayleightInScatter * AllScatter;
                    Mie += OneMieInScatter * AllScatter;
                    SamplePoint += SampleStepDir;
                }
                float PointToCameraCos = -dot(CameraViewDir, _WorldSpaceLightPos0);
                float3 FinalColor = _LightColor0 * _Brightness * ((Rayleight * RayleightPhase(PointToCameraCos) * _RayleightColor)
                                                + (Mie * MiePhase(PointToCameraCos) * _MieColor * _MieAdjust));
                //或许需要开一个HDR
                return float4(min(FinalColor, float3(1, 1, 1)), 1);
                //return float4(FinalColor, 1);
                //return float4(0, 0, 0, 1);
            }
            ENDCG
        }
    }
    FallBack "Diffuse"
}

本次实现的效果没有优化，是实时计算，所以当采样数量很大的时候性能会明显不足，可以去别家大佬的文章中学习一下查找表，拟合函数等优化

参考了很多大佬的现成教程，仅有一些自己的思考，写到此处方便以后复习

参考资料

大气密度与海拔高度的关系 - 知乎 (zhihu.com)

游戏魔法编程：unity实现完整大气散射 - 知乎 (zhihu.com)

在Unity中制作Procedural Sky (Part 1-Atmospheric Scattering) - 知乎 (zhihu.com)