联级阴影贴图CSM(Cascaded shadow map)原理与实现

CSM是利用分层的ShadowMap技术，实现大场景的阴影算法。示意图如下图：

我们通过给眼视锥分片，为每个分片生成一个相同分辨率的光空间深度图。利用眼睛看到的片段，根据其所在分片位置，转换为相应光空间深度，和光空间深度图比对，若深度大于深度图内的深度，则处于阴影。

该算法分为4个步骤：

将眼视锥分成多个深度分片。
将光的视锥体分成多个较小的视锥体，每个视锥体(也可以称为包围盒)都覆盖一个分片。
为每个光视锥体渲染一个阴影贴图。
渲染整个场景的阴影。

根据以上步骤，我们首先构建类：

#define CSM_MAX_SPLITS 4class ShadowMap{
public:ShadowMap();~ShadowMap();void init(Camera* camera);//生成CSM所需的资源(分片光空间变换矩阵，眼空间到光空间变换矩阵)void pre_depth_write(Camera* camera, const glm::vec3& lightdir);glm::mat4 crop_matrix(int t_split_index);glm::mat4 projection_matrix(int t_split_index);glm::mat4 modelview_matrix();int num_splits() const;int depth_tex_size() const;GLuint fbo() const;GLuint texture() const;glm::vec4 far_bounds();float* texture_matrices();private:void create_fbo();void create_texture();void update_split_distances(Camera* camera);//计算世界空间中的相机视锥切片边界点void update_split_frustum_points(Camera* camera);//生成光视锥切片的光空间变换矩阵void generate_crop_matrices(const glm::mat4& t_modelview);//更新每个分片far边界void update_far_bounds(const glm::mat4& projection);void update_texture_matrices(const glm::mat4& projection, const glm::mat4& view_inverse);GLuint m_fbo;//CSM的fboGLuint m_texture_array;//texture数组int m_num_splits;//视锥分片数int m_depth_tex_size;//深度图尺寸float m_split_weight;//划分权重(lambda)float m_far_bounds[CSM_MAX_SPLITS];//相机空间的深度切面Frustum m_frustums[CSM_MAX_SPLITS];glm::mat4 m_bias;glm::mat4 m_modelview;//光空间view矩阵glm::mat4 m_crop_matrices[CSM_MAX_SPLITS];//分片光空间变换矩阵(投影*观察)glm::mat4 m_projection_matrices[CSM_MAX_SPLITS];//光空间投影矩阵glm::mat4 m_texture_matrices[CSM_MAX_SPLITS];//眼空间转换为光空间矩阵
};

初始化函数，非配摄像机要素，以及为眼视锥分片：

void ShadowMap::init(Camera* camera) {float camera_fov = camera->Zoom;float width = camera->Width;float height = camera->Height;float ratio = width / height;//略微添加偏移(这里为0.2)，避免在边框处有伪像for(int i = 0 ; i < m_num_splits ; i++) {m_frustums[i].fov(camera_fov + 0.2f);m_frustums[i].ratio(ratio);}//光空间的深度需要m_bias将[-1，1]变换为[0，1]，从而在texture中使用。m_bias = glm::mat4(0.5f, 0.0f, 0.0f, 0.0f,0.0f, 0.5f, 0.0f, 0.0f,0.0f, 0.0f, 0.5f, 0.0f,0.5f, 0.5f, 0.5f, 1.0f);//计算摄像机视线空间中每个平截头体切片的近距离和远距离update_split_distances(camera);create_fbo();create_texture();
}

眼视锥分片

“Parallel-Split Shadow Maps for Large-scale Virtual Environments”这篇文章给出了一个非常好的眼视锥分片方法，称为PSSM，即平行分片阴影贴图。为了理解其公式由来，我们看下图：

