全局光照算法：reflective shadow maps

1. 技术理解

RSM的全称是reflective shadow maps，受到Instant Radiosity这个离线技术的启发，其思想和ShadowMap的思想近似。在正式介绍和了解这个技术之前，我需要确定RSM用处何在？我想，《RTR4》中对它的分类很正确——Dynamic Diffuse Global Illumination，这是一个处理动态全局漫反射的技术：

• GI（全局光照）：用于二次及以上bounce造成的间接光。
• Dynamic（动态）：可以实时更新，可作用于动态物体。
• Diffuse（漫反射）：更加细致的说，这个技术考虑的是间接光照中的漫反射部分。

RSM和立即辐射度方法一样，都是在直接光照亮的区域，选择采样点作为发光物（虚拟点光源VPL），来计算间接光照。主要分为两个Pass：

1. Pass 1：从光源的角度，对整个场景进行一次渲染。这个过程在引擎中可以直接和ShadowMap的生成放在一块。不过不同的是，除了存储深度之外（这个感觉就是直接使用ShadowMap的结果），我们还需要存储世界空间位置、法线、辐射通量。
2. Pass 2：在正常的lighting pass中，考虑基于RSM的间接光。

1.1 推导

在RSM中，每个像素都被解释为一个间接照亮场景的像素光。通量 ϕ p \phi_p ϕp​定义了其亮度。

辐射通量 ϕ p \phi_p ϕp​：radiant flux，单位是W（瓦特）。描述的是光源的总体能量——每秒的发出的辐射能量。而辐照度E的单位是 W / m 2 W/m^2 W/m2，辐射强度I的单位是 W / s r W/sr W/sr。我们在光照积分中比较常见的是：辐射率L，单位是 W / ( m 2 s r ) W/(m^2sr) W/(m2sr)
注意：当光源的面积无限小时，基本可以将辐射强度I等同于辐射度L。

忽略可见性，我们的光照积分应该时如下形式：
L o ( p , w o ) = ∫ Ω p a c t h L i ( p , w i ) V ( p , w i ) f r 1 ( p , w i , w o ) cos ⁡ θ p d w i = ∫ A p a c t h L i ( p , w i ) f r 1 ( p , w i , w o ) cos ⁡ θ p cos ⁡ θ q ∣ ∣ q − p ∣ ∣ 2 d A L_o(p,w_o)=\int_{\Omega_{pacth}}{L_i(p,w_i)V(p,w_i)f_{r1}(p,w_i,w_o)\cos{\theta_p}dw_i} \\ =\int _{A_{pacth}}L_i(p,w_i)f_{r1}(p,w_i,w_o)\frac{\cos{\theta_p} \cos{\theta_q}}{||q-p||^2} dA \\ Lo​(p,wo​)=∫Ωpacth​​Li​(p,wi​)V(p,wi​)fr1​(p,wi​,wo​)cosθp​dwi​=∫Apacth​​Li​(p,wi​)fr1​(p,wi​,wo​)∣∣q−p∣∣2cosθp​cosθq​​dA

参考上图，作者假设光源是无限小的，而且 d w = d A cos ⁡ θ q ∣ ∣ x / − x ∣ ∣ 2 dw=\frac{dA\cos{\theta_q}}{||x^/-x||^2} dw=∣∣x/−x∣∣2dAcosθq​​，我们可以做出如下推导：法线为n的表面点p因像素光q而产生的辐照度为：
d E ( p ) = L i ( p , w i ) cos ⁡ θ p d w i E ( p ) = ∫ L i cos ⁡ θ p d w i = ∫ A p a t c h L i cos ⁡ θ p cos ⁡ θ q ∣ ∣ q − p ∣ ∣ 2 d A dE(p)=L_i(p,w_i)\cos{\theta_p}dw_i \\ \\ E(p)=\int L_i\cos{\theta_p}dw_i=\int_{A_{patch}}L_i\frac{\cos{\theta_p} \cos{\theta_q}}{||q-p||^2} dA dE(p)=Li​(p,wi​)cosθp​dwi​E(p)=∫Li​cosθp​dwi​=∫Apatch​​Li​∣∣q−p∣∣2cosθp​cosθq​​dA
又因为对于一个diffuse patch来说，所有方向的出射光都是相同的——可以有 L i = f r ⋅ Φ d A L_i=f_r \cdot \frac{\Phi}{dA} Li​=fr​⋅dAΦ​（这个公式如何理解？），所以：
E ( p ) = cos ⁡ θ p cos ⁡ θ q ∣ ∣ q − p ∣ ∣ 2 Φ p E(p)=\frac{\cos{\theta_p} \cos{\theta_q}}{||q-p||^2}\Phi_p E(p)=∣∣q−p∣∣2cosθp​cosθq​​Φp​

实际上： Φ p = Φ l i g h t ⋅ f r q \Phi_p=\Phi_{light}\cdot f_{r_q} Φp​=Φlight​⋅frq​​。最终推导为：

