2016年8月20日 星期六

暑假作業3:簡單的filter實作以及自由發揮

接下來就是正式進入影像處理的部分,影像簡單的來看可以看成是一個矩陣,矩陣中的數字就是每個點的顏色,我們最基礎的影像處理方式是濾波器,也就是對矩陣做矩陣乘法,接下來就來看一個最簡單的濾波器。





1.mean filter

均值濾波器,這個濾波效果可以讓雜訊消去,至於為什麼可以把雜訊消去呢?先不要看影像我們來看一個簡單的波,下面是一個有雜訊的波以及去除雜訊處理後的波



他是怎麼做到的呢?基本上雜訊的頻率都很高,基本會會都高於原本正常的波,這種高頻雜訊,我們可以利用這個特性,對波上的每一個點取周遭點做平均,那麼高頻率的雜訊被影響就非常的大,數字就被周圍的點改變非常多,但是相對低頻率的正常的波不會被影響那麼多(還是會),如此以來就可以去除雜訊變成下面正常的波形了,這個在聲音處理非常有用,因為聲音是個低頻率的波,但是周遭的雜訊(收音儀器電磁波干擾等等)頻率幾乎是聲音的數百倍,所以做平均化的濾波處理去除雜訊的效果會非常好,所以雜訊與原波的頻率差越多,均值濾波的效果會更好。

那回到影像處理的部分,我們只要改成對一個像素取周圍半徑R裡面的點加起來取平均賦值就是影像中的mean filter在做的事了,均值的意思就是每個點的權重都一樣,影像的問題就在這了,如果雜訊像聲音這種東西一樣頻率是數百倍的差距就好了,但是影像本身的雜訊點出現的頻率通常不會很高,所以做mean filter的缺點就是我們影響到原本的影像程度也會比較高,這個呈現的結果就會是,影像會變得模糊,下圖是高頻率的雜訊做mean filter的效果會非常好。



基本的理論部分,我們已經懂了,接下來就是要怎麼變成程式碼。

最簡單的,我們可以用基本巢狀迴圈

filterSize代表filter矩陣的大小,取半徑Diameter


Mat mat = imread("image.jpg", CV_LOAD_IMAGE_COLOR);   // Read the file
int filterSize = 5;
int Diameter = filterSize/2;
Mat final_mat;
final_mat.create(height, width, CV_8UC3);
 //cin >> Diameter;
 for (int i = Diameter; i < height - Diameter; i++)
 {
  for (int j = Diameter; j < width - Diameter; j++)
  {
   Vec3f sum = (0.0f,0.0f,0.0f);
   for (int p = i - Diameter; p <= i + Diameter; p++)
   {
    for (int q = j - Diameter; q <= j + Diameter; q++)
    {
     sum += mat.at(p, q);
    }
   }
   sum /= FilterSize*FilterSize;
   final_mat.at(i, j) = sum;
  }
 }


mat.at(p,g)是8UC3的未處理原影像,每次loop後的sum結果會存在final_mat中,這個例子是5X5的結果,所以就是對點(i,j)周圍半徑內一共25個點加起來除以25做為新的值。
結果的影像如下,上方為原圖,下面就是mean filter的結果


最後來提一下影像的頻率低是什麼意思,以1000hz的波來看的話,代表一張圖至少要2000*2000的像素才夠顯示一個1000Hz的雜訊,這就是有名的奈奎斯特頻率(Nyquist frequence),例如,CD音訊訊號的採樣頻率為44100Hz,那麼它的奈奎斯特頻率就是22050 Hz,這是CD音訊數據所能表現的最高頻率。

一般的影像幾乎很少有這麼高的像素,所以也就很少看得到用mean filter之後既能去除雜訊又不會讓基本影像被模糊化。

既然有這個缺點,mean filter有一個改良版本,就是接下來要看的gaussian filter


2.gaussian filter
高斯濾波器不是高斯發明的,它的名稱由來是因為他用到了一個最重要的分布曲線—-高斯分布,也就是所謂的常態分佈
他的數學公式就是自然常數e的指數函數型態
之前提到,影像被濾波過後會出現嚴重的模糊現象,為了改善這個問題,想要做到「模糊又不會太模糊」,基本還是要從取周遭的點來做加總平均,只是呢,這次的平均我們讓原本的像素權重高一點來保留原本的樣子,減少模糊的程度。
而這個權重,剛好我們可以使用高斯函數作為係數,可以看到在中心點的係數特別高,遠方的點系數幾乎等於0,這樣不就達成了去除雜訊又能保留原本影像輪廓的目的了!