我们的关键问题是如何指定分割位置，即计算Ci。在上图中，从表面(橙色线段)投射到屏幕上的观察光束dp的大小约为ndy/z。φ和θ分别表示曲面法线与屏幕和阴影贴图平面之间的角度。由于dy =dzcosφ/cosθ，因此阴影图锯齿误差dp/ds(阴影贴图的尺寸为ds×ds)可以表示为:

当dp大于屏幕的像素大小时，会出现阴影贴图欠采样，这可能在透视误差dz/zds或投影误差cosφ/cosθ变大时发生。文章中提出了三种分割方案，即计算Ci的方案。

对数拆分方案

理想情况下，透视误差的最佳分布可使dp/ds在整个深度范围内保持恒定。即dz/zds =ρ(忽略投影误差)，其中ρ为常数。有：

在这里，有一个重要的注意事项：阴影贴图s∈[0，1]的局部参数化意味着一个假设“阴影贴图准确地覆盖了视锥，并且在场景的不可见部分上没有浪费任何分辨率”。这意味着上述参数化只能在理论上实现。根据s∈[0，1]的假设，很容易得到ρ= ln(f / n)。因此，理论上最佳的透视阴影图参数化s(z)应满足以下对数函数：

但是，当前硬件支持的唯一非线性变换是透视变换(s = A/z+B)。为了近似对数函数，必须在几个深度层将其离散化。我们将眼视锥划分为更多的层，以使该函数具有更好的近似值。这是该文章中提出的“分片”思想的基本动机。

上述方程可以用 $z=\left \{ C_{i}^{log} \right \}$ 进行离散化:

即：

因为此分片方案旨在产生透视误差理论上均匀的分布，所以分配给[Ci，Ci-1]的分辨率应为整体纹理分辨率的i/m。将si = i/m代入上述公式可得出：

该分割方案的主要缺点是观看者附近的分割部分的长度太小，因此这些分割部分中可以包含的对象很少。这是由于理论上的最佳参数设置是假设阴影贴图准确地覆盖了视锥，并且没有任何分辨率浪费在场景的不可见部分上。换句话说，此假设要求每个z∈[n，f]必须在规范化纹理空间中映射到唯一的s∈[0,1]。但是，这种假设在实践中无法满足，以至于靠近观看者的部分发生过采样，而远离观看者的部分发生欠采样。因此，对数拆分方案很难在实践中应用。

均匀拆分方案

最简单的分割方案是沿z轴均匀放置分割平面,其分割点为：

忽略投影误差，我们获得下述阴影图锯齿误差：

可以看出当z减小时，阴影误差急剧增大。

加权拆分方案

文章提出的PSSM方法为将对数和均匀相加乘上二分之一，去获取较平衡的位置，而现如今的拆分多用对数和均匀的加权值，即：

引出了 $\lambda$ 去人为调节想要的权重分配，从而获得误差减小和分辨率浪费之间的平衡。根据以上描述，我们的分片点选取代码如下：

//计算摄像机视线空间中每个平截头体切片的近距离和远距离
void ShadowMap::update_split_distances(Camera* camera) {float nd = camera->near;float fd = camera->far;float lambda = m_split_weight;float ratio = fd / nd;m_frustums[0].near(nd);//使用文章”Parallel-Split Shadow Maps for Large-scale Virtual Environments“的切片方法(PSSM)//文章使用均匀法和对数法的中值计算，而我们现在可以引用lambda去调节它们的比率。for(int i = 1 ; i < m_num_splits ; i++) {float si = i / (float)m_num_splits;float t_near = lambda * (nd * powf(ratio, si)) + (1 - lambda) * (nd + (fd - nd) * si);float t_far = t_near * 1.005f;//略微增加重合，避免断裂m_frustums[i].near(t_near);m_frustums[i-1].far(t_far);}m_frustums[m_num_splits-1].far(fd);
}

值得注意的是我们避免阴影断裂，适当的延长每个视锥切片的远平面，视锥类如下：

