opencv阈值法分割图像

opencv阈值法分割图像

from：

http://www.cnblogs.com/skyseraph/archive/2010/12/21/1913058.html

1。手动设定阈值

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 );

2。自适应法

（１）将图像平均分为若干行列的子图像，行数和列数可以根据应用需要进行选取；　　　　

（２）对每一个子图像，利用全局阈值化方法求最优阈值，所得阈值用图像表示；　　　　

（３）将表示阈值的图像进行必要的延拓后进行滤波；　　　　

（４）将滤波结果作为子图像的新阈值对每个子图像进行阈值化处理

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);

3. 最大熵阀值分割法

http://www.lw23.com/pdf_2db12bdb-1f91-4d63-92d7-30232df7e382/lunwen.pdf

http://read.pudn.com/downloads62/sourcecode/others/214566/%E6%9C%80%E5%A4%A7%E7%86%B5%E5%8E%9F%E7%90%86/%E4%BD%8E%E5%AF%B9%E6%AF%94%E5%BA%A6%E5%9B%BE%E5%83%8F%E4%B8%AD%E6%94%B9%E8%BF%9B%E7%9A%84%E4%BA%8C%E7%BB%B4%E7%86%B5%E9%98%88%E5%80%BC%E5%88%86%E5%89%B2%E6%B3%95.PDF

图像的信息熵反映了图像的总体概貌。若图像中包含目标，则在目标与背景可分割的交界处信息量（即熵）最大。

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 );

///

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 }

//

4. 基本全局阀值法

http://www.lw23.com/pdf_ed5d04b5-a4e3-4d33-b96b-a33c1e7ed707/lunwen.pdf

p(x,y)代表点(x,y)处的灰度值，q(x,y)代表点邻域的某种局部特性。

如果阈值的选取只与p(x,y)有关，则是全局阈值，全局阈值是利用全局信息（例如整个图像的灰度直方图）得到的阈值，它仅与全图各像素的本事性质有关，对全图使用；如果阈值的选取与p(x,y),q(x,y)有关，则是局部阈值。

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);

//

   
   

    
     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 
    
    }
   
   
   
   

    
    ///
   
   
   
   

    
    

   
   
   
   

    
    5.
    
    
    
    

     
     
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);
    
    
   
   
   
   

    
    
    
    

     
     


    
    
   
   
   
   

    
    
    
    

     
     

    
    
   
   
   
   

    
    
    
    

     
     
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 } 
    
    
   
   
   
   

    
    
    
    

     
     

    
    
   
   
   
   

    
    
    
    

     
     


    
    
   
   
   
   

    
    
    
    

     
     


    
    
   
   
   
   

    
    
    
    

     
     
6.
    
    
    
    
    
    

     
     
/*上下阀值法：利用正态分布求可信区间*/
    
    
   
   
   
   

    
    
    
    

     
     


    
    
   
   

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);

/

7. 迭代法

迭代法是基于逼近的思想，其步骤如下： 
1． 求出图象的最大灰度值和最小灰度值，分别记为ZMAX和ZMIN，令初始阈值T0=(ZMAX+ZMIN)/2； 
2． 根据阈值TK将图象分割为前景和背景，分别求出两者的平均灰度值ZO和ZB 
3． 求出新阈值TK+1=(ZO+ZB)/2； 
4． 若TK=TK+1，则所得即为阈值；否则转2，迭代计算。

当TK=TK+1时，说明阈值已经达到稳定了，不需在进行迭代，当TK！=TK+1时，说明阈值还没有稳定，需要用迭代再进一步来确定。

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

///