本篇文章谢绝转载，也禁止用于任何商业目的。文章中的源码实现的下载地址参见https://github.com/twinklingstar20/hdr_demo

1. 色阶重建

色阶重建（Tone Mapping）是什么样的一个技术呢，它与高动态光照渲染（High Dynamic Range, HDR）又有什么关系呢？^[1]，且听我娓娓道来。

假设你是一名摄影爱好者，一定听过HDR这个词，它的专业称法为高动态范围成像（High Dynamic Range Image, HDRI）^{[2, 3]}，是用来实现比普通图像（每个颜色通道占8位）表示的曝光范围更大的一种技术，高动态范围成像的目的就是要正确地表示真实世界中从太阳光直射到最暗的阴影这样大的范围亮度。举个例子来解释，大晴天出门去摄（zhuang）影（bi），拍摄了图1所示的一张图，一定会被人喷不专业，因为你这是想证明你对自然景色更感兴趣呢，还是对中国的文化遗产更感兴趣呢？照顾了天空的曝光，就使得暗部欠缺曝光。

图1. 想更准确的拍摄天空

你为了证明你是爱中国的文化遗产，你想重点关注建筑物，于是你拍摄了另外一张图，如图2所示，确实证明了你的爱，但是背景又是什么鬼，那是霾吗，大北京这么难能可贵的大晴天，却拍出了雾霾的效果。

分析下造成这个问题的原因，简单来说，照相机拍摄的亮度（Luminance）范围是有限的，如果想要亮部更加清晰，特别暗的区域就显得比较黑；相反，如果想要暗部更加清晰，特别亮的区域就会被截断（Clamp to white）成白色，显得特别的亮。我们想要的结果如图3所示，不仅能清楚的显现出建筑的模样，同时还有一个比较清晰的背景，这就是色阶重建技术达到的效果，它把一个高动态范围的亮度，变换为低动态范围（Low Dynamic Range）的亮度，使仅用低动态范围的亮度能清晰、有效地还原场景的效果，PERFECT！！！

图3. 高动态范围成像技术达到的效果图

再举个例子，如图4所示^[2]，我们可以用不同的曝光度，拍摄一组图片，由暗到亮，然后采用色阶重建技术把它变换为一组明暗协调的图片，如图5所示。色阶重建技术，又可以细分为全局色阶重建和局部色阶重建，想要了解各种色阶重建技术可以参见Reinhard etc^[4]。

图4. 不同曝光度的一组图像^[2]

图5. 全局和局部色阶重建后的图像^[2]

至此，介绍了高动态范围成像（HDRI）和色阶重建之间的关系，可以简单把色阶重建技术划分为全局色阶重建和局部色阶重建。接下来，又提出一个问题，在实时渲染中，特别在多光源的光照计算中，为什么存在无法渲染高动态范围的问题呢？

在计算机的颜色缓冲区中^[5]，存储的颜色范围是$[0,1]$，当多光源的光照计算，由于亮度的叠加，容易导致亮度超过1.0，此时，硬件设备（GPU）会把它截断至亮度1.0，大部分情况下做这样一种截断处理是没有问题的。但是如果一个场景中绝大部分的片断（像素）经过多光源的光照叠加，它的亮度都超过1.0，简单的把颜色截断至$[0,1]$范围内，那么整个场景有会显得曝光过度，物体的细节无法有效的得到呈现，如图6所示，整块区域就显得白花花的。可以模仿摄影中的色阶重建技术，同时渲染多张不同曝光度的高动态范围的图像，再通过色阶重建技术得到低动态范围的图像。那么问题来了，在实时渲染中，颜色缓冲区会把亮度超过1.0的图像自动截断至1.0，换句话说，就是GL_RGB类型的颜色缓冲区是无法存储高动态范围的帧，此时又如何存储渲染出来的高动态范围的帧呢？

图6. 曝光过度的渲染效果

这就有了新的图像存储格式的诞生，最早由Greg Ward在1985年提出了RGBE 32位的图像存储格式，RGB表示颜色值，E（Exponent）表示RBG颜色的指数系数，通过指定不同的E，就能有效的存储不同高动态范围的图像。后续又出现了多种不同的HDR格式的标准，例如Pixar Log Encoding（TIFF）、Radiance RGBE Encoding (HDR)、SGI LogLuv (TIFF)、ILM OpenEXR (EXR)^[9]、Microsoft/HP scRGB Encoding。支持HDR的图像存储格式，在许多地方都有所应用，基于物理的全局光照技术（Global Illumination Technique, i.e., Physically-Based Rendering），数字电影等。举个实时渲染的例子，在DX安装包DirectX Sample Browser中提供了一个示例，采用HDR格式存储的图像实现环境贴图，只需要额外对HDR格式图像的解析消耗，就可以实现高动态范围的渲染效果，如图7所示，图中的夜壶还经过了Bloom后处理技术（略有光晕的效果），会不会觉得效果很赞？！当然，HDR图像存储格式和Bloom后处理又是与高动态范围渲紧密相关的两个主题，本篇文章不作详细讨论。

