http://www.cnblogs.com/carekee/articles/3643394.htmlphp
圖像二值化的目的是最大限度的將圖象中感興趣的部分保留下來,在不少狀況下,也是進行圖像分析、特徵提取與模式識別以前的必要的圖像預處理過程。這個看似簡單的問題,在過去的四十年裏受到國內外學者的普遍關注,產生了數以百計的閾值選取方法,但如同其餘圖像分割算法同樣,沒有一個現有方法對各類各樣的圖像都能獲得使人滿意的結果。html
在這些龐大的分類方法中,基於直方圖的全局二值算法佔有了絕對的市場份額,這些算法都從不一樣的科學層次提出了各自的實施方案,而且這類方法都有着一些共同的特色:算法
一、簡單;網絡
二、算法容易實現;app
三、執行速度快。jsp
本文摘取了若干種這類方法進行了介紹。ide
一:灰度平局值值法:函數
一、描述:即便用整幅圖像的灰度平均值做爲二值化的閾值,通常該方法可做爲其餘方法的初始猜測值。post
二、原理: this
三、實現代碼:
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; }
2、百分比閾值(P-Tile法)
一、描述
Doyle於1962年提出的P-Tile (即P分位數法)能夠說是最古老的一種閾值選取方法。該方法根據先驗機率來設定閾值,使得二值化後的目標或背景像素比例等於先驗機率,該方法簡單高效,可是對於先驗機率難於估計的圖像卻無能爲力。
二、該原理比較簡單,直接以代碼實現。
/// <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; }
3、基於谷底最小值的閾值
一、描述:
此方法實用於具備明顯雙峯直方圖的圖像,其尋找雙峯的谷底做爲閾值,可是該方法不必定能得到閾值,對於那些具備平坦的直方圖或單峯圖像,該方法不合適。
二、實現過程:
該函數的實現是一個迭代的過程,每次處理前對直方圖數據進行判斷,看其是否已是一個雙峯的直方圖,若是不是,則對直方圖數據進行半徑爲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.
三、實現代碼:
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、基於雙峯平均值的閾值
一、描述:
該算法和基於谷底最小值的閾值方法相似,只是最後一步不是取得雙峯之間的谷底值,而是取雙峯的平均值做爲閾值。
二、參考代碼:
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); }
三、效果:
原圖 二值圖 原始直方圖 平滑後的直方圖
5、迭代最佳閾值
一、描述:
該算法先假定一個閾值,而後計算在該閾值下的前景和背景的中心值,當前景和背景中心值得平均值和假定的閾值相同時,則迭代停止,並以此值爲閾值進行二值化。
二、實現過程:
(1)求出圖象的最大灰度值和最小灰度值,分別記爲gl和gu,令初始閾值爲:
(2) 根據閾值T0將圖象分割爲前景和背景,分別求出二者的平均灰度值Ab和Af:
(3) 令
若是Tk=Tk+1,則取Tk爲所求得的閾值,不然,轉2繼續迭代。
三、參考代碼:
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; }
四、效果:
原圖 二值圖 直方圖
6、OSTU大律法
一、描述:
該算法是1979年由日本大津提出的,主要是思想是取某個閾值,使得前景和背景兩類的類間方差最大,matlab中的graythresh便是以該算法爲原理執行的。
二、原理:
關於該算法的原理,網絡上有不少,這裏爲了篇幅有限,不加以贅述。
三、參考代碼:
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; }
四、效果:
該算法對於那些具備平坦的直方圖的圖像具備必定的適應能力。
七、一維最大熵
一、描述:
該算法把信息論中熵的概念引入到圖像中,經過計算閾值分割後兩部分熵的和來判斷閾值是否爲最佳閾值。
二、算法原理
這方面的文章也比較多,留給讀者自行去查找相關資料。
三、參考代碼:
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; }
8、力矩保持法
一、描述:
該算法經過選擇恰當的閾值從而使得二值後的圖像和原始的灰度圖像具備三個相同的初始力矩值。
二、原理:
參考論文:W. Tsai, 「Moment-preserving thresholding: a new approach,」 Comput.Vision Graphics Image Process., vol. 29, pp. 377-393, 1985.
因爲沒法下載到該論文(收費的),僅僅給出從其餘一些資料中找到的公式共享一下。
其中的A\B\C的函數可見代碼部分。
三、參考代碼:
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; }
對於不少圖像,該算法頁能取得比較滿意的結果。
9、基於模糊集理論的閾值
該算法的具體分析可見:基於模糊集理論的一種圖像二值化算法的原理、實現效果及代碼
此法也借用香農熵的概念,該算法通常都能得到較爲理想的分割效果,不論是對雙峯的仍是單峯的圖像。
10、Kittler最小錯誤分類法
因爲精力有限,如下幾種算法僅僅給出算法的論文及相關的代碼。
該算法具體的分析見:
參考代碼:
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法)
參考論文:
參考代碼(未作整理):
// 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; }
12、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; }
十3、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
//圖像的熵=========================================================================================
private void Menu_Entropy_Click(object sender, EventArgs e)
{
if (curBitmap != null)
{
//計算熵
double entropy = GetEntropy(curBitmap, curBitmap.Width, curBitmap.Height);
MessageBox.Show("圖像:"+curFileName +"\n\n"+"圖像熵爲 H = " + entropy);
}
}
//計算熵的函數-----------------------------------------------
public double GetEntropy(Bitmap bmp, int w, int h)
{
int g;
double H = 0.0;
int[] pix = new int[256];
Array.Clear(pix, 0, 256);
//獲取各灰度值的像素個數
for (int i = 0; i < h; i++)
{
for (int j = 0; j < w; j++)
{
Color c = bmp.GetPixel(j, i);
g = (int)(0.299*c.R + 0.587*c.B + 0.114*c.G);
pix[g]++;
}
}
//計算熵
for (int i = 0; i < 256; i++)
if (pix[i] > 0)
H = H + pix[i] * Math.Log10(pix[i]);
H = Math.Log10(w * h) / Math.Log10(2) - H / (w * h * Math.Log10(2));
return H;
}
//灰度平均值=======================================================================================
private void Menu_GrayAverage_Click(object sender, EventArgs e)
{
if (curBitmap != null)
{
double aver = GetGrayAverage(curBitmap, curBitmap.Width, curBitmap.Height);
MessageBox.Show("圖像:" + curFileName +"\n\n"+ "灰度平均值爲 A = " + aver);
}
}
//計算灰度平均值的函數---------------------------------------
public double GetGrayAverage(Bitmap bmp, int w, int h)
{
double sum = 0;
//計算平均值
for (int i = 0; i < w; i++)
for (int j = 0; j < h; j++)
{
Color c = bmp.GetPixel(i, j);
sum = sum + (int)(0.299 * c.R + 0.587 * c.B + 0.114 * c.G);
}
double aver = sum / (w * h);
return aver;
}
//灰度中值=========================================================================================
private void Menu_GrayMid_Click(object sender, EventArgs e)
{
if (curBitmap != null)
{
Bitmap bmp = curBitmap;
//由5 X 5的矩形獲取值
int[] pix = new int[25];
for (int j = 0; j < 5; j++)
{
for (int i = 0; i < 5; i++)
{
pix[i + j * 5] = bmp.GetPixel(100 + i, 100 + j).B;
}
}
pix = MedianSorter(pix,25);
MessageBox.Show("圖像:" + curFileName + "\n" +
"左上角座標爲(100,100)的5X5矩形塊\n灰度中值爲" +
" M = " + pix[12]);
}
}
//排序-------------------------------------------------------
public int[] MedianSorter(int[] data, int m)
{
int tem;
for (int i = m - 1; i >= 1; i--)
for (int j = 1; j <= i; j++)
if (data[j - 1] > data[j])
{
tem = data[j];
data[j] = data[j - 1];
data[j - 1] = tem;
}
return data;
}
//方差=============================================================================================
private void Menu_Variance_Click(object sender, EventArgs e)
{
if (curBitmap != null)
{
double vari = GetVariance(curBitmap, curBitmap.Width, curBitmap.Height);
MessageBox.Show("圖像:" + curFileName + "\n\n" + "方差值 S = " + vari);
}
}
//計算方差值的函數-------------------------------------------
public double GetVariance(Bitmap bmp, int w, int h)
{
double aver, vari =0;
int data;
//求平均值
aver = GetGrayAverage(bmp ,w , h);
//求方差
for (int i = 0; i < w; i++)
{
for (int j = 0; j < h; j++)
{
Color c = bmp.GetPixel(i, j);
data = (int)(0.299 * c.R + 0.587 * c.B + 0.114 * c.G);
vari = vari + (data - aver) * (data - aver);
}
}
vari = vari / (w * h);
return vari;
}