unity程序化天空盒--实现一次单次大气散射

一.前言

//采用的URP管线;

大气散射这一部分还是涉及比较多的知识点,因此决定单独写一篇文章，总结一下大气散射的相关知识点。

首先安利几篇文章,写的真的不错，个人感觉看这几篇足矣嘿嘿:

从零实现一套完整单次大气散射

https://zhuanlan.zhihu.com/p/238048524

文章一共四篇，看完应该对大气散射就会有一个很好的了解了。

由于这里我们要实现的是偏向于卡通化的渲染，因此我们做了两个修改:

1. 仅计算Mie散射的效果;

2. 光学距离部分只计算了D(PA¯)，省去了计算D(CP¯)；

关于卡通渲染的大气散射的实现,我们还可以参考这篇文章:

https://blog.csdn.net/qq_41835314/article/details/128409989

二.一些思路分析

2.1(总思路)主要是以我们的视角看向天空,即为人眼与天空相交的过程,如果射线被物体挡到，则返回为黑色,如果没有被挡到，则要返回一个单次大气散射的值;

2.2在计算dpa以及dpc的时候（这里我们忽略了dpc的计算），我们会在AB（看下图）之间取多个点，对于每一个点分别值并进行累加，就可以得出最终的dpa值;

2.3

下面开始代码的分析

三.代码框架

代码部分主要围绕下面两个方程：

3.1构建坐标系

首先我们会先构造坐标系,如下图，由于我们是以地球的顶部中心处作为我们的(0,0,0)的位置,因此我们不用进行世界坐标的调整转换，直接用Unity的世界坐标就可以得到正确的结果;

对于上图的解释，人眼就是在如图的A点处，我们的射线会与大气的球体相交于B点处,AB是射线的方向;

在主函数的计算之前,我们会先求出在不同位置处的大气散射系数函数,相位函数，以及在不同高度的大气密度系数函数:

3.2射线与大气球体的相交

首先是计算射线与球体的相交，思路主要是：

1.联立视线与球体的方程，得到一个二次函数的表达式，如果这一个二次函数的表达式有解，则说明存在交点，若不存在则返回-1;

2.这里我们返回的是射线方程中的+-t的值，之后我们用rayOrigin + t * rayDir(射线方程)就可以得到两个交点的三维坐标

关于数学的分析可以参考以下的文章:

射线与球体/三角面片求交、重心坐标、插值

https://blog.csdn.net/xiji333/article/details/108581804#:~:text=%E7%90%83%E4%BD%93%E6%B1%82%E4%BA%A4%20%E4%B8%8D%E5%A6%A8%E8%AE%BE%E7%90%83%E4%BD%93%E7%9A%84%E6%96%B9%E7%A8%8B%E4%B8%BA%20%28X%20%E2%88%92%20C%292%20%3D%20R2%20%EF%BC%8C%E5%85%B6%E4%B8%AD,%E2%88%97P%20%2BP%202%20%2BC%202%20%E2%88%92R2%20%3D%200

代码如下：

//--RaySphereIntersection --首先是两个点的方向t值的计算;
//返回的是±t，rayOrigin+t*rayDir就能得到两个交点,这两个交点就是光线从与大气相交的位置传入到我们的眼睛位置的两个点;float2 RaySphereIntersection(float3 rayOrigin,float3 rayDir,float3 sphereCenter,float sphereRadius){rayOrigin -= sphereCenter;//因为一开始传入的是一个点,因此这里我们减去原点的坐标得到一个向量;float a = dot(rayDir,rayDir);//这里是二次方程的a;float b = 2.0 * dot(rayOrigin,rayDir);//这里是二次方程的b;float c = dot(rayOrigin,rayOrigin)-(sphereRadius * sphereRadius);//下面是二次函数的求解;float d = b*b -4*a*c;if(d<0){return -1;//如果b^2 - 4ac小于01则返回,此时无解;}else{d = sqrt(d);return float2(-b-d,-b +d)/(2*a);//此处最后返回的的+-t,之后我们仅需要rayOrigin+t*rayDir就能得到两个交点;}}

接下来是大气散射的密度系数函数以及相位函数的方程;

3.3大气密度函数

主要是根据公式exp(-height / H);//其中H为海平面的大气密度,height为高度为height处的大气密度，这里返回的是一个比值;

代码如下:

//1.首先我们先计算大气密度--GetAtmosphereDensity--;//思路:这里我们会根据position以及地表中心距离减去地球的半径来求得高度,最后我们会输出dpc以及dpa,这里我们省去了dpc;void GetAtmosphereDensity(float3 position,float3 plantCenter,float3 lightDir,out float dpa,out float dpc){//计算公式为dpa = exp{-h/H};float height = length(position - plantCenter) - _PlanetRadius;//获得高度;dpa = exp(-height / _DensityScaleHeight);//求得此处的大气密度;dpc = 0;//dpc我们暂时省略;}

3.4相位函数

接下来是相位函数的计算;

相位函数代表的是在所有的散射方向中，某一个方向的比例,由于我们这里仅会计算Mie散射，因此先给出Mie散射的相位函数公式:

代码部分我们仅需要求出这个公式的值并与传入的scattering相乘即为得出最终在此方向上的sacttering；

//2.下面是相位函数的计算;--ApplyPhaseFunction--;//这里我们仅引入Mie相位函数!!;void AppltPhaseFunction(inout float scatterMie,float cosAngle){float g = _MieG;//注意这里的g在传入的时候要将值设定为0.76,这是一个固定的值!;float phase = (1.0 / (4.0 * PI)) * ((3.0 * (1.0 - g*g)) / (2.0 * (2.0 + g*g))) * ((1 + cosAngle * cosAngle) / (pow((1 + g*g - 2 * g * cosAngle), 3.0 / 2.0)));scatterMie *= phase;}

