【GAMES202】作业1—实时阴影

作业1的内容主要是实现实时渲染管线中的阴影计算，在作业1的assignment中可以看到一共包含如下三种阴影的实现：

硬阴影

Two Pass Shadow Map 方法
软阴影

PCF (Percentage Closer Filter)

PCSS (Percentage Closer Soft Shadow)

$L_0(p,\omega_o) \approx \frac{\int_{\Omega^+}V(p,\omega_i)d\omega_i}{\int_{\Omega^+}d\omega_i}\cdot \int_{\Omega^+}{L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta_id\omega_i}$

如上式所示，后半部分的积分是作业 0 中完成的 Blinn-Phong 模型，前半部分是要实现的阴影方法。这里作业框架提供了blinnPhong() 函数，在实现了阴影算法，得到可见性信息后，将两部分相乘得到最终的着色结果。

代码中对应部分如下：

//visibility = useShadowMap(uShadowMap, vec4(shadowCoord, 1.0));
//visibility = PCF(uShadowMap, vec4(shadowCoord, 1.0));
//visibility = PCSS(uShadowMap, vec4(shadowCoord, 1.0));

vec3 phongColor = blinnPhong();
gl_FragColor = vec4(phongColor * visibility, 1.0);

一、原理介绍

1、Shadow Mapping

用于渲染一个点光源在场景中投射的阴影。

Pro：

只需要Shadow Map而不需要实际场景中的几何，就可以得到阴影

Con：

自遮挡现象
锯齿

【一些问题】

（1）自遮挡

精度问题，产生自遮挡，因为Shadow Map 记录的深度不连续。

在光源和地板垂直时没有这个问题，在光源和地板平行时（夕阳西下）容易出现这个问题。

【解决方法】

加一个较小的 Bias，通过偏移值bias，将被离散的距离遮挡部分剔除。

但此时又出现了新的问题：

存在阴影失真、阴影悬浮现象。针对该现象，学术界提出以下解决方法：

对于每个Shadow Map，不仅存最小深度，还存次小深度，用它们的中间深度来进行比较。

【缺点】

要求物体都是watertight，有正面一定有反面，比如地板就不是watertight

（2）锯齿 Aliasing

通过PCF解决

2、硬阴影Two Pass Shadow Map

渲染场景两次：

第一次从光源所在方向看向这个场景，生成Shadow Map，只记录最浅的深度

第二次使用上面的pass渲染出的纹理，判断是否处于阴影中

当前点到light的距离和上一个pass中的距离一致，则可见，若Shadow Map上的值更小，则被遮挡

3、一个重要约等式

实时渲染领域的一个重要约等式：

$\int _\Omega f(x)g(x)dx \approx \frac{\int_\Omega f(x)dx}{\int_\Omega dx}\cdot \int_{\Omega}g(x)dx$

于是就可以将渲染方程进行一些变化：

$L_o(p,\omega_o) = \int_\Omega L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta_i V(p,\omega_i)d\omega_i \tag{原方程}$ $L_o(p,\omega_o)\approx \frac{\int_{\Omega^+}V(p,\omega_i)d\omega_i}{\int_{\Omega^+}d\omega_i} \int_\Omega L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)cos\theta_i d\omega_i \tag{近似后}$

对于点光源和方向光源，该变换较为准确

4、软阴影PCF

对阴影边缘做 anti-aliasing

PCF一开始是用来做抗锯齿的
对阴影的比较结果进行一个 filtering，可以理解为求平均

并不是对 Shadow进行filtering

Texture filtering仅仅是对颜色的比较进行filtering，一开始会得到一个较为模糊的Shadow Map
对深度值进行均值化，然后进行比较，可见性的依然是二值（非1即0）

5、软阴影PCSS

黄色线段代表光源
绿色线段代表遮挡物（或者说是物体自身）
蓝色线段代表阴影的接收方

图中存在一对相似三角形，d_Blocker 的长度会影响 W 的长度，d_Blocker越长， W越小，通俗的讲就是物体离光源越近，半影范围越大，filter越大，阴影看起来就更“软”。

【翻译成如下公式】

$W_{Penumbra} = W_{light} * \frac{(d_{Receiver} - d_{Blocker})}{d_{Blocker}}$

