【图像处理与分析】图像分割方法总结

一、灰度阈值法

1、全局阈值法

整幅图像使用同一个阈值做分割处理，适用于目标物和背景有明显对比的图像。这种方法是只考虑像素值不考虑空间特征的，例如峰谷法、最大熵分割法。

如下是一维最大熵分割法实现。

最大熵原理：只掌握关于未知分布的部分知识时，应选取符合这些知识的但是熵值最大的概率分布。因为熵越大，随机变量的不确定性越大，做准确预测的难度越大。

int maxEntropy(IplImage *src);
float calEntropy(CvHistogram *hist,int begin,int end);

int main()
{
//Mat img =imread("E:/b.jpg",1);

IplImage *src=cvLoadImage("G:\\4.bmp",CV_LOAD_IMAGE_GRAYSCALE);
IplImage *dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
int threshold=maxEntropy(src);
cvThreshold(src, dst, threshold, 255, CV_THRESH_BINARY);
//imshow("一维最大熵分割",img);
cvNamedWindow("一维最大熵分割");
cvShowImage("一维最大熵分割", dst);
cvWaitKey(0);
return 0;
}

int maxEntropy(IplImage *src)
{
assert(src);
assert(src->depth == 8);
assert(src->nChannels == 1);

int sizes=256;
float range[2]={0,255};
float *ranges[1]={&range[0]};
CvHistogram *hist = cvCreateHist(1, &sizes, CV_HIST_ARRAY, ranges, 1);
cvCalcHist(&src,hist,0,0);
cvNormalizeHist(hist,0);

int threthold=0;
int Entropy_max=0;
for(int i=0;i<sizes;i++)
{
float entropy=calEntropy(hist,0,i)+calEntropy(hist,i,sizes);
if(entropy>=Entropy_max)
Entropy_max=entropy;
}
return Entropy_max;
}

float calEntropy(CvHistogram *hist,int begin,int end)
{
float total=0;//总概率
for(int i=begin;i<end;i++)
{
total+=cvQueryHistValue_1D(hist,i);
}
float entropy=0;
for(int i=begin;i<end;i++)
{
float probability=cvQueryHistValue_1D(hist,i);
if(probability==0)
continue;
probability/=total;
entropy-=probability*log(probability);
}

return entropy;
}

2、自适应阈值法

在许多情况下，目标和背景的对比度在图像中的各处是不一样的，实际处理时，需要按照具体情况将图像分成若干子区域分别选择阈值，或者动态地根据一定的邻域范围选择每点处的阈值，进行图像分割。

OTSU法，它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。

以前在遥感影像灰度阈值分割上尝试了一些自适应阈值法，影像背景地物比较复杂，效果都不太理想。

3、人工确定

二、基于边缘的图像分割

1、差分法

图像差分，就是把两幅图像的对应像素值相减，以削弱图像的相似部分，突出显示图像的变化部分。【原理都是共通的，这个也是比较简单的变化检测方法】

2、梯度法

Canny\Sobel算子等，他们是里面包含有计算梯度幅值和方向的步骤，也是一种梯度算子吧。

Canny 步骤：首先消除噪声，一般使用高斯平滑滤波；然后计算梯度幅值和方向；然后非极大值抑制，也就是说选择局部范围内的极大值，抑制局部范围内的非极大值；最后是双阈值算法检测和连接，这一步是为了保留主要边缘，剔除不连续的、次要的边缘。

这里原理解释主要参考：

http://blog.csdn.net/likezhaobin/article/details/6892176

代码根据该博主提供的代码进行了重构

http://blog.csdn.net/likezhaobin/article/details/6892629

