图像分割--分水岭分割法

一：介绍

            分水岭法将一幅图像看成一个拓扑地形图，灰度值被认为是地形高度值。高灰度值对应着山峰，低灰度值对应着山谷。水从高处流下时，会流到地势低的地方，直到某一局部低洼处才停下来，这个低洼处被称为吸水盆地。最终所有的水会分局在不同的吸水盆地，吸水盆地之间的山脊被称为分水岭，水会从分水岭流下去，朝不同的吸水盆地的流去的可能性是等同的。

        对于图像分割，分水岭法就是在灰度图像中找出不同的吸水盆地和分水岭，不同的吸水盆地代表不同的区域，分水岭代表边缘，分割的主要目的是找到分水线。

         分水岭的计算过程是一个迭代标注过程。分水岭比较经典的计算方法是L. Vincent提出的。在该算法中，分水岭计算分两个步骤，一个是排序过程，一个是淹没过程。首先对每个像素的灰度级进行从低到高排序，然后在从低到高实现淹没过程中，对每一个局部极小值在h阶高度的影响域采用先进先出(FIFO)结构进行判断及标注。

         分水岭变换得到的是输入图像的集水盆图像，集水盆之间的边界点，即为分水岭。显然，分水岭表示的是输入图像极大值点。因此，为得到图像的边缘信息，通常把梯度图像作为输入图像，即

                            g(x,y)=grad(f(x,y))={[f(x,y)-f(x-1,y)]2[f(x,y)-f(x,y-1)]2}0.5

                                                 式中，f(x,y)表示原始图像，grad{.}表示梯度运算。

分水岭算法对微弱边缘具有良好的响应，图像中的噪声、物体表面细微的灰度变化，都会产生过度分割的现象。但同时应当看出，分水岭算法对微弱边缘具有良好的响应，是得到封闭连续边缘的保证的。另外，分水岭算法所得到的封闭的集水盆，为分析图像的区域特征提供了可能。

为消除分水岭算法产生的过度分割，通常可以采用两种处理方法，一是利用先验知识去除无关边缘信息。二是修改梯度函数使得集水盆只响应想要探测的目标。

为降低分水岭算法产生的过度分割，通常要对梯度函数进行修改，一个简单的方法是对梯度图像进行阈值处理，以消除灰度的微小变化产生的过度分割。即

                                          g(x,y)=max(grad(f(x,y)),gθ)

                                          式中，gθ表示阈值。

程序可采用方法：用阈值限制梯度图像以达到消除灰度值的微小变化产生的过度分割，获得适量的区域，再对这些区域的边缘点的灰度级进行从低到高排序，然后在从低到高实现淹没的过程，梯度图像用Sobel算子计算获得。对梯度图像进行阈值处理时，选取合适的阈值对最终分割的图像有很大影响，因此阈值的选取是图像分割效果好坏的一个关键。缺点：实际图像中可能含有微弱的边缘，灰度变化的数值差别不是特别明显，选取阈值过大可能会消去这些微弱边缘。

二：算法步骤

三：实例

#include <iostream>  
#include<opencv2\highgui\highgui.hpp>
#include<opencv2\opencv.hpp>
  
using namespace std;    
using namespace cv;

IplImage* marker_mask = 0;  
IplImage* markers = 0;  
IplImage* img0 = 0, *img = 0, *img_gray = 0, *wshed = 0;  
CvPoint prev_pt = {-1,-1};  

void on_mouse( int event, int x, int y, int flags, void* param );//opencv 会自动给函数传入合适的值  
static void ShowHelpText() ;

