---

time:2015年10月04日 星期日 00時00分27秒

opencv筆記4：模板運算和常見濾波操作

---

這一篇主要是學習模板運算，了解各種模板運算的運算過程和分類，理論方面主要參考《圖像工程——圖像處理》（章毓晉）一書第3章，空域增強：模板操作。同時也有個疑問：此書第四章，頻域圖像增強，講了低通濾波和高通濾波，然而這些東西和模板運算中的平滑、銳化操作有什么區別？。。。
以下是正文：

模板運算

首先我們把所有圖像看作矩陣。
模板一般是nxn（n通常是3、5、7、9等很小的奇數）的矩陣。模板運算基本思路：將原圖像中某個像素的值，作為它本身灰度值和其相鄰像素灰度值的函數。模板中有一個錨點（anchor point），通常是矩陣中心點，和原圖像中待計算點對應；整個模板對應的區域，就是原圖像中像素點的相鄰區域。模板也稱為核（kernel）。

前面的解釋翻譯成公式就是：

g(x,y)=function(f(x,y), template)

常見的function操作有卷積和排序兩種。卷積可以立即為一個map-reduce過程：元素對應相乘（mapper），乘積累加（reducer）。顯然，卷積是一個線性操作。
排序操作也不難理解：模板的錨點和待計算點綁定后，鄰域內所有點進行排序操作，將排序結果中符合策略規定的作為結果。一般的排序算法是O(n log n)的，不知道是不是因此有人認為模板排序運算不是線性的。其實通常處理的圖像像像素值都是unsigned char類型的，是[0,255]之間的非負整數，顯然用桶排序是可以O(n)復雜度內完成排序的，依我看也是一種線性運算。
如果模板排序前，需要對應元素和模板元素相乘，然后將乘積排序，那么這時候乘積可能是浮點數，排序就基本上是O(n log n)了，這確實是非線性操作了。不過目前我沒有見到類似的操作，也覺得沒有什么實際的用處。

濾波

模板運算的效果，可能讓圖像變好，也可能讓圖像變壞。我們當然需要好的那種模板運算了：）利用像素本身及其鄰域像素的灰度關系進行增強的方法，被稱為濾波，濾波使用到的模板就是濾波器。（注意：濾波器是一個模板矩陣，也就是核kernel，而具體的卷積操作還是排序操作，不是濾波器）

濾波和卷積的區別

卷積是濾波的一種實現方式。卷積是一種具體的運算，雖然它其實也是有點一種抽象的表述；而濾波則是比卷積要抽象的描述。

高頻分量和低頻分量

先看看頻率的本意：(狹義概念)頻率是單位時間內完成周期性變化的次數。推廣開來，（廣義概念）頻率就是指一定時間內的變化次數。
頻率在信號處理領域大量使用。信號處理中的函數自變量是時間；數字圖像處理被看作類似信號處理，只不過這里的函數自變量不再是時間，而是換成了圖像矩陣的像素灰度值。
原來在信號處理中，從前一秒到后一秒，信號周期性變化的次數，就是頻率；相應地，在數字圖像處理中，從一個像素點到相鄰的一個像素點，灰度值變化的多少，就是頻率。
所謂高頻分量，就是頻率值高，就是像素之間灰度變化大，這通常對應着圖像區域邊緣等；而低頻分量，就是頻率值低，就是像素灰度之間灰度變化小，這通常是圖像中穩定的區域，是在一個object的內部，同屬於一個superpixel...

總之，這樣的理解下，高頻分量對應圖像邊緣等像素變化大的像素點；低頻分量對應着圖像中穩定的區域。

平滑濾波和銳化濾波

平滑濾波能去除高頻分量，而銳化濾波能去除低頻分量。這么說還是抽象，具體講是：平滑濾波去處噪聲，銳化濾波強化邊緣、細節與周圍的對比度。
平滑濾波主要包括：線性平滑濾波（方框濾波、均值濾波、高斯濾波等）、非線性平滑濾波（中值濾波、序統計濾波）。opencv中對應boxblur、blur、gaussianblur函數。
銳化濾波主要包括：線性銳化濾波（拉普拉斯算子、高頻提升濾波）、非線性銳化濾波（基於梯度的銳化濾波、最大-最小銳化變換等）