/*****
int* M ;                       //梯度幅值
unsigned char* N;            //非极大值抑制结果
double* Theta ;             //梯度方向
******/
void Gaussian(unsigned char* nImageData,unsigned char* pCanny,int nWidth,int nHeight);
void nonMaximumSuppression(unsigned char *pCanny,unsigned char* N,int* M,double* Theta,int nWidth,int nHeight);//极大值抑制
void doubleThresholdDetection(unsigned char* N,int* M,int nWidth,int nHeight);//双阈值检测
void TraceEdge(int y, int x, int nThrLow, unsigned char* pResult, int *pMag, int nWidth,int nHeight);
void main()
{
/******彩色图灰度化********/
IplImage *colorImage=cvLoadImage("G:\\5.jpg",-1);//读图，获取彩色图指针

int width=colorImage->width;
int height=colorImage->height;
int sizeofPixel=width*height;//像素数
IplImage *opencvGrayImage=cvCreateImage(cvSize(width,height),IPL_DEPTH_8U,1);//定义变换后的灰度图指针
char *gray=(char*)malloc(sizeof(char
4000
)*width*height);//存放转换后的灰度像素值
for(int i=0;i<sizeofPixel;i++)
{
uchar b=(uchar)colorImage->imageData[3*i+0];
uchar g = (uchar)colorImage->imageData[3 * i + 1];
uchar r = (uchar)colorImage->imageData[3 * i + 2];
gray[i] = r * 0.299 + g * 0.587 + b * 0.114;
}
opencvGrayImage->imageData=gray;

unsigned char* nImageData = new unsigned char[width*height];   //保存图像中的数据
for(int j=0; j<height; j++)                                     //获取数据
{
for(int i=0; i<width; i++)
nImageData[j*width+i] = (unsigned char)opencvGrayImage->imageData[j*width+i];
}

unsigned char*pCanny = new unsigned char[width*height];        //为平滑后的图像数据分配内存
Gaussian(nImageData,pCanny,width,height);//高斯模糊，不太理解

int* M = new int[width*height];                       //梯度幅值
unsigned char* N =new unsigned char[width*height];  //非极大值抑制结果
double* Theta = new double[width*height];             //梯度方向
nonMaximumSuppression(pCanny,N,M,Theta,width,height);//非极大值抑制结果
doubleThresholdDetection(N,M,width,height);//双阈值检测
opencvGrayImage->imageData=(char*)N;
cvNamedWindow("GrayImage",CV_WINDOW_AUTOSIZE);
cvShowImage("GrayImage",opencvGrayImage);               //显示灰度图
cvWaitKey(0);
cvDestroyWindow("GrayImage");
}

void Gaussian(unsigned char* nImageData,unsigned char* pCanny,int nWidth,int nHeight)
{
/****高斯滤波****/
double nSigma = 0.4;                            //定义高斯函数的标准差
int nWidowSize = 1+2*ceil(3*nSigma);            //定义滤波窗口的大小
int nCenter = (nWidowSize)/2;                   //定义滤波窗口中心的索引

double* nData = new double[nWidth*nHeight];                      //两次平滑的中间数据

//////////////////////生成一维高斯滤波系数//////////////////////////////////
double* pdKernal_2 = new double[nWidowSize*nWidowSize]; //定义一维高斯核数组
double  dSum_2 = 0.0;                                   //求和，进行归一化
///////////////////////二维高斯函数公式////////////////////////////////////
////                         x*x+y*y                        ///////////////
////                   -1*--------------                ///////////////
////         1             2*Sigma*Sigma                ///////////////
////   ---------------- e                                   ///////////////
////   2*pi*Sigma*Sigma                                     ///////////////
///////////////////////////////////////////////////////////////////////////
for(int i=0; i<nWidowSize; i++)
{
for(int j=0; j<nWidowSize; j++)
{
int nDis_x = i-nCenter;
int nDis_y = j-nCenter;
pdKernal_2[i+j*nWidowSize]=exp(-(1/2)*(nDis_x*nDis_x+nDis_y*nDis_y)
/(nSigma*nSigma))/(2*3.1415926*nSigma*nSigma);
dSum_2 += pdKernal_2[i+j*nWidowSize];
}
}
for(int i=0; i<nWidowSize; i++)
{
for(int j=0; j<nWidowSize; j++)                 //进行归一化
{
pdKernal_2[i+j*nWidowSize] /= dSum_2;
}
}
int x;
int y;
for(int i=0; i<nHeight; i++)
{
for(int j=0; j<nWidth; j++)
{
double dFilter=0.0;
double dSum = 0.0;
for(x=(-nCenter); x<=nCenter; x++)                     //行
{
for(y=(-nCenter); y<=nCenter; y++)             //列
{
if( (j+x)>=0 && (j+x)<nWidth && (i+y)>=0 && (i+y)<nHeight) //判断边缘
{
dFilter += (double)nImageData [(i+y)*nWidth + (j+x)]
* pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];
dSum += pdKernal_2[(y+nCenter)*nWidowSize+(x+nCenter)];
}
}
}
pCanny[i*nWidth+j] = (unsigned char)dFilter/dSum;
}
}
}

