Move to totally PBR

Date: 2018-11-05 13:15
Categories: 渲染

草稿，未完…

Validating

首先要说的是，做基于物理渲染（着色）Validation非常重要，因为稍有不同的代码和数据可能在某些情形下不会产生十分引人注目的走样甚至有可能看起来还挺正确，但这个看起来还满正确的结果却有可能是完全错误的，他在其他某些情形下可能完全不work.
In our approach we tried to measure and verify actual data, such as light intensity and falloff, sky brightness and camera effects. But all these steps are time consuming and not always easy to set up.
测量brdfMERL数据库包含了至少100种材质的brdf数据，可以作为Ref来验证Brdf模型。
Dice将自己的模型渲染结果与mitsuba的离线渲染结果对比，他们写了一个插件来导出引擎内的灯光和材质参数配置一遍渲染等价的离线结果。同时因为每次为了快速迭代他们还实现了一个引擎内的暴力采样渲染器，一遍渲染引擎内reference。
【此部分工作尚未开始】我准备选择的方案是，将离线渲染器与实时渲染集成到一起，首先是离线渲染器与Mitsuba进行validation，然后使用离线渲染器与实时渲染器做validation。
注意：作为Reference的必须是没有近似的方法，方程和数据，比如菲涅尔不能使用Schilick近似之类的，如果使用的Reference本身就是近似方法，那么结果将不一定可信。但总而言是，唯一可以被完全信任的Reference只有真实世界。

理解漫反射和高光反射

首先要说的是，漫反射是由于光线透过到材质内部多次反射后出射形成的，也就是说漫反射考虑次表面的作用，是相当复杂的过程，但在数学建模中由于其相对均匀的散射分布一般使用比较简单的模型；而高光模型相对来说比较简单，可以比较精确地建模（参考理解微表面）。

材质与着色

能量守恒

无论你选择什么 f d f_d fd​(Diffuse term)和 f r f_r fr​(Specular term)，都必须保证能量守恒，也就是要hemispherical-directional reflectance<=1：

ρ h d = ∫ Ω ( f r + f d ) m a x ( n ⋅ l , 0 ) d l ≤ 1 \rho_{hd}=\int_{\Omega}(f_r+f_d)max(n·l,0)dl \le 1 ρhd​=∫Ω​(fr​+fd​)max(n⋅l,0)dl≤1

选择lambertian+微表面就要引入能量守恒系数 k s k_s ks​和 k d k_d kd​:

L o ( p , ω o ) = ∫ Ω ( k d c π + k s D F G 4 ( ω o ⋅ n ) ( ω i ⋅ n ) ) L i ( p , ω i ) n ⋅ ω i d ω i L_o(p,\omega_o) = \int\limits_{\Omega} (k_d\frac{c}{\pi} + k_s\frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p,\omega_i) n \cdot \omega_i d\omega_i Lo​(p,ωo​)=Ω∫​(kd​πc​+ks​4(ωo​⋅n)(ωi​⋅n)DFG​)Li​(p,ωi​)n⋅ωi​dωi​

这里 k s k_s ks​直接选择为菲涅尔系数Fr——菲涅尔可以代表specular反射了多少光，而1-Fr近似代表1-Fr的光线透射后与次表面作用产生了漫反射，而 k d = 1 − k s k_d=1-k_s kd​=1−ks​，上述方程直接变为：

L o ( p , ω o ) = ∫ Ω ( ( 1 − F ) c π + D F G 4 ( ω o ⋅ n ) ( ω i ⋅ n ) ) L i ( p , ω i ) n ⋅ ω i d ω i L_o(p,\omega_o) = \int\limits_{\Omega} ((1-F)\frac{c}{\pi} + \frac{DFG}{4(\omega_o \cdot n)(\omega_i \cdot n)}) L_i(p,\omega_i) n \cdot \omega_i d\omega_i Lo​(p,ωo​)=Ω∫​((1−F)πc​+4(ωo​⋅n)(ωi​⋅n)DFG​)Li​(p,ωi​)n⋅ωi​dωi​

