基於模糊集理論的一種影象二值化演算法的原理、實現效果及程式碼

Imageshop發表於2013-09-05
演算法

　　這是篇很古老的論文中的演算法，發表與1994年，是清華大學黃良凱（Liang-kai Huang) 所寫，因此國外一些論文裡和程式碼裡稱之為Huang's fuzzy thresholding method。雖然古老也很簡單，但是其演算法的原理還是值得學習的。

     該論文的原文可從此處下載： Image thresholding by minimizing the measure of fuzziness。

     該論文結合了當時處於研究熱潮的模糊集理論，提出了一種具有較好效果的影象二值化演算法，本文主要是對其進行簡單的翻譯和註釋，並提供了測試程式碼。

     一、模糊集及其隸屬度函式

     首先，我們假定X代表一副大小為M×N的具有L個色階的灰度影象，而xmn代表影象X中點（m,n)處的畫素灰度值，定義μx(xmn)表示該點具有某種屬性的隸屬度值，也就是說我們定義了一個從影象X對映到[0,1]區間的模糊子集，用專業的模糊集表達，即有：

                                       

       其中0≤μx(xmn)≤1，m=0,1,...M-1，n=0,1,...N-1。對於二值化來說，每個畫素對於其所屬的類別（前景或背景）都應該有很相近的關係，因此，我們可以這種關係來表示μx(xmn)的值。

       定義h(g)表示影象中具有灰度級g的畫素的個數，對於一個給定的閾值t,背景和前景各自色階值的平均值μ0和μ1可用下式表示：

                                

                                         

      上述μ0和μ1，可以看成是指定閾值t所對應的前景和背景的目標值，而影象X中某一點和其所述的區域之間的關係，在直覺上應該和該點的色階值與所屬區域的目標值之間的差異相關。因此，對於點(m,n)，我們提出如下的隸屬度定義函式：

         

　　其中C是一個常數，該常數使得0.5≤μx(xmn)≤1。因此，對於一個給定的閾值t，影象中任何一個畫素要麼屬於背景，要麼屬於前景，因此，每個畫素的隸屬度不應小於0.5。

     C值在實際的程式設計中，可以用影象的最大灰度值減去最小灰度值來表達，即 C=gmax-gmin;

　　二、模糊度的度量及取閾值的原則

     模糊度表示了一個模糊集的模糊程度，有好幾種度量方式已經被提及了，本文僅僅使用了夏農熵函式來度量模糊度。

     基於夏農熵函式，一個模糊集A的熵定義為：

      

      其中夏農函式：

     

      擴充套件到2維的影象，影象X的熵可以表達為：

    

　　因為灰度影象至多隻有L個色階，因此使用直方圖式（7）可進一步寫成：

        

    可以證明式（6）在區間[0,0.5]之間是單調遞增而在[0.5,1]之間是單調遞減的，並且E(X)具有以下屬性：

     （1）0≤E(X)≤1 ;

      (2)如果μx(xmn)=0或者μx(xmn)=1時，E(X)具有最小值0，在本文中μx(xmn)只可能為1，此時分類具有最好的明確性。

     （3）當μx(xmn)=0.5，E(X)獲得最大值1，此時的分類具有最大的不明確性。

     那麼對於影象X，我們確定最好的閾值t的原則就是：對於所有的可能的閾值t,取夏農熵值最小時的那個t為最終的分割閾值。

     三、程式設計中的技巧

     有了上述原理，其實程式設計也是件很容易的事情了，你可以按照你的想法去做，不過作者論文中的闡述會讓程式碼寫起來更清晰、更有效。

     首先，為了表達方便，我們定義如下一些表示式：

       

     根據上述表示式，可以知道S(L-1)及W(L-1)對於一副影象來說是個常量，其中S(L-1)明顯就是畫素的總個數。

　　我們的演算法步驟如下：

     （1）、計算S(L-1)、W(L-1)，設定初始閾值t=gmin,令S(t-1)=0、W(t-1)=0;

      (2)、 計算下面算式:

               

     　稍微有點數學基礎的人都應該能看懂上述算式的推導原理。

        根據式（2）和式（3），可以知道背景和前景的區域的平均灰度值為：

      　　　

　　 上式中int表示取整操作。

       （3）根據式（4）及式（11）計算影象的模糊度;

       （4）令t=t+1，然後重新執行步驟2，直到t=gmax-1;

        (5)找到整個過程中的最小模糊度值對應的閾值t，並作為最佳的分割閾值。

　　   為了稍微加快點速度，上述式4中的計算可以在步驟1中用一查詢表實現。

    四、參考程式碼：

public static int GetHuangFuzzyThreshold(int[] HistGram)
{
    int X, Y;
    int First, Last;
    int Threshold = -1;
    double BestEntropy = Double.MaxValue, Entropy;
    //   找到第一個和最後一個非0的色階值
    for (First = 0; First < HistGram.Length && HistGram[First] == 0; First++) ;
    for (Last = HistGram.Length - 1; Last > First && HistGram[Last] == 0; Last--) ;
    if (First == Last) return First;                // 影象中只有一個顏色
    if (First + 1 == Last) return First;            // 影象中只有二個顏色

    // 計算累計直方圖以及對應的帶權重的累計直方圖
    int[] S = new int[Last + 1];
    int[] W = new int[Last + 1];            // 對於特大圖，此陣列的儲存資料可能會超出int的表示範圍，可以考慮用long型別來代替
    S[0] = HistGram[0];
    for (Y = First > 1 ? First : 1; Y <= Last; Y++)
    {
        S[Y] = S[Y - 1] + HistGram[Y];
        W[Y] = W[Y - 1] + Y * HistGram[Y];
    }

    // 建立公式（4）及（6）所用的查詢表
    double[] Smu = new double[Last + 1 - First];
    for (Y = 1; Y < Smu.Length; Y++)
    {
        double mu = 1 / (1 + (double)Y / (Last - First));               // 公式（4）
        Smu[Y] = -mu * Math.Log(mu) - (1 - mu) * Math.Log(1 - mu);      // 公式（6）
    }

    // 迭代計算最佳閾值
    for (Y = First; Y <= Last; Y++)
    {
        Entropy = 0;
        int mu = (int)Math.Round((double)W[Y] / S[Y]);             // 公式17
        for (X = First; X <= Y; X++)
            Entropy += Smu[Math.Abs(X - mu)] * HistGram[X];
        mu = (int)Math.Round((double)(W[Last] - W[Y]) / (S[Last] - S[Y]));  // 公式18
        for (X = Y + 1; X <= Last; X++)
            Entropy += Smu[Math.Abs(X - mu)] * HistGram[X];       // 公式8
        if (BestEntropy > Entropy)
        {
            BestEntropy = Entropy;      // 取最小熵處為最佳閾值
            Threshold = Y;
        }
    }
    return Threshold;
}

   程式碼其實還是很簡單的。

   五、效果：

　針對一些影象，我們做了如下測試：

  　

　　　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　二值圖，閾值=175

   上圖使用OSTU等經典演算法都無法獲得上圖的理想效果。

       

   　　　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　   二值圖，閾值=67

    上圖和其他一些二值演算法的效果也是非常類似的。

 

 

*********************************作者： laviewpbt   時間： 2013.9.5    聯絡QQ:  33184777  轉載請保留本行資訊************************

相關文章