大气散射

资料

原理

总览：
06.游戏中地形大气和云的渲染(下) | GAMES104-现代游戏引擎：从入门到实践
一个简短的公式推导：
[Rendering] 基于物理的大气渲染
一个详细的推导和最暴力且直观的实现：
【实战】从零实现一套完整单次大气散射
$I=I_S\beta(\lambda)\gamma(\theta)\sum_{P\in\overline{AB}}T(\overline{CP})T(\overline{PA})\rho(h)ds$
$T(\overline{CP})T(\overline{PA})=exp\lbrace-\beta(\lambda)(D(\overline{CP})+D(\overline{PA}))\rbrace$
这篇文章也不错：
基于物理的大气体渲染从理论到实践（一）
天空为什么白天是蓝的，黄昏是红的：
更通俗易懂之天空为啥那么蓝——瑞利散射
一个有趣的实验：
麻省理工教授：答应我，如果不是为瑞利散射，永远不要同时抽四根烟！！！
一篇详尽过头的文章，我的建议是当小说看：
游戏魔法编程：unity 实现完整大气散射

Accurate Atmospheric Scattering

只考虑单次散射，预计算并拟合了任意两点间的光学距离：
GPU Gems 2 Chapter 16. Accurate Atmospheric Scattering
这也是 Unity 的做法：
在 Unity 中制作 Procedural Sky
解除了 Sean O’Neil 实现中标准海拔、大气高度和地球半径之间比例的限制。
Flexible Physical Accurate Atmosphere Scattering

Precomputed Atmospheric Scattering

预计算了多重散射，考虑了相机在太空的情况，考虑了臭氧层对光线的吸收，考虑了太阳对着色点 Direct Irradiance 的贡献，考虑了大气对着色点 Indirect Irradiance 的贡献，考虑了看向着色点时视线光路上发生的散射，考虑了位于阴影中的大气部分产生的 LightShafts:
Precomputed Atmospheric Scattering
Precomputed Atmospheric Scattering
Precomputed Atmosphere Scaterring
【Unreal 从 0 到 1】【第四章：天气系统】4.2，从零开始的物理大气散射
 解读 GLSL 的大气散射函数学习笔记
 Precomputed Atmospheric Scattering Shader 代码笔记
 CPU 精确计算与 GPU 近似的对比
 对不同近似的对比的描述 Test cases 部分

如何在 win + VS2022 编译该工程

ebruneton/precomputed_atmospheric_scattering
文档中的 make demo 是 linux 的编译方式，windows 的编译方式藏在：
precomputed_atmospheric_scattering/platform/windows/。
download_build_run.bat 会依次运行 download_dependencies.bat，generate_project.bat，build.bat，run.bat。中间如果出错会直接关闭命令提示符窗口，手动一个个试才能知道具体报错。
generate_project.bat
报错：

CMake Error at CMakeLists.txt:78 (project):
Failed to run MSBuild command:
    MSBuild.exe
to get the value of VCTargetsPath:
-- Configuring incomplete, errors occurred!

将 MSBuild.exe 加入环境变量即可，我的在 C:\Program Files\Microsoft Visual Studio\2022\Community\Msbuild\Current\Bin。
报错：

1
2
3

无法找到 Visual Studio 2017 的生成工具(平台工具集 =“v141”)。若要使用 v141 生成工具进行生成，
请安装 Visual Studio 2017 生成工具。或者，可以升级到当前 Visual Studio 工具，方式是通过选择
“项目”菜单或右键单击该解决方案，然后选择“重定解决方案目标”。

VisualStudio Installer，修改，单个组件，搜索 2017，勾选 MSVC v141 - VS2017 C++ x64/x86 生成工具(v14.16)。
报错：

1	`无法启动程序 /ALL_BUILD 拒绝访问`

右键项目设为启动项目。

关于几种新做法的系列文章：
大气散射相关

这篇是 Eric Bruneton 那篇的进一步改进：
Rendering Parametrizable Planetary Atmospheres with Multiple Scattering in Real-Time
【译】【大气散射】[Elek09] 实时渲染参数化的、有多次散射的行星大气

