图形学论文阅读 - 使用盒型滤波图像金字塔实时近似大卷积核

778 阅读4分钟

本文内容皆引用于知乎用户 Monica的小甜甜

The Power of Box Filters: Real-time Approximation to Large Convolution Kernel by Box-filtered Image Pyramid

SIGGRAPH Asia 2019 Technical Briefs

介绍

在空间范围内进行滤波本质是对一个着色点与他周围着色点进行一些处理,最关键的就是如何获取该着色点周围着色点的信息。和传统的高斯滤波直接对周围点进行采样不同,本文提出了一个新方法,即从 mipmapmipmap 里获取周围点信息用于滤波,显然采样率就大大减少了。

例如,对于一个 1024×10241024\times1024 分辨率的纹理进行采样,mipmapmipmap 的层级仅有 1010 级,即需要 1010 次采样。而 mipmapmipmap 的生成可以直接通过硬件加速生成:在一些图形 API 中,图形驱动程序已经实现了盒型滤波器的 mipmapmipmap 生成,并做了专门的优化,因此盒型滤波器的实现成本很低。

方法

采样公式

高斯滤波采样公式

对于一个着色点 xx 想要进行高斯滤波,首先要取该着色点周围 R×RR\times R 方形区域里面的像素然后在乘以高斯滤波函数。

p(x)=R×RG(x,y)t(x,y)dxdy(1)p(x)=\int_{R \times R} G(x, y) t(x, y) d x d y \quad (1)

G(x,y)G(x,y) :二维的高斯滤波函数。由于高斯滤波有个很好的特性:可以先在 xx 方向进行滤波,然后再在 yy 方向滤波。因此实际过程中通常用于一维的高斯滤波函数。

t(x,y)t(x,y) :对着色点进行采样的函数。采样时只是单纯的后处理,即对一张纹理进行处理,因此没有z方向。

盒型滤波采样公式

p(x)=Lmw(L)pdown (L)Lmw(L)(2)p(x)=\frac{\sum_{L}^{m} w(L) p_{\text {down }}(L)}{\sum_{L}^{m} w(L)} \quad (2)

pdown(L)p_{down} (L) :着色点 xx 在第 LL 层(当前层)的 mipmapmipmap 纹理颜色,该纹理采样是对等待滤波的纹理进行采样。

w(L)w(L) :盒型滤波的权重,该值是根据高斯滤波所推导得出,在后面介绍如何推导。

mmmipmapmipmap 的最大层数理论上为无穷大。

论文中对上述公式进行了变形,给出了递归求法:

