大气散射学习和实践

原理

学习了大气散射的原理，看了不少文章，讲得不错的就是这几篇
[Rendering] 基于物理的大气渲染
从零实现一套完整单次大气散射
基本的原理明白了。

大气散射求解的是从视线方向AB看去应该是什么颜色，那么应该是所有太阳光经过AB上的点散射到A点方向上的颜色和（这只考虑单次散射）。
解决这个问题，可以先考虑传输方程，光每经过一段距离，就会衰减，假如每单位距离的衰减系数是delta，那么可以列出方程。
$$dL=-\delta Lds$$
这就是经典的存款利息问题，解这个常微分方程就是
$$L=L_0{e}^{-\delta s}$$
如果delta是s的函数，那么上面的方程应该写为
$$L=L_0{e}^{{\int }_A^B-\delta (s)ds}$$
现在的关键问题是在大气中这个detla由哪些因素影响，可以想象，肯定受大气密度影响，其次光波长本身也会有影响（不同的波长被散射的比重不同）。因此detla又可以写为
$$\delta (s)=\beta (\lambda )\rho (h)$$
所以公式可以写为
$$\begin{gathered} L=L_0{e}^{{\int }_A^B-\beta (\lambda )\rho (h)ds}，其中\beta (\lambda )是波长的函数，可以拿出来L=L_0{e}^{-\beta (\lambda ){\int }_A^B\rho (h)ds} \\ {\int }_A^B\rho (h)ds称为光学距离，它只和介质的也密度有关，描述的是介质本身的属性和它通过什么光无关。 \end{gathered}$$

$$而T(AB)=e^{-\beta (\lambda ){\int }_A^B\rho (h)ds}被称为传输方程，表明波长\lambda 的光经过AB剩余的能量$$
$$再来看\rho (h)，空气密度只和海拔高度有关，也是一个负指数形式的公式，\rho (h)=e^{-h/H}$$
意思是每上升单位高度，密度变化量也是当前密度的一个比例。又是经典的复利问题。
如果只是计算直射的光的传输方程，那么上面的方程就够了，但是我们要考虑的是单次散射的行为，因此需要有一个分布函数来描述，光在指定入射方向上散射后，能量在出射方向上的球面分布。那么这个函数就是散射系数，它肯定是光波长，大气密度的函数，同时因为又是一个球面分布，可以想象和球面坐标相关
$$写作S(\lambda ,\rho ,\theta ,\varphi )$$
但我猜可能这个这个函数以入射光方向为x轴的话，是x轴是径向对称的，所以只与天顶角有关与方位角无关，所以函数可以写为
$$写作S(\lambda ,\rho ,\theta )$$
这个函数表示了散射出去的能量，因此对它进行球面积分的话，就应该得到上面的detla函数，刚好，球面积分会消去 theta.

$${\int }_{\Omega }S(\lambda ,\rho ,\theta )dw=I_s(\lambda ,\rho )=\delta (s)=\beta (\lambda )\rho (h)$$

$$S(\lambda ,\rho ,\theta )的积分可以用\beta (\lambda )\rho (h)表示，为了方便，我们把S本身也用\beta 和\rho 来表示$$

$$因此就写成S(\lambda ,\rho ,\theta )=\beta (\lambda )\rho (h)P(\theta )$$

$$其中P(\theta )=\frac{S(\lambda ,\rho ,\theta )}{{\int }_{\Omega }S(\lambda ,\rho ,\theta )dw}=\frac{S(\lambda ,\rho ,\theta )}{\beta (\lambda )\rho (h)}，只和角度有关，被称为相位函数$$

公式拆解完毕，那么光从C到P到A的方程为

这里引用第一篇文章的图，最后的D(PA)应该为D(CP)
最终把视线方向上所有的点加起来，就得到求解颜色。

实践

先给camera挂上脚本如下:

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEngine.Rendering;

