Shadow Mapping 基本实现

Shadow Mapping 基本原理：

1. 第一个Pass：
   • 额外设置一个摄像机在光源位置（Light Camera，光源摄像机），并且朝光照方向看去。
   • 用一张 Texture（称为 阴影贴图 Shadow Map）来记录 Light Camera 所看到的像素深度（每个像素位置只记录所见最近深度，而不用做别的 shading 计算）来作为遮挡深度。

// shadowVertex.glsl
// ...
void main(void) {
  vNormal = aNormalPosition;
  vTextureCoord = aTextureCoord;
  gl_Position = uLightMVP * vec4(aVertexPosition, 1.0);
}

// shadowFragment.glsl
// ...
void main(){
  gl_FragColor = pack(gl_FragCoord.z);
}

如图，Shadow Map 记录了 Light Camera 所看到的最近深度图，颜色越深，离摄像机越近：

2. 第二个Pass：
   • 主摄像机需要渲染屏幕每个像素时，该像素对应的世界坐标进行 Light Camera 的MVP变换后能得到在 Light Camera 屏幕空间中的对应位置 \(shadowCoord = (x‘,y‘,z‘)\)。
   • Shadow Map 里用\((x‘,y‘)\)采样得到的遮挡深度 \(depth\) 与深度值 \(z‘\) 做比较： 若 \(depth < z‘\)（意味着该像素的光被遮挡），这时就可以对该像素降低可见度（Visibility）。

// phongVertex.glsl
// ...
void main(void) {
  vNormal = (uModelMatrix * vec4(aNormalPosition, 0.0)).xyz;
  vTextureCoord = aTextureCoord;
  vFragPos = (uModelMatrix * vec4(aVertexPosition, 1.0)).xyz;
  vPositionFromLight = uLightMVP * vec4(aVertexPosition, 1.0);
  gl_Position = uProjectionMatrix * uViewMatrix * uModelMatrix * vec4(aVertexPosition, 1.0);
}

// phongFragment.glsl 
// ...
void main(){
    // 归一化坐标
    vec3 projCoords = vPositionFromLight.xyz / vPositionFromLight.w;
	vec3 shadowCoord = projCoords * 0.5 + 0.5;
    // Shadow
	float visibility = 1.0;
	float depthInShadowmap = unpack(texture2D(shadowMap,shadowCoord.xy).rgba); //将rgba四通道（32位）的值unpack成float类型的数值
	if(depthInShadowmap < shadowCoord.z){
    	visibility = 0.0;
	}
    // blinnPhong光照着色
    vec3 color = blinnPhong();
    
    gl_FragColor = vec4(color * visibility,1.0);
}

如图为主摄像机每个像素经过变换后比较深度的结果，其中绿色点意味着深度 \(depth \approx z’\) （没有遮挡光照），非绿色点意味着 \(depth < z‘\)（被遮挡了光照）:

Shadow Bias

直接使用Shadow Map可能会在不应该出现阴影的位置出现一些黑白条纹相间的现象（称为 Shadow Acne）：

其本质原因在于，Shadow Map 是一个二维数组，离散的存储方式很难完全表示实际的几何信息。尤其当光照方向不垂直于平面时，遮挡深度的采样会和实际深度产生偏差（如图一个不受遮挡的几何平面，但黑色加粗部分却被Shadow Mapping方法认为是被遮挡的）：

解决方法：

• 直接给采样阴影深度加一个 偏移量 Bias（相当于把阴影深度往远处加，从而更不容易产生遮挡）。

// phongFragment.glsl
//...
void main(){
    // 归一化坐标
    vec3 projCoords = vPositionFromLight.xyz / vPositionFromLight.w;
	vec3 shadowCoord = projCoords * 0.5 + 0.5;
    // Shadow Bias
	const float BIAS = 0.005;
    // Shadow
	float visibility = 1.0;
	float depthInShadowmap = unpack(texture2D(shadowMap,shadowCoord.xy).rgba); //将rgba四通道（32位）的值unpack成float类型的数值
	if(depthInShadowmap + BIAS < shadowCoord.z){
    	visibility = 0.0;
	}
    // blinnPhong光照着色
    vec3 color = blinnPhong();
    
    gl_FragColor = vec4(color * visibility,1.0);
}

