图像分割算法：概述及常用边缘检测算法

图像分割算法概述
用计算机进行数字图像处理的目的有两个，一是产生更适合人类视觉观察和识别的图像，二是希望计算机能够自动进行识别和理解图像。无论是为了何种目的，图像处理的关键一步是对包含有大量各式各样景物信息的图像进行分解。分解的最终结果就是图像被分成一些具有各种特征的最小成分，这些成分就称为图像的基元。产生这些基元的过程就是图像分割的过程。图像分割作为图像处理领域中极为重要的内容之一，是实现图像分析与理解的基础。从概念上来说，所谓图像分割就是按照一定的原则将一幅图像或景物分为若干个部分或子集的过程。目前图像处理系统中我们只能得到二维图像信息，因此只能进行图像分割而不是景物分割（景物是三维信息）；图像分割也可以理解为将图像中有意义的特征区域或者需要应用的特征区域提取出来，这些特征区域可以是像素的灰度值、物体轮廓曲线、纹理特性等，也可以是空间频谱或直方图特征等。在图像中用来表示某一物体的区域，其特征都是相近或相同的，但是不同物体的区域之间，特征就会急剧变化。目前已经提出的图像分割方法很多，从分割依据的角度来看，图像的分割方法可以分为相似性分割和非连续性分割。相似性分割就是将具有同一灰度级或相同组织结构的像素聚集在一起，形成图像的不同区域；非连续性分割就是首先检测局部不连续性，然后将它们连接在一起形成边界，这些边界将图像分成不同的区域。由于不同种类的图像，不同的应用场合，需要提取的图像特征是不同的，当然对应的图像特征提取方法也就不同，因此并不存在一种普遍适应的最优方法。

　　图像分割方法又可分为结构分割方法和非结构分割方法两大类。结构分割方法是根据图像的局部区域象素的特征来实现图像分割，如阈值分割、区域生长、边缘检测、纹理分析等，这些方法假定事先知道这些区域的特性，或者在处理过程中能够求得这些特性，从而能够寻找各种形态或研究各像素群。非结构分割法包括统计模式识别、神经网络方法或其它利用景物的先验知识实现的方法等等。总之，图像分割可以分为图像的边缘提取和图像的二值化二部分内容，下面我们首先来讨论一下各种常用的图像边缘提取的方法。

 图像边缘检测
　　数字图像的边缘检测是图像分割、目标区域的识别、区域形状提取等图像分析领域十分重要的基础，是图像识别中提取图像特征的一个重要属性，图像理解和分析的第一步往往就是边缘检测，目前它以成为机器视觉研究领域最活跃的课题之一，在工程应用中占有十分重要的地位。物体的边缘是以图像的局部特征不连续的形式出现的，也就是指图像局部亮度变化最显著的部分，例如灰度值的突变、颜色的突变、纹理结构的突变等，同时物体的边缘也是不同区域的分界处。图像边缘有方向和幅度两个特性，通常沿边缘的走向灰度变化平缓，垂直于边缘走向的像素灰度变换剧烈，根据灰度变化的特点，可分为阶跃型、房顶型和凸缘型，

　　由于边缘是图像上灰度变化最剧烈的地方，传统的边缘检测就是利用了这个特点，对图像各个像素点进行微分或求二阶微分来确定边缘像素点。一阶微分图像的峰值处对应着图像的边缘点；二阶微分图像的过零点处对应着图像的边缘点。根据数字图像的特点，处理图像过程中常采用差分来代替导数运算，对于图像的简单一阶导数运算，图像边缘提取的常用梯度算子有Robert算子、Sobel算子、Prewitt算子、Krisch算子等。

下面以边缘检测Sobel算子为例来讲述数字图像处理中边缘检测的实现：对于数字图像，可以用一阶差分代替一阶微分；

　　　   　△xf(x,y)=f(x,y)-f(x-1,y);

              　△yf(x,y)=f(x,y)-f(x,y-1);

　　求梯度时对于平方和运算及开方运算，可以用两个分量的绝对值之和表示，Sobel梯度算子是先做成加权平均，再微分，然后求梯度，即：　

              △xf(x,y)= f(x-1,y+1) + 2f(x,y+1) + f(x+1,y+1)-f(x-1,y-1) - 2f(x,y-1) - f(x+1,y-1);

　　　　△yf(x,y)= f(x-1,y-1) + 2f(x-1,y) +f(x-1,y+1)- f(x+1,y-1) - 2f(x+1,y) - f(x+1,y+1);

　　　　G[f(x,y)]=|△xf(x,y)|+|△yf(x,y)|;

　　

常用的检测实现公式列出如下（检测模版可以从相应的算法很容易的得到）：

　　Roberts算子：G[i,i]=|f[i,j]-f[i+1,j+1]|+|f[i+1,j]-f[i,j+1]|;

　　Sobe算子：G[i,i]=|f[i-1,j+1]+2f[i,j+1]+f[i+1,j+1]-f[i-1,j-1]-2f[i,j-1]-f[i+1,j-1]|

