常用圖像閾值分割算法

 閾值分割法可以說是圖像分割中的經典方法，它利用圖像中要提取的目標與背景在灰度上的差異，通過設置閾值來把像素級分成若干類，從而實現目標與背景的分離。

    一般流程：通過判斷圖像中每一個像素點的特徵屬性是否滿足閾值的要求，來確定圖像中的該像素點是屬於目標區域還是背景區域，從而將一幅灰度圖像轉換成二值圖像。

用數學表達式來表示，則可設原始圖像f(x,y),T為閾值，分割圖像時則滿足下式：

                         

按照閾值確定的來源，可以分成：

人工經驗選擇法；直方圖方法；類間方差法；自適應閾值等。

按照閾值確定中的運算範圍，可以分為全局閾值法，局部閾值法。

下面對各種方法進行混合展示：

第一類：全局閾值處理

圖像閾值化分割是一種傳統的最常用的圖像分割方法，因其實現簡單、計算量小、性能較穩定而成為圖像分割中最基本和應用最廣泛的分割技術。它特別適用於目標和背景佔據不同灰度級範圍的圖像。難點在於如何選擇一個合適的閾值實現較好的分割。

0.直接按照經驗給定閾值

我們自己根據需要處理的圖像的先驗知識，對圖像中的目標與背景進行分析。

*PDB文件問題

程序所依賴的所有動態連結庫（dll 文件）也會被編譯，編譯過程中每個 dll 都會產生一個pdb文件，又稱為「符號文件」，是一個存儲數據的信息文件，其包含 dll 庫在編譯過程的某些調試信息，例如程序中所用到的全局變量、局部變量、函數名以及他們的入口地址等。

http://c.biancheng.net/view/474.html

1.基本全局閾值

這個也是直方圖的方法：利用直方圖進行分析，並根據直方圖的波峰和波谷之間的關係，選擇出一個較好的閾值。這樣方法，準確性較高，但是只對於存在一個目標和一個背景的，且兩者對比明顯的圖像，且直方圖是雙峰的那種最有價值。

2. 邊緣改進全局閾值

3.用拉普拉斯邊緣信息改進全局閾值處理

4. 最大類間方差法（OTSU）

    OTSU是一種使用最大類間方差的自動確定閾值的方法。是一種基於全局的二值化算法，它是根據圖像的灰度特性,將圖像分為前景和背景兩個部分。當取最佳閾值時，兩部分之間的差別應該是最大的，在OTSU算法中所採用的衡量差別的標準就是較為常見的最大類間方差。前景和背景之間的類間方差如果越大，就說明構成圖像的兩個部分之間的差別越大，當部分目標被錯分為背景或部分背景被錯分為目標，都會導致兩部分差別變小，當所取閾值的分割使類間方差最大時就意味著錯分概率最小。

 

    記T為前景與背景的分割閾值，前景點數佔圖像比例為w0，平均灰度為u0；背景點數佔圖像比例為w1，平均灰度為u1，圖像的總平均灰度為u，前景和背景圖象的方差g，則有：

                                              

 

聯立上式得：

                                                             

或：

                                                                       

 

        當方差g最大時，可以認為此時前景和背景差異最大，此時的灰度T是最佳閾值。類間方差法對噪聲以及目標大小十分敏感，它僅對類間方差為單峰的圖像產生較好的分割效果。當目標與背景的大小比例懸殊時（例如受光照不均、反光或背景複雜等因素影響），類間方差準則函數可能呈現雙峰或多峰，此時效果不好。

5.最大熵閾值分割法（KSW熵算法）

最大熵閾值分割法和OTSU算法類似，假設將圖像分為背景和前景兩個部分。熵代表信息量，圖像信息量越大，熵就越大，最大熵算法就是找出一個最佳閾值使得背景與前景兩個部分熵之和最大。

 熵的公式：

                                

 

第二類：自適應閾值分割

在分割過程中對圖像上的每個像素都使用了相等的閾值。但在實際情況中，當照明不均勻、有突發噪聲或者背景變化較大時，整幅圖像分割時將沒有合適的單一閾值，如果仍採用單一的閾值去處理每一個像素，可能會將目標和背景區域錯誤劃分。而自適應閾值分割的思想，將圖像中每個像素設置可能不一樣的閾值。

如何確定局部閾值呢？可以計算某個鄰域(局部)的均值、中值、高斯加權平均(高斯濾波)來確定閾值。值得說明的是：如果用局部的均值作為局部的閾值，就是常說的移動平均法

基本原理：

    一種較為簡單的自適應閾值選取方法是對每個像素確定以其自身為中心的一個領域窗口，尋找窗口內像素的最大值與最小值，並取二者的平均值作為閾值，或者將窗口內所有像素的平均值作為閾值，亦或者將窗口內的所有像素的高斯卷積作為閾值。

當背景不同的時候，局部閾值的還是效果好一些。

下面對幾種代表性的方法進行對比：

