十三種基於直方圖的影象全域性二值化演算法原理、實現、程式碼及效果。

Imageshop發表於2013-09-07

影象二值化的目的是最大限度的將圖象中感興趣的部分保留下來，在很多情況下，也是進行影象分析、特徵提取與模式識別之前的必要的影象預處理過程。這個看似簡單的問題，在過去的四十年裡受到國內外學者的廣泛關注，產生了數以百計的閾值選取方法，但如同其他影象分割演算法一樣，沒有一個現有方法對各種各樣的影象都能得到令人滿意的結果。

在這些龐大的分類方法中，基於直方圖的全域性二值演算法佔有了絕對的市場份額，這些演算法都從不同的科學層次提出了各自的實施方案，並且這類方法都有著一些共同的特點：

　　1、簡單；

2、演算法容易實現；

3、執行速度快。

本文摘取了若干種這類方法進行了介紹。

一：灰度平局值值法：

　　1、描述:即使用整幅影象的灰度平均值作為二值化的閾值，一般該方法可作為其他方法的初始猜想值。

　　2、原理:　　　

3、實現程式碼：

 public static int GetMeanThreshold(int[] HistGram)
    {
        int Sum = 0, Amount = 0;
        for (int Y = 0; Y < 256; Y++)
        {
            Amount += HistGram[Y];
            Sum += Y * HistGram[Y];
        }
        return Sum / Amount;
    }

二、百分比閾值（P-Tile法）

1、描述

Doyle於1962年提出的P-Tile (即P分位數法)可以說是最古老的一種閾值選取方法。該方法根據先驗概率來設定閾值，使得二值化後的目標或背景畫素比例等於先驗概率，該方法簡單高效，但是對於先驗概率難於估計的影象卻無能為力。

　　2、該原理比較簡單，直接以程式碼實現。

    /// <summary>
    /// 百分比閾值
    /// </summary>
    /// <param name="HistGram">灰度影象的直方圖</param>
    /// <param name="Tile">背景在影象中所佔的面積百分比</param>
    /// <returns></returns>
    public static int GetPTileThreshold(int[] HistGram, int Tile = 50)
    {
        int Y, Amount = 0, Sum = 0;
        for (Y = 0; Y < 256; Y++) Amount += HistGram[Y];        //  畫素總數
         for (Y = 0; Y < 256; Y++)
        {
            Sum = Sum + HistGram[Y];
            if (Sum >= Amount * Tile / 100) return Y;
        }
        return -1;
    }

三、基於谷底最小值的閾值

1、描述：

　　此方法實用於具有明顯雙峰直方圖的影象，其尋找雙峰的谷底作為閾值，但是該方法不一定能獲得閾值，對於那些具有平坦的直方圖或單峰影象，該方法不合適。

　　2、實現過程：

　　該函式的實現是一個迭代的過程，每次處理前對直方圖資料進行判斷，看其是否已經是一個雙峰的直方圖，如果不是，則對直方圖資料進行半徑為1（視窗大小為3）的平滑，如果迭代了一定的數量比如1000次後仍未獲得一個雙峰的直方圖，則函式執行失敗，如成功獲得，則最終閾值取兩個雙峰之間的谷底值作為閾值。

注意在編碼過程中，平滑的處理需要當前畫素之前的資訊，因此需要對平滑前的資料進行一個備份。另外，首資料型別精度限制，不應用整形的直方圖資料，必須轉換為浮點型別資料來進行處理，否則得不到正確的結果。

該演算法相關參考論文如下：

　 J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," innnals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
C. A. Glasbey, "An analysis of histogram-based thresholding algorithms," CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.