class Frustum {
public:glm::vec3 m_points[8];//平截头体8个点Frustum();void set(float t_fov, float t_ratio, float t_near, float t_far);float fov() const;float ratio() const;float near() const;float far() const;void near(float t_near);void far(float t_far);void fov(float t_fov);void ratio(float t_ratio);private:float m_fov;float m_ratio;//宽高比float m_near;float m_far;
};
Frustum::Frustum() :m_fov(45.0),m_ratio(0.5),m_near(1.0),m_far(200.0) {
}
void Frustum::set(float t_fov, float t_ratio, float t_near, float t_far) {m_fov = t_fov;m_ratio = t_ratio;m_near = t_near;m_far = t_far;
}
float Frustum::fov() const {return m_fov;
}
float Frustum::ratio() const {return m_ratio;
}
float Frustum::near() const {return m_near;
}
float Frustum::far() const {return m_far;
}
void Frustum::near(float t_near) {m_near = t_near;
}
void Frustum::far(float t_far) {m_far = t_far;
}
void Frustum::fov(float t_fov) {m_fov = t_fov;
}
void Frustum::ratio(float t_ratio) {m_ratio = t_ratio;
}

构建光视锥体

完成了分割点的选取后，我们需要在光空间生成投影和观察矩阵，以便之后在光空间下的深度比对和渲染深度图。该部分主要分为以下四个小步骤：

更新摄像机视锥分块的世界空间位置(8个顶点)
生成光空间下的摄像机分块包围盒投影矩阵
更新分片位置(相机空间，标准化为[0,1])
更新眼空间到光空间变换矩阵

我们构建光视锥体代码如下：

//生成CSM所需的资源(分片光空间变换矩阵，眼空间到光空间变换矩阵(负灯光向量))
void ShadowMap::pre_depth_write(Camera* camera, const glm::vec3& lightdir) {glm::mat4 t_modelview = glm::lookAt(glm::vec3(0.0, 0.0, 0.0),glm::vec3(-lightdir.x, -lightdir.y, -lightdir.z),glm::vec3(-1.0f, 0.0f, 0.0f));//更新摄像机视锥分块的世界空间位置(8个顶点)update_split_frustum_points(camera);//生成光空间的摄像机分片包围盒投影矩阵generate_crop_matrices(t_modelview);m_modelview = t_modelview;glm::mat4 t_view = camera->GetViewMatrix();glm::mat4 t_view_inverse = glm::inverse(t_view);glm::mat4 t_projection = camera->GetPerspectiveMatrix();//更新分片位置(相机空间，标准化为[0,1])update_far_bounds(t_projection);//更新眼空间到光空间变换矩阵update_texture_matrices(t_projection, t_view_inverse);
}

1.更新摄像机视锥分块的世界空间位置(8个顶点)：

因为我们获得了分割点的位置，但是这个点位于眼睛空间，为了计算光投影锥体，我们获得其在世界空间的位置，并计算出每盒分块眼睛视锥体的8个顶点，用以之后的光投影锥体的构建。我们计算分块锥体8个顶点的代码如下：