在DICE的Notes中，他们使用了disney的漫反射brdf，高光依然是微表面，但是Disney的漫反射能量不守恒（为了艺术家创作目的刻意的），DICE于是配合他们使用的高光brdf专门对这个漫反射brdf做了renormalization，以保证总体能量守恒。

BRDF选择

首先必须要注意，漫反射看起来使用的模型很简单，但实际上却恰恰是因为漫反射建模远比高光复杂得多，有时候需要考虑次表面散射。

而对于金属，金属不会折射光线它只吸收光线，所以一般他的介电系数是复数，假设在二氧化硅中参入铝，那么复数的实部表示折射，而虚部表示金属对光的吸收，所以对于金属，考虑为仅仅反射光线。

也因为纯金属吸收所有的折射光，所以他没有漫反射.

漫反射

Lambert:
f r , d = k d c d / π f_{r,d} = k_{d}c_{d}/\pi fr,d​=kd​cd​/π

Disney:
f d = b a s e C o l o r π ( 1 + ( F D 90 − 1 ) ( 1 − c o s θ l ) 5 ) ( 1 + ( F D 90 − 1 ) ( 1 − c o s θ v ) 5 ) f_{d}=\frac{baseColor}{\pi}(1+(F_{D90}-1)(1-cos\theta_{l})^5)(1+(F_{D90}-1)(1-cos\theta_{v})^5) fd​=πbaseColor​(1+(FD90​−1)(1−cosθl​)5)(1+(FD90​−1)(1−cosθv​)5)
where
F D 90 = 0.5 + 2 ∗ r o u g h n e s s ∗ c o s 2 θ d F_{D90} = 0.5+2*roughness*cos^2\theta_{d} FD90​=0.5+2∗roughness∗cos2θd​

为了简单起见这里暂时选择了Lambert。

高光

法线分布GGX：
N D F G G X T R ( n , h , α ) = α 2 π ( ( n ⋅ h ) 2 ( α 2 − 1 ) + 1 ) 2 NDF_{GGX TR}(n, h, \alpha) = \frac{\alpha^2}{\pi((n \cdot h)^2 (\alpha^2 - 1) + 1)^2} NDFGGXTR​(n,h,α)=π((n⋅h)2(α2−1)+1)2α2​

几何Term（G的意义）Smith’s method with Schlick-GGX：
G S c h l i c k G G X ( n , v , k ) = n ⋅ v ( n ⋅ v ) ( 1 − k ) + k G_{SchlickGGX}(n, v, k) = \frac{n \cdot v} {(n \cdot v)(1 - k) + k } GSchlickGGX​(n,v,k)=(n⋅v)(1−k)+kn⋅v​
For direct lighting:
k d i r e c t = ( α + 1 ) 2 8 k_{direct} = \frac{(\alpha + 1)^2}{8} kdirect​=8(α+1)2​
For IBL lighting:
k I B L = α 2 2 k_{IBL} = \frac{\alpha^2}{2} kIBL​=2α2​

注意：此处的alpha取决于粗糙度，但是要看引擎如何转化粗糙度。一般对于DG，roughness平方是一个相对真实的选择。

根据Smith提出的方法，可以试用

G ( n , v , l , k ) = G s u b ( n , v , k ) G s u b ( n , l , k ) G(n, v, l, k) = G_{sub}(n, v, k) G_{sub}(n, l, k) G(n,v,l,k)=Gsub​(n,v,k)Gsub​(n,l,k)

来作为G.

菲涅尔Schlick approximation：
F S c h l i c k ( n , v , F 0 ) = F 0 + ( 1 − F 0 ) ( 1 − ( n ⋅ v ) ) 5 F_{Schlick}(n, v, F_0) = F_0 + (1 - F_0) ( 1 - (n \cdot v))^5 FSchlick​(n,v,F0​)=F0​+(1−F0​)(1−(n⋅v))5
where
F 0 F_0 F0​是垂直看向表面时的reflectivity 。