數學公式知道了,要怎麼轉化成程式碼呢?基本上還是用基本巢狀迴圈,接下來就是發揮這個演算法的精髓,重點只有一個,中心點的權重越大越好,高斯分布有一個離散的逼近版本,那就是接下來要使用的二項式濾波器
我們都知道二項式係數的特徵-巴斯卡三角形
我們取一個奇數行的係數來看
就用第五行的【1    4     6     4       1】我們把它放到一個5X5矩陣中,先放在第一列中
1
4
6
4
1





















接下來第二列也是【1    4     6     4       1】但是要乘上【1    4     6     4       1】中第二個數字的「4」
1
4
6
4
1
4
16
24
16
4
















接下來第三列也是【1    4     6     4       1】但是要乘上【1    4     6     4       1】中的第三個數字「6」
1
4
6
4
1
4
16
24
16
4
6
24
36
24
6











接下來第四列也是【1    4     6     4       1】但是要乘上【1    4     6     4       1】中第四個數字的「4」
1
4
6
4
1
4
16
24
16
4
6
24
36
24
6
4
16
24
16
4





接下來第五列也是【1    4     6     4       1】但是要乘上【1    4     6     4       1】中第五個數字的「1」
1
4
6
4
1
4
16
24
16
4
6
24
36
24
6
4
16
24
16
4
1
4
6
4
1
有發現到雖然是用二項式生成的矩陣,但是中心點的係數是最大的,更重要的是越邊緣的點與中心點的差距會更大,感覺得到與鐘形分布的高斯分布有那麼幾分相似對吧,用這個方法就不用去套用那複雜的公式,就可以很簡單的生成一個足夠逼近高斯矩陣效果的二項式濾波器了,想到這個方法的人真的是很聰明。所以接下來我們不是用gaussian filter而是較為簡化版本的二項式濾波器--binomial filter。

以這個規則來看7X7的binomail filter就會是這樣
      1,    6,      15,    20,     15,     6,     1,
      6,   36,     90,  120,     90,   36,     6,
    15,   90,   225,  300,   225,   90,   15,
    20,  120,  300,  400,   300, 120,   20,
    15,   90,   225,  300,   225,   90,   15,
      6,   36,     90,  120,     90,   36,     6,
      1,    6,      15,    20,     15,     6,     1,
所以生成更大的filter係數只要照這個規則算就可以了。

只是要注意,這裡的係數並沒有「單位化」,簡單的說我們要再除以總和256,這樣系數才是我們要的

矩陣生成完了,那麼我們的程式碼只要就簡單地改用這個矩陣來乘,就完成了binomail filter的濾波效果了,非常簡單。



Mat mat = imread("image.jpg", CV_LOAD_IMAGE_COLOR); // Read the file

int filterSize = 5;

int Diameter = filterSize/2;

Mat final_mat;
final_mat.create(height, width, CV_8UC3);

int  mask[5][5]={

      1,4,6,4,1,

      4,16,24,16,4,

      6,24,36,24,6,

      4,16,24,16,4,

      1,4,6,4,1

};


for (int i = Diameter; i < height - Diameter; i++)

 {

  for (int j = Diameter; j < width - Diameter; j++)

  {

   int m = 0;

   Vec3f sum = (0.0f,0.0f,0.0f);

   for (int p = i - Diameter; p <= i + Diameter; p++)

   {

    int n = 0;

    for (int q = j - Diameter; q <= j + Diameter; q++)

    {

     sum = sum + mat.at<Vec3b>(p, q) *  mask[m][n];

     n++;

    }

    m++;

   }

   final_mat.at<Vec3b>(i, j) = sum/256;

  }

 }

 



結果如下圖,基本上就是不怎麼模糊的影像
這樣看可能看不太出來,所以我們把原圖還有mean filter一起拿來比較
由上至下:原圖、mean filter、binomial filter

雖然binomail filter解決了mean filter的問題,但是基本上他還是在作模糊,只是我們利用高斯分布來設定權重讓模糊的程度沒那麼大而已,最讓我們感受強烈的地方就是邊緣部分,像是鬍鬚和背景交會處,他們不應該模糊在一起的,可以發現就算用了binomial filter鬍鬚還是幾乎因為模糊效果幾乎不見了,我們想要在模糊的同時保留這些邊緣處的像素,為了解決這個問題,我們接下來要介紹一個更加進階的filter--bilateral filter

