灰度影象閾值化分割常見方法總結及VC實現
在影象處理領域，二值影象運算量小，並且能夠體現影象的關鍵特徵，因此被廣泛使用。將灰度影象變為二值影象的常用方法是選定閾值，然後將待處理影象的每個畫素點進行單點處理，即將其灰度值與所設定的門限進行比對，從而得到二值化的黑白圖。這樣一種方式因為其直觀性以及易於實現，已經在影象分割領域處於中心地位。本文主要對最近一段時間作者所學習的閾值化影象分割演算法進行總結，全文描述了作者對每種演算法的理解，並基於OpenCV和VC6.0對這些演算法進行了實現。最終將原始碼公開，希望大家一起進步。（本文的程式碼暫時沒有考慮執行效率問題）
首先給出待分割的影象如下：
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的類間灰度方差最大的灰度值。在程式中需要對每個出現的灰度值據此進行尋優。
本人的VC實現程式碼如下。閾值分割結果如下圖，求解所得的閾值為116./***************************************************************************** * * \函式名稱： * OneDimentionOtsu() * * \輸入引數: * pGrayMat: 二值影象資料 * width: 圖形尺寸寬度 * height: 圖形尺寸高度 * nTlreshold: 經過演算法處理得到的二值化分割閾值 * \返回值: * 無 * \函式說明:實現灰度圖的二值化分割——最大類間方差法（Otsu演算法，俗稱大津演算法） * ****************************************************************************/ void CBinarizationDlg::OneDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold) { 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(i=0; i<height; i++) { for(int j=0; j<width; j++) { unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, 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; for(i=0; i<256; i++) //計算全域性平均灰度 { nHistogram[i] /= N; W0 += i*nHistogram[i]; } for(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(i=0; i<256; i++) { if(dVariance[i]>temp) { temp = dVariance[i]; nThreshold = i; } } }
2、一維交叉熵值法
這種方法與類間最大方差很相似，是由Li和Lee應用了資訊理論中熵理論發展而來。首先簡要介紹交叉熵的概念。
對於兩個分佈P和Q，定義其資訊交叉熵D如下：
這代表的物理意義是兩個分佈之間資訊理論距離，另外一種理解是，將分佈P變為Q後所帶來的資訊變化。那麼對於影象分割來說，如果要用分割影象來替換原來的影象，最優的分割依據應該就是使得兩幅影象之間的交叉熵最小。以下對最小交叉熵法的過程進行簡要總結。
可以假設上文的P為源影象的灰度分佈，Q為所得到的分割影象的灰度分佈，其中：
上式中H為統計直方圖；
N為影象總的畫素點數；
L為源影象總的灰度級數；
P代表源影象，其每個元素代表每個灰度級上的灰度分佈（平均灰度值）；
Q為分割後的二值影象，兩個u分別代表兩個分割後的區域的平均灰度值，其中t為分割影象所採用的閾值。
根據以上定義，以每個灰度級上的灰度和為計算量，可以很容易根據交叉熵的公式，推匯出P和Q之間的交叉熵定量表達式：
根據上文所述思路，使得D最小的t即為最小交叉熵意義下的最優閾值。
作者VC實現程式碼如下。/*****************************************************************************
*
* \函式名稱：
* MiniCross()
*
* \輸入引數:
* pGrayMat: 二值影象資料
* width: 圖形尺寸寬度
* height: 圖形尺寸高度
* nTlreshold: 經過演算法處理得到的二值化分割閾值
* \返回值:
* 無
* \函式說明:實現灰度圖的二值化分割——最小交叉熵演算法
*
****************************************************************************/
void CBinarizationDlg::MiniCross(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)
{
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(i=0; i<height; i++)
{
for(int j=0; j<width; j++)
{
unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, 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; //中間值
for(i=0; i<256; i++) //計算全域性平均灰度
{
dHistogram[i] /= N;
P0 += i*dHistogram[i];
}
for(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(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(i=1; i<256; i++)
{
if(dEntropy[i]<temp)
{
temp = dEntropy[i];
nThreshold = i;
}
}
}閾值分割結果如下圖，求解所得的閾值為106.
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也可以，因為從求解結果可以看出，這這個數值往往很接近。
具體的實現程式碼如下：
/*****************************************************************************
*
* \函式名稱：
* TwoDimentionOtsu()
*
* \輸入引數:
* pGrayMat: 二值影象資料
* width: 圖形尺寸寬度
* height: 圖形尺寸高度
* nTlreshold: 經過演算法處理得到的二值化分割閾值
* \返回值:
* 無
* \函式說明:實現灰度圖的二值化分割——最大類間方差法（二維Otsu演算法）
* \備註：在構建二維直方圖的時候，採用灰度點的3*3鄰域均值
******************************************************************************/
void CBinarizationDlg::TwoDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)
{
double dHistogram[256][256]; //建立二維灰度直方圖
double dTrMatrix = 0.0; //離散矩陣的跡
int N = height*width; //總畫素數
for(int i=0; i<256; i++)
{
for(int j=0; j<256; j++)
dHistogram[i][j] = 0.0; //初始化變數
}
for(i=0; i<height; i++)
{
for(int j=0; j<width; j++)
{
unsigned char nData1 = (unsigned char)cvmGet(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)cvmGet(pGrayMat, m, n); //當前的灰度值
}
}
nData2 = (unsigned char)(nData3/9); //對於越界的索引值進行補零,鄰域均值
dHistogram[nData1][nData2]++;
}
}
for(i=0; i<256; i++)
for(int 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;
double temp = 0.0; //尋求最大值
for(i=0; i<256; i++)
{
for(int j=0; j<256; j++)
{
Pti += i*dHistogram[i][j];
Ptj += j*dHistogram[i][j];
}
}
for(i=0; i<256; i++)
{
for(int 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;
/* W1=dData3=dData4=0.0; //對內迴圈的資料進行初始化
for(int s=i+1; s<256; s++)
{
for(int t=j+1; t<256; t++)
{
W1 += dHistogram[s][t];
dData3 += s*dHistogram[s][t]; //方法2
dData4 += t*dHistogram[s][t]; //也可以增加迴圈進行計算
}
}*/
Pai = dData1/W0;
Paj = dData2/W0;
Pbi = dData3/W1;
Pbj = dData4/W1; // 得到兩個均值向量，用4個分量表示
dTrMatrix = ((W0*Pti-dData1)*(W0*Pti-dData1)+(W0*Ptj-dData1)*(W0*Ptj-dData2))/(W0*W1);
if(dTrMatrix > temp)
{
temp = dTrMatrix;
nThreshold_s = i;
nThreshold_t = j;
}
}
}
nThreshold = nThreshold_t; //返回結果中的灰度值
//nThreshold = 100;
}閾值分割結果如下圖，求解所得的閾值為114.（s=114，t=117）
參考文獻
[1] Nobuyuki Otsu. A Threshold SelectionMethod from Gray-Level Histograms.
[2] C. H. Li and C. K. Lee. Minimum CrossEntropy Thresholding
[3] Jian Gong, LiYuan Liand WeiNan Chen. Fast Recursive Algorithms For Two-Dimensional Thresholding.
由於網速不給力，這篇文章公式以及圖片上傳好慢，用了作者好長時間，不容易啊。因此如果各位看官要轉載，麻煩請註明轉載地址，以算是對作者辛苦的回報。同時希望能夠多多留言，給作者積分。