一维OTSU法、最小交叉熵法、二维OTSU法及C++源码
1.最大類間方差法(Otsu法)? ? ? ?
? ? ? ? 該算法是日本人Otsu提出的一種動態閾值分割算法。它的主要思想是按照灰度特性將圖像劃分為背景和目標2部分,劃分依據為選取門限值,使得背景和目標之間的方差最大。(背景和目標之間的類間方差越大,說明這兩部分的差別越大,當部分目標被錯劃分為背景或部分背景錯劃分為目標都會導致這兩部分差別變小。因此,使用類間方差最大的分割意味著錯分概率最小。)這是該方法的主要思路。其主要的實現原理為如下:
?? ? ?1)建立圖像灰度直方圖(共有L個灰度級,每個出現概率為p)
?? ? ?2)計算背景和目標的出現概率,計算方法如下:
?? ? ?上式中假設t為所選定的閾值,A代表背景(灰度級為0~N),根據直方圖中的元素可知,Pa為背景出現的概率,同理B為目標,Pb為目標出現的概率。
?? ? ?3)計算A和B兩個區域的類間方差如下:
?? ?? ?第一個表達式分別計算A和B區域的平均灰度值;
?? ? ?第二個表達式計算灰度圖像全局的灰度平均值;
?? ? ?第三個表達式計算A、B兩個區域的類間方差。?? ? ?4)以上幾個步驟計算出了單個灰度值上的類間方差,因此最佳分割門限值應該是圖像中能夠使得A與B的類間灰度方差最大的灰度值。在程序中需要對每個出現的灰度值據此進行尋優。
? ? ? ? 以下對兩個不同類型的圖像采用最大類間方差法尋找最優閾值進行分割:
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
? ? ? ? C++源碼:
#include "stdio.h" #include "cv.h" #include "highgui.h" #include "Math.h" #include <tchar.h> int Otsu(IplImage* src);int _tmain(int argc, _TCHAR* argv[]) {IplImage* img = cvLoadImage("E:\\11.jpg", 0);IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);int threshold = Otsu(img);cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);cvNamedWindow("img", 1);cvShowImage("img", dst);cvWaitKey(-1);cvReleaseImage(&img);cvReleaseImage(&dst);cvDestroyWindow("dst");return 0; } int Otsu(IplImage* src) {int height = src->height;int width = src->width;double nHistogram[256]; //灰度直方圖 double dVariance[256]; //類間方差 int N = height*width; //總像素數 for (int i = 0; i<256; i++){nHistogram[i] = 0.0;dVariance[i] = 0.0;}for (int i = 0; i<height; i++){unsigned char* nData = (unsigned char*)src->imageData + src->widthStep * i;for (int j = 0; j < width; j++){nHistogram[*nData++]++; //建立直方圖}}double Pa = 0.0; //背景出現概率 double Pb = 0.0; //目標出現概率 double Wa = 0.0; //背景平均灰度值 double Wb = 0.0; //目標平均灰度值 double W0 = 0.0; //全局平均灰度值 double dData1 = 0.0, dData2 = 0.0;int nThreshold;for (int i = 0; i<256; i++) //計算全局平均灰度 {nHistogram[i] /= N;W0 += i*nHistogram[i];}for (int i = 0; i<256; i++) //對每個灰度值計算類間方差 {Pa += nHistogram[i];Pb = 1 - Pa;dData1 += i*nHistogram[i];dData2 = W0 - dData1;Wa = dData1 / Pa;Wb = dData2 / Pb;dVariance[i] = (Pa*Pb* pow((Wb - Wa), 2));} //遍歷每個方差,求取類間最大方差所對應的灰度值 double temp = 0.0;for (int i = 0; i<256; i++){if (dVariance[i]>temp){temp = dVariance[i];nThreshold = i;}}printf("threshold = %d\n", nThreshold);return nThreshold;} ? ? ? ? ? ? ? ? ? ? ? ??
? ? ? ? 以上是 采用 最大類間方差法進行最優閾值分割的結果,可以看出原圖采用不同尺寸的彩色圖像和灰度圖像,處理之后的效果有顯著的差別。 ? ? ? ? ? ? ? ? ? ? ? ? ??
2.最小交叉熵法
? ? ? 這種方法與類間最大方差很相似,是由Li和Lee應用了信息論中熵理論發展而來。首先簡要介紹交叉熵的概念。
?? ? ?對于兩個分布P和Q,定義其信息交叉熵D如下:
?? ? ?這代表的物理意義是兩個分布之間信息理論距離,另外一種理解是,將分布P變為Q后所帶來的信息變化。那么對于圖像分割來說,如果要用分割圖像來替換原來的圖像,最優的分割依據應該就是使得兩幅圖像之間的交叉熵最小。以下對最小交叉熵法的過程進行簡要總結。
?? ? ?可以假設上文的P為源圖像的灰度分布,Q為所得到的分割圖像的灰度分布,其中:
?? ? ?
?? ? ?上式中H為統計直方圖;
?? ? ?N為圖像總的像素點數;
?? ? ?L為源圖像總的灰度級數;
?? ? ?P代表源圖像,其每個元素代表每個灰度級上的灰度分布(平均灰度值);
?? ? ?Q為分割后的二值圖像,兩個u分別代表兩個分割后的區域的平均灰度值,其中t為分割圖像所采用的閾值。
?? ? ?根據以上定義,以每個灰度級上的灰度和為計算量,可以很容易根據交叉熵的公式,推導出P和Q之間的交叉熵定量表達式:
?? ? ?根據上文所述思路,使得D最小的t即為最小交叉熵意義下的最優閾值。
? ? ? ? ? ? ? ? ???
? ? ? ? 以下采用最小交叉熵法對彩色圖像和灰度圖進行閾值化分割:
? ? ? ? 源碼:
int MiniCross(IplImage* src) {int height = src->height;int width = src->width;double dHistogram[256]; //灰度直方圖 double dEntropy[256]; //每個像素的交叉熵 int N = height*width; //總像素數 for (int i = 0; i<256; i++){dHistogram[i] = 0.0;dEntropy[i] = 0.0;}for (int i = 0; i<height; i++){unsigned char* nData = (unsigned char*)src->imageData + src->widthStep * i;for (int j = 0; j < width; j++){dHistogram[*nData++]++; //建立直方圖}}double Pa = 0.0; //區域1平均灰度值 double Pb = 0.0; //區域2平均灰度值 double P0 = 0.0; //全局平均灰度值 double Wa = 0.0; //第一部分熵 double Wb = 0.0; //第二部分的熵 double dData1 = 0.0, dData2 = 0.0; //中間值 double dData3 = 0.0, dData4 = 0.0; //中間值 int nThreshold;for (int i = 0; i<256; i++) //計算全局平均灰度 {dHistogram[i] /= N;P0 += i*dHistogram[i];}for (int i = 0; i<256; i++){Wa = Wb = dData1 = dData2 = dData3 = dData4 = Pa = Pb = 0.0;for (int j = 0; j<256; j++){if (j <= i){dData1 += dHistogram[j];dData2 += j*dHistogram[j];}else{dData3 += dHistogram[j];dData4 += j*dHistogram[j];}}Pa = dData2 / dData1;Pb = dData4 / dData3;for (int j = 0; j<256; j++){if (j <= i){if ((Pa != 0) && (dHistogram[j] != 0)){double d1 = log(dHistogram[j] / Pa);Wa += j*dHistogram[j] * d1 / log(2);}}else{if ((Pb != 0) && (dHistogram[j] != 0)){double d2 = log(dHistogram[j] / Pb);Wb += j*dHistogram[j] * d2 / log(2);}}}dEntropy[i] = Wa + Wb;}//遍歷熵值,求取最小交叉熵所對應的灰度值 double temp = dEntropy[0];for (int i = 1; i<256; i++){if (dEntropy[i]<temp){temp = dEntropy[i];nThreshold = i;}}printf("threshold = %d\n", nThreshold);return nThreshold; }? ? ? ? ? ? ? ? ? ?
? ? ? ? 以上是最小交叉熵法尋找最優閾值后處理得到的分割圖像,可以看出,彩色圖像處理后輪廓信息基本已經消失,灰度圖像經最小交叉熵處理后的效果與采用最大類間方差處理后得到的分割圖像效果相似。
3.二維OTSU法
? ? ? ? 這種方法是對類間最大方差法的擴展,將其從求兩個一維分布最大類間方差擴充為求解類間離散度矩陣的跡的最大值,考慮像素點灰度級的基礎上增加了對像素點鄰域平均像素值的考慮。
?? ? ?以下按照本人的理解對該方法的思路以及推倒過程進行分析:
?? ? ?1)首先需要建立二維的灰度統計直方圖P(f, g);
?? ? ?圖像的灰度級為L級,那么其每個像素點的8鄰域灰度平均值的灰度級也為L級,據此來構建直方圖P。二維統計直方圖的橫軸為每個像素點的灰度值f(i, j),縱坐標為同一個點對應的鄰域平均值g(i, j) 其中(0≤i<height, 0≤j<width, 0≤f(i, j)<L),而所對應的P(f,g)整幅圖像中灰度值為f,鄰域灰度均值為g的點的統計值占總像素點的比例(即為灰度值出現的聯合概率密度)。其中P的每個元素滿足如下公式:
?? ? ?n為整幅圖像中灰度值為f,鄰域灰度均值為g的點的統計值;
?? ? ?N為圖像總的像素點個數;
?? ? ?2)對于下圖所示的二維統計直方圖,t代表橫坐標(灰度值),s代表縱坐標(像素點鄰域的灰度均值)
?? ? ?對已圖像中的閾值點(t,s)來說,其灰度值t和其鄰域內的灰度均值s不應該相差太多,如果t比s大很多(點位于上圖中的II區域),說明像素的灰度值遠遠大于其臨域的灰度均值,故而該點很可能是噪聲點,反之如果t比s小很多(點位于途中的IV區域),即該點的像素值比其臨域均值小很多,則說明是一個邊緣點。據此我們在進行背景前景分割的時候忽略這些干擾因素,認為這兩個區域內Pi,j=0。剩下的I區域和III區域則分別代表了前景和背景。以下據此來推導對于選定的閾值(t, s),進行離散度判據的最優推導。
?? ? ?3)推導閾值(t, s)點處的離散度矩陣判據
?? ? ?根據上文分析可知,由閾值(t, s)所分割的前景和背景出現的概率如下:
?? ? ?定義兩個中間變量,方便下面推導:
?? ? ?據此,這兩部分的灰度均值向量可以推導如下(兩個分量分別根據灰度值以及每個點的灰度均值計算):
?? ? ?整幅圖像的灰度均值向量為:
?? ? ?與一維的大津法一樣的思路,推導類間方差,這里是二維因此要用矩陣形式。參考一維法,可同樣定義類間“方差”矩陣如下:
?? ? ?為了在實現的時候,容易判斷出這樣一個矩陣的“最大值”,因此數學中采用矩陣的跡(對角線之和)來衡量矩陣的“大小”。因此以該矩陣的跡作為離散度測度,推導如下:
?? ? ?這樣就可以通過求解使得這個參數最大時的(t,s)即為所求得的最佳閾值組合。
?? ? ?以下為具體算法實現過程:
?? ? ?1)建立二維直方圖
?? ? ?2)對直方圖進行遍歷,計算每個(t,s)組合所得到的矩陣離散度,也就是一維大津法中所謂的最大類間方差。
?? ? ?3)求得使“類間方差”最大的(t,s),由于t代表灰度值,s代表改點在其鄰域內的灰度均值,因此本人認為在選擇閾值時可以選擇s為最佳,當然用t也可以,因為從求解結果可以看出,這這個數值往往很接近。
?? ? ?具體的實現代碼如下:
/***************************************************************************** * 函數名稱: cvOtsu2D() * 函數參數: CvMat* pGrayMat:灰度圖形相對應的矩陣 * 返回值: int nThreshold * 函數說明:實現灰度圖的二值化分割——最大類間方差法(二維Otsu算法) * 備注:在構建二維直方圖的時候,采用灰度點的3*3鄰域均值 ******************************************************************************/#include <opencv2\opencv.hpp> #include <highgui.h> #include <cv.h> #include <iostream> #include <stdio.h> #include<opencv\highgui.h> #include<opencv2\core\core.hpp> #include "math.h" #include <tchar.h> int cvOtsu2D(CvMat *pGrayMat) {double dHistogram[256][256]; //建立二維灰度直方圖 double dTrMatrix = 0.0; //離散矩陣的跡int height = pGrayMat->rows;int width = pGrayMat->cols;int N = height*width; //總像素數int i, j;for (i = 0; i < 256; i++){for (j = 0; j < 256; j++)dHistogram[i][j] = 0.0; //初始化變量 }for (i = 0; i < height; i++){for (j = 0; j < width; j++){unsigned char nData1 = (unsigned char)cvGetReal2D(pGrayMat, i, j);//當前的灰度值 unsigned char nData2 = 0;int nData3 = 0;//注意9個值相加可能超過一個字節 for (int m = i - 1; m <= i + 1; m++){for (int n = j - 1; n <= j + 1; n++){if ((m >= 0) && (m < height) && (n >= 0) && (n < width))nData3 += (unsigned char)cvGetReal2D(pGrayMat, m, n); //當前的灰度值 }}nData2 = (unsigned char)(nData3 / 9); //對于越界的索引值進行補零,鄰域均值 dHistogram[nData1][nData2]++;}}for (i = 0; i < 256; i++)for (j = 0; j < 256; j++)dHistogram[i][j] /= N; //得到歸一化的概率分布 double Pai = 0.0; //目標區均值矢量i分量 double Paj = 0.0; //目標區均值矢量j分量 double Pbi = 0.0; //背景區均值矢量i分量 double Pbj = 0.0; //背景區均值矢量j分量 double Pti = 0.0; //全局均值矢量i分量 double Ptj = 0.0; //全局均值矢量j分量 double W0 = 0.0; //目標區的聯合概率密度 double W1 = 0.0; //背景區的聯合概率密度 double dData1 = 0.0;double dData2 = 0.0;double dData3 = 0.0;double dData4 = 0.0; //中間變量 int nThreshold_s = 0;int nThreshold_t = 0;int nThreshold;double temp = 0.0; //尋求最大值 for (i = 0; i < 256; i++){for (j = 0; j < 256; j++){Pti += i*dHistogram[i][j];Ptj += j*dHistogram[i][j];}}for (i = 0; i < 256; i++){for (j = 0; j < 256; j++){W0 += dHistogram[i][j];dData1 += i*dHistogram[i][j];dData2 += j*dHistogram[i][j];W1 = 1 - W0;dData3 = Pti - dData1;dData4 = Ptj - dData2;Pai = dData1 / W0;Paj = dData2 / W0;Pbi = dData3 / W1;Pbj = dData4 / W1; // 得到兩個均值向量,用4個分量表示 dTrMatrix = ((W0 * Pti - dData1) * (W0 * Pti - dData1) + (W0 * Ptj - dData2) * (W0 * Ptj - dData2)) / (W0 * W1);if (dTrMatrix > temp){temp = dTrMatrix;nThreshold_s = i;nThreshold_t = j;}}}nThreshold = (nThreshold_s + nThreshold_t) / 2;//返回閾值,有多種形式,可以單獨某一個,也可 //是兩者的均值 printf("nThreshold= %d\n", nThreshold);return nThreshold; } int main() {IplImage *img = cvLoadImage("E:\\12.jpg", 0);//加載圖片cvNamedWindow("Miss", CV_WINDOW_AUTOSIZE);//創建顯示窗口cvShowImage("Miss", img);//顯示圖片cvNamedWindow("二維Otsu", CV_WINDOW_AUTOSIZE);CvMat *imgMat = cvCreateMat(img->height, img->width, CV_8UC1);//創建矩陣,其type為CV_8UC1cvConvert(img, imgMat);//將圖像轉換成矩陣形式int t = cvOtsu2D(imgMat);//求二維閾值分割的閾值cvThreshold(img, img, t, 255, CV_THRESH_BINARY);//用求得的閾值分割圖像cvShowImage("二維Otsu", img);cvWaitKey(0);//保持圖像的顯示窗口cvReleaseImage(&img);//釋放圖像資源cvReleaseMat(&imgMat);//釋放矩陣資源cvDestroyAllWindows();//銷毀窗口return 0; }? ? ? ? ? ? ? ? ? ? ??
? ? ? ?
? ? ? ? ? 此圖的效果看不出一維Otsu法和二維Otsu法的區別。
本文理論部分參考自:http://blog.csdn.net/likezhaobin/article/details/6915755
總結
以上是生活随笔為你收集整理的一维OTSU法、最小交叉熵法、二维OTSU法及C++源码的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 美国股市泡沫有多严重
- 下一篇: memset()函数用法