争议： ∣ ∣ x − x p ∣ ∣ ||x-x_p|| ∣∣x−xp​∣∣的上标是4，还是2，有所争议。
参考：

2. Pass 1：Generation

2.1 Data

由果推因，最后的计算公式需要世界空间位置、法线、辐射通量，我们就存储它们。但现在，依然存在一个问题，我们存取的通量怎么获得？

如果说，平行光的通量是 Φ \Phi Φ，那么照明这个像素块之后，出射辐射率是： L o = f r ⋅ Φ d A L_o=f_r\cdot \frac{\Phi}{dA} Lo​=fr​⋅dAΦ​。而这个像素块此时的辐射通量是：
Φ p = ∫ A r e a ∫ Ω ( L o ) d A d w o = ∫ Ω ( f r ⋅ Φ ) d w i \Phi_p=\int_{Area}\int_{\Omega}(L_o)dAdw_o=\int_{\Omega}(f_r\cdot \Phi) dw_i Φp​=∫Area​∫Ω​(Lo​)dAdwo​=∫Ω​(fr​⋅Φ)dwi​
由于这个patch是漫反射的，所以 Φ \Phi Φ和 f r f_r fr​都是常量。而 f r = ρ / π f_r=\rho/\pi fr​=ρ/π和 ∫ ( 1 ) d w i = π \int(1) dw_i=\pi ∫(1)dwi​=π，最终：
Φ p = ρ ⋅ Φ \Phi_p=\rho\cdot \Phi Φp​=ρ⋅Φ

2.2 实现

所以，最终，我们在第一个Pass中，这样做：

代码来自：https://github.com/AngelMonica126/GraphicAlgorithm/blob/master/001_Reflective%20shadow%20map/RSMBuffer_FS.glsl

void main()
{
	vec3 TexelColor = texture(u_DiffuseTexture, v2f_TexCoords).rgb;
	//TexelColor = pow(TexelColor, vec3(2.2f));
	vec3 VPLFlux = u_LightColor * TexelColor;
	Flux_ = VPLFlux;
	Normal_ = v2f_Normal;
	Position_ = v2f_FragPosInViewSpace;
}

3.3 点光源

之前我们考虑的都是平行光，如果是点光源，我们或许应该在这里考虑一下光线衰减和余弦问题：
Φ p = ρ ⋅ Φ ⋅ d o t ( x L − x p , n p ) / ( ∣ ∣ x L − x p ∣ ∣ 2 ) \Phi_p=\rho\cdot \Phi \cdot dot(x_L-x_p,n_p)/(||x_L-x_p||^2) Φp​=ρ⋅Φ⋅dot(xL​−xp​,np​)/(∣∣xL​−xp​∣∣2)

3 Pass 2：Lighting

3.1 基础实现

主要流程，读取上一个pass存的数据，利用下面的公式，计算间接照明。

float3 indirectIllumination = float3(0, 0, 0);
//最远采样半径
float rMax = rsmRMax;

// rsmSampleCount = hight * width（etc. 512*512）
for (uint i = 0; i < rsmSampleCount; ++i)
{
	// 这里就是随机值
	float2 rnd = rsmSamples[i].xy;
	float2 coords = textureSpacePosition.xy + rMax * rnd;
	// 依次读取位置、法线、通量
	float3 vplPositionWS = g_rsmPositionWsMap.Sample(g_clampedSampler, coords.xy).xyz;
	float3 vplNormalWS = g_rsmNormalWsMap.Sample(g_clampedSampler, coords.xy).xyz;
	float3 flux = g_rsmFluxMap.Sample(g_clampedSampler, coords.xy).xyz;
	// 计算当前像素在此RSM像素灯光的影响下，导致的辐照度E
	float3 result = flux
	* ((max(0, dot(vplNormalWS, P – vplPositionWS))
	* max(0, dot(N, vplPositionWS – P)))
	/ pow(length(P – vplPositionWS), 4));
	indirectIllumination += result;
}
indirectIllumination = result / rsmSampleCount;
return saturate(indirectIllumination * rsmIntensity);

3.2 改进方法

对于一个典型的阴影图来说，像素的数量是很大的（ 512 × 512 512\times 512 512×512），所以上述sum计算是非常昂贵的，在实时情况下不实用。相反，作者必须将总和减少到有限的光源数量，例如400个。作者使用重要性驱动的方法来做到这一点，试图将采样集中到相关像素灯上。

一般来说，可以说x和阴影图中的像素光 x p x_p xp​之间的距离是它们在世界空间中的距离的合理近似值。如果相对于光源的深度值差别很大，世界空间的距离就会大得多，会高估其影响。然而，重要的间接光源总是很接近，这些光源在阴影图中也必须是接近的。
所以作者决定按以下方式获得像素光的样本：

• 首先，作者将x投影到阴影图 ( → ( s , t ) ) (→(s,t)) (→(s,t))中。