函數如下:

// 大津閾值

Mat OtsuAlgThreshold(Mat &image)

{

if (image.channels() != 1)

{

cout << "Please input Gray-image!" << endl;

}

int T = 0; //Otsu算法閾值  

double varValue = 0; //類間方差中間值保存

double w0 = 0; //前景像素點數所佔比例  

double w1 = 0; //背景像素點數所佔比例  

double u0 = 0; //前景平均灰度  

double u1 = 0; //背景平均灰度  

double Histogram[256] = { 0 }; //灰度直方圖，下標是灰度值，保存內容是灰度值對應的像素點總數  

uchar *data = image.data;

double totalNum = image.rows*image.cols; //像素總數

for (int i = 0; i < image.rows; i++)

{

for (int j = 0; j < image.cols; j++)

{

if (image.at<uchar>(i, j) != 0) Histogram[data[i*image.step + j]]++;

}

}

int minpos, maxpos;

for (int i = 0; i < 255; i++)

{

if (Histogram[i] != 0)

{

minpos = i;

break;

}

}

for (int i = 255; i > 0; i--)

{

if (Histogram[i] != 0)

{

maxpos = i;

break;

}

}

for (int i = minpos; i <= maxpos; i++)

{

//每次遍歷之前初始化各變量  

w1 = 0;       u1 = 0;       w0 = 0;       u0 = 0;

//***********背景各分量值計算**************************  

for (int j = 0; j <= i; j++) //背景部分各值計算  

{

w1 += Histogram[j];   //背景部分像素點總數  

u1 += j * Histogram[j]; //背景部分像素總灰度和  

}

if (w1 == 0) //背景部分像素點數為0時退出  

{

break;

}

u1 = u1 / w1; //背景像素平均灰度  

w1 = w1 / totalNum; // 背景部分像素點數所佔比例

//***********背景各分量值計算**************************  

//***********前景各分量值計算**************************  

for (int k = i + 1; k < 255; k++)

{

w0 += Histogram[k];  //前景部分像素點總數  

u0 += k * Histogram[k]; //前景部分像素總灰度和  

}

if (w0 == 0) //前景部分像素點數為0時退出  

{

break;

}

u0 = u0 / w0; //前景像素平均灰度  

w0 = w0 / totalNum; // 前景部分像素點數所佔比例  

//***********前景各分量值計算**************************  

//***********類間方差計算******************************  

double varValueI = w0 * w1*(u1 - u0)*(u1 - u0); //當前類間方差計算  

if (varValue < varValueI)

{

varValue = varValueI;

T = i;

}

}

Mat dst;

threshold(image, dst, T, 255, CV_THRESH_OTSU);

return dst;

}

// 自適應閾值

void myadaptive(InputArray _src, OutputArray _dst, double maxValue,

int method, int type, int blockSize, double delta)

{

Mat src = _src.getMat();

CV_Assert(src.type() == CV_8UC1);

CV_Assert(blockSize % 2 == 1 && blockSize > 1);

Size size = src.size();

_dst.create(size, src.type());

Mat dst = _dst.getMat();

if (maxValue < 0)

{

dst = Scalar(0);

return;

}

Mat mean;

if (src.data != dst.data)

mean = dst;

if (method == ADAPTIVE_THRESH_GAUSSIAN_C)

{

GaussianBlur(src, mean, Size(blockSize, blockSize), 0, 0, BORDER_REPLICATE);

}

else if (method == ADAPTIVE_THRESH_MEAN_C)

{

boxFilter(src, mean, src.type(), Size(blockSize, blockSize),

Point(-1, -1), true, BORDER_REPLICATE);

}

else

{

CV_Error(CV_StsBadFlag, "Unknown/unsupported adaptive threshold method");

}

int i, j;

uchar imaxval = saturate_cast<uchar>(maxValue);

int idelta = type == THRESH_BINARY ? cvCeil(delta) : cvFloor(delta);

uchar tab[768];

if (type == CV_THRESH_BINARY)

for (i = 0; i < 768; i++)

tab[i] = (uchar)(i - 255 > -idelta ? imaxval : 0);

else if (type == CV_THRESH_BINARY_INV)

for (i = 0; i < 768; i++)

tab[i] = (uchar)(i - 255 <= -idelta ? imaxval : 0);

else

{

CV_Error(CV_StsBadFlag, "Unknown/unsupported threshold type");

}

if (src.isContinuous() && mean.isContinuous() && dst.isContinuous())

{

size.width *= size.height;

size.height = 1;

}

for (i = 0; i < size.height; i++)

{

const uchar* sdata = src.data + src.step*i;

const uchar* mdata = mean.data + mean.step*i;

uchar* ddata = dst.data + dst.step*i;

for (j = 0; j < size.width; j++)

// 將[-255, 255] 映射到[0, 510]然後查表

ddata[j] = tab[sdata[j] - mdata[j] + 255];

}

}

// 最大熵閾值

Mat EntropySeg(Mat src)