3、實現程式碼：

    public static int GetMinimumThreshold(int[] HistGram)
    {
        int Y, Iter = 0;
        double[] HistGramC = new double[256];           // 基於精度問題，一定要用浮點數來處理，否則得不到正確的結果
        double[] HistGramCC = new double[256];          // 求均值的過程會破壞前面的資料，因此需要兩份資料
        for (Y = 0; Y < 256; Y++)
        {
            HistGramC[Y] = HistGram[Y];
            HistGramCC[Y] = HistGram[Y];
        }

        // 通過三點求均值來平滑直方圖
        while (IsDimodal(HistGramCC) == false)                                        // 判斷是否已經是雙峰的影象了      
        {
            HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3;                 // 第一點
            for (Y = 1; Y < 255; Y++)
                HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3;     // 中間的點
            HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3;         // 最後一點
            System.Buffer.BlockCopy(HistGramCC, 0, HistGramC, 0, 256 * sizeof(double));
            Iter++;
            if (Iter >= 1000) return -1;                                                   // 直方圖無法平滑為雙峰的，返回錯誤程式碼
        }
       // 閾值極為兩峰之間的最小值 
        bool Peakfound = false;
        for (Y = 1; Y < 255; Y++)
        {
            if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peakfound = true;
            if (Peakfound == true && HistGramCC[Y - 1] >= HistGramCC[Y] && HistGramCC[Y + 1] >= HistGramCC[Y])
                return Y - 1;
        }
        return -1;
    }

　　其中IsDimodal函式為判斷直方圖是否是雙峰的函式，程式碼如下：

    private static bool IsDimodal(double[] HistGram)       // 檢測直方圖是否為雙峰的
    {
        // 對直方圖的峰進行計數，只有峰數位2才為雙峰 
        int Count = 0;
        for (int Y = 1; Y < 255; Y++)
        {
            if (HistGram[Y - 1] < HistGram[Y] && HistGram[Y + 1] < HistGram[Y])
            {
                Count++;
                if (Count > 2) return false;
            }
        }
        if (Count == 2)
            return true;
        else
            return false;
    }

　　4、效果：

　　　　　　　　原圖二值圖原始直方圖　　　　　　　　　　　　　　　　平滑後的直方圖

　　對於這種有較明顯的雙峰的影象，該演算法還是能取得不錯的效果的。

四、基於雙峰平均值的閾值

1、描述：

　　該演算法和基於谷底最小值的閾值方法類似，只是最後一步不是取得雙峰之間的谷底值，而是取雙峰的平均值作為閾值。

2、參考程式碼：

 public static int GetIntermodesThreshold(int[] HistGram)
    {
        int Y, Iter = 0, Index;
        double[] HistGramC = new double[256];           // 基於精度問題，一定要用浮點數來處理，否則得不到正確的結果
        double[] HistGramCC = new double[256];          // 求均值的過程會破壞前面的資料，因此需要兩份資料
        for (Y = 0; Y < 256; Y++)
        {
            HistGramC[Y] = HistGram[Y];
            HistGramCC[Y] = HistGram[Y];
        }
        // 通過三點求均值來平滑直方圖
        while (IsDimodal(HistGramCC) == false)                                                  // 判斷是否已經是雙峰的影象了      
        {
            HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3;                   // 第一點
            for (Y = 1; Y < 255; Y++)
                HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3;       // 中間的點
            HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3;           // 最後一點
            System.Buffer.BlockCopy(HistGramCC, 0, HistGramC, 0, 256 * sizeof(double));         // 備份資料，為下一次迭代做準備
            Iter++;
            if (Iter >= 10000) return -1;                                                       // 似乎直方圖無法平滑為雙峰的，返回錯誤程式碼
        }
// 閾值為兩峰值的平均值
        int[] Peak = new int[2];
        for (Y = 1, Index = 0; Y < 255; Y++)
            if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peak[Index++] = Y - 1;
        return ((Peak[0] + Peak[1]) / 2);
    }

3、效果：

　　　　　　　原圖二值圖原始直方圖　　　　　　　　　　　　　　　　平滑後的直方圖

五、迭代最佳閾值

　　1、描述：

該演算法先假定一個閾值，然後計算在該閾值下的前景和背景的中心值，當前景和背景中心值得平均值和假定的閾值相同時，則迭代中止，並以此值為閾值進行二值化。

2、實現過程：