int main(int argc, char** argv)
{
	CvMemStorage* storage=cvCreateMemStorage(0);
	CvRNG rng=cvRNG(-1);   //初始化随机数生成器并返回其状态
	ShowHelpText();
	if (argc==2 && (img0=cvLoadImage(argv[1],-1))!=0)
	{
		return -1;
	}
	else
	{
		img0=cvLoadImage("0.jpg",1);
	}
	cvNamedWindow("image",1);
	cvNamedWindow("watershed transform",1);
	img=cvCloneImage(img0);
	img_gray=cvCloneImage(img0);
	wshed=cvCloneImage(img0);
	marker_mask=cvCreateImage(cvGetSize(img),8,1);
	markers=cvCreateImage(cvGetSize(img),IPL_DEPTH_32S,1);
	cvCvtColor(img,marker_mask,CV_BGR2GRAY);//彩色空间转换
	cvCvtColor(marker_mask,img_gray,CV_GRAY2BGR);
	cvZero(marker_mask);
	cvZero(wshed);
	cvShowImage("image",img);
	cvShowImage("watershed transform",wshed);
	cvSetMouseCallback("image",on_mouse,0);
	for (;; )
	{
		int c=waitKey(0);
		if ((char)c==27)
		{
			break;
		}
		if ((char)c=='r')
		{
			cvZero(marker_mask);
			cvCopy(img0,img);
			cvShowImage("image",img);
		}
		if ((char)c=='w' || (char)c==' ')
		{
			CvSeq* contours=0;
			CvMat* color_tab=0;
			int i ,j ,comp_count=0;

			//下面选将标记的图像取得其轮廓, 将每种轮廓用不同的整数表示  
			//不同的整数使用分水岭算法时，就成为不同的种子点  
			//算法本来就是以各个不同的种子点为中心扩张  
			cvClearMemStorage(storage);  
			cvFindContours( marker_mask, storage, &contours, sizeof(CvContour),  CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE );  
			cvZero( markers );  
			for( ; contours != 0; contours = contours->h_next, comp_count++ )  
			{  
				cvDrawContours(markers, contours, cvScalarAll(comp_count+1),  
					cvScalarAll(comp_count+1), -1, -1, 8, cvPoint(0,0) );  
			}  
			//cvShowImage("image",markers);  
			if( comp_count == 0 )  
				continue;  
			color_tab = cvCreateMat( 1, comp_count, CV_8UC3 );//创建随机颜色列表  
			//不同的整数标记
			for( i = 0; i < comp_count; i++ )  
			{  
				uchar* ptr = color_tab->data.ptr + i*3;  
				ptr[0] = (uchar)(cvRandInt(&rng)%180 + 50);  
				ptr[1] = (uchar)(cvRandInt(&rng)%180 + 50);  
				ptr[2] = (uchar)(cvRandInt(&rng)%180 + 50);  
			}  
			{  
				double t = (double)cvGetTickCount();  
				cvWatershed( img0, markers );  
				cvSave("img0.xml",markers);  
				t = (double)cvGetTickCount() - t;  
				printf( "exec time = %gms\n", t/(cvGetTickFrequency()*1000.) );  
			}  
			// paint the watershed image  
			for( i = 0; i < markers->height; i++ )  
				for( j = 0; j < markers->width; j++ )  
				{  
					int idx = CV_IMAGE_ELEM( markers, int, i, j );//markers的数据类型为IPL_DEPTH_32S  
					uchar* dst = &CV_IMAGE_ELEM( wshed, uchar, i, j*3 );//BGR三个通道的数是一起的,故要j*3  
					if( idx == -1 ) //输出时若为-1，表示各个部分的边界  
						dst[0] = dst[1] = dst[2] = (uchar)255;  
					else if( idx <= 0 || idx > comp_count )  //异常情况  
						dst[0] = dst[1] = dst[2] = (uchar)0; // should not get here  
					else //正常情况  
					{  
						uchar* ptr = color_tab->data.ptr + (idx-1)*3;  
						dst[0] = ptr[0]; dst[1] = ptr[1]; dst[2] = ptr[2];  
					}  
				}  
				cvAddWeighted( wshed, 0.5, img_gray, 0.5, 0, wshed );//wshed.x.y=0.5*wshed.x.y+0.5*img_gray+0加权融合图像  
				cvShowImage( "watershed transform", wshed );  
				cvReleaseMat( &color_tab );  
		}  
	}  
	return 1;  
}

static void ShowHelpText()    
{    
    //输出一些帮助信息    
	cout<<endl<<"图像分割--分水岭分割示例程序~"<<endl<<endl;
	cout<<"当前使用的OpenCV版本为 OpenCV "<<CV_VERSION<<endl;
	cout<<"按 'ESC' 退出程序！ "<<endl<<endl;
	cout<<"按 'r' 重新显示原图像！ "<<endl<<endl;
	cout<<"按 'w' 或者 ' ' 运算该算法！ "<<endl<<endl;
	cout<<endl<<"键盘按键任意键- 退出程序"<<endl;
} 

void on_mouse( int event, int x, int y, int flags, void* param )//opencv 会自动给函数传入合适的值  
{  
    if( !img )  
        return;  
    if( event == CV_EVENT_LBUTTONUP || !(flags & CV_EVENT_FLAG_LBUTTON) )  
        prev_pt = cvPoint(-1,-1);  
    else if( event == CV_EVENT_LBUTTONDOWN )  
        prev_pt = cvPoint(x,y);  
    else if( event == CV_EVENT_MOUSEMOVE && (flags & CV_EVENT_FLAG_LBUTTON) )  
    {  
        CvPoint pt = cvPoint(x,y);  
        if( prev_pt.x < 0 )  
            prev_pt = pt;  
        cvLine( marker_mask, prev_pt, pt, cvScalarAll(255), 5, 8, 0 );//CvScalar 成员：double val[4] RGBA值A=alpha  
        cvLine( img, prev_pt, pt, cvScalarAll(255), 5, 8, 0 );  
        prev_pt = pt;  
        cvShowImage( "image", img);  
    }  
}