【图像算法】七种常见阈值分割代码(Otsu、最大熵、迭代法、自适应阀值、手动、迭代法、基本全局阈值法)

图像算法：图像阈值分割

SkySeraph Dec 21st 2010  HQU

Email：zgzhaobo@gmail.com    QQ：452728574

Latest Modified Date：Dec.21st 2010 HQU

一、工具：VC+OpenCV

二、语言：C++

三、原理(略)

四、程序

主程序(核心部分) 

代码

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: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); // 查询直方块的值 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   }