　　（1）求出圖象的最大灰度值和最小灰度值，分別記為g_l和g_u，令初始閾值為：

(2) 根據閾值T⁰將圖象分割為前景和背景，分別求出兩者的平均灰度值A_b和A_f:

(3) 令

　　　　如果T^k=T^k+1,則取T^k為所求得的閾值，否則，轉2繼續迭代。

3、參考程式碼：

    public static int GetIterativeBestThreshold(int[] HistGram)
    {
        int X, Iter = 0;
        int MeanValueOne, MeanValueTwo, SumOne, SumTwo, SumIntegralOne, SumIntegralTwo;
        int MinValue, MaxValue;
        int Threshold, NewThreshold;

        for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
        for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;

        if (MaxValue == MinValue) return MaxValue;          // 影象中只有一個顏色             
        if (MinValue + 1 == MaxValue) return MinValue;      // 影象中只有二個顏色

        Threshold = MinValue;
        NewThreshold = (MaxValue + MinValue) >> 1;
        while (Threshold != NewThreshold)    // 當前後兩次迭代的獲得閾值相同時，結束迭代    
        {
            SumOne = 0; SumIntegralOne = 0;
            SumTwo = 0; SumIntegralTwo = 0;
            Threshold = NewThreshold;
            for (X = MinValue; X <= Threshold; X++)         //根據閾值將影象分割成目標和背景兩部分，求出兩部分的平均灰度值      
            {
                SumIntegralOne += HistGram[X] * X;
                SumOne += HistGram[X];
            }
            MeanValueOne = SumIntegralOne / SumOne;
            for (X = Threshold + 1; X <= MaxValue; X++)
            {
                SumIntegralTwo += HistGram[X] * X;
                SumTwo += HistGram[X];
            }
            MeanValueTwo = SumIntegralTwo / SumTwo;
            NewThreshold = (MeanValueOne + MeanValueTwo) >> 1;       //求出新的閾值
            Iter++;
            if (Iter >= 1000) return -1;
        }
        return Threshold;
    }

　　4、效果：

　　　　原圖二值圖直方圖　

六、OSTU大律法

　　1、描述：

該演算法是1979年由日本大津提出的，主要是思想是取某個閾值，使得前景和背景兩類的類間方差最大，matlab中的graythresh即是以該演算法為原理執行的。

2、原理：

關於該演算法的原理，網路上有很多，這裡為了篇幅有限，不加以贅述。

3、參考程式碼：

    public static int GetOSTUThreshold(int[] HistGram)
    {
        int X, Y, Amount = 0;
        int PixelBack = 0, PixelFore = 0, PixelIntegralBack = 0, PixelIntegralFore = 0, PixelIntegral = 0;
        double OmegaBack, OmegaFore, MicroBack, MicroFore, SigmaB, Sigma;              // 類間方差;
        int MinValue, MaxValue;
        int Threshold = 0;

        for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
        for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
        if (MaxValue == MinValue) return MaxValue;          // 影象中只有一個顏色             
        if (MinValue + 1 == MaxValue) return MinValue;      // 影象中只有二個顏色

        for (Y = MinValue; Y <= MaxValue; Y++) Amount += HistGram[Y];        //  畫素總數

        PixelIntegral = 0;
        for (Y = MinValue; Y <= MaxValue; Y++) PixelIntegral += HistGram[Y] * Y;
        SigmaB = -1;
        for (Y = MinValue; Y < MaxValue; Y++)
        {
            PixelBack = PixelBack + HistGram[Y];
            PixelFore = Amount - PixelBack;
            OmegaBack = (double)PixelBack / Amount;
            OmegaFore = (double)PixelFore / Amount;
            PixelIntegralBack += HistGram[Y] * Y;
            PixelIntegralFore = PixelIntegral - PixelIntegralBack;
            MicroBack = (double)PixelIntegralBack / PixelBack;
            MicroFore = (double)PixelIntegralFore / PixelFore;
            Sigma = OmegaBack * OmegaFore * (MicroBack - MicroFore) * (MicroBack - MicroFore);
            if (Sigma > SigmaB)
            {
                SigmaB = Sigma;
                Threshold = Y;
            }
        }
        return Threshold;
    }