图7. 采用HDR图像格式的图像实现的环境贴图

而在现代的硬件设备中，为了支持HDR，提供了浮点缓冲区（Floating Point Buffer），规定颜色缓冲区的内部格式（Internal Format）为：GL_RGB16F、GL_RGBA16F、GL_RGB32F、GL_RGBA32F，这样就支持亮度超过1.0的颜色存储，这种方案就能完美的解决高动范围的渲染存储问题。当然，这就要求比普通的RGB颜色缓冲区更大的存储消耗。在GL中，创建一个浮点数缓冲区，用于高动态范围的渲染的实现代码如下所示

	// Set up floating point framebuffer to render scene
	GLuint hdrFBO;
	glGenFramebuffers(1, &hdrFBO);

	// - Create floating point color buffer
	GLuint colorBuffer;
	glGenTextures(1, &colorBuffer);
	glBindTexture(GL_TEXTURE_2D, colorBuffer);
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA16F, SCR_WIDTH, SCR_HEIGHT, 0, GL_RGBA, GL_FLOAT, NULL);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

	// - Create depth buffer (renderbuffer)
	GLuint rboDepth;
	glGenRenderbuffers(1, &rboDepth);
	glBindRenderbuffer(GL_RENDERBUFFER, rboDepth);
	glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT, SCR_WIDTH, SCR_HEIGHT);

	// - Attach buffers
	glBindFramebuffer(GL_FRAMEBUFFER, hdrFBO);
	glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, colorBuffer, 0);
	glFramebufferRenderbuffer(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, rboDepth);
	if (glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE){
		printf("Framebuffer not complete!\n");
		exit(0);
	}

创建完浮点数缓冲区后，渲染高动态范围的帧：

	// Render scene into floating point framebuffer
	glBindFramebuffer(GL_FRAMEBUFFER, hdrFBO);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
	// Render the scene.
	renderScene();
	
	glBindFramebuffer(GL_FRAMEBUFFER, 0);

最后，通过色阶重建算子（Tone Mapping Operator），把浮点缓冲区中的帧数据映射至普通的RGB颜色缓冲区，输出到显示设备上，这就是最简单的色阶重建技术在GL中的实现流程：

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	// Tone Mapping.
	toneMappingShader.use();
	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, colorBuffer);
	renderQuad();

现在考虑色阶重建算子的实现，介绍一种Reinhard在2002年^[11]提出的全局色阶重建算子，该算法的核心代码如等式（1）所示：

$\overline{{{L}_{w}}}=\frac{1}{N}\exp \left( \sum\limits_{x,y}{\log \left( \delta +{{L}_{w}}\left( x,y \right) \right)} \right) \tag{1}$

其中，${{L}_{w}}\left( x,y \right)$表示图像像素$\left( x,y \right)$的亮度，$N$表示图像的像素个数，$\delta$是一个很小的数值，为了避免对数函数的真数为0。图像上每个像素的亮度就通过计算出来的亮度平均值$\overline{{{L}_{w}}}$进行缩放，缩放函数为：

$L\left( x,y \right)=\frac{a}{\overline{{{L}_{w}}}}{{L}_{w}}\left( x,y \right) \tag{2}$

其中，$a$是一个介于0与1之间的缩放系数。$a$值较大，则可以模拟白天的场景渲染；$a$值较小，则可以模拟傍晚的场景渲染。

另外，还有一个更简单粗暴的色阶重建算子，如等式（3）所示，这个等式可以保证高动态范围的亮度在低动态范围内有效表示出来，${{L}_{x}}\left( x,y \right)$越大，计算出来的${{L}_{d}}\left( x,y \right)$越接近1.0。

${{L}_{d}}\left( x,y \right)=\frac{{{L}_{w}}\left( x,y \right)}{1+{{L}_{w}}\left( x,y \right)} \tag{3}$

扩展等式（3），规定一个场景中能表示的最大的亮度值，设为${{L}_{white}}$，所有大于${{L}_{white}}$的亮度可以被截断至1.0，得到等式（4），它的函数函数曲线如图8所示。