注意 實際上對圖像進行二維傅立葉變換得到頻譜圖，就是圖像梯度的分布圖，當然頻譜圖上的各點與圖像上各點並不存在一一對應的關系，即使在不移頻的情況下也是沒有。傅立葉頻譜圖上我們看到的明暗不一的亮點，實際上圖像上某一點與鄰域點差異的強弱，即梯度的大小，也即該點的頻率的大小（可以這么理解，圖像中的低頻部分指低梯度的點，高頻部分相反）。一般來講，梯度大則該點的亮度強，否則該點亮度弱。這樣通過觀察傅立葉變換后的頻譜圖，也叫功率圖，我們首先就可以看出，圖像的能量分布，如果頻譜圖中暗的點數更多，那么實際圖像是比較柔和的（因為各點與鄰域差異都不大，梯度相對較小），反之，如果頻譜圖中亮的點數多，那么實際圖像一定是尖銳的，邊界分明且邊界兩邊像素差異較大的。

平滑濾波

線性平滑濾波

opencv現在有3個線性平滑濾波器：方框濾波、均值濾波、高斯濾波

先說均值濾波。均值濾波就是用指定大小的、元素全為1的模糊核，對原圖進行卷積操作（其實，就是原圖像中當前位置對應的核大小的區域，各個元素相加），然后除以核的元素個數。也可以理解為，其核為：元素全為1、系數為元素個數的矩陣。

然后是方框濾波。是均值濾波的推廣：元素全為1、系數為alpha的矩陣。alpha等於核的個數時，就是均值濾波；否則一般取1（哦，為什么要搞這么大？難道不會超出255麻？）

再看高斯濾波。前面兩個濾波的核（也叫模板），元素值都只有一種。如果模板的元素不只一種，就是加權線性濾波了。
高斯濾波是加權線性濾波的一種，准確說是：模板元素的分布符合二次高斯分布。高斯分布其實就是正態分布。因為一般認為噪聲的分布都符合高斯分布，那么去噪也用符合高斯分布的模板，效果會比較好。
通過為GaussianBlur函數傳入sigmaX,sigmaY,size等參數，函數能生成相應的符合高斯分布的模板。然后對原圖像和模板進行卷積操作，就得到濾波后的圖像。

非線性平滑濾波

opencv現在有2個非線性平滑濾波：中值濾波和雙邊濾波

中值濾波：模板限定區域內，取像素灰度的中值（我理解為中位數），作為計算結果。中值濾波的效果是，讓與周圍像素灰度值的差比較大的像素改取與周圍像素值接近的值，消除了信號序列（這里是模板框定范圍內像素點灰度值）的孤立點。因為不是簡單的取均值，產生的模糊更少些，通常能比均值濾波更好地保持圖像的細節。

測試發現，在有白色噪聲的圖像上，中值濾波和均值濾波的對比效果很明顯：

#include <iostream>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

int main(){
    //【1】定義原圖像
    //【2】調用模糊函數
    //【3】顯示結果

    //【1】定義原圖像
    Mat srcImage = imread("/home/chris/workspace/clion/blur_img1.png");
    Mat meanBlurImage, medianBlurImage;

    //【2】調用模糊函數
    blur(srcImage, meanBlurImage, Size(3,3));
    medianBlur(srcImage, medianBlurImage, 3);


    //【3】顯示結果
    imshow("原圖", srcImage);
    imshow("均值濾波", meanBlurImage);
    imshow("中值濾波", medianBlurImage);
    while(waitKey(1)!='q'){}
    destroyAllWindows();

    return 0;
}

雙邊濾波：
雙邊濾波是一種簡單的、非迭代的保邊平滑過濾器：能夠去除圖像噪聲，同時很好地保持邊界。缺點是比其他過濾器慢。
查看opencv官方文檔，目前（opencv3.0）中，雙邊濾波器函數的實現依然有問題，“This filter does not work inplace.”
看起來有點沮喪，不管了，先了解下原理。

雙邊濾波同時考慮了空間域和值域的差別：空間域給人一種“出身”的感覺，模板框定了你周圍的像素點，這些點不管它們灰度值是多少，你總要按相應權重對待它們（按模板中對應元素值來處理），比如均值濾波是“一視同仁”，高斯濾波是“像沖擊波一樣從自身衰減”。空間域濾波器的效果是，能去除噪聲。
而從值域的角度看，給人一種“看后天努力程度”的感覺：對於模板框定的周圍像素點，考慮它們的灰度值，而不去官它們當中的“老幼尊卑”。這方面的代表是alpha-截尾均值濾波器。值域濾波器的效果是，能保留邊界效果。

