使用C++、opencv进行分水岭分割图像

分水岭概念是以对图像进行三维可视化处理为基础的：其中两个是坐标，另一个是灰度级。基于“地形学”的这种解释，我们考虑三类点：

a.属于局部性最小值的点，也可能存在一个最小值面，该平面内的都是最小值点

b.当一滴水放在某点的位置上的时候，水一定会下落到一个单一的最小值点

c.当水处在某个点的位置上时，水会等概率地流向不止一个这样的最小值点

对一个特定的区域最小值，满足条件(b)的点的集合称为这个最小值的“汇水盆地”或“分水岭”。满足条件(c)的点的集合组成地形表面的峰线，称做“分割线”或“分水线”。

分水岭分割方法，是一种基于拓扑理论的数学形态学的分割方法，目前较著名且使用较多的有2种算法:

(1) 自下而上的模拟泛洪的算法 (2) 自上而下的模拟降水的算法

这里介绍泛洪算法的过程。

算法主要思想：

我们把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，模拟泛洪算法的基本思想是：假设在每个区域最小值的位置上打一个洞并且让水以均匀的上升速率从洞中涌出，从低到高淹没整个地形。当处在不同的汇聚盆地中的水将要聚合在一起时，修建的大坝将阻止聚合。水将达到在水线上只能见到各个水坝的顶部这样一个程度。这些大坝的边界对应于分水岭的分割线。所以，它们是由分水岭算法提取出来的(连续的)边界线。

原图像： 地形俯视图：

?

原图像显示了一个简单的灰度级图像，其中“山峰”的高度与输入图像的灰度级值成比例。为了阻止上升的水从这些结构的边缘溢出，我们想像将整幅地形图的周围用比最高山峰还高的大坝包围起来。最高山峰的值是由输入图像灰度级具有的最大值决定的。

?

?

?

图一被水淹没的第一个阶段，这里水用浅灰色表示，覆盖了对应于图中深色背景的区域。在图二和三中，我们看到水分别在第一和第二汇水盆地中上升。由于水持续上升，最终水将从一个汇水盆地中溢出到另一个之中。

?

?

左图中显示了溢出的第一个征兆。这里，水确实从左边的盆地溢出到右边的盆地，并且两者之间有一个短“坝”(由单像素构成)阻止这一水位的水聚合在一起。随着水位不断上升，如右图所显示的那样。这幅图中在两个汇水盆地之间显示了一条更长的坝，另一条水坝在右上角。这条水坝阻止了盆地中的水和对应于背景的水的聚合。

这个过程不断延续直到到达水位的最大值(对应于图像中灰度级的最大值)。水坝最后剩下的部分对应于分水线，这条线就是要得到的分割结果。

?

对于这个例子，分水线在图中显示为叠加到原图上的一个像素宽的深色路径。注意一条重要的性质就是分水线组成一条连通的路径，由此给出了区域之间的连续的边界。

动图演示了整个分水岭算法的过程：

?

算法实现：

?

?

算法应用：

分水岭算法对噪声等影响非常敏感。所以在真实图像中，由于噪声点或者其它干扰因素的存在，使用分水岭算法常常存在过度分割的现象，这是因为很多很小的局部极值点的存在，比如下面的图像，这样的分割效果是毫无用处的。

?

?

为了解决过度分割的问题，可以使用基于标记(mark)图像的分水岭算法，就是通过先验知识，来指导分水岭算法，以便获得更好的图像分段效果。通常的mark图像，都是在某个区域定义了一些灰度层级，在这个区域的洪水淹没过程中，水平面都是从定义的高度开始的，这样可以避免一些很小的噪声极值区域的分割。下面的动图很好的演示了基于mark的分水岭算法过程：

?

上面的过度分割图像，我们通过指定mark区域，可以得到很好的分段效果：

?

?

以上参考：冈萨雷斯《数字图象处理（第三版）》和https://www.cnblogs.com/mikewolf2002/p/3304118.html

相关API：

void setMousecallback(const string& winname, MouseCallback onMouse, void* userdata=0)

winname:窗口的名字
onMouse:鼠标响应函数，回调函数。指定窗口里每次鼠标时间发生的时候，被调用的函数指针。 这个函数的原型应该为void on_Mouse(int event, int x, int y, int flags, void* param);
userdate:传给回调函数的参数

void on_Mouse(int event, int x, int y, int flags, void* param)

event： CV_EVENT_*变量之一
x和y：鼠标指针在图像坐标系的坐标（不是窗口坐标系）
flags：CV_EVENT_FLAG的组合， param是用户定义的传递到setMouseCallback函数调用的参数。