{

int tbHist[256] = { 0 };

int index = 0;

double Property = 0.0;

double maxEntropy = -1.0;

double frontEntropy = 0.0;

double backEntropy = 0.0;

int TotalPixel = 0;

int nCol = src.cols*src.channels();

for (int i = 0; i < src.rows; i++)

{

uchar* pData = src.ptr<uchar>(i);

for (int j = 0; j < nCol; j++)

{

++TotalPixel;

tbHist[pData[j]] += 1;

}

}

for (int i = 0; i < 256; i++)

{

double backTotal = 0;

for (int j = 0; j < i; j++)

{

backTotal += tbHist[j];

}

for (int j = 0; j < i; j++)

{

if (tbHist[j] != 0)

{

Property = tbHist[j] / backTotal;

backEntropy += -Property * logf((float)Property);

}

}

for (int k = i; k < 256; k++)

{

if (tbHist[k] != 0)

{

Property = tbHist[k] / (TotalPixel - backTotal);

frontEntropy += -Property * logf((float)Property);

}

}

if (frontEntropy + backEntropy > maxEntropy)

{

maxEntropy = frontEntropy + backEntropy;

index = i;

}

frontEntropy = 0.0;

backEntropy = 0.0;

}

Mat dst;

threshold(src, dst, index, 255, 0);

return dst;

}

//迭代閾值

Mat IterationThreshold(Mat src)

{

int width = src.cols;

int height = src.rows;

int hisData[256] = { 0 };

for (int j = 0; j < height; j++)

{

uchar* data = src.ptr<uchar>(j);

for (int i = 0; i < width; i++)

hisData[data[i]]++;

}

int T0 = 0;

for (int i = 0; i < 256; i++)

{

T0 += i * hisData[i];

}

T0 /= width * height;

int T1 = 0, T2 = 0;

int num1 = 0, num2 = 0;

int T = 0;

while (1)

{

for (int i = 0; i < T0 + 1; i++)

{

T1 += i * hisData[i];

num1 += hisData[i];

}

if (num1 == 0)

continue;

for (int i = T0 + 1; i < 256; i++)

{

T2 += i * hisData[i];

num2 += hisData[i];

}

if (num2 == 0)

continue;

T = (T1 / num1 + T2 / num2) / 2;

if (T == T0)

break;

else

T0 = T;

}

Mat dst;

threshold(src, dst, T, 255, 0);

return dst;

}

調試過程中：

*PDB文件問題

程序所依賴的所有動態連結庫（dll 文件）也會被編譯，編譯過程中每個 dll 都會產生一個pdb文件，又稱為「符號文件」，是一個存儲數據的信息文件，其包含 dll 庫在編譯過程的某些調試信息，例如程序中所用到的全局變量、局部變量、函數名以及他們的入口地址等。

http://c.biancheng.net/view/474.html

* 內存溢出問題通過修改絕對路徑或者將cpp與圖片放在同一層解決。

可見廣泛使用的大津法其實效果不是太好， 二值化後一些細節部分難以通過前景和背景的差異表現出來。最大熵法通過信息量最大化進行分割，有良好的的效果，但是在青蛙的帽子上面有誤差。局部自適應閾值法對深色的蝴蝶翅膀提取效果不好。

這兩個都效果還好，最大熵使得全局信息門限值最大化，但是畢竟是全局閾值的方法；自適應方法很明顯窗口大小太重要了，平緩背景下分割效果不好，對窗口大小進行調整：

OpenCV提供的API，對窗口大小進行調整後進一步實驗：

說明下各參數：

InputArray src：源圖像

OutputArray dst：輸出圖像，與源圖像大小一致

int adaptiveMethod：在一個鄰域內計算閾值所採用的算法，有兩個取值，分別為 ADAPTIVE_THRESH_MEAN_C 和 ADAPTIVE_THRESH_GAUSSIAN_C 。

ADAPTIVE_THRESH_MEAN_C的計算方法是計算出領域的平均值再減去第七個參數double C的值。

ADAPTIVE_THRESH_GAUSSIAN_C的計算方法是計算出領域的高斯均值再減去第七個參數double C的值。

int thresholdType：這是閾值類型，只有兩個取值，分別為 THRESH_BINARY 和THRESH_BINARY_INV  具體的請看官方的說明，這裡不多做解釋。

int blockSize：adaptiveThreshold的計算單位是像素的鄰域塊，這是局部鄰域大小，3、5、7等。

double C：這個參數實際上是一個偏移值調整量，用均值和高斯計算閾值後，再減或加這個值就是最終閾值

將閾值從15（左圖）改成35（右圖）後好多了，但是又有點過了.怎麼樣才能自動確定出最佳窗口大小呢，如果 運行成本低又適應不同粒度區域的方法的話就很好了..

* 不同粒度程度用不同的窗口，是不是會出現邊緣提取不連接的情況？林洪磊師兄的那個長方形窗口好處是啥/