Peter Panning 问题 & 简单 Trick

然而由于增加了Bias，可能会导致 Peter Panning 现象：往往在物体缝隙间发生漏光。

解决方法：

• 避免使用单薄的几何体（例如薄墙、薄地面）；只要几何体厚度大于Bias，影子边界便会产生在几何体内部，从而不易看见影子与几何体的分离现象。

有一种有别于Bias的方法（但实际上也是殊途同归）：

• 不使用Bias

• 第一个Pass（Light Camera记录深度的那个）设置成仅渲染背面（正面剔除）

这样可以让一些具有厚度的几何体背面作为深度记录，从而部分避免了几何体正面的 Shadow Acne现象。实际上这个跟使用了Bias+加厚几何体思想是差不多的，区别只不过在于：前者是低门限加一个偏移，后者则是直接给出高门限

Slope Scale Based Depth Bias

通过上面知道，Bias 过小时可能不能解决 Shadow Acne 现象，Bias 过大时又可能导致严重的 Peter Panning问题。

Slope Scale Based Depth Bias ：为了尽可能减少由于 Bias 过大过小引起的问题，采取了根据平面倾角的一种自适应 Bias（例如：当光线与平面垂直时，Bias应该为0；当光线与平面的夹角越小，则Bias应越大）。

float bias = max(0.05 * (1.0 - dot(normal, lightDir)), 0.005);

Percentage Closer Filtering（PCF）

Shadow Mapping 还存在 阴影锯齿（Shadow Aliasing） 问题：

Percentage Closer Filtering（PCF）正是解决阴影锯齿的方案，它的核心想法是计算阴影时不是考虑单个采样点，而是在一定范围内进行多重采样，这样可以让阴影的边缘不那么锯齿，因为 Visibility 不再是非0即1，而是带有渐变的取值。

分布采样函数

vec2 disk[NUM_SAMPLES]; // 经过分布采样函数运算后得到NUM_SAMPLES个采样坐标

在对周围一定范围内若干个坐标进行采样的时候，可以通过分布采样函数来确定 NUM_SAMPLES 个采样位置，为了让阴影边缘更加柔和，我们可以用一些较好的分布采样函数。

均匀圆盘分布采样（Uniform-Disk Sample）：圆范围内随机取一系列坐标作为采样点；看上去比较杂乱无章，采样效果的 noise 比较严重。

泊松圆盘分布采样（Poisson-Disk Sample）：圆范围内随机取一系列坐标作为采样点，但是这些坐标还需要满足一定约束，即坐标与坐标之间至少有一定距离间隔。

// 均匀圆盘分布
void uniformDiskSamples( const in vec2 randomSeed ) {
  // 随机种子
  float randNum = rand_2to1(randomSeed);
  // 随机取一个角度
  float sampleX = rand_1to1( randNum ) ;
  float angle = sampleX * PI2;
  // 随机取一个半径
  float sampleY = rand_1to1( sampleX ) ;
  float radius = sqrt(sampleY);
  for( int i = 0; i < NUM_SAMPLES; i ++ ) {
    disk[i] = vec2(radius * cos(angle) , radius * sin(angle));
    // 继续随机取一个半径
    sampleX = rand_1to1( sampleY ) ;
    radius = sqrt(sampleY);
    // 继续随机取一个角度
    sampleY = rand_1to1( sampleX ) ;
    angle = sampleX * PI2;
  }
}

// 泊松圆盘分布
void poissonDiskSamples( const in vec2 randomSeed ) {
  // 初始弧度
  float angle = rand_2to1( randomSeed ) * PI2;
  // 初始半径
  float INV_NUM_SAMPLES = 1.0 / float( NUM_SAMPLES );
  float radius = INV_NUM_SAMPLES;
  // 一步的弧度
  float ANGLE_STEP = 3.883222077450933;// (sqrt(5)-1)/2 *2PI
  // 一步的半径
  float radiusStep = radius;

  for( int i = 0; i < NUM_SAMPLES; i ++ ) {
    disk[i] = vec2(cos(angle),sin(angle)) * pow( radius, 0.75 );
    radius += radiusStep;
    angle += ANGLE_STEP;
  }
}

PCF 算法过程

Percentage Closer Filtering（PCF） 的算法过程：

