一種快速自適應的影象二值化方法介紹 (Wellner 1993)

NO IMAGE

在手機模式識別的時候, 我們首先viewfinder裡面拿到的frame通常是RGB的或者YUV的, 如果我們需要用來做模式識別的話, 通常需要首先把彩色圖首先轉化成灰度圖. 對於RGB影象而言, 網上有充足的公式, 比如Y = 0.299R 0.587G 0.114B 等等. 如果是YUV的話, 直接用Y就是灰度圖了. 順帶說一句, 這種灰度圖通常我們用.raw檔案來表示, 用photoshop或者irfanview是可以直接開啟看效果的. 比如說這裡就有一個灰度圖的例子

data matrix

這個圖就是現在很流行的所謂Data Matrix的sample, 我們用手機的照相機拿到的灰度圖. 現在我們要把它變化成為黑白圖(二值圖). 在網上廣為流傳著很多辦法. 什麼雙峰法, P引數法等等. 今天的辦法和這些都不相同. 這個方法就是被稱之為Quick Adaptive Thresholding algorithm, 提出這個觀點的人名字叫做Pierre D. Wellner.  這裡的網頁上就有這個演算法的說明:

http://www.xrce.xerox.com/Publications/Attachments/1993-110/EPC-1993-110.pdf

這個演算法的基本思想要確定一個畫素的黑或者白, 用他周圍, 或者掃描順序上的其他點的一些平均值來評估閥值就可以了.用閥值和畫素值比較即可. 我們現在定義出這樣的模型, 比方說我們用P(n)來表示第n個點的灰度值. T(n)來表示二值化後的值

用f­s (n) 來表示第n個點之前s個點的灰度值的和, 就是

公式1

 

用這個s和另一個變數t就可以簡單的說明P(n)應該是0還是1了, 這個公式就是

公式二

而且根據經驗值來看, 這裡的s和t最佳的取值範圍是s= image.width/8, 而t=15的時候效果最好.

好的, 到這裡為止, 我們的理解就是一個點t(n), 他是0還是1取決於什麼呢? 就是前面s個點的和除以s (就是前s個點的平均值)*0.85 如果這個點的灰度值<前面的值, 那麼就是1黑色, 如果大於就是0白色. 是不是非常簡單? 至少到現在為止確實是的.但是這個演算法有個問題, 我們忽略了一個問題, 就是我們現在定義T(n)的時候, 用的是平均值, 也就是說之前掃描過的若干點對於當前點的影響或者說權重是一樣的, 也就是說當前點1個畫素距離的畫素和s-1個畫素點的距離的畫素的灰度值對當前點的影響是一樣的. 顯然根據我們直觀的理解來看,
應該是離當前點越近的畫素對當前點的影響越大, 越遠則越小. 所以演算法的作者發明了這個個更合適更高效的替代值gs (n). 這個值的意義就是:

公式三

可以看到, 這裡的gs (n) 和f­s (n) 的區別在於f­s (n) 直接是不做任何修正的s個灰度值的和, 而gs (n)則是一定比例的灰度值的和, 可以看到, 離這個n越近的畫素的比重越高, 越遠越低. 顯然這樣描述對把握畫素的顏色更為準確.
而且這裡的 gs (n)和 gs (n-1)通過加法和乘法就可以遞迴得到, 計算效率是比較高的.

 

即使到了這一步了, 還有一個問題存在, 就是我現在的顏色計算其實依賴於我的掃描順序, 也就是說P(n)的這個序列的定義就是我的掃描順序(一般都是水平掃描的). 這樣的話, 我的畫素值實際上取決於我水平位置上的鄰接點的灰度值, 可是豎直方向的畫素如何關聯起來呢? 這裡也有一個說明, 我們可以維護前面依次水平掃描產生的g_prev(n)序列, 在某個g(n)被使用之前, 我們可以讓他和前一個g_prev(n)取一個平均值, 這樣的話, 這個最終的值就更有說服力了.

公式四

