图像分割技术（1）

图像分割综述：

针对单色的图像分割通常基于处理灰度值的两类特性：不连续性和相似性。

第一类中，方法是以灰度突变为基础分割一幅图像，比如图像的边缘。

第二类中，主要方法是根据一组预定义的准则将一幅图像分割为相似的区域，如与之处理、区域生长、区域分裂和区域聚合。

     分割是将图像细分为子区域，这些子区域互不重叠，并集为初始图像，每个子区域内的像素分布符合预定规则。用数学描述可表示为：

一般有4种分割思路：基于点线、边缘的分割；基于阈值的分割（二值化）；基于区域的分割；基于分水岭的分割

1. 基于点线、边缘的分割

这种分割的目的是 找出来，这也是大部分机器视觉应用的必要步骤。

无论是点、线还是边缘的分割，都是基于图像微分（《图像局部不变性特征与描述》阅读笔记（3）-- 点与边缘检测）

由(LOG)Laplacian of Guassian & (DOH)Determinant of Hessian 斑点检测对高斯一阶和二阶微分的讨论，可以知道：

1）一阶导数得到的边缘比二阶要粗

2）二阶导数对细线、噪声更敏感

3）二阶导数在灰度过渡处会产生双边缘，此外通过二阶导数的符号可以推断边缘的过渡是从暗到亮还是相反

基于以上讨论，孤立点的检测一般用二阶导数（laplace），如Harris角点检测，除此之外还有很多别的角点检测算法，如SUSAN，SIFT，SURF

线的检测主要针对空间上连续的线段，一般使用对方向敏感的模板对线做方向筛选

图1-1. 各向互异的Laplace线检测算子

边缘检测比点线检测要复杂，因为边缘一般是弯曲的，因此需要首先检测点或小线段，然后做拟合，其大致步骤是：

1）平滑图像以过滤噪声

2）使用微分获得边缘点

3）基于一定规则对边缘点进行筛选并连接

简单的方法一般有Roberts、Prewitt、Sobel等线检测算子，之后做点膨胀连接下断点也就可以了

更有效的方法是Canny算法，效果很好

然后再补充下边缘连接的方法吧，虽然觉得canny已经足够了，而这些方法往往会增加很多计算量：

1）局部处理方式：对每个边缘点A，考虑其邻域的点B，如果两个点无论梯度幅值还是梯度角度都在误差范围内，那么认为B在A代表的边缘上

2）多边形拟合：找到一个多边形，尽可能的反映当前离散点的分布

3）全局处理方式：比较有针对性的场合，比如需要找直线，找圆，用Hough

2. 基于阈值的分割

常用的是单阈值分割，这样得到的图像就是二值的，当然也有多阈值分割以区分多种区域。

而基于阈值的作用域又可分为全局阈值分割和局部阈值分割。

阈值分割的主要参考是图像直方图，阈值选择在波谷以区分不同区域

2.1 全局阈值分割

主要介绍几种主流的全局阈值分割方法（单阈值，多阈值在这个基础上做二次迭代就好了）：最大熵、基本全局、迭代法、OTSU（以上4种思路及代码主要参考SkySeraph的文章《图像算法：图像阈值分割》，手动选择阈值的方法就不说了）

a) 最大熵阈值

依次选择每个灰度级（0~255），将直方图分成两部分，分别计算两部分的熵再相加，取使得熵为最大值的那个灰度级为阈值。

令Ai为直方图A部分第i级的值，At为A部分所有灰度级的和，那么A部分的熵为：

[cpp]  view plain  copy

