openCv学习笔记（九）-自己实现分水岭算法并且用openCv测试

   为了验证上一篇博文中的算法，自己稍作修改，即利用openCv里的图像结构。为了简单，直接采用Mat.at<uchar>(i,j)进行图像操作，效率肯定低，如果感兴趣可参考http://blog.csdn.net/caiye917015406/article/details/7791815改进算法。不过这里只是测试效果，就将就吧。。。

  在算法中要手动选取阀值进行图像的二值化，这给分水岭算法的效果有很大影响。也曾采用openCv里的cvAdaptiveThreshold自动选取 寻求,但效果还不如不好，期待改进。。。

 以下是测试代码

#include<iostream>
#include <queue>
#include <vector>

#include<cv.h>
#include<highgui.h>
using namespace std;
using namespace cv;
void Watershed_1(const Mat OriginalImage,Mat SeedImage,Mat LabelImage);
/*====================================================================
函数名：            Watershed
功能：              用标记-分水岭算法对输入图像进行分割
算法实现：          无
输入参数说明：      OriginalImage --输入图像（灰度图，0~255）
                    SeedImage       --标记图像（二值图，0-非标记，1-标记）
                    LabelImage      --输出图像（1-第一个分割区域，2-第二个分割区域，...）
返回值说明：        无      
====================================================================*/
int main()
{
	//声明变量
	IplImage *OriginalImage,*SeedImage,*LabelImage,*LabelImage2;
	OriginalImage=cvLoadImage("D:\\openCV\\openCVProject\\img\\stuff.jpg",CV_LOAD_IMAGE_GRAYSCALE);
	//OriginalImage_gray =cvCreateImage(cvGetSize(OriginalImage), 8,1); 
//	cvtColor( OriginalImage, OriginalImage, CV_BGR2GRAY );//转化为灰度图
	cvNamedWindow("sorc");
	cvNamedWindow("seed");
	cvNamedWindow("result1");
	cvNamedWindow("result2");
    //创建操作图像
	LabelImage = cvCreateImage(cvGetSize(OriginalImage), 8,1);
	SeedImage = cvCreateImage(cvGetSize(OriginalImage), 8,1);
	LabelImage2= cvCreateImage(cvGetSize(OriginalImage), 8,3);
    //图像二值化，要手动取阀值
	cvThreshold(OriginalImage,SeedImage,195,1,CV_THRESH_BINARY);
    //cvAdaptiveThreshold(OriginalImage,SeedImage,100);
	
	cvShowImage("sorc",OriginalImage);
	cvShowImage("seed",SeedImage);
    
	Watershed_1( OriginalImage,SeedImage, LabelImage);
    cvShowImage("result1",LabelImage);
    //找出处理后的轮廓
	CvMemStorage* storage = cvCreateMemStorage();
	CvSeq* first_contour =NULL;
	int Nc = cvFindContours(LabelImage,storage,&first_contour,sizeof(CvContour),CV_RETR_CCOMP);
	int n=0;
	printf("Toatal Contours Detected: %d\n",Nc);
	cvZero(LabelImage2);
	//画出处理后的轮廓
	for(CvSeq* c=first_contour;c!=NULL;c=c->h_next)
	{
		//cvCvtColor(img_8uc1,img_8uc3,CV_GRAY2BGR);
		cvDrawContours(LabelImage2,c,CV_RGB(0,200,0), CV_RGB(255, 0, 0),2,2,8);
		//printf("Contour #%d\n",n);
		//cvShowImage("result",LabelImage);
		//printf("%d elements:\n",c->total);
		//for(int i=0;i<c->total;++i)
		//{
			//CvPoint* p = CV_GET_SEQ_ELEM(CvPoint,c,i);
			//printf("(%d,%d)\n",p->x,p->y );
		
		//}
		//cvWaitKey(0);
		n++;
	}

	cvShowImage("result2",LabelImage2);
	//cvShowImage("result2",markers);
	cvWaitKey(0);
	cvReleaseImage(&OriginalImage);
	cvDestroyWindow("sorc");
	cvReleaseImage(&SeedImage);
	cvDestroyWindow("seed");
	cvReleaseImage(&LabelImage);
	cvDestroyWindow("result1");
	cvReleaseImage(&LabelImage2);
	cvDestroyWindow("result2");
	return 0;
}
void Watershed_1(const Mat OriginalImage,Mat SeedImage,Mat LabelImage)
{
	int row,col;
	row = OriginalImage.rows;
	col = OriginalImage.cols;
	//标记区域标识号，从1开始
	int Num=0;
	int i,j;

	//保存每个队列种子个数的数组
	vector<int*> SeedCounts;
	//临时种子队列
	queue<POINT> quetem;
	//保存所有标记区域种子队列的数组，里面放的是种子队列的指针
	vector<queue<POINT>*> vque;

	int* array;
	//指向种子队列的指针.at<uchar>(i,j)
	queue<POINT> *pque;
	POINT temp;

	for(i=0;i<row;i++)
	{
		for(j=0;j<col;j++)
			LabelImage.at<uchar>(i,j)=0;
	}


	int m,n,k=0;
	BOOL up,down,right,left,upleft,upright,downleft,downright;//8 directions...

	//预处理,提取区分每个标记区域，并初始化每个标记的种子队列
	//种子是指标记区域边缘的点，他们可以在水位上升时向外淹没（或者说生长）
	//pan's words：我的理解是梯度值较小的象素点，或者是极小灰度值的点。
	for(i=0;i<row;i++)
	{
		for(j=0;j<col;j++)
		{
			//如果找到一个标记区域
			if(SeedImage.at<uchar>(i,j)==1)
			{
				//区域的标识号加一
				Num++;
				//分配数组并初始化为零，表示可有256个灰阶
				array=new int[256];
				ZeroMemory(array,256*sizeof(int));
				//种子个数数组进vector，每次扫描则生成一个数组，并用区域标识号来做第一维。灰度级做第二维。
				//表示某个盆地区域中某灰阶所对应的点的数目。
				SeedCounts.push_back(array);
				//分配本标记号的优先队列，256个种子队列，
				//表示对应一个灰阶有一个队列，并且每个队列可以存储一个集合的点信息
				pque=new queue<POINT>[256];
				//加入到队列数组中，对应的是本标记号Num的
				vque.push_back(pque);
				//当前点放入本标记区域的临时种子队列中
				temp.x=i;
				temp.y=j;
				quetem.push(temp);
				//当前点标记为已处理
				LabelImage.at<uchar>(i,j)=Num;
				SeedImage.at<uchar>(i,j)=127;//表示已经处理过

				//让临时种子队列中的种子进行生长直到所有的种子都生长完毕
				//生长完毕后的队列信息保存在vque中，包括区域号和灰阶，对应点数存储在seedcounts中
				while(!quetem.empty())
				{
					up=down=right=left=FALSE;
					upleft=upright=downleft=downright=FALSE;
					//队列中取出一个种子
					temp=quetem.front();
					m=temp.x;
					n=temp.y;
					quetem.pop();
					//注意到127对扫描过程的影响，影响下面的比较，但是不影响while语句中的扫描

					if(m>0)
					{
						//上方若为可生长点则加为新种子
						if(SeedImage.at<uchar>(m-1,n)==1)
						{
							temp.x=m-1;
							temp.y=n;
							quetem.push(temp);//如果这样的话，那么这些标记过的区域将再次在while循环中被扫描到，不会，因为值是127
							//新种子点标记为已淹没区域，而且是当前区域，并记录区域号到labelImage
							LabelImage.at<uchar>(m-1,n)=Num;
							SeedImage.at<uchar>(m-1,n)=127;
						}
						else//否则上方为不可生长
						{
							up=TRUE;
						}
					}
					if(m>0&&n>0)
					{
						if(SeedImage.at<uchar>(m-1,n-1)==1)//左上方若为可生长点则加为新种子
						{
							temp.x=m-1;
							temp.y=n-1;
							quetem.push(temp);
							//新种子点标记为已淹没区域，即下一个循环中以127来标识不再扫描，而且是当前区域
							LabelImage.at<uchar>(m-1,n-1)=Num;
							SeedImage.at<uchar>(m-1,n-1)=127;
						}
						else//否则左上方为不可生长
						{
							upleft=TRUE;
						}
					}

					if(m<row-1)
					{
						if(SeedImage.at<uchar>(m+1,n)==1)//下方若为可生长点则加为新种子
						{
							temp.x=m+1;
							temp.y=n;
							quetem.push(temp);
							//新种子点标记为已淹没区域，而且是当前区域
							LabelImage.at<uchar>(m+1,n)=Num;
							SeedImage.at<uchar>(m+1,n)=127;
						}
						else//否则下方为不可生长
						{
							down=TRUE;
						}
					}
					if(m<(row-1)&&n<(col-1))
					{
						if(SeedImage.at<uchar>(m+1,n+1)==1)//下方若为可生长点则加为新种子
						{
							temp.x=m+1;
							temp.y=n+1;
							quetem.push(temp);
							//新种子点标记为已淹没区域，而且是当前区域
							LabelImage.at<uchar>(m+1,n+1)=Num;
							SeedImage.at<uchar>(m+1,n+1)=127;
						}
						else//否则下方为不可生长
						{
							downright=TRUE;
						}
					}

					if(n<col-1)
					{
						if(SeedImage.at<uchar>(m,n+1)==1)//右方若为可生长点则加为新种子
						{
							temp.x=m;
							temp.y=n+1;
							quetem.push(temp);
							//新种子点标记为已淹没区域，而且是当前区域
							LabelImage.at<uchar>(m,n+1)=Num;
							SeedImage.at<uchar>(m,n+1)=127;
						}
						else//否则右方为不可生长
						{
							right=TRUE;
						}
					}
					if(m>0&&n<(col-1))
					{
						if(SeedImage.at<uchar>(m-1,n+1)==1)//右上方若为可生长点则加为新种子
						{
							temp.x=m-1;
							temp.y=n+1;
							quetem.push(temp);
							//新种子点标记为已淹没区域，而且是当前区域
							LabelImage.at<uchar>(m-1,n+1)=Num;
							SeedImage.at<uchar>(m-1,n+1)=127;
						}
						else//否则右上方为不可生长
						{
							upright=TRUE;
						}
					}

					if(n>0)
					{
						if(SeedImage.at<uchar>(m,n-1)==1)//左方若为可生长点则加为新种子
						{
							temp.x=m;
							temp.y=n-1;
							quetem.push(temp);
							//新种子点标记为已淹没区域
							LabelImage.at<uchar>(m,n-1)=Num;
							SeedImage.at<uchar>(m,n-1)=127;
						}
						else//否则左方为不可生长
						{
							left=TRUE;
						}
					}
					if(m<(row-1)&&n>0)
					{
						if(SeedImage.at<uchar>(m+1,n-1)==1)//左下方若为可生长点则加为新种子
						{
							temp.x=m+1;
							temp.y=n-1;
							quetem.push(temp);
							//新种子点标记为已淹没区域
							LabelImage.at<uchar>(m+1,n-1)=Num;
							SeedImage.at<uchar>(m+1,n-1)=127;
						}
						else//否则左方为不可生长
						{
							downleft=TRUE;
						}
					}

					//上下左右只要有一点不可生长，那么本点为初始种子队列中的一个
					//这里可否生长是由seedimage中的值来决定的。
					if(up||down||right||left||
						upleft||downleft||upright||downright)
					{
						temp.x=m;
						temp.y=n;
						//下面这个矢量数组：第一维是标记号；第二维是该图像点的灰度级
						//m,n点对应的是while循环里面扫描的像素点。
						//Num是当前的区域号
						//这样这个二维信息就表示了，某个区域中对应某个灰度级对应的成员点的集合与个数
						//分别由下面两个量来表达
						vque[Num-1][OriginalImage.at<uchar>(m,n)].push(temp);//这两句中我把Num-1改成了Num...pan's codes...
						SeedCounts[Num-1][OriginalImage.at<uchar>(m,n)]++;
					}     
				}//while结束，扫描到quetem为空而止。也就是对应所有的节点都得到不可生长为止（或者是周围的点要么不可生长，要么已生长）
			}//if结束
		}
	}
	//在上述过程中，如果标记的点为0则表示，没有扫描到的点，或者表明不是输入的种子点
	//这里相当于是找seedimage传过来的初始区域的分水岭界线的所有的点；并且用标号记录每个区域，同时集水盆的边缘点进入队列。
	//上面是找集水盆的程序。同时也是连通区域。
漫水//
	bool actives;//在某一水位处，所有标记的种子生长完的标志
	int WaterLevel;

	//淹没过程开始，水位从零开始上升，水位对应灰度级，采用四连通法
	for(WaterLevel=0;WaterLevel<180;WaterLevel++)//第二维。。。
	{
		actives=true;
		while(actives)
		{
			actives=false;
			//依次处理每个标记号所对应的区域，且这个标记号对应的区域的点的个数在SeedCounts里面
			for(i=0;i<Num;i++)//第一维。。。
			{
				if(!vque[i][WaterLevel].empty())//对应的分水岭不为空集,i表示区域号，waterlevel表示灰阶
				{
					actives=true;
					while(SeedCounts[i][WaterLevel]>0)
					{
						SeedCounts[i][WaterLevel]--;//取出一个点，个数少一
						temp=vque[i][WaterLevel].front();//取出该区域的一个点，清空这个边缘点，表示当前
						//灰度级该像素已经处理掉了。
						vque[i][WaterLevel].pop();
						m = temp.x;
						n = temp.y;//当前种子的坐标
						if(m>0)
						{
							if(!LabelImage.at<uchar>(m-1,n))//上方若未处理，表示没有标号，应该在输入前已经作过初始化为0
								//本函数中在开头也作过初始化
							{
								temp.x=m-1;
								temp.y=n;
								LabelImage.at<uchar>(m-1,n)=i+1;//上方点标记为已淹没区域
								//注意到这个标记是与扫描点的区域号相同，一定在这个标号所属的区域吗？是的
								//这样在下一轮至少会扫描到这个点，确保不遗漏，但是下一轮的处理会使它合理
								//归类吗？问题还有这样标记并没有一定将它加入到种子队列。也就是说它
								//只是被淹没而不能向上淹没。只有满足下述可生长条件才行。
								if(OriginalImage.at<uchar>(m-1,n)<=WaterLevel)//上方若为可生长点则加入当前队列，当前高度的队列
								{
									vque[i][WaterLevel].push(temp);
								}
								else//否则加入OriginalImage[m-1][n]对应的灰度级的队列，为什么？
								{
									vque[i][OriginalImage.at<uchar>(m-1,n)].push(temp);
									SeedCounts[i][OriginalImage.at<uchar>(m-1,n)]++;
								}
							}
						}

						if(m<row-1)
						{
							if(!LabelImage.at<uchar>(m+1,n))//下方若未处理
							{
								temp.x=m+1;
								temp.y=n;
								LabelImage.at<uchar>(m+1,n)=i+1;//下方点标记为已淹没区域

								if(OriginalImage.at<uchar>(m+1,n)<=WaterLevel)//下方若为可生长点则加入当前队列
								{
									vque[i][WaterLevel].push(temp);
								}
								else//否则加入OriginalImage[m+1][n]级队列
								{
									vque[i][OriginalImage.at<uchar>(m+1,n)].push(temp);
									SeedCounts[i][OriginalImage.at<uchar>(m+1,n)]++;
								}
							}
						}

						if(n<col-1)
						{
							if(!LabelImage.at<uchar>(m,n+1))//右边若未处理
							{
								temp.x=m;
								temp.y=n+1;
								LabelImage.at<uchar>(m,n+1)=i+1;//右边点标记为已淹没区域

								if(OriginalImage.at<uchar>(m,n+1)<=WaterLevel)//右边若为可生长点则加入当前队列
								{
									vque[i][WaterLevel].push(temp);
								}
								else//否则加入OriginalImage[m][n+1]级队列
								{
									vque[i][OriginalImage.at<uchar>(m,n+1)].push(temp);
									SeedCounts[i][OriginalImage.at<uchar>(m,n+1)]++;
								}
							}
						}

						if(n>0)
						{
							if(!LabelImage.at<uchar>(m,n-1))//左边若未处理
							{
								temp.x=m;
								temp.y=n-1;
								LabelImage.at<uchar>(m,n-1)=i+1;//左边点标记为已淹没区域

								if(OriginalImage.at<uchar>(m,n-1)<=WaterLevel)//左边若为可生长点则加入当前队列
								{
									vque[i][WaterLevel].push(temp);
								}
								else//否则加入OriginalImage[m][n-1]级队列
								{
									vque[i][OriginalImage.at<uchar>(m,n-1)].push(temp);
									SeedCounts[i][OriginalImage.at<uchar>(m,n-1)]++;
								}
							}
						}
					}//while循环结束
					SeedCounts[i][WaterLevel]=vque[i][WaterLevel].size();
				}//if结束
			}//for循环结束
		}//while循环结束
	}//for循环结束
	while(!vque.empty())
	{
		pque=vque.back();
		delete[] pque;
		vque.pop_back();
	}
	while(!SeedCounts.empty())
	{
		array=SeedCounts.back();
		delete[] array;
		SeedCounts.pop_back();
	}

}