void nonMaximumSuppression(unsigned char *pCanny,unsigned char* N,int* M,double* Theta,int nWidth,int nHeight)
{
//////////////////同样可以用不同的检测器/////////////////////////
/////    P[i,j]=(S[i,j+1]-S[i,j]+S[i+1,j+1]-S[i+1,j])/2     /////
/////    Q[i,j]=(S[i,j]-S[i+1,j]+S[i,j+1]-S[i+1,j+1])/2     /////
/////////////////////////////////////////////////////////////////
double* P = new double[nWidth*nHeight];                 //x向偏导数
double* Q = new double[nWidth*nHeight];                 //y向偏导数

//计算x,y方向的偏导数
for(int i=0; i<(nHeight-1); i++)
{
for(int j=0; j<(nWidth-1); j++)
{
P[i*nWidth+j] = (double)(pCanny[i*nWidth + min(j+1, nWidth-1)] - pCanny[i*nWidth+j] + pCanny[min(i+1, nHeight-1)*nWidth+min(j+1, nWidth-1)] - pCanny[min(i+1, nHeight-1)*nWidth+j])/2;
Q[i*nWidth+j] = (double)(pCanny[i*nWidth+j] - pCanny[min(i+1, nHeight-1)*nWidth+j] + pCanny[i*nWidth+min(j+1, nWidth-1)] - pCanny[min(i+1, nHeight-1)*nWidth+min(j+1, nWidth-1)])/2;
}
}
//计算梯度幅值和梯度的方向
for(int i=0; i<nHeight; i++)
{
for(int j=0; j<nWidth; j++)
{
M[i*nWidth+j] = (int)(sqrt(P[i*nWidth+j]*P[i*nWidth+j] + Q[i*nWidth+j]*Q[i*nWidth+j])+0.5);
Theta[i*nWidth+j] = atan2(Q[i*nWidth+j], P[i*nWidth+j]) * 57.3;
if(Theta[i*nWidth+j] < 0)
Theta[i*nWidth+j] += 360;              //将这个角度转换到0~360范围
}
}

//非极大值抑制

int g1=0, g2=0, g3=0, g4=0;                            //用于进行插值，得到亚像素点坐标值
double dTmp1=0.0, dTmp2=0.0;                           //保存两个亚像素点插值得到的灰度数据
double dWeight=0.0;                                    //插值的权重

//对边界进行初始化
for(int i=0; i<nWidth; i++)
{
N[i] = 0;
N[(nHeight-1)*nWidth+i] = 0;
}
for(int j=0; j<nHeight; j++)
{
N[j*nWidth] = 0;
N[j*nWidth+(nWidth-1)] = 0;
}
//局部最大值查找
for(int i=1; i<(nWidth-1); i++)
{
for(int j=1; j<(nHeight-1); j++)
{
int nPointIdx = i+j*nWidth;       //当前点在图像数组中的索引值
if(M[nPointIdx] == 0)
N[nPointIdx] = 0;         //如果当前梯度幅值为0，则不是局部最大对该点赋为0
else
{
////////首先判断属于那种情况，然后根据情况插值///////
////////////////////第一种情况///////////////////////
/////////       g1  g2                  /////////////
/////////           C                   /////////////
/////////           g3  g4              /////////////
/////////////////////////////////////////////////////
if( ((Theta[nPointIdx]>=90)&&(Theta[nPointIdx]<135)) ||
((Theta[nPointIdx]>=270)&&(Theta[nPointIdx]<315)))
{
//////根据斜率和四个中间值进行插值求解
g1 = M[nPointIdx-nWidth-1];
g2 = M[nPointIdx-nWidth];
g3 = M[nPointIdx+nWidth];
g4 = M[nPointIdx+nWidth+1];
dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]);   //反正切
dTmp1 = g1*dWeight+g2*(1-dWeight);
dTmp2 = g4*dWeight+g3*(1-dWeight);
}
////////////////////第二种情况///////////////////////
/////////       g1                      /////////////
/////////       g2  C   g3              /////////////
/////////               g4              /////////////
/////////////////////////////////////////////////////
else if( ((Theta[nPointIdx]>=135)&&(Theta[nPointIdx]<180)) ||
((Theta[nPointIdx]>=315)&&(Theta[nPointIdx]<360)))
{
g1 = M[nPointIdx-nWidth-1];
g2 = M[nPointIdx-1];
g3 = M[nPointIdx+1];
g4 = M[nPointIdx+nWidth+1];
dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]);   //正切
dTmp1 = g2*dWeight+g1*(1-dWeight);
dTmp2 = g4*dWeight+g3*(1-dWeight);
}
////////////////////第三种情况///////////////////////
/////////           g1  g2              /////////////
/////////           C                   /////////////
/////////       g4  g3                  /////////////
/////////////////////////////////////////////////////
else if( ((Theta[nPointIdx]>=45)&&(Theta[nPointIdx]<90)) ||
((Theta[nPointIdx]>=225)&&(Theta[nPointIdx]<270)))
{
g1 = M[nPointIdx-nWidth];
g2 = M[nPointIdx-nWidth+1];
g3 = M[nPointIdx+nWidth];
g4 = M[nPointIdx+nWidth-1];
dWeight = fabs(P[nPointIdx])/fabs(Q[nPointIdx]);   //反正切
dTmp1 = g2*dWeight+g1*(1-dWeight);
dTmp2 = g3*dWeight+g4*(1-dWeight);
}
////////////////////第四种情况///////////////////////
/////////               g1              /////////////
/////////       g4  C   g2              /////////////
/////////       g3                      /////////////
/////////////////////////////////////////////////////
else if( ((Theta[nPointIdx]>=0)&&(Theta[nPointIdx]<45)) ||
((Theta[nPointIdx]>=180)&&(Theta[nPointIdx]<225)))
{
g1 = M[nPointIdx-nWidth+1];
g2 = M[nPointIdx+1];
g3 = M[nPointIdx+nWidth-1];
g4 = M[nPointIdx-1];
dWeight = fabs(Q[nPointIdx])/fabs(P[nPointIdx]);   //正切
dTmp1 = g1*dWeight+g2*(1-dWeight);
dTmp2 = g3*dWeight+g4*(1-dWeight);
}
}
//////////进行局部最大值判断，并写入检测结果////////////////
if((M[nPointIdx]>=dTmp1) && (M[nPointIdx]>=dTmp2))
N[nPointIdx] = 128;
else
N[nPointIdx] = 0;
}
}
}