好了, 到現在為止, 我們描述了整個演算法的全過程, 在加上我們定義的初始g(n)值127*s(127表示0-255之間的中間值)就可以開始實現演算法了

  1. void quickAdaptiveThreshold(unsigned char* grayscale, unsigned char*& thres, int width, int height )  
  2. {  
  3.       
  4.     /**           / 
  5.     *            | FOREGROUND, if pn < ((gs(n)   gs(n-w)) / (2*s)) * 
  6.     * color(n) = |                     ((100-t)/100) 
  7.     *            | BACKGROUND_QR, otherwise 
  8.     *            / 
  9.     * where pn = gray value of current pixel, 
  10.     *        s = width of moving average, and 
  11.     *        t = threshold percentage of brightness range 
  12.     *    gs(n) = gs(n-1) * (1-1/s)   pn 
  13.     *    gs(n-w) = gs-value of pixel above current pixel 
  14.     * 
  15.     */  
  16.     int t = 15;   
  17.     int s = width >> 3; // s: number of pixels in the moving average (w = image width)  
  18.     const int S = 9; // integer shift, needed to avoid floating point operations  
  19.     const int power2S = 1 << S;  
  20.     // for speedup: multiply all values by 2^s, and use integers instead of floats  
  21.     int factor = power2S * (100-t) / (100*s); // multiplicand for threshold  
  22.     int gn = 127 * s; // initial value of the moving average (127 = average gray value)  
  23.     int q = power2S – power2S / s; // constant needed for average computation  
  24.     int pn, hn;  
  25.     unsigned char *scanline = NULL;  
  26.     int *prev_gn = NULL;  
  27.           
  28.     prev_gn = new int[width];  
  29.     for (int i = 0; i < width; i ) {  
  30.         prev_gn[i] = gn;  
  31.     }  
  32.     thres = new unsigned char[width*height];  
  33.     for (int y = 0; y < height; y   )  
  34.     {  
  35.         int yh = y * width;  
  36.         scanline = grayscale   y * width;  
  37.         for ( int x = 0; x <width; x   )  
  38.         {  
  39.             pn = scanline[x] ;  
  40.             gn = ((gn * q) >> S)   pn;   
  41.             hn = (gn   prev_gn[x]) >> 1;  
  42.             prev_gn[x] = gn;          
  43.             pn < (hn*factor) >> S ? thres[yh x] = 0 : thres[yh x] = 0xff;  
  44.         }  
  45.         y   ;  
  46.         if ( y == height)  
  47.             break;  
  48.         yh = y * width;  
  49.         scanline = grayscale   y * width;  
  50.         for ( int x = width-1; x >= 0; x –)  
  51.         {  
  52.             pn = scanline[x] ;  
  53.             gn = ((gn * q) >> S)   pn;   
  54.             hn = (gn   prev_gn[x]) >> 1;  
  55.             prev_gn[x] = gn;          
  56.             pn < (hn*factor) >> S ? thres[yh x] = 0 : thres[yh x] = 0xff;  
  57.         }  
  58.     }  
  59.     delete prev_gn;  
  60. }  

這個演算法也不是我發明創造的, 這個演算法從http://mikie.iki.fi/lxr/source/VisualCodeSystem/src/RecognitionAlgorithm.cpp?v=v3 這個網址上看過來.
不過是去除了一些Symbian的痕跡, 還有有的細節上做了一些改進, 讓程式碼更加合理了些. 經過這個演算法, 我們可以來看看效果了

原圖1:                                                                       二值圖1:

原圖1 效果圖1

 

原圖2:                                                                           二值圖2:

原圖2 二值圖2

我對這個效果還是比較滿意的. 哈哈. 一直要寫這個文章, 今天終於寫完了, 心裡真是痛快了.

 

Wellner
1993快速自適應的影象二值化方法的提高 (Derek Bradley and Gerhard Roth 2007)

前面一種方案實際上還是存在一定的問題的, 就是這個避重就輕的初始g(n)值127*s(127表示0-255之間的中間值), 這個東西帶來的最直接的問題就是邊緣的效果在這個演算法下是不咋地的。 其實從這個所謂的”Wellner 1993″, 後人又做了很多的改進, 使之效率更高, 效果更好。比方說這個Derek Bradley和Gerhard Roth搞的這個所謂 Adaptive Thresholding Using the Integral Image 在這個網頁

http://www.scs.carleton.ca/~roth/iit-publications-iti/docs/gerh-50002.pdf 可以看到一些他的蹤跡。

 

這個演算法的基本思想是這樣的,為了打破原來演算法的初始值問題以及掃描順序的問題, 這裡的畫素二值化的時候, 直接使用周圍矩形畫素的顏色作比較,這樣來判斷畫素值更科學。我們對演算法的介紹從求和麵積表(Summed-Area Table)開始. 這個求和麵積表簡單點說就是維護一張表, 表中的元素值就是它左上位置的所有畫素的畫素值和。(數學公式在這裡編輯簡直是噩夢!只能放圖了無圖無真相:))

示意圖

