完全理解數字圖像處理中的卷積-以Sobel算子爲例

完全理解數字圖像處理中的卷積-以Sobel算子爲例

• 做者:FreeBlues
• 修訂記錄
  • 2016.08.04 初稿完成

概述

卷積在信號處理領域有極其普遍的應用, 也有嚴格的物理和數學定義. 本文只討論卷積在數字圖像處理中的應用.

在數字圖像處理中, 有一種基本的處理方法:線性濾波. 待處理的平面數字圖像可被看作一個大矩陣, 圖像的每一個像素對應着矩陣的每一個元素, 假設咱們平面的分辨率是 1024*768, 那麼對應的大矩陣的行數= 1024, 列數=768.

用於濾波的是一個濾波器小矩陣(也叫卷積核), 濾波器小矩陣通常是個方陣, 也就是 行數 和 列數 相同, 好比常見的用於邊緣檢測的 Sobel 算子 就是兩個 3*3 的小矩陣.

進行濾波就是對於大矩陣中的每一個像素, 計算它周圍像素和濾波器矩陣對應位置元素的乘積, 而後把結果相加到一塊兒, 最終獲得的值就做爲該像素的新值, 這樣就完成了一次濾波.

上面的處理過程能夠參考這個示意圖:

圖像卷積計算示意圖:

對圖像大矩陣和濾波小矩陣對應位置元素相乘再求和的操做就叫卷積(Convolution)或協相關(Correlation).

協相關(Correlation)和卷積(Convolution)很相似, 二者惟一的差異就是卷積在計算前須要翻轉卷積核, 而協相關則不須要翻轉.

以 Sobel 算子爲例

Sobel 算子 也叫 Sobel 濾波, 是兩個 3*3 的矩陣, 主要用來計算圖像中某一點在橫向/縱向上的梯度, 看了很多網絡上講解 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 來用動畫過程模擬圖像卷積的計算過程.

參考

圖像卷積與濾波的一些知識點
Sobel Derivatives
Wiki:Kernel (image processing)