注意：Fresnel-Schlick approximation仅仅用于绝缘体，此时 F 0 F_0 F0​代表表面的reflectivity ，使用折射率计算 F 0 = ( n 1 − n 2 n 1 + n 2 ) 2 F_0=(\frac{n_1-n_2}{n_1+n_2})^2 F0​=(n1​+n2​n1​−n2​​)2，通常对于空气近似有 n 1 = 1 n_1=1 n1​=1。对于导体为了方便，通常预计算 F 0 F_0 F0​，然后依然可以使用这个公式。

常见的 F 0 F_0 F0​:

|Material|F0 (Linear)|F0 (sRGB)|
|Water	|(0.02, 0.02, 0.02)|	 (0.15, 0.15, 0.15)  |
|Plastic / Glass (Low)|	(0.03, 0.03, 0.03)|	(0.21, 0.21, 0.21)	|
|Plastic High	|(0.05, 0.05, 0.05)	|(0.24, 0.24, 0.24)	|
|Glass (high) / Ruby|	(0.08, 0.08, 0.08)|	(0.31, 0.31, 0.31)|	
|Diamond	|(0.17, 0.17, 0.17)	|(0.45, 0.45, 0.45)	|
|Iron	|(0.56, 0.57, 0.58)	|(0.77, 0.78, 0.78)	|
|Copper|	(0.95, 0.64, 0.54)|	(0.98, 0.82, 0.76)|	
|Gold	|(1.00, 0.71, 0.29)	|(1.00, 0.86, 0.57)	|
|Aluminium	|(0.91, 0.92, 0.92)|	(0.96, 0.96, 0.97)|	
|Silver|	(0.95, 0.93, 0.88)|	(0.98, 0.97, 0.95)|	

注意绝缘体的 F 0 F_0 F0​通常都是白色的(也就是说是一个标量)，而且通常都很低，不会大于0.17，而导体一般都在0.5到1.0之间变化，而且是逐波长变化的，所以直接观察金属表面是有色的。
金属和绝缘体的性质差异使得我们使用一个叫做金属工作流的工作流程。我们将使用一个额外的参数金属性来描述表面是金属还是非金属。
很显然理论上，表面要么是金属要么是非金属，如果是金属他将不会呈现出漫反射颜色（金属吸收所有折射光线，他的粗糙只由高光部分反映出来），而高光则完全贴近albedo，所以可以直接使用表面纹理颜色来作为他的F0。
而0.04对于大多数绝缘体来说都是一个基础的F0，所以我们假设非金属F0总是0.04.
但实际上往往把金属性线性分布到0到1之间，以此可以表现出金属表面的尘埃，刮痕，以及相互掺杂等效果。
所以：
F0 = lerp(0.04,albedo,metallic)
这表明：随着金属性的增强，垂直高光颜色将不断从0.04贴近albedo.

Albedo纹理对于绝缘体应该保存表面颜色，对于金属则应该保存 F 0 F_0 F0​(base reflectivity).Albedo纹理和漫反射纹理不同的地方就在于金属和漫反射纹理有时候会画一些AO之类的光照信息在上面。在PBR工作流中，AO将被单独放到一张纹理中。

但是要注意：金属表面因为没有漫反射，所以要讲金属性为1的材质漫反射term完全归零：

ks=Fr
kd=(1-ks)(1-metallic)

灯光

IBL