　　4、效果：

　　該演算法對於那些具有平坦的直方圖的影象具有一定的適應能力。

七、一維最大熵

1、描述：

該演算法把資訊理論中熵的概念引入到影象中，通過計算閾值分割後兩部分熵的和來判斷閾值是否為最佳閾值。

2、演算法原理

這方面的文章也比較多，留給讀者自行去查詢相關資料。

3、參考程式碼：

public static int Get1DMaxEntropyThreshold(int[] HistGram)
    {
        int  X, Y,Amount=0;
        double[] HistGramD = new double[256];
        double SumIntegral, EntropyBack, EntropyFore, MaxEntropy;
        int MinValue = 255, MaxValue = 0;
        int Threshold = 0;

        for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
        for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
        if (MaxValue == MinValue) return MaxValue;          // 影象中只有一個顏色             
        if (MinValue + 1 == MaxValue) return MinValue;      // 影象中只有二個顏色

        for (Y = MinValue; Y <= MaxValue; Y++) Amount += HistGram[Y];        //  畫素總數

        for (Y = MinValue; Y <= MaxValue; Y++)   HistGramD[Y] = (double)HistGram[Y] / Amount+1e-17;

        MaxEntropy = double.MinValue; ;
        for (Y = MinValue + 1; Y < MaxValue; Y++)
        {
            SumIntegral = 0;
            for (X = MinValue; X <= Y; X++) SumIntegral += HistGramD[X];
            EntropyBack = 0;
            for (X = MinValue; X <= Y; X++) EntropyBack += (-HistGramD[X] / SumIntegral * Math.Log(HistGramD[X] / SumIntegral));
            EntropyFore = 0;
            for (X = Y + 1; X <= MaxValue; X++) EntropyFore += (-HistGramD[X] / (1 - SumIntegral) * Math.Log(HistGramD[X] / (1 - SumIntegral)));
            if (MaxEntropy < EntropyBack + EntropyFore)
            {
                Threshold = Y;
                MaxEntropy = EntropyBack + EntropyFore;
            }
        }
        return Threshold;
    }

八、力矩保持法
1、描述:

　　該演算法通過選擇恰當的閾值從而使得二值後的影象和原始的灰度影象具有三個相同的初始力矩值。

　　2、原理：

　參考論文：W. Tsai, “Moment-preserving thresholding: a new approach,” Comput.Vision Graphics Image Process., vol. 29, pp. 377-393, 1985.

由於無法下載到該論文（收費的），僅僅給出從其他一些資料中找到的公式共享一下。

其中的A\B\C的函式可見程式碼部分。

3、參考程式碼：

public static byte GetMomentPreservingThreshold(int[] HistGram)
    {
        int X, Y, Index = 0, Amount=0;
        double[] Avec = new double[256];
        double X2, X1, X0, Min;

        for (Y = 0; Y <= 255; Y++) Amount += HistGram[Y];        //  畫素總數
        for (Y = 0; Y < 256; Y++) Avec[Y] = (double)A(HistGram, Y) / Amount;       // The threshold is chosen such that A(y,t)/A(y,n) is closest to x0.

        // The following finds x0.

        X2 = (double)(B(HistGram, 255) * C(HistGram, 255) - A(HistGram, 255) * D(HistGram, 255)) / (double)(A(HistGram, 255) * C(HistGram, 255) - B(HistGram, 255) * B(HistGram, 255));
        X1 = (double)(B(HistGram, 255) * D(HistGram, 255) - C(HistGram, 255) * C(HistGram, 255)) / (double)(A(HistGram, 255) * C(HistGram, 255) - B(HistGram, 255) * B(HistGram, 255));
        X0 = 0.5 - (B(HistGram, 255) / A(HistGram, 255) + X2 / 2) / Math.Sqrt(X2 * X2 - 4 * X1);

        for (Y = 0, Min = double.MaxValue; Y < 256; Y++)
        {
            if (Math.Abs(Avec[Y] - X0) < Min)
            {
                Min = Math.Abs(Avec[Y] - X0);
                Index = Y;
            }
        }
        return (byte)Index;
    }

    private static double A(int[] HistGram, int Index)
    {
        double Sum = 0;
        for (int Y = 0; Y <= Index; Y++)
            Sum += HistGram[Y];
        return Sum;
    }

    private static double B(int[] HistGram, int Index)
    {
        double Sum = 0;
        for (int Y = 0; Y <= Index; Y++)
            Sum += (double)Y * HistGram[Y];
        return Sum;
    }

    private static double C(int[] HistGram, int Index)
    {
        double Sum = 0;
        for (int Y = 0; Y <= Index; Y++)
            Sum += (double)Y * Y * HistGram[Y];
        return Sum;
    }

    private static double D(int[] HistGram, int Index)
    {
        double Sum = 0;
        for (int Y = 0; Y <= Index; Y++)
            Sum += (double)Y * Y * Y * HistGram[Y];
        return Sum;
    }

