七种常见阈值分割代码(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  }