寒霜引擎的实现，Oskar Elek 那篇进一步改进：
Physically Based Sky, Atmosphere & Cloud Rendering in Frostbite 2016 那篇
【译】【寒霜引擎】【大气渲染】基于物理的天空、大气和云的渲染（上）

英特尔的实现，包括完整的大气+地形+阴影+后处理。
Outdoor Light Scattering Update
GameTechDev/OutdoorLightScattering

纯拟合：
An Analytic Model for Full Spectral Sky-Dome Radiance 最后一篇

实现

Single Scattering

$I=I_S\beta(\lambda)\gamma(\theta)\sum_{P\in\overline{AB}}T(\overline{CP})T(\overline{PA})\rho(h)ds$
$T(\overline{CP})T(\overline{PA})=exp\lbrace-\beta(\lambda)(D(\overline{CP})+D(\overline{PA}))\rbrace$
单个 Fragment Shader 即可实现，LightDir 和 CameraPos 由 uniform 传入，viewDir 由 skyBox 的坐标得到。

$input v_worldPos

uniform vec4 u_cameraPos[1];
uniform vec4 u_LightDir[1];

// 单位：km
#define _PlanetRadius 6371.0
#define _AtmosphereHeight 80.0
#define _DensityScaleHeight vec2(7.994, 1.2)
#define _PlanetCenter vec3(0.0, -_PlanetRadius, 0.0)

#define _ExtinctionR vec3(0.0058, 0.0135, 0.0331)
#define _ExtinctionM vec3(0.0200, 0.0200, 0.0200)

#define _MieG 0.98
#define _IncomingLight vec3(12.0, 12.0, 12.0)
#define _PI 3.1415926

float GetAltitude(vec3 position){
    return abs(length(position - _PlanetCenter) - _PlanetRadius);
}

// 射线球面相交函数
vec2 RaySphereIntersection(vec3 rayOrigin, vec3 rayDir,
    vec3 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.0 * a * c;
    if (d < 0.0) {
    	return vec2(-1.0, -1.0);
    }
    else {
    	d = sqrt(d);
    	return vec2(-b - d, -b + d) / (2.0 * a);
    }
}

vec2 RayEarthIntersection(vec3 rayOrigin, vec3 rayDir) {
    return RaySphereIntersection(rayOrigin, rayDir, _PlanetCenter, _PlanetRadius);
}

vec2 RayAtmosphereIntersection(vec3 rayOrigin, vec3 rayDir) {
    return RaySphereIntersection(rayOrigin, rayDir, _PlanetCenter,  
        _PlanetRadius + _AtmosphereHeight);
}

// 采样得到两点之间的光学深度。
vec2 DensitySampling(vec3 startPos, vec3 endPos) {
	vec2 opticalDepth = vec2_splat(0.0);
	
	const int stepCount = 16;
	vec3 step = (endPos - startPos) / stepCount;
	float stepSize = length(step);
	
	for(int s = 0; s <= stepCount; ++s) {
		vec3 position = startPos + step * s;
		float height = GetAltitude(position);
		vec2 localDensity = exp(-vec2_splat(height) / _DensityScaleHeight.xy);
		
		// 梯形法则
		float weight = (s == 0 || s == stepCount) ? 0.5 : 1.0;
		opticalDepth += localDensity * stepSize * weight;
	}
	
	return opticalDepth;
}

// 射线起点到射线与大气顶部交点之间的光学深度。
vec2 LightDensitySampleing(vec3 position, vec3 lightDir) {
    // 在这里假设射线一定不与地表相交。
	vec2 atmIntersection = RayAtmosphereIntersection(position, lightDir);
	vec3 rayEnd = position + lightDir * atmIntersection.y;
	return DensitySampling(position, rayEnd);
}

// T(PC) * T(PA) * rho(h)
void ComputeLocalInscattering(vec2 rho, vec2 densityCP, vec2 densityPA,
    out vec3 localInscatterR, out vec3 localInscatterM) {
	vec2 densityCPA = densityCP + densityPA;
	
	vec3 Tr = densityCPA.x * _ExtinctionR;
	vec3 Tm = densityCPA.y * _ExtinctionM;

	vec3 extinction = exp(-(Tr + Tm));
		
	localInscatterR = extinction * rho.x;
	localInscatterM = extinction * rho.y;
}