• 然后，作者选择 ( s , t ) (s,t) (s,t)周围的像素光，样本密度随着与 ( s , t ) (s,t) (s,t)距离的平方而减少。这可以通过选择相对于 ( s , t ) (s,t) (s,t)的极坐标样本轻松实现，也就是说，如果 ξ 1 ξ_1 ξ1​和 ξ 2 ξ_2 ξ2​是均匀分布的随机数，作者选择位置：
  

• 然后，必须通过用 ξ 1 2 ξ^2_1 ξ12​对样本进行加权，来补偿不同的采样密度（以及最后的归一化）。下图显示了一个采样模式的例子。

• 
  实际实现过程中，我们在CPU端通过低差异序列，生成随机数 ( ξ 1 , ξ 2 ) (\xi_1,\xi_2) (ξ1​,ξ2​)。s和t不用管，就是GPU端像素的UV坐标，我们只需要计算： r m a x ξ 1 sin ⁡ ( 2 π ξ 2 ) r_{max}\xi_1\sin{(2\pi\xi_2)} rmax​ξ1​sin(2πξ2​)、 r m a x ξ 1 cos ⁡ ( 2 π ξ 2 ) r_{max}\xi_1\cos{(2\pi\xi_2)} rmax​ξ1​cos(2πξ2​)、 ξ 1 2 \xi_1^2 ξ12​。将这三个数据存储四维向量数组，作为uniform data传入GPU：

此代码非原创，来自：https://github.com/AngelMonica126/GraphicAlgorithm/blob/master/001_Reflective%20shadow%20map/ShadingWithRSMPass.cpp

std::default_random_engine e;
std::uniform_real_distribution<float> u(0, 1);
for (int i = 0; i < m_VPLNum; ++i)
{
	float xi1 = u(e);
	float xi2 = u(e);
	m_VPLsSampleCoordsAndWeights.push_back({ xi1 * sin(2 * ElayGraphics::PI * xi2), xi1 * cos(2 * ElayGraphics::PI * xi2), xi1 * xi1, 0 });
}

genBuffer(GL_UNIFORM_BUFFER, m_VPLsSampleCoordsAndWeights.size() * 4 * sizeof(GL_FLOAT), m_VPLsSampleCoordsAndWeights.data(), GL_STATIC_DRAW, 1);

然后，在GPU端主要加入的就是这个权重：

for (int i = 0; i < u_NumSamples; i++)
{
    vec3 offset    = texelFetch(s_Samples, ivec2(i, 0), 0).rgb;
    vec2 tex_coord = light_coord.xy + offset.xy * u_SampleRadius + (((offset.xy * u_SampleRadius) / 2.0) * dither_offset);

    vec3 vpl_pos    = texture(s_RSMWorldPos, tex_coord).rgb;
    vec3 vpl_normal = normalize(texture(s_RSMNormals, tex_coord).rgb);
    vec3 vpl_flux   = texture(s_RSMFlux, tex_coord).rgb;

    vec3 result = light_attenuation(vpl_pos) * vpl_flux * ((max(0.0, dot(vpl_normal, (P - vpl_pos))) * max(0.0, dot(N, (vpl_pos - P)))) / pow(length(P - vpl_pos), 4.0));
	
	// 权重
    result *= offset.z * offset.z;
    indirect += result;
}

3.3 关于最后的结果是否要乘上albedo

鄙人认为，最后得到的 E p E_p Ep​，不是当前像素本身产生辐照度（对眼睛生效），而是其他像素灯对此像素产生的辐照度，或者推导上可以看看：

L o ( p , w o ) = ∫ A p a c t h L i ( p , w i ) f r 1 ( p , w i , w o ) cos ⁡ θ p cos ⁡ θ q ∣ ∣ q − p ∣ ∣ 2 d A L_o(p,w_o) =\int _{A_{pacth}}L_i(p,w_i)f_{r1}(p,w_i,w_o)\frac{\cos{\theta_p} \cos{\theta_q}}{||q-p||^2} dA \\ Lo​(p,wo​)=∫Apacth​​Li​(p,wi​)fr1​(p,wi​,wo​)∣∣q−p∣∣2cosθp​cosθq​​dA

E ( p ) = ∫ L i cos ⁡ θ p d w i = ∫ A p a t c h L i cos ⁡ θ p cos ⁡ θ q ∣ ∣ q − p ∣ ∣ 2 d A E(p)=\int L_i\cos{\theta_p}dw_i=\int_{A_{patch}}L_i\frac{\cos{\theta_p} \cos{\theta_q}}{||q-p||^2} dA E(p)=∫Li​cosθp​dwi​=∫Apatch​​Li​∣∣q−p∣∣2cosθp​cosθq​​dA

而 f r f_r fr​又和面积没有关系，所以：
L o ( p , w o ) = E ( p ) ∗ f r 1 L_o(p,w_o)=E(p)*f_{r1} Lo​(p,wo​)=E(p)∗fr1​
也就是说，我认为，最终的间接光照结果，应该是：

indirect = indirect * albedo / PI;
indirect = indirect / VPL_NUM;
light_result = directLight + indirect;

4 其他

ToDo