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

图像算法：图像阈值分割

SkySeraph Dec 21st 2010  HQU

Email：zgzhaobo@gmail.com    QQ：452728574

Latest Modified Date：Dec.21st 2010 HQU

更多精彩请直接访问SkySeraph个人站点：www.skyseraph.com 

 图像算法系列： http://skyseraph.com/2011/08/27/CV/图像算法专题/ 

一、工具：VC+OpenCV

二、语言：C++

三、原理(略)

四、程序

主程序(核心部分) 

/*===============================图像分割=====================================*/
 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);

模块程序

迭代法

/*======================================================================*/
 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代码一 

/*======================================================================*/
 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代码二 

/*======================================================================*/
 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 }

 

最大熵阀值 

/*============================================================================
 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 }

基本全局阀值法 

/*============================================================================
 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 }

五 效果(略)

 

Author:         SKySeraph

Email/GTalk: zgzhaobo@gmail.com    QQ：452728574

From：         http://www.cnblogs.com/skyseraph/

本文版权归作者和博客园共有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文连接，请尊重作者的劳动成果。