徹底理解數字圖像處理中的卷積-以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 來用動畫過程模擬圖像卷積的計算過程.