面光源无法生成Shadow Map，常用点光源的方式去生成Shadow Map

【步骤】

找到shading point 在Shadow Map上的点，然后再Shadow Map上确定一个范围，把这个范围内每一个像素的深度信息与shading point的实际深度进行比较，判断shading point是否在阴影中，若在，记录所有blocker的深度，累加取平均
知道了blocker的距离以后，就可以用上面的公式算出filter的大小
再次使用PCF

一个新的问题，该用哪个区域进行blocker的查询呢？

估计上图红色区域内，所有能挡住shading point的点，计算它们对应深度的平均值。使用自适应的采样半径，将采样半径与光源宽度、视锥体宽度关联起来，从shading point连向光源，计算在Shadow Map上覆盖的区域，顶点距离光源越远，红色区域就越大，应使用更大的采样半径，反之则使用更小的采样半径。

二、具体实现

1、光源的变换矩阵实现

我们需要以光源的位置对场景进行渲染，得到一张Shadow Map，相当于将一个摄像机摆在光源的位置，朝向光源的方向，根据需求确定渲染的范围，然后对场景进行渲染，得到在该位置的每个像素上的最近点的深度值。在 ShadowMaterial.js 中需要向 Shader 传递正确的 uLightMVP 矩阵，其调用了 light 中的 CalcLightMVP 函数。