${{L}_{d}}\left( x,y \right)=\frac{{{L}_{w}}\left( x,y \right)\left( 1+\frac{L\left( x,y \right)}{L_{white}^{2}} \right)}{1+L\left( x,y \right)} \tag{4}$

图8. 等式（4）的函数曲线

同时，在Reinhard色阶重建算子中添加伽马校正能达到更好的效果，相应的片断着色器代码如下所示，渲染出来的效果如图9所示。

#version 330 core
out vec4 color;
in vec2 TexCoords;

uniform sampler2D hdrBuffer;
uniform bool hdr;

void main(){             
    const float gamma = 2.2;
    vec3 hdrColor = texture(hdrBuffer, TexCoords).rgb;
    if(hdr){
        // reinhard
         vec3 result = hdrColor / (hdrColor + vec3(1.0));
        // also gamma correct while we're at it       
        result = pow(result, vec3(1.0 / gamma));
        color = vec4(result, 1.0f);
    }
    else{
        vec3 result = pow(hdrColor, vec3(1.0 / gamma));
        color = vec4(result, 1.0);
    }
}

图9. Reinhard色阶重建算子的渲染效果图

在真实的世界中，最亮的光照和最暗的光照间的亮度范围非常的大，它们的比值可以达到${{10}^{12}}:1$，人眼为了能看清眼前的场景，会根据进入眼中的光线强弱动态调整曝光度，简单说，人眼就是一个能够动态调整曝光度的照相机。例如，当人从黑暗的隧道走出^[6]，第一反应是感觉阳光很刺眼，但是慢慢地，会适应当前的亮度，逐渐看清眼前的场景，这个过程称为光适应（Light Adaptation），光适应的过程就是人眼动态调整曝光度的过程。计算机，也可以通过调整曝光度参数来控制场景的亮度，这样的色阶重建算子可以表示为：

${{L}_{d}}\left( x,y \right)=1-\exp \left( -L\left( x,y \right)\cdot \operatorname{exposure} \right) \tag{5}$

其中，$\operatorname{exposure}$表示曝光度，通过对曝光度的控制，来调整场景的整体亮度。

采用曝光度参数的片断着色器代码如下所示，渲染出来的效果如图10所示。

#version 330 core
out vec4 color;
in vec2 TexCoords;
uniform sampler2D hdrBuffer;
uniform bool hdr;
uniform float exposure;

void main(){             
    const float gamma = 2.2;
    vec3 hdrColor = texture(hdrBuffer, TexCoords).rgb;
    if(hdr){
        // exposure
        vec3 result = vec3(1.0) - exp(-hdrColor * exposure);
        // also gamma correct while we're at it       
        result = pow(result, vec3(1.0 / gamma));
        color = vec4(result, 1.0f);
    }
    else{
        vec3 result = pow(hdrColor, vec3(1.0 / gamma));
        color = vec4(result, 1.0);
    }
}

       图10. 不同曝光度的渲染效果

根据时间戳动态调整曝光度，可以模拟人眼由黑色的隧道走出看到场景的效果，起始渲染时，采用比较大的曝光度，然后逐渐降低曝光度至正常的水平，这也是一个很有意思的应用实例。上述介绍的色阶重建的完整实现代码参见[11]，但是它并没有考虑整个屏幕的光照效果，只是独立的对某些像素进行处理。

还有一种常规的算法，根据整个屏幕的光照强度，计算出一个光照的均值，利用该均值实现HDR效果，它的基本流程包括以下几个步骤：

1. 渲染初始纹理帧

采用浮点缓冲区，渲染高动光场景，得到初始帧纹理$fram{{e}_{original}}$。可以采用Phong光照渲染场景，得到帧纹理，如图11所示，顶点着色器和片断着色器的代码如下所示：

图11. 帧纹理$fram{{e}_{original}}$

#version 330 core
layout (location = 0) in vec3 position;
layout (location = 1) in vec3 normal;
layout (location = 2) in vec2 texCoords;
out VS_OUT{
	vec3 fragPos;
	vec3 normal;
	vec2 texCoords;
} vs_out;

uniform mat4 projMatrix;
uniform mat4 viewMatrix;
uniform mat4 modelMatrix;

void main(){
	gl_Position = projMatrix * viewMatrix * modelMatrix * vec4(position, 1.0f);
	vs_out.fragPos = vec3(modelMatrix * vec4(position, 1.0));   
	vs_out.texCoords = texCoords;
	vs_out.normal = normalize(mat3(modelMatrix) * normal);
}

