二维高斯模糊和可分离核形式的快速实现

来源：互联网发布：js解析xml文件的作用编辑：程序博客网时间：2024/06/03 21:22

高斯模糊原理

基本概念

二维高斯模糊，或者说高斯滤波，是图像处理中非常常见的操作。操作的核心是使用一个从高斯分布中采样得到的掩膜，或者叫核，和输入图片中的每个像素及其邻域进行计算，结果保存到输出图片中。假设高斯核窗口尺寸为 (2w+1)×(2w+1)，高斯分布的标准差为 σ，则高斯核可以表示为矩阵的形式

G = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ G (- w, - w) ⋮ G (0, - w) ⋮ G (w, - w) \dots \dots \dots G (- w, 0) ⋮ G (0, 0) ⋮ G (w, 0) \dots \dots \dots G (- w, w) ⋮ G (0, w) ⋮ G (w, w) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

其中

G (u, v) = 1 S exp (- u 2 2 σ 2 - v 2 2 σ 2)

其中

u 表示行，

v 表示列，

u,v∈{−w,−w+1,…,w−1,w}，

S 是归一化常数

S = \sum u = - w w \sum v = - w w exp (- u 2 2 σ 2 - v 2 2 σ 2)

由于高斯分布的概率密度函数的非零值区间主要集中在

(−3σ,3σ) 内，所以为了保证选取的高斯核的完整性，一般取

w≈3σ。

说完了高斯核，该说高斯模糊的表达式了。设输入图片为 X，输出图片为 Y，第 i 行第 j 列的数据表示为 X(i,j) 和 Y(i,j)，则使用窗口大小为 (2w+1)×(2w+1)，标准差为 σ 的高斯核计算后的结果为

Y (i, j) = \sum u = - w w \sum v = - w w X (i + u, j + v) G (u, v) (1)

根据这个表达式，为了得到位置

(i,j) 上的输出，需要将高斯核的中心置于输入图片的位置

(i,j) 处，让高斯核的每一个值和输入图片对应位置上的值相乘，进行

(2w+1)×(2w+1) 次乘法计算，然后再进行

(2w+1)×(2w+1)−1 次加法计算，所以时间复杂度是

O(w2) 的。

可分离核形式实现

但是，注意到，高斯核的表达式是可分离的。下面为了表示方便，令

g (x) = exp (- x 2 2 σ 2)

则有

G (u, v) = 1 S g (u) g (v)

那么高斯核矩阵又可以改写成归一化常数乘以一个行向量乘以一个列向量的形式，如下

G = 1 S ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ g (- w) g (- w) ⋮ g (0) g (- w) ⋮ g (w) g (- w) \dots \dots \dots g (- w) g (0) ⋮ g (0) g (0) ⋮ g (w) g (0) \dots \dots \dots g (- w) g (w) ⋮ g (0) g (w) ⋮ g (w) g (w) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = 1 S ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ g (- w) ⋮ g (0) ⋮ g (w) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ \times [g (- w) \dots g (0) \dots g (w)]

另外

S = \sum u = - w w \sum v = - w w g (u) g (v) = \sum u = - w w [\sum v = - w w g (v)] g (u) = [\sum u = - w w g (u)] \times [\sum v = - w w g (v)] = S' \times S'

其中

S′ 可以认为是一维高斯核的归一化系数的倒数。所以有

G = G 1 \times G 2

其中

G 1 G 2 = 1 S ' [g (- w) \dots g (0) \dots g (w)] T = 1 S ' [g (- w) \dots g (0) \dots g (w)]

由此可见，

G 可以分离成两个向量的乘积的形式。

下面对 (1) 式进行改写