该函数如下:

 CalcLightMVP(translate, scale) 
{
   	let lightMVP = mat4.create();
   	let modelMatrix = mat4.create();
   	let viewMatrix = mat4.create();
   	let projectionMatrix = mat4.creat	
    
   	// Model transform
    
   	// View transform
    
   	// Projection transform
    
   	mat4.multiply(lightMVP, projectionMatrix, viewMatrix);
   	mat4.multiply(lightMVP, lightMVP, modelMatr	
   	return lightMVP;
}

我们需要完成这三个变换矩阵

（1）Model transform

模型矩阵的构建先调用位移方法，再调用缩放方法，没有调用旋转方法

// Model transform
mat4.identity(modelMatrix);
mat4.translate(modelMatrix, modelMatrix, translate);
mat4.scale(modelMatrix, modelMatrix, scale);

（2）View transform

对于一个光源，其包含 lightPos、focalPoint、lightUp 三个变量，这里就可以直接调用 glMatrix 库提供的接口 mat4.lookAt()

lookAt 方法可以根据上述三个信息构建出从当前坐标系转换到观察坐标系的矩阵，其构建过程可以参考LearnOpenGL网站中的方法：摄像机/观察空间](https://learnopengl-cn.github.io/01 Getting started/09 Camera/))

1
2
3

// View transform
mat4.identity(viewMatrix);
mat4.lookAt(viewMatrix, this.lightPos, this.focalPoint, this.lightUp);

(3) Projection transform

作业指导中在此处推荐使用正交投影

这里依旧可以直接调用 glMatrix 库提供的接口 mat4.ortho()

1 2	// Generates a orthogonal projection matrix with the given bounds (static) ortho(out, left, right, bottom, top, near, far) → {mat4}

1
2
3

// Projection transform
mat4.identity(projectionMatrix);
mat4.ortho(projectionMatrix, -150, 150, -150, 150, 1e-2, 400);

2、坐标空间转换

由作业提示可知，首先需要将 vPositionFromLight 的坐标范围映射到 0 ~ 1 之间，因为纹理采样的坐标范围是 0 ~ 1 ，Shadowmap 中存储的纹理深度范围也是 0 ~ 1

1 2	vec3 shadowCoord = vPositionFromLight.xyz / vPositionFromLight.w; shadowCoord.xyz = (shadowCoord.xyz + 1.0) / 2.0;

然后调用 useShadowMap() 函数

1	visibility = useShadowMap(uShadowMap, vec4(shadowCoord, 1.0));

3、硬阴影

由作业提示可知，硬阴影主要需要实现phongFragment.glsl 中的 useShadowMap 函数

该函数具体实现如下:

float useShadowMap(sampler2D shadowMap, vec4 shadowCoord)
{
  // 获取距离最近的深度值
  float shadowDepth = unpack(texture2D(shadowMap,shadowCoord.xy));
  // 获取当前深度值
  float currentDepth = shadowCoord.z;
  // 将二者进行比较
  return shadowDepth + EPS <= currentDepth ? 0.0 : 1.0;
}

其中作业框架以及提供了 unpack 函数

float unpack(vec4 rgbaDepth) 
{
    const vec4 bitShift = vec4(1.0, 1.0/256.0, 1.0/(256.0*256.0), 1.0/(256.0*256.0*256.0));
    
    return dot(rgbaDepth, bitShift);
}

运行结果如下：

可以观察到，当光线方向与顶点法线方向夹角较大时容易出现阴影瑕疵（shadow acne），这是由于 ShadowMap 的精度问题（采样率低）而产生的自遮挡（self occlusion）现象，在Lecture 03中有相关的论述。

我们可以引入 bias 来对阴影纹理的采样值进行一定程度的偏移来处理自遮挡问题，但这种处理方式同样会产生漏光现象（一些本应处于阴影中的位置被意外地照亮了）

此处可参考自适应Shadow Bias算法](https://zhuanlan.zhihu.com/p/370951892))

可以在phongVertex.glsl里令 gl_Position = vPositionFromLight，就可以可视化出shadow map渲染的范围

4、软阴影 PCF

PCF的实现需要完善代码框架中的 PCF 函数

![下载](C:/Users/EKKO/Desktop/%E4%B8%8B%E8%BD%BD.png)float PCF(sampler2D shadowMap, vec4 coords) 
{
  return 1.0;
}

对采样点周围的一片区域进行随机采样

$被遮挡的程度 = \frac{被采样点遮挡的顶点数量}{采样点的总数量}$

关于随机采样的方法，框架中提供了泊松圆盘采样和均匀圆盘采样两种采样函数，参考链接](https://codepen.io/arkhamwjz/pen/MWbqJNG?editors=1010))中还给出了两种采样方式的可视化展示

【泊松圆盘采样 poisson disk】

【均匀圆盘采样 uniform disk】

两种采样方法都在一个单位圆内随机地生成二维坐标向量，并根据 NUM_SAMPLES 的值确定采样点的数量，最后将生成的采样坐标存储在 poissonDisk 数组中，用于后续在PCF中对着色点对应的SM坐标加上该数组中的值来得到采样点的SM坐标。

上述两种采样方式对应的代码如下:

vec2 poissonDisk[NUM_SAMPLES];

void poissonDiskSamples( const in vec2 randomSeed ) 
{
  float ANGLE_STEP = PI2 * float( NUM_RINGS ) / float( NUM_SAMPLES );
  float INV_NUM_SAMPLES = 1.0 / float( NUM_SAMPLES );

  float angle = rand_2to1( randomSeed ) * PI2;
  float radius = INV_NUM_SAMPLES;
  float radiusStep = radius;

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

void uniformDiskSamples( const in vec2 randomSeed ) 
{
  float randNum = rand_2to1(randomSeed);
  float sampleX = rand_1to1( randNum ) ;
  float sampleY = rand_1to1( sampleX ) ;

  float angle = sampleX * PI2;
  float radius = sqrt(sampleY);

  for( int i = 0; i < NUM_SAMPLES; i ++ ) {
    poissonDisk[i] = vec2( radius * cos(angle) , radius * sin(angle)  );

    sampleX = rand_1to1( sampleY ) ;
    sampleY = rand_1to1( sampleX ) ;

    angle = sampleX * PI2;
    radius = sqrt(sampleY);
  }
}

根据随机采样的得到的uv偏移值，继续完善 PCF 函数，遍历 poissonDisk 数组，为每一个uv采样点的坐标加上偏移值，然后将采样得到的深度值与当前顶点在 light space 的深度值比较，若判断顶点存在于阴影中，则将 blocker 的值+1。

最后可见程度 = 1 - 存在于阴影中的采样点 / 采样点总数

float PCF(sampler2D shadowMap, vec4 coords) 
{
    float blocker = 0.0;  // 统计被遮挡的数量
    float stride = 5.0;
    float textureSize = 2048;          // 2048为生成的纹理的分辨率
    poissonDiskSamples(coords.xy);

    for(int i = 0; i < NUM_SAMPLES; i++)
    {
        vec2 uvBias = poissonDisk[i] * stride / textureSize;
        float shadowDepth = unpack(texture2D(shadowMap, coords.xy + uvBias));
        if( coords.z > shadowDepth + EPS)
        {
            blocker++;
        }            
    }
    float visibility = 1.0 - blocker / float(NUM_SAMPLES);
    return visibility;
}

效果图如下:

【stride = 5 ， poissonDiskSamples】

【stride = 5 ， uniformDiskSamples】

【NUM_SAMPLES = 200， stride = 5 ，poissonDiskSamples】

5、软阴影 PCSS

PCSS 依据现实当中观测到的物体产生的阴影离物体越近，就会感觉越阴影越”浓”，离物体越远，就会感觉阴影越”淡“，PCSS就是用来还原这一现象的，具体原理以及在上文提过，此处不再赘述。

有作业说明可知，PCSS要求完成 phongFragment.glsl 中的 findBlocker 和 PCSS 函数

float findBlocker( sampler2D shadowMap,  vec2 uv, float zReceiver ) 
{
	return 1.0;
}

float PCSS(sampler2D shadowMap, vec4 coords)
{
  // STEP 1: avgblocker depth，遮挡点以及其周围的遮挡点的平均深度

  // STEP 2: penumbra size，利用相似三角形的原理计算新的半影直径

  // STEP 3: filtering
  
  return 1.0;
}

【findBlocker】

该函数需要完成对遮挡物平均深度的计算，同PCF一样，首先是用采样函数生成随机偏移数组，对于每一个偏移的 uv 坐标采样得到深度，将产生了遮挡的部分进行累加，然后计算得到被遮挡的平均深度。此处uv的偏移对结果影响较大，若选择的stride的值较小，则对于距离光源较远的部分而言，打到该部分上的光更容易被遮挡，从而影响平均遮挡深度的计算。若选择的stride的值较大，会造成顶点周围的平均遮挡深度相似，无法表现出距离越远阴影越模糊的效果。

float findBlocker( sampler2D shadowMap,  vec2 uv, float zReceiver ) 
{
  // 定义相关参数
  float blocker = 0.0;
  float blockDepth = 0.0;
  float textureSize = 2048.0;
  float stride = 5.0;

  // 泊松采样
  poissonDiskSamples(uv);

  // 对被遮挡的部分进行累加
  for(int i = 0 ;i < NUM_SAMPLES; i++)
  {
     vec2 uvBias = poissonDisk[i] * stride / textureSize;
     float shadowDepth = unpack( texture2D(shadowMap, uv + uvBias));
     if( shadowDepth + 0.01 <= zReceiver)
     {
        // 在阴影中
        blocker++;
        blockDepth += shadowDepth;
     }
  }

  if( blocker < 0.1 )
  {
    return 1.0;
  }

  float avgblocker = blockDepth / blocker;
  return avgblocker;
}

【PCSS】

由相似三角形计算半影范围的公式:

$W_{Penumbra} = W_{light} * \frac{(d_{Receiver} - d_{Blocker})}{d_{Blocker}}$

具体代码如下:

float PCSS(sampler2D shadowMap, vec4 coords)
{

  float w_Light = 1.0;
  float d_Receiver = coords.z;

  // STEP 1: avgblocker depth
  float d_Blocker = findBlocker(shadowMap, coords.xy, coords.z);

  // STEP 2: penumbra size 可根据公式计算得到
  float w_penumbra = w_Light * ( d_Receiver - d_Blocker ) / d_Blocker;

  // STEP 3: filtering 
  // 对上面的 PCF 进行一些变化: PCF 的uvBias不再固定，而是根据上面计算出来的半影直径动态地调整
  // 采样点距离遮挡点越远，其偏移半径越大
  float visibility = 0.0;
  for( int i = 0; i < NUM_SAMPLES ; i++)
  {
      vec2 uvBias = poissonDisk[i] * 10.0 / 2048.0 * w_penumbra;
      float shadowDepth = unpack( texture2D(shadowMap, coords.xy + uvBias ));
      if( coords.z < (shadowDepth + 0.01) )
      {
         visibility += 1.0;
      }
  }
  return visibility / float(NUM_SAMPLES);
}

效果图如下：

采样数为20，stride为10，可以看到噪声明显

采样数为20，stride为10，可以看到噪声情况有明显改善

经过一番调参，在采样数为200的情况下，stride为25时效果较好