#version 330 core
out vec4 fragColor;
in VS_OUT {
    vec3 fragPos;
    vec3 normal;
    vec2 texCoords;
} fs_in;
struct Light {
    vec3 position;
    vec3 color;
};

const int LIGHT_NUMS = 2;
uniform Light lights[2];
uniform float phongCoeff;
uniform float phongExp;
uniform float diffCoeff;
uniform sampler2D diffuseTexture;

void main(){           
    vec3 color = texture(diffuseTexture, fs_in.texCoords).rgb;
    vec3 n = normalize(fs_in.normal);
    // Ambient
    float amb = 0.05;
    // Lighting
    vec3 lighting = vec3(0.0f);
    for(int i = 0; i < 2; i++){
        // Diffuse
        vec3 l = normalize(lights[i].position - fs_in.fragPos);
        // because of the floating frame buffer, clamp() hear.
        float diff = max(dot(l, n), 0.0);
        float spec = max(dot(reflect(-l, n), n), 0.0);
        vec3 res = lights[i].color * (diff + amb + phongCoeff * pow(spec, phongExp)) * color;     
        // Attenuation (use quadratic as we have gamma correction)
        float d = length(fs_in.fragPos - lights[i].position);
        res *= 1.0 / (d * d);
        lighting += res;
    }
    //lighting = vec3(1.0f, 0.0f, 0.0f);
    fragColor = vec4(lighting, 1.0f);
}

2. 缩放初始帧纹理

将初始帧纹理缩放至1/4，得到缩放帧纹理$fram{{e}_{scale}}$，缩放之后的帧纹理如图12所示，片断着色器代码如下所示：

图12. 帧纹理$fram{{e}_{scale}}$

#version 330 core
out vec4 color;
in vec2 TexCoords;

uniform sampler2D hdrBuffer;
uniform bool hdr;
uniform float exposure;

void main(){             
    const float gamma = 2.2;
    vec3 hdrColor = texture(hdrBuffer, TexCoords).rgb;
    if(hdr){
        // exposure
        vec3 result = vec3(1.0) - exp(-hdrColor * exposure);
        // also gamma correct while we're at it       
        result = pow(result, vec3(1.0 / gamma));
        color = vec4(result, 1.0f);
    }
    else{
        vec3 result = pow(hdrColor, vec3(1.0 / gamma));
        color = vec4(result, 1.0);
    }
}

3. 提取高光纹理

从缩放帧纹理中提取高光部分，得到高光纹理$fram{{e}_{hdr}}$，高光纹理如图13所示，片断着色器代码如下所示：

图13. 帧纹理$fram{{e}_{hdr}}$

#version 330 core
out vec4 color;
in vec2 TexCoords;

// The per-color weighting to be used for luminance calculations in RGB order.
const vec3 BRIGHT_PASS_THRESHOLD  	= vec3(6.0, 6.0, 6.0);
const vec3 BRIGHT_PASS_OFFSET 		= vec3(10.0, 10.0, 10.0);
uniform sampler2D quadBuffer;

void main(){
	vec3 sample = texture2D(quadBuffer, TexCoords).rgb;
	sample -= BRIGHT_PASS_THRESHOLD;
	// Clamp to 0
	sample = max(sample, 0.0f);
	// Map the resulting value into the 0 to 1 range. Higher values for
	// BRIGHT_PASS_OFFSET will isolate lights from illuminated scene 
	// objects.
	color = vec4(sample / (BRIGHT_PASS_OFFSET + sample), 1.0);
}

4. 计算平均亮度值

可以采用公式（6）将颜色值转化为相应的亮度，再根据公式（1）计算出平均亮度，得到初始的亮度纹理，如图14所示，（注意：此时不做取对数的运算），片断着色器代码如下所示：

$Lum\left( R,G,B \right)=0.2125R+0.7154G+0.0721B \tag{6}$

图14. 帧纹理$fram{{e}_{lum}}$

#version 330 core
out vec4 color;
in vec2 TexCoords;

// The per-color weighting to be used for luminance calculations in RGB order.
const vec3 LUMINANCE_VECTOR  = vec3(0.2125f, 0.7154f, 0.0721f);
uniform sampler2D quadBuffer;
uniform int width;
uniform int height;