Y (i, j) = \sum u = - w w \sum v = - w w X (i + u, j + v) G (u, v) = \sum u = - w w \sum v = - w w X (i + u, j + v) 1 S g (u) g (v) = \sum u = - w w \sum v = - w w X (i + u, j + v) 1 S ' g (u) 1 S ' g (v) = \sum u = - w w [\sum v = - w w X (i + u, j + v) 1 S ' g (v)] 1 S ' g (u) = \sum u = - w w Z (i + u) 1 S ' g (u)

上面的式子表明，为了获得最终的高斯滤波的结果，可以先用横向一维高斯核

G2 与输入图片

X 进行计算，得到中间结果

Z。再用纵向一维高斯核

G1 与中间结果

Z 进行计算，得到输出

Y。利用

X (i + u, j - w), \dots, X (i + u, j), \dots, X (i + u, j + w)

计算得到

Z(i+u,j) 需要进行

(2w+1) 次乘法计算和

2w 次加法计算。利用

Z (i - w, j), \dots, Z (i, j), \dots, Z (i + w, j)

计算得到

Y(i,j) 需要进行

(2w+1) 次乘法计算和

2w 次加法计算。总的来说，计算出

Y(i,j) 的值需要进行

(4w+2) 次乘法计算和

4w 次加法计算，时间复杂度仅为

O(w)，比直接采用

(1) 式的计算方法快了很多。但是该算法需要使用和输入图片尺寸一致的内存保存中间结果。

实例分析

基本实现

假设我们有如下 8×8 的数据

X = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0816243240485619172533414957210182634425058311192735435159412202836445260513212937455361614223038465462715233139475563 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

我们要对这些数据进行窗口大小为

5×5，标准差为

σ=1.0 的高斯模糊。这时我们先计算纵向高斯核

G1 和横向高斯核

G2 的数值，得到

G 1 G 2 = [0.054 0.244 0.403 0.244 0.054] T = [0.054 0.244 0.403 0.244 0.054]

下面我们用

G2 对

X 的每一行进行模糊计算。以计算第 0 行为例。由于高斯核的长度等于 5，所以计算前要对第 0 行的数据进行扩展，这里采用镜像对称的方式进行扩展，扩展之后的数据为

[2 - 1 - 01234567 6 - 5 -]

其中加了下划线的数据就是扩展的数据。下图展示了高斯核从左向右移动并计算的过程，向下的箭头下面表示的是计算结果，也就是

Z 的第 0 行

2 0.054 1 0.244 0.054 0 0.403 0.244 0.054 0.706 1 0.244 0.403 0.244 1.109 2 0.054 0.244 0.403 2.000 3 0.054 0.244 3.000 4 0.054 \dots 0.054 ↓ 4.000 5 0.244 0.054 5.000 6 0.403 0.244 5.891 7 0.244 0.403 6.294 6 0.054 0.244 5 0.054

对

X 的每一行都这么处理，就可以得到中间结果

Z = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 0.706 8.706 16.706 24.706 32.706 40.706 48.706 56.706 1.109 9.109 17.109 25.109 33.109 41.109 49.109 57.109 2.000 10.000 18.000 26.000 34.000 42.000 50.000 58.000 3.000 11.000 19.000 27.000 35.000 43.000 51.000 59.000 4.000 12.000 20.000 28.000 36.000 44.000 52.000 60.000 5.000 13.000 21.000 29.000 37.000 45.000 53.000 61.000 5.891 13.891 21.891 29.891 37.891 45.891 53.891 61.891 6.294 14.294 22.294 30.294 38.294 46.294 54.294 62.294 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

下面我们用

G1 对

Z 的每一列进行模糊计算。以计算第 0 列为例。第 0 列采用镜像对称的方式进行扩展，扩展之后的数据为

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 16.706 - - - - - 8.706 - - - - - 0.706 8.706 16.706 24.706 32.706 40.706 48.706 56.706 48.706 - - - - - 40.706 - - - - - ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

其中加了下划线的数据就是扩展的数据。下图展示了高斯核从上向下移动并计算的过程，第一个向右的箭头右边表示的是浮点数计算结果，第二个向右箭头的右边表示四舍五入后成为整数的计算结果

16.706 8.706 0.706 8.706 16.706 24.706 32.706 40.706 48.706 56.706 48.706 40.706 0.054 0.244 0.403 0.244 0.054 0.054 0.244 0.403 0.244 0.054 0.054 0.244 0.403 0.244 0.054 \dots 0.054 0.244 0.403 0.244 0.054 0.054 0.244 0.403 0.244 0.054 \to 6.357 9.578 16.706 24.706 32.706 40.706 47.835 51.056 \to 610172533414851

对

Z 的每一列都这么处理，就可以得到模糊结果

Y = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 6101725334148517101725334148518111826344249529121927354350531013202836445154111421293745525512152230384653561215223038465357 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

符合局部性原则的内存访问加速

下面来考虑上述方法在内存访问效率方面的问题。利用 G2 和 X 计算 Z 的过程中，内存访问都是连续的，都是从左到右的形式。但是在利用 G1 和 Z 计算 Y 的过程中，取出每一列中的相邻数据，需要跨行。如果需要处理的图片宽度比较大，跨行访问数据可能会导致 Cache Miss，这是违反了内存访问局部性原则的。为了解决这一问题，利用 G1 和 Z 计算 Y 的方法需要调整。

实际上，利用 G1 和 Z 计算 Y 同样可以按行的方式计算。为了表述方便，以计算 Y 的第 2 行（下标从 0 开始）Y(2,⋅)为例，

Y (2, \cdot) = G 1 (0) Z (0, \cdot) + G 1 (1) Z (1, \cdot) + G 1 (2) Z (2, \cdot) + G 1 (3) Z (3, \cdot) + G 1 (4) Z (4, \cdot)

其中

G1(i) 表示

G1 的第

i 个元素，

Z(i,⋅) 表示

Z 的第

i 行。

在代码实现的时候，为了计算 Y(2,⋅)，初始化一个长度为 8 的浮点数行向量 T，令里面的值全等于零，然后用遍历行元素的方式进行如下计算

T T T T T = T + G 1 (0) Z (0, \cdot) = T + G 1 (1) Z (1, \cdot) = T + G 1 (2) Z (2, \cdot) = T + G 1 (3) Z (3, \cdot) = T + G 1 (4) Z (4, \cdot)

最后将

T 中的浮点数的值四舍五入赋值给

Y(2,⋅)。这样就避免了内存访问跨行的问题。注意，为了满足内存访问的局部性，增加了内存使用量，多用了

T。

对于边界行，按照镜像对称的方式选取相应行进行计算。比如，为了计算 Y(0,⋅)，初始化一个长度为 8 的浮点数行向量 T，令里面的值全等于零，然后用遍历行元素的方式进行如下计算

T T T T T = T + G 1 (0) Z (2, \cdot) = T + G 1 (1) Z (1, \cdot) = T + G 1 (2) Z (0, \cdot) = T + G 1 (3) Z (1, \cdot) = T + G 1 (4) Z (2, \cdot)

最后将

T 中的浮点数的值四舍五入赋值给

Y(0,⋅)。

扩展与总结

本文中所讲述的高斯模糊的计算方法，可以扩展到任意尺寸可分离核的滤波的实现。

设输入数据为 X，hX 行 wX 列，滤波核为 K，(2hK+1) 行 (2wK+1) 列，使用 K 对 X 进行二维滤波的结果是 Y。而直接采用二维循环的原始计算方法，需要进行 (2hK+1)×(2wK+1) 次乘法计算和 (2hK+1)×(2wK+1)−1 次加法计算。计算的时间复杂度是 O(wKhK) 的。

如果 K 是可分离核，可以写成列向量 Kvertical 和行向量 Kvertical 相乘的形式，即 K=Kvertical×Khorizontal。那么在计算滤波结果 Y 的时候，可以先用 Khorizontal 对 X 进行行滤波计算，将计算结果保存到 Z 中，计算 Z 中的每一个数值需要 (2wK+1) 次乘法计算和 2wK 次加法计算。再使用 Kvertical 对 Z 进行列滤波计算，得到最终结果 Y。在 Z 的基础上计算 Y 中的每一个数值需要 (2hK+1) 次乘法计算和 2hK 次加法计算。总的来说，根据 X 计算 Y 中的一个数值，需要进行 (2hK+2wK+2) 次乘法计算和 2hK+2wK 次加法计算。计算的时间复杂度从 O(wKhK) 降至 O(wK+hK)。

列滤波的过程还可以考虑内存访问的局部性原则，以提高程序的运行效率。

可分离核的实现方法和列滤波的内存访问加速的实现方法，都需要消耗额外的内存，用空间复杂度的提高换取时间复杂度和效率的改进。

0 0