p(x,L)={(1α(L))p(x,L+1)+α(L)pdown (x,L)LLmax p(x,L)=pdown (x,L)L=Lmax (3)p(x, L)=\left\{\begin{array}{ll} (1-\alpha(L)) p(x, L+1)+\alpha(L) p_{\text {down }}(x, L) & L \neq L_{\text {max }} \\ p(x, L)=p_{\text {down }}(x, L) & L=L_{\text {max }} \end{array}\right. \quad (3)
α(L)=w(L)lw(L)dL\alpha(L)=\frac{w(L)}{\int_{l}^{\infty} w(L)dL}

p(x,L)p(x,L) :着色点 xx 在第 LL 层的纹理颜色。

α(L)\alpha(L) :每层采样过程的混合权重 。

image.png

权重推导

对于一个着色点 pp ,对其应用盒型滤波的结果应该与高斯滤波是相同的。但对于盒型滤波,由于 mipmapmipmap 的金字塔式生成,不同层级的贡献肯定是不同的。

p(x)=R×RG(x,y)t(x,y)dxdy=Lw(l)tbox(l)dl(4)p(x)=\int_{R \times R} G(x, y) t(x, y) d x d y =\int_{L}^{\infty} w(l) t_{b o x}(l) d l \quad (4)
tbox(l)=LB(l,x,y)t(x,y)dxdy(5)t_{b o x}(l)=\int_{L}^{\infty} B(l, x, y) t(x, y) d x d y \quad (5)

tbox(l)t_{box}(l) 为采样第 llmipmapmipmap 的加权采样函数。

B(l,x,y)B(l,x,y) :第 llmipmapmipmap 的权重。

又因为 mipmapmipmap 由每层的四个像素求平均值以生成下一层的像素,因此每层的权重都是上一层的四分之一。然后将边界条件转为为一个同时包含 x,y,lx,y,l 的边界条件,虽然该边界条件更加准确的做法应该根据式 (6)(6) 的边界条件所围成的面积来计算,但为了方便后续的计算贴合一维高斯滤波函数,可以近似成式 (7)(7) 的形势,该近似对结果的影响非常小:

B(l,x,y)={1N=4l2x2l&&2y2l0 otherwise (6)B(l, x, y)=\left\{\begin{array}{ll} \frac{1}{N}=4^{-l} & 2|x| \leq 2^{l} \& \& 2|y| \leq 2^{l} \\ 0 & \text { otherwise } \end{array}\right. \quad (6)
x2+y24lπ(7)x^{2}+y^{2} \leq \frac{4^{l}}{\pi} \quad (7)

NN :盒子的尺寸,由 mipmapmipmap 的层级 ll 所决定的,即 N=2l2lN=2^l * 2^l

由上式可以看出采样点的坐标 (x,y)(x,y) 与采样的层级 ll 存在的不等式关系。

将以上公式代回式 (4)(4) 的盒型滤波采样公式:

p(L)=R×Rt(x,y)Lw(l)4ldldxdy(8)p(L)=\int_{R \times R} t(x, y) \int_{L}^{\infty} w(l) 4^{-l} d l d x d y \quad (8)

和高斯滤波采样公式相比,显然可以看出盒型滤波函数和高斯滤波函数的关系。再将式 (7)(7) 求得的边界条件带入界,即可得到和 mipmapmipmap 的层级 LL 有关的高斯滤波函数,同时可以看出对边界条件的近似的目的就是为了方便高斯函数的计算:

G(x,y)=G(L)=Lw(l)4ldl(9)G(x,y)=12πσ2e(x2+y22σ2)G(L)=12πσ2e(4L2πσ2)\begin{array}{c} G(x, y)=G(L)=\int_{L}^{\infty} w(l) 4^{-l} d l \quad (9)\\ G(x, y)=\frac{1}{2 \pi \sigma^{2}} e^{\left(-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right)} \\ G(L)=\frac{1}{2 \pi \sigma^{2}} e^{\left(-\frac{4^{L}}{2 \pi \sigma^{2}}\right)} \end{array}

假设 w(l)4lw(l)4^{-l} 的原函数为 g(l)g(l) ,即 g(l)=w(l)4lg'(l)=w(l)4^{-l} ,由牛顿-莱布尼兹公式和式 (9)(9) 可得 :

G(L)=Lw(l)4ldl=limlg(l)g(L)(10)G(L)=\int_{L}^{\infty} w(l) 4^{-l} d l =\lim _{l \rightarrow \infty} g(l)-g(L) \quad (10)

由于层级越大,盒型滤波的权重越小。故式 (10)(10) 的极限结果应为 00。因此可得高斯滤波函数与盒型滤波函数的关系为:

G(L)=g(L)(11)G(L)=-g(L) \quad (11)

g(l)=w(l)4lg'(l)=w(l)4^{-l} 或对式 (9)(9) 求导皆可得出权重 w(l)w(l) 与层级 ll 的关系式及其原函数,再由权重 w(l)w(l) 及其原函数可以得出第 LL 层的混合权重 α(L)\alpha(L)

w(l)=4lg(l)=4lG(l)=16lln44π2σ4e(4l2πσ2)\begin{aligned} w(l) &=4^{l} g^{\prime}(l) \\ &=-4^{l} G^{\prime}(l) \\ &=\frac{16^{l} \ln 4}{4 \pi^{2} \sigma^{4}} e^{\left(-\frac{4^{l}}{2 \pi \sigma^{2}}\right)} \end{aligned}
W(l)=12πσ2e(4l2πσ2)W(l)=\frac{1}{2 \pi \sigma^{2}}e ^ {\left(-\frac{4^{l}}{2 \pi \sigma^{2}}\right)}
α(L)=w(L)lw(L)dL=16Lln42πσ2(4L+2πσ2)\alpha(L) =\frac{w(L)}{\int_{l}^{\infty} w(L)dL} =\frac{16^{L} \ln 4}{2 \pi \sigma^{2}\left(4^{L}+2 \pi \sigma^{2}\right)}

代码实现

GLSL代码实现

#version 430 core
#define PI 3.141592654
in  vec2 v2f_TexCoords;
out vec4 Color_;

uniform sampler2D u_SourceTexture;
uniform sampler2D u_CoarserTexture;
uniform vec2	g_focus = vec2(0);
uniform	float	g_sigma = 5.0;
uniform int	g_level = 0;

求混合参数:float MipGaussianBlendWeight(vec2 tex)

  1. 首先依据纹理坐标生成高斯函数的标准差 σ\sigma 和方差 σ2\sigma^2
float sigma = g_sigma;
if (uint(g_focus.x) != 0xffffffff)
{
    const vec2 r = (2.0 * tex - 1.0) - g_focus; //(0,1)^2 to (-1,1)^2
    sigma *= dot(r, r);//length of vector r
}
const float sigma2 = sigma * sigma;
  1. α(L)\alpha (L) 的公式先求出分子部分: 16Lln4{16^{L} \ln 4}
const float numerator = (1 << (g_level << 2)) * log(4.0);
  1. 求出 α(L)\alpha (L) 公式的分母部分: 2πσ2(4L+2πσ2)2 \pi \sigma^{2}\left(4^{L}+2 \pi \sigma^{2}\right)
const float c = 2.0 * PI * sigma2;
const float denominator = c * ((1 << (g_level << 1)) + c);
  1. 返回混合权重 α(L)\alpha (L) ,并确保 α(L)\alpha (L) 的范围在 (0,1)(0,1)
return clamp(numerator / denominator,0,1);

主函数:void main()

  1. 只生成到第 10 层,到第10层后结束递归。
if(g_level == 10)
{
    Color_ = textureLod(u_CoarserTexture, v2f_TexCoords,10).rgba;
    return ;
}
  1. 获取第 LL 层的混合权重 α(L)\alpha (L)
float alpha = MipGaussianBlendWeight(v2f_TexCoords);
  1. 获取当前层和下一层的 mipmapmipmap
vec3 color_source = textureLod(u_SourceTexture, v2f_TexCoords,g_level + 1).rgb;
vec3 color_coarser = textureLod(u_CoarserTexture, v2f_TexCoords,g_level).rgb;
  1. 由公式 (1α(L))p(x,L+1)+α(L)pdown (x,L)(1-\alpha(L)) p(x, L+1)+\alpha(L) p_{\text {down }}(x, L) 求出纹理颜色值。
Color_ = vec4((1 - alpha) * color_source + alpha * color_coarser, 1.0f);