void main(){
	vec3 sample;
	float logLumSum = 0.0;
	vec2 offset;
	int x, y;
	float tu = 1.0 / float(3 * width), tv = 1.0 / float(3 * height);
	for(x = -1; x <= 1 ; x++){
		for(y = -1; y <= 1; y++){
			offset = vec2(float(x) * tu, float(y) * tv);
			sample = texture2D(quadBuffer, TexCoords + offset).rgb;
			logLumSum += log(dot(sample, LUMINANCE_VECTOR) + 0.0001f);
		}
	}
	logLumSum /= 9;
	color = vec4(logLumSum, logLumSum, logLumSum, 1.0);
}

接着，把亮度纹理帧$fram{{e}_{lum}}$逐步缩放至$2\times 2$大小，如图15所示（为了方便观察，把图像放大了三倍），片断着色器代码如下所示：

图15. 对亮度纹理帧的缩放

#version 330 core
out vec4 color;
in vec2 TexCoords;

uniform sampler2D quadBuffer;
uniform int width;
uniform int height;

void main(){
	vec3 sample;
	float logLumSum = 0.0;
	vec2 offset;
	int x, y;
	float tu = 1.0 / float(width), tv = 1.0 / float(height);
	for(x = 0; x < 4 ; x++){
		for(y = 0; y < 4; y++){
			offset = vec2(float(x - 1.5) * tu, float(y - 1.5) * tv);
			logLumSum += texture2D(quadBuffer, TexCoords + offset).r;
		}
	}
	logLumSum /= 16;
	color = vec4(logLumSum, logLumSum, logLumSum, 1.0);
}

最后，把$2\times 2$的纹理再缩放至$1\times 1$的纹理，进行对数运算，计算出平均亮度，此时，$1\times 1$的纹理$fram{{e}_{ave}}$就存储着整个场景的平均亮度，片断着色器代码如下所示：

#version 330 core
out vec4 color;
in vec2 TexCoords;

uniform sampler2D quadBuffer;
uniform int width;
uniform int height;

void main(){
	float logLumSum = 0.0;
	vec2 offset;
	int x, y;
	float tu = 1.0 / float(width), tv = 1.0 / float(height);
	for(x = 0; x < 4 ; x++){
		for(y = 0; y < 4; y++){
			offset = vec2(float(x - 1.5) * tu, float(y - 1.5) * tv);
			logLumSum += texture2D(quadBuffer, TexCoords + offset).r;
		}
	}
	logLumSum = exp(logLumSum / 16);
	color = vec4(logLumSum, logLumSum, logLumSum, 1.0);
}

5. 利用平均亮度，重新计算初始帧纹理$fram{{e}_{original}}$上每个像素的颜色

利用公式（2）（3）计算初始帧纹理上每个像素的RGB值，得到最终的渲染效果，如图16所示，图（1）与初始帧存在差异，是因为它采用了伽马校正。片断着色器代码如下所示：

图16. 色阶重建前后的对比效果

#version 330 core
out vec4 color;
in vec2 TexCoords;

uniform sampler2D sceneBuffer;
uniform sampler2D bloomBuffer;
uniform sampler2D hdrBuffer;

void main(){
	vec3 sample = texture2D(sceneBuffer, TexCoords).rgb;
	vec3 bloom = texture2D(bloomBuffer, TexCoords).rgb;
	float adapt = texture2D(hdrBuffer, vec2(0.5, 0.5)).r;
	sample /= adapt;
	sample /= (vec3(1.0) + sample);
	sample *= 0.8;
	color = vec4(sample, 1.0);
}

接下来，阐述Bloom光晕效果的实现原理和实现，首先介绍Bloom技术需要的高斯滤波的原理。

2. 高斯滤波

高斯滤波器（Gaussian Smoothing）是一个二维的卷积算子（Convolution Operator），用于模糊化图像，移除图像中的细节和噪音，这与均值滤波器（Mean Filter）有类似的效果，只是它采用了不同的核函数。这是图像处理的范畴，就不得不提冈萨雷斯^[12]的《数字图像处理》这本巨作，能力有限只能摘取其中的部分内容进行阐述。

17. 图像（左）和卷积核（右）

所谓的卷积核^[13]是什么东西，每个图像都可以看成是一个像素构成的矩阵，设该图像有$N$行$M$列，如图17（左）所示，卷积核可以看成比该图像小的数值矩阵，设有$n$行$m$列，如图17（右）所示。对该图像进行卷积运算，输出的图像就有$N-n+1$行$M-m+1$列，计算公式表示为：

$F\left( i,j \right)=\sum\limits_{k=1}^{m}{\sum\limits_{l=1}^{n}{{{I}_{i+k-1,j+l-1}}\cdot {{K}_{k,l}}}} \tag{7}$