雙邊濾波綜合考慮了空間域和值域，其計算公式中的權重系統，是定義域核與值域核的乘積。

銳化濾波

主要包括：線性銳化濾波（拉普拉斯算子、高頻提升濾波）、非線性銳化濾波（基於梯度的銳化濾波、最大-最小銳化變換等）
線性銳化濾波可以借助模板卷積實現。對應積分運算的模板卷積可以平滑圖像，反過來對應微分運算的模板卷積可以銳化圖像。銳化模板系數的取值，應該在中心為正而周圍遠離中心處為負。

當然，圖像銳化還可以用高通濾波法來做，不過不屬於模板操作的范圍，現在還不懂，以后再說。

非線性銳化濾波

這次先來看看非線性銳化濾波。所謂線性還是非線性，是從最后的結果來看，計算步驟是不是線性的：雖然拉普拉斯算子是二階差分得到的，但是結果上開它等同於做線性模板卷積運算；而Sobel算子等一階差分方法，因為要分別考慮x、y兩個方向然后再合並，整個步驟沒法簡化，所以是非線性的。從推導的角度看，要先看一階差分操作，也就是非線性的幾個算子。

梯度銳化

梯度銳化依然是一種模板算法卷積算法，經過一系列推導，並整理出對應的運算模板，就是我們最終需要的。
圖像一般是二維矩陣，因此梯度銳化法一般在x、y方向分別計算出梯度幅值Gx、Gy，然后再合並。
用差分來近似微分.比如水平垂直差分:

f'(x) = f(i,j)-f(i+1,j)
f'(y) = f(i,j)-f(i,j+1)

或者交叉差分：

f'(x) = f(i,j)-f(i+1, j+1)
f'(y) = f(i+1,j ) - f(i, j+1) 

x、y兩個方向上差分的結果，可以通過距離公式和在一起，用來表示最后的計算結果，比如采用水平垂直差分+曼哈頓距離公式，有：

g(i,j)=|f(i,j)-f(i+1,j)|+|f(i,j)-f(i,j+1)|

而使用交叉差分+歐幾里得距離公式，有：

