彻底理解数字图像处理中的卷积-以Sobel算子为例

转自：http://www.cnblogs.com/freeblues/p/5738987.html

概述

卷积在信号处理领域有极其广泛的应用, 也有严格的物理和数学定义. 本文只讨论卷积在数字图像处理中的应用.

在数字图像处理中, 有一种基本的处理方法:线性滤波. 待处理的平面数字图像可被看做一个大矩阵, 图像的每个像素对应着矩阵的每个元素, 假设我们平面的分辨率是 1024*768, 那么对应的大矩阵的行数= 1024, 列数=768.

用于滤波的是一个滤波器小矩阵(也叫卷积核), 滤波器小矩阵一般是个方阵, 也就是 行数 和 列数 相同, 比如常见的用于边缘检测的 Sobel 算子 就是两个 3*3 的小矩阵.

进行滤波就是对于大矩阵中的每个像素, 计算它周围像素和滤波器矩阵对应位置元素的乘积, 然后把结果相加到一起, 最终得到的值就作为该像素的新值, 这样就完成了一次滤波.

上面的处理过程可以参考这个示意图:

图像卷积计算示意图:

对图像大矩阵和滤波小矩阵对应位置元素相乘再求和的操作就叫卷积(Convolution)或协相关(Correlation).

协相关(Correlation)和卷积(Convolution)很类似, 两者唯一的差别就是卷积在计算前需要翻转卷积核, 而协相关则不需要翻转.

以 Sobel 算子为例

Sobel 算子 也叫 Sobel 滤波, 是两个 3*3 的矩阵, 主要用来计算图像中某一点在横向/纵向上的梯度, 看了不少网络上讲解 Sobel 算子 的文章, 发现人们常常把它的横向梯度矩阵和纵向梯度矩阵混淆. 这可能与 Sobel 算子 在它的两个主要应用场景中的不同用法有关.

Sobel 算子的两个梯度矩阵: Gx 和 Gy

这里以 Wiki 资料为准, Sobel 算子 有两个滤波矩阵: Gx 和 Gy, Gx 用来计算横向的梯度, Gy 用来计算纵向的梯度, 下图就是具体的滤波器:

• 注意:这里列出的这两个梯度矩阵对应于横向从左到右, 纵向从上到下的坐标轴, 也就是这种:

原点
O ------->  x轴
|
|
|
V  y轴

Sobel 算子的用途

它可以用来对图像进行边缘检测, 或者用来计算某个像素点的法线向量. 这里需要注意的是:

• 边缘检测时: Gx 用于检测纵向边缘, Gy 用于检测横向边缘.
• 计算法线时: Gx 用于计算法线的横向偏移, Gy 用于计算法线的纵向偏移.

计算展开

假设待处理图像的某个像素点周围的像素如下:

左上	上	右上
左	中心像素	右
左下	下	右下

那么用 Gx 计算展开为:

横向新值 = (-1)*[左上] + (-2)*[左] + (-1)*[左下] + 1*[右上] + 2*[右] + 1*[右下]

用 Gy 计算展开为:

纵向新值 = (-1)*[左上] + (-2)*[上] + (-1)*[右] + 1*[左下] + 2*[下] + 1*[右下]

前面说过, 做图像卷积时需要翻转卷积核, 但是我们上面的计算过程没有显式翻转, 这是因为 Sobel 算子 绕中心元素旋转 180 度后跟原来一样. 不过有些 卷积核 翻转后就变了, 下面我们详细说明如何翻转卷积核.

卷积核翻转

前面说过, 图像卷积计算, 需要先翻转卷积核, 也就是绕卷积核中心旋转 180度, 也可以分别沿两条对角线翻转两次, 还可以同时翻转行和列, 这3种处理都可以得到同样的结果.

对于第一种卷积核翻转方法, 一个简单的演示方法是把卷积核写在一张纸上, 用笔尖固定住中心元素, 旋转 180 度, 就看到翻转后的卷积核了.

下面演示后两种翻转方法, 示例如下:

假设原始卷积核为:

a	b	c
d	e	f
g	h	i

方法2:沿两条对角线分别翻转两次

先沿左下角到右上角的对角线翻转, 也就是 a和i, b和f, d和h交换位置, 结果为:

i	f	c
h	e	b
g	d	a

再沿左上角到右下角的对角线翻转, 最终用于计算的卷积核为:

i	h	g
f	e	d
c	b	a

方法3:同时翻转行和列

在 Wiki 中对这种翻转的描述:

convolution is the process of flipping both the rows and columns of the kernel and then multiplying locationally similar entries and summing.

也是把卷积核的行列同时翻转, 我们可以先翻转行, 把 a b c跟 g h i 互换位置, 结果为:

g	h	i
d	e	f
a	b	c

再翻转列, 把 g d a 和 i f c 互换位置, 结果为:

i	h	g
f	e	d
c	b	a