//更新摄像机视锥分块的世界空间位置(8个顶点)
void ShadowMap::update_split_frustum_points(Camera* camera) {glm::vec3 center = camera->Position;glm::vec3 view_dir = -camera->Front;glm::vec3 up(0.0f, 1.0f, 0.0f);glm::vec3 right = glm::cross(view_dir, up);for(int i = 0 ; i < m_num_splits ; i++) {Frustum& t_frustum = m_frustums[i];glm::vec3 fc = center + view_dir * t_frustum.far();glm::vec3 nc = center + view_dir * t_frustum.near();right = glm::normalize(right);up = glm::normalize(glm::cross(right, view_dir));// 计算当前分片的近平面和远平面宽高的一半float near_height = tan(t_frustum.fov() / 2.0f) * t_frustum.near();float near_width = near_height * t_frustum.ratio();float far_height = tan(t_frustum.fov() / 2.0f) * t_frustum.far();float far_width = far_height * t_frustum.ratio();//记录眼视锥8个顶点t_frustum.m_points[0] = nc - up * near_height - right * near_width;t_frustum.m_points[1] = nc + up * near_height - right * near_width;t_frustum.m_points[2] = nc + up * near_height + right * near_width;t_frustum.m_points[3] = nc - up * near_height + right * near_width;t_frustum.m_points[4] = fc - up * far_height - right * far_width;t_frustum.m_points[5] = fc + up * far_height - right * far_width;t_frustum.m_points[6] = fc + up * far_height + right * far_width;t_frustum.m_points[7] = fc - up * far_height + right * far_width;}
}

我们使用摄像机的位置和方向，并利用简单的三角函数就能计算出每个视锥分块需要的8个顶点。

2.生成光空间的摄像机分块包围盒投影矩阵：

接下来非常重要的一步就是为每个眼视锥分块找到一个AABB，该包围盒位于光空间下。我们为每个AABB构造一个投影观察矩阵，之后的每一个眼视锥分块的深度图都是用与其对应的光空间投影观察矩阵生成。变换矩阵生成代码如下：

//生成光视锥切片的光空间变换矩阵(t_modelview为光观察矩阵)
void ShadowMap::generate_crop_matrices(const glm::mat4& t_modelview) {glm::mat4 t_projection;for(int i = 0 ; i < m_num_splits ; i++) {Frustum& t_frustum = m_frustums[i];glm::vec3 tmax(-1000.0f, -1000.0f, 0.0f);glm::vec3 tmin(1000.0f, 1000.0f, 0.0f);//寻找光空间八个顶点的最大最小z值glm::vec4 t_transf = t_modelview * glm::vec4(t_frustum.m_points[0], 1.0f);tmin.z = t_transf.z;tmax.z = t_transf.z;for(int j = 1 ; j < 8 ; j++) {t_transf = t_modelview * glm::vec4(t_frustum.m_points[j], 1.0f);if(t_transf.z > tmax.z) { tmax.z = t_transf.z; }if(t_transf.z < tmin.z) { tmin.z = t_transf.z; }}//可能在距离光源较近处存在遮挡物(不知道给啥值，给个大一点的吧)tmax.z += 10;//生成光的正交投影矩阵(长宽用单位1，每个视锥分片的光正交投影矩阵可以通过该单位宽高投影矩阵缩放平移后得到)glm::mat4 t_ortho = glm::ortho(-1.0f, 1.0f, -1.0f, 1.0f, -tmax.z, -tmin.z);glm::mat4 t_shad_mvp = t_ortho * t_modelview;//找到在光空间下视锥切片的(x,y)范围for(int j = 0 ; j < 8 ; j++) {t_transf = t_shad_mvp * glm::vec4(t_frustum.m_points[j], 1.0f);t_transf.x /= t_transf.w;t_transf.y /= t_transf.w;if(t_transf.x > tmax.x) { tmax.x = t_transf.x; }if(t_transf.x < tmin.x) { tmin.x = t_transf.x; }if(t_transf.y > tmax.y) { tmax.y = t_transf.y; }if(t_transf.y < tmin.y) { tmin.y = t_transf.y; }}glm::vec2 tscale(2.0f / (tmax.x - tmin.x), 2.0f / (tmax.y - tmin.y));glm::vec2 toffset(-0.5f * (tmax.x + tmin.x) * tscale.x, -0.5f * (tmax.y + tmin.y) * tscale.y);glm::mat4 t_shad_crop=glm::mat4(1.0);t_shad_crop[0][0] = tscale.x;t_shad_crop[1][1] = tscale.y;t_shad_crop[0][3] = toffset.x;t_shad_crop[1][3] = toffset.y;t_shad_crop = glm::transpose(t_shad_crop);//注意glm按列储存，实际矩阵要转置t_projection = t_shad_crop * t_ortho;//保存视锥切片的光投影矩阵m_projection_matrices[i] = t_projection;//保存世界坐标到光空间变换的矩阵(每个视锥切片保存一个自己的)m_crop_matrices[i] = t_projection * t_modelview;}
}

这里第一个需要注意的是tmax.z+=10。这是由于以下情况：

因此该值也可以定义一个设置接口，根据不同的场景，可随意调整该值。

在构造正交投影矩阵时，我们使用：

glm::mat4 t_ortho = glm::ortho(-1.0f, 1.0f, -1.0f, 1.0f, -tmax.z, -tmin.z);

创建了一个x,y∈[-1,1]，z∈[-tmax.z,-tmin.z]的正交投影矩阵。首先，使用x,y∈[-1,1]，是因为我们的每个分块矩阵都可以使用单位x,y缩放平移后获得。而我们设定z∈[-tmax.z,-tmin.z]，是因为眼空间指向-Z方向，而glm::ortho传入的是近平面和远平面(指向+Z)。

平移缩放矩阵可以根据构造正交投影矩阵的原理得出，公式如下：

更详细的可以看我介绍投影矩阵的博客：https://blog.csdn.net/qq_39300235/article/details/90670282

计算完所有的分块光空间变换矩阵后，我们将其保存。

3.更新分片位置(相机空间，标准化为[0,1])

我们依旧可以使用投影矩阵的特点，直接计算标准化后的z，并变换至[0,1]，详细介绍也在上述博客中。根据投影公式：

眼空间中w为1。所以我们更新分片在相机空间的位置代码如下：

//更新分片位置(相机空间，标准化为[0,1])
void ShadowMap::update_far_bounds(const glm::mat4& projection) {for(int i = 0 ; i < m_num_splits ; i++) {//将分片的位置标准化到(0,1)，需要注意我们初识给的视锥分片是正值，然而利用投影矩阵标准化需要在眼空间上执行，//眼空间到Z方向为指向负方向，所以我们使用视锥分片位置需要乘-1。而w分量正好是-z，所以不需要乘m_far_bounds[i] =0.5f*((-1.0f*m_frustums[i].far()*projection[2][2]+projection[3][2])/m_frustums[i].far())+0.5f;}
}

4.更新眼空间到光空间变换矩阵

为了方便起见，我们还直接保存从眼空间到光空间的变换矩阵：


//更新眼空间到光空间的变换矩阵
void ShadowMap::update_texture_matrices(const glm::mat4& projection, const glm::mat4& view_inverse) {for(int i = 0 ; i < m_num_splits ; i++) {//为了比对眼空间的深度和光空间的深度，需要m_bias将[-1，1]变换为[0，1]m_texture_matrices[i] = m_bias * m_crop_matrices[i] * view_inverse;}
}

为每个光视锥体渲染一个阴影贴图

完成所有矩阵的获取后，我们可以为每一个光视锥体分块渲染一个深度图。这里我们的fbo和texture构造代码如下：

void ShadowMap::create_fbo() {glGenFramebuffers(1, &m_fbo);glBindFramebuffer(GL_FRAMEBUFFER, m_fbo);glDrawBuffer(GL_NONE);glReadBuffer(GL_NONE);glBindFramebuffer(GL_FRAMEBUFFER, 0);
}void ShadowMap::create_texture() {if(m_texture_array) {glDeleteTextures(1, &m_texture_array);}glGenTextures(1, &m_texture_array);glBindTexture(GL_TEXTURE_2D_ARRAY, m_texture_array);glTexImage3D(GL_TEXTURE_2D_ARRAY, 0, GL_DEPTH_COMPONENT, depth_tex_size(), depth_tex_size(),CSM_MAX_SPLITS, 0, GL_DEPTH_COMPONENT, GL_FLOAT, nullptr);glTexParameteri(GL_TEXTURE_2D_ARRAY, GL_TEXTURE_MIN_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D_ARRAY, GL_TEXTURE_MAG_FILTER, GL_LINEAR);glTexParameteri(GL_TEXTURE_2D_ARRAY, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);glTexParameteri(GL_TEXTURE_2D_ARRAY, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);glBindTexture(GL_TEXTURE_2D_ARRAY, 0);
}

写入深度的着色器很简单，我们不需要做其他的操作，只需要把世界空间转换为光视锥分块空间，代码如下：

顶点着色器：

#version 410 corelayout (location = 0) in vec3 aPos;
layout (location = 1) in vec2 aTexCoord;uniform mat4 lightSpaceMatrix;
uniform mat4 view;
uniform mat4 model;void main(){gl_Position=lightSpaceMatrix*view*model*vec4(aPos,1.0f);
}

我们为每个光包围盒投影矩阵传入各自的lightSpaceMatrixd。

片元着色器：

#version 410 coreout vec4 FragColor;void main(){}

完成后，我们可以可视化我们写入的深度图，效果如下：

渲染整个场景的阴影

写入每个分块的深度图后，我们可以为整个场景渲染阴影。我们判断每个摄像机内的片元，根据其位于哪个分块内，将其转换为对应光空间，并与深度图的深度比对，如果大于深度图的深度，则处于阴影。我们的着色器代码如下：

顶点着色器：

#version 410 core
layout (location = 0) in vec3 position;
layout (location = 1) in vec2 aTexCoord;out vec2 texCoord;
out vec4 viewPos;
out vec4 farBounds;
out mat4 FragPosLightSpace[4];uniform mat4 model;layout (std140) uniform Matrices{vec4 afarBounds;mat4 projection;mat4 view;mat4 aFragPosLightSpace[4];
};void main()
{gl_Position=projection*view*model*vec4(position, 1.0f);texCoord=aTexCoord;farBounds=afarBounds;for (int i=0; i<4; i++) {FragPosLightSpace[i]=aFragPosLightSpace[i];}viewPos=view*model*vec4(position, 1.0f);
}

片元着色器：

#version 410 corein vec2 texCoord;
in vec4 viewPos;
in vec4 farBounds;
in mat4 FragPosLightSpace[4];out vec4 FragColor;uniform sampler2DArray shadowMap;
uniform vec2 texSize;float shadowCal(vec4 viewPos);
void shadowJud();
void main()
{float shadow=shadowCal(viewPos);vec3 lighting=(1.0f-shadow)*vec3(0.5,0.5,0.5);FragColor = vec4(lighting,1.0f);
}void shadowJud(){}float shadowCal(vec4 viewPos){int aIndex=3;if (gl_FragCoord.z<farBounds.x) {aIndex=0;}else if(gl_FragCoord.z<farBounds.y) {aIndex=1;}else if(gl_FragCoord.z<farBounds.z) {aIndex=2;}vec4 projCoords=FragPosLightSpace[aIndex]*viewPos;float closestDepth=texture(shadowMap,vec3(projCoords.xy,float(aIndex))).r;float currentDepth=projCoords.z;float bias=0.001;//阴影偏移，可以使所有片元在表面之上，消除失真float shadow=0.0f;
//    shadow=currentDepth-bias>closestDepth?1:0;vec2 texelSize=1.0/texSize;for(int x=-1;x<=1;x++){for(int y=-1;y<=1;y++){float pcfDepth=texture(shadowMap,vec3(projCoords.xy+vec2(x,y)*texelSize,float(aIndex))).r;shadow+=currentDepth-bias>pcfDepth?0.7f:0.0f;}}shadow/=9.0f;if(projCoords.z>1.0)//超出视界的z都大于1，大于深度贴图中的1.0，所以返回黑色，因此我们要加限制shadow=0.0f;return shadow;
}

我们这里做了简单的3*3的PCF，完成效果如下图：

我们可视化摄像机空间的分块效果，如下图：

我们实验中将视锥分为4片，而市面上许多游戏分三层就能达到非常好的阴影效果了。