1. 计算 Visibility 时，原本对 Shadow Map 的一次坐标采样换成对周围一定范围内若干个坐标进行采样。
2. 各个采样结果同样用来与 \(z‘\) 做比较，最后取比较结果的平均作为 Visibility。

float visibility_PCF(sampler2D shadowMap, vec4 coords) {
  const float bias = 0.005;
  float sum = 0.0;
  // 初始化泊松分布
  poissonDiskSamples(coords.xy);
  // 采样
  for(int i = 0;i<NUM_SAMPLES;++i){
    float depthInShadowmap = unpack(texture2D(shadowMap,coords.xy+poissonDisk[i]*0.001).rgba);
    sum += ((depthInShadowmap + bias)< coords.z?0.0:1.0);
  }
  // 返还平均采样结果
  return sum/float(NUM_SAMPLES);
}

Percentage Closer Soft Shadows（PCSS）

Shadow Mapping 还存在硬阴影（Hard Shadow）的问题，因为现实世界的影子往往是软阴影（Soft Shadow）。

一个现实观察是，当投影物与阴影之间的距离越远，则阴影越软（如下图：笔尖阴影由于与笔尖的距离较近，因此阴影边缘较为锐利；而远处笔身阴影则因与笔身距离较远，阴影边缘较为发散且模糊）。

这是因为较大的光源面会有一些区域被遮蔽一部分光又接受一部分光，从而产生半影（Penumbra），直观看就是没那么暗的边缘处阴影。

Penumbra Size

用二维平面的图去描述，实际上就是光源段 \(w_{Light}\) 两端与遮挡物连直线后打在被投影物上的即是 半影段 \(w_{Penumbra}\) ，也就是说这段半影需要有渐变的阴影效果。假如我们用 PCF 算法中的圆盘半径大小等同于这个半影段的尺寸 \(w_{Penumbra}\)，就能实现这段的渐变阴影效果（可以想想为什么）。

现在，由下图的几何关系容易推出：

\(w_{\text {Penumbra }}=\left(d_{\text {Receiver }}-d_{\text {Blocker }}\right) \cdot w_{\text {Light }} / d_{\text {Blocker }}\)

其中，\(w_{Light}\) 是光源面积尺寸，\(d_{Blocker}\) 是遮挡物的深度，\(d_{Receiver}\) 是被投影物（实际上就是shading point）的深度。

但是 PCF 算法的圆盘半径大小是固定的，因此处处的边缘看起来都带有相同的渐变范围，这和我们看到的笔尖阴影现象不符合（近处边缘渐变应该更少些，远处边缘渐变应该多些），所以我们可以只要根据不同位置动态地修改圆盘半径大小（实际上就是动态地计算 \(w_{Penumbra}\) ），这个也就是PCSS的核心部分。

Blocker Search

我们不能简单把一个投影点变换成Shadow Map的坐标后，直接拿单个坐标采样 ShadowMap 的深度来作为 \(d_{Blocker}\) 。这是因为投影点的单次采样实际上就是单一直线连向了光源面的中心，而这条直线要是没有碰到遮挡物（即 \(d_{Blocker}=d_{Receiver}\) ），从而得出该投影点为全亮的结论。

但实际很多场景中（如下图），投影点和光源面处处连线后会发现有相当一部分光线会碰到遮挡物，因此该投影点应该属于半影范围内。

为此，我们可以对 ShadowMap 的一定范围内进行多重采样，每次采样得到的深度若小于 \(d_{Receiver}\) 则认为遇到遮挡物并算入平均遮挡深度的贡献，这样多重采样之后得到的平均遮挡深度就作为 \(d_{Blocker}\)。

如何确定采样的范围半径呢？两个参数决定：\(w_{Light}\) 的尺寸、投影点与光源的距离（可以结合上图推理一下为什么）

\(SampleSize=w_{Light}\cdot z_{Receiver} \cdot c\)

这样，计算 Blocker 平均遮挡深度的整个过程为：