[ExecuteInEditMode]
#if UNITY_5_4_OR_NEWER
[ImageEffectAllowedInSceneView]
#endif
public class SkyScatter : MonoBehaviour {


    public Material skyMat;
	protected Mesh mesh;
    protected Camera m_camera;
    protected CommandBuffer skyBoxCmdBuffer;
    public float rDensityScale = 7994.0f;

    public float mDensityScale = 1200;

    public Vector3 rCoef = new Vector3(5.8f, 13.5f, 33.1f);
    public float rScatterStrength = 1f;
    public float rExtinctionStrength = 1f;

    public Vector3 mCoef = new Vector3(2.0f, 2.0f, 2.0f);
    public float mScatterStrength = 1f;
    public float mExtinctionStrength = 1f;
    public float mieG = 0.625f;

    public Vector3 lightDir = new Vector3(-1, -1, 0);

    public Color lightFromOuterSpace = Color.white;

    public float _PlanetRadius = 6357000.0f;
    public float _AtmosphereHeight = 12000.0f;
    void CreateMesh()
    {
        mesh = new Mesh();
        Vector3[] v = new Vector3[]
        {
            new Vector3(-1,-1,1),
            new Vector3(1,-1,1),
            new Vector3(-1,1,1),
            new Vector3(1,1,1),
        };

        int[] triangles = new int[]
        {
            0,1,2,1,3,2
        };

        Vector2[] uv = new Vector2[]
        {
            new Vector2(0, 0),
            new Vector2(1,0),
            new Vector2(0,1),
            new Vector2(1,1),
        };
        mesh.vertices = v;
        mesh.triangles = triangles;
        mesh.uv = uv;
        mesh.UploadMeshData(true);
    }

    void Start()
    {
        m_camera = gameObject.GetComponent<Camera>();

        _SetShaderParam();
        CreateMesh();
        skyBoxCmdBuffer = new CommandBuffer();
        skyBoxCmdBuffer.name = "skyBoxCmdBuffer";
        skyBoxCmdBuffer.DrawMesh(mesh, Matrix4x4.identity, skyMat);

        m_camera.AddCommandBuffer(CameraEvent.BeforeForwardAlpha,  skyBoxCmdBuffer);
    }

    void _SetShaderParam()
    {
        skyMat.SetFloat("_PlanetRadius", _PlanetRadius);
        skyMat.SetFloat("_AtmosphereHeight", _AtmosphereHeight);
        skyMat.SetVector("_DensityScaleHeight", new Vector4(rDensityScale, mDensityScale));
        var rCoefReal = this.rCoef * 0.000001f;
        var mCoefReal = this.mCoef * 0.00001f;
        skyMat.SetVector("_ScatteringR", rCoefReal * rScatterStrength);
        skyMat.SetVector("_ScatteringM", mCoefReal * mScatterStrength);
        skyMat.SetVector("_ExtinctionR", rCoefReal * rExtinctionStrength);
        skyMat.SetVector("_ExtinctionM", mCoefReal * mExtinctionStrength);
        skyMat.SetFloat("_MieG", mieG);
        skyMat.SetVector("_IncomingLight", lightFromOuterSpace);

        skyMat.SetVector("_LightDir", lightDir.normalized);

        var projectionMatrix = GL.GetGPUProjectionMatrix(m_camera.projectionMatrix, false);

        skyMat.SetMatrix("_InverseViewMatrix", m_camera.worldToCameraMatrix.inverse);
        skyMat.SetMatrix("_InverseProjectionMatrix", projectionMatrix.inverse);
    }

    void Update () {

        _SetShaderParam();
    }
}

然后创建一个shader和一个cginc赋给skyMat

Shader "Unlit/SkyBox"
{
	Properties
	{
	}
	SubShader
	{
		Tags { "RenderType"="Opaque" }

		Pass
		{
			CGPROGRAM
			#pragma vertex vert
			#pragma fragment frag
			
			#include "UnityCG.cginc"
			#include "Scatter.cginc"
			
			ENDCG
		}
	}
}

