|
|
前言：
本文旨在对于大气散射中的原理推导与代码实现，在阅读本文之前，笔者希望你有对于微积分的一定知识，限度大概在于大学高数水平，其余所有前置知识理论，笔者会完全帮助你理解并进行理论推导。本文分为两个大的模块，原理和实现，在原理篇中会对大气散射中的散射理论：瑞利散射与米氏散射，同时还有学习散射的前置理论知识，例如麦克斯韦方程组，伴随勒让德方程多项式等等进行推导，帮助读者完整了解大气散射理论。而在实现篇中，本文会使用unity引擎进行大气散射的具体实现，同时介绍了ComputerShader，也说明了Lut图的计算原理。同时在文章最后，会给出大气散射的开源工程。让我们开始这一场冒险吧
另外，笔者其实本文删减了有一半的内容，因为想说的两部分内容：全动态实时天空和大气散射中的数值积分方式的对比。写了有本文的一半以上内容，发现会让本文变成裹脚布，也不利于读者吸收，于是后续会开新篇章进行说明。希望各位阅读后能够留下宝贵的意见，十分感谢
理解大气散射原理
大气散射实际是一个对外散射和对内散射的积分贡献问题。
简要观察大气散射过程：
在大部分的渲染过程中都是对于物体表面的着色描述，但是对于一束光线进入物体内部进行反复偏转得到一系列的结果的过程就很少有描述，但是大气散射就是这样一种着色过程。从这个描述中就能体会到SDF等技术为何不适用大气散射。在看大气散射之前先来看一下散射：
大气散射中的散射：
大气散射中的散射根据实际情况可人为分成两类，一类是向外散射，一类是向内散射。后文称为外散和内散。
img
img
这里提一点：内散射允许我们观察本来不能观察的光线，而实际的物理效果就是光晕现象。因此实际散射能接触的光线是下图这样的：
img
实际散射过程：
让我先掏出这一张图片，这是我们切入大气散射的钥匙，我们将通过这一张图片深入大气散射，并且掌握大气散射。
img
我需要声明地是相机位于A点处，现在在p点发生了一次大气单散射，我需要得到p的情况。让我们简化这个过程：
1.相机的方向从A到B
2.尊重物理,我们将要考虑每个点的散射真实情况，也就是上文所说的外散和内散的结合
3.光源只有太阳
4.光的p量会在穿过大气层的时候受到散射的影响也就是cp的阶段
5.当光在p点时候又会受到内散射向相机进行传输
6.并不是全部光到了A，在AP的过程中也存在散射偏转
这是对于这一个散射过程的语言描述，接下来让我们进入一个有趣的旅程，用数学来描述这个过程
数学描述
img
在物理世界当中，我们会发现太阳的传输大致有两个阶段，一个是宇宙阶段，一个是大气阶段，宇宙阶段我们假设是真空条件，并不会减少光的传输量，那么现在我们有一个光照的最初值Light_InT_C;
很明显在c到p的过程中是会发生光照衰减的，会减少光照量。定义一个透射率：cp两点的光照量比值
img
我们知道这个cp的过程也是存在散射的，让我们先过滤这个事实，让问题简化（因为在物理上C到P的过程和P到A的过程是一样的），先来到p点的光照量。
后文会对这一点进行说明。
有了光照比值也就是透射率，那么就能得到p点的光照量：
img
散射函数：
回忆一下我们的语言描述，得到了p点的光照的时候，是不是有通过内散得到转向A的光照。所以我们这里要先定义散射函数。这也是十分头疼的一部，不过笔者会提供推导过程。真实世界中为什么天空是蓝色的，夕阳的时候是红色的，我们都知道是散射，这个散射有具体的专业名词来描述这个现象，瑞利散射。
img
在其他关于大气散射的文章当中你会看到这样的一个式子：
img
实际上这就是计算的是p点的光照在p点进行了内散然后发生了衰减到达了A点的过程。而具体的散射函数S（xxx）这一堆是什么东西呢？笔者会进行推导，进行到这里其实对于大气散射就应该有一个大致的概念了，但是存在一个问题，这个过程并不是正确的。
问题：
在PA的过程中，不是一个P点，在PA段中也有无数的P点进行着同样的过程，而A点得到的结果是无数点发生大气散射后的积分求和。这个过程我们称之为：数值积分
相信你，读到这里的时候，依然是理解笔者在说什么的，你也和我想到了同一个问题：光照呢，每一个点的光照呢？
光照假定都是一样的，因为实际物理世界的太阳光可以用平行光进行近似，为了计算简单，就不考虑太阳光的偏转的，都设定为一样的光照方向。
img
那么这个时候我们的A点的光照实际是每一个点的求和
img
或许你会有一个小问题，为何有ds。这实际上是数学上的一种近似处理方式，在计算每一个点的时候，他们各自的光照贡献度是不一样的，通过乘以一个ds达到一个连续现象的近似结果。如果不能理解的话，暂且记住：ds是为了让结果更加准确
然后让我们回忆起前文的公式
img
现在来看这个公式不是那么可怕了，我们至少知道它在干什么。那么结合上面所说的，一个大气散射完整的过程A点的光照贡献总和就应当是下文的这种形式
img
其实到这里，你已经对于大气散射了解了一个完整的简单过程。接下来，让我们更加深入，来到每一步的具体实现部分，那么首先便是散射函数
散射函数定义和推导：
散射是有两种方式的，按照结果为导向就是前文所说。但很明显，我并不是在说这个，而是散射函数的不同，目前大多数的光学效应有两种：瑞利散射和米氏散射。后文会进行推导说明两种散射函数。
瑞利散射的定义：
一束光在空旷的空间中，撞击到了一个粒子，那么光发生的运动就能通过瑞利散射来预测。
img
这个过程说明了在瑞利散射中常见的现象：基本上不会发生偏转了90度的光而更多的光照偏转是由下图的方向：
img
那么什么光会倾向于这样的主方向偏转呢，仔细想一下，光中有七种不同波长的光，不同的波长决定了不同的偏转倾向：
img
img
这里大致就能说明为啥看到的天空是蓝色的。让我们回到整体瑞利散射函数上：
瑞利散射原始函数，笔者搜寻众多，没有找到相关的推导。感觉是一个经验公式。
没错，前文所使用的散射函数就是瑞利散射函数。看着似乎有点复杂，整不明白这是什么玩意，说明一下具体的字面意义：
其中，λ是光的波长，n是粒子折射率，N是海平面处的大气密度，ρ(h)是高度h*处的相对大气密度（即相对于海平面的密度，可以理解成h处真正的大气密度与海平面处大气密度的比值，因此它在海平面处值为1，随着h增加不断减小）
推导瑞利散射：
瑞利散射函数分为了两部分，瑞利散射系数和瑞利散射相位函数，也就是这样：
img
前置推导：密度系数
这一步的推导我想放在前面，因为涉及到了瑞利散射的推导。所以放在前面，会在推导瑞利散射函数的更加觉得亲切。
这个里面有一个东西是R（h）/N的东西。虽然有文字说明，但想必也觉得十分痛苦，这是大气的密度比。
从逻辑的角度来看，散射强度与大气密度成正比是对的，每平方米更多的份子意味着光子被散射的机会更多。挑战在于大气的组成十分复杂，由于具有不同的压力，密度和温度的几层组成，幸运的是大部分的瑞利散射发生在大气层的前60公里，也就是对流层，在对流层中，温度是线性下降的，压力则是指数级下降。
img
而上面的R（h）来源于这里：
img
但是或许你也意识到了，这个变化过程不是如此线性的，而是指数级的。遵循指数衰变，实际上的R（h）是这样的
img
而其中的H0是比例高度的比例因子。对于底层大气的瑞利散射，密度比的图是下图这样的：
img
指数衰减简要说明：
在进行一次的内散之后，这个光还会进行内散吗？是的，它会。令人感到复杂的一个过程，但是不要急，让我们完成这个过程，这个过程中也会发生散射，意思就是发生衰减，那么这个衰减的系数，设为 β
那么衰减一次：
img
img
但实际过程呢？是发生无数次的衰减：
img
幸运地是，这个东西是指数函数的定义：
img
每一个点在P到A之后传输的过程也被我们成功的描述出来了。
这里再回忆前文，C到P的过程和P到A的过程是不是一样的。那么同样的描述是不是也能用于C到P。那么我们前文用的简陋不堪的透射率就能不使用了，而是使用这个指数函数衰减来描述C点光照到达P点后的光照量
瑞利散射系数推导：
请注意到这样的一个事实：瑞利散射是在三维空间中进行的：
img
img
img
img
img
瑞利相位函数推导：
瑞利散射相位函数描述了角度的依赖性。
我们得到了最终使用的瑞利散射的系数的推导过程。那么为了得到完整的瑞利散射，接下来要推导瑞利散射相位函数：
img
如果到这里你依然没有记住也没有理解我在说什么，那么请先记住下面的公式：
瑞利散射函数：
img
瑞利散射系数函数：
img
瑞利散射相位函数：
img
密度比：
img
其中H=8500m时候称之为刻高高度
海平面的瑞利散射系数：
img
完善大气散射：
在这里其实我们已经基本说完了大气散射，但是，介于为了让效果更加真实，我们引入了无数的点在P到A的描述和平行光的描述。无数的点的描述，通过指数衰减已经可以得到一个衰减过程，不需要再去考虑。但是平行光的过程呢：
很明显这个衰减也不完全是指数的，也存在大气密度的影响，为了更加真实地描述大气散射，通过把PC阶段分成两个阶段（举例说明，其实不止两个），先计算Q点的衰减，再计算P点的衰减
img
img
img
如果是均匀长度切分，那么就能如上面的数值积分方法一样后面乘上一个ds达到近似地效果。
img
再简化一下写法：
img
那么我们就得到了平行光的描述，大气散射函数也到此为止。
前文简陋说明的替换：
从c到p的透射率替换
img
瑞利散射中实际使用的衰减因子的替换：
img
由于许多的衰减因子是常数，所以可以分解出来，就是下面的样子
img
用于求和的量称之为光学深度D（cp）。这是shader里面要实际计算的东西，其余都是乘法系数，计算一次就行了。
附加内容：米氏散射定义和推导
前文的叙述已经是一个较为完善的大气渲染思路，为了保证读者的思考连续性所以没有说明米氏散射，在进入shader编程之前，让我们停下脚步休息一下，来体会一下米氏散射的力。米氏散射不仅可以用于大气的渲染，也能做地表的雾气。一个真实且美丽的大气系统，是离不开瑞利散射和米氏散射的共同作用的。
米氏散射和瑞利散射的不同之处：
作为一个知识补充了解，这边还会说到一个散射：拉曼散射。但是不会进行讨论和推导，仅仅是作为对比使用。
瑞利散射：
粒子半径小于十分之一光波长，比如空气的分子与可见光。
瑞利散射光的强度和入射光波长的四次方呈现反比。因此紫光比红光更容易发生瑞利散射。但是天空为什么不是紫色的，其实是的，但原因在于眼睛对于黄色和绿色的敏感度更高，往两边呈现钟形分布，因此看蓝紫色只是更感觉蓝色。上面也给出了瑞利散射的系列知识，因此不再讨论。
拉曼散射：
当光子打到直径大于自己的粒子的时候，会与其发生碰撞，导致方向发生偏转。但是其中多数的光子都是发生的弹性碰撞，故而散射的光子和碰撞之前的光子在波长，频率与能量相同（批注：此处说明的是光子）。那么自然会想到瑞利散射，拉曼散射和瑞利散射有何不同呢？在碰撞过程中，有一小部分的光子（大约千万分之一）和介质分子之间发生的是非弹性碰撞，存在能量交换，因而在散射后的波长，频率，能量都会发生变换，称之为拉曼散射。
img
两者在物理层面的区别：
分子的外层电子在辐射能的照射下,吸收能量使电子激发至基态中较高的振动能级,在10-12s左右跃回原能级并产生光辐射,这种发光现象称为瑞利散射.
分子的外层电子在辐射能的照射下,吸收能量使电子激发至基态中较高的振动能级,在10-12s左右跃回原能级附近的能级并产生光辐射,这种发光现象称为拉曼散射。
米氏散射：
粒子比光的波长相似或者大的时候，发生的散射现象。比如烟雾比可见光。
当微粒半径的大小接近或者大于入射光线的波长的时候，大部分的入射光线会沿着光前进的方向进行散射，这种现象被称为米氏散射。这种大的微粒包括灰尘，水滴，来自污染颗粒物质，烟雾等等。形成所谓的丁达尔效应。
因此，光的通路其实更是米氏散射的作用，因此实现体积光，也能采用米氏散射的方法进行制作，但本文并不会讨论体积光的米氏散射实现，笔者会在体积光的文章当中说明如何实现米氏散射，同时对比raymatching的方法。
米氏散射定义和推导：
米氏散射的历史来由：
在看了这么久的文章之后，让我们短暂休息一下。来到散射的历史。在1908年的时候，Gustav Mie通过对麦克斯韦方程组求解，得到了均匀介质球体对于弹性散射的严格解。而这个可以计算任何大小，任意材料的球形颗粒的散射，就是Mie（米氏）散射理论。
前置知识：
麦克斯韦方程组：
麦克斯韦方程组（Maxwell's equations)是英国物理学家詹姆斯 克拉克 麦克斯韦在19世纪建立的一套用于描述电场，磁场与电荷密度，电流密度之间关系的偏微分方程。总共由四个方程组成：
描述电荷如何产生电场的高斯定律
论述磁单极子不存在的高斯磁定律
描述电流和时变电场怎么产生磁场的麦克斯韦-安培定律
描述时变磁场如何产生电场的法拉第感应定律
麦克斯韦的方程推导并不是很难，因此本文也简要推导一下四大定律。让我们重新再出发，短暂来到物理的海洋。
高斯定律：
高斯定理（Gauss' law）也称为高斯通量理论（Gauss' flux theorem），或称作散度定理、高斯散度定理、高斯－奥斯特罗格拉德斯基公式、奥氏定理或高－奥公式（通常情况的高斯定理都是指该定理，也有其它同名定理）
在静电学中，表明在闭合曲面内的电荷之和与产生的电场在该闭合曲面上的电通量积分之间的关系。 高斯定律（Gauss' law）表明在闭合曲面内的电荷分布与产生的电场之间的关系。高斯定律在静电场情况下类比于应用在磁场学的安培定律，而二者都被集中在麦克斯韦方程组中。因为数学上的相似性，高斯定律也可以应用于其它由平方反比律决定的物理量，例如引力或者辐照度（来自百度百科）
img
观察这一副图，高斯公式实际就是在描述这样一个过程，输入等于输出。想象一个水流从A到B，这个过程当中，输入的水流始终等于输出的水流。本质就是守恒。用作数学语言表达就是下面这个方程：
img
其中三个cos部分，为外部法线向量的方向余弦
也就是说矢量通过任意闭合曲面的通量等于矢量的散度对闭合面所包围的体积的积分。这个公式乍一看十分夸张，有体积分，有偏导，还有曲面积分。所以慢慢拆分研究。
高斯公式推导理解：
等号的左边是一个体积分（也就是微积分中的三重积分，类比面积分来说，就是拿了xyz三个变量做了三重积分实现了计算体积质量的计算）。体积分的标准公式如下：
img
介于笔者的思想：让读者从零完整理解知识体系，不会默认读者明白什么理论，因而会短暂说一下体积分。这一部分有相关知识的朋友都能跳过。
用定积分的思想：可以延长三个方向把几何体划分成许多小块（i,j,k)，其中每一个小块的长度是
ΔXi，ΔYj，ΔZk，其中Xi，Yj，Zk分别是ΔXi，ΔYj，ΔZk中的任意一点，物体的总质量可以用极限表示为：
img
用三重积分可以把这个三重极限表示为：
img
让我们回到等式的左边：
img
所以这个方程实际表示的是从左边进入的所有的流体的量，其中具体的表示是，通过P对x求得微分，Q对y求得微分，R对z求得微分，得到一个微分的立方体。这个说法并不准确，因为实际上只是xyz三个方向的变化率而已，但是请允许我这么说，让这个公式更加容易理解，通过上面关于体积分的表示，可以知道现在得到了一个微分立方体，这个立方体其实是xyz三方向微分总和，那么通过一个三重积分就能得到整个体积的质量（在这里就是流体的量），后面的dV实际上就是dx,dy,dz。
高斯公式左边推导：
img
先是得到在水平方向的变化率：
img
这个时候变化率乘以长度，就能得到速度：
img
再乘以时间就能得到位移，所以在水平方向的位移表示为：
img
这是移动的距离，那么加上原来的距离，就能得到在水平方向下的微元团：
img
那么同理可得，在y方向和z方向的微元团也是同样的思路：
img
img
所以在流动一个微分下的体积变为了：
img
这是原来的体积加上新的微元体积，那么要减去原来的体积买得到相对的体积变化量：
img
那么再进行上面的逆过程，得到相对体积的变化率就应该是：
img
化简得到：
img
因此再进行对于dV的三重积分，就能得到进入的流量。
也就是：
img
高斯公式右边推导：
img
所以根据题图可以得出目前三维条件下，流量的计算公式为：
img
积分后，便可得到通过边界
的总流量：
img
那么很自然地就可以得到高斯公式：
img
可喜可贺，在这里我们成功推出了麦克斯韦方程的第一个定律，高斯定律。高斯定律实际表示的就是一个守恒的关系。
高斯磁定律：
高斯的磁定律对于理解电磁波和导电介质之间的相互作用极为重要。但是这个定律可以被写作一个很简单的形式：
img
本文不会讨论过多单理论的物理含义，毕竟本义在于实现一个最为完整的大气散射系统
毕奥萨戈尔定律（The Biot-Savart law)
要推导高斯磁定律，那么首先就要了解磁场是如何产生的。本文中的讨论不会涉及微观尺度下的表现，仅仅是在宏观条件下，因此不会涉及到统计力学的相关知识，请读者放心食用。
电流的公式（高中知识）：
img
电磁感应强度公式：
img
电磁感应强度的微分形式：
img
其中r是与电流源（导线）的距离。我们现在通过考虑导线dl的每个元素的贡献来建立一个总磁场。为此，我们考虑到每一点的磁场方向，这是由导线的切向分量和位置矢量之间的叉积给出的。
而
img
则是r-r'指向的单位矢量
然后把这个方程两边进行在导线长度上的积分，可以得到磁场的表达式：
img
高斯磁定律推导：
有了对于磁场的了解之后，我们就能来看是如何得到高斯磁定律的。说明一下较为重要的矢量微积分特性，之前关于微分的比值其实可以称之为散度。那么在空间中的每一点的散度会发生什么呢：
img
然后对于这三个式子进行推导可以得出：
img
img
然后进行到这里的时候，我们距离高斯磁定律就不远了，但是还得说明一个特性：标量三重积。也就是如下的式子
img
三重积 ，又称 混合积 ，是三个向量相乘的结果
结合标量三重积，再进行推导，可以得到：
img
这里需要注意的是：
img
因为这是取了一个向量和自身的叉积。
同时对于以上式子的研究可以发现：
img
所以最终可以得出这两个式子是完全成立的：
img
也就能得到磁场散度是：
img
最后说一下高斯磁定律的微分形式：
img
法拉第定律：
幸运地是在后面两个定律，我们在高中的时候，基本都有接触，那么现在让我们推导完整个麦克斯韦方程组吧。
定义电动势：电路中单位电荷对于所受合力绕闭合电路一圈的路径积分：
img
那么这个时候E的路径积分是0.不过前提是E全是由静电荷产生的，这也是平方反比的直接结果。
根据高中所学的物理知识对于磁通量变化的方程有：
img
然后根据前面的电动势和磁通量的定义可以得到：
img
这个时候通过用stokes'Theorem转化可以得到：
img
安培定律：
回到之前提到的安培定律的微分形式，对该等式两边同时取散度, 数学立刻要求等式左侧为0 （divergence of curl must be 0)，等式右侧为电流的散度在恒定电流下为0 （之前给出证明），但在交变电流下却不是，比如一个电容器连接一个余玄交变电压，在两极板间建立一个封闭圆形曲线C ，高斯定律指出由于C包含的电流I为0，周围就没有磁场，但实验表明电容器周围存在磁场！这样安培定律就不适用于非恒定电流， 不过这也是意料之中，毕竟安培定律来自适用于恒定电流的Biot-Savart 定律。直到Maxwell 对之前提到的局部电流守恒公式换了个造型:
,如果在安培定律右侧减去这一项不就让散度恒为0了吗? 于是Maxwell 方程的第四个就这出来了,也就是：
img
麦克斯韦方程组总结：
img
看到这里，估计你也已经有了对于麦克斯韦方程组的基本认识，推导过程其实相对精简，笔者不想在物理这里停留太多的时间，所以只是大致提供了一定的思路，那么接下来让我们进入重头戏：推导米氏散射。
米氏散射的正式推导：
请回忆起来，米氏散射是对于与自身波长差不多的粒子的碰撞发生的散射现象，形成了包括丁达尔现象在内的许多有趣的光学现象。而米氏理论为我们提出了一组关于散射电磁场的基本方程，可以很容易地在计算机上实现、并且通过计算机模拟得到光是如何从这个粒子（后文称之为小球体）上散射的，以及散射如何取决于光的颜色的，也就是光的频率。
为了更高效地推导米氏散射，会对于一些方程进行约定。开始回忆米氏散射的问题：在空气（或者真空）中的电磁辐射平面波撞击半径为a的小介点球体，例如水分子。为了找到这个散射结果，首先就得得到一个坐标系来描述散射情况：
坐标系的确定：
img
在这个坐标系的描述中，原点位于小球体的中心，其中补充一些基本知识：在笛卡尔的坐标系中的单位基向量表示为ex，ey和ez；球坐标中的单位基向量表示为
img
。为了得到这个散射的基本情况，需要找到麦克斯韦方程组的解，通过引入一些矢量的球谐函数，满足赫姆霍兹方程的矢量场，可以更容易找到描述辐射的解、他们是球谐函数的类似。
球谐函数推导和说明：
球谐函数（SPherical Harmonics)，是一组基函数，所谓的基，便要结合上文关于坐标系的确立。对于基函数，其实我们都是很了解的，只是以前不会这么称呼，列如多项式基函数：
img
再比如三角函数基函数：
img
基函数中的基到底代表什么呢？在向量的时候，我们都会说基底，意思是通过得到一个xyz的基底，任意的向量都能分解到基底上，意思就是通过基底表达出来。那么基函数的作用呢？自然也同理，基函数确立之后，任何一个函数都能描述成几个基函数的加权和。
例如：
img
那么球体的基函数是什么呢？
其中 P_{l}^{m} 是伴随勒让德多项式，而后面的 K_{l}^{m} 则是一个缩放系数，用于实现归一化。
伴随勒让德多项式： 勒让德方程式是物理和工程领域里面经常遇见的一类常微分方程，同时也是拉普拉斯微分方程的一种变形，当试图在求坐标当中得到三维拉普拉斯方程（其余偏微分方程也一样）的解的时候，问题就会编程对勒让德方程进行求解，而此时方程满足 |x|<1 时，可以得到有界解（解级数收敛）。并且当 nZ^{*} 时，X等于正负一处也有边界解。在这种情况下，方程的解随着n值变化而变化，构成的一组由正交多项式组成的多项式序列，称为勒让德多项式:
此时相当于就是在球坐标系下对于偏微分的方程进行求解，因此会诞生勒让德多项式。
产生的伴随勒让德多项式有两个参数构成 l,m ,并且定义可以通过普通的勒让德多项式进行定义：
对麦克斯韦方程求解：
在推导过程中，我们使用麦克斯韦方程的宏观版本：
img
其中关系为：
D=∈E\\B=uH
由于前文假设的球体，是建立在水分子的基础上进行散射，那么就可以假设: p=0,J=0
对于上文的问题再进行细分说明就应该是：
形式为:
的入射平面波撞击了一个水滴，然后开始求解散射问题。因为在米氏散射中粒子可以大于投射的波长的，因此通过这样的说法，帮助读者减小想象力上的困难，同时帮助理解阳光照射到大气中的雨滴发生的一系列散射自然现象。
我们通过选择E和H作为变量，并且结合上文此时的:
E(x, y, z, t) = E(x, y, z)e iωt
H(x, y, z, t) = H(x, y, z)e iωt.
通过对材质u和内部，我们假设在半径为a的球体（理解成液滴）内部取得2 和 2，那么对应的在球体外部取值就应该是1 和 1.故而要求解的方程就变成了：
img
结合以上方程可以发现E和H必须要满足向量赫姆霍兹方程：
img
其中式子的说明： k 2 = ω2 △ = ·
球坐标中的微分算子：
我们将要讨论球坐标中的微分算子，并且说明如何使用变量分离的方式找到赫姆霍兹方程的解，同时为了解决后续的水滴散射问题，我们将在这一小节中了解球面贝塞尔函数和相关勒让德函数的深入讨论。
如果f是一个可微的函数，那么f的梯度就可以在球坐标系中通过下面的式子计算得出：
img
那么一个可微的向量场的散度x就能通过下面的式子得出：
img
同时由由于可微向量场X的旋度是向量场，其分量是球体坐标中下面式子可以计算得出的：
img
然后这个时候，神奇的事情就会发生，可以得到平滑函数f的拉普拉斯算子，如下：
img
为了求解赫姆霍兹方程的平滑函数，我们将使用分离变量的方法，并且设置：
img
然后代入赫姆霍兹的球体坐标方程当中就能找到解：
img
如果我们乘以整个式子以
img
，就会发现包含
img
在内，取决于
img
而没有其他的变量，也正是因此，必须有一个常数-m2，我们将其写为：
img
那么现在公式就会变为：
img
同理可得，可以得到其他的常数c，例如：
img
img
那么现在φ方程就十分容易求解，可以取
img
实际上，只有当分离常数的形式的时候，才会存在周期解，对于常数c，需要做更多的事情，对于某些非负整数来说，常数C需要采取 l(l+1) 的形式才行，并且m可以取到负一到一之间范围的值。
对于径向方程，通过输入x=kr可以得到其方程：
img
而这个解就是球面的贝塞尔函数 ω_{l}(x) ,存在两种形式：正则函数，记作 j_{i}(x) 和不规则函数，记作 y_{1}(x)
然后对于角度方程，通过输入
img
可以发现满足方程，进而得到：
img
而这个函数就是常规下的伴随勒让德函数plm。
任何在原点处最多不规则的赫姆霍兹方程的解f都可以用以上的函数进行展开。由于特殊函数：伴随勒让德函数和球面贝塞尔函数在后续部分还有作用，所以我们将在后文进行探讨。
另外我们经常会在一些较为紧凑的空间上遇到正交的函数集，也就是关于紧凑空间上的积分问题，由于所考虑的空间是十分紧凑的，我们只能找到可数的正交函数集，而在这种情况下，克罗内克函数
img
的使用就无处不在了，如果
img
或者
img
，
img
同时
img
，通过使用克罗内克的符号就能得到一组正交函数的特征：
img
如果此时Na=1，那么该集合就是正交的集合。
矢量球谐函数：
没错，我们又回到了球谐函数上，前文对于球谐函数的说明，只是大致的说法和了解。我们首先考虑在满足赫姆霍兹方程的基础上，去构造一个向量赫姆霍兹方程。基本的思想就是通过将两种向量场与f进行关联，用L和K表示，其中r表示为径向矢量场
img
然后定义：
img
通过运算，很容易就能得出：
img
，其中的恒等式正式麦克斯韦方程组为E和H规定的恒等式。一切又绕了回来了，也在这里，我们阐明了为何麦克斯韦方程组会用于解决米氏理论中的散射问题。
根据前文所说，满足赫姆霍兹方程的复函数f都可以用以上的函数进行展开，其中plm是伴随勒让德函数，还有球面贝塞尔函数。
那么我们就能解得：
以及在这一步通过
img
与之结合就能发现：
img
特殊函数的性质：
这里我将说明一些函数的性质，球面贝塞尔函数和伴随勒让德函数的一些性质： 我们在上文定义勒让德多项式：
img
然后又知道勒让德多项式是微分方程的正则解：
img
就能得到勒让德多项式是正交的：
img
在其他关于勒让德多项式的文章当中，也能发现这一性质，但是我们还能找到一个佐证。那就是当勒让德多项式从P0到P1跨越多项式空间直到1次，那么就能让勒让德多项式从P0到Pn进行扩展，因此如果要保证在0到1的扩展就必须要有：
img
因此得到一个有趣的积分：
img
这一步也能通过前文的函数进行证明，笔者不做过多讨论了。
勒让德多项式还满足微分关系：
img
满足递推关系：
img
回想我们定义的伴随勒让德函数：
img
那么显然就能有：
img
故而可以证明伴随勒让德函数是微分方程的正则解：
img
使用这个微分方程，我们就可以通过部分积分，来证明对于固定的m存在：
img
通过以上的描述，plm，sin可以在球体上构成一组正交函数，那么球面上连续函数都可以根据这些基函数进行展开。通过正交性，这样的扩展就是独一无二的，如果fg分别是球面的两个函数，现在就能写出球面的内积fg：
img
球面贝塞尔函数wl的两种形式：
在上文以及说明了两种形式，现在给出定义：
img
img
球面贝塞尔函数就是微分方程的解：
img
另外一个相当有用的定义是：
img
，当找到了出射辐射的良好扩展的时候，就能得到：
img
另外再补充两个有用但是不需要提的等式：
img
平面波展开：
img
从一个众所周知的结果出发：
img
这个等式是成立的。
然后进行展开：
由于eikz在求解赫姆霍兹方程并且在原点处是正则的，所以可以在基本函数中展开为
img
同时没有φ的依赖性，因此我们有：
img
和之前一样，通过正价关系，可以得到：
img
推导主要部分：
根据矢量谐波扩展E：
img
我们选择wl=jl，因为E在原点处是绝对的正则，通过取内积：
img
得到了只有m=1的项式可以贡献的，排除了其余的项，因此继续推导：
img
为了求解以上的函数，注意到：
img
组成了两组正交函数，我们具有：
img
此处可以通过部分积分证明。
联系方程
img
就能得到一个展开式：
重新观察微分关系，可以得到：
img
从之前的推导式子中可以得到：
img
并可以进而推导得到：
img
因此存在关系：
img
那么对于任意系数c1都可以进行重写了：
img
也就是说现在c1=
img
但是在过程中我们还必须要解决：
img
对于较为低的p和k的值，可以知道以下关系：
img
代入上面的式子可以发现必须要解决：
img
通过比较不同球面贝塞尔函数的系数，可以得出要求解的方程的解是：
img
那么现在把入射辐射的电场和磁场写为：
img
检查部分：
检查径向方程：
img
通过使用jl（x）进行检查，获得方程：
img
img
结合前文所说，再进行推导可得：
img
用 π_{l}^{+} 进行检查：
通过前文式子得出：
img
使用类似的代数步骤可以发现：
img
也就能称c1的式子为：
img
同时也能表示 π_{l}^{+} 为：
img
使用和上面一样的方法就可以得到：
img
有了这些，就能推导出：
img
通过计算重新排列得到：
img
另一个方面， π_{l}^{+} 有：
img
img
再比较不同球面前的贝塞尔函数前的系数为
img
img
实现接口条件：
对于散射电磁场，还需要知识处理表面电荷和电流的水滴的界面条件，为了处理这些表面上的条件设为：
img
所以可以重写矢量谐波了：
img
同时，整数M只取非负数，有了这个限制，可以得到：
img
使用这些关系联系上文的正交性：
img
那么在球体内部上内积的正交关系就应该是：
img
而其他部分的内积就是0：
img
另外由于我们首先关注的是连续性，那么还得考虑向量积才行：
img
同时也由于正交特效，可以推导得出：
img
最后推导：
对于内部的辐射进行一般化分析：
img
结合前文：
img
然后得到内积，并且采用m=1的时候取适当的方程：
img
化简得到：
img
然后为了得到更好看的形式，引入折射率：
img
同时定义：
img
那么就能化简得到：
img
我们分散场和液滴内部场的扩展到现在终于变为了：
img
而这个便是我们最终的表达式
原理到实现的过渡（附加内容）：非预计算实现大气散射
本段说明：
本段属于附加内容，不会作为对于读者的要求，可以把这一部分当做“小说”来看。因为自从在08年之后，基本不再使用这一种方法，但是在最早期（上个世纪九十年代）的时候，本文所说的方法却是十分重要，也为后续的预计算提供了一定基础。本段的引入在于完善一个历史的理论发展进程。读者需要注意的是：本段仅仅是原理部分，后续代码会在实现篇给出。
求解散射方程：
从上文可以发现，在实际的大气散射的计算机计算当中，实则在进行一个无法解析求解的嵌套积分。虽然可以用梯形法则使用数值积分来做，但也会让积分的计算成本更高。因此为了尽量减小计算的压力，实则也是假设相机始终在地面上或非常靠近地面。因为这可以让大气在所有高度都具有恒定的密度。
回顾一下前文：笔者曾说过一句话，“大气散射实则是对内散射和外散射的贡献积分问题”。
渲染模型：
img
很明显，我们需要地是两个大步骤：第一确定积分区域。第二，开始积分。首先说前者，射线从A到B的过程中，是从进入大气层到打到地面的过程。需要做一次射线检测来确定A和B两点的位置。而当我们有了AB的长度之后，想要近似地去描述沿其长度的大气散射的积分就是我们的第二个目标。命名其中的五个点作为研究对象。例如P5这个点，阳光到达p5的时候会发生散射，同时在传输到相机的过程中也会有衰减，这都是前文说过的事情。
散射过程：
上文我们推过散射方程了，读者应该对此深有记忆了。但是这里使用的是Henyey-Greenstein函数的改编。该相位函数为：
img
这实际上是上文瑞利散射相位函数的变式而已，而且是通过改变数值的方式来近似模拟米氏散射相位函数。
瑞利散射的近似：把g设置为0
米氏散射的近似：把g设置为-0.75~-0.999之间。切记不要设置为1和-1，因为会让方程变成0，最终缺失瑞利或者米氏
外散射和内散射：
前文已经进行了相关说明，但是没有具体点明，现在正式点明
外散射方程：
img
仔细观察，你想到了什么？没错，就是前文所说的衰减函数。实际过程中，外散射就是内积分，也就是从A到B的衰减积分。
K是散射常数，取决于大气密度，变量h是采样点的高度。在我的实现中，高度是按比例缩放的，0 代表海平面，1 代表大气层的顶部。理论上，大气层没有固定的顶部，但出于实际目的，我们必须选择一些高度来渲染天穹。H0是标高，这是找到大气平均密度的高度。
积分的部分决定了光学深度，具体的计算是从A到B的大气密度乘以射线的长度。可以视为基于光线沿光线路径中的空气粒子数的加权因子。方程分为两部分，一部分是4πXk（）。这个说明的是光线散射出的光量。而后面就是前文说过的指数衰减积分（说法是指数衰减积分，但是在实现中的做法是做一个密度比得到衰减的因子所以也能说是大气密度），把每一个点的衰减加起来乘以长度得到总衰减再乘以光量。那么梳理到这里相对简单了，但是具体的做法还需要声明：米氏散射不依赖于波长
内散射方程：
img
内散射方程描述了由于来自太阳的光散射，有多少光被添加到穿过大气层的光线中。对于沿从pa到Pb的射线的每个点P，PP c是从该点到太阳的射线，而PPa是从采样点到相机的射线。通过前文的描述，读者应该明白此处内散射方程在做什么了。不再进行过多说明
额外散射方程：表面散射方程
前面只是说：我们会散射太阳光，但是反射光呢？太阳照射到其他星球例如月球，然后反射到大气层的光照呢？我们如何进行散射，因此引入了表面散射方程
I&#39;v () = Iv () + I e() x exp (–t (Pa Pb , ))
表面散射方程中描述：Ie（）是反射的光量，会被后面的外散射因子衰减。如果全部衰减掉了呢？散射的光应该就是从太阳光来的光的散射量，也就是Iv。
实际上到这里，以前的人使用的方法的数学原理大致没了，实现会在后文给出步骤，欢迎移步观看。
原理到实现的过渡（正式内容）：预计算大气散射
本段说明：
纵然我们通过前面的原理推导，学会了大气散射的系统理论，涵盖了瑞利散射和米氏散射的计算过程。但很明显的是，这样的计算过程对于计算机来说也是庞大的过程。如何加速大气散射的过程，保证实时的帧率呢？在论文《precomputed Atmosphereic Scattering》中提出了解决方案——预计算大气散射。另外前提说明：本段的实现过程不包含臭氧层的渲染，后续会提及。
该方法的特色：
考虑了大气散射中最为重要和主要的散射：瑞利散射和米氏散射，同时实现了多种散射。
实现图片：
img
img
其中最为主要的是，基于光传输方程的公式，是实现了所有视点和视图方向和太阳方向进行预计算。并且所提出的方法可以在GPU运算，在几秒钟内得出预计算结果。这种预计算后的数据可以让我们在恒定时间内评估光传输方程，无需任何采样，同时考虑了阴影和光轴的地面。可谓是十分俱全。这个方法是第一个实现所有视点，所有视图，太阳方向，多重散射的方法。
后文描述所用符号：
img
渲染模型：
前面只是说了散射的过程和计算公式，但实际散射的空间呢？哪一些空间会发生散射，哪一些不会。如何决定这一件事情，考虑到渲染大气的影响要素：局部介质属性的物理模型（大气本身的空间），和对观察者眼睛的全局照明的散射的模拟（散射过程）这两个要素，这其中包含了与地面的交互，因为地面是不会进行大气渲染的 ，那么这个空间可以视为具有反射率a（与前文一致的定义说法）,法线n（表示高度场）的一个郎伯曲面。
郎伯曲面的定义：
用于反射的郎伯曲面，是从所有视角方向均呈现均匀明亮并反射整个反射光的曲面。其中伴随的是性质是郎伯反射率，郎伯反射率的提出是为了描述反射过程，进而对哑光的表面和漫反射表面进行描述。
例如各种表面的反射使用各自的BRDF（双向反射分布函数）的时候，如果是郎伯曲面，BRDF的值就是常数值。郎伯曲面在计算机图形学中算一种理想扩散表面。
CG中的物理模型说明：
CG中最常使用的物理模型就是基于空气分子和气溶胶颗粒两种成分的模型，大气的密度在6330千米到6420千米这个空气密度之间递减的薄球形层，这也是我们将要采用大气模型的特质。
例如下面这张图：
img
需要说明的是，其中a和b是指的是单散射和多散射
这个穹顶到X0的位置，按照道理来说，就是我们需要进行积分的区域来计算散射衰减。但是为了减少计算量，我们采取的近似方法是忽视遮挡（shadow volume)的区域积分，如下图的处理：
img
这个时候，我们只积分X到Xs的区域，那么减少很多计算，包括积分界限的区分，都能做到明确的处理。当然会损失一些效果。但影响很小。
回忆前文内容：
瑞利散射取决于空气分子的密度和入射光的波长，而米氏散射系数仅仅取决于气溶胶的密度，而不取决于波长。为了帮助读者回忆前文的内容，这里再拖沓地拿出公式。
瑞利散射系数定义：
img
再顺带说一下其中的物理量：n是空气的折射率，Ne是海平面空气的分子密度，然后 λ 是入射光的波长。如果有读者看到这里说：你写错了！前文推导的结果和这个函数不一致，那么请仔细观察上面的函数，其实两者是等价的。
米氏散射系数定义：
笔者在前文完整地推导了米氏散射系数，但正如所见，其复杂性决定了无法用于目前的计算机实时渲染上，但是我们有更便宜的做法，公式如下：
米氏散射相位函数：
上文纵然说了一下HG函数（Henyey-Greenstein)函数，但是只是应用了一下，这里简单说明一下由来。读者可以看到笔者推导得到的米氏散射相位函数的解十分复杂，难以用到计算机当中。最常用的近似方法就是HG函数，HG函数的首次提出在1941年，形式为：
img
稍微细心一点的朋友会发现这和上文说的HG函数不一样，上面的HG函数是本函数的改编，而这里笔者将采用的是上面给出的形式，该改编在1992年得出，参数g通常称之为首选散射方向，也有一些晦涩难懂的论文喜欢叫均值余弦，为什么通过改变g值可以近似瑞利散射的相位函数呢？这是因为g值描述的是米氏散射的各向异性，换一句话说，g值控制散射是、几何的不对称性，正值会增加向前的散射光量，负值会表示向后的散射光量。
img
瑞利散射相位函数：
img
与HG函数相比，瑞利相位函数就简单很多了。
密度函数：
使用前两个函数，已经可以计算海平面上单一点的散射。该密度函数指的是空气中分子和气溶胶密度的变化，在前文实则也已经有了说明，这里不做过多讲述。
img
为何是指数的形式，前文也已经有了十分详尽的推导过程。这里也不做过多说明。
渲染解方程：
散射过程：
img
通过上文的学习，可以得到Is的函数方程为：
img
实际要解决的也是这样的函数方程，其中 各个函数是什么就不再单独说明了。
衰减函数：
上面方程可以得到定义在某一个方向上散射的光量，但是我们还需要表示这种光是如何衰减的。也就是指数衰减的部分，相信读者这里有点迷惑了，我在说密度函数的时候也在说指数衰减，笔者很有必要地说的是：前文说的密度函数的部分形式是指数衰减的形式，但是实际并不是在说如何描述光的衰减的。还记得我们最开始有一个透射率函数吗？我们也成功推导出来了，那才是我们的衰减函数，函数方程为：
img
相信读者不理解Be是什么东西，请查看笔者最开始给出的系数表，其中有说明Be的意义，Be是消光系数。当分子中只是存在散射现象的时候， β e R = βR ,但是作为气溶胶的颗粒是存在散射现象和吸收光线的现象的，此时的消光系数 β e M = βM + β a M 。实际并不是很难理解，所以光在传输过程中的实际情况应该是：
img
前文实则对这一部分说了不少了，实则到目前为止我们已经可以去做出一个逼真的天空了，但是笔者不想在此处停留下来。在衰减的过程中，我们似乎只有大气和地球，我们忘记了什么。在真实的物理世界当中是存在臭氧层的，那么臭氧层又如何通过对衰减函数的操作来控制呢？
臭氧层的特点在于不会散射光，只是会吸收光，因此只需要在衰减方程（透射率方程）中考虑臭氧就可以了。
img
其中O是臭氧的消光系数，臭氧的密度函数为： ρO = 6e7 ρR ，臭氧的消光系数可以在物理世界中得出，请参照下表：
img
我们完成了一个有臭氧层的大气散射。而这是08年论文中所没有完成的事情。
单散射：
和前文不一样的是，这里不再会说外散射函数和内散射函数，因为比起上面的非预计算的方法，我们这里更为复杂，同时最重要的是在单散射方程中就蕴含了外散射和内散射，因此上文对于外散射和内散射的说法只是为了让大家更加轻松地去实现大气散射。让我们开始单散射部分，单散射的定义在上文也给出了，但是这里再重申一次：在到达观察者之前经历了一次且仅是一次散射事件的单次散射光，结合上面的散射过程函数和透射率函数可以得到单散射公式：
img
img
然后把米氏散射和瑞利散射的结果相加就能得到散射结果：
img
多重散射：
令人激动的标题，在08年后才有了一个完整的多重散射的方法，虽然可以通过嵌套多维积分获得高阶散射阶数，实现多重散射，但是这太慢了，我们会使用一个新的函数来帮助我们实现多重散射：
img
该函数定义了到达某一个点p的散射k阶的光量和在方向v上的散射。例如一阶的聚焦函数是将所有方向上的单散射进行积分并使用相位函数。k阶的聚焦函数表示聚焦在p点的光，该点的光是恰好经历了k个散射事件
img
使用k阶函数让我们对上面的单散射函数的最终公式变形得到我们的多重散射最终公式：
img
img
然后总强度就是把瑞利和米氏的结果相加，最后再计算k阶散射后到达观察者的所有光，需要做的是就是把单次散射和更高阶散射的结果相加：
img
预计算
在进行这一段的冒险过程中，我会在后文给出一段实现无需用表来实现大气散射。可以来体会一下两者的区别。
在这一段中，笔者会介绍预计算的原理同时在最后给出优化方法。
前述
为每一个观察者位置P（XYZ），每一个视图方向（XYZ），每一个光源方向（XYZ）提供散射信息就需要一个具有9个维度的查找表，这不仅会让纹理变得十分巨大，而且会让散射变得十分漫长。所以为了简化问题，采取以下假设：
1.地球可以被看成一个完美的球形
2.空气中的颗粒物密度只随高度变化和地理位置无关
3.由于地球与太阳之间的距离很远，所有到达大气层的光都可以视为平行光
通过以上的假设，问题可以简化维度到四个自由度：高度，日天顶角，视天顶角和日视方位角。然后进一步简化，可以去除太阳视角的方位角，就只需要三个维度。所以我们最终查找表的三个参数为：高度，视角-天空顶角，太阳-天空顶角。
那么瑞利散射的相位函数就会变为：
img
另外需要对瑞利散射相位函数做一个处理，由于去除了太阳视角方位角（影响地球对大气的阴影）纵然这种现象由于多次散射之后几乎看不见了，并且时常会被地形遮挡。因此需要去除太阳视角方位角，也正因此需要调节一下瑞利散射相位函数，不然的话，最终会在太阳的米氏散射上有块状的阴影，解决办法就是把相位函数的评估推迟到片元shader上。
img
纹理化：
为了能够在3d纹理中获取到值，需要把三个参数转化到【0,1】的范围当中
img
这种转化函数最直接的实现就是使用线性参数化：
只要纹理的分辨率足够就可以存下，但后续有人发现可以通过线性参数化方程的改变让精度更一步提高，之所以需要提高的原因在于：在地平线附近的视角处的色调大多不一样，因此参数化必须在这些区域提供更高的精度，但是上面的方法是均匀的，无法做到，后来由Bruneton和Neyret两人提出了线性参数化的重映射方法。
后文部分用 cs代替cos（s）,cv代替cos（v）方便说明
线性参数化的升级：
高度重映射：
由于大气密度随着海拔高度呈现指数下降，因此天空颜色的差异在低海拔的地区最大，这意味着在低海拔地区应该有更多地参数化集中，重映射方法修改为：
img
img
视角-天空重映射：
和上文一样，应该在低海拔地区集中参数化。 重映射方程为：
img
img
太阳-天空重映射：
虽然太阳方向的重要性在白天不会发生剧烈的变化，但是我们不需要再表格中包含低于某个阈值的太阳角度，因为夜间几乎不存在来自太阳的非散射光
重映射方程为：
img
img
代码实现大气散射：
大气散射的实现（不用LUT表）
在前文说明实现的时候，关于预计算的部分，笔者曾说这是解决传统大气渲染模型的很棒的方法，但是以前的人关于大气散射怎么做的呢？他们是怎么做到不用LUT表，完成那样复杂的计算的呢？这其中又蕴含了多少前人的智慧。因此在本段，笔者会带领大家来实现一个不用LUT表的大气散射。
img
说明:
本段的内容参考来源GPUGEMS2的16章节关于大气散射部分。本身属于附加内容，因为作为单一技术来说价值已经不大，但是对于了解上个世纪九十年代的做法非常有利，去看看。介于本文到目前为止的复杂性，会在最开始加入推荐导读部分。
本文章的算法支持在DirectX shader Model2.0的GPU上实时运行，但是笔者打算用unity来做。原理部分参考上文，另外这一段的原理和实现其实算附加内容，读者可以带着轻松心情来阅读。
实现效果图：
img
img
img
相位函数实现：
根据上面的原理说明，瑞利散射相位函数和米氏的相位函数都是通过一个方程近似完成的。帮助读者回忆一下近似用的公式：
img
瑞利散射相位函数代码：
实现采用上文给出的公式，并且请读者回忆前文的说法：近似瑞利散射需要g取得0。所以代码就应该如下所示：
float Equation_ray(float cc) {
return (0.75) * (1.0 + cc);
}但这里的相位函数的散射强度不是很够，所以效果如下：
img
所以这里要把瑞利散射相位函数的取值变小，笔者采用的方法是GPUGEM2中的代码方法，最终效果表现就在上面。代码如下：
float Equation_ray(float cc) {
return (0.1875 / PI) * (1.0 + cc);
}米氏散射相位函数的代码：
Mie H-G相位函数是一个很简单近似米氏散射相位函数的方法，甚至比上面给出的方程更简单，H-G函数可以很好地体现mie散射的向前峰值的主要特质。但也有一个问题不能很好地正确模拟向后散射，为了克服这个问题，就能使用上面给出的方程双H-G函数来模拟米氏散射相位函数
//Mie H-G相位函数
// 1 - g^2
//F = ----------------------------------
// 4pi * ( 1 + g^2 - 2g * cos)^(3/2)
//这个函数非常简单易用， 他可以很好地体现mie散射的向前峰值的主要特征， 但是不能正确模拟向后散射，为了克服这个缺陷，就用了 双 H-G相位函数
// 改编自 Henyey-Greenstein 函数， 双 Henyey-Greenstein相位函数 的 单参数版 (原双相位函数拥有3个参数, 要确定3个参数非常复杂)
// g : ( -0.75, -0.999 )
// 3 * ( 1 - g^2 ) 1 + cos^2
// F = ----------------- * -------------------------------
// <4pi> * 2 * ( 2 + g^2 ) ( 1 + g^2 - 2 * g * cos )^(3/2)
// 4 * PI
#define MA_VERSION_0
float Equation_mie(float g, float cos, float cosSquare) {
float gg = g * g;
#if defined(MA_VERSION_2)
//H-G相位函数
float M = (1.0 - gg) / (4.0 * PI * ( 1 + gg - 2 * g * cos));
M = pow(M, 1.5);
return M;
#elif defined(MA_VERSION_1)
//双相位函数 外网大佬的做法，效果上相差不大
float FMA = (1.0 - gg) * (1.0 + cosSquare);
float FMB_L = 2.0 + gg;
float FMB_R = 1.0 + gg - 2.0 * g * cos;
FMB_R *= sqrt(FMB_R) * FMB_L;
return (3.0 / 8.0 / PI) * FMA / FMB_R;
#else
//双相位函数 标准实现
float MA = (1.0 - gg) * (1 + cosSquare);
float MB_L = 2.0 + gg;
float MB_R = 1.0 + gg - 2.0 * g * cos;
MB_R = pow(MB_R, 1.5) * MB_L;
return (3.0 / 8.0 / PI ) * (MA / MB_R);
#endif
}代码部分实在过于简单，没有太多好说的。读者请详看注释和方程进行对应理解
密度函数部分：
这密度函数部分是一个较为单独的过程，所以单独写出来了。代码如下：
float Density(float3 p, float ph) {
return exp(-max(length(p) - r_inner, 0.0) / ph);
}
//说是指数衰减，实则做的是光学路径上的密度积分的过程。当然是从不同的理解层面去说表面散射部分：
这部分在上面说了，遗忘了方程的话，往上翻阅即可。具体实现也较为简单，是不同内散射和外散射的额外部分。代码如下：
float OutScatter(float3 pSource, float3 pTarget, float ph) {
float sum = 0.0;
float3 step = (pTarget - pSource) / float(OUT_SCATTER_COUNT);
float3 pCurrent = pSource;
for (int i = 0; i < OUT_SCATTER_COUNT; i++) {
sum += Density(pCurrent, ph);
pCurrent += step;
}
sum *= length(step);
return sum;
}
//表面散射部分球体相交检测部分：
还记得最开始我们说大气渲染模型的时候吗，需要做一个相交的检测。具体实现逻辑如下：
float2 RaySphereIntersection(float3 rayOrigin, float3 rayDir, float3 sphereCenter, float sphereRadius) {
rayOrigin -= sphereCenter;
float a = dot(rayDir, rayDir);
float b = 2.0 * dot(rayOrigin, rayDir);
float c = dot(rayOrigin, rayOrigin) - (sphereRadius * sphereRadius);
float d = b * b - 4 * a * c;
if (d < 0){
return -1;
}else{
d = sqrt(d);
return float2(-b - d, -b + d) / (2 * a);
}
}
//球体相交检测，得出积分区域（大气层部分）外散射和内散射部分：
过程也相对简单。不再多说，上代码：
float3 InnerScatter(float3 eyePos, float3 viewDir, float2 e, float3 lightDir) {
float3 sum_ray = float3(0.0, 0.0, 0.0);
float3 sum_mie = float3(0.0, 0.0, 0.0);
float n_ray0 = 0.0;
float n_mie0 = 0.0;
float stepLen = (e.y - e.x) / float(IN_SCATTER_COUNT);
float3 step = viewDir * stepLen;
float3 currentPoint = eyePos + viewDir * (e.x);
for (int i = 0; i < IN_SCATTER_COUNT; i++, currentPoint += step) {
float currentStep_rayValue = Density(currentPoint, slope_ray) * stepLen;
float currentStep_mieValue = Density(currentPoint, slope_mie) * stepLen;
float2 pIntersection = RaySphereIntersection(currentPoint, lightDir, ObjectPosition, R);
float3 point_PlanetEdgeinLightDirfromCurrentPoint = currentPoint + lightDir * pIntersection.y; // 当前点 沿着光线方向 与 大气层最边缘 的交点
//内散射
n_ray0 += currentStep_rayValue;
n_mie0 += currentStep_mieValue;
//外散射
float n_ray1 = OutScatter(currentPoint, point_PlanetEdgeinLightDirfromCurrentPoint, slope_ray);
float n_mie1 = OutScatter(currentPoint, point_PlanetEdgeinLightDirfromCurrentPoint, slope_mie);
float3 attribute = exp(-(n_ray0 + n_ray1) * ray_ColorAdjust - ((n_mie0 + n_mie1) * mie_ColorAdjust) * mie_Adjust);
sum_ray += currentStep_rayValue * attribute;
sum_mie += currentStep_mieValue * attribute;
}
float cos = dot(-lightDir, viewDir); //相位函数描述了&#34;角度&#34;向相机方向散射多少光
float cosSquare = cos * cos;
float3 res = sum_ray * ray_ColorAdjust * Equation_ray(cosSquare) + sum_mie * mie_ColorAdjust * Equation_mie(-0.78, cos, cosSquare);
return res * 10;
}说是实现外散射和内散射，其实把之前的过程走了一遍。实现方面能说的不多，要深刻理解地话，请仔细参考前文的原理篇章和代码，都能看懂。
工程开源链接：
结束总结：
那么不使用预计算的方法所实现的大气散射，就到此结束了。其实相对后面两种，简单很多。所以笔者更愿意把这一段的原理讲解和工程实现说明当做给读者的休息时间。实现代码也不足两百行，但却实实在在地走了一次大气散射的过程。休息完毕，让我们开始本文的第二部分吧。
大气散射实现说明：
从笔者对于本文的安排来说，这一部分实际才是本文的第二部分：实现篇。前面的非预计算实现大气散射是作为对于读者关于90年代的做法的一个了解。
预计算实现大气散射：
在上面进行了关于大气散射的非预计算实现，那么接下来自然就是进行关于大气散射的预计算实现，而这一部分内容，则是需要有一定的编程基础和图形基础，笔者也有意识到本文到目前为止的整体长度过大了，所以会说一下核心代码方面的写法，至于原理请参照上文。
项目初始化：
依然是先从需求入手，既然要做预计算就得把Lut图的RT进行创建，同时在在ComputerShader中进行计算，然后把结果存成LUT图，用于skybox。但是为了方便后期的检查，以及运行的时候的测试来判断我们写的大气散射结果是否正确，先写一个控制光源的脚本，因为后期的大气散射计算也是基于相机数据的。定义一下该脚本的功能：可以控制光源的方向，可以控制相机的旋转和升降。
Scripts：Suncontrol
public float Height = 5;//起始相机高度
public Transform DirLightTransform;//光源的位置
public bool Help = true;//用于GUI，表明是否显示脚本功能
private Vector3 prevMousePos;//鼠标位置
然后自然就是在update函数当中对输入的数据进行判定，首先判定是否要更改光源或者相机的位置：
if (Input.GetMouseButtonDown(0) || Input.GetMouseButtonDown(1))//如果输入的是0代表鼠标左键，1代表右键
{
prevMousePos = Input.mousePosition;//更新位置
}
然后就是对单独的左右键的行为进行规定，例如左键，就是更改输入的DirLight的位置属性，右键同理，这部分比较简单就直接上代码了：
if (Input.GetMouseButton(0))//控制光源
{
DirLightTransform.Rotate(0, mouseDelta.x * 0.1f, 0, Space.World);
DirLightTransform.Rotate(mouseDelta.y * 0.1f, 0, 0, Space.Self);
}
if (Input.GetMouseButton(1))//控制相机
{
transform.Rotate(0, mouseDelta.x * 0.1f, 0, Space.World);
transform.Rotate(-mouseDelta.y * 0.1f, 0, 0, Space.Self);
}
在完成了左右键控制光源和相机之后，自然是希望相机更加自由，因此通过AD键控制其高度移动移动，而这部分本质就是更新上面的高度变量。
if (Input.GetKey(KeyCode.A))//控制移动
Height += 2000 * Time.deltaTime;
if (Input.GetKey(KeyCode.D))
Height -= 2000 * Time.deltaTime;
再然后就是GUI，用来表明该脚本的功能：
void OnGUI()//显示GUI
{
if (Help)
{
GUILayout.Label(&#34;功能&#34;);
GUILayout.BeginHorizontal();
GUILayout.Label(&#34;相机高度&#34;, GUILayout.ExpandWidth(false));
Height = GUILayout.HorizontalSlider(Height, 10, 40000, GUILayout.Width(400));
GUILayout.EndHorizontal();
GUILayout.Label(&#34;鼠标左键 - 旋转光源&#34;);
GUILayout.Label(&#34;鼠标右键 - 旋转相机&#34;);
GUILayout.Label(&#34;A/D - 控制相机高度&#34;);
}
}
然后就完成了对于相机和光源的控制：
这部分内容完成之后，便能正式进入大气散射的编写了。
预计算实现大气散射：
由于本文的shader编写是使用的ComputerShader，所以笔者会稍微介绍一些ComputerShader的基础。
ComputerShader基础：
// Each #kernel tells which function to compile; you can have many kernels
#pragma kernel CSMain
// Create a RenderTexture with enableRandomWrite flag and set it
// with cs.SetTexture
RWTexture2D<float4> Result;
[numthreads(8,8,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
// TODO: insert actual code here!
Result[id.xy] = float4(id.x & id.y, (id.x & 15)/15.0, (id.y & 15)/15.0, 0.0);
}这是最简单的一个ComputerShader代码段，众多教学博客都是基于这一段进行说明，本文也不例外。将会基于这段代码说明ComputerShader的作用机制
Kernel:
正如英文中所说的一样，“kernel tells which function to complie”，kernel会告诉编译器哪一些函数会被编译。其格式为：
pragma kernel 函数名
RWTexture2D:
R=read,W=write。所以RWTexture说明的一张可读写的纹理。可在脚本中对其进行读写的操作。其格式为：
RWTexture2D纹理名
numthreads:
numthreads实则是number of threads。线程的数量，所以在进行这一部分的说明之前，让我们了解一下什么是线程。在计算机当中，唯有多核心的cpu可以进行真正的多任务，例如播放音乐算一个任务，运行Steam也算一个任务。这一个一个任务对于计算机来说，就是进程（process），但是播放音乐这个任务当中，还能进行更多的事情，例如点击暂停，停止播放音乐等，这些事情则是线程（Thread）
在了解了什么是线程之后，就得来说一下怎么为每一个像素来指定线程了。首先说id，即是标识符而在ComputerShader当中的id则是分配的线程：
id :SV_DispatchThreadID指明了此时的id代表的是线程的id，而我们就要线程如何去计算像素，假如我们分配的id是（1，1），然后除以屏幕宽长，那么就会只会渲染最左上角的像素。所以我们要让表示每一个像素是干什么的，就得用整个xy去除以屏幕宽长，为每一个像素分配一个线程。这里可以参考笔者的另一篇文章《RayTracing:unity中实现光线追踪》：
CSMain:
这部分自然不用说了，就是对上文要编译的主函数的撰写。
对LUT的认识：
回忆一下前文实现大气散射的时候，要计算进行LUT图，而LUT图到底是什么，让我们来看一下，也只有明白Lut图是什么才能继续在Scripts中进行编写。
LUT英文的全称是Look Up Table,其简称就是Lut表，本质就是一张颜色查找表，是一种降低GPU运算量的技术，通过事先把颜色值存储在一张缓存表中当需要运算的时候，直接从这张表去拿索引对应的颜色值，本质就是拿存储空间换运算时间的算法应用，常用的领域其实图像调色的领域，而在计算图形的领域也有不小的应用，例如之前米哈游就公开的卡通渲染技术当中就用到了Lut技术。
Lut被分为1DLut和3DLut。本质区别在于索引输出值所需要的索引数，1D的自然就是用单个变量进行索引，而我们的3DLut自然就是用三个变量进行索引，大气散射计算的Lut自然也应该是3D的。
这部分的内容笔者思来想去，还是决定放在这个实现部分进行Lut理论的撰写，由于前文的理论太多了，再去补充，很容易就让读者失去阅读兴趣，关于Lut的原理部分，笔者是参考的论文：
《An Efficient 3D Color LUT Compression Algorithm Based on a Multi-Scale Anisotropic Diffusion Scheme》一文。下面开始讲述Lut的原理
Lut的原理：
如果令RGB为连续的颜色域[0,255]，那么就能用一个立方体来表示这些数据：
其映射函数，大致可以写为：
如果上面这张图不是很能明白笔者在说什么的话，可以参考下面这张图，3DLut本质就是一个立体的颜色隐射表：
然后经常在软件中看到的cube文件或者3d1文件到底是什么？本质就是三维的数据二维平面化展开：
所以3dLut的平面化其实是这样的一张图：
所以可以发现3dlut的数据处理就是一个三维rgb的映射，不过这个说法也不是很准确。但很明显的是如果颜色给与256个取值，那么lut的文件大小就应该是257 257 257的数据大小。那么如果拿到了根据三变量计算得到的RGB按照这样的对应关系进行3dlut的建立，就能得到其3dlut图，想拿回计算结果，就直接根据输入的变量去查找就行。
但是可以发现的一点在于：这样的数据存储消耗内存实在太大了，有没有办法去减少数据量，一般的处理就是降低采样来减少数据量。
用代码的形式来表述，就是：
def gen_3dlut(src, target, lut_dim, use_diff):
&#34;&#34;&#34;
:param src: n*3 array, float, 0-1
:param target: n*3 array, float, 0-1
:return: 3dlut
&#34;&#34;&#34;
lut = np.zeros([lut_dim*lut_dim*lut_dim, 3], dtype=np.float64)
lut_weight = np.zeros([lut_dim*lut_dim*lut_dim, 3], dtype=np.float64)
print(&#39;lut len:&#39;, len(lut))
if len(src.shape) != 2 or src.shape[1] != 3:
print(&#34;error input !&#34;)
return -1
#src = np.clip(src, 0, 0.999999)
# K near
bin_size = 1.000000001 / (lut_dim - 1)
src_ind000 = np.floor(src / bin_size).astype(np.int32)
src_ind100 = src_ind000 + np.array([1, 0, 0]).reshape(-1, 3)
src_ind010 = src_ind000 + np.array([0, 1, 0]).reshape(-1, 3)
src_ind001 = src_ind000 + np.array([0, 0, 1]).reshape(-1, 3)
src_ind110 = src_ind000 + np.array([1, 1, 0]).reshape(-1, 3)
src_ind101 = src_ind000 + np.array([1, 0, 1]).reshape(-1, 3)
src_ind011 = src_ind000 + np.array([0, 1, 1]).reshape(-1, 3)
src_ind111 = src_ind000 + np.array([1, 1, 1]).reshape(-1, 3)
print(src_ind000.max(), src_ind100.max(), src_ind010.max(), src_ind001.max(), src_ind101.max(), src_ind111.max(), lut_dim)
src_pos000 = src_ind000[:, 0] + src_ind000[:, 1] * lut_dim + src_ind000[:, 2] * lut_dim * lut_dim
src_pos100 = src_ind100[:, 0] + src_ind100[:, 1] * lut_dim + src_ind100[:, 2] * lut_dim * lut_dim
src_pos010 = src_ind010[:, 0] + src_ind010[:, 1] * lut_dim + src_ind010[:, 2] * lut_dim * lut_dim
src_pos001 = src_ind001[:, 0] + src_ind001[:, 1] * lut_dim + src_ind001[:, 2] * lut_dim * lut_dim
src_pos110 = src_ind110[:, 0] + src_ind110[:, 1] * lut_dim + src_ind110[:, 2] * lut_dim * lut_dim
src_pos101 = src_ind101[:, 0] + src_ind101[:, 1] * lut_dim + src_ind101[:, 2] * lut_dim * lut_dim
src_pos011 = src_ind011[:, 0] + src_ind011[:, 1] * lut_dim + src_ind011[:, 2] * lut_dim * lut_dim
src_pos111 = src_ind111[:, 0] + src_ind111[:, 1] * lut_dim + src_ind111[:, 2] * lut_dim * lut_dim
print(src_pos000.max(), src_pos001.max(), src_pos111.max(), src_pos101.max(), src_pos010.max(), src_pos011.max())
print(src_ind100.shape, src_ind100.dtype, )
print(lut[src_pos000].shape, (target-src).shape, lut[src_pos000][0], (target-src)[0],np.maximum(lut[src_pos000], (target-src)))
lut[src_pos000] = np.maximum(lut[src_pos000], (target-src))
lut[src_pos100] = np.maximum(lut[src_pos100], (target-src))
lut[src_pos010] = np.maximum(lut[src_pos010], (target-src))
lut[src_pos001] = np.maximum(lut[src_pos001], (target-src))
lut[src_pos110] = np.maximum(lut[src_pos110], (target-src))
lut[src_pos101] = np.maximum(lut[src_pos101], (target-src))
lut[src_pos011] = np.maximum(lut[src_pos011], (target-src))
lut[src_pos111] = np.maximum(lut[src_pos111], (target-src))
#
# src_weight000 = cal_invdis(src, src_ind000 * bin_size)
# src_weight100 = cal_invdis(src, src_ind100 * bin_size)
# src_weight010 = cal_invdis(src, src_ind010 * bin_size)
# src_weight001 = cal_invdis(src, src_ind001 * bin_size)
# src_weight110 = cal_invdis(src, src_ind110 * bin_size)
# src_weight101 = cal_invdis(src, src_ind101 * bin_size)
# src_weight011 = cal_invdis(src, src_ind011 * bin_size)
# src_weight111 = cal_invdis(src, src_ind111 * bin_size)
# src_weight = np.hstack((src_weight000.reshape(-1, 1),
# src_weight100.reshape(-1, 1),
# src_weight010.reshape(-1, 1),
# src_weight001.reshape(-1, 1),
# src_weight110.reshape(-1, 1),
# src_weight101.reshape(-1, 1),
# src_weight011.reshape(-1, 1),
# src_weight111.reshape(-1, 1)))
# src_weight = src_weight
# print(np.hstack((src - src_ind000 * bin_size, src_weight000[..., None])), src_weight000.shape, src_weight000.shape, src_weight000.dtype, src_weight000, )
# #use_diff = 0
# print(&#39;dd : &#39;, src.shape, src_ind000.reshape(-1, 3).shape, src_pos000.reshape(-1, 1).shape, src_weight.shape)
#
# np.set_printoptions(precision=3, suppress=True)
# # print(tt)
# tt1 = np.hstack((src, src_ind000.reshape(-1, 3), src_pos000.reshape(-1, 1), src_weight))
# np.savetxt(&#39;tt1.txt&#39;, tt1, fmt=&#39;%.5f&#39;, delimiter=&#39; &#39;)
# if use_diff:
# lut[src_pos000] += (target-src) * src_weight000[..., None]
# lut[src_pos100] += (target-src) * src_weight100[..., None]
# lut[src_pos010] += (target-src) * src_weight010[..., None]
# lut[src_pos001] += (target-src) * src_weight001[..., None]
# lut[src_pos110] += (target-src) * src_weight110[..., None]
# lut[src_pos101] += (target-src) * src_weight101[..., None]
# lut[src_pos011] += (target-src) * src_weight011[..., None]
# lut[src_pos111] += (target-src) * src_weight111[..., None]
#
# else:
# lut[src_pos000] += target * src_weight000[..., None]
# lut[src_pos100] += target * src_weight100[..., None]
# lut[src_pos010] += target * src_weight010[..., None]
# lut[src_pos001] += target * src_weight001[..., None]
# lut[src_pos110] += target * src_weight110[..., None]
# lut[src_pos101] += target * src_weight101[..., None]
# lut[src_pos011] += target * src_weight011[..., None]
# lut[src_pos111] += target * src_weight111[..., None]
#
# lut_weight[src_pos000] += src_weight000[..., None]
# lut_weight[src_pos100] += src_weight100[..., None]
# lut_weight[src_pos010] += src_weight010[..., None]
# lut_weight[src_pos001] += src_weight001[..., None]
# lut_weight[src_pos110] += src_weight110[..., None]
# lut_weight[src_pos101] += src_weight101[..., None]
# lut_weight[src_pos011] += src_weight011[..., None]
# lut_weight[src_pos111] += src_weight111[..., None]
# #print(lut[src_pos000].shape, lut_weight[src_pos111].shape)
# lut[src_pos000] /= lut_weight[src_pos000]
# lut[src_pos100] /= lut_weight[src_pos100]
# lut[src_pos010] /= lut_weight[src_pos010]
# lut[src_pos001] /= lut_weight[src_pos001]
# lut[src_pos110] /= lut_weight[src_pos110]
# lut[src_pos101] /= lut_weight[src_pos101]
# lut[src_pos011] /= lut_weight[src_pos011]
# lut[src_pos111] /= lut_weight[src_pos111]
#
# weight = np.hstack((src_weight000.reshape(-1, 1) / lut_weight[src_pos000], src_weight100.reshape(-1, 1) / lut_weight[src_pos100]))
# print(src.shape, src_ind000.shape, src_pos000.shape, weight.shape)
# tt = np.hstack((src, src_ind000, src_pos000.reshape(-1, 1), weight))
# print(&#39;tt, src, ind, pos, weight :&#39;)
#
# np.set_printoptions(precision=3, suppress=True)
# #print(tt)
# np.savetxt(&#39;tt.txt&#39;, tt, fmt=&#39;%.4f&#39;, delimiter=&#39; &#39;)
# fig = plt.figure()
# ax1 = fig.add_subplot(1, 3, 1)
# ax1.imshow(src_pos000.reshape(40, 60)[..., None])
# ax2 = fig.add_subplot(1, 3, 2)
# ax2.imshow(src_weight000.reshape(40, 60)[..., None])
# ax3 = fig.add_subplot(1, 3, 3)
# ax3.imshow((target - src).reshape(40, 60, 3))
# plt.show()
return lut
def cal_invdis(src, src_ind, method=&#34;inv&#34;):
&#34;&#34;&#34;
1. method = gauss:
distance = sqrt(x*x + y*y + z*z)
inv_dis = 1 / std / sqrt(2*pi) * exp(-(distance**2)/(2*std*std))
2. medhod = inv:
distance = sqrt(x*x + y*y + z*z)
inv_dis = 1 / (distance+eps)
:param src:
:param src_ind:
:param method:
:return:
&#34;&#34;&#34;
if method == &#39;gauss&#39;:
std = 1
distance2 = np.sum((src - src_ind)**2, axis=-1)
inv_dis = 1 / (std * np.sqrt(2*np.pi)) * np.exp(-distance2/(2*std*std))
if method == &#39;inv&#39;:
distance2 = np.sum((src - src_ind) ** 2, axis=-1)
inv_dis = 1 / (np.sqrt(distance2) + 0.0000001)
#inv_dis2 = np.ones_like(inv_dis)
return inv_dis
//代码来自于CSDN博主：tony365
在大致明白了Lut图的作用之后，回到大气散射的Lut计算。
大气散射Lut预计算：
首先自然就是在c#脚本当中去声明这些Lut参数：
private RenderTexture _particleDensityLUT = null;
private Texture2D _randomVectorsLUT = null;
private RenderTexture _skyboxLUT;
private RenderTexture _skyboxLUT2;
private RenderTexture _inscatteringLUT;
private RenderTexture _extinctionLUT;
private Color[] _directionalLightLUT;
private Color[] _ambientLightLUT;
初始化散射的lut函数：
private void InitializeInscatteringLUT()
{
_inscatteringLUT = new RenderTexture((int)_inscatteringLUTSize.x, (int)_inscatteringLUTSize.y, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Linear);
_inscatteringLUT.volumeDepth = (int)_inscatteringLUTSize.z;
_inscatteringLUT.isVolume = true;
_inscatteringLUT.enableRandomWrite = true;
_inscatteringLUT.name = &#34;InscatteringLUT&#34;;
_inscatteringLUT.Create();
_extinctionLUT = new RenderTexture((int)_inscatteringLUTSize.x, (int)_inscatteringLUTSize.y, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Linear);
_extinctionLUT.volumeDepth = (int)_inscatteringLUTSize.z;
_extinctionLUT.isVolume = true;
_extinctionLUT.enableRandomWrite = true;
_extinctionLUT.name = &#34;ExtinctionLUT&#34;;
_extinctionLUT.Create();
}
然后在CompuerShader计算完之后自然就要创建Lut图作用于Skybox上：
private void PrecomputeSkyboxLUT()//与预计算Lut作用于skybox
{
if (_skyboxLUT == null)
{
_skyboxLUT = new RenderTexture((int)_skyboxLUTSize.x, (int)_skyboxLUTSize.y, 0, RenderTextureFormat.ARGBHalf, RenderTextureReadWrite.Linear);
_skyboxLUT.volumeDepth = (int)_skyboxLUTSize.z;
_skyboxLUT.isVolume = true;
_skyboxLUT.enableRandomWrite = true;
_skyboxLUT.name = &#34;SkyboxLUT&#34;;
_skyboxLUT.Create();
}
#if HIGH_QUALITY
if (_skyboxLUT2 == null)
{
_skyboxLUT2 = new RenderTexture((int)_skyboxLUTSize.x, (int)_skyboxLUTSize.y, 0, RenderTextureFormat.RGHalf, RenderTextureReadWrite.Linear);
_skyboxLUT2.volumeDepth = (int)_skyboxLUTSize.z;
_skyboxLUT2.isVolume = true;
_skyboxLUT2.enableRandomWrite = true;
_skyboxLUT2.name = &#34;SkyboxLUT2&#34;;
_skyboxLUT2.Create();
}
#endif
int kernel = ScatteringComputeShader.FindKernel(&#34;SkyboxLUT&#34;);
ScatteringComputeShader.SetTexture(kernel, &#34;_SkyboxLUT&#34;, _skyboxLUT);
#if HIGH_QUALITY
ScatteringComputeShader.SetTexture(kernel, &#34;_SkyboxLUT2&#34;, _skyboxLUT2);
#endif
UpdateCommonComputeShaderParameters(kernel);
ScatteringComputeShader.Dispatch(kernel, (int)_skyboxLUTSize.x, (int)_skyboxLUTSize.y, (int)_skyboxLUTSize.z);
}
前面进行完lut的创建之后，自然就是要把CompueterShader要用的数据从scripts从脚本传递给ComputerShader。撰写ComputerShaderParameters函数：
private void ComputeShaderParameters(int kernel)
{
ScatteringComputeShader.SetTexture(kernel, &#34;_ParticleDensityLUT&#34;, _particleDensityLUT);
ScatteringComputeShader.SetFloat(&#34;_AtmosphereHeight&#34;, AtmosphereHeight);
ScatteringComputeShader.SetFloat(&#34;_PlanetRadius&#34;, PlanetRadius);
ScatteringComputeShader.SetVector(&#34;_DensityScaleHeight&#34;, DensityScale);
ScatteringComputeShader.SetVector(&#34;_ScatteringR&#34;, RayleighSct * RayleighScatterCoef);
ScatteringComputeShader.SetVector(&#34;_ScatteringM&#34;, MieSct * MieScatterCoef);
ScatteringComputeShader.SetVector(&#34;_ExtinctionR&#34;, RayleighSct * RayleighExtinctionCoef);
ScatteringComputeShader.SetVector(&#34;_ExtinctionM&#34;, MieSct * MieExtinctionCoef);
ScatteringComputeShader.SetVector(&#34;_IncomingLight&#34;, IncomingLight);
ScatteringComputeShader.SetFloat(&#34;_MieG&#34;, MieG);
}
这部分要传递的参数，根据上文的原理对应就明白应该传递什么参数了。
最后就要进行Lut图的更新，说是更新Lut图，其实是在更新ComputerShader中的参数：
private void UpdateInscatteringLUT()//更新Lut图中的参数
{
int kernel = ScatteringComputeShader.FindKernel(&#34;InscatteringLUT&#34;);
ScatteringComputeShader.SetTexture(kernel, &#34;_InscatteringLUT&#34;, _inscatteringLUT);
ScatteringComputeShader.SetTexture(kernel, &#34;_ExtinctionLUT&#34;, _extinctionLUT);
ScatteringComputeShader.SetVector(&#34;_InscatteringLUTSize&#34;, _inscatteringLUTSize);
ScatteringComputeShader.SetVector(&#34;_BottomLeftCorner&#34;, _FrustumCorners[0]);
ScatteringComputeShader.SetVector(&#34;_TopLeftCorner&#34;, _FrustumCorners[1]);
ScatteringComputeShader.SetVector(&#34;_TopRightCorner&#34;, _FrustumCorners[2]);
ScatteringComputeShader.SetVector(&#34;_BottomRightCorner&#34;, _FrustumCorners[3]);
ScatteringComputeShader.SetVector(&#34;_CameraPos&#34;, transform.position);
ScatteringComputeShader.SetVector(&#34;_LightDir&#34;, Sun.transform.forward);
ScatteringComputeShader.SetFloat(&#34;_DistanceScale&#34;, DistanceScale);
UpdateCommonComputeShaderParameters(kernel);
ScatteringComputeShader.Dispatch(kernel, (int)_inscatteringLUTSize.x, (int)_inscatteringLUTSize.y, 1);
}
球体相交：
不知道读者，是否还记得做大气散射，需要先进行球体相交，当然从这里开始的函数都是在ComputerShader中撰写：
float2 RaySphereIntersection(float3 rayOrigin, float3 rayDir, float3 sphereCenter, float sphereRadius)
{
rayOrigin -= sphereCenter;
float a = dot(rayDir, rayDir);
float b = 2.0 * dot(rayOrigin, rayDir);
float c = dot(rayOrigin, rayOrigin) - (sphereRadius * sphereRadius);
float d = b * b - 4 * a * c;
if (d < 0)
{
return -1;
}
else
{
d = sqrt(d);
return float2(-b - d, -b + d) / (2 * a);
}
}//球体相交函数计算大气强度：
计算大气密度就是在做指数衰减。
void GetAtmosphereDensity(float3 position, float3 planetCenter, float3 lightDir, out float2 localDensity, out float2 densityToAtmTop)
{
float height = length(position - planetCenter) - _PlanetRadius;
localDensity = exp(-height.xx / _DensityScaleHeight.xy);
float cosAngle = dot(normalize(position - planetCenter), -lightDir.xyz);
//densityToAtmTop = _ParticleDensityLUT.SampleLevel(PointClampSampler, float2(cosAngle * 0.5 + 0.5, (height / _AtmosphereHeight)), 1.0).xy;
densityToAtmTop = _ParticleDensityLUT.SampleLevel(LinearClampSampler, float2(cosAngle * 0.5 + 0.5, (height / _AtmosphereHeight)), 0.0).xy;
}大气层散射计算：
在获取到了大气层的过程以及得到了大气密度之后，就要进行大气的散射和衰减计算：
void ComputeLocalInscattering(float2 localDensity, float2 densityPA, float2 densityCP, out float3 localInscatterR, out float3 localInscatterM)
{
float2 densityCPA = densityCP + densityPA;
float3 Tr = densityCPA.x * _ExtinctionR;
float3 Tm = densityCPA.y * _ExtinctionM;
float3 extinction = exp(-(Tr + Tm));
localInscatterR = localDensity.x * extinction;
localInscatterM = localDensity.y * extinction;
}相位函数：
计算得到的值只要加上了相位函数的计算，才是真的完成了大气散射的瑞丽和米氏计算。
void ApplyPhaseFunction(inout float3 scatterR, inout float3 scatterM, float cosAngle)
{
// r
float phase = (3.0 / (16.0 * PI)) * (1 + (cosAngle * cosAngle));
scatterR *= phase;
// m
float g = _MieG;
float g2 = g * g;
phase = (1.0 / (4.0 * PI)) * ((3.0 * (1.0 - g2)) / (2.0 * (2.0 + g2))) * ((1 + cosAngle * cosAngle) / (pow((1 + g2 - 2 * g*cosAngle), 3.0 / 2.0)));
scatterM *= phase;
}//使用相位函数计算总散射结果：
计算完各自的散射之后，自然就是要得到整体的散射结果
ScatteringOutput IntegrateInscattering(float3 rayStart, float3 rayDir, float rayLength, float3 planetCenter, float3 lightDir)
{
float sampleCount = 64;
float3 step = rayDir * (rayLength / sampleCount);
float stepSize = length(step);
float2 densityCP = 0;
float3 scatterR = 0;
float3 scatterM = 0;
float2 localDensity;
float2 densityPA;
float2 prevLocalDensity;
float3 prevLocalInscatterR, prevLocalInscatterM;
GetAtmosphereDensity(rayStart, planetCenter, lightDir, prevLocalDensity, densityPA);
ComputeLocalInscattering(prevLocalDensity, densityPA, densityCP, prevLocalInscatterR, prevLocalInscatterM);
// P - current integration point
// C - camera position
// A - top of the atmosphere
[loop]
for (float s = 1.0; s < sampleCount; s += 1)
{
float3 p = rayStart + step * s;
GetAtmosphereDensity(p, planetCenter, lightDir, localDensity, densityPA);
densityCP += (localDensity + prevLocalDensity) * (stepSize / 2.0);
prevLocalDensity = localDensity;
float3 localInscatterR, localInscatterM;
ComputeLocalInscattering(localDensity, densityPA, densityCP, localInscatterR, localInscatterM);
scatterR += (localInscatterR + prevLocalInscatterR) * (stepSize / 2.0);
scatterM += (localInscatterM + prevLocalInscatterM) * (stepSize / 2.0);
prevLocalInscatterR = localInscatterR;
prevLocalInscatterM = localInscatterM;
}
ScatteringOutput output;
output.rayleigh = scatterR;
output.mie = scatterM;
return output;
}skyboxLut计算：
散射已经计算完了，自然就要计算skybox了：
void SkyboxLUT(uint3 id : SV_DispatchThreadID)
{
float w, h, d;
_SkyboxLUT.GetDimensions(w, h, d);
float3 coords = float3(id.x / (w - 1), id.y / (h - 1), id.z / (d - 1));
float height = coords.x * coords.x * _AtmosphereHeight;
float ch = -(sqrt(height * (2 * _PlanetRadius + height)) / (_PlanetRadius + height));
float viewZenithAngle = coords.y;
if (viewZenithAngle > 0.5)
{
viewZenithAngle = ch + pow((viewZenithAngle-0.5)*2, 5) * (1 - ch);
}
else
{
viewZenithAngle = ch - pow(viewZenithAngle*2, 5) * (1 + ch);
}
float sunZenithAngle = (tan((2 * coords.z - 1 + 0.26)*0.75)) / (tan(1.26 * 0.75));// coords.z * 2.0 - 1.0;
float3 planetCenter = float3(0, -_PlanetRadius, 0);
float3 rayStart = float3(0, height, 0);
float3 rayDir = float3(sqrt(saturate(1 - viewZenithAngle * viewZenithAngle)), viewZenithAngle, 0);
float3 lightDir = -float3(sqrt(saturate(1 - sunZenithAngle * sunZenithAngle)), sunZenithAngle, 0);
float rayLength = 1e20;
float2 intersection = RaySphereIntersection(rayStart, rayDir, planetCenter, _PlanetRadius + _AtmosphereHeight);
rayLength = intersection.y;
intersection = RaySphereIntersection(rayStart, rayDir, planetCenter, _PlanetRadius);
if (intersection.x > 0)
rayLength = min(rayLength, intersection.x);
ScatteringOutput scattering = IntegrateInscattering(rayStart, rayDir, rayLength, planetCenter, lightDir);
_SkyboxLUT[id.xyz] = float4(scattering.rayleigh.xyz, scattering.mie.x);
#ifdef HIGH_QUALITY
_SkyboxLUT2[id.xyz] = scattering.mie.yz;
#endif
}计算大气散射：
相当于我们现在完成了对于光照的衰减过程计算，散射过程计算，计算结果之后的总体散射计算，散射函数的计算，以及最终的skybox计算。但是我们还没有进行对于上面过程总和计算，表达一个大气散射过程的整体描述，那么最为核心的散射函数为：
void PrecomputeLightScattering(float3 rayStart, float3 rayDir, float rayLength, float3 planetCenter, float3 lightDir, uint3 coords, uint sampleCount)
{
float3 step = rayDir * (rayLength / (float)(sampleCount-1));
float stepSize = length(step) * _DistanceScale;
float2 densityCP = 0;
float3 scatterR = 0;
float3 scatterM = 0;
float2 localDensity;
float2 densityPA;
float2 prevLocalDensity;
float3 prevLocalInscatterR, prevLocalInscatterM;
GetAtmosphereDensity(rayStart, planetCenter, lightDir, prevLocalDensity, densityPA);
ComputeLocalInscattering(prevLocalDensity, densityPA, densityCP, prevLocalInscatterR, prevLocalInscatterM);
_InscatteringLUT[coords] = float4(0, 0, 0, 1);
_ExtinctionLUT[coords] = float4(1, 1, 1, 1);
// P - current integration point
// C - camera position
// A - top of the atmosphere
[loop]
for (coords.z = 1; coords.z < sampleCount; coords.z += 1)
{
float3 p = rayStart + step * coords.z;
GetAtmosphereDensity(p, planetCenter, lightDir, localDensity, densityPA);
densityCP += (localDensity + prevLocalDensity) * (stepSize / 2.0);
prevLocalDensity = localDensity;
float3 localInscatterR, localInscatterM;
ComputeLocalInscattering(localDensity, densityPA, densityCP, localInscatterR, localInscatterM);
scatterR += (localInscatterR + prevLocalInscatterR) * (stepSize / 2.0);
scatterM += (localInscatterM + prevLocalInscatterM) * (stepSize / 2.0);
prevLocalInscatterR = localInscatterR;
prevLocalInscatterM = localInscatterM;
float3 currentScatterR = scatterR;
float3 currentScatterM = scatterM;
ApplyPhaseFunction(currentScatterR, currentScatterM, dot(rayDir, -lightDir.xyz));
float3 lightInscatter = (currentScatterR * _ScatteringR + currentScatterM * _ScatteringM) * _IncomingLight.xyz;
float3 lightExtinction = exp(-(densityCP.x * _ExtinctionR + densityCP.y * _ExtinctionM));
_InscatteringLUT[coords] = float4(lightInscatter, 1);
_ExtinctionLUT[coords] = float4(lightExtinction, 1);
}
}预计算收尾：
接下来在一切计算过程，表述完整之后，自然就要用要编译的主函数去进行大气散射的预计算：
void InscatteringLUT(uint3 id : SV_DispatchThreadID)
{
float w, h, d;
_InscatteringLUT.GetDimensions(w, h, d);
float2 coords = float2(id.x / (w - 1), id.y / (h - 1));
float3 v1 = lerp(_BottomLeftCorner.xyz, _BottomRightCorner.xyz, coords.x);
float3 v2 = lerp(_TopLeftCorner.xyz, _TopRightCorner.xyz, coords.x);
float3 rayEnd = lerp(v1, v2, coords.y);
float3 rayStart = _CameraPos.xyz;
float3 rayDir = rayEnd - rayStart;
float rayLength = length(rayDir);
rayDir /= rayLength;
float3 planetCenter = float3(0, -_PlanetRadius, 0);
PrecomputeLightScattering(rayStart, rayDir, rayLength, planetCenter, normalize(_LightDir), id, d);
}结尾：
最后你会得到这样的一张图：
也许你会说这里为什么和我上面的图不一样，这是因为没有进行tonemapping。而关于tonemapping的内容就在其余文章讲述了。
如果你在加上lightshift效果，就能得到很好看的光照效果，但是需要说明的是，关于这里的内容由于写作本文耗费的精力实在有点大，所以笔者会根据读者的要求来看是不是要进行这部分内容的讲述和迭代，开源工程将会在笔者完成了对于代码的梳理和注释的编写后上传github，到时候会迭代在本文当中。最后就用笔者实现的大气散射效果图镇文吧
如果对于本文感觉阅读有一定压力，可以告诉笔者，笔者会进行迭代。祝愿读者在读完本文之后，收获了自己的大气散射，那么最后我是在野，感谢你的阅读。
参考链接：
https://blog.csdn.net/weixin_43803133/article/details/116354462
https://zhuanlan.zhihu.com/p/463551881
https://zhuanlan.zhihu.com/p/259316527
https://www.zhihu.com/search?type=content&q=%E9%AB%98%E6%96%AF%E5%AE%9A%E5%BE%8B%E7%9A%84%E6%8E%A8%E5%AF%BC
https://wuli.wiki/online/IntN.html
https://www.zhihu.com/question/326568092/answer/860302644
https://zhuanlan.zhihu.com/p/421185092
https://zhuanlan.zhihu.com/p/43020220
https://zhuanlan.zhihu.com/p/123835582
https://zhuanlan.zhihu.com/p/351289217
https://www.azooptics.com/Article.aspx?ArticleID=790
https://zhuanlan.zhihu.com/p/363866197
http://download.nvidia.com/developer/SDK/Individual_Samples/DEMOS/Direct3D9/src/HLSL_RainbowFogbow/docs/RainbowFogbow.pdf
http://nis-lab.is.su-tokyo.ac.jp/~nis/cdrom/pg/pg01_lightning.pdf
http://www.markmark.net/PDFs/RTClouds_HarrisEG2001.pdf
http://www.ati.com/developer/dx9/ATI-LightScattering.pdf
http://nis-lab.is.su-tokyo.ac.jp/~nis/cdrom/sig93_nis.pdf
http://www.gamedev.net/columns/hardcore/atmscattering/
https://hello-david.github.io/archives/f116047e.html
https://www.zhangxinxu.com/wordpress/2020/02/3d-lut-principle/
[https://blog.csdn.net/tywwwww/article/details/122338432
参考论文：
01.《MieTheory》-米氏散射的推导和研究
02.《Beitrge zur Optik trüber Medien, speziell kolloidaler Metallsungen》-Gustay Mie1908年米氏散射理论提出
03.《Precomputed Atmospheric Scattering》-Eric Brunetion,Fabrice Neyret
04.《An Efficient 3D Color LUT Compression Algorithm Based on a Multi-Scale Anisotropic Diffusion Scheme》
参考书籍：
《GPU GEM2》-16章：精确大气散射
致谢：
感谢所有参考链接的作者，文章的题图和知识离不开各位的付出。本文旨在对于知识的学习和传播，不作任何商业用途，希望发展一个良好的社区学习氛围。完善一个系统的大气散射过程，也是我对于游戏世界真实大气的追求。 |
本帖子中包含更多资源
您需要 登录 才可以下载或查看，没有账号？立即注册
×
|
窥视卡
雷达卡
发表于 2022-8-2 10:06
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
千斤顶
显身卡