举个例子，对图17（左）的图像进行处理，新图像的第1行第1列的像素的计算方法为：

${{F}_{11}}={{I}_{11}}\cdot {{K}_{11}}+{{I}_{12}}\cdot {{K}_{12}}+{{I}_{13}}\cdot {{K}_{13}}+{{I}_{21}}\cdot {{K}_{21}}+{{I}_{22}}\cdot {{K}_{22}}+{{I}_{23}}\cdot {{K}_{23}}+{{I}_{31}}\cdot {{K}_{31}}+{{I}_{32}}\cdot {{K}_{32}}+{{I}_{33}}\cdot {{K}_{33}}$

而均值滤波和高斯滤波的区别就在于采用的卷积核不同，高斯滤波采用的卷积核是一个二维的高斯函数：

$G\left( x,y \right)=\frac{1}{2\pi {{\sigma }^{2}}}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2{{\sigma }^{2}}}}} \tag{8}$

对于标准方差$\sigma =1.0$的高斯滤波函数的函数分布图，如图18所示。对图像的滤波处理，其实是一个二维核函数的积分过程，如等式（9）所示：

$L\left( x,y \right)=\int\limits_{-\infty }^{+\infty }{\int\limits_{-\infty }^{+\infty }{I\left( X-x,Y-y \right)\cdot G\left( x,y \right)}}dxdy \tag{9}$

其中，$I\left( x,y \right)$表示图像信息，$G\left( x,y \right)$表示核函数，$L\left( x,y \right)$表示处理后的图像信息。

图18. 二维高斯函数的函数分布图

把等式（8）代入等式（9）中，得到：

$L\left( x,y \right)=\int\limits_{-\infty }^{+\infty }{\int\limits_{-\infty }^{+\infty }{I\left( X-x,Y-y \right)\cdot \frac{1}{2\pi {{\sigma }^{2}}}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2{{\sigma }^{2}}}}}}}dxdy \tag{10}$

该积分方程可以表示为有限数的和，如等式（11）所示，这就是二维的高斯滤波函数的方程式，每个滤波生成一个像素，都需要$O\left( {{n}^{2}} \right)$的运算规模。

$L\left( x,y \right)=\sum\limits_{y=-n/2}^{n/2}{\sum\limits_{x=-m/2}^{m/2}{I\left( X-x,Y-y \right)\cdot \frac{1}{2\pi {{\sigma }^{2}}}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2{{\sigma }^{2}}}}}}} \tag{11}$

其中，$\sum\limits_{y=-n/2}^{n/2}{\sum\limits_{x=-m/2}^{m/2}{\frac{1}{2\pi {{\sigma }^{2}}}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2{{\sigma }^{2}}}}}}}=1$。

由于高斯卷积的可分离性，可以把运算规模降低到$O\left( n \right)$。首先介绍下可分离卷积^[15]（Separable Convolution）是什么，可以这样定义^[16]：

如果一个二维滤波核能用两个向量的积来表示，那么它就是可分离性。

举个例子，设矩阵$M$为：

$$M = \left[ {\begin{array}{*{20}{c}}{ - 1}&0&1\\{ - 2}&0&2\\{ - 1}&0&1\end{array}} \right]$$

它可以通过一个列向量和一个行向量的积表示：

$$M = V \cdot H = \left[ {\begin{array}{*{20}{c}}1\\2\\1\end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}}{ - 1}&0&1\end{array}} \right]$$

我们就可以称$M$为可分离性卷积，它有一个重要的性质，如等式（12）所示，它的证明参见Songho^[17]：

$S\cdot M=S\cdot \left( V\cdot H \right)=\left( S\cdot V \right)\cdot H \tag{12}$

高斯滤波核由于它的可分离性，二维的高斯滤波方程式可以表示为：

$L\left( x,y \right)=\left[ \sum\limits_{y=-n/2}^{n/2}{I\left( X,Y-y \right)\cdot \frac{1}{\sqrt{2\pi }\sigma }{{e}^{-\frac{{{y}^{2}}}{2{{\sigma }^{2}}}}}} \right]\cdot \left[ \sum\limits_{x=-m/2}^{m/2}{I\left( X-x,Y \right)\cdot \frac{1}{\sqrt{2\pi }\sigma }{{e}^{-\frac{{{x}^{2}}}{2{{\sigma }^{2}}}}}} \right] \tag{13}$

以$5\times 5$的高斯滤波核为例，如下所示：

图19. 二维高斯滤波系数