struct appdata
{
	float4 vertex : POSITION;
	float2 uv : TEXCOORD0;
};

struct v2f
{
	float4 vertex : SV_POSITION;
	float2 uv : TEXCOORD0;
};



v2f vert (appdata v)
{
	v2f o;
	o.vertex = v.vertex;
	o.uv = v.uv;

#if UNITY_UV_STARTS_AT_TOP
    o.uv.y = 1 - o.uv.y;
#else
#endif
	return o;
}

#define PI 3.14159265359

float _AtmosphereHeight;
float _PlanetRadius;
float2 _DensityScaleHeight;

float3 _ScatteringR;
float3 _ScatteringM;
float3 _ExtinctionR;
float3 _ExtinctionM;

float4 _IncomingLight;
float _MieG;

float _SunIntensity;
float _DistanceScale;

float3 _LightDir;
sampler2D_float _CameraDepthTexture;
float4 _CameraDepthTexture_ST;

float4x4 _InverseViewMatrix;
float4x4 _InverseProjectionMatrix;

float3 GetWorldSpacePosition(float2 i_UV)
{
	float depth = SAMPLE_DEPTH_TEXTURE(_CameraDepthTexture, i_UV);

	float4 positionViewSpace = mul(_InverseProjectionMatrix, float4(2.0 * i_UV - 1.0, depth, 1.0));
	positionViewSpace /= positionViewSpace.w;


	float3 positionWorldSpace = mul(_InverseViewMatrix, float4(positionViewSpace.xyz, 1.0)).xyz;
	return positionWorldSpace;
}


//-----------------------------------------------------------------------------------------
// Helper Funcs : RaySphereIntersection
//-----------------------------------------------------------------------------------------
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);
	}
}


//----- Input
// position			视线采样点P
// lightDir			光照方向

//----- Output : 
// opticalDepthCP:	dcp
bool lightSampleing(
	float3 position,							// Current point within the atmospheric sphere
	float3 lightDir,							// Direction towards the sun
	out float2 opticalDepthCP)
{
	opticalDepthCP = 0;


	float3 rayStart = position;
	float3 rayDir = -lightDir;

	float3 planetCenter = float3(0, -_PlanetRadius, 0);
	float2 intersection = RaySphereIntersection(rayStart, rayDir, planetCenter, _PlanetRadius + _AtmosphereHeight);
	float3 rayEnd = rayStart + rayDir * intersection.y;

	// compute density along the ray
	float stepCount = 50;// 250;
	float3 step = (rayEnd - rayStart) / stepCount;
	float stepSize = length(step);
	float2 density = 0;

	for (float s = 0.5; s < stepCount; s += 1.0)
	{
		float3 position = rayStart + step * s;
		float height = abs(length(position - planetCenter) - _PlanetRadius);
		float2 localDensity = exp(-(height.xx / _DensityScaleHeight));

		density += localDensity * stepSize;
	}

	opticalDepthCP = density;

	return true;
}

//----- Input
// position			视线采样点P
// lightDir			光照方向

//----- Output : 
//dpa
//dcp
bool GetAtmosphereDensityRealtime(float3 position, float3 planetCenter, float3 lightDir, out float2 dpa, out float2 dpc)
{
	float height = length(position - planetCenter) - _PlanetRadius;
	dpa = exp(-height.xx / _DensityScaleHeight.xy);

	bool bOverGround = lightSampleing(position, lightDir, dpc);
	return bOverGround;
}

//----- Input
// localDensity			rho(h)
// densityPA
// densityCP

//----- Output : 
// localInscatterR 
// localInscatterM
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;
}

//----- Input
// cosAngle			散射角

//----- Output : 
// scatterR 
// scatterM
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;
}