void doubleThresholdDetection(unsigned char* N,int* M,int nWidth,int nHeight)
{
//双阈值检测,最大的梯度幅值为360
int nHist[1024];
int nEdgeNum;             //可能边界数
int nMaxMag = 0;          //最大梯度数
int nHighCount;
//统计直方图
for(int i=0;i<1024;i++)
nHist[i] = 0;
for(int i=0; i<nHeight; i++)
{
for(int j=0; j<nWidth; j++)
{
if(N[i*nWidth+j]==128)
nHist[M[i*nWidth+j]]++;
}
}
//获取最大幅值及潜在边缘个数
nEdgeNum = nHist[0];
nMaxMag = 0;                    //获取最大的梯度值
for(int i=1; i<1024; i++)           //统计经过“非最大值抑制”后有多少像素
{
if(nHist[i] != 0)       //梯度为0的点是不可能为边界点的
{
nMaxMag = i;
}
nEdgeNum += nHist[i];   //经过non-maximum suppression后有多少像素
}
//这段代码的意思是，按照灰度值从低到高的顺序，选取前79%个灰度值中的最大的灰度值为高阈值，低阈值大约为高阈值的一半,这是根据经验数据的来的.
double  dRatHigh = 0.79;
double  dThrHigh;
double  dThrLow;
double  dRatLow = 0.5;
nHighCount = (int)(dRatHigh * nEdgeNum + 0.5);
int j=1;
nEdgeNum = nHist[1];
while((j<(nMaxMag-1)) && (nEdgeNum < nHighCount))
{
j++;
nEdgeNum += nHist[j];
}
dThrHigh = j;                                   //高阈值
dThrLow = (int)((dThrHigh) * dRatLow + 0.5);    //低阈值
for(int i=0; i<nHeight; i++)
{
for(int j=0; j<nWidth; j++)
{
if((N[i*nWidth+j]==128) && (M[i*nWidth+j] >= dThrHigh))
{
N[i*nWidth+j] = 255;
TraceEdge(i, j, dThrLow, N, M, nWidth,nHeight);
}
}
}
}