3.bilateral filter
這個filter,中文叫做雙邊濾波器,顧名思義他有兩種不同的濾波組合而成,一樣是使用權重來達到濾波效果,只是他的重點在,邊緣處的權重要特別大來保留下來,而且非邊緣的地方還是要作模糊化,所以才會使用兩種濾波的組合。
接下來一樣從數學原理開始看起,最後在利用他的數學定義的精神來轉化為程式碼實作。
先從一般的係數矩陣的數學形式來看
x是要處理的點座標,h()是output,f()是input,c()是係數functuion,希臘字母Epsilon代表一個range內的維度,這裡就可以看成「選定一個點x,他會等於周圍範圍內的所有點epsilon與係數函數相乘後加起來的結果」,當然這裡的點是連續的,就會是積分。
那Kd是什麼呢?
可以看到它是係數function的總和,所以前面乘了一個Kd^-1代表我們單位化
這些是不是都很眼熟
c()就是這個矩陣
1
4
6
4
1
4
16
24
16
4
6
24
36
24
6
4
16
24
16
4
1
4
6
4
1
Kd就是上面的總和256。
f()就是一張未處理的圖,h()就是結果圖。

這些數學公式有了意義,我們也看得懂了,所以接下就進入正題
考慮了上述的數學式,我們有一個s()作為similarity function相似度函數
有以下的關係

Kr則一樣為s()的積分



好像根本沒有不同的地方,但是注意s()中不是x而是一個新的算式—-f(epsilon)-f(x),這代表,我們這次不是用位置(epsilon-x)來算系數矩陣了,而是用像素的差異(f(epsilon)-f(x))來計算系數矩陣,這就是這個演算法的精髓,為什麼要用像素之間的差異來計算系數矩陣呢?
回到一開始的目標吧,「保留邊緣」,那第一步對於程式碼來說,甚麼是邊緣呢?馬上就想到「顏色差異大的地方」,例如白色牆上的人影,邊緣處的顏色差異最大,於是我們馬上就有了計算s()函數的概念了,與周遭的像素差異越大代表這個點屬於邊緣,那麼給他大一點的權重保留他,反之與周遭像素差異並不大的點,不屬於邊緣,給他小一點的值,如此一來邊緣的權重矩陣就生成了。
融合以上兩者c()與s()就是這次介紹的bilateral filter
我們只要把他們用乘法結合在一起就可以了,所以單位化用的k就長這個樣子

每個係數相乘之後的總和,讓我們在計算完係數矩陣後能夠單位化矩陣。
知道了數學公式的背後意義,那麼我們就來進行程式的實作。
c()非常簡單,我們用之前生成的binomial filter就可以了,代表了模糊化用的矩陣,接下是重點s()相似度矩陣,生成的方法就用直觀的方式--巢狀迴圈
每次對點(i,j)往上下左右一個矩陣半徑的長度範圍內的所有周圍點(p,q)取RGB的相差
RGB的差,我們用簡單的歐式距離,空間中的兩點距離等於兩點的(x,y,z)相減平方後加起來開根號來計算。
for (int i = Diameter; i < height - Diameter; i++) { for (int j = Diameter; j < width - Diameter; j++) { int m = 0; for (int p = i - Diameter; p <= i + Diameter; p++) { int n = 0; for (int q = j - Diameter; q <= j + Diameter; q++) { float diff = 0.0f; diff += pow((mat.at<Vec3b>(p, q)[0] - mat.at<Vec3b>(i, j)[0]),2); diff += pow((mat.at<Vec3b>(p, q)[1] - mat.at<Vec3b>(i, j)[1]), 2); diff += pow((mat.at<Vec3b>(p, q)[2] - mat.at<Vec3b>(i, j)[2]), 2); diff = sqrt(o_diff); s[m][n] = diff; n++; } m++; } } }
如此一來我們就可以生成s矩陣了,但是這裡有一個大陷阱,那就是RGB的差越大權重越大這點,我們來看看這個s裡面放的會是甚麼東西,s[m][n]代表著選定某個像素(i,j)他的周圍點(m,n)與他的顏色差異,那麼有一個點一定會是0對吧?沒錯,正是中心點,因為她就是本身那個點,所以(i,j)跟(i,j)的顏色差異當然會等於0,那麼不管c矩陣數字是什麼,在矩陣相乘的結果下本身的點權重變成了0,那麼二項式濾波的模糊效果不就完全不是我們要的嗎?這讓我們有個疑問,我們不是要算顏色差異最大的權重越高嗎?怎麼反而讓本身的像素消失了,所以有兩種方法,一個是全部都加上1,這個缺點就在於會讓差距的倍率變小可能邊緣保留效果會變很差,第二個就是取一個門檻值,差異度高於這個門檻的權重設為2,低於它的就設為1。


