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

一、工具：VC+OpenCV

二、語言：C++

三、原理

     

      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演算法之前，需要計算該映像的灰階長條圖。

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

映像閾值分割---迭代演算法

   1 .處理流程：

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

       適合映像長條圖有明顯波穀

四、程式

主程式(核心部分) 

閾值分割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_C14  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)的值，存於長條圖數組F24 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 }68Otsu代碼一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 image7 rows, cols --- size of image8 x0, y0, dx, dy --- region of vector used for computing threshold9 vvv --- debug option, is 0, no debug information outputed10 */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 counters26 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 everything47 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 method64 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 value91 92 // debug code to display thresholding values93 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 counters15 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 everything36 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 method54 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:crond123@163.com 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);//查詢直方塊的值 P30427 }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 }