什么是大气散射
大气散射模拟太阳光在大气层穿行时发生的一系列复杂的光学现象，进而模拟出大自然中令人惊艳的视觉效果。UE提供了一个现成的大气散射组件，大致效果如下所示：
传统大气的渲染是使用类似 phong 模型的经验公式去逼近，虽然结果还行，但参数对于艺术家不好理解，且不是基于物理的
为什么需要大气散射
很直接的一个原因：为了好康！上图可以看出加了大气散射使得整体氛围拉满
理论
其中的推导部分可以看大哥的文章
散射、衰减、吸收
在谈论物理渲染部分都绕不开这三货：散射、衰减、吸收
散射：和次表面散射的散射是同一个意思——光线打到物体表面，会出现反射、折射，且带有能量损失
衰减：光线在空气中行进也会伴有能量损失
吸收：光线打到物体表面时，部分散射，部分被吸收
大气中根据颗粒的直径大小，主要存在两种散射——瑞利散射和米氏散射
- 瑞利散射（颗粒直径远小于波长）：不同波长的光损失的能量各有不同，但是散射向各个方向的能量相差不大(
各向同性)。瑞利散射逃逸较多的蓝光，这是天空呈现蓝色的原因 - 米氏散射（粒子直径和波长相当）：不同波长的光损失的能量相同，但是不同的散射方向出射的能量有着巨大的差异(
各向异性)。米氏散射方向性较强，是太阳光晕的成因
大气
- 大气的组成：大气里包含许多不同的颗粒，这些颗粒有大有小，当光线打到不同的颗粒时，会发生不同的散射。当打到小颗粒时，会发生Rayleigh 散射；当打到大颗粒时，会发生Mei散射
- 大气的密度：大气的密度并不是均匀的，
密度大小与海拔有关，密度越大散射强度越强
散射方程
瑞利散射和米氏散射都有着不同的散射方程，用以描述入射光散射到特定方向上的能量百分比
- 主要形式：S(λ, θ, h).除此之外皆为常量
-
定义：波长λ的光，在海拔h处，散射到θ角度的能量百分比。散射方程真正表示的是散射到相机的能量大小（因为散射会朝各个方向，但渲染时只关心散射到相机方向的百分比）
先着手相对简单的瑞利散射，一步步地来解读散射方程，搞明白后再看更加复杂的米氏散射
瑞利散射
瑞利散射方程：S(λ, θ, h) = \frac{Π^2(n^2 - 1)^2}{2} \frac{ρ(h)}{N} \frac{1}{λ^4} (1 + cos^2θ)
散射系数
- 定义：散射到某方向的百分比，更准确地说散射系数描述了光散射后的能量损失百分比
简单来说，散射方程描述了入射光散射到特定方向上的能量百分比，这个百分比就由散射系数决定
-
主要形式：β = (λ, h)
散射系数只和波长、海拔高度有关，与方向无关
-
瑞利散射系数：β = (λ, h) = \frac{8Π^3 (n^2 - 1)^2}{3} \frac{ρ_r(h)}{N} \frac{1}{λ^4}
这里先展示瑞利散射系数的主要形式，细节等后面再介绍，不然人要晕的
-
简单计算
瑞利散射虽然比米氏散射更简洁，但计算量依然不小，在实时渲染中一股脑的计算是不可取的，需要一定的手段来简化
通常在简化时，需要把散射系数分为β(λ) ρ(h)，其中β(λ)表示海拔高度为0的散射系数，即\frac{8Π^3 (n^2 - 1)^2}{3} \frac{1}{N} \frac{1}{λ^4}
在实际计算时，大多设定红绿蓝的波长为440、550、680nm，即β(440,550,680) = (33.1, 13.5, 5.8) * 10^{-6}。实际上，这个值是散射系数的最大值，当海拔高度越高，散射系数越小
为什么要拆分开呢？因为可以看出当β(λ)与海拔高度无关时，仅仅是一个常数，那么可以提前在cpu端将β(λ)传给gpu，以减少gpu压力
剩下的只需在GPU中计算ρ(h)，但ρ(h)的计算并不简单，因为它代表着海拔高度对散射的影响，需要对他计算大量积分
相位函数
目前已经知道散射方程中的波长与海拔高度了，只剩下计算指定方向占散射方向的百分比，而相位函数就在做这个工作
- 主要形式：γ(θ)
-
瑞利散射相位函数：γ(θ) = \frac{3}{16Π}(1 + cos^2θ)
根据该函数不难看出，光线很难沿着90°方向散射，很容易沿着原方向或反方向散射
大气密度比例函数ρ(h)
ρ(h)意为大气密度比例方程，表示海拔为h处的大气密度和海平面处的大气密度的比值
- 函数形式：ρ(h) = \frac{density(h)}{density(0)}
density(0)可以测算出来，但density(h)需要计算且难度不小，因为大气的组成非常复杂，由不同压力、密度、温度的几层组成。好在瑞丽散射和米氏散射发生在较低的大气层，在这里大气密度随着高度指数衰减，可以使用一条曲线近似模拟density(h)
但在实际计算时，不会用到上述函数，而是使用近似——ρ(h) = exp(-\frac{h}{H})，其中H是基准高度（固定值），瑞利散射为8500米
能量衰减
已经知道散射系数了，那么假设光线初始强度为I_0，它在一处散射系数为β的均匀介质行驶x距离后，剩余的能量I_1应为多少？
如果只发生一次散射，那I_1肯定为I_1 = I0(1 - β)
但在大气中是存在许多许多颗粒的，大概率是会打中多个颗粒，导致多次散射的，因此需要一个公式来表示多次散射的结果
能量衰减公式：I_1 = I_0 \space exp(-β(λ,h)x)
大气透射率
上处能量衰减是建立在均匀介质的基础上计算的，但现实中大气的介质并非均匀的
如下图所示，光线从c点进入大气，此处能量为I_c，一直行进到点P，在途中发生了多次散射，需要求点P的能量I_p的值
- 为了解决该问题需要引出大气透射率方程，它描述了光线在大气中行进一段距离后， 剩余的能量百分比
-
公式：T_{CP} = \frac{I_P}{I_C}
-
推导过程
- 根据能量衰减公式，不难得出：
I_Q = Ic * exp(-β(λ,h_0) * Length(CQ))
I_P = I_Q * exp(-β(λ,h_1) * Length(QP))
- 那么：I_P = [Ic * exp(-β(λ,h_0) * Length(CQ))] * [exp(-β(λ,h_1) * Length(QP))]
-
但在实际计算中，会存在多个CQ、QP这样的距离，且它们的距离值都是固定的，那么就可以进行积分推广：I_P = I_c * exp(-β(λ, h_0) - β(λ, h_1) ds)
$I_P = I_c exp(- \sum_{Q∈CP} β(λ, h_Q) ds)$
其中，点Q看作CQ上许许多多用于积分的点
- 那么T_{CP} = exp(- \sum_{Q∈CP} β(λ, h_Q) ds)
- 根据能量衰减公式，不难得出：
- 优化
- 把散射系数公式代入T_{CP}，可得：T_{CP} = exp(- \sum_{Q∈CP} \frac{8Π^3 (n^2 - 1)^2}{3} \frac{ρ(h_Q)}{N} \frac{1}{λ^4} ds)
-
移动一下：T_{CP} = exp(- \frac{8Π^3 (n^2 - 1)^2}{3}\frac{1}{λ^4} \frac{1}{N} \sum_{Q∈CP} ρ(h_Q) ds)
可以看到前面的是常量，尾部的\sum_{Q∈CP} ρ(h_Q) ds即为积分项，计算时使用Raymarch即可
为了简便，这里将\sum_{Q∈CP} ρ(h_Q) ds看作D(CP)（译作OpticalDepth 光学深度，意为一段距离， 大气密度比例ρ(h)的累加），可得T_{CP} = exp(-β(λ) D(CP))
米氏散射
米氏散射方程：S(λ, θ, h) = Π^2(n^2 - 1)^2 \frac{ρ(h)}{N} \frac{1 - g^2}{2 + g^2} \frac{1 + cos^2θ}{(1 - g^2 - 2gcosθ)^{\frac{3}{2}}}
相位函数
γ(θ) = \frac{1 - g^2}{4Π(1 + g^2 - 2gcos(θ))^{\frac{3}{2}}}
大气密度函数
ρ(h) = exp(-\frac{h}{H}), H = 1200.0f
大气散射模型
事实上光线进入大气后并不会以一条直线一直向前行进，而是碰到颗粒，向多个方向散射多条光线，最终进入相机
这就存在单次散射和多次散射的情况
- 单次散射：光线进入大气后，第一次打中颗粒导致散射，最终进入相机
- 多次散射：光线进入大气后，两次及以上打中颗粒导致散射，最终进入相机
因为多次散射非常复杂，对于实现渲染来说是不可能的，且单次散射的渲染效果已经不错了，所以这里只介绍单次散射方程的推导
单次散射方程
根据先前提到的散射方程，不难推出：
I_{PA} = I_P S(λ, θ, h) T(PA)
I_{PA} = I_C \space T(CP) \space S(λ, θ, h) \space T(PA)
I_{PA} = I_C \space S(λ, θ, h) \space T(CP) \space T(PA)
但AB上不止存在一个点P使光线发生散射，极有可能有多个，那么就需要积分了：
I_A = I_{PA} \sum_{P∈AB} ds
最终得到单次散射模型的渲染方程：I_A = I_C \space S(λ, θ, h) \space T(CP) \space T(PA) \sum_{P∈AB} ds
优化
散射方程
S(λ,θ,h) = β(λ,h) γ(θ) = β(λ) ρ(h) γ(θ)
因为β(λ)表示散射系数，与海拔高度无关，是个常量，且γ(θ)表示相位函数，是个常量。所以可以将β(λ) γ(θ)这两货从散射方程中提出来，可得：I_A = I_C \space β(λ) γ(θ) \space T(CP) \space T(PA) \sum_{P∈AB} ρ(h) ds
大气透射率
接下来优化T(CP) T(PA)
因为大气散射率公式为T_{CP} = exp(-β(λ) D(CP))
那么可拆分为：
exp(-β(λ) D(CP)) \space exp(-β(λ) D(PA))
exp(-β(λ) (D(CP) + D(PA)))
其中OpticalDepth D(CP)=\sum_{Q∈CP} ρ(h_Q) ds = \sum_{Q∈CP} exp(-\frac{h_Q}{H}) ds
积分
目前得到的渲染方程如下
I_A = I_C \space β(λ) γ(θ) \space T(CP) \space T(PA) \sum_{P∈AB} ρ(h) ds
T(CP) T(PA) = exp(-β(λ) (D(CP) + D(PA)))，其中光学深度D(CP) = \sum_{Q∈CP} exp(-\frac{h_Q}{H}) ds
根据这个方程可以将其分为两部分，一部分是常量，由CPU端提前传入；另一部分是需要实时计算的
真正需要优化的是需要实时计算的部分，目前存在三个积分项
- 第一个积分：AB上存在许多颗粒产生散射
- 第二个积分：大气透射率，OpticalDepth CP
- 第三个积分：大气透射率，OpticalDepth PA
说人话就是取AB上的一个颗粒点P，计算CP的大气透射率进并累加，再计算PA的大气透射率并累加，并计算AB上其他所有的颗粒
这里还阔以进一步优化为两个积分，也就是优化D(PA)，因为点P是在AB上的，那么上一次积分得到的D(P_{i-1}A)，可以直接用于下次积分D(P_iA)
实现
思路
主要思路：相机看向天空，若射线被物体遮挡，则返回0；若未被挡到，则返回单次大气散射的值
后续实现主要围绕渲染方程展开，具体形式如下：
I_A = I_C \space β(λ) γ(θ) \space T(CP) \space T(PA) \sum_{P∈AB} ρ(h) ds
T(CP) T(PA) = exp(-β(λ) (D(CP) + D(PA)))，其中光学深度D(CP) = \sum_{Q∈CP} exp(-\frac{h_Q}{H}) ds
这里为了方便使用URP的Render Feature实现，需要注意大气散射的pass在渲染管线中的位置是位于半透明之前的
重构世界坐标
这个无需多腌，很简单
// positionVP : uv
// depth : deviceDepth
float3 ReBuildPosWS(float2 positionVP, float depth)
{
////////////////////
// 变换到 NDC space
////////////////////
float3 positionNDC = float3(positionVP * 2.f - 1.f, depth);
#if defined (UNITY_UV_STARTS_AT_TOP)
positionNDC.y = - positionNDC.y;
#endif
////////////////////
// 变换到 world space
////////////////////
float4 positionWS = mul(UNITY_MATRIX_I_VP, float4(positionNDC, 1.f));
positionWS.xyz /= positionWS.w;
return positionWS;
}
在render feature中的RenderPassEvent设为RenderPassEvent.BeforeRenderingTransparents
Debug可以得到如下结果：
确定坐标系
为了配合其他渲染效果，如雾效，需要确定一个适合的坐标系，避免额外的计算
因为大气层是包围地球的，所以将地球顶部作为原点——使用unity的世界坐标即可
射线和大气层求交
假设相机在A点，它发出的射线与大气层的交点为C点。已知A点，需要求交点C，使用一元二次方程即可
/// ------------------ In
// rayOrigin : 相机位置
// rayDir : 相机朝向
// sphereCenter : 球体中心点
// sphereRadius : 球体半径
/// ------------------ Out
// return : rayLength
float2 RaySphereIntersection(float3 rayOrigin, float3 rayDir, float3 sphereCenter, float3 sphereRadius)
{
float3 oc = rayOrigin - sphereCenter;
float a = dot(rayDir, rayDir);
float b = 2.f * dot(rayOrigin, oc);
float c = dot(oc, oc) - Pow2(sphereRadius);
float delta = Pow2(b) - 4 * a * c;
if(delta < 0.f)
{
return -1.f;
}
else
{
return float2(-b - sqrt(delta), -b + sqrt(delta)) * rcp(2 * a);
}
}
求解T(CP) T(PA) \sum_{P∈AB} ρ(h) ds
接下来一步一步地求解大气渲染方程：I_A = I_C \space β(λ) γ(θ) \space T(CP) \space T(PA) \sum_{P∈AB} ρ(h) ds
首先来解T(CP) T(PA) = exp(-β(λ) (D(CP) + D(PA)))中的D(CP) + D(PA)
因为D(CP) = \sum_{Q∈CP} exp(-\frac{h_Q}{H}) ds，主要是求解h_Q
/// ------------------ In
// currPos : 采样点P
// lightDir : 光源方向
/// ------------------ Out
// opticalDepthCP : CP的光学深度
void LightSampling(float3 currPos, float3 lightDir, float3 planetCenter,
inout float2 opticalDepthCP)
{
float3 rayOrigin = currPos;
float3 rayDir = -lightDir;
float2 intersection = RaySphereIntersection(rayOrigin, rayDir, planetCenter, _PlanetRadius + _SkyAtmosphereHeight);
float3 rayEnd = rayOrigin + rayDir * intersection.y;
float3 step = (rayEnd - rayOrigin) / _SampleCounts;
float stepLength = length(step);
float2 density = 0;
for(float s = 0.5; s < _SampleCounts; ++s)
{
float3 pos = rayOrigin + step * s;
float height = abs(length(pos - planetCenter) - _PlanetRadius);
float2 currDensity = exp(-(height.xx / float2(_SeaLevelHeight_R, _SeaLevelHeight_M))); // 瑞丽和米氏都实现
density += currDensity * stepLength;
}
opticalDepthCP = density;
}
然后再求DCP + DPA
/// ------------------ In
// currPos : 采样点P
// lightDir : 光源方向
/// ------------------ Out
// opticalDepth CP & opticalDepth PA
void GetAtmosphereOpticalDepth(float3 currPos, float3 lightDir, float3 planetCenter,
inout float2 dpa, inout float2 dpc)
{
float height = length(currPos - planetCenter) - _PlanetRadius;
dpa = exp(- height.xx / float2(_SeaLevelHeight_R, _SeaLevelHeight_M));
LightSampling(currPos, lightDir, planetCenter, dpc);
}
因为求解DPA时可以做优化，累加之前的DPA（DPA_{i-1}），所以每次只需计算当前的DPA
接下来就可以解T(CP) \space T(PA) \sum_{P∈AB} ρ(h) ds这个方程了
/// ------------------ In
// localDensity : rho(h)
/// ------------------ Out
// localInscatter_R : 瑞利散射相关
// localInscatter_M : 米氏散射相关
void CalcLocalInscatter(float2 localDensity, float2 D_PA, float2 D_CP,
inout float3 localInscatter_R, inout float3 localInscatter_M)
{
float2 D_CPA = D_PA + D_CP;
// _Extinction_R:瑞丽散射系数
// _Extinction_M:米氏散射系数
float3 T_R = D_CPA.x * _Extinction_R;
float3 T_M = D_CPA.y * _Extinction_M;
float3 extinction = exp(-(T_M + T_R));
localInscatter_R = localDensity.x * extinction;
localInscatter_M = localDensity.y * extinction;
}
最后只需求解前三个系数I_C \space β(λ) γ(θ)
I_C即光源能量，而β(λ)是散射系数，一个常量。因此只需求解相位函数γ(θ)
/// ------------------ In
// cosAngle : 散射角
/// ------------------ Out
// scatter_R : 乘上瑞丽相位函数
// scatter_M : 乘上米氏相位函数
void GetPhaseFun(float cosAngle,
inout float3 scatter_R, inout float3 scatter_M)
{
float phase = 3.f * rcp(16.f * PI) * (1.f + Pow2(cosAngle));
scatter_R *= phase;
float g = _MieG;
float g2 = Pow2(g);
phase = rcp(4.f * PI) * 3.f * (1.f - g2) * rcp(2.f * (2.f + g2)) * (1.f + Pow2(cosAngle)) * rcp(pow(1.f + g2 - 2 * g * cosAngle, 1.5f));
scatter_M *= phase;
}
最后整合这些函数
/// ------------------ In
// rayOrigin : 相机位置
// rayDir : 相机朝向
// rayLength : AB长度
// planetCenter : 地球中心位置
// lightDir : 光源方向
// sampleCounts : AB采样次数
// distanceScale : 世界坐标尺寸
/// ------------------ Out
// intensity : 散射后的能量
// extinction : 透射率
void CalcSkyAtmosphereIntensity(float3 rayOrigin, float3 rayDir, float rayLength, float3 planetCenter, float3 lightDir, float sampleCounts, float distanceScale,
inout float4 intensity, inout float4 extinction)
{
float3 step = rayDir * (rayLength / sampleCounts);
float stepLength = length(step) * distanceScale;
float2 D_PA = 0;
float3 scatter_R = 0;
float3 scatter_M = 0;
float2 localDensity = 0;
float2 D_CP = 0;
// 存储上一个DPA
float2 preLocalDensity = 0;
float3 preLocalInscatter_R, preLocalInscatter_M = 0;
GetAtmosphereOpticalDepth(rayOrigin, lightDir, planetCenter, preLocalDensity, D_CP);
CalcLocalInscatter(preLocalDensity, D_PA, D_CP, preLocalInscatter_R, preLocalInscatter_M);
[loop]
for(float i = 1; i < sampleCounts; ++i)
{
float3 currPos = rayOrigin + i * step;
GetAtmosphereOpticalDepth(currPos, lightDir, planetCenter, localDensity, D_CP);
D_PA += (localDensity + preLocalDensity) * (stepLength / 2.f);
float3 localInscatter_R, localInscatter_M = 0;
CalcLocalInscatter(localDensity, D_PA, D_CP, localInscatter_R, localInscatter_M);
scatter_R += (localInscatter_R + preLocalInscatter_R) * (stepLength / 2.f);
scatter_M += (localInscatter_M + preLocalInscatter_M) * (stepLength / 2.f);
preLocalInscatter_R = localInscatter_R;
preLocalInscatter_M = localInscatter_M;
preLocalDensity = localDensity;
}
GetPhaseFun(dot(rayDir, -lightDir), scatter_R, scatter_M);
intensity = float4((scatter_R * _Scattering_R + scatter_M * _Scattering_M) * _LightColor, 1.f);
extinction = exp(-(D_CP.x * _Extinction_R + D_CP.y * _Extinction_M));
extinction.w = 0.f;
}
注意事项
为了使得大气散射不出错，需要保证当物体遮挡skybox，该处采样值返回0，没有被遮挡则计算大气散射值
PSOutput SkyAtmosphere_PS(PSInput i)
{
PSOutput o = (PSOutput)0;
// 重建世界坐标
float2 uv = (i.positionCS.xy - 0.5f) * _RTSize.zw;
float deviceDepth = _CameraDepthTexture.SampleLevel(sampler_LinearClamp, uv, 0).r;
float3 posWS = ReBuildPosWS(uv, deviceDepth);
// 相机相关
// A点即为相机起点
float3 rayOrigin = GetCameraPositionWS();
float3 rayDir = posWS - rayOrigin;
float rayLength = length(rayDir);
rayDir /= rayLength;
float3 planetCenter = float3(0, -_PlanetRadius, 0);
// 被遮挡返回黑色
if(deviceDepth > 0.000001f)
{
rayLength = 1e20;
}
float2 intersection = RaySphereIntersection(rayOrigin, rayDir, planetCenter, _PlanetRadius + _SkyAtmosphereHeight);
rayLength = min(rayLength, intersection.y);
float4 skyAtmosphereExtinction = 0;
float4 skyAtmosphereIntensity = 0;
if(deviceDepth < 0.000001f)
{
CalcSkyAtmosphereIntensity(rayOrigin, rayDir, rayLength, planetCenter, _MainLightDir, _SampleCounts, _DistanceScale, skyAtmosphereIntensity, skyAtmosphereExtinction);
}
o.color = skyAtmosphereIntensity + _Source_RT.SampleLevel(sampler_LinearClamp, uv, 0);
return o;
}
效果
优化
其实每次求DPA时，都会采样很多段CP，开销挺大的，笔者测试了下，性能耗时得有0.8ms。若不是基于物理的大气散射，可以把对CP的积分砍掉，只计算DPA
void GetAtmosphereOpticalDepth(float3 currPos, float3 lightDir, float3 planetCenter,
inout float2 dpa, inout float2 dpc)
{
float height = length(currPos - planetCenter) - _PlanetRadius;
dpa = exp(- height.xx / float2(_SeaLevelHeight_R, _SeaLevelHeight_M));
//LightSampling(currPos, lightDir, planetCenter, dpc);
}
这样性能耗时仅有0.1ms，也可以得到不错的效果
代码仓库
https://github.com/chenglixue/unity-SkyAtmosphere/tree/main
Comments | NOTHING