對於很多影象，該演算法也能取得比較滿意的結果。

九、基於模糊集理論的閾值

該演算法的具體分析可見：基於模糊集理論的一種影象二值化演算法的原理、實現效果及程式碼

此法也借用夏農熵的概念，該演算法一般都能獲得較為理想的分割效果，不管是對雙峰的還是單峰的影象。

十、Kittler最小錯誤分類法

由於精力有限，以下幾種演算法僅僅給出演算法的論文及相關的程式碼。

該演算法具體的分析見：

Kittler, J & Illingworth, J (1986), "Minimum error thresholding", Pattern Recognition 19: 41-47

　　參考程式碼：

 public static int GetKittlerMinError(int[] HistGram)
    {
        int X, Y;
        int MinValue, MaxValue;
        int Threshold ;
        int PixelBack, PixelFore;
        double OmegaBack, OmegaFore, MinSigma, Sigma, SigmaBack, SigmaFore;
        for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
        for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
        if (MaxValue == MinValue) return MaxValue;          // 影象中只有一個顏色             
        if (MinValue + 1 == MaxValue) return MinValue;      // 影象中只有二個顏色
        Threshold = -1;
        MinSigma = 1E+20;
        for (Y = MinValue; Y < MaxValue; Y++)
        {
            PixelBack = 0; PixelFore = 0;
            OmegaBack = 0; OmegaFore = 0;
            for (X = MinValue; X <= Y; X++)
            {
                PixelBack += HistGram[X];
                OmegaBack = OmegaBack + X * HistGram[X];
            }
            for (X = Y + 1; X <= MaxValue; X++)
            {
                PixelFore += HistGram[X];
                OmegaFore = OmegaFore + X * HistGram[X];
            }
            OmegaBack = OmegaBack / PixelBack;
            OmegaFore = OmegaFore / PixelFore;
            SigmaBack = 0; SigmaFore = 0;
            for (X = MinValue; X <= Y; X++) SigmaBack = SigmaBack + (X - OmegaBack) * (X - OmegaBack) * HistGram[X];
            for (X = Y + 1; X <= MaxValue; X++) SigmaFore = SigmaFore + (X - OmegaFore) * (X - OmegaFore) * HistGram[X];
            if (SigmaBack == 0 || SigmaFore == 0)
            {
                if (Threshold == -1)
                    Threshold = Y;
            }
            else
            {
                SigmaBack = Math.Sqrt(SigmaBack / PixelBack);
                SigmaFore = Math.Sqrt(SigmaFore / PixelFore);
                Sigma = 1 + 2 * (PixelBack * Math.Log(SigmaBack / PixelBack) + PixelFore * Math.Log(SigmaFore / PixelFore));
                if (Sigma < MinSigma)
                {
                    MinSigma = Sigma;
                    Threshold = Y;
                }
            }
        }
        return Threshold;
    }

　　從實際的執行效果看，該演算法並不很好。