3.5大气散射主函数

到此处，前置的一些方程都已经完成，接下来就是主函数的计算，在这一个函数中我们会返回最中对应的颜色值;

代码如下:

//3.到这里先前工作都完成,接下来我们要进行的是大气散射的主函数!!!!--IntegrateInscattering--;//此函数计算返回值就是在我们的A点得到的总能量大小!!(用返回的颜色来进行表示);//rayStart为摄像机的起点,也就是理论图中的A点;//rayDir为射线的方向;//raylength为AB的长度;half4 IntegrateInscattering(float3 rayStart, float3 rayDir,float rayLength,float3 plantCenter,float3 lightDir,float samplerCount){//参考Mie散射的光的波长来进行赋值; float4 MieSct = float4(2.0f, 2.0f, 2.0f, 0.0f) * 0.00001f;//光波对应的颜色值，与散射系数的乘积，这里会在最后进行的I与散射系数的相乘中用到;_ExtinctionM = MieSct * MieExtinctionCoef;//Mie对应波长乘以散射系数;_ScatteringM = MieSct * MieScatterCoef;float3 stepVector = rayDir * (rayLength / samplerCount);float stepSize = length(stepVector);//这里会根据samplerCount分成几个点,此处是每两个点之间的长度;//下面是一些代传入参数的创建;float prevDPA =0;//在循环前算的的Density大气密度对应值;float prevTransmittance = 0;//透射率方程的参数,用于存储循环前的最后投射率值;float densityPA = 0;//PA的大气密度;float densityCP = 0;//CP的大气密度;float localDPA = 0;//地平线处的大气密度;float scatterMie = 0;//下面先求T(CP) * T(PA)的部分,具体看原理,主要就是先求出两个光学深度的和再乘上地平线处的大气密度;GetAtmosphereDensity(rayStart,plantCenter,lightDir,localDPA,densityCP);//最后我们得到的是摄像机位置的大气密度localDensity;//下面是PA的大气密度,这里进行加和;densityPA +=stepSize * localDPA;prevDPA = localDPA;//下面求T(CP) * T(PA)透射率方程的部分;float Transmittance = exp(-(densityPA + densityCP)*_ExtinctionM) * localDPA;//这里获取到透射率方程部分_ExtinctionM用于调节;prevTransmittance = Transmittance;//获取循环计算前的透射率比例值;//下面进行主函数的计算;求出的是BA上无数个(samplerCount个)点经过一次散射之后的光线能量,依次加入到最后的结果;for(int i = 0;i<samplerCount;i+=1.0){float p = rayStart + stepVector * i;//获取到每一个分点的位置;GetAtmosphereDensity(p,plantCenter,lightDir,localDPA,densityCP);//这里求出该点位置的大气密度;densityPA +=(localDPA + prevDPA) *stepSize/2.0;//加上两点间的密度值的累加,这里取平均;Transmittance = exp(-(densityPA+densityCP)* _ExtinctionM) * localDPA;scatterMie +=  (prevTransmittance +Transmittance) * stepSize /2.0;//同样进行一个累加,表示该点到A的透射率;//准备下一次的循环;prevTransmittance = Transmittance;prevDPA = densityPA;}//下面结束循环,还差三个量，相位值，散射系数，以及I;I与散射系数的乘积已经表示出来了;;AppltPhaseFunction(scatterMie,dot(rayDir,-lightDir.xyz));//角度为摄像机方向和光照反方向间的夹角;half3 lightInscatter = _ScatteringM * scatterMie;//乘上散射值与光强I的乘积，这一部分在前面有提及;;return half4(lightInscatter,1);//最后算出颜色值;}

关于累加项的dpa计算，我们可以在AB上取不同个数的相等距离的点p，之后对PA,PB分别求dpa dpb再相加，然后跳到下一个点上进行处理;

如图：

3.6片元着色器

最后我们在片元着色器里面传入参数即可(原文是通过cpp传入参数，这里参数我都直接放在Shader里里面)

 //Mie scatteringfloat3 scatteringColor = 0;float3 rayStart = float3(0,10,0);//人眼观察的位置;rayStart.y = saturate(rayStart.y);float3 rayDir = normalize(i.uv.xyz);float3 planetCenter = float3(0, -_PlanetRadius, 0);//地球的中心;float2 intersection = RaySphereIntersection(rayStart,rayDir,planetCenter,_PlanetRadius + _AtmosphereHeight);//获取两个t值;float rayLength = intersection.y;//获取正的t,也就是rayLength;intersection = RaySphereIntersection(rayStart, rayDir, planetCenter, _PlanetRadius);//这里判断遮挡是利用射线的方向与地球来进行相交，得到两个大于0的t值;
if (intersection.x > 0){//没有被遮挡则可以进行光线长度的计算;rayLength = min(rayLength, intersection.x*100);}//下面传入参数进行散射的计算;float4 inscattering = IntegrateInscattering(rayStart, rayDir, rayLength,planetCenter, -main_light.direction.xyz, 16);scatteringColor = _MieColor*_MieStrength * inscattering;//加入到最终的颜色里面;finalSkyColor +=scatteringColor;

3.7参数的调节

美术比较垃调得不太行。。。。可以参考下自己调节一下。。。。

跳出来的效果大概酱紫:

好了大概就是酱紫，基于物理的大气散射等做完天空盒再来尝试以下嘿嘿嘿......