g(i,j)=sqrt( (f(i,j)-f(i+1,j+1))^2 + (f(i+1,j)-f(i,j+1)^2 )

或者使用交叉差分+曼哈頓距離，有：

g(i,j)=|(f(i,j)-f(i+1,j+1)| + |(f(i+1,j)-f(i,j+1)|

上式就是Roberts算子的梯度幅值公式，對應的模板為：

Sx=[1  0]
   [0 -1]

Sy=[0 -1]
   [1  0]
   
Gx=Sx * A   #卷積操作
Gy=Sy * A   #卷積操作

梯度銳化的改進 上述算法計算出來的，是一個考慮了周遭像素的、差分、然后相加的結果。這個結果拿來取代原來的像素值，不一定合適，我們可以對此加以判斷，如果它落在某個閾值范圍內（比如大於閾值T），才算作有效，否則仍然取原來的f(i,j)灰度值。

梯度銳化的不足 考慮低頻區域中的一個點f(i,j)，其g(i,j)通常是0；再噪聲點f(i,j)，顯然它對應的g(i,j)不為0，而且還比較大。往往噪聲點在梯度銳化算法中，被增強的效果比普通點更大。因此使用此算法前盡量去除噪聲。

梯隊銳化對應的模板？ 顯然上面提到的公式，對應的模板是2x2的。這和通常使用的奇階方陣不一樣，是不實用的。無論是x方向還是y方向，2x2的方陣都僅僅是考慮了半鄰域，而不是整個鄰域。不過這個思路是可取的。Sobel算子就是考慮了整個鄰域的一階算子，當然還可以用拉普拉斯這樣的二階算子。

Prewitt算子

最簡單的、考慮了整個鄰域的算子，是Prewitt算子，其卷積模板為：

   [-1 0 1]
Sx=[-1 0 1]
   [-1 0 1]
   
   [ 1  1  1]
Sy=[ 0  0  0]
   [-1 -1 -1]

顯然，它僅僅在一個方向上考慮了權值，而另一個方向上則沒有考慮。Sobel算法則是其的進一步優化。

Soble和Canny

Sobel算子用來計算圖像的導數，目的是獲得圖像邊緣。因此，它常常被用於邊緣檢測。Canny算法是比Sobel更完整的邊緣檢測算法，包括了預處理（高斯濾波）、后處理（閾值法去除非邊緣點）。
為什么用導數檢測邊緣？考慮一維圖像f(x)，其邊緣是一個點。作為邊緣點，其左右兩側像素灰度的變化濾肯定不一樣（否則，就是真的平滑圖像了，哪還有什么邊界）。如果計算f(x)的導數，f(x)一定是導數中的極值。推廣到二維圖像這也是成立的。只不過求導公式使用差分公式來近似代替了。Sobel和前面的銳化梯度的區別在於，它在x、y方向分別使用一個模板進行卷積運算，得到Gx、Gy兩個分量，然后用G=sqrt(Gx^2 + Gy^2)得到梯度幅值。

Sobel算子，一般取模板:

   [-1 0 1]
Sx=[-2 0 2]
   [-1 0 1]
   
   [ 1  2  1]
Sy=[ 0  0  0]
   [-1 -2 -1]

Canny是指Canny邊緣檢測算法，其步驟包括：

1. 濾波：讀取灰度圖像后，使用高斯濾波做平滑處理（去噪）
2. 增強：用一階偏導的有限差分來計算梯度的幅值和方向（銳化）。這一步可以使用Roberts、Sobel、Prewitt等算子,或者說，使用類似Sobel濾波器的濾波步驟。
3. 檢測：使用閾值法將非邊緣點去除，獲得真正的邊緣點
   如此看來，Canny邊緣檢測是模板操作的綜合應用了，既有平滑處理，也有銳化處理。

在Sobel算子中，容易計算：

Gx(x,y)=Sx * A   #卷積操作
Gy(x,y)=Sy * A   #卷積操作
G(x,y)=sqrt(Gx^2+Gy^2)   #梯度幅值
theta=arctan(Gy(x,y)/Gx(x,y))   #梯度幅角

其中，梯度幅角在Canny最后一步的檢測判斷中使用到。

P.S.:opencv中的Sobel函數，如果指定的模板規格為3x3，按照前面的模板矩陣，會產生明顯的不精確的結果。通過調用Scharr算子的模板，效果會好些：

   [-3  0  3]
Sx=[-10 0 10]
   [-3  0  3]
   
   [-3 -10 -3]
Sy=[ 0  0   0]
   [-3 -10 -3]

線性銳化濾波

前面使用的是一階差分，現在使用二階差分，同樣能得到用來銳化濾波的操作。只不過很巧妙的是，拉普拉斯算子的推導過程到最后發現，能夠等價於一個權值模板卷積操作，所以是一階操作。

拉普拉斯算子

拉普拉斯算子是二階微分算子，也用於線性銳化濾波:

g(i,j)=f''(x)+f''(y)

仍然使用差分的方式，容易得到：

f''(x) = 2f(i,j)-f(i+1,j)-f(i-1,j)
f''(y) = 2f(i,j)-f(i,j+1)-f(i,j-1)
g(i,j)=4f(i,j)-f(i+1,j)-f(i-1,j)-f(i,j+1)-f(i,j-1)

對應的模板是：

 0 -1  0
-1  4 -1
 0 -1  0

這是只考慮4鄰域的情況。如果考慮8鄰域，對應的模板是：

-1 -1 -1
-1  8 -1
-1 -1 -1

高頻提升濾波

用原始圖像見去平滑或模糊圖像能得到非銳化掩模，將非銳化掩模加到原始圖像上能得到銳化圖像。更進一步，可以把原始圖像放大A倍后再減去平滑圖像：

h(x,y) = A*f(x,y)-g(x,y) = (A-1)*f(x,y)+(f(x,y)-g(x,y)) = (A-1)*f(x,y) + mask(x,y)

其實就是稍微復雜點的線性組合了。

自定義濾波

如果想到了什么新的算法，或者純粹想試一試碰碰運氣，可以自定義模板矩陣，扔給opencv的filter2d函數，就可以看到效果了。

ref

RachealZhang:雙邊濾波器的原理及實現
淺墨:【OpenCV入門教程之九】 非線性濾波專場：中值濾波、雙邊濾波
MrMystery:圖像增強-圖像銳化
奮斗斌斌的專欄:Canny邊緣檢測算法原理及其VC實現詳解(一)