　　十一：ISODATA(也叫做intermeans法）

　　參考論文：

Ridler, TW & Calvard, S (1978), "Picture thresholding using an iterative selection method", IEEE Transactions on Systems, Man and Cybernetics 8: 630-632, <http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4310039>

參考程式碼（未做整理)：

 // Also called intermeans
    // Iterative procedure based on the isodata algorithm [T.W. Ridler, S. Calvard, Picture 
    // thresholding using an iterative selection method, IEEE Trans. System, Man and 
    // Cybernetics, SMC-8 (1978) 630-632.] 
    // The procedure divides the image into objects and background by taking an initial threshold,
    // then the averages of the pixels at or below the threshold and pixels above are computed. 
    // The averages of those two values are computed, the threshold is incremented and the 
    // process is repeated until the threshold is larger than the composite average. That is,
    //  threshold = (average background + average objects)/2
    // The code in ImageJ that implements this function is the getAutoThreshold() method in the ImageProcessor class. 
    //
    // From: Tim Morris (dtm@ap.co.umist.ac.uk)
    // Subject: Re: Thresholding method?
    // posted to sci.image.processing on 1996/06/24
    // The algorithm implemented in NIH Image sets the threshold as that grey
    // value, G, for which the average of the averages of the grey values
    // below and above G is equal to G. It does this by initialising G to the
    // lowest sensible value and iterating:

    // L = the average grey value of pixels with intensities < G
    // H = the average grey value of pixels with intensities > G
    // is G = (L + H)/2?
    // yes => exit
    // no => increment G and repeat
    //
    // There is a discrepancy with IJ because they are slightly different methods

    public static int GetIsoDataThreshold(int[] HistGram)
    {
        int i, l, toth, totl, h, g = 0;
        for (i = 1; i < HistGram.Length; i++)
        {
            if (HistGram[i] > 0)
            {
                g = i + 1;
                break;
            }
        }
        while (true)
        {
            l = 0;
            totl = 0;
            for (i = 0; i < g; i++)
            {
                totl = totl + HistGram[i];
                l = l + (HistGram[i] * i);
            }
            h = 0;
            toth = 0;
            for (i = g + 1; i < HistGram.Length; i++)
            {
                toth += HistGram[i];
                h += (HistGram[i] * i);
            }
            if (totl > 0 && toth > 0)
            {
                l /= totl;
                h /= toth;
                if (g == (int)Math.Round((l + h) / 2.0))
                    break;
            }
            g++;
            if (g > HistGram.Length - 2)
            {
                return 0;
            }
        }
        return g;
    }

十二、Shanbhag 法

參考論文：

Shanbhag, Abhijit G. (1994), "Utilization of information measure as a means of image thresholding", Graph. Models Image Process. (Academic Press, Inc.) 56 (5): 414--419, ISSN 1049-9652, DOI 10.1006/cgip.1994.1037

參考程式碼（未整理）：

