图像分割技术（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部分的熵为：

// 计算当前位置的能量熵 double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,int state) {int start,end;int total =0;double cur_entropy =0.0;if(state == 1)           // 1是计算当前点之前的能量熵{start =0;end = cur_threshold; }else {start = cur_threshold;end =256; } for(int i=start;i<end;i++)  total += (int)cvQueryHistValue_1D(Histogram1,i);//查询直方块的值 P304for(int j=start;j<end;j++){if((int)cvQueryHistValue_1D(Histogram1,j)==0)continue;double percentage = cvQueryHistValue_1D(Histogram1,j)/total;/*熵的定义公式*/cur_entropy +=-percentage*logf(percentage);}return cur_entropy; }  //寻找最大熵阈值并分割 int MaxEntropy(IplImage *src) {int HistogramBins=256;CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,0,1);//创建一个指定尺寸的直方图//参数含义：直方图包含的维数、直方图维数尺寸的数组、直方图的表示格式、方块范围数组、归一化标志cvCalcHist(&src,hist);//计算直方图double maxentropy =-1.0;int max_index =-1;// 循环测试每个分割点，寻找到最大的阈值分割点for(int i=0;i<HistogramBins;i++) {double cur_entropy = caculateCurrentEntropy(hist,i,0)+caculateCurrentEntropy(hist,i,1);if(cur_entropy>maxentropy){maxentropy = cur_entropy;max_index = i;}}cvReleaseHist(&hist);return max_index; }

b) 基本全局阈值

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

/*=============================================================================  代码内容：基本全局阈值法                              ==============================================================================*/int BasicGlobalThreshold(int*pg,int start,int end){                                               int  i,t,t1,t2,k1,k2;    double u,u1,u2;        t=0;         u=0;    for (i=start;i<end;i++)     {        t+=pg[i];                u+=i*pg[i];    }    k2=(int) (u/t);                          //  计算此范围灰度的平均值        do     {        k1=k2;        t1=0;            u1=0;for (i=start;i<=k1;i++)         {            //  计算低灰度组的累加和            t1+=pg[i];                u1+=i*pg[i];        }        t2=t-t1;        u2=u-u1;if (t1)             u1=u1/t1;                     //  计算低灰度组的平均值else             u1=0;if (t2)             u2=u2/t2;                     //  计算高灰度组的平均值else             u2=0;        k2=(int) ((u1+u2)/2);                 //  得到新的阈值估计值    }    while(k1!=k2);                           //  数据未稳定，继续   //cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<k1<<endl;   return(k1);                              //  返回阈值}

c) 迭代法

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

/*======================================================================*//* 迭代法*//*======================================================================*/// nMaxIter：最大迭代次数；nDiffRec：使用给定阀值确定的亮区与暗区平均灰度差异值int DetectThreshold(IplImage*img, int nMaxIter, int& iDiffRec) {     //图像信息    int height = img->height;    int width = img->width;    int step = img->widthStep/sizeof(uchar);    uchar *data = (uchar*)img->imageData;    iDiffRec =0;    int F[256]={ 0 }; //直方图数组    int iTotalGray=0;//灰度值和    int iTotalPixel =0;//像素数和    byte bt;//某点的像素值    uchar iThrehold,iNewThrehold;//阀值、新阀值    uchar iMaxGrayValue=0,iMinGrayValue=255;//原图像中的最大灰度值和最小灰度值    uchar iMeanGrayValue1,iMeanGrayValue2;    //获取(i,j)的值，存于直方图数组F    for(int i=0;i<width;i++)    {for(int j=0;j<height;j++)        {            bt = data[i*step+j];            if(bt<iMinGrayValue)                iMinGrayValue = bt;    if(bt>iMaxGrayValue)                iMaxGrayValue = bt;            F[bt]++;        }    }    iThrehold =0;//    iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始阀值    iDiffRec = iMaxGrayValue - iMinGrayValue;    for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中止条件    {iThrehold = iNewThrehold;//小于当前阀值部分的平均灰度值for(int i=iMinGrayValue;i<iThrehold;i++){iTotalGray += F[i]*i;//F[]存储图像信息iTotalPixel += F[i];}iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel);//大于当前阀值部分的平均灰度值iTotalPixel =0;iTotalGray =0;for(int j=iThrehold+1;j<iMaxGrayValue;j++){iTotalGray += F[j]*j;//F[]存储图像信息iTotalPixel += F[j];    }iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel);iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2;        //新阀值iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1);     }//cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl;return iThrehold;}

d) OTSU法

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

/* OTSU 算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。 下面的代码最早由 Ryan Dibble提供，此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改，补正。 算法对输入的灰度图像的直方图进行分析，将直方图分成两个部分，使得两部分之间的距离最大。 划分点就是求得的阈值。 */ int otsu (IplImage *image) { int thresholdValue=1;              // 阈值                   int i,k;                           // various counters int n, n1, n2;     <span style="white-space:pre"></span> double m1, m2, sum, csum, fmax, sb; // 生成直方图 int HistogramBins=256; CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,0,1);//创建一个指定尺寸的直方图 cvCalcHist(&image,hist);//计算直方图 // set up everything sum = csum =0.0; n =0; for (k =0; k <256; k++)  {i=cvQueryHistValue_1D(hist,k);sum += k * i; /* x*f(x) 质量矩*/n += i; /* f(x) 质量 */ }if (!n) { // if n has no value, there is problems... thresholdValue =160;}    <span style="white-space:pre"></span>// do the otsu global thresholding methodfmax =-1.0;n1 =0;for (k =0; k <255; k++) {i=cvQueryHistValue_1D(hist,k);n1 += i;if (!n1) { continue; }n2 = n - n1;if (n2 ==0) { break; }csum += k *i;m1 = csum / n1;m2 = (sum - csum) / n2;sb = n1 * n2 *(m1 - m2) * (m1 - m2);/* bbg: note: can be optimized. */if (sb > fmax){fmax = sb;thresholdValue = k;}}cvReleaseHist(&hist);return(thresholdValue); }

下篇图像分割技术（2）