1. // 计算当前位置的能量熵  
2.  double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,int state)  
3.  {  
4.     int start,end;  
5.     int total =0;  
6.     double cur_entropy =0.0;  
7.     if(state == 1)           // 1是计算当前点之前的能量熵  
8.     {  
9.         start =0;  
10.         end = cur_threshold;   
11.     }  
12.     else   
13.     {  
14.         start = cur_threshold;  
15.         end =256;   
16.     }   
17.     for(int i=start;i<end;i++)   
18.          total += (int)cvQueryHistValue_1D(Histogram1,i);//查询直方块的值 P304  
19.     for(int j=start;j<end;j++)  
20.     {  
21.         if((int)cvQueryHistValue_1D(Histogram1,j)==0)  
22.         continue;  
23.         double percentage = cvQueryHistValue_1D(Histogram1,j)/total;  
24.         /*熵的定义公式*/  
25.         cur_entropy +=-percentage*logf(percentage);  
26.     }  
27.     return cur_entropy;  
28.   
29.  }  
30.    
31.  //寻找最大熵阈值并分割  
32.  int MaxEntropy(IplImage *src)  
33.  {  
34.     int HistogramBins=256;  
35.     CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,0,1);//创建一个指定尺寸的直方图  
36.     //参数含义：直方图包含的维数、直方图维数尺寸的数组、直方图的表示格式、方块范围数组、归一化标志  
37.     cvCalcHist(&src,hist);//计算直方图  
38.     double maxentropy =-1.0;  
39.     int max_index =-1;  
40.     // 循环测试每个分割点，寻找到最大的阈值分割点  
41.     for(int i=0;i<HistogramBins;i++)   
42.     {  
43.         double cur_entropy = caculateCurrentEntropy(hist,i,0)+caculateCurrentEntropy(hist,i,1);  
44.         if(cur_entropy>maxentropy)  
45.         {  
46.             maxentropy = cur_entropy;  
47.             max_index = i;  
48.         }  
49.     }  
50.     cvReleaseHist(&hist);  
51.     return max_index;  
52.  }  

b) 基本全局阈值

基本思路是首先设定一个阈值，然后通过计算质量矩调整阈值，直到阈值稳定

[cpp]  view plain  copy

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. }  

c) 迭代法

基本思路和全局阈值法差不多，也是基于一个初始阈值不停迭代直到稳定，只是衡量准则有点不同

[cpp]  view plain  copy

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. }  

d) OTSU法

这应该是阈值分割里最有名的方法了.OTSU(大津算法)将类间方差作为区域分割的衡量标准，认为方差越大分割越精确

[cpp]  view plain  copy

1. /* 
2.  OTSU 算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。 
3.  下面的代码最早由 Ryan Dibble提供，此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改，补正。 
4.  算法对输入的灰度图像的直方图进行分析，将直方图分成两个部分，使得两部分之间的距离最大。 
5.  划分点就是求得的阈值。 
6.  */  
7.  int otsu (IplImage *image)  
8.  {  
9.      int thresholdValue=1;              // 阈值                    
10.      int i,k;                           // various counters  
11.      int n, n1, n2;  
12.      <span style="white-space:pre">   </span> double m1, m2, sum, csum, fmax, sb;  
13.      // 生成直方图  
14.      int HistogramBins=256;  
15.      CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,0,1);//创建一个指定尺寸的直方图  
16.      cvCalcHist(&image,hist);//计算直方图  
17.      // set up everything  
18.      sum = csum =0.0;  
19.      n =0;  
20.      for (k =0; k <256; k++)   
21.      {  
22.         i=cvQueryHistValue_1D(hist,k);  
23.         sum += k * i; /* x*f(x) 质量矩*/  
24.         n += i; /* f(x) 质量 */  
25.      }  
26.     if (!n)   
27.     {  
28.      // if n has no value, there is problems...  
29.      thresholdValue =160;  
30.     }  
31.     <span style="white-space:pre">    </span>// do the otsu global thresholding method  
32.     fmax =-1.0;  
33.     n1 =0;  
34.     for (k =0; k <255; k++)   
35.     {  
36.         i=cvQueryHistValue_1D(hist,k);  
37.         n1 += i;  
38.         if (!n1) { continue; }  
39.         n2 = n - n1;  
40.         if (n2 ==0) { break; }  
41.         csum += k *i;  
42.         m1 = csum / n1;  
43.         m2 = (sum - csum) / n2;  
44.         sb = n1 * n2 *(m1 - m2) * (m1 - m2);  
45.         /* bbg: note: can be optimized. */  
46.         if (sb > fmax)  
47.         {  
48.             fmax = sb;  
49.             thresholdValue = k;  
50.         }  
51.     }  
52.     cvReleaseHist(&hist);  
53.     return(thresholdValue);  
54.  }