At normal incident direction, the shape of a BRDF is isotropic. At grazing angles
the shape of a BRDF is anisotropic. Removing the view dependency for pre-integrating Equation 46
would make the assumption that the BRDF shape is isotropic at all view angles. This leads to key
visual differences, preventing stretched reflections. This approximation can be quite noticeable on flat
surfaces as shown on Figure 55 but less on curvy surfaces .
V=R=N

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7phQJJOG-1605237060413)(https://github.com/wubugui/FXXKTracer/raw/master/pic/newblog/atm/dome1.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-uotcwwLu-1605237060415)(https://github.com/wubugui/FXXKTracer/raw/master/pic/newblog/atm/demo2.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gT3pLB3o-1605237060416)(https://github.com/wubugui/FXXKTracer/raw/master/pic/newblog/atm/demo3.png)]

实现

Split Sum导出

实现摘要

注意正确地做gamma矫正，再线性空间中做光照计算，输出的时候也要正确矫正。
一种简单的tonemapping方式和gamma矫正：

 color = color / (color+1.f);
 color = pow(color,1/2.2f);

预过滤漫反射envmap

基本的公式：

L o ( p , ω o ) = k d c π ∫ Ω L i ( p , ω i ) n ω ˙ i d ω i L_o(p,\omega_o)=k_d\frac{c}{\pi}\int_{\Omega}Li(p,\omega_i)n \dot \omega_i d\omega_i Lo​(p,ωo​)=kd​πc​∫Ω​Li(p,ωi​)nω˙i​dωi​

一种方法是使用黎曼和直接展开积分，对于cubemap的每一个纹素都计算一次。这种方法使用很多采样来离散化积分，缺点是会产生很多banding。

L o ( p , θ o , ϕ o ) = c π ∫ Φ ∫ Θ L i ( p , θ i , ϕ i ) c o s ( θ i ) s i n ( θ i ) d θ i d ϕ i L_o(p,\theta_o,\phi_o) = \frac{c}{ \pi } \int\limits_{\Phi} \int\limits_{\Theta} L_i(p,\theta_i,\phi_i) cos(\theta_i) sin(\theta_i) d\theta_i d\phi_i Lo​(p,θo​,ϕo​)=πc​Φ∫​Θ∫​Li​(p,θi​,ϕi​)cos(θi​)sin(θi​)dθi​dϕi​

对这个式子的两个积分变量进行均匀采样（其实会发现此时的雷曼和和MC积分形式是一致的），得到：

L o ( p , θ o , ϕ o ) = c π 2 π N 1 π 2 N 2 ∑ N 1 ∑ N 2 L i ( p , θ i , ϕ i ) c o s ( θ i ) s i n ( θ i ) L_o(p,\theta_o,\phi_o) = \frac{c}{ \pi }\frac{2 \pi}{ N_1 } \frac{\pi}{ 2 N_2 } \sum_{}^{N_1} \sum_{}^{N_2} L_i(p,\theta_i,\phi_i) cos(\theta_i) sin(\theta_i) Lo​(p,θo​,ϕo​)=πc​N1​2π​2N2​π​∑N1​​∑N2​​Li​(p,θi​,ϕi​)cos(θi​)sin(θi​)
L o ( p , θ o , ϕ o ) = π c N 1 N 2 ∑ N 1 ∑ N 2 L i ( p , θ i , ϕ i ) c o s ( θ i ) s i n ( θ i ) L_o(p,\theta_o,\phi_o) = \frac{\pi c}{ N_1 N_2 } \sum_{}^{N_1} \sum_{}^{N_2} L_i(p,\theta_i,\phi_i) cos(\theta_i) sin(\theta_i) Lo​(p,θo​,ϕo​)=N1​N2​πc​∑N1​​∑N2​​Li​(p,θi​,ϕi​)cos(θi​)sin(θi​)

但我这里采用的并非此法，我这里因为使用的离线工具在CPU计算积分，所以我们将更加容易地利用MIS（甚至可以使用更多更好的方法）来减少方差加快收敛，这里对于计算这积分，漫反射因为brdf本身是均匀的，只有envmap不是，所以直接重要性采样envmap即可（注意IS的时候要注意只要有值无论多小都要求pdf不可以为0！）

下面是使用MIS采样brdf和envmap的代码，去掉brdf采样部分就可以得到envmap only IS.

float3 DiffuseIBLWeightMis(const float3& normal, const YRInfiniteAreaLight& light,YRIndependentSampler& sp)
{
	//sample cosine hemsphere pdf just cos\theta/PI
	//so the estimator is : L_o = c/N\sum L_i
	//we pre-integrate : W_o = 1/N\sum L_i
	const uint sample_count_diffuse = SN;
	RBColorf weight = 0.0;
	for (uint i = 0; i < sample_count_diffuse; ++i)
	{
		float2 rr = sp.get2d();
		float2 rr1 = sp.get2d();
		RBVector3 light_wi;
		float light_pdf;
		float brdf_pdf;
		RBColorf light_val = light.sample_l(RBVector3::zero_vector, 0.001f, 0, &light_wi, &light_pdf, nullptr, rr1);
		if (!light_val.is_black())
		{
			brdf_pdf = RBVector3::dot_product(light_wi, normal) / PI;
			if (brdf_pdf > 0)
			{
				float weight_light =  power_mis(light_pdf, brdf_pdf);
				weight += light_val*RBVector3::dot_product(light_wi, normal) / light_pdf*weight_light;
			}
		}
		
		float3 L = hemisphereSample_cos(rr.x, rr.y, normal);
		float NoL = RBMath::clamp(float3::dot_product(normal, L), 0.f, 1.f);
		light_val = light.l(L);
		light_pdf = light.pdf(RBVector3(0.f), L);
		brdf_pdf = NoL / PI;
		float weight_brdf = power_mis(brdf_pdf, light_pdf);
		weight += light_val*PI*weight_brdf;
		
	}
	weight /= sample_count_diffuse;
	return weight;
}

Environment Probe

Sky Lighting

首先，大气要作为光源照亮场景，如果是离线渲染，当然是毫不犹豫的积分。但在实时中是不可取的。为了解决这个问题，当前有一个初步想法是，将着色分为两个部分，漫反射和高光，漫反射使用irradiance map，高光用light probe。但是这种清醒下我们有个问题，好处

https://software.intel.com/sites/default/files/managed/b9/1d/gdc2013-lightscattering-final.pdf

回顾

我之前参考【http://www.scratchapixel.com/lessons/procedural-generation-vritual-worlds/simulating-sky】做了一次大气，大体的记录如下：

大气强度衰减公式：
d e n s i t y ( h ) = d e n s i t y ( 0 ) e − h H density(h)=density(0) e^{-\frac{h}{H}} density(h)=density(0)e−Hh​

大气散射粒子分为气体原子散射和尘埃闪射。前者叫做Rayleigh scattering后者叫做Mie scattering。前者负责呈现天空的蓝色和日落日出时的红色，后者负责灰白色（大气污染那种颜色）。

注意单位差距，粒子的单位在毫米微米级，而星球的半径则在千米级。

Rayleigh Scattering

在体渲染中，衰减系数决定散射量，相位函数决定散射方向。其中衰减系数为散射系数和吸收系数之和。
散射系数：
β R s ( h , λ ) = 8 π 3 ( n 2 − 1 ) 2 3 N λ 4 e − h H R \beta_R^s (h, \lambda) = \frac{8\pi^3(n^2-1)^2}{3N\lambda^4} e^{-\frac{h}{H_R}} βRs​(h,λ)=3Nλ48π3(n2−1)2​e−HR​h​
h h h是海拔高度， N N N是海平面气体原子密度， λ \lambda λ是波长， n n n是气体折射率， H R H_R HR​是scale height。我们使用 H R = 8 k m H_R=8km HR​=8km.
注意：我们完全可以使用一些测量值来计算N，n等等参数，但是我们这里将会使用预计算的Rayleigh散射系数（ 33.1 e − 6 m − 1 33.1 \mathrm{e^{-6}} \mathrm{m^{-1}} 33.1e−6m−1, 13.5 e − 6 m − 1 13.5 \mathrm{e^{-6}} \mathrm{m^{-1}} 13.5e−6m−1, 5.8 e − 6 m − 1 5.8 \mathrm{e^{-6}} \mathrm{m^{-1}} 5.8e−6m−1 for wavelengths 440, 550 and 680 respectively.）
在大气渲染中，吸收可以忽略不计所以：
β R e = β R s \beta_R^e = \beta_R^s βRe​=βRs​
相位函数为：
P R ( θ ) = 3 16 π ( 1 + θ 2 ) P_R(\theta)=\frac{3}{16\pi}(1+\theta^2) PR​(θ)=16π3​(1+θ2)
其中 θ \theta θ是入射光线与视线的夹角。

Mie Scattering

散射系数：
β M s ( h , λ ) = β M s ( 0 , λ ) e − h H M \beta_M^s(h,\lambda)=\beta_M^s(0,\lambda) e^{-\frac{h}{H_M}} βMs​(h,λ)=βMs​(0,λ)e−HM​h​
H M H_M HM​是scale height 通常为1.2km.我们使用预计算值 β S = 210 e − 5 m − 1 \beta_S = 210 \mathrm{e^{-5}} \mathrm{m^{-1}} βS​=210e−5m−1。衰减系数是散射系数的约1.1倍。
相位函数：
P M ( μ ) = 3 8 π ( 1 − g 2 ) ( 1 + μ 2 ) ( 2 + g 2 ) ( 1 + g 2 − 2 g μ ) 3 2 P_M(\mu)=\frac{3}{8\pi}\frac{(1-g^2)(1+\mu^2)}{(2+g^2)(1+g^2-2g\mu)^{\frac{3}{2}}} PM​(μ)=8π3​(2+g2)(1+g2−2gμ)23​(1−g2)(1+μ2)​
g g g用来决定各向异性，前向为[0,1],后向为[-1,0].我们取0.76，0代表各向同。

Ray March

体渲染方程：
L ( P c , P a ) = ∫ P c P a T ( P c , X ) L s u n ( X ) d s L(P_c, P_a)=\int_{P_c}^{P_a} T(P_c,X)L_{sun} (X)ds L(Pc​,Pa​)=∫Pc​Pa​​T(Pc​,X)Lsun​(X)ds
T ( P a , P b ) = L a L b = e x p ( − ∑ P a P b β e ( h ) d s ) T(P_a,P_b)=\frac{L_a}{L_b}=exp(-\sum_{P_a}^{P_b} \beta_e(h)ds) T(Pa​,Pb​)=Lb​La​​=exp(−Pa​∑Pb​​βe​(h)ds)

β s ( h ) = β s ( 0 ) e x p ( − h H ) \beta_s(h) = \beta_s(0)exp(-\frac {h}{H}) βs​(h)=βs​(0)exp(−Hh​)
β e ( h ) = β s ( h ) + β a ( h ) = β s ( h ) + 0 = β s ( h ) \beta_e(h)=\beta_s(h)+\beta_a(h)=\beta_s(h)+0=\beta_s(h) βe​(h)=βs​(h)+βa​(h)=βs​(h)+0=βs​(h)
T ( P a , P b ) = e x p ( − β e ( 0 ) ∑ P a P b e x p ( − h H ) d s ) T(P_a, P_b)=exp(-\beta_e(0) \sum_{P_a}^{P_b} exp(-\frac{h}{H})ds) T(Pa​,Pb​)=exp(−βe​(0)Pa​∑Pb​​exp(−Hh​)ds)

L s u n ( X ) = S u n I n t e n s i t y T ( X , P s ) P ( V , L ) ∗ β S ( h ) L_{sun}(X)=Sun Intensity T(X,P_s) P(V,L)*\beta_S(h) Lsun​(X)=SunIntensityT(X,Ps​)P(V,L)∗βS​(h)

最终：

S k y   C o l o r ( P c , P a ) = ∫ P c P a T ( P c , X ) S u n   I n t e n s i t y P ( V , L ) T ( X , P s ) β s ( h ) d s Sky \: Color(P_c, P_a) = \int_{P_c}^{P_a} T(P_c, X) Sun \: Intensity P(V,L) T(X,P_s) \beta_s(h)ds SkyColor(Pc​,Pa​)=∫Pc​Pa​​T(Pc​,X)SunIntensityP(V,L)T(X,Ps​)βs​(h)ds

对于方程中的 S u n I n t e n s i t y SunIntensity SunIntensity和 P ( V , L ) P(V,L) P(V,L)都是常量可以提出来，两个 T T T合在一起：

T ( P c , X ) T ( X , P s ) = e − β e 0 e − β e 1 = e − ( β e 0 + β e 1 ) T(P_c,X) T(X,P_s)=e^{ -\beta_{e0} } e^{ -\beta_{e1} }=e^{ -(\beta_{e0} + \beta_{e1}) } T(Pc​,X)T(X,Ps​)=e−βe0​e−βe1​=e−(βe0​+βe1​)

最后我们将分别计算Rayleigh和Mie各一次加到一起。

GPU初步实现

我是用了一个全屏pass渲染的，基本就是一个后处理。但是渲染在frame的最前面。PS传入摄像机参数生成每个像素的View方向。
需要注意的地方是数据精度问题，宇宙级的半径和原子半径差距太大，再加上积分操作，很容易造成精度丢失。我最初的时候把单位统一为cm，最后输出了一片黑，大概算了一下，发现光学深度这个数据非常的小，因为拿取做指数衰减的项非常大，宇宙级的。
仔细观察不难发现那些很小的部分都是线性乘上去的，前面的计算和mm级的数据没有交集，所以前面的计算可以使用千米为单位，在某些预估可能会比较大的地方转化成m。在整个计算过程中根据计算项的意义不断地调整单位，最后可以把单位全部统一到m。
这个统一到m的量就可以使用一个scale来统一缩放了。
还有一点是，在计算view方向的积分的时候，最好是把有效的采样范围缩到最小，这样可以提高采样率，否则当摄像机在大气外时，地球半径级别的距离会导致采样率非常低。

struct VS_Output
{  
    float4 Pos : SV_POSITION;              
    float2 Tex : TEXCOORD0;
};

VS_Output main_vs(uint id : SV_VertexID)
{
    VS_Output Output;
    Output.Tex = float2((id << 1) & 2, id & 2);
    Output.Pos = float4(Output.Tex * float2(2,-2) + float2(-1,1), 0, 1);
    return Output;
}

#define M_PI 3.141592653

int isc_sphere(float3 ro,float3 rd,float3 pos,float r,inout float t1,inout float t2)
{
	float t;
	float3 temp = ro - pos;
	float a = dot(rd, rd);
	float b = 2.f*dot(temp, rd);
	float c = dot(temp, temp) - r*r;
	float ds = b*b - 4 * a*c;

		
	int isc_n = 0;
	if (ds < 0.f)
	{
		return isc_n;
	}
	else
	{
		float e = sqrt(ds);
		float de = 2.f * a;
		t = (-b - e) / de;

		if (t > 0.00001)
		{
			t1 = t;
			isc_n++;
			//return isc_n;
		}

		t = (-b + e) / de;

		if (t>0.00001)
		{
			t2 = t;
			isc_n++;
			//return isc_n;
		}
	}
	return isc_n;
}

cbuffer CB
{
	//km
	float3 sphere_pos_v;
	float view_dist;
	float3 sun_dir;
	float sacle_factor;
	float tan_fov_x;
	float tan_fov_y;
	int n_sample;
	int n_sample_light;
	matrix mv_mat;
};

float4 main_ps(VS_Output pin): SV_TARGET0
{
	sphere_pos_v = mul(sphere_pos_v*1000,mv_mat);
	sphere_pos_v/=1000;
	sun_dir = mul(sun_dir,mv_mat);
	
	//km
	float er = 6360;
	float ar = 6420;
	//m
	float hr = 7994;
	float hm = 1200;
	
	//mm
	float3 betar= float3(0.0058,0.0135,0.0331);
	float3 betam = float3(0.021,0.021,0.021);
	
	float tx = 2*pin.Tex.x - 1;
	float ty = 2*(1-pin.Tex.y) - 1;
	float vx = tx*view_dist* tan_fov_x;
	float vy = ty*view_dist*tan_fov_y;
	
	//cm
	float3 view_dir = float3(vx,vy,view_dist);
	view_dir = normalize(view_dir);
	
	float t1 = -1;
	float t2 = -1;
	float t11 = -1;
	float t12 = -1;
	int ik = isc_sphere(float3(0,0,0),view_dir,sphere_pos_v,ar,t1,t2);
	//检测内球交点，提升采样率
	int ic = isc_sphere(float3(0,0,0),view_dir,sphere_pos_v,er,t11,t12);
	float3 pa,pb;
	if(ik==0)
		discard;
		//return float4(1,0,0,1);
	if(ik==1)
	{
		pa = float3(0,0,0);
		pb = view_dir*t2;
		//return float4(0,0,1,1);
	}
	else if(ik==2)
	{
		
		if(ic<=1)
		{
			pa = view_dir*t1;
			pb = view_dir*t2;
		}
		else if(ic==2 )
		{
			pa = view_dir*t1;
			pb = view_dir*t11;
		}
		//return float4(0,1,0,1);
	}
	//km
	float seg_len = length(pb-pa)/n_sample;
	float3 sumr=0,summ=0;
	float optiacal_depth_r=0,optical_depth_m=0;
	float theta = dot(view_dir,sun_dir);
	float phaser = 3.f / (16.f * M_PI) * (1 + theta * theta);
	float g = 0.76;
	float phasem = 3.f / (8.f * M_PI) * ((1.f - g * g) * (1.f + theta * theta)) / ((2.f + g * g) * pow(1.f + g * g - 2.f * g * theta, 1.5f));
		
	for(uint i=0;i<n_sample;++i)
	{
		//km
		float3 sample_pos = pa + seg_len*(i+0.5)*view_dir;
		float samlpe_height = length(sample_pos - sphere_pos_v) - er;
		//km to m
		samlpe_height*=1000;
		float h_r = exp(-samlpe_height/hr)*seg_len;
		float h_m = exp(-samlpe_height/hm)*seg_len;
		//to km
		optiacal_depth_r += h_r;
		optical_depth_m += h_m;
		float t0Light;
		float t1Light;
		int sn = isc_sphere(sample_pos,sun_dir,sphere_pos_v,ar,t0Light,t1Light);

		//km
		float seg_len_sun = t1Light/n_sample_light;
		float optiacal_depth_sun_r=0,optical_depth_sun_m=0;
		uint j = 0;
		for(j=0;j<n_sample_light;++j)
		{
			float3 sample_pos_sun = sample_pos + seg_len_sun*(j+0.5)*sun_dir;
			float sample_height_sun = length(sample_pos_sun - sphere_pos_v) - er;
			//km to m
			sample_height_sun*=1000;
			//km
			optiacal_depth_sun_r += exp(-sample_height_sun/hr)*seg_len_sun;
			optical_depth_sun_m += exp(-sample_height_sun/hm)*seg_len_sun;			
		}
		{
			//mm*km=m^2
			float3 tau = betar*(optiacal_depth_r+optical_depth_m) + betam*1.1f*(optiacal_depth_sun_r+optical_depth_sun_m);
			float3 attenuation = float3(exp(-tau.x),exp(-tau.y),exp(-tau.z));
			sumr += attenuation*h_r;
			summ += attenuation*h_m;
		}
	}
	
	float3 sky_color = clamp((sumr*betar*phaser + summ*betam*phasem)*sacle_factor,0.0001,100000000);
	return float4(sky_color,1);
}

大致实现的结果是这样的：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VlB6vw7w-1605237060418)(https://github.com/wubugui/FXXKTracer/raw/master/pic/newblog/atm/reviewed.png)]

这个结果一般，效率捉急，而且仅仅是渲染了大气，却并没有成为光源。为了在实时中取得更高的效率和更好的效果，我们必须要做一定的预计算。

To be continued…