float findBlocker( sampler2D shadowMap, vec2 uv, float zReceiver ) {
  float dBlocker = zReceiver * 0.01;
  const float wLight = 0.006;
  const float c = 100.0;	
  float sampleSize = wLight * zReceiver * c;
  float sum = 0.01;	// 取0.01一是为了避免出现0除问题，二是当多重采样没有贡献时的dBlocker/sum将等于zReceiver
  for(int i = 0;i<BLOCKER_SEARCH_NUM_SAMPLES;++i){
    float depthInShadowmap = unpack(texture2D(shadowMap,uv+poissonDisk[i]*sampleSize).rgba);
    if(depthInShadowmap < zReceiver){
      dBlocker += depthInShadowmap;
      sum += 1.0;
    }
  }
  return dBlocker/float(sum);
}

PCSS 算法过程

Percentage Closer Soft Shadows（PCSS） 的算法过程：

1. Blocker Search：通过多重采样，计算出平均遮挡深度 \(d_{Blocker}\)

2. Penumbra Size：计算圆盘半径大小 \(w_{\text {Penumbra }}=\left(d_{\text {Receiver }}-d_{\text {Blocker }}\right) \cdot w_{\text {Light }} / d_{\text {Blocker }}\)

3. Filtering：通过多重采样，计算出平均 Visibility（实际上就是调用PCF算法）

float visibility_PCSS(sampler2D shadowMap, vec4 coords){
  poissonDiskSamples(coords.xy);
  // STEP 1: avgblocker depth
  float dBlocker = findBlocker(shadowMap,coords.xy,coords.z);
  // STEP 2: penumbra size
  const float wLight = 0.006;
  float penumbra = (coords.z-dBlocker)/dBlocker * wLight;
  // STEP 3: filtering
  const float bias = 0.005;
  float sum = 0.0;
  for(int i = 0;i<PCF_NUM_SAMPLES;++i){
    float depthInShadowmap = unpack(texture2D(shadowMap,coords.xy+disk[i]*penumbra).rgba);
    sum += ((depthInShadowmap + bias)< coords.z?0.0:1.0);
  }
  return sum/float(PCF_NUM_SAMPLES);
}

PCF算法效果图：

PCSS算法效果图：

Variance Soft Shadow Mapping（VSSM）

PCSS、PCF 的算法都需要多重采样，尤其 PCSS 需要两个多重采样（第一步的Blocker Search和第三步的PCF），这使得算法速度较慢。

为了避免多重采样的计算，Variance Soft Shadow Mapping（VSSM） 假定一定范围内的深度的分布符合 正态分布（Normal Distribution） ，那么只要知道该段范围的 均值（实际上就是期望值）E 、方差 Var，就能先得到该范围的正态分布模型（即知道对应的 概率密度函数 PDF）。

\(PDF(x) = \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{(x-\mu)^{2}}{2 \sigma^{2}}\right)\)

其中，\(\mu = E\)，\(\sigma^2 = Var\)。

接着可以通过该正态分布模型的 累计分布函数（即 CDF），就能快速推算出该范围内有多少比例的 x 大于（或小于）给定的某个值。

\(CDF(x) = \int^x_{-\infin} PDF(t) \mathrm{d}t\)

Variance Soft Shadow Mapping（VSSM） ：简单来说，VSSM 算法就是依据 ShadowMap 的深度符合正态分布的假设来快速完成 PCSS 中的第一步（Blocker Search）和第三步（PCF算法）的一种阴影算法。

VSSM效果图：

计算平均值 & 方差

为了快速查询得到某段范围的均值、方差，我们可以先选以下一种数据结构来快速查询 Shadow Map 某段范围的均值（期望值）\(E(X)\) 。

• 硬件 Mipmap：当 Shadow Map 更新时，需要重新生成 Mipmap，不过GPU硬件实现的 Mipmap 算法非常快的开销非常小；查询某段方形范围时，需要根据方形中心所在的位置（相对于周围四个纹素的坐标）、上下层级做三线性插值（Trillinear interpolation），得到的结果即是近似的均值（期望值）。

• 前缀和数组（Summed Area Tables/SAT）：当 Shadow Map 更新时，需要重新进行二维前缀和计算；需要编写 Compute Shader 实现该算法，比Mipmap方法更慢一些，但百分百精准；查询某段方形范围时，就可以通过如下图方法快速查询得到某段范围的总和，除掉范围面积就能得到均值（期望值）。

我们需要存储 \(E(X)\)、\(E(X^2)\)? ，这样就能计算某段范围的平均值、方差：

1. 平均值 \(E(X)\)