通过将二维的高斯模糊处理可分解为对图像分别进行垂直和水平两通道的模糊处理，每个像素需要25次的乘法运算，减低为10次，把运算规模由$O\left( {{n}^{2}} \right)$降低为$O\left( n \right)$，更重要的是它降低了对GPU的访问次数，计算机中GPU访问纹素的时间损耗比计算损耗还大。

最后，提供一个GLSL实现的高斯模糊的代码示例，在片断着色器里面采用实时的计算高斯系数，在实际的应用中，可以把数值预先计算好存储起来，不需要每次计算。采用常规的高斯滤波方法和行/列双通道的高斯滤波能达到的效果是一样，效果如图20所示。

图20. 高斯滤波效果图

#version 330 core

out vec4 color;
in vec2 TexCoords;

uniform sampler2D sampleBuffer;
uniform int width;
uniform int height;
uniform int sampleWidth;
uniform int sampleHeight;
uniform int sampleType;
const float PI = 3.141592653589;

float gaussianDistribution( float x, float y, float rho){
	float g = 1.0f / sqrt( 2.0f * PI * rho * rho );
	g *= exp( -( x * x + y * y ) / ( 2 * rho * rho));
	return g;
}

vec4 gaussBlurOrg(int smW, int smH, int bW, int bH){
	int i, j;
	float sum = 0.0;
	float weight;
	vec2 offset;
	int w = smW / 2, h = smH / 2;
	vec4 sample = vec4(0.0);
	float tu = 1 / float(bW), tv = 1 / float(bH);
	for(i = - w; i <= w ; i ++)
	{
		for(j = -h; j <= h; j ++)
		{
			offset = vec2(tu * i, tv * j);
			weight = gaussianDistribution(i, j, 1.0);
			sample += weight * texture2D(sampleBuffer, TexCoords + offset);
			sum += weight;
		}
	}
	sample /= sum;
	return sample;
}

vec4 gaussBlurVert(int smH, int bH){
	int i;
	float sum = 0.0;
	float weight;
	vec2 offset;
	int h = smH / 2;
	vec4 sample = vec4(0.0);
	float tv = 1 / float(bH);
	for(i = - h; i <= h ; i ++){
		offset = vec2(0, tv * i);
		weight = gaussianDistribution(0, i, 1.0);
		sample += weight * texture2D(sampleBuffer, TexCoords + offset);
		sum += weight;
	}
	sample /= sum;

	return sample;
}

vec4 gaussBlurHori(int smW, int bW){
	int i;
	float sum = 0.0;
	float weight;
	vec2 offset;
	int h = smW / 2;
	vec4 sample = vec4(0.0);
	float tu = 1 / float(bW);
	for(i = - h; i <= h ; i ++){
		offset = vec2(tu * i, 0);
		weight = gaussianDistribution(i, 0, 1.0);
		sample += weight * texture2D(sampleBuffer, TexCoords + offset);
		sum += weight;
	}
	sample /= sum;
	return sample;
}

void main(){
	if(sampleType == 0)
		color = gaussBlurOrg(sampleWidth, sampleHeight, width, height);
	else if(sampleType == 1)
		color = gaussBlurHori(sampleWidth, width);
	else if(sampleType == 2)
		color = gaussBlurVert(sampleHeight, height);
}

3. Bloom

本节主要介绍基于高动态光照实现的光晕效果，与前面介绍的色阶重建技术相比，它相当于把第3阶段得到的高光纹理做模糊处理，然后组合到最终的纹理帧上，如图21所示。

图21. Bloom技术

对高光纹理进行高斯模糊处理，可以采用行/列双通道的高斯模糊处理，片断着色器代码如下所示：

#version 330 core

out vec4 color;
in vec2 TexCoords;

uniform sampler2D sampleBuffer;
uniform int width;
uniform int height;
uniform int sampleType;

const float GAUSSIAN[5] = float[5](0.054489, 0.244201, 0.402620, 0.244201, 0.054489);

vec4 gaussBlurHori(){
	vec4 sample = vec4(0.0);
	float tu = 1 / float(width);
	sample += GAUSSIAN[0] * texture2D(sampleBuffer, TexCoords + vec2(-2 * tu, 0));
	sample += GAUSSIAN[1] * texture2D(sampleBuffer, TexCoords + vec2(-tu, 0));
	sample += GAUSSIAN[2] * texture2D(sampleBuffer, TexCoords + vec2(0, 0));
	sample += GAUSSIAN[3] * texture2D(sampleBuffer, TexCoords + vec2(tu, 0));
	sample += GAUSSIAN[4] * texture2D(sampleBuffer, TexCoords + vec2( 2 * tu, 0));

	return sample;
}