for (int i = Diameter; i < height - Diameter; i++)
{
 for (int j = Diameter; j < width - Diameter; j++)
 {
  int m = 0;
  for (int p = i - Diameter; p <= i + Diameter; p++)
  {
   int n = 0;
   for (int q = j - Diameter; q <= j + Diameter; q++)
   {
    float diff = 0.0f;
    diff += pow((mat.at<Vec3b>(p, q)[0] - mat.at<Vec3b>(i, j)[0]),2);
    diff += pow((mat.at<Vec3b>(p, q)[1] - mat.at<Vec3b>(i, j)[1]), 2);
    diff += pow((mat.at<Vec3b>(p, q)[2] - mat.at<Vec3b>(i, j)[2]), 2);
    
                                diff = sqrt(o_diff);
                                diff = diff > 0.25? 2:0;
    s[m][n] = diff;
    n++;
   }
  m++;
  }
 } 
} 
我們加入diff = diff > 0.25? 2:0; 來取得非0的係數。
這樣的矩陣加乘後別忘了還要在做一次binomail filter去平滑模糊之後用這兩個矩陣系數的總和單位化,如此一來就完成了baliteral filter
我們來看看結果
影像根本沒有產生變化,原因很簡單,不要忘了每一個像素的值都是0~255,所以平方後開根號基本上會是大於零的整數,不會有點小於0.25--中心點例外,所以只有中心點是2其他都是0那乘上二項式係數之後也只有中心點有值所以沒有模糊的效果在。
所以我們調整一下算式就可以了。

//max,min
diff = diff > max/2? 2:1;

我們去找到這個矩陣中與中心點顏色差異最大的值取一半當作門檻,大於這個門檻就等於2,小於這個門檻就等於1。
結果就會是下面這樣

好像成功了,但是這裡要特別注意到邊緣的線似乎被放大了,仔細看到鬍鬚跟瞳孔邊緣,似乎有兩條線,這原因是因為我們忘記考慮到binomail filter的影響了,原本binomial filter是對每個像素都一視同仁的,只會根據與中心點的遠近給予不同的權重,但是在這之前我們透過similarity matrix相似度矩陣s給了邊緣點上的像素2的值,非邊緣點則是1,看起來沒甚麼,別忘了我們有單位化,所以他的差距不是1而是「兩倍」,這個兩倍會讓邊緣上的像素可以說是100%的去影響周圍的點,因為他的權重太高了,導致結果圖上看到的,邊緣都會變粗,這個問題其實非常好解決,應該說其實重點就在這裡,來看看wiki怎麼定義顏色差異函數的生成

這不是高斯分布的函數嗎?
只是sigma中變成了像素差異f(epsilon)-f(x)而已,原來我們定義的顏色差異是錯誤的,應該是中心點(顏色差異最小)的權重最高,接下來顏色差異越大的權重要越低,那麼不禁會疑問,這樣難道不會反而讓邊緣模糊嗎?其實這個擔心是多餘的,因為要記得一件事,這是高斯分布的「函數」,如果這個5X5範圍內每個點都與中心點顏色差異等於0,會發生甚麼事呢?沒錯s會等於{1,1,1,1,1,1,1,1,1,1,1,1,1...},那拿去給binomail filter乘不就回到原本的模糊化,反而如果周圍都是顏色差異非常之大的點,那麼會長成這樣{0,0,0,0,0,0,0........0,1,0,0,0,0.......}只有中心點會是1,這樣不就把所謂的邊緣(中心)給保留下來了嗎?所以這是數學定義上與程式碼中我們忽略的重點,我們的目標維度是哪一個才是重點,但是我們忽略了這個,就會往反方向反而弱化了邊緣的權重,所以具體而微地來看,數學公式變成程式碼其實還是有不一樣的地方,好了,那我們修正錯誤,反過來讓相似度越高的點權重越大,那麼結果果然是正確的
這樣看可能沒有甚麼感覺所以來比較一下原圖,左邊清晰的是原圖,右邊是bilateral filter

在來看看binomail filter,左邊是binomial filter,右邊是bilateral filter,可以看到左邊的鬍鬚幾乎不見了,但右邊的還有保留住,還有一些白色的毛都被保留住清晰的影像了

但是這個演算法有個嚴重的缺點,那就是非常的花費時間,基本上I7的CPU來跑一個760*600的小圖片就要花一分鐘。
我們會在下一章節,以程式設計的角度去拆解所有filter演算法,利用數學性質以及演算法的理論去探討,然後加速它們。

那麼filter的介紹就到這裡結束了。

參考資料:wiki、google圖片
參考網址:Bilateral Filtering for Gray and Color Images雙邊濾波器 (Bilateral Filter)

沒有留言:

張貼留言