2. 方差 \(Var(X)=E(X^2)-E^2(X)\)

\(E(X^2)\) 即 ShadowMap 每个纹素再求个平方后作为额外的ShadowMap，然后再生成 Mipmap 或 SAT。

计算累计分布函数（CDF）

有了上面的期望值与方差，我们就能确定一个正态分布。但是它对应的 CDF 函数是没有解析解的，而有数值解（称为 Error Function），但是计算比较繁琐。

**切克比夫不等式（Chebyshev’s Inequality） **：\(P(x>t) \leq \frac{\sigma^{2}}{\sigma^{2}+(t-\mu)^{2}}\)

实际上这个切克比夫不等式不仅可用在正态分布，其它的很多分布也是可以套用这个不等式的。

将这个不等式改造一下，就成了一个大胆的近似公式：

\(P(x>t) \approx \frac{\sigma^{2}}{\sigma^{2}+(t-\mu)^{2}}\)

注意：这里求的是 \(x>t\) 的部分，即 \(P(t)=1-CDF(t)\)。

当然这个近似公式肯定不是精确的，但是计算开销非常小，也就被用在 VSSM 算法中。

加速 Blocker Search 算法

PCSS 算法中的 Blocker Search 步骤：在一定范围内多重采样，每次采样得到的深度若小于 \(d_{Receiver}\) 则认为遇到遮挡物并算入平均遮挡深度的贡献，这样多重采样之后得到的平均遮挡深度 \(z_{occ}\) 就作为 \(d_{Blocker}\)。

如下图5X5的采样结果若设 \(d_{Receiver}\) 为7，那么平均遮挡深度 \(z_{occ}\) 则为红色部分的平均值。

设该采样范围的面积为 \(N\)，无遮挡的面积占有 \(N_1\)，有遮挡的面积则占有 \(N_2\) ，则有：

\(\frac{N_{1}}{N} z_{\text {unocc }}+\frac{N_{2}}{N} z_{o c c}=z_{A v g}\)

我们做出两个假设：

• \(\frac{N_1}{N} = P(x>d_{Receiver})\) ，\(\frac{N_2}{N} = 1-P(x>d_{Receiver})\) ；这个假设基于认为深度分布为正态分布，通过切克比夫不等式获得近似解（即上面两节的内容）。

• \(z_{unocc} = d_{Receiver}\) ；这个假设基于认为绝大部分没被遮挡的情况都属于同一个深度（相当于在同一个垂直于光方向的平面），即可认为均为深度 \(d_{Receiver}\)。

那么 VSSM 加速该算法的公式表示为：

\(d_{Blocker} = z_{occ} = \frac{N\cdot z_{Avg} - N_1 \cdot z_{unocc}}{N_2} = \frac{E(x)-P(x>d_{Receiver})\cdot d_{Receiver}}{1-P(x>d_{Receiver})}\)

加速 PCF 算法

PCF 算法中的多重采样：每次采样得到的遮挡物深度用来与 \(z‘\) 做大小比较（小于 \(z‘\) 则视为被遮挡，大于 \(z‘\) 则视为全亮），最后取比较结果的平均作为 Visibility。

我们做出一个假设：

• \(不被遮挡的概率 = P(x > z‘)\) ；这个假设基于认为深度分布为正态分布，通过切克比夫不等式获得近似解。

那么 VSSM 加速该算法的公式表示为：

\(Visibility = P(x>z‘) \cdot 1 + (1-P(x>z‘)) \cdot 0 = P(x>z‘)\)

VSSM 的缺陷

VSSM 的主要缺陷表现：

• 并不是任何深度的分布都是符合正态分布模型的，例如对于图右的简单几何体反而用正态分布表示会很不适合。

• 漏光（Light Leaking）现象，在一些应当被阴影完全遮蔽的内部有可能仍产生亮度。

• 在加速 Blocker Search 算法中的假设 \(z_{unocc} = d_{Receiver}\) 基于认为绝大部分没被遮挡的情况都属于同一个深度，但实际上有些不被遮挡的地方深度并不等于 \(d_{Reveiver}\) 。

Moment Shadow Mapping

