方差阴影贴图与切比雪夫不等式

Tags:

要理解方差阴影贴图的来龙去脉,必须先深刻理解概率论中的几个概念:

  • 矩(Moment)
  • 数学期望(Mean)
  • 方差(Variance)
  • 马可夫不等式 (Markov's Inequality)
  • 切比雪夫不等式 (Chebyshev's inequality)
  • 切比雪夫不等式的one-tailed版本 (one-tailed version of Chebyshev's inequality)

矩(Moment)

https://en.wikipedia.org/wiki/Moment_(mathematics)

给定关于实变量x、常数c的实值连续函数f(x),它的n阶矩(n-th moment)的公式是:

\[ \mu ^{n} = \int _{-\infty }^{ +\infty } (x - c)^{n} f(x) dx \]

数学期望(Mean)

当c = 0,n = 1时,上述公式变成:

\[ \mu = \int _{-\infty }^{ +\infty } x f(x) dx \]

这也就是数学期望(Mean)的积分公式。

方差(Variance)

当\(c = \mu\)时,n阶矩可称为n阶中心矩;当\(c = \mu,n = 2\)时,2阶中心矩的公式为:

\[ \mu ^{2} = \int _{-\infty }^{ +\infty } (x - \mu )^{ 2 } f(x) dx \]

这其实就是方差(Variance)的积分公式。下面作简单推导。

方差的定义式为:

\[ Var(X) = E[(X - \mu)^{2}] = E[(X - E[X])^{2}] = \sigma ^{2} \]

可以推出:

\[ Var(X) = E[X^{2} - 2 X E[X] + E[X]^{2} ] \]

\[ = E[X^{2}] - 2 E[X] E[X] + E[X]^{2} \]

\[ = E[X^{2}] - E[X]^{2} \]

而2阶中心矩公式可以推出:

\[ \int (x - \mu )^{ 2 } f(x) dx \]

\[ = \int x^{ 2 } f(x) dx - 2\mu \int x f(x) dx + \int \mu ^{ 2 } f(x) dx \]

\[ = \int x^{ 2 } f(x) dx - 2\mu \cdot \mu + \mu ^{ 2 } \]

\[ = \int x^{ 2 } f(x) dx - \mu ^{ 2 } \]

\[ =E[X^{ 2 }] - E[X] ^{ 2 } \]

马可夫不等式 (Markov's Inequality)

设X是非负的随机变量,且有a > 0,那么X大于等于a的概率不超过X的数学期望除以a:

\[ P _{ X \geq a } \leq \frac { E[X] }{ a } \]

(Note:这里的P是指概率)

证明:

证明前需要先理解一个概念:示性函数(Indicator)。对于任意事件e,当e发生时,\( I _{e} = 1\), 当E没发生时,\( I _{e} = 0\)。

那么把\( X \geq a \)当作一个事件e,当e发生时,有:

\[ I _{ X \geq a } = 1 \]

\[ a I _{ X \geq a } \leq X \]

两边同时变成数学期望,不等式依然成立:

\[ E[a I _{ X \geq a }] \leq E[X] \]

又因为数学期望的线性关系,有:

\[ E[a I _{ X \geq a }] = a \cdot E[I _{ X \geq a }] \]

又因为函数\( I _{ X \geq a } \)的取值只有2种,所以可直接得到:

\[ a \cdot E[I _{ X \geq a }] = a \cdot ( 1\cdot P _{ X \geq a } + 0\cdot P _{ X \lt a } ) \]

\[ = a \cdot P _{ X \geq a } \]

综上,就得到了:

\[ a \cdot P _{ X \geq a } \leq E[X] \]

\[ P _{ X \geq a } \leq \frac { E[X] }{ a } \]

切比雪夫不等式 (Chebyshev's inequality)

设有随机变量X以及它的数学期望\(\mu \) 、有限且不等于0的方差\( \sigma ^{2} \),对于任意>0的实数k,以下不等式成立:

\[ P _{ | X - \mu | \geq k\sigma } \leq \frac { 1 }{ k^{2} } \]

(Note:这里的P是指概率)

这就是切比雪夫不等式。其中,因为概率P永远小于等于1,所以k值要大于1这个不等式才有意义。

证明:

设有随机变量\( Y = (X - \mu )^{2} \) 以及 \( a = (k\sigma )^{2} \),代入马可夫不等式后:

\[ P _{ Y \geq a } \leq \frac { E[Y] }{ a } \]

\[ P _{ (X - \mu )^{2} \geq (k\sigma )^{2} } \leq \frac { E[(X - \mu )^{2}] }{ (k\sigma )^{2} } \]

回顾下方差公式:

\[ Var(X) = E[(X - \mu)^{2}] = \sigma ^{2} \]

显然有:

\[ P _{ (X - \mu )^{2} \geq (k\sigma )^{2} } \leq \frac { E[(X - \mu )^{2}] }{ (k\sigma )^{2} } = \frac { \sigma ^{2} }{ (k\sigma )^{2} } = \frac { 1 }{ k^{2} } \]

左边的式子可以进一步简化:

\[ (X - \mu )^{2} \geq (k\sigma )^{2} \]

\[ |X - \mu | \geq k\sigma \]

(右边没有绝对值是因为有前提条件k>0,即使标准差\(\sigma < 0 \)该等式依然成立 )

于是切比雪夫不等式成立:

\[ P _{ | X - \mu | \geq k\sigma } \leq \frac { 1 }{ k^{2} } \]

切比雪夫不等式的one-tailed版本

切比雪夫不等式的one-tailed版本其实就是坎泰利不等式Cantelli's inequality。坎泰利不等式公式如下:

1.png

(from wiki,Pr等价于上文的P)

而切比雪夫不等式的one-tailed版本如下:

\[ P _{ X - \mu \geq t } \leq \frac { \sigma ^{2} }{ \sigma ^{2} + t ^{2} } ,t > 0 \]

一模一样的。

证明:

要证明one-tailed公式,要用到马可夫不等式。首先定义\(Y = X - \mu \),那么就有\(E[Y] = E[X - \mu ] = E[X]- E[\mu ] = \mu - \mu = 0\),以及:\( Var[Y] = Var[X - \mu ] = Var[X] - Var[\mu ] = Var[X] - 0 = Var[X] \)。

于是:

\[ P _{ Y \geq t } = P _{ Y + \mu \geq t + \mu } \leq P _{ (Y + \mu)^{2} \geq (t + \mu)^{2} } \]

这时候用上马可夫不等式 \( P _{ X \geq a } \leq \frac { E[X] }{ a } \),得到:

\[ P _{ (Y + \mu)^{2} \geq (t + \mu)^{2} } \leq \frac { E[(Y + \mu)^{2}] }{ (t + \mu)^{2} } \]

再用上方差公式:

\[ E[(Y + \mu)^{2}] = Var[Y + \mu] + E[Y + \mu]^{2} \]

\[ = Var[Y] + Var[\mu] + E[Y + \mu]^{2} \]

\[ = \sigma ^{2} + 0 + (E[Y] + E[\mu])^{2} \]

\[ = \sigma ^{2} + (0 + E[\mu])^{2} \]

\[ = \sigma ^{2} + \mu ^{2} \]

所以有:

\[ P _{ Y \geq t } \leq \frac { \sigma ^{2} + \mu ^{2} }{ (t + \mu)^{2} } \]

接着令 \( \phi(\mu ) = \frac { \sigma ^{2} + \mu ^{2} }{ (t + \mu)^{2} } \),求导\( \phi '(\mu ) = 0\)时的\( \mu \)值。这个导数算起来比较复杂,我找了个在线导数计算工具来辅助下(这不是广告)。

先进入http://zh.numberempire.com/derivativecalculator.php,输入: (a^2+x^2)/((b+x)^2),得到导数公式: (2*b*x-2*a^2)/(x^3+3*b*x^2+3*b^2*x+b^3)。

a就是\( \sigma \),b就是t,x就是\(\mu \),把这个式子弄成latex:

\[ \phi '(\mu ) = \frac { 2bx - 2a^{2} }{ x^{3} + 3bx^{2} + 3b^{2}x + b^{3} } = 0 \]

这个方程也是复杂,继续用工具来算就好了。

进入http://zh.numberempire.com/equationsolver.php

输入:(2*b*x-2*a^2)/(x^3+3*b*x^2+3*b^2*x+b^3) = 0,

得到:\( x = \frac { a^{2} } { b } \) ,即:

\[ \mu = \frac { \sigma ^{2} } { t } \]

也就是说当\( \mu = \frac { \sigma ^{2} } { t } \)时,\( \phi(\mu ) \)取得最小值。而又因为对任意的\( \mu \),概率\( P _{ Y \geq t } \) 都不超过\( \phi(\mu ) \),所以原先的不等式可以进一步简化成:

\[ P _{ Y \geq t } \leq \frac { \sigma ^{2} + \mu_{ * } ^{2} }{ (t + \mu _{ * } )^{2} } = \frac { \sigma ^{2} + (\frac { \sigma ^{2} } { t }) ^{2} }{ (t + \frac { \sigma ^{2} } { t } )^{2} } \]

右边的式子继续简化;

\[ \frac { \sigma ^{2} + (\frac { \sigma ^{2} } { t }) ^{2} }{ (t + \frac { \sigma ^{2} } { t } )^{2} } = \frac { \frac { \sigma ^{2}t^{2} +\sigma ^{4} } { t^{2} } } { \frac { (t^{2} + \sigma ^{2})^{2} } { t^{2} } } = \frac { \sigma ^{2}t^{2} +\sigma ^{4} } { (t^{2} + \sigma ^{2})^{2} } = \frac { \sigma ^{2} } { t^{2} + \sigma ^{2} } \]

所以:

\[ P _{ Y \geq t } = P _{ X - \mu \geq t } \leq \frac { \sigma ^{2} } { t^{2} + \sigma ^{2} } \]

得证。

方差阴影贴图(VSM)

生成VSM,相比SSM(Standard Shadow Map),除了把深度d存到深度图,还要多存一个d*d。似乎看起来有点蠢,明明有d了,要用d平方不就是运算的时候d*d就行了嘛。其实VSM的原理还得结合一些硬件知识来理解。相比SSM是单通道贴图,VSM则是双通道贴图,并且VSM在用于光照计算前,可以做纹理过滤或抗锯齿处理,纹理过滤,相当于纹理被模糊化(VSM能实现软阴影的关键之处)。

SSM并不能被模糊,而VSM能模糊化的好处就先不说了。模糊化后,d和d平方的值就在一定范围内做了偏移。那么我们可以设:

  • 模糊前的深度d和d平方分别为x和\( x^{2} \)

  • 模糊后的深度d和d平方分别为一阶矩M1(x的数学期望)和一阶矩M2(\( x^{2} \)的数学期望)

再用上前面给出的方差公式:

\[ Var(X) = E[X^{2}] - E[X]^{2} \]

所以对于模糊前的x(真实深度),我们可以求出它的数学期望和方差:

\[ \mu = E(x) = M1 \]

\[ \sigma ^{2} = Var(x) = E[x^{2}] - E[x]^{2} = M2 - M1^{2} \]

到了这里,VSM的原理已经相当清楚了:先是存了原始的深度x和x平方,然后又对x和x平方做了一些坏事(对我们来说是好事):把x和x平方稍微搅浑了,搅浑了后想恢复原始的x是不可能了,但是我们可以利用搅浑后的x和x平方(一、二阶矩),算出x的方差,再结合搅浑后的x(x的数学期望,一阶矩),就能对原始的x做概率估计。

(深度值越大就越深,越接近1,即越白,所以靠近摄像机的深度偏白色,远处的偏黑色)

但还没完,VSM还和切比雪夫不等式的one-tailed版本有关系。首先拿出该公式:

\[ P _{ X - \mu \geq t } \leq \frac { \sigma ^{2} }{ \sigma ^{2} + t ^{2} } ,t > 0 \]

设\( t' = t - \mu ,t > \mu > 0 \),代入上式:

\[ P _{ X - \mu \geq t' } \leq \frac { \sigma ^{2} }{ \sigma ^{2} + t' ^{2} } \]

\[ P _{ X - \mu \geq t - \mu } \leq \frac { \sigma ^{2} }{ \sigma ^{2} + (t - \mu) ^{2} } \]

\[ P _{ X \geq t } \leq \frac { \sigma ^{2} }{ \sigma ^{2} + (t - \mu) ^{2} } \]

(在Variance Shadow Maps 论文中就会看到这个不等式)

下面的就按照论文的思路讲。假定一个3d场景里有一块挡板(Occluder),产生了阴影(VSM),然后在fragment shader里对某一个fragment做阴影计算,这时候fragment可能被遮挡,也可能没被遮挡(薛定谔的猫?)。

首先把fragPos转换到light space下,从而取出过过滤后的VSM里的深度,称为\(d_{1}\)。而fragPos当前在light place的深度则称为\(d_{2} \)。 \( P _{ X \geq t } \)里的X则是VSM里未过滤前的深度。

把t当成\(d_{2} \),再加上用前面的公式算出来的数学期望\(\mu \)和方差\(\sigma ^{2} \),就可以确定,任意X超过\(d_{2} \)的概率,不大于\( \frac { \sigma ^{2} }{ \sigma ^{2} + (t - \mu) ^{2} } \)。不过,这个只是我们需要的概率的上界, \( P _{ X \geq t } \)并不能被准确算出来。

再设这个表示frag没有被挡板遮挡的实际概率\( P_{ X \geq t } \)为p,那么frag的d(即X)的数学期望为:

\[ \mu = E(x) = p d_{2} + (1 - p)d_{1} \]

(这里可以这样理解:当p=1时,期望值为\(d_{2}\),即fragPosLightSpace.z,表示frag在light space深度和在VSM的深度完全无关,即完全没有被遮挡,0%的阴影;当p=0时,期望值为\(d_{1}\),即在light space深度绝对等于在VSM的深度,表示被完全遮挡,100%的阴影)

\(d^{2}\)的数学期望为:

\[ E(x^{2}) = p d_{2}^{2} + (1 - p)d_{1}^{2} \]

d的方差为:

\[ \sigma ^{2} = E[x^{2}] - E[x]^{2} = p d_{2}^{2} + (1 - p)d_{1}^{2} - (p d_{2} + (1 - p)d_{1} )^{2} \]

\[ = p d_{2}^{2} + (1 - p)d_{1}^{2} - p^{2} d_{2}^{2} - 2p(1 - p)d_{1} d_{2} - (1 - p)^{2} d_{1} ^{2} \]

\[ = (p - p^{2}) d_{2}^{2} + (1 - p)(d_{1}^{2} - 2pd_{1} d_{2} - (1 - p)d_{1}^{2} ) \]

\[ = (p - p^{2}) d_{2}^{2} + (1 - p)(d_{1}^{2} - 2pd_{1} d_{2} - d_{1}^{2} + pd_{1}^{2} ) \]

\[ = (p - p^{2}) d_{2}^{2} + (1 - p)( - 2pd_{1} d_{2} + pd_{1}^{2} ) \]

\[ = (p - p^{2}) d_{2}^{2} + (p - p^{2})( - 2d_{1} d_{2} + d_{1}^{2} ) \]

\[ = (p - p^{2})(d_{2}^{2} - 2d_{1} d_{2} + d_{1}^{2} ) \]

\[ = (p - p^{2})(d_{2} - d_{1})^{2} \]

把这些参数都代入到上面的不等式,求出\(p_{max}(d_{2}) \):

\[ P _{max}(d_{2}) = \frac { \sigma ^{2} }{ \sigma ^{2} + (d_{2} - \mu) ^{2} } \]

\[ = \frac { (p - p^{2})(d_{2} - d_{1})^{2} }{ (p - p^{2})(d_{2} - d_{1})^{2} + ( p d_{2} + (1 - p)d_{1} - d_{2} ) ^{2} } \]

\[ = \frac { (p - p^{2})(d_{2} - d_{1})^{2} }{ (p - p^{2})(d_{2} - d_{1})^{2} +(1 - p) ^{2} (d_{2} - d_{1} ) ^{2} } \]

\[ = \frac { p - p^{2} }{ (p - p^{2}) +(1 - p) ^{2} } \]

\[ = \frac { p - p^{2} }{ (p - p^{2}) +(1 - p) ^{2} } \]

\[ = \frac { p - p^{2} }{ 1 - p } \]

\[ = p \]

也就是说,\( P_{ max }( X \geq d_{2}) \) 这个最大值就刚好是我们想要求的\( P( X \geq d_{2}) \) !

\[ P _{ X \geq t } \leq \frac { \sigma ^{2} }{ \sigma ^{2} + (t - \mu) ^{2} } \]

这里的\( \leq \) 刚好是\( = \),不等式变成了等式。Amazing!

理论到此为止,剩下的就是工程问题了。

工程实现

如果你的工程里已经实现了Standard Shadow Map,那么只需要几个步骤即可实现基本的VSM:

  • 把Shadow Map Buffer改为输出depth和depth平方到颜色buffer。可能需要先修改RT的创建代码,下面以OpenGL 4.x为例:
// 32位高精度双通道纹理 用来存moment1和moment2
glTexImage2D(GL_TEXTURE_2D, 0, GL_RG32F, width, height, 0, GL_RG, GL_FLOAT, NULL);

然后修改shadow map着色器:

vs:

#version 410 core

layout (location = 0) in vec3 position;

uniform mat4 lightPV;
uniform mat4 model;

void main()
{
    gl_Position = lightPV * model * vec4(position, 1.0);
}

fs:

#version 410 core

in vec4 v_position;
out vec2 outColor;

void main()
{             
    float depth = v_position.z / v_position.w;
    depth = depth * 0.5 + 0.5;// TO NDC [0, 1]

    float moment1 = depth; // 一阶矩
    float moment2 = depth * depth; // 二阶矩

    float dx = dFdx(depth);
    float dy = dFdy(depth);
    moment2 += 0.25 * (dx * dx + dy * dy); // 解决acne问题
    outColor = vec2(moment1, moment2);
}

可以渲染下Shadow Map是否正常:

rg通道一起输出:

2.png

只有r通道(moment1):

3.png

只有g通道(moment2):

4.png

然后是应用的问题,需要修改lighting shader里的阴影计算代码:

// 只贴下核心代码:
uniform sampler2D shadowMap;

// 这个函数是用来算出上文说的Pmax
float chebyshevUpperBound(sampler2D shadowMap, float d, vec2 coord)
{
    vec2 moments = texture(shadowMap, coord).rg;
    // Surface is fully lit. as the current fragment is before the light occluder
    if (d <= moments.x)
        return 1.0;

    // The fragment is either in shadow or penumbra. We now use chebyshev's upperBound to check
    // How likely this pixel is to be lit (p_max)
    float variance = moments.y - (moments.x * moments.x);
    //variance = max(variance, 0.000002);
    variance = max(variance, 0.00002);

    float d_minus_mean = d - moments.x;
    float p_max = variance / (variance + d_minus_mean * d_minus_mean);

    return p_max;
}

// 返回阴影百分比[0,1], 然后拿去乘以光照颜色即可
float ShadowCalculation_Dir(vec3 fragPos, Light light) {

    vec4 fragPosLightSpace = light.lightPV * vec4(fragPos, 1.0);
    // perform perspective divide
    vec3 projCoords = fragPosLightSpace.xyz / fragPosLightSpace.w;
    // transform to [0,1] range
    projCoords = projCoords * 0.5 + 0.5;
    float currentDepth = projCoords.z;
    shadow = 1 - chebyshevUpperBound(shadowMap, currentDepth, projCoords.xy);
    return shadow;
}

最终效果:

5.png

这里可以看到阴影边缘是变模糊了的,不过还是有锯齿。可以进一步对shadow map做blur来实现软化。

(未经授权禁止转载)
Written on November 17, 2017

博主将十分感谢对本文章的任意金额的打赏^_^