// 相函数。
void ApplyPhaseFunction(inout vec3 scatterR, inout vec3 scatterM, float cosAngle) {
	float cosAngle2  = cosAngle * cosAngle;
	float g = _MieG;
	float g2 = g * g;
	
	float phase = (3.0 / (16.0 * _PI)) * (1.0 + cosAngle2);
	scatterR *= phase;

	phase = (1.0 / (4.0 * _PI)) *
        ((3.0 * (1.0 - g2)) / (2.0 * (2.0 + g2))) *
        ((1.0 + cosAngle2) / (pow(abs((1.0 + g2 - 2.0 * g * cosAngle)), 1.5)));
	scatterM *= phase;
}

vec3 IntegrateInscattering(vec3 rayStart, vec3 rayDir, float rayLength, vec3 lightDir) {
	// 在循环中累加
	vec3 scatterR  = vec3_splat(0.0);
	vec3 scatterM  = vec3_splat(0.0);
	vec2 densityPA = vec2_splat(0.0);
	vec2 densityCP = vec2_splat(0.0);
	
	// P - 当前发生散射的位置
	// A - 相机位置
	// C - 与大气顶部的交点
	const int sampleCount = 256;
	vec3 step = rayDir * (rayLength / sampleCount);
	float stepSize = length(step);
	for (int s = 0; s <= sampleCount; ++s) {
		vec3 samplePos = rayStart + step * s;
		
        // 这里只考虑单次散射，所以如果一个采样点对太阳不可见
        // 他不会收到任何 Radiance 自然也不会散射出任何 Radiance。
		vec2 crtEarthIntersection = RayEarthIntersection(samplePos, lightDir);
		bool vis = (crtEarthIntersection.x > 0.0 || crtEarthIntersection.y > 0.0)
            ? false : true;
		
		if(vis) {
			densityCP = LightDensitySampleing(samplePos, lightDir);
			
            // DensityPA 可以在每一次采样点 P 步进的时候累加，这样可以省下 PA 之间的积分。
			float height = GetAltitude(samplePos);
			vec2 localDensityPA = exp(-vec2_splat(height) / _DensityScaleHeight.xy);
			densityPA += localDensityPA * stepSize;
			
			// localDensityPA 便是 rho(h)。
			vec3 localInscatterR;
			vec3 localInscatterM;
			ComputeLocalInscattering(localDensityPA, densityCP, densityPA,
            localInscatterR, localInscatterM);
			
			// 梯形法则
			float weight = (s == 0 || s == sampleCount) ? 0.5 : 1.0;
			scatterR += localInscatterR * stepSize * weight;
			scatterM += localInscatterM * stepSize * weight;
		}
	}
	ApplyPhaseFunction(scatterR, scatterM, dot(rayDir, lightDir));
	return (scatterR * _ExtinctionR + scatterM * _ExtinctionM) * _IncomingLight.xyz;
}

void main()
{
	vec3 rayStart = u_cameraPos[0].xyz;
	vec3 rayDir = normalize(v_worldPos.xyz);
	
    // 当相机位于太空时，沿着视线方向将相机移动至大气顶部。
	vec2 atmIntersection = RayAtmosphereIntersection(rayStart, rayDir);
	if(atmIntersection.x > 0.0) {
		rayStart +=  rayDir * atmIntersection.x;
		atmIntersection = RayAtmosphereIntersection(rayStart, rayDir);
	}
	float rayLength = atmIntersection.y;
	
	// 如果射线与地表相交，更新射线的长度。
	vec2 earthIntersection = RayEarthIntersection(rayStart, rayDir);
	if (earthIntersection.x > 0.0 || earthIntersection.y > 0.0) {
		rayLength = earthIntersection.x;
	}
	
	vec3 lightDir = -normalize(u_LightDir[0].xyz);
	vec3 radiance = IntegrateInscattering(rayStart, rayDir, rayLength, lightDir);
	
	gl_FragColor = vec4(radiance, 1.0);
}

Accurate Atmospheric Scattering