//----- Input
// rayStart			视线起点 A
// rayDir			视线方向
// rayLength		AB 长度
// planetCenter		地球中心坐标
// distanceScale	世界坐标的尺寸
// lightdir			太阳光方向
// sampleCount		AB 采样次数

//----- Output : 
// extinction       T(PA)
// inscattering:	Inscatering
float4 IntegrateInscatteringRealtime(float3 rayStart, float3 rayDir, float rayLength, float3 planetCenter, float distanceScale, float3 lightDir, float sampleCount, out float4 extinction)
{
	float3 step = rayDir * (rayLength / sampleCount);
	float stepSize = length(step) * distanceScale;

	float2 densityPA = 0;
	float3 scatterR = 0;
	float3 scatterM = 0;

	float2 localDensity;
	float2 densityCP;

	float2 prevLocalDensity;
	float3 prevLocalInscatterR, prevLocalInscatterM;
	GetAtmosphereDensityRealtime(rayStart, planetCenter, lightDir, prevLocalDensity, densityCP);

	ComputeLocalInscattering(prevLocalDensity, densityCP, densityPA, prevLocalInscatterR, prevLocalInscatterM);

	// P - current integration point
	// A - camera position
	// C - top of the atmosphere
	[loop]
	for (float s = 1.0; s < sampleCount; s += 1)
	{
		float3 p = rayStart + step * s;

         // 下面的densityCP和densityPA，其实是指光学厚度
		GetAtmosphereDensityRealtime(p, planetCenter, lightDir, localDensity, densityCP);
		densityPA += (localDensity + prevLocalDensity) * (stepSize / 2.0);
		float3 localInscatterR, localInscatterM;
		ComputeLocalInscattering(localDensity, densityCP, densityPA, localInscatterR, localInscatterM);

		scatterR += (localInscatterR + prevLocalInscatterR) * (stepSize / 2.0);
		scatterM += (localInscatterM + prevLocalInscatterM) * (stepSize / 2.0);

		prevLocalInscatterR = localInscatterR;
		prevLocalInscatterM = localInscatterM;

		prevLocalDensity = localDensity;
	}

	float3 m = scatterM;
	// phase function
	ApplyPhaseFunction(scatterR, scatterM, dot(rayDir, -lightDir.xyz));
	float3 lightInscatter = (scatterR * _ScatteringR  + scatterM * _ScatteringM) * _IncomingLight.xyz;
	float3 lightExtinction = exp(-(densityCP.x * _ExtinctionR + densityCP.y * _ExtinctionM));

	extinction = float4(lightExtinction, 0);
	return float4(lightInscatter, 1);
}



			
float4 frag(v2f i) : SV_Target
{
	float deviceZ = SAMPLE_DEPTH_TEXTURE(_CameraDepthTexture, i.uv);

	float3 positionWorldSpace = GetWorldSpacePosition(i.uv);

	float3 rayStart = _WorldSpaceCameraPos;
	float3 rayDir = positionWorldSpace - _WorldSpaceCameraPos;
	float rayLength = length(rayDir);
	rayDir /= rayLength;

	if (deviceZ < 0.000001)
	{
		rayLength = 1e20;
	}


	float3 planetCenter = float3(0, -_PlanetRadius, 0);
	float2 intersection = RaySphereIntersection(rayStart, rayDir, planetCenter, _PlanetRadius + _AtmosphereHeight);

	rayLength = min(intersection.y, rayLength);
	float4 extinction;
	if (deviceZ < 0.000001)
	{
		float4 inscattering = IntegrateInscatteringRealtime(rayStart, rayDir, rayLength, planetCenter, 1, _LightDir, 16, extinction);
		return inscattering;
	}
	else
		return 0;
		
}

其中用到的beta在三个颜色上的系数是

在公式中有两处都用到了beta，但在代码中并不完全一致，前面的beta在shader中用_ScatteringR和_ScatteringM表示，后面的beta在shader中用_ExtinctionR和_ExtinctionM表示，是为了加些控制参数进去