左邊就是原始畫素值, 右邊的就是累加得到的表, 比方說這個表裡面的(2,2)位置的8就是通過2 3 3 0得到的, 而這個最大值28就是所有畫素的累加和。得到這個和和我們的二值化有什麼關聯呢?前面我們提到了在新的這個演算法裡面畫素的值以來於周圍畫素的顏色, 那周圍畫素的顏色如何表示呢? 我們可以通過這個表輕鬆獲得, 且看下面一張圖:

示意圖2

這裡的UL, LL, UR, LR表示的就是前面這個求和表裡面的值, 如果我們要判斷綠色區域中這個 號位置的值, 我們就要計算整個綠色區域的平均畫素值, 如何計算呢? 有了新的表就方便了,右邊其實給出了這個公式,這裡的LR-UR-LL UL就是整個綠色區域的畫素值和。這個什麼道理其實已經自己可以推斷出來了, 如果還嫌這裡不清楚的話,我們就給個更清楚的圖:

示意圖2

這個圖和前面一樣,但是如果還是用LR-UR-LL UL來表示的話,這裡就可以寫成:

LR-UR-LL UL = (A B C D)-(A B)-(A C) A = D, 這樣就清楚很多了吧。 得到的這個值D就是D這個區域的畫素值和, 那D中最中心的畫素的顏色就可以用D/(widith*height)來做比較了。 所以演算法的流程就是首先得到這個求和麵積表, 其次遍歷所有的畫素, 然後以這些畫素為中心點, 計算S*S大小的矩形的平均顏色, 用來和當前畫素比較即可。這個流程可以說是相當精煉啊!這裡依然用到了原來的S, T, 還保持了一致S是寬度的八分之一, 而T則是15,下面有一段我改過的實現程式碼:

  1. void adaptiveThreshold(unsigned char* input, unsigned char*& bin, int width, int height)  
  2. {  
  3.     int S = width >> 3;  
  4.     int T = 15;  
  5.       
  6.     unsigned long* integralImg = 0;  
  7.     int i, j;  
  8.     long sum=0;  
  9.     int count=0;  
  10.     int index;  
  11.     int x1, y1, x2, y2;  
  12.     int s2 = S/2;  
  13.       
  14.     bin = new unsigned char[width*height];  
  15.     // create the integral image  
  16.     integralImg = (unsigned long*)malloc(width*height*sizeof(unsigned long*));  
  17.     for (i=0; i<width; i )  
  18.     {  
  19.         // reset this column sum  
  20.         sum = 0;  
  21.         for (j=0; j<height; j )  
  22.         {  
  23.             index = j*width i;  
  24.             sum  = input[index];  
  25.             if (i==0)  
  26.                 integralImg[index] = sum;  
  27.             else  
  28.                 integralImg[index] = integralImg[index-1]   sum;  
  29.         }  
  30.     }  
  31.     // perform thresholding  
  32.     for (i=0; i<width; i )  
  33.     {  
  34.         for (j=0; j<height; j )  
  35.         {  
  36.             index = j*width i;  
  37.             // set the SxS region  
  38.             x1=i-s2; x2=i s2;  
  39.             y1=j-s2; y2=j s2;  
  40.             // check the border  
  41.             if (x1 < 0) x1 = 0;  
  42.             if (x2 >= width) x2 = width-1;  
  43.             if (y1 < 0) y1 = 0;  
  44.             if (y2 >= height) y2 = height-1;  
  45.             count = (x2-x1)*(y2-y1);  
  46.             // I(x,y)=s(x2,y2)-s(x1,y2)-s(x2,y1) s(x1,x1)  
  47.             sum = integralImg[y2*width x2] –  
  48.                 integralImg[y1*width x2] –  
  49.                 integralImg[y2*width x1]    
  50.                 integralImg[y1*width x1];  
  51.             if ((long)(input[index]*count) < (long)(sum*(100-T)/100))  
  52.                 bin[index] = 0;  
  53.             else  
  54.                 bin[index] = 255;  
  55.         }  
  56.     }  
  57.     free (integralImg);  
  58. }  

這裡也有一點效果圖可以看看, 同時有和前面一個演算法的比較:

 

原始1                                      wellnar演算法                            最新

原始圖1 wellnar 最新dm

 

 

還有一組:

ez_raw

wellnar:

wellnar

最新演算法:

 

new

 

 

這些個貼圖其實還不是特別的具體, 其實這個演算法特別適用於光照強度變化很大的畫素, 這裡有些網頁也給出了鮮明的對比:http://www.derekbradley.ca/AdaptiveThresholding/index.html 效果的差距還是很明顯的。
總的來說這個演算法實現簡單, 效率很高,確實是不錯的選擇。 而且還很新!在07年的雜誌上發表的,現在記錄下來與君共勉之!