Unity3D_实现大气散射

发表于2016-06-11
评论4 1.84w浏览

一、前言

  本文旨在与大家一起探讨学习新知识,如有疏漏或者谬误,请大家不吝指出,以下内容参考了GPU精粹2中第16章关于大气散射的知识,先来看下完成后的效果:








二、原理

天空为什么是蓝色的?赤橙黄绿青蓝紫,谁持彩练舞当空。毛大大的一句诗告诉我们,太阳的白光含有七种颜色。其中各个颜色的波长按顺序递减。红色波长最长,紫色波长最短。当阳光穿过大气层时,会发生散射,因为大气中存在着尺度比光的波长还小的粒子,这些粒子吸收长波长的光,反射短波长的光,所以大气层呈现蓝色。当然有人可能会问,紫光不是波长更短吗,为什么我们看到的不是紫光?之所以会这样有两个原因,一是人的眼睛对紫光不敏感,二是阳光光谱的峰值位于绿色与蓝色之间,就是说同一束阳光中蓝光的能量更多。

上述散射现象指的是小粒子的散射,称为Rayleigh散射(Rayleigh是个人名,他发现了这种现象),大气中还存在大粒子,大粒子对各个波长的散射比较均匀,所以散射出来的光是灰白色的。这个散射现象叫做Mie散射(还是一个人名)。当空气中大粒子比较多时,天空就是灰蒙蒙的,比如乌云还有雾霾。如果大粒子较少,Mie散射少,那么天空就是纯净的蓝色。当然要注意一点,在高空Rayleigh散射更强烈,低空则是Mie散射更多。

 

三、散射模型

   

1


    如图1,我们用两个球来表示地球和大气层,小的那个表示地球,大的表示大气层的最外层。当相机A在大气层中时,它观察到大气层的B点,那么我们就要来计算AB段上大气的散射结果。

下面直接给出大气散射的完整公式,我们不去细究公式怎么得来的,只要搞清楚每个参数的意义就OK了。

这里的公式是直接套用的GPU精粹2中的公式。

内向散射方程:

内向散射方程表示的是有多少散射光汇聚到了AB线段上并射向相机。因为大气中散射是无处不在的,而且散射的方向是四面八方的。就像图2所画的那样(图丑不要介意)。

2


可以看出方程里有积分,在shader中,积分一般都是通过raymaching方式来计算,即把AB段分成很多个小线段,然后通过累加各个线段的积分结果来近似求解方程。

内向散射中,是阳光的波长;表示的是阳光的光强,你可以用一个float的类型定义它,也可以用float3是散射常数,一般来说Rayleigh散射的散射常数跟波长有关,通常要除以波长的四次方,而Mie散射的散射常数可以认为与波长无关;是相位函数,表示光朝某个方向散射的概率,是指光与某个方向的夹角,在这里指图三中的夹角,g是常数,Rayleigh散射的常数g可以等于0,而Mie散射的常数g一般在-0.75~-0.99之间,g不能等于1-1hP3点的高度,在这里我们把高度归一化,即认为在地平线上高度为0,在大气层最外层高度为1H0是标尺高度,表示在多高的地方大气的密度是平均密度;t函数是外向散射方程。

 

外向散射方程:

外向散射方程表示光在线段中散射出去多少。各个参数在内向散射中都有说明到。

然后我们考虑怎么在GPU中实现上述方程的求解。

1


上面有说过方程中有积分项的求解方法。那么如图1所示,我们把AB分成5段,然后考虑在P3点处怎么计算内向散射方程的结果(图画的有点不正确,P1P5指的是各个小线段的中点)。首先计算t函数,即外向散射方程,这个方程中又有积分,那么我们就再对AP3P3C进行分段求解;然后求出两个t函数的解带入内向散射方程中,得到s3,将s3乘以小线段的长度,这就近似表示这一个小线段的内向散射结果,其他线段依次求解,将各个线段的结果累加就得到了方程的近似解。

其他参数都很容易求得,有一个比较复杂的是C点的求解,当然通过列出球的公式以及射线的公式联立求解也很简单。

  ——球的方程

(a, b, c)是球心坐标,r是球半径。

   ——射线方程

其中P是射线起点,即相机的位置,d是距离,Ray是射线方向。

联立两个方程可得一个关于d的一元二次方程,求解公式相信大家都很清楚。


四、实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
//P and raydefine the line; O is the sphere's center point
float3 CalculateIntersection(float3P, float3 ray, float3 O, float r)
{
         float3PO = O - P;
         floatlPO = length(PO);
         floatB = 2.0*dot(PO, ray);
  
         floatC = B*B - 4*(lPO*lPO-r*r);
         floatdet = sqrt( max(0, C) );
  
         //elsed: .5*(B+det)
         floatd = 0.5*(B-det);
         float3intersect = P+d*ray;
  
         returnintersect;
}
//g=0
floatPhaseRayleigh(float costheta)
{
         return0.75*(1+costheta*costheta);
}
floatPhaseMie(float costheta, float g, float g2)
{
         return1.5*( (1.0-g2) / (2.0+g2) ) *
(1.0+costheta*costheta) / pow(1.0+g2-2.0*g*costheta, 1.5);
}
  
