1. 程式人生 > >七種常見閾值分割程式碼(Otsu、最大熵、迭代法、自適應閥值、手動、迭代法、基本全域性閾值法)

七種常見閾值分割程式碼(Otsu、最大熵、迭代法、自適應閥值、手動、迭代法、基本全域性閾值法)

整理了一些主要的分割方法,以後用省的再查,其中大部分為轉載資料,轉載連結見資料;

一、工具:VC+OpenCV

二、語言:C++

三、原理

   (1) otsu法(最大類間方差法,有時也稱之為大津演算法)使用的是聚類的思想,把影象的灰度數按灰度級分成2個部分,使得兩個部分之間的灰度值差異最大,每個部分之間的灰度差異最小,通過方差的計算來尋找一個合適的灰度級別 來劃分。 所以 可以在二值化的時候 採用otsu演算法來自動選取閾值進行二值化。otsu演算法被認為是影象分割中閾值選取的最佳演算法,計算簡單,不受影象亮度和對比度的影響。因此,使類間方差最大的分割意味著錯分概率最小。

設t為設定的閾值。

wo: 分開後 前景畫素點數佔影象的比例

uo: 分開後 前景畫素點的平均灰度

w1:分開後 被景畫素點數佔影象的比例

u1: 分開後 被景畫素點的平均灰度

u=w0*u0 + w1*u1 :影象總平均灰度

從L個灰度級遍歷t,使得t為某個值的時候,前景和背景的方差最大, 則 這個 t 值便是我們要求得的閾值。

其中,方差的計算公式如下:

g=wo * (uo - u) * (uo - u) + w1 * (u1 - u) * (u1 - u)

[ 此公式計算量較大,可以採用: g = wo * w1 * (uo - u1) * (uo - u1) ]

由於otsu演算法是對影象的灰度級進行聚類,so 在執行otsu演算法之前,需要計算該影象的灰度直方圖。

(2)迭代法原理:迭代選擇法是首先猜測一個初始閾值,然後再通過對影象的多趟計算對閾值進行改進的過程。重複地對影象進行閾值操作,將影象分割為物件類和背景類,然後來利用每一個類中的灰階級別對閾值進行改進。

影象閾值分割---迭代演算法

   1 .處理流程:

       1.為全域性閾值選擇一個初始估計值T(影象的平均灰度)。        2.用T分割影象。產生兩組畫素:G1有灰度值大於T的畫素組成,G2有小於等於T畫素組成。        3.計算G1和G2畫素的平均灰度值m1和m2;        4.計算一個新的閾值:T = (m1 + m2) / 2;        5.重複步驟2和4,直到連續迭代中的T值間的差小於一個預定義引數為止。

       適合影象直方圖有明顯波谷

(3)

直方圖雙峰法(mode 法)

  Prewitt 等人於六十年代中期提出的直方圖雙峰法(也稱 mode 法) 是典型的全域性單閾值分割方法。該方法的基本思想是:假設影象中有明顯的目標和背景,則其灰度直方圖呈雙峰分佈,當灰度級直方圖具有雙峰特性時,選取兩峰之間的谷對應的灰度級作為閾值。如果背景的灰度值在整個影象中可以合理地看作為恆定,而且所有物體與背景都具有幾乎相同的對比度,那麼,選擇一個正確的、固定的全域性閾值會有較好的效果。例如圖3.1所示:

圖3.1原始灰度影象

圖3.2灰度直方圖

選定閾值M為100

演算法實現:找到第一個峰值和第二個峰值, 再找到第一和第二個峰值之間的谷值,谷值就是那個閥值了。

四、程式

主程式(核心部分) 