public static int GetShanbhagThreshold(int[] HistGram)
    {
        int threshold;
        int ih, it;
        int first_bin;
        int last_bin;
        double term;
        double tot_ent;  /* total entropy */
        double min_ent;  /* max entropy */
        double ent_back; /* entropy of the background pixels at a given threshold */
        double ent_obj;  /* entropy of the object pixels at a given threshold */
        double[] norm_histo = new double[HistGram.Length]; /* normalized histogram */
        double[] P1 = new double[HistGram.Length]; /* cumulative normalized histogram */
        double[] P2 = new double[HistGram.Length];

        int total = 0;
        for (ih = 0; ih < HistGram.Length; ih++)
            total += HistGram[ih];

        for (ih = 0; ih < HistGram.Length; ih++)
            norm_histo[ih] = (double)HistGram[ih] / total;

        P1[0] = norm_histo[0];
        P2[0] = 1.0 - P1[0];
        for (ih = 1; ih < HistGram.Length; ih++)
        {
            P1[ih] = P1[ih - 1] + norm_histo[ih];
            P2[ih] = 1.0 - P1[ih];
        }

        /* Determine the first non-zero bin */
        first_bin = 0;
        for (ih = 0; ih < HistGram.Length; ih++)
        {
            if (!(Math.Abs(P1[ih]) < 2.220446049250313E-16))
            {
                first_bin = ih;
                break;
            }
        }

        /* Determine the last non-zero bin */
        last_bin = HistGram.Length - 1;
        for (ih = HistGram.Length - 1; ih >= first_bin; ih--)
        {
            if (!(Math.Abs(P2[ih]) < 2.220446049250313E-16))
            {
                last_bin = ih;
                break;
            }
        }

        // Calculate the total entropy each gray-level
        // and find the threshold that maximizes it 
        threshold = -1;
        min_ent = Double.MaxValue;

        for (it = first_bin; it <= last_bin; it++)
        {
            /* Entropy of the background pixels */
            ent_back = 0.0;
            term = 0.5 / P1[it];
            for (ih = 1; ih <= it; ih++)
            { //0+1?
                ent_back -= norm_histo[ih] * Math.Log(1.0 - term * P1[ih - 1]);
            }
            ent_back *= term;

            /* Entropy of the object pixels */
            ent_obj = 0.0;
            term = 0.5 / P2[it];
            for (ih = it + 1; ih < HistGram.Length; ih++)
            {
                ent_obj -= norm_histo[ih] * Math.Log(1.0 - term * P2[ih]);
            }
            ent_obj *= term;

            /* Total entropy */
            tot_ent = Math.Abs(ent_back - ent_obj);

            if (tot_ent < min_ent)
            {
                min_ent = tot_ent;
                threshold = it;
            }
        }
        return threshold;
    }

十三、Yen法

　　參考論文：

1) Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion for Automatic Multilevel Thresholding" IEEE Trans. on Image Processing, 4(3): 370-378
2) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding Techniques and Quantitative Performance Evaluation" Journal of Electronic Imaging, 13(1): 146-165

參考程式碼（未整理）：

 // M. Emre Celebi
    // 06.15.2007
    // Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
    public static int GetYenThreshold(int[] HistGram)
    {
        int threshold;
        int ih, it;
        double crit;
        double max_crit;
        double[] norm_histo = new double[HistGram.Length]; /* normalized histogram */
        double[] P1 = new double[HistGram.Length]; /* cumulative normalized histogram */
        double[] P1_sq = new double[HistGram.Length];
        double[] P2_sq = new double[HistGram.Length];

        int total = 0;
        for (ih = 0; ih < HistGram.Length; ih++)
            total += HistGram[ih];

        for (ih = 0; ih < HistGram.Length; ih++)
            norm_histo[ih] = (double)HistGram[ih] / total;

        P1[0] = norm_histo[0];
        for (ih = 1; ih < HistGram.Length; ih++)
            P1[ih] = P1[ih - 1] + norm_histo[ih];

        P1_sq[0] = norm_histo[0] * norm_histo[0];
        for (ih = 1; ih < HistGram.Length; ih++)
            P1_sq[ih] = P1_sq[ih - 1] + norm_histo[ih] * norm_histo[ih];

        P2_sq[HistGram.Length - 1] = 0.0;
        for (ih = HistGram.Length - 2; ih >= 0; ih--)
            P2_sq[ih] = P2_sq[ih + 1] + norm_histo[ih + 1] * norm_histo[ih + 1];

        /* Find the threshold that maximizes the criterion */
        threshold = -1;
        max_crit = Double.MinValue;
        for (it = 0; it < HistGram.Length; it++)
        {
            crit = -1.0 * ((P1_sq[it] * P2_sq[it]) > 0.0 ? Math.Log(P1_sq[it] * P2_sq[it]) : 0.0) + 2 * ((P1[it] * (1.0 - P1[it])) > 0.0 ? Math.Log(P1[it] * (1.0 - P1[it])) : 0.0);
            if (crit > max_crit)
            {
                max_crit = crit;
                threshold = it;
            }
        }
        return threshold;
    }