v2f vert(appdata v)
{
         v2fo;
         o.vertex= mul(UNITY_MATRIX_MVP, v.vertex);
         o.uv= TRANSFORM_TEX(v.uv, _MainTex);
         //
         float3lightDir = normalize(_WorldSpaceLightPos0.xyz);
         float3worldPos = mul(_Object2World, v.vertex);
         float3ray = worldPos - _WorldSpaceCameraPos.xyz;
         floatdepth = length(ray);
         ray= ray/depth;
         o.viewDir= _WorldSpaceCameraPos.xyz-worldPos;
                                     
         //
float sampleLength = depth/_SampleNum;
         float3sampleRay = ray*sampleLength;
         float3samplePoint = _WorldSpaceCameraPos.xyz + sampleRay*0.5;
         float3inScattering = float3(0,0,0);
         float3outScatteringAB = float3(0,0,0);
         floatinvOuterH = 1.0/(_OuterRadius-_Sealevel);
                                     
         for(inti=0; i<_SampleNum; i++)
         {
                   float3crossPoint = CalculateIntersection(samplePoint, -lightDir, _PlanetPos,_OuterRadius);
                   //APout scattering
                   float3APRay = (samplePoint-_WorldSpaceCameraPos.xyz)/_SampleNum;
                   floatAPRayLength = length(APRay);
                   float3sampleAP = _WorldSpaceCameraPos.xyz + APRay*0.5;
                   float3outScatteringAP = float3(0,0,0);
  
                   for(intm=0; m<_SampleNum; m++)
                   {
                                                                                                                outScatteringAP+=exp(-(length(sampleAP-_PlanetPos)-_Sealevel)*invOuterH/_H0)*APRayLength*invOuterH;
                            sampleAP+= APRay;
                   }
                   outScatteringAP= 4*3.14*(_Km+_Kr*_Wave)*outScatteringAP;
  
                   //PCout scattering
                   float3PCRay = (crossPoint - samplePoint)/_SampleNum;
                   floatPCRayLength = length(PCRay);
                   float3samplePC = samplePoint + PCRay*0.5;
                   float3outScatteringPC = float3(0,0,0);
  
                   for(intn=0; n < _SampleNum; n++)
                   {
                            outScatteringPC+=exp(-(length(samplePC-_PlanetPos)-_Sealevel)*invOuterH/_H0)*PCRayLength*invOuterH;
                            samplePC+= PCRay;
                   }
                   outScatteringPC= 4*3.14*(_Km+_Kr*_Wave)*outScatteringPC;
  
                   inScattering+= exp(-(length(samplePoint-_PlanetPos)-_Sealevel)*invOuterH/_H0)*sampleLength*invOuterH*exp(-outScatteringAP-outScatteringPC);
                   outScatteringAB= outScatteringAP;
                   samplePoint+=sampleRay;
         }
  
         inScattering= inScattering*_ESun*_Km;
         o.scattering= inScattering;
         o.attenuate= inScattering*_ESun*_Kr*_Wave;
  
         UNITY_TRANSFER_FOG(o,o.vertex);
         returno;
}
                            
float4 frag (v2fi) : SV_Target
{
         //sample the texture
         float4col = tex2D(_MainTex, i.uv);
         float3lightDir = normalize(_WorldSpaceLightPos0.xyz);
         floatcostheta = dot(normalize(i.viewDir), lightDir);
         floatphaseMie = PhaseMie(costheta, _G, _G2);
         floatphaseRayleigh = PhaseRayleigh(costheta);
                                     
         col.xyz= i.scattering*phaseMie + i.attenuate*phaseRayleigh;
         col.a= col.b;
         //apply fog
         UNITY_APPLY_FOG(i.fogCoord,col);                                 
         returncol;
}


这里给出的是大气层的shader,代码稍微有点复杂,不过对照上述第三小节来看应该还是能看懂的。这里的代码没有做任何优化,为的就是希望大家能看懂。地面的散射代码跟大气层的稍微有点不一样,地面反射的光在射向相机时也会发生散射,所以需要计算外向散射。

即:

其中表示地面的反射的光的总量。


五、优化

GPU精粹2中,有对优化进行了大量的描述,并且给出了优化后的演示代码。其主要思路就是预计算外向散射,然后得出张查找表,再根据需要从查找表中获取。这是第一步优化。更深入的优化是根据查找表画出值的变化曲线,然后使用方程来拟合曲线,这样就用方程来替代了查找表。更详细的内容请大家自行查看相关资料,由于对曲线拟合这块不熟悉,这里就不做细讲。 

完整代码大家可以下载附件进行查看。

 

如社区发表内容存在侵权行为,您可以点击这里查看侵权投诉指引