GPU Gems 2 Chapter 16. Accurate Atmospheric Scattering
提出了预计算光学深度的做法，高度视线天顶角夹角作为参数得到一张二维的 LUT，代表了任意一点到大气顶部的光学深度。如果要计算 AB 两点之间的光学深度，假设射线 AB 与大气顶部的交点为 C，Density(AB) 即 Density(AC) - Density(BC)。如果射线 AB 指向地表，因为 Density(AB) == Density(BA)，这时计算 BA 两点之间的光学深度即可。
然后作者又将该 LUT 进行了拟合，公式如下：

float GetDensityFit(float height, float cosine) {
    float heightTerm = exp(-4.0 * height);
    float c = 1.0 - cosine;
    float cosineTerm =
        0.25 * exp(-0.00287 + c * (0.459 + c * (3.38 + c * (-6.8 + c * 5.25))));
    return heightTerm * cosineTerm;
}

Precomputed Atmospheric Scattering

内容比较多，建议直接看 Precomputed Atmospheric Scattering，这里只有一些笔记。
如果仅仅要实现大气散射，这篇文档包含了不少冗余的内容：

额外提供了全光谱的渲染模式，而非只使用 RGB 对应的三种波长。
为防止不合理的类型转换，自定义了所有物理量，但是 glsl 没有静态类型检查，所以所有物理量实际上都是 float 或 vec3。
为验证各个 trick 的效果，额外提供了 CPU 版本的渲染模式，通过 #include "atmosphere/functions.glsl" 和一些宏定义，使得 CPU、GPU 共用同一份代码。

所以 demo 和 reference 下的文件是基本不用看的。
还有一些：

没有实际的 Shader 文件，所有 Shader 代码是在 model 类里用 string 拼出来的。
model 类的初始化包括了最重要的 AtmosphereParameters 结构体的初始化。
roeas/PrecomputedAtmosphereTexture 是一个用来验证预计算贴图可行性的小工程，只计算了 Transmittance 并将其保存为图片，其中还包括了如何生成定义 AtmosphereParameters 对象的 GLSL 代码。
Sharer 里同一张 Texture 名有时代表的是最终存储的 LUT，有时代表的是临时的 delta LUT，具体要看 CUP 端绑的是哪张。为防止混淆我做了一点整理：

存储并使用的贴图

transmittance_texture_
scattering_texture_
irradiance_texture_
optional_single_mie_scattering_texture_（可选项）

临时贴图

delta_irradiance_texture
delta_rayleigh_scattering_texture（仅用于计算二阶散射）
delta_mie_scattering_texture（仅用于计算二阶散射）
delta_scattering_density_texture
delta_multiple_scattering_texture = delta_rayleigh_scattering_texture（复用，仅用于计算三阶及以上的散射）

Transmittance

输入

输出

transmittance_texture_

direct Irradiance

输入

transmittance_texture_

输出

delta_irradiance_texture
irradiance_texture_ = 0.0

single Scattering

输入

transmittance_texture_

输出

delta_rayleigh_scattering_texture
delta_mie_scattering_texture
scattering_texture_ = rayleigh.rgb + mie.r
optional_single_mie_scattering_texture_ = mie（可选项）

Scattering Density

输入

transmittance_texture_
delta_rayleigh_scattering_texture
delta_mie_scattering_texture
delta_multiple_scattering_texture
delta_irradiance_texture

输出

delta_scattering_density_texture

indirect Irradiance

输入

delta_rayleigh_scattering_texture
delta_mie_scattering_texture
delta_multiple_scattering_texture

输出

delta_irradiance_texture
irradiance_texture_

multiple Scattering

输入

transmittance_texture_
delta_scattering_density_texture

输出

delta_multiple_scattering_texture
scattering_texture_

#AtmosphericScattering #GlobalIllumination #Math #Rendering #VolumeRendering

Assimp 上一篇

【Archive】笔记 E 星计划下一篇

大气散射

资料

原理

Accurate Atmospheric Scattering

Precomputed Atmospheric Scattering

如何在 win + VS2022 编译该工程

A Scalable and Production Ready Sky and Atmosphere Rendering Technique

大镖客中的天气

其他资料

实现

Single Scattering

Accurate Atmospheric Scattering

Precomputed Atmospheric Scattering

存储并使用的贴图

临时贴图

Transmittance

direct Irradiance

single Scattering

Scattering Density

indirect Irradiance

multiple Scattering