一行很多程式碼是摘自開源軟體ImageJ的資料，讀者也可以參考：http://fiji.sc/wiki/index.php/Auto_Threshold 這裡獲得更多的資訊。

最後，我對這些演算法的做了簡單的UI介面，供有興趣的讀者參考。

工程程式碼下載：http://files.cnblogs.com/Imageshop/HistgramBinaryzation.rar

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

基於模糊集理論的一種影象二值化演算法的原理、實現效果及程式碼
2013-09-05
演算法
直方圖均衡化原理與實現
2023-12-01
直方圖
【影像處理】基於OpenCV實現影像直方圖的原理
2021-10-08
OpenCV直方圖
OpenCV之影象直方圖均衡化
2019-01-04
OpenCV直方圖
一種強化的基於區域性直方圖裁剪均衡化的對比度調節演算法。
2017-12-07
直方圖演算法
灰度影像直方圖均衡化公式及實現
2018-01-28
直方圖公式
灰度直方圖均衡化及其實現
2020-11-01
直方圖
基於C++和OpenCv的SIFT_影象區域性特徵檢測演算法程式碼的實現
2017-02-27
C++OpenCV特徵演算法
【數字影像處理】直方圖均衡化詳解及程式設計實現
2011-10-28
直方圖程式設計
實現基於React的全域性提示元件Toast
2019-03-04
React元件AST
收集直方圖及檢視直方圖資訊
2015-12-05
直方圖
OpenCV成長之路(4)：影象直方圖
2015-02-10
OpenCV直方圖
【16位RAW影像處理三】直方圖均衡化及區域性直方圖均衡用於16點陣圖像的細節增強。
2021-07-23
直方圖
[Python影象處理] 十一.灰度直方圖概念及OpenCV繪製直方圖
2018-11-06
Python直方圖OpenCV
直方圖、基數、選擇性、群集因子
2017-03-19
直方圖
直方圖均衡化
2023-01-08
直方圖
OpenCV成長之路(5)：影象直方圖的應用
2015-02-10
OpenCV直方圖
考拉Android全域性滑動返回及聯動效果的實現
2018-02-06
Android
cbo心得(選擇率,基數,直方圖)(二)
2008-08-10
直方圖
Python-OpenCV 處理影象（四）：影象直方圖和反向投影
2016-05-30
PythonOpenCV直方圖
DES原理及程式碼實現
2020-04-20
影象縮放的雙線性內插值演算法的原理解析
2013-03-19
演算法
iOS基於灰度的均值雜湊演算法實現影象匹配
2016-08-22
iOS演算法
常見排序演算法原理及JS程式碼實現
2020-08-07
排序演算法JS
【效能優化】直方圖
2010-09-03
優化直方圖
微信小程式實現全域性搜尋程式碼高亮
2018-03-29
微信小程式
UIButton實現各種圖文結合的效果以及原理
2018-02-24
UI
影像演算法之直方圖均衡化（灰度影像）
2020-12-30
演算法直方圖
淺析MySQL 8.0直方圖原理
2024-05-27
MySql直方圖
改進的二值影像畫素標記演算法及程式實現
2014-08-19
演算法
HC(Histogram-based Contrast) 基於直方圖對比度的顯著性
2021-10-13
HistogramAST直方圖
Android影象處理 - 高斯模糊的原理及實現
2018-02-11
Android
[Python影象處理] 七.影象閾值化處理及演算法對比
2018-10-30
Python演算法
四種比較簡單的影象顯著性區域特徵提取方法原理及實現-----> AC/HC/LC/FT。
2014-08-04
特徵
UITableView的原理——探究及重新實現程式碼
2024-08-06
UIView
影象二值化(python+opencv)
2018-12-04
PythonOpenCV
MATLAB實現頻數直方圖——hist的使用
2013-01-20
Matlab直方圖
轉貼oracle直方圖histogram(二)
2012-11-24
Oracle直方圖Histogram

十三種基於直方圖的影象全域性二值化演算法原理、實現、程式碼及效果。

相關文章