2008 年在一个 PS 讨论群里,有网友不解 Photoshop 的高斯模糊中的半径是什么含义,因此当时我写了这篇文章:
在那篇文章中,主要讲解了高斯模糊中的半径的含义,是二维正态分布的方差的平方根,并且给出了算法的理论描述。现在我又打算把该算法用 c++ 实现出来,于是有了下面的这个 DEMO。
起初我是按照算法理论直接实现,即使用了二维高斯模板,结果发现处理时间很长,对一个图片竟然能达到大约数分钟之久。这样肯定是不对的,所以我百度了一下,发现这个问题应该采用分别进行两次一维高斯模糊就可以了[1],这样算法的时间复杂度的一个系数,就从 O ( σ ^2 ) 降低到了 O ( σ )。这样算法于是速度提高到了毫秒级。下表给出分别用二维模糊的原始方法,和两次一维模糊累加的方法的算法成本比较:
算法 | 时间复杂度 | 空间复杂度 |
---|---|---|
(1) 二维高斯模糊 | O(σ ^ 2) * O(n) (慢) | O(σ ^ 2) (较小) |
(2) 两次一维高斯模糊的累加 | O(σ) * O(n) (快) | O(n) + O(σ) ≈ O(n) (较大) |
其中:σ :方差平方根(Photoshop 中的高斯模糊半径);n = w * h (图片的像素数量)。具体时间和图片大小和高斯半径的大小有关,一个粗略的大概情况为,算法(1)的耗时为分钟级,算法(2)的耗时为毫秒到秒级。可见算法(2)比算法(1)速度更快,但相比算法(1)来说算法(2)具有较高的空间需求。
注:当然上面的空间复杂度并不是绝对的,例如,可以通过对图像进行串行的切片处理,既可减小算法(2)的空间需求。
两种算法在高斯半径为常数条件下,都是关于图片大小的线性算法,区别在于常数系数的大小不同,前者是高斯半径(模板尺寸)的平方级,后者是高斯半径(模板尺寸)的线性级别。这个改进,非常类似于我此前有一篇博客中给出的,对一个油画效果滤镜的算法改进,也是通过把常数系数,从模板尺寸的平方级别降低到线性级别,使算法速度获得提高的。
在理论上,高斯模板是无边界和无限扩展的一个二维曲面,在实现时,就必须对这个曲面截断为有限大的二维模板。那么在哪里截断呢?根据下图所示的一维正态分布贡献:
图1. 正态分布的贡献比
此图来自参考资料 [1],根据资料文中叙述,此图实际来源于(Maybe blocked by the GFW):
http://zh.wikipedia.org/wiki/File:Standard_deviation_diagram.svg。
从图 1 中可以看到,在 3σ 以外的贡献比例非常小,为 0.1 %,因此我们截断模板时,对模板边界定义为 3σ ;
int r = ( int ) ( sigma * 3 + 0.5 ); // 完整模板的逻辑尺寸:( 2 * r + 1 ) * ( 2 * r + 1 );
二维高斯模板的计算公式是:
下图给出了二维模板的可视化结果。采用的可视化方法是,根据上面的公式和模板边界,生成二维高斯模板,然后取一个缩放因子 f = 255 / 模板中心点的数据,以此缩放因子把模板数据等比缩放,然后绘制成灰度图片,这样中心点的亮度就被提高到最亮。可视化效果中,每个单元格对应着一个模板数据,单元格大小为 8 * 8 或者 16 * 16 像素。
图 2. 二维模板的可视化结果
图 2 中,左侧是人们常见的 3 x 3 模板(σ ≈ 0.849),围绕中心点的 3 x 3 的浮点数据为:
sigma = 0.849: 0.055 0.110 0.055 0.110 0.221 0.110 0.055 0.110 0.055
采用算法(2),要完成高斯模糊,对图片分别进行两个方向的一维高斯模糊即可。例如,先对图片进行水平方向的模糊,得到中间结果,然后再对这个中间结果进行垂直方向的模糊,即得到最终结果。下图是一个演示图,给出了原图在两个方向上分别单独进行一维高斯模糊的结果,以及最终的结果:
图3. 算法(2)中一维模糊的中间结果
仅在这个图片的例子中,我把我写的算法的处理结果,在 Photoshop 中打开和 Photoshop 自带的高斯模糊的处理结果做差值对比,发现两者是相同的。
我实现的 DEMO 程序(Windows 平台)的界面如下所示:
图4. DEMO 程序的主窗口 UI
通过点击菜单 - 可视化 - 二维高斯模板可以在右侧的视图中生成一个灰度图片,即二维高斯模板的可视化结果。
在程序界面的客户区下方有一个控制面板,可以选择高斯模糊的算法参数,高斯半径的意义和 Photoshop 中的高斯模糊半径的意义相同,都是算法中的 σ。
算法参数中:
(1)支持多线程处理。根据我的观察,线程数设置为和 CPU 核心数相同是比较合适的。线程数比 CPU 核心数更多,也是没有什么意义的,因为算法执行时,CPU 已经满负荷运转了。开启更多线程,也不能再提高速度了。
假设 CPU 核数为 p,开启的多线程数量 >= p,则算法速度大约为单线程处理的 p 倍。(当 CPU 满负荷时,线程数量取得更大,也没有提高速度的意义了)
注:此处的 CPU 核心数应该为 CPU 的物理核心数,而非模拟出来的多核数目。
(2)浮点类型:支持 float 和 double。它是高斯模板的数据的类型,也是进行像素加权累加时的数据类型,根据我的观察,float 和 double 的速度相差不大。基本相同。
(3)高斯半径:即 σ。算法的常数系数为 O(σ)。很显然,σ 的值取得越大,算法耗时将会越长。在 DEMO 中,其允许范围和 Photoshop 的要求一致,是 [0.1, 250]。
在实现算法时,我也尝试了对 255 个灰度值 * 模板数据的结果进行缓存和查表处理,但是发现不能有效提高速度,所以最终我放弃了这种方法。这可能是因为,算法的计算只是一个浮点乘法,对数据的读取动作,并不能做到比浮点乘法更快。所以这里采用缓存也就显得没有必要了。
在本 DEMO 中,滤镜处理是放在 UI 线程中进行的,这使得在滤镜处理时间较长时(例如高斯半径取值很大,图片也很大),界面会有些卡,可以把滤镜处理动作放在一个新建的后台线程中执行。这是比较容易实现的。
对算法的使用方法:
在 C++ 程序中,使用我写的这个算法是非常简单的,例如:
#include "GaussBlurFilter.hpp" CGaussBlurFilter<double> _filter; _filter.SetSigma(3.5); // 设置高斯半径 _filter.SetMultiThreads(true, 4); // 开启多线程,用户建议的线程数为 4; // lpSrcBits / lpDestBits: 像素数据的起始地址,必须以 4 bytes 对齐,
// 注意:不论高度为正或者负,lpBits 都必须为所有像素中地址值最低的那个像素的地址。 // bmWidth, bmHeight: 图像宽度和高度(像素),高度允许为负值; // bpp: 位深度,支持 8(灰度), 24(真彩色), 32 _filter.Filter(lpSrcBits, lpDestBits, bmWidth, bmHeight, bpp);
需要注意的是,在多线程处理中,我使用了 Windows API (例如 CreateThread)等,这使得 GaussBlurFilter.hpp 目前只能用在 Windows 平台,如果要在其他平台使用,应当修改和多线程有关的 API 函数调用。
高度值可以为正也可以为负,但像素数据的地址 lpBits 都必须是所有像素中,地址值最小的那个像素的地址。即,假设图片左上角点的坐标为原点,如果图片高度为正数(bottom - up),则 lpBits 是左下角像素 (col = 0,row = height - 1)的地址。如果图片高度为负数(top-down),则 lpBits 是左上角像素(col = 0,row = 0) 的地址。图像数据的扫描行宽度必须以 4 Bytes 对齐,即通过下面的公式计算扫描行宽度:
int stride = ( bmWidth * bpp + 31 ) / 32 * 4; //扫描行宽度,对齐到 4 Bytes
(上式为编程语言表达,非数学表达,即利用了整数除法对小数部分的截断性。)
bpp:图像的像素位深度。只支持 8 (灰度索引图像),24,32 这几个值。对于 32 bpp 的图像来说,最后一个像素通道是表征像素的不透明度,也就是 alpha,对于 alpha 如何参与到算法中,我想了下,有多种处理方法,但都好像没有什么容易理解的物理意义,所以在代码里我忽略了 alpha 通道。
【相关下载】:
(1)Demo 可执行文件(包含 GaussBlurFilter.hpp):GaussBlurDemo_Bin.zip
(2)Demo 完整源码(包含 GaussBlurFilter.hpp 和 可执行文件):GaussBlurDemo_Src.zip
【参考资料】
[1]. 高斯模糊算法的实现和优化;
[注] 文中的公式,采用如下网址生成:http://www.codecogs.com/latex/eqneditor.php。
参考自:博客中插入公式——之在线数学公式生成;