　　　　　　　　　　　+|f[i-1,j-1]+2f[i-1,j]+f[i-1,j+1]-f[i+1,j-1]-2f[i+1,j]-f[i+1,j+1]|;

　　其中G[i,j]表示处理后(i,j)点的灰度值，f[i,j]表示处理前该点的灰度值。

在视图类中定义了响应菜单命令的边缘检测Sobel算子实现灰度图像边缘检测的函数：

void CDibView::OnMENUSobel()

{

　CClientDCpDC(this);

HDChDC=pDC.GetSafeHdc();//获取当前设备上下文的句柄；

SetStretchBltMode(hDC,COLORONCOLOR);

HANDLEdata1handle;

　LPDIBHDRTMAPINFOHEADERlpDIBHdr;

　CDibDoc*pDoc=GetDocument();

　HDIBhdib;

　unsignedchar *lpDIBBits;

　unsignedchar *data;

　hdib=pDoc->m_hDIB;//得到图象数据；

　lpDIBHdr=(LPDIBHDRTMAPINFOHEADER)GlobalLock((HGLOBAL)hdib);

　lpDIBBits=lpDIBHdr +* (LPDWORD)lpDIBHdr + 256*sizeof(RGBQUAD);

　//得到指向位图像素值的指针;

data1handle=GlobalAlloc(GMEM_SHARE,WIDTHBYTES(lpDIBHdr->biWidth*8)*lpDIBHdr->biHeight);　//申请存放处理后的像素值的缓冲区

data=(unsignedchar*)GlobalLock((HGLOBAL)data1handle);

　AfxGetApp()->BeginWaitCursor();

　inti,j,buf,buf1,buf2;

　for(j=0; jbiHeight; j++)//以下循环求(x,y)位置的灰度值

　　for(i=0; ibiWidth; i++)

　{

　　　if(((i-1)>=0)&&((i+1)biWidth)&&((j-1)>=0)&&((j+1)biHeight))

　　　{//对于图像四周边界处的向素点不处理

　　　buf1=(int)*(lpDIBBits+(i+1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j-1))

　　　　　+2*(int)*(lpDIBBits+(i+1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j))

　　　　　+(int)(int)*(lpDIBBits+(i+1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j+1));

　　　　buf1=buf1-(int)(int)*(lpDIBBits+(i-1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j-1))

　　　　　-2*(int)(int)*(lpDIBBits+(i-1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j))

　　　　　-(int)(int)*(lpDIBBits+(i-1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j+1));

　　　　//以上是对图像进行水平（x）方向的加权微分

　　　　buf2=(int)(int)*(lpDIBBits+(i-1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j+1))

　　　　　+2*(int)(int)*(lpDIBBits+(i)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j+1))

　　　　　+(int)(int)*(lpDIBBits+(i+1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j+1));　

　　　　buf2=buf2-(int)(int)*(lpDIBBits+(i-1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j-1))

　　　　　-2*(int)(int)*(lpDIBBits+(i)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j-1))

　　　　　-(int)(int)*(lpDIBBits+(i+1)*WIDTHBYTES(lpDIBHdr->biWidth*8)+(j-1));

　

　　　　//以上是对图像进行垂直(y)方向加权微分

　　　　buf=abs(buf1)+abs(buf2);//求梯度

　　　if(buf>255) buf=255;

　　　　　if(buf<0)buf=0;

　　　*(data+i*WIDTHBYTES(lpDIBHdr->biWidth*8)+j)=(BYTE)buf;

　　　}

　　　else *(data+i*lpDIBHdr->biWidth+j)=(BYTE)0;

　　}

　for(j=0; jbiHeight; j++)

　for(i=0; ibiWidth; i++)

　*(lpDIBBits+i*WIDTHBYTES(lpDIBHdr->biWidth*8)+j)=*(data+i*WIDTHBYTES(lpDIBHdr->biWidth*8)+j);　　　　　 //处理后的数据写回原缓冲区

StretchDIBits(hDC,0,0,lpDIBHdr->biWidth,lpDIBHdr->biHeight,0,0,

lpDIBHdr->biWidth,lpDIBHdr->biHeight,

lpDIBBits,(LPDIBHDRTMAPINFO)lpDIBHdr,

DIB_RGB_COLORS,

SRCCOPY);　

} 

对于边缘检测，只要知道有若干个检测模板（既边缘检测矩阵），就可以直接实现检测功能了，常用的检测实现公式列出如下（检测模版可以从相应的算法很容易的得到）：

　　Roberts算子：G[i,i]=|f[i,j]-f[i+1,j+1]|+|f[i+1,j]-f[i,j+1]|;

　　Sobe算子：G[i,i]=|f[i-1,j+1]+2f[i,j+1]+f[i+1,j+1]-f[i-1,j-1]-2f[i,j-1]-f[i+1,j-1]|+|f[i-1,j-1]+2f[i-1,j]+f[i-1,j+1]-f[i+1,j-1]-2f[i+1,j]-f[i+1,j+1]|;