void TraceEdge(int y, int x, int nThrLow, unsigned char* pResult, int *pMag, int nWidth,int nHeight)
{
//对8邻域像素进行查询
int xNum[8] = {1,1,0,-1,-1,-1,0,1};
int yNum[8] = {0,1,1,1,0,-1,-1,-1};
long yy,xx,k;
for(k=0;k<8;k++)
{
yy = y+yNum[k];
xx = x+xNum[k];
if(pResult[yy*nWidth+xx]==128 && pMag[yy*nWidth+xx]>=nThrLow )
{
//该点设为边界点
pResult[yy*nWidth+xx] = 255;
//以该点为中心再进行跟踪
TraceEdge(yy,xx,nThrLow,pResult,pMag,nWidth,nHeight);
}
}
}

3、高通滤波

图像的边缘、细节主要是集中在高频部分，为了突出边缘，则采用高通滤波器让高频通过，而削弱低频部分，对条带状噪声的影像做去噪工作也是有奇效，但是，高通滤波容易产生振铃现象或者伪现象。

常用滤波器有二维理想高通滤波器、巴特沃斯高通滤波器、指数高通滤波器、梯形高通滤波器。

步骤：

1、首先对一副图像进行如下二维傅里叶变换。



然后F[0,0]是直流分量



直流分量表示图像较为平滑没有太大变化的部分,其他为交流分量,所以理想的情况应该是保持直流分量而对其他部分进行增幅，但是滤波器是会衰减直流分量的，所以会产生振铃效应，导致图像盲恢复会产生图像信息的损失。

2、然后选择高通滤波器，如

二维理想高通滤波器：



D0是滤波器的阻带半径，而D(u,v)是点到滤波器中央的距离。理想高通的滤波器的振幅特性如下所示



巴特沃斯高通滤波器：



巴特沃斯高通滤波器也可以通过改变次数n，对过度特性进行调整。过大的n会造成振铃现象。



高斯高通滤波器



高斯滤波器的过度特性很好，所以不会发生振铃现象。



参考资料：

https://www.cnblogs.com/fuhaots2009/p/3465149.html

4、模板匹配

模板匹配与图像分割交集比较少，目前看的基本是基于形状模型的图像分割，模板匹配对其有一定帮助，局限性挺大的。我同学的毕设老师说模板匹配太低级了，不给用哈哈哈哈

模板匹配是一种最原始、最基本的模式识别方法，研究某一特定对象物的图案位于图像的什么地方，进而识别对象物，这就是一个匹配问题。它是图像处理中最基本、最常用的匹配方法。模板匹配具有自身的局限性，主要表现在它只能进行平行移动，若原图像中的匹配目标发生旋转或大小变化，该算法无效。

在 OpenCV 中，提供了相应的函数完成这个操作。

matchTemplate 函数：在模板和输入图像之间寻找匹配,获得匹配结果图像

minMaxLoc 函数：在给定的矩阵中寻找最大和最小值，并给出它们的位置

其中OpenCV 提供了 6 种计算两幅图像相似度的方法。

差值平方和匹配 CV_TM_SQDIFF

标准化差值平方和匹配 CV_TM_SQDIFF_NORMED

相关匹配 CV_TM_CCORR

标准相关匹配 CV_TM_CCORR_NORMED

相关匹配 CV_TM_CCOEFF

标准相关匹配 CV_TM_CCOEFF_NORMED

具体方法公式可以参考以下链接。

参考资料http://blog.csdn.net/zhi_neng_zhi_fu/article/details/51029864

http://blog.csdn.net/liyuanbhu/article/details/49837661

三、基于区域的图像分割

1、区域生长法

首先对每个需要分割的区域找出一个种子像素作为生长的七点，然后将种子像素周围邻域中与种子有相同或相似性质的像素（根据事先确定的生长或相似准则来确定）合并到种子像素所在的区域中。而新的像素继续作为种子向四周生长，直到再没有满足条件的像素可以包括进来，一个区域就生长而成了。

这个方法有两个关键：一个是种子点怎么选，一个是阈值的限定。原理代码以前po过。

http://blog.csdn.net/weixin_37051000/article/details/60142281

2、分水岭算法

这个怕是一时半会讲不清，地学领域很爱用。

这里又发展出了多种面向对象的分割方法

四、分割后处理

1、数学形态学：膨胀、腐蚀、开运算、闭运算

2、图像平滑、锐化