在 Wiki 中有一个计算展开式, 也说明了这种翻转:

• 注意:这里要跟矩阵乘法区分开, 这里只是借用了矩阵符号, 实际做的是对应项相乘, 再求和.

图像边缘像素的处理

以上都默认待处理的像素点周围都有像素, 但是实际上图像边缘的像素点周围的像素就不完整, 比如顶部的像素在它上方就没有像素点了, 而图像的四个角的像素点的相邻像素更少, 我们以一个图像矩阵为例:

左上角	...		...	右上角
...	...	...	...	...
左侧	...	...	...	右侧
...	...	...	...	...
左下角	...		...	右下角

位于左上角的像素点的周围就只有右侧和下方有相邻像素, 遇到这种情况, 就需要补全它所缺少的相邻像素, 具体补全方法请参考下一节的代码.

用GPU进行图像卷积

如果在 CPU 上实现图像卷积算法需要进行4重循环, 效率比较差, 所以我们试着把这些卷积计算放到 GPU 上, 用 shader 实现, 结果发现性能相当好, 而且因为顶点着色器和片段着色器 本质就是一个循环结构, 我们甚至不需要显式的循环, 代码也清晰了很多.

图像卷积在代码中的实际应用, 下面是一个 GLSL 形式的着色器, 它可以根据纹理贴图生成对应的法线图:

-- 用 sobel 算子生成法线图    generate normal map with sobel operator
genNormal1 = {
vertexShader = [[
attribute vec4 position;
attribute vec4 color;
attribute vec2 texCoord;

varying vec2 vTexCoord;
varying vec4 vColor;
varying vec4 vPosition;

uniform mat4 modelViewProjection;

void main()
{
    vColor = color;
    vTexCoord = texCoord;
    vPosition = position;
    gl_Position = modelViewProjection * position;
}
]],

fragmentShader = [[
precision highp float;

varying vec2 vTexCoord;
varying vec4 vColor;
varying vec4 vPosition;

// 纹理贴图
uniform sampler2D tex;
uniform sampler2D texture;

//图像横向长度-宽度, 图像纵向长度-高度
uniform float w;
uniform float h;

float clamp1(float, float);
float intensity(vec4);

float clamp1(float pX, float pMax) {
    if (pX > pMax) 
        return pMax;
    else if (pX < 0.0)
        return 0.0;
    else
        return pX;   
}

float intensity(vec4 col) {
    // 计算像素点的灰度值
    return 0.3*col.x + 0.59*col.y + 0.11*col.z;
}

void main() {
    // 横向步长-每像素点宽度，纵向步长-每像素点高度
    float ws = 1.0/w ;
    float hs = 1.0/h ;
    float c[10];
    vec2 p = vTexCoord;
    lowp vec4 col = texture2D( texture, p );

    // sobel operator
    // position.      Gx.            Gy
    // 1 2 3     |-1. 0. 1.|   |-1. -2. -1.|
    // 4 5 6     |-2. 0. 2.|   | 0.  0.  0.|
    // 7 8 9     |-1. 0. 1.|   | 1.  2.  1.|
    // 右上角，右，右下角
c[3] = intensity(texture2D( texture, vec2(clamp(p.x+ws,0.,w), clamp(p.y+hs,0.,h) )));
c[6] = intensity(texture2D( texture, vec2(clamp1(p.x+ws,w), clamp1(p.y,h))));
c[9] = intensity(texture2D( texture, vec2(clamp1(p.x+ws,w), clamp1(p.y-hs,h))));

// 上, 下
c[2] = intensity(texture2D( texture, vec2(clamp1(p.x,w), clamp1(p.y+hs,h))));
c[8] = intensity(texture2D( texture, vec2(clamp1(p.x,w), clamp1(p.y-hs,h))));

// 左上角, 左, 左下角
c[1] = intensity(texture2D( texture, vec2(clamp1(p.x-ws,w), clamp1(p.y+hs,h))));
c[4] = intensity(texture2D( texture, vec2(clamp1(p.x-ws,w), clamp1(p.y,h)))); 
c[7] = intensity(texture2D( texture, vec2(clamp1(p.x-ws,w), clamp1(p.y-hs,h))));

    // 先进行 sobel 滤波, 再把范围从 [-1,1] 调整到 [0,1]
    // 注意: 比较方向要跟坐标轴方向一致, 横向从左到右, 纵向从下到上
    float dx = (c[3]+2.*c[6]+c[9]-(c[1]+2.*c[4]+c[7]) + 1.0) / 2.0;
    float dy = (c[7]+2.*c[8]+c[9]-(c[1]+2.*c[2]+c[3]) + 1.0) / 2.0;
    float dz = (1.0 + 1.0) / 2.0;

    gl_FragColor = vec4(vec3(dx,dy,dz), col.a);

}
]]
}

后续有时间的话考虑写一个 APP 来用动画过程模拟图像卷积的计算过程.