　　其中G[i,j]表示处理后(i,j)点的灰度值，f[i,j]表示处理前该点的灰度值。

Kirsch算子实现起来相对来说稍微麻烦一些，它采用8个模板对图像上的每一个像素点进行卷积求导数，这8个模板代表8个方向，对图像上的8个特定边缘方向作出最大响应，运算中取最大值作为图像的边缘输出（上述算法中用到的8个模板在下面的实现代码中给出）。为了便于读者理解该算法的实现，这里我们给出实现该算法的函数代码，读者可以稍加改动应用到自己的项目中去。

BOOL Kirsch(BYTE *pData,int Width,int Height)

{//定义实现Kirsch算法的8个模板；

inti,j,s,t,k,max,sum[8];

statica[3][3]={{+5,+5,+5},{-3,0,-3},{-3,-3,-3}};

statica1[3][3]={{-3,+5,+5},{-3,0,+5},{-3,-3,-3}};

statica2[3][3]={{-3,-3,+5},{-3,0,+5},{-3,-3,+5}};

statica3[3][3]={{-3,-3,-3},{-3,0,+5},{-3,+5,+5}};

statica4[3][3]={{-3,-3,-3},{-3,0,-3},{+5,+5,+5}};

statica5[3][3]={{-3,-3,-3},{+5,0,-3},{+5,+5,-3}};

statica6[3][3]={{+5,-3,-3},{+5,0,-3},{+5,-3,-3}};

statica7[3][3]={{+5,+5,-3},{+5,0,-3},{-3,-3,-3}};

BYTE*pData1;

if(pData==NULL)

{

AfxMessageBox("图像数据为空!");

returnFALSE;

}

pData1=(BYTE*)newchar[Width*Height];

if(pData1==NULL)

{

AfxMessageBox("图像缓冲数据区申请失败!");

returnFALSE ;

}

memcpy(pData1,pData,Width*8*Height);

//kirsch算子处理，对每一像素点求取八个方向的导数；；

for(i=1;i<Height-1;i++)

for(j=1;j<Width-1;j++)

{ 

sum[1]=sum[2]=sum[3]=sum[4]=sum[5]=sum[6]=sum[7]=sum[8]=0;

for(t=-1;t<2;t++)

{

for(s=-1;s<2;s++)

{sum[1]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a[1+t][1+s];

sum[2]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a1[1+t][1+s];sum[3]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a2[1+t][1+s];sum[4]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a3[1+t][1+s];sum[5]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a4[1+t][1+s];sum[6]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a5[1+t][1+s];sum[7]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a6[1+t][1+s];sum[8]+=*(pData+WIDTHBYTES(Width*8)*(i+t)+j+s)*a7[1+t][1+s];

}

}

//取最大方向的导数；

for(k=0;k<8;k++)

{ 

max=0; 

if(max<sum[k])

max=sum[k];

}

if(max<0)

max=0;

if(max>255)

max=255;

*(pData1+Width*8*i+j)=max;

}

memcpy(pData,pData1,Width*8*Height);

deletepData1;

returnTRUE;

} 

　　另外还有一种称为拉普拉斯的算子是不依赖于边缘方向的二阶微分算子，对于数字图像来说拉普拉斯算子可以简单表示为：G[I,j]=|f[i+1,j]+f[i-1,j]+f(i,j+1)+f[i,j-1]-4f[i,j]|;它是一个标量而不是向量，具有旋转不变,既各向同性的性质，它常常用在图像处理的过程中。梯度算子和拉普拉斯算子对噪声敏感，它们都使噪声成份加强，因此在处理含有较大噪声的图像时，常常先对图像进行平滑操作，然后再进行二阶微分，这就产生了所谓的LOG(又称为Marr方法)边缘检测方法。它先用高斯函数对图像进行平滑，然后再用拉普拉斯算子进行运算。总的来说，传统的边缘检测算子的噪声平滑能力和边缘定位能力是矛盾的，为了克服这个不足，正确地得到图像的边缘信息，人们提出了很多方法，如多尺度空间滤波、Facet模型检测边缘、模板匹配、Hough变换、小波变换、人工神经网络、模糊推理等算法。但这些方法绝大多数没有经典的算法精简，要么难以获得合理的计算复杂度，要么需要人为的调节各种参数，有时甚至难以实时运行。因为传统边缘的定义为图像中灰度的突变，所以这样定义边缘既失去了边缘的部分信息，又把噪声的影响包含在了边缘中。其实，边缘往往具有以下特征：

1）灰度突变；

2）是不同区域的边界；

3）具有方向性；

　　根据边缘的这三个特征，可以判断所关心的区域其特征是否存在差异来判断是否存在边缘的可能性。如果特征没有差异，则认为是平滑区；如果特征有差异，则判断为边缘点。

参考文献：

1.    Visual C++数字图像识别技术典型案例

      2.   冈萨雷斯 数字图像处理