vec4 gaussBlurVert(){
	vec4 sample = vec4(0.0);
	float tv = 1 / float(height);
	sample += GAUSSIAN[0] * texture2D(sampleBuffer, TexCoords + vec2(0, -2 * tv));
	sample += GAUSSIAN[1] * texture2D(sampleBuffer, TexCoords + vec2(0, -tv));
	sample += GAUSSIAN[2] * texture2D(sampleBuffer, TexCoords + vec2(0, 0));
	sample += GAUSSIAN[3] * texture2D(sampleBuffer, TexCoords + vec2(0, tv));
	sample += GAUSSIAN[4] * texture2D(sampleBuffer, TexCoords + vec2(0, 2 * tv));

	return sample;
}

void main(){
	if(sampleType == 1)
		color = gaussBlurHori();
	else if(sampleType == 2)
		color = gaussBlurVert();
}

在得到模糊化的高光纹理后，将它加到最终的渲染帧上，片断着色器代码如下所示：

#version 330 core

out vec4 color;
in vec2 TexCoords;

uniform sampler2D sceneBuffer;
uniform sampler2D bloomBuffer;
uniform sampler2D hdrBuffer;
uniform int useHdr;
uniform int useBloom;

void main(){
	vec3 sample = texture2D(sceneBuffer, TexCoords).rgb;
	vec3 bloom = texture2D(bloomBuffer, TexCoords).rgb;
	if(useHdr == 1){
		float adapt = texture2D(hdrBuffer, vec2(0.5, 0.5)).r;
		sample /= adapt;
		sample /= (vec3(1.0) + sample);
		
	}
	sample *= 0.8;
	if(useBloom == 1)
	{
		sample += bloom;
	}
	color = vec4(sample, 1.0);
}

 

参考

[1] Matt Grum. "What is tone mapping? How does it relate to HDR?." website < http://photo.stackexchange.com/questions/7630/what-is-tone-mapping-how-does-it-relate-to-hdr>.
[2] Wikipedia. "High-dynamic-range imaging." website<https://en.wikipedia.org/wiki/High-dynamic-range_imaging>.
[3] Wikipedia. "高动态范围成像."website<https://zh.wikipedia.org/wiki/%E9%AB%98%E5%8A%A8%E6%80%81%E8%8C%83%E5%9B%B4%E6%88%90%E5%83%8F>.
[4] Erik Reinhard, et al. High dynamic range imaging: acquisition, display, and image-based lighting. Morgan Kaufmann, 2010.
[5] Joey De Vries. "High-dynamic-range imaging." website<http://learnopengl.com/#!Advanced-Lighting/HDR >.
[6] Tomas Akenine-Möller, Eric Haines, and Naty Hoffman. Real-time rendering. CRC Press, 2008.
[7] NBRANCACCIO. "HDR in WebGL." website<http://www.floored.com/blog/202015webgl-hdr/>.
[8] Greg Ward. "High Dynamic Range Image Encodings." website<http://www.anyhere.com/gward/hdrenc/hdr_encodings.html>, anyhere software.
[9] Florian Kainz, Rod Bogart, and Drew Hess. "The OpenEXR image file format." GPU Gems, 2004. website< http://http.developer.nvidia.com/GPUGems/gpugems_ch26.html>.
[10] Erik Reinhard, et al. "Photographic tone reproduction for digital images." ACM Transactions on Graphics (TOG), vol.21, no. 3, ACM, 2002.
[11] Joey De Vries. "LearnOpenGL." website<https://github.com/JoeyDeVries/LearnOpenGL>.
[12] Rafael C. Gonzalez, and Richard E. Woods. Digital image processing. Nueva Jersey, 2008.
[13] R. Fisher, S. Perkins, A. Walker and E. Wolfart. "Convolution." website<http://homepages.inf.ed.ac.uk/rbf/HIPR2/convolve.htm>, 2003.
[14] R. Fisher, S. Perkins, A. Walker and E. Wolfart. "Gaussian Smoothing." website<http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm>, 2003.
[15] songho. "Convolution." website<http://www.songho.ca/dsp/convolution/convolution.html>.
[16] Steve Eddins . "Separable convolution." website <http://blogs.mathworks.com/steve/2006/10/04/separable-convolution/>.
[17] Songho. "Proof of Separable Convolution 2D. " website<http://www.songho.ca/dsp/convolution/convolution2d_separable.html>.