閾值分割
1 /*===============================影象分割=====================================*/
2 /*---------------------------------------------------------------------------*/
3 /*手動設定閥值*/
4 IplImage* binaryImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);
5 cvThreshold(smoothImgGauss,binaryImg,71,255,CV_THRESH_BINARY); 
6 cvNamedWindow("cvThreshold", CV_WINDOW_AUTOSIZE );
7 cvShowImage( "cvThreshold", binaryImg );
8 //cvReleaseImage(&binaryImg); 
9  /*---------------------------------------------------------------------------*/
10 /*自適應閥值 //計算像域鄰域的平均灰度,來決定二值化的值*/
11 IplImage* adThresImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);
12 double max_value=255;
13 int adpative_method=CV_ADAPTIVE_THRESH_GAUSSIAN_C;//CV_ADAPTIVE_THRESH_MEAN_C
14  int threshold_type=CV_THRESH_BINARY;
15 int block_size=3;//閾值的象素鄰域大小
16  int offset=5;//視窗尺寸
17   cvAdaptiveThreshold(smoothImgGauss,adThresImg,max_value,adpative_method,threshold_type,block_size,offset);
18 cvNamedWindow("cvAdaptiveThreshold", CV_WINDOW_AUTOSIZE );
19 cvShowImage( "cvAdaptiveThreshold", adThresImg );
20 cvReleaseImage(&adThresImg);
21 /*---------------------------------------------------------------------------*/
22 /*最大熵閥值分割法*/ 
23 IplImage* imgMaxEntropy = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
24 MaxEntropy(smoothImgGauss,imgMaxEntropy);
25 cvNamedWindow("MaxEntroyThreshold", CV_WINDOW_AUTOSIZE );
26 cvShowImage( "MaxEntroyThreshold", imgMaxEntropy );//顯示影象
27   cvReleaseImage(&imgMaxEntropy ); 
28 /*---------------------------------------------------------------------------*/
29 /*基本全域性閥值法*/
30 IplImage* imgBasicGlobalThreshold = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
31 cvCopyImage(srcImgGrey,imgBasicGlobalThreshold);
32 int pg[256],i,thre; 
33 for (i=0;i<256;i++) pg[i]=0;
34 for (i=0;i<imgBasicGlobalThreshold->imageSize;i++) // 直方圖統計
35   pg[(BYTE)imgBasicGlobalThreshold->imageData[i]]++; 
36 thre = BasicGlobalThreshold(pg,0,256); // 確定閾值
37   cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<thre<<endl;//輸出顯示閥值
38   cvThreshold(imgBasicGlobalThreshold,imgBasicGlobalThreshold,thre,255,CV_THRESH_BINARY); // 二值化 
39   cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE );
40 cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//顯示影象
41   cvReleaseImage(&imgBasicGlobalThreshold);
42 /*---------------------------------------------------------------------------*/
43 /*OTSU*/
44 IplImage* imgOtsu = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
45 cvCopyImage(srcImgGrey,imgOtsu);
46 int thre2;
47 thre2 = otsu2(imgOtsu);
48 cout<<"The Threshold of this Image in Otsu is:"<<thre2<<endl;//輸出顯示閥值
49 cvThreshold(imgOtsu,imgOtsu,thre2,255,CV_THRESH_BINARY); // 二值化 
50 cvNamedWindow("imgOtsu", CV_WINDOW_AUTOSIZE );
51 cvShowImage( "imgOtsu", imgOtsu);//顯示影象 
52 cvReleaseImage(&imgOtsu);
53 /*---------------------------------------------------------------------------*/
54 /*上下閥值法:利用正態分佈求可信區間*/
55 IplImage* imgTopDown = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
56 cvCopyImage(srcImgGrey,imgTopDown);
57 CvScalar mean ,std_dev;//平均值、 標準差
58 double u_threshold,d_threshold;
59 cvAvgSdv(imgTopDown,&mean,&std_dev,NULL); 
60 u_threshold = mean.val[0] +2.5* std_dev.val[0];//上閥值
61 d_threshold = mean.val[0] -2.5* std_dev.val[0];//下閥值
62 //u_threshold = mean + 2.5 * std_dev; //錯誤
63 //d_threshold = mean - 2.5 * std_dev;
64 cout<<"The TopThreshold of this Image in TopDown is:"<<d_threshold<<endl;//輸出顯示閥值
65 cout<<"The DownThreshold of this Image in TopDown is:"<<u_threshold<<endl;
66 cvThreshold(imgTopDown,imgTopDown,d_threshold,u_threshold,CV_THRESH_BINARY_INV);//上下閥值
67 cvNamedWindow("imgTopDown", CV_WINDOW_AUTOSIZE );
68 cvShowImage( "imgTopDown", imgTopDown);//顯示影象 
69 cvReleaseImage(&imgTopDown);
70 /*---------------------------------------------------------------------------*/
71 /*迭代法*/
72 IplImage* imgIteration = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
73 cvCopyImage(srcImgGrey,imgIteration);
74 int thre3,nDiffRec;
75 thre3 =DetectThreshold(imgIteration, 100, nDiffRec);
76 cout<<"The Threshold of this Image in imgIteration is:"<<thre3<<endl;//輸出顯示閥值
77 cvThreshold(imgIteration,imgIteration,thre3,255,CV_THRESH_BINARY_INV);//上下閥值
78 cvNamedWindow("imgIteration", CV_WINDOW_AUTOSIZE );
79 cvShowImage( "imgIteration", imgIteration);
80 cvReleaseImage(&imgIteration);
迭代
1 /*======================================================================*/
2 /* 迭代法*/
3 /*======================================================================*/
4 // nMaxIter:最大迭代次數;nDiffRec:使用給定閥值確定的亮區與暗區平均灰度差異值
5 int DetectThreshold(IplImage*img, int nMaxIter, int& iDiffRec) //閥值分割:迭代法
6 {
7 //影象資訊
8 int height = img->height;
9 int width = img->width;
10 int step = img->widthStep/sizeof(uchar);
11 uchar *data = (uchar*)img->imageData;
12 
13 iDiffRec =0;
14 int F[256]={ 0 }; //直方圖陣列
15 int iTotalGray=0;//灰度值和
16 int iTotalPixel =0;//畫素數和
17 byte bt;//某點的畫素值
18 
19 uchar iThrehold,iNewThrehold;//閥值、新閥值
20 uchar iMaxGrayValue=0,iMinGrayValue=255;//原影象中的最大灰度值和最小灰度值
21 uchar iMeanGrayValue1,iMeanGrayValue2;
22 
23 //獲取(i,j)的值,存於直方圖陣列F
24 for(int i=0;i<width;i++)
25 {
26 for(int j=0;j<height;j++)
27 {
28 bt = data[i*step+j];
29 if(bt<iMinGrayValue)
30 iMinGrayValue = bt;
31 if(bt>iMaxGrayValue)
32 iMaxGrayValue = bt;
33 F[bt]++;
34 }
35 }
36 
37 iThrehold =0;//
38 iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始閥值
39 iDiffRec = iMaxGrayValue - iMinGrayValue;
40 
41 for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中止條件
42 {
43 iThrehold = iNewThrehold;
44 //小於當前閥值部分的平均灰度值
45 for(int i=iMinGrayValue;i<iThrehold;i++)
46 {
47 iTotalGray += F[i]*i;//F[]儲存影象資訊
48 iTotalPixel += F[i];
49 }
50 iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel);
51 //大於當前閥值部分的平均灰度值
52 iTotalPixel =0;
53 iTotalGray =0;
54 for(int j=iThrehold+1;j<iMaxGrayValue;j++)
55 {
56 iTotalGray += F[j]*j;//F[]儲存影象資訊
57 iTotalPixel += F[j]; 
58 }
59 iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel);
60 
61 iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2; //新閥值
62 iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1);
63 }
64 
65 //cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl;
66 return iThrehold;
67 }
68
Otsu程式碼一
1 /*======================================================================*/
2 /* OTSU global thresholding routine */
3 /* takes a 2D unsigned char array pointer, number of rows, and */
4 /* number of cols in the array. returns the value of the threshold */
5 /*parameter: 
6 *image --- buffer for image
7 rows, cols --- size of image
8 x0, y0, dx, dy --- region of vector used for computing threshold
9 vvv --- debug option, is 0, no debug information outputed
10 */
11 /*
12 OTSU 演算法可以說是自適應計算單閾值(用來轉換灰度影象為二值影象)的簡單高效方法。
13 下面的程式碼最早由 Ryan Dibble提供,此後經過多人Joerg.Schulenburg, R.Z.Liu 等修改,補正。
14 演算法對輸入的灰度影象的直方圖進行分析,將直方圖分成兩個部分,使得兩部分之間的距離最大。
15 劃分點就是求得的閾值。
16 */
17 /*======================================================================*/
18 int otsu (unsigned char*image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
19 {
20 
21 unsigned char*np; // 影象指標
22 int thresholdValue=1; // 閾值
23 int ihist[256]; // 影象直方圖,256個點
24 
25 int i, j, k; // various counters
26 int n, n1, n2, gmin, gmax;
27 double m1, m2, sum, csum, fmax, sb;
28 
29 // 對直方圖置零
30 memset(ihist, 0, sizeof(ihist));
31 
32 gmin=255; gmax=0;
33 // 生成直方圖
34 for (i = y0 +1; i < y0 + dy -1; i++) 
35 {
36 np = (unsigned char*)image[i*cols+x0+1];
37 for (j = x0 +1; j < x0 + dx -1; j++)
38 {
39 ihist[*np]++;
40 if(*np > gmax) gmax=*np;
41 if(*np < gmin) gmin=*np;
42 np++; /* next pixel */
43 }
44 }
45 
46 // set up everything
47 sum = csum =0.0;
48 n =0;
49 
50 for (k =0; k <=255; k++) 
51 {
52 sum += (double) k * (double) ihist[k]; /* x*f(x) 質量矩*/
53 n += ihist[k]; /* f(x) 質量 */
54 }
55 
56 if (!n) 
57 {
58 // if n has no value, there is problems...
59 fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
60 return (160);
61 }
62 
63 // do the otsu global thresholding method
64 fmax =-1.0;
65 n1 =0;
66 for (k =0; k <255; k++)
67 {
68 n1 += ihist[k];
69 if (!n1) 
70 { 
71 continue; 
72 }
73 n2 = n - n1;
74 if (n2 ==0)
75 { 
76 break; 
77 }
78 csum += (double) k *ihist[k];
79 m1 = csum / n1;
80 m2 = (sum - csum) / n2;
81 sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
82 /* bbg: note: can be optimized. */
83 if (sb > fmax) 
84 {
85 fmax = sb;
86 thresholdValue = k;
87 }
88 }
89 
90 // at this point we have our thresholding value
91 
92 // debug code to display thresholding values
93 if ( vvv &1 )
94 fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
95 thresholdValue, gmin, gmax);
96 
97 return(thresholdValue);
98 } 
Otsu程式碼二
1 /*======================================================================*/
2 /* OTSU global thresholding routine */
3 /*======================================================================*/
4 int otsu2 (IplImage *image)
5 {
6 int w = image->width;
7 int h = image->height;
8 
9 unsigned char*np; // 影象指標
10 unsigned char pixel;
11 int thresholdValue=1; // 閾值
12 int ihist[256]; // 影象直方圖,256個點
13 
14 int i, j, k; // various counters
15 int n, n1, n2, gmin, gmax;
16 double m1, m2, sum, csum, fmax, sb;
17 
18 // 對直方圖置零...
19 memset(ihist, 0, sizeof(ihist));
20 
21 gmin=255; gmax=0;
22 // 生成直方圖
23 for (i =0; i < h; i++) 
24 {
25 np = (unsigned char*)(image->imageData + image->widthStep*i);
26 for (j =0; j < w; j++) 
27 {
28 pixel = np[j];
29 ihist[ pixel]++;
30 if(pixel > gmax) gmax= pixel;
31 if(pixel < gmin) gmin= pixel;
32 }
33 }
34 
35 // set up everything
36 sum = csum =0.0;
37 n =0;
38 
39 for (k =0; k <=255; k++) 
40 {
41 sum += k * ihist[k]; /* x*f(x) 質量矩*/
42 n += ihist[k]; /* f(x) 質量 */
43 }
44 
45 if (!n) 
46 {
47 // if n has no value, there is problems...
48 //fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
49 thresholdValue =160;
50 goto L;
51 }
52 
53 // do the otsu global thresholding method
54 fmax =-1.0;
55 n1 =0;
56 for (k =0; k <255; k++) 
57 {
58 n1 += ihist[k];
59 if (!n1) { continue; }
60 n2 = n - n1;
61 if (n2 ==0) { break; }
62 csum += k *ihist[k];
63 m1 = csum / n1;
64 m2 = (sum - csum) / n2;
65 sb = n1 * n2 *(m1 - m2) * (m1 - m2);
66 /* bbg: note: can be optimized. */
67 if (sb > fmax)
68 {
69 fmax = sb;
70 thresholdValue = k;
71 }
72 }
73 
74 L:
75 for (i =0; i < h; i++) 
76 {
77 np = (unsigned char*)(image->imageData + image->widthStep*i);
78 for (j =0; j < w; j++) 
79 {
80 if(np[j] >= thresholdValue)
81 np[j] =255;
82 else np[j] =0;
83 }
84 }
85 
86 //cout<<"The Threshold of this Image in Otsu is:"<<thresholdValue<<endl;
87 return(thresholdValue);
88 }
最大熵閥值
1 /*============================================================================
2 = 程式碼內容:最大熵閾值分割 
3 = 修改日期:2009-3-3 
4 = 作者:crond123 
5 = 部落格:http://blog.csdn.net/crond123/
6 = E_Mail:[email protected] 
7 ===============================================================================*/
8 // 計算當前位置的能量熵
9 double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,entropy_state state)
10 {
11 int start,end;
12 int total =0;
13 double cur_entropy =0.0;
14 if(state == back) 
15 {
16 start =0;
17 end = cur_threshold; 
18 }
19 else 
20 {
21 start = cur_threshold;
22 end =256; 
23 } 
24 for(int i=start;i<end;i++) 
25 {
26 total += (int)cvQueryHistValue_1D(Histogram1,i);//查詢直方塊的值 P304
27 }
28 for(int j=start;j<end;j++)
29 {
30 if((int)cvQueryHistValue_1D(Histogram1,j)==0)
31 continue;
32 double percentage = cvQueryHistValue_1D(Histogram1,j)/total;
33 /*熵的定義公式*/
34 cur_entropy +=-percentage*logf(percentage);
35 /*根據泰勒展式去掉高次項得到的熵的近似計算公式
36 cur_entropy += percentage*percentage;*/ 
37 }
38 return cur_entropy;
39 // return (1-cur_entropy);
40 }
41 
42 //尋找最大熵閾值並分割
43 void MaxEntropy(IplImage *src,IplImage *dst)
44 {
45 assert(src != NULL);
46 assert(src->depth ==8&& dst->depth ==8);
47 assert(src->nChannels ==1);
48 CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);//建立一個指定尺寸的直方圖
49 //引數含義:直方圖包含的維數、直方圖維數尺寸的陣列、直方圖的表示格式、方塊範圍陣列、歸一化標誌
50 cvCalcHist(&src,hist);//計算直方圖
51 double maxentropy =-1.0;
52 int max_index =-1;
53 // 迴圈測試每個分割點,尋找到最大的閾值分割點
54 for(int i=0;i<HistogramBins;i++) 
55 {
56 double cur_entropy = caculateCurrentEntropy(hist,i,object)+caculateCurrentEntropy(hist,i,back);
57 if(cur_entropy>maxentropy)
58 {
59 maxentropy = cur_entropy;
60 max_index = i;
61 }
62 }
63 cout<<"The Threshold of this Image in MaxEntropy is:"<<max_index<<endl;
64 cvThreshold(src, dst, (double)max_index,255, CV_THRESH_BINARY);
65 cvReleaseHist(&hist);
66 }
基本全域性閥值法
1 /*============================================================================
2 = 程式碼內容:基本全域性閾值法 
3 ==============================================================================*/
4 int BasicGlobalThreshold(int*pg,int start,int end)
5 { // 基本全域性閾值法
6 int i,t,t1,t2,k1,k2;
7 double u,u1,u2; 
8 t=0; 
9 u=0;
10 for (i=start;i<end;i++) 
11 {
12 t+=pg[i]; 
13 u+=i*pg[i];
14 }
15 k2=(int) (u/t); // 計算此範圍灰度的平均值 
16 do 
17 {
18 k1=k2;
19 t1=0; 
20 u1=0;
21 for (i=start;i<=k1;i++) 
22 { // 計算低灰度組的累加和
23 t1+=pg[i]; 
24 u1+=i*pg[i];
25 }
26 t2=t-t1;
27 u2=u-u1;
28 if (t1) 
29 u1=u1/t1; // 計算低灰度組的平均值
30 else 
31 u1=0;
32 if (t2) 
33 u2=u2/t2; // 計算高灰度組的平均值
34 else 
35 u2=0;
36 k2=(int) ((u1+u2)/2); // 得到新的閾值估計值
37 }
38 while(k1!=k2); // 資料未穩定,繼續
39 //cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<k1<<endl;
40 return(k1); // 返回閾值
41 }
 

---------------------  資料:

https://blog.csdn.net/xw20084898/article/details/17564957