全局光照算法：reflective shadow maps

1. 技术理解

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

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

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

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

1.1 推导

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

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

忽略可见性，我们的光照积分应该时如下形式：
Lo(p,wo)=∫ΩpacthLi(p,wi)V(p,wi)fr1(p,wi,wo)cos⁡θpdwi=∫ApacthLi(p,wi)fr1(p,wi,wo)cos⁡θpcos⁡θq∣∣q−p∣∣2dAL_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)=∫ΩpacthLi(p,wi)V(p,wi)fr1(p,wi,wo)cosθpdwi=∫ApacthLi(p,wi)fr1(p,wi,wo)∣∣q−p∣∣2cosθpcosθqdA

参考上图，作者假设光源是无限小的，而且dw=dAcos⁡θq∣∣x/−x∣∣2dw=\frac{dA\cos{\theta_q}}{||x^/-x||^2}dw=∣∣x/−x∣∣2dAcosθq，我们可以做出如下推导：法线为n的表面点p因像素光q而产生的辐照度为：
dE(p)=Li(p,wi)cos⁡θpdwiE(p)=∫Licos⁡θpdwi=∫ApatchLicos⁡θpcos⁡θq∣∣q−p∣∣2dAdE(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θpdwiE(p)=∫Licosθpdwi=∫ApatchLi∣∣q−p∣∣2cosθpcosθqdA
又因为对于一个diffuse patch来说，所有方向的出射光都是相同的——可以有Li=fr⋅ΦdAL_i=f_r \cdot \frac{\Phi}{dA}Li=fr⋅dAΦ（这个公式如何理解？），所以：
E(p)=cos⁡θpcos⁡θq∣∣q−p∣∣2ΦpE(p)=\frac{\cos{\theta_p} \cos{\theta_q}}{||q-p||^2}\Phi_p E(p)=∣∣q−p∣∣2cosθpcosθqΦp

实际上：Φp=Φlight⋅frq\Phi_p=\Phi_{light}\cdot f_{r_q}Φp=Φlight⋅frq。最终推导为：

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

2. Pass 1：Generation

2.1 Data

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

如果说，平行光的通量是Φ\PhiΦ，那么照明这个像素块之后，出射辐射率是：Lo=fr⋅ΦdAL_o=f_r\cdot \frac{\Phi}{dA}Lo=fr⋅dAΦ。而这个像素块此时的辐射通量是：
Φp=∫Area∫Ω(Lo)dAdwo=∫Ω(fr⋅Φ)dwi\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Φ和frf_rfr都是常量。而fr=ρ/πf_r=\rho/\pifr=ρ/π和∫(1)dwi=π\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=ρ⋅Φ⋅dot(xL−xp,np)/(∣∣xL−xp∣∣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像素灯光的影响下，导致的辐照度Efloat3 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×512512\times 512512×512），所以上述sum计算是非常昂贵的，在实时情况下不实用。相反，作者必须将总和减少到有限的光源数量，例如400个。作者使用重要性驱动的方法来做到这一点，试图将采样集中到相关像素灯上。

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

首先，作者将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是均匀分布的随机数，作者选择位置：
然后，必须通过用ξ12ξ^2_1ξ12对样本进行加权，来补偿不同的采样密度（以及最后的归一化）。下图显示了一个采样模式的例子。
实际实现过程中，我们在CPU端通过低差异序列，生成随机数(ξ1,ξ2)(\xi_1,\xi_2)(ξ1,ξ2)。s和t不用管，就是GPU端像素的UV坐标，我们只需要计算：rmaxξ1sin⁡(2πξ2)r_{max}\xi_1\sin{(2\pi\xi_2)}rmaxξ1sin(2πξ2)、rmaxξ1cos⁡(2πξ2)r_{max}\xi_1\cos{(2\pi\xi_2)}rmaxξ1cos(2πξ2)、ξ12\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

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

Lo(p,wo)=∫ApacthLi(p,wi)fr1(p,wi,wo)cos⁡θpcos⁡θq∣∣q−p∣∣2dAL_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)=∫ApacthLi(p,wi)fr1(p,wi,wo)∣∣q−p∣∣2cosθpcosθqdA

E(p)=∫Licos⁡θpdwi=∫ApatchLicos⁡θpcos⁡θq∣∣q−p∣∣2dAE(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)=∫Licosθpdwi=∫ApatchLi∣∣q−p∣∣2cosθpcosθqdA

而frf_rfr又和面积没有关系，所以：
Lo(p,wo)=E(p)∗fr1L_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