Moment Shadow Mapping 正是为了解决 VSSM 缺陷的一种算法，它主要想法是：使用高阶的矩去描述一个分布的 CDF。这样就能通过记录 m 阶的矩，就能复原成足够接近实际 CDF 函数的效果，从而能适应不同的深度分布模型（有些地方可能接近正态分布，有些地方可能奇奇怪怪的分布）。

Moment Shadow Mapping将使用最简单的形式来标识矩：\(z,z^2,z^3,z^4,...\)

实际上，VSSM 本质便是记录 2 阶的矩来复原 CDF 函数，而 Moment Shadow Mapping 一般使用4阶的矩就已经足够接近实际 CDF 了。

虽然 Moment Shadow Mapping 效果相当不错，很好的解决了 VSSM 绝大部分缺陷，但是它仍需要相当的额外空间开销和重建矩的额外性能开销。

Distance Field Soft Shadows

Distance Field Soft Shadows 是与 Shadow Mapping 系列技术（PCF、PCSS、VSSM、Moment Shadow Mapping）截然不同的阴影技术路线，它主要想法是：

将点 \(o\)（Shading Point）与光源面中心点 \(p_{light}\) 相连形成一条方向为 \(l\) 的中心线段，而这条中心线上各个点 \(p_i\) 都可以通过 SDF 查得与其最近几何物体的距离并且推算出安全角度（点\(o\) 能打到光源面的直线与中心线的最大夹角）为 \(\theta_i = arcsin \frac{\operatorname{SDF}(p_i)}{p_i-o}\)

SDF 相关可以看几何（Geometry）部分，这里假定已经对场景生成了 SDF 信息。

那么所有这些点中对应的安全角度之中取最小的安全角度 \(\theta = min\{\theta_i\}\) ，这个安全角度与最大角度的比例决定了光源面的光照覆盖率，也就决定了点 \(o\) 的Visibility。

使用 Distance Field Soft Shadows 的好处很多：

• 计算阴影很快（假设已经生成了SDF的情况下，比传统Shadow Mapping类技术是要快的多）
• 阴影质量很高，而且完美解决 Shadow Ance / Peter Panning / 采样噪声等传统Shadow Mapping会出现的问题

然而代价是：

• SDF 需要预计算，这就意味着场景物体需要是静态的，当然也可以使用一些算法使能和动态物体相结合，尽量减少重新生成SDF的成本。
• SDF 需要较大的存储空间（一般采用三维数组表示空间各个网格的SDF值，但是可以使用八叉树等空间数据结构或者其它方法做进一步优化）。

计算安全角度

计算某个点 \(p_i\) 的安全角度时，直观的几何关系便是：

\(\theta_i = \arcsin \frac{\operatorname{SDF}(p_i)}{p_i-o}\)

而在实践中，往往会使用：

\(\theta_i = \min \left\{\frac{k \cdot \operatorname{SDF}(p_i)}{p_i-o}, 1.0\right\}\)

这样的近似公式实际效果相当接近原几何关系，而且也能减少复杂的 arcsin 运算开销，最后它还能通过 \(k\) 这个参数来调整阴影的硬软程度。

如下图分别为 \(k=32\) 、\(k=8\)、\(k=2\) 的效果：

Distance Field Soft Shadows 算法过程

具体算法过程：

1. 将 \(o\) 点（shading point）设为第一个步进点，即 \(p_0 = o\)

2. 每次算出下一个步进点 \(p_{i+1} = p_{i} + l \cdot SDF(p_{i})\) 并记录安全角度 \(\theta_i = \min \left\{\frac{k \cdot \operatorname{SDF}(p_i)}{p_i-o}, 1.0\right\}\)

3. 重复 "步骤2"，直到满足 \(l \cdot (p_{i+1}-p_{light}) < 0\) （即意味着已经步进到光源点背面了）

4. 取所有次步进的最小安全角度 \(\theta = min\{\theta_i\}\) ，则可见度则为 \(Visibility = \frac{\theta}{c}\) （其中 \(c\) 为点 \(o\) 与光源面连接的最大角度）

参考

• [1] GAMES202-高质量实时渲染-闫令琪
• [2] Tutorial 16 : Shadow mapping | opengl-tutorial
• [3] Chapter 8. Summed-Area Variance Shadow Maps | NVIDIA Developer
• [4] Inigo Quilez :: Distance Field Soft Shadows (iquilezles.org)

实时阴影技术（Real-time Shadows）