附常用的event：CV_EVENT_MOUSEMOVE、CV_EVENT_LBUTTONDOWN 、CV_EVENT_RBUTTONDOWN、 CV_EVENT_LBUTTONUP 、 CV_EVENT_RBUTTONUP

和标志位flags有关的：CV_EVENT_FLAG_LBUTTON

C++: void watershed(InputArray image,InputoutputArray markers)

第一个参数，InputArray类型的src，输入图像，即源图像，填Mat类的对象即可，且需为8位三通道的彩色图像。

第二个参数，InputOutput Array类型的markers，函数调用后的运算结果存在这里，输入/输出32位单通道图像的标记结果。即这个参数用于存放函数调后的输出结果，需和源图片有一样的尺寸和类型。

代码实现：

#include "stdafx.h"
#include "opencv2/imgproc/imgproc.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include <iostream>  
#include <fstream>  

using namespace cv;
using namespace std;

#define WINDOW_NAME1 "【程序窗口1】"        //为窗口标题定义的宏   
#define WINDOW_NAME2 "【分水岭算法效果图】"        //为窗口标题定义的宏  

//描述：全局变量的声明  
Mat g_maskImage, g_srcImage;
Point prevPt(-1, -1);
//描述：全局函数的声明  
static void ShowHelpText();
static void on_Mouse(int event, int x, int y, int flags, void*);

int main()
{
	//【0】改变console字体颜色  
	system("color 02");

	//【1】载入原图并显示，初始化掩膜和灰度图
	g_srcImage = imread("D:\\pic-sam\\哀.JPG", 1);
	namedWindow(WINDOW_NAME1, WINDOW_NORMAL);
	imshow(WINDOW_NAME1, g_srcImage);
	Mat srcImage, grayImage;
	g_srcImage.copyTo(srcImage);
	cvtColor(g_srcImage, g_maskImage, COLOR_BGR2GRAY);
	cvtColor(g_maskImage, grayImage, COLOR_GRAY2BGR);
	g_maskImage = Scalar::all(0);
	//【2】设置鼠标回调函数
	setMouseCallback(WINDOW_NAME1, on_Mouse, 0);

	//【3】轮询按键，进行处理
	while (1)
	{
		//获取键值
		int c = waitKey(0);

		//若按键键值为ESC时，退出
		if ((char)c == 27)
			break;

		//按键键值为2时，恢复源图
		if ((char)c == '2')
		{
			g_maskImage = Scalar::all(0);
			srcImage.copyTo(g_srcImage);
			imshow("image", g_srcImage);
		}

		//若检测到按键值为1或者空格，则进行处理
		if ((char)c == '1' || (char)c == ' ')
		{
			//定义一些参数
			int i, j, compCount = 0;
			vector<vector<Point> > contours;
			vector<Vec4i> hierarchy;

			//寻找轮廓
			findContours(g_maskImage, contours, hierarchy, RETR_CCOMP, CHAIN_APPROX_SIMPLE);

			//轮廓为空时的处理
			if (contours.empty())
				continue;

			//拷贝掩膜
			Mat maskImage(g_maskImage.size(), CV_32S);
			maskImage = Scalar::all(0);

			//循环绘制出轮廓
			for (int index = 0; index >= 0; index = hierarchy[index][0], compCount++)
				drawContours(maskImage, contours, index, Scalar::all(compCount + 1), -1, 8, hierarchy, INT_MAX);

			//compCount为零时的处理
			if (compCount == 0)
				continue;

			//生成随机颜色
			/*vector<Vec3b> colorTab;
			for (i = 0; i < compCount; i++)
			{
				int b = theRNG().uniform(0, 255);
				int g = theRNG().uniform(0, 255);
				int r = theRNG().uniform(0, 255);

				colorTab.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));
			}*/

			//计算处理时间并输出到窗口中
			double dTime = (double)getTickCount();
			watershed(srcImage, maskImage);
			dTime = (double)getTickCount() - dTime;
			printf("\t处理时间 = %gms\n", dTime*1000. / getTickFrequency());

			//双层循环，将分水岭图像遍历存入watershedImage中
			Mat watershedImage(maskImage.size(), CV_8UC3);
			int index1 = 0;
			for (i = 0; i < maskImage.rows; i++)
				for (j = 0; j < maskImage.cols; j++)
				{
					if(maskImage.at<int>(i, j)>index1)
					index1 = maskImage.at<int>(i, j);
				}
			for (i = 0; i < maskImage.rows; i++)
				for (j = 0; j < maskImage.cols; j++)
				{
					int index = maskImage.at<int>(i, j);
					//对watershed函数生成的index的规律不是很清楚，经测试，并不是按照标记顺序给出index的
					//具体每一块的index是怎么给出的还需要研究源码
					if (index == -1)
						watershedImage.at<Vec3b>(i, j) = Vec3b(255, 255, 255);
					else if (index <= 0 || index > compCount)
						watershedImage.at<Vec3b>(i, j) = Vec3b(0, 0, 0);
					else if (index ==index1)
						watershedImage.at<Vec3b>(i, j) = Vec3b(255, 255, 255);
					else
						watershedImage.at<Vec3b>(i, j) = Vec3b(index*10, 0, 0);//这里想给不同的物体标记为不同程度的颜色
																				//方便后面去除背景，显示目标物体
				}

			//混合灰度图和分水岭效果图并显示最终的窗口
			//watershedImage = watershedImage*0.5 + grayImage*0.5;
			imshow(WINDOW_NAME2, watershedImage);//直接显示分水岭的效果图
			//这里想直接根据index，将背景显示为黑色，需要分割出来的目标物体直接显示
			//但对index生成的规律还未搞清楚，结果可能不是很稳定
			Mat src = imread("D:\\pic-sam\\哀.JPG", 1);
			for (int i = 0; i < src.rows; i++)
				for (int j = 0; j < src.cols; j++)
				{
					int a = abs(watershedImage.at<Vec3b>(i, j)[0] - 250) / 150;
					src.at<Vec3b>(i, j)[0] *= a;
					src.at<Vec3b>(i, j)[1] *= a;
					src.at<Vec3b>(i, j)[2] *= a;
				}
			namedWindow("dst", WINDOW_NORMAL);
			imshow("dst", src);
		}
	}	
	return 0;
}

//鼠标消息回调函数  
static void on_Mouse(int event, int x, int y, int flags, void*)
{
	//处理鼠标不在窗口中的情况  
	if (x < 0 || x >= g_srcImage.cols || y < 0 || y >= g_srcImage.rows)
		return;

	//处理鼠标左键相关消息  
	if (event == CV_EVENT_LBUTTONUP || !(flags & CV_EVENT_FLAG_LBUTTON))
		prevPt = Point(-1, -1);
	else if (event == CV_EVENT_LBUTTONDOWN)
		prevPt = Point(x, y);

	//鼠标左键按下并移动，绘制出线条  
	else if (event == CV_EVENT_MOUSEMOVE && (flags & CV_EVENT_FLAG_LBUTTON))
	{
		Point pt(x, y);
		if (prevPt.x < 0)
			prevPt = pt;
		line(g_maskImage, prevPt, pt, Scalar::all(255), 4, 8, 0);
		line(g_srcImage, prevPt, pt, Scalar::all(255), 4, 8, 0);
		prevPt = pt;
		imshow(WINDOW_NAME1, g_srcImage);
	}
}

//      描述：输出一些帮助信息    
static void ShowHelpText()
{
	printf("\n\n\t\t\t   当前使用的OpenCV版本为：" CV_VERSION);
	printf("\n\n  ----------------------------------------------------------------------------\n");
	//输出一些帮助信息    
	printf("\n\n\n\t欢迎来到【分水岭算法】示例程序~\n\n");
	printf("\t请先用鼠标在图片窗口中标记出大致的区域，\n\n\t然后再按键【1】或者【SPACE】启动算法。"
		"\n\n\t按键操作说明: \n\n"
		"\t\t键盘按键【1】或者【SPACE】- 运行的分水岭分割算法\n"
		"\t\t键盘按键【2】- 恢复原始图片\n"
		"\t\t键盘按键【ESC】- 退出程序\n\n\n");
}

源图像：

?

进行标记的图像：

?

分水岭算法得到的图像：

?

分割后图像：

?

代码的第108-122行是对opencv分水岭算法生成的结果图进行分析，目前对watershed函数生成的index的规律不是很清楚，经测试，并不是按照标记顺序给出index的，具体每一块的index是怎么给出的还需要研究源码

代码第130-138行，目的是想直接根据分水岭算法生成的图像中的index，将背景显示为黑色，需要分割出来的目标物体直接显示，但对index生成的规律还未搞清楚，结果可能不是很稳定

以上部分参考： 毛星云 《OpenCV3编程入门》

-----------------------------------------------------

2019年4月19日增加：

查阅到opencv分水岭算法中，在“循环绘制出轮廓”时用到一个参数compCount，这个参数并不是记录轮廓数目的，它的作用是把每个轮廓设为同一像素值，而maskImage中的像素值就是用1-compcount 的像素值标注的，这样问题又转化为不清楚在查找轮廓时，算法是按照什么样的顺序找出轮廓放入vector中的。