分水岭算法（watershed）C++版

#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/objdetect/objdetect.hpp>
#include <opencv2/highgui/highgui.hpp>
#include<vector>
#include<iostream>
#include<queue>
#include<fstream>
cv::Mat marker_mask;
cv::Mat g_markers;
cv::Mat img0, img, img_gray, wshed;
cv::Point_<int> prev_pt(-1, -1);
using std::vector;
using std::queue;
static void my_watershed(cv::Mat img, cv::Mat& markers, int comp_count);
static void my_watershed(cv::Mat img0, cv::Mat& markers, int comp_count)
{
	cv::Mat gray = cv::Mat(cv::Size(img0.rows, img0.cols), 8, 1);
	cv::cvtColor(img0, gray, CV_BGR2GRAY);
	cv::Mat imge = cv::Mat(cv::Size(img0.rows, img0.cols), 8, 1);
	cv::Sobel(gray, imge, CV_8U, 1, 1);
	vector<queue<cv::Point2i>*>Labeleddata;//图像中各连通区域的点
	queue<cv::Point2i>* pque;//某连通区域已包含的点
	queue<cv::Point2i> quetem; //用于提取某连通区域中输入种子点中的初始种子点
	vector<int*> SeedCounts;
	int* Array;
	cv::Point2i temp;
	int row = imge.rows, col = imge.cols;
	cv::Mat marker_saved = markers.clone();
	bool up, down, right, left, uplef, uprig, downlef, downrig;
	int m, n;
	for (int i = 0; i < comp_count; i++)
	{
		Array = new int[256];
		SeedCounts.push_back(Array);//统计某waterlevel的各个连通区域中种子点的个数
		pque = new queue<cv::Point2i>[256];
		Labeleddata.push_back(pque);//存储该连通区域中种子生长所得的点		
	}
	for (int i = 0; i < row; i++)
		for (int j = 0; j < col; j++)
		{
			if (markers.at<int>(i, j) > 0)
			{
				temp.x = i;
				temp.y = j;
				quetem.push(temp);
				int num = markers.at<int>(i, j);
				markers.at<int>(i, j) = -1;//该点已处理,其他种子点生长时将绕过该点
				while (!quetem.empty())
				{
					up = down = right = left = uplef = uprig = downlef = downrig = false;
					temp = quetem.front(); //提取出一个点，在该点的八连通区域内寻找可生长点
					m = temp.x;
					n = temp.y;
					quetem.pop();

					if (m - 1 >= 0)//若上方可生长则添加为新种子
					{
						if (markers.at<int>(m - 1, n) == num)
						{
							temp.x = m - 1;
							temp.y = n;
							quetem.push(temp);
							markers.at<int>(m - 1, n) = -1;
						}
						else {
							up = true;
						}
					}
					if (m - 1 >= 0 && n - 1 >= 0)
					{
						if (markers.at<int>(m - 1, n - 1) == num)
						{
							temp.x = m - 1;
							temp.y = n - 1;
							quetem.push(temp);
							markers.at<int>(m - 1, n - 1) = -1;
						}
						else {
							uplef = true;
						}
					}
					if (m + 1 <= row - 1)
					{
						if (markers.at<int>(m + 1, n) == num)
						{
							temp.x = m + 1;
							temp.y = n;
							quetem.push(temp);
							markers.at<int>(m + 1, n) = -1;
						}
						else {
							down = true;
						}
					}
					if (m + 1 <= row - 1 && n + 1 <= col - 1)
					{
						if (markers.at<int>(m + 1, n + 1) == num)
						{
							temp.x = m + 1;
							temp.y = n + 1;
							quetem.push(temp);
							markers.at<int>(m + 1, n + 1) = -1;
						}
						else {
							downrig = true;
						}
					}
					if (n + 1 <= col - 1)
					{
						if (markers.at<int>(m, n + 1) == num)
						{
							temp.x = m;
							temp.y = n + 1;
							quetem.push(temp);
							markers.at<int>(m, n + 1) = -1;
						}
						else {
							right = true;
						}
					}
					if (m - 1 >= 0 && n + 1 <= col - 1)
					{
						if (markers.at<int>(m - 1, n + 1) == num)
						{
							temp.x = m - 1;
							temp.y = n + 1;
							quetem.push(temp);
							markers.at<int>(m - 1, n + 1) = -1;
						}
						else {
							uprig = true;
						}
					}
					if (n - 1 >= 0)
					{
						if (markers.at<int>(m, n - 1) == num)
						{
							temp.x = m;
							temp.y = n - 1;
							quetem.push(temp);
							markers.at<int>(m, n - 1) = -1;
						}
						else {
							left = true;
						}
					}
					if (m + 1 <= row - 1 && n - 1 >= 0)
					{
						if (markers.at<int>(m + 1, n - 1) == num)
						{
							temp.x = m + 1;
							temp.y = n - 1;
							quetem.push(temp);
							markers.at<int>(m + 1, n - 1) = -1;
						}
						else {
							downlef = true;
						}
					}
					//八连通区域中有未标记点，则该点属于初始种子点
					if (up || down || right || left || uplef || downlef || uprig || downrig)
					{
						temp.x = m;
						temp.y = n;
						Labeleddata[comp_count - 1][imge.at<uchar>(m, n)].push(temp);
						SeedCounts[comp_count - 1][imge.at<uchar>(m, n)]++;
					}
				}
			}
		}
	bool active;
	int waterlevel;
	for (waterlevel = 0; waterlevel < 180; waterlevel++)
	{
		active = true;
		while (active) //当1-count_com个连通区域都无可生长点时结束循环
		{
			active = false;
			for (int i = 0; i < comp_count; i++)//将区域i中将waterlevel梯度以下的点用于区域增长
			{
				//区域增长，经过多次迭代，直至该区域，该waterlevel无可生长点。
				if (!Labeleddata[i][waterlevel].empty())
				{
					active = true;
					//SeedCount中保留了前一轮生长后各区域，各waterlevel的种子点个数，
					//本轮生长结束后，将根据Labeleddata中的元素个数更新
					while (SeedCounts[i][waterlevel] > 0)
					{
						SeedCounts[i][waterlevel]--;
						temp = Labeleddata[i][waterlevel].front();
						Labeleddata[i][waterlevel].pop();
						m = temp.x;
						n = temp.y;
						int num = marker_saved.at<int>(m, n);
						if (m - 1 >= 0)
						{
							if (!marker_saved.at<int>(m - 1, n))//上方点未处理过
							{
								temp.x = m - 1;
								temp.y = n;
								marker_saved.at<int>(m - 1, n) = num;
								if (imge.at<uchar>(m - 1, n) <= waterlevel)
									Labeleddata[i][waterlevel].push(temp);
								else {
									//本次生长不处理，可能在waterlevel变化到某值时再用于生长
									Labeleddata[i][imge.at<uchar>(m - 1, n)].push(temp);
									SeedCounts[i][imge.at<uchar>(m - 1, n)]++;
								}
							}
						}
						if (m + 1 <= row - 1)
						{
							if (!marker_saved.at<int>(m + 1, n))
							{
								temp.x = m + 1;
								temp.y = n;
								marker_saved.at<int>(m + 1, n) = num;
								if (imge.at<uchar>(m + 1, n) <= waterlevel)
									Labeleddata[i][waterlevel].push(temp);
								else {
									Labeleddata[i][imge.at<uchar>(m + 1, n)].push(temp);
									SeedCounts[i][imge.at<uchar>(m + 1, n)]++;
								}
							}
						}
						if (n + 1 <= col - 1)
						{
							if (!marker_saved.at<int>(m, n + 1))
							{
								temp.x = m;
								temp.y = n + 1;
								marker_saved.at<int>(m, n + 1) = num;
								if (imge.at<uchar>(m, n + 1) <= waterlevel)
									Labeleddata[i][waterlevel].push(temp);
								else {
									Labeleddata[i][imge.at<uchar>(m, n + 1)].push(temp);
									SeedCounts[i][imge.at<uchar>(m, n + 1)]++;
								}
							}
						}
						if (n - 1 >= 0)
						{
							if (!marker_saved.at<int>(m, n - 1))
							{
								temp.x = m;
								temp.y = n - 1;
								marker_saved.at<int>(m, n - 1) = num;
								if (imge.at<uchar>(m, n - 1) <= waterlevel)
									Labeleddata[i][waterlevel].push(temp);
								else
								{
									Labeleddata[i][imge.at<uchar>(m, n - 1)].push(temp);
									SeedCounts[i][imge.at<uchar>(m, n - 1)]++;
								}
							}
						}
					}
					SeedCounts[i][waterlevel] = Labeleddata[i][waterlevel].size();
				}

			}
		}
	}
	markers = marker_saved.clone();
}


static void mouse_event(int event, int x, int y, int flags, void*)
{
	if (img.rows == 0)
		return;
	if (event == CV_EVENT_LBUTTONUP || !(flags&CV_EVENT_FLAG_LBUTTON))
		prev_pt = cv::Point_<int>(-1, -1);
	else if (event == CV_EVENT_LBUTTONDOWN)
		prev_pt = cv::Point2i(x, y);
	else if (event == CV_EVENT_MOUSEMOVE && (flags&CV_EVENT_FLAG_LBUTTON))
	{
		cv::Point2i pt(x, y);
		if (prev_pt.x < 0)
			prev_pt = pt;
		cv::line(marker_mask, prev_pt, pt, cv::Scalar(255, 255, 255), 1, 8, 0);
		cv::line(img, prev_pt, pt, cv::Scalar(255, 255, 255), 1, 8, 0);
		prev_pt = pt;
		cv::imshow("image", img);
	}
}
int main()
{
	img0 = cv::imread("D:/0.png", 1);
	img = img0.clone();
	CvRNG rng = cvRNG(-1);
	img_gray = img0.clone();
	wshed = img0.clone();
	marker_mask = cv::Mat(cv::Size(img0.cols, img0.rows), 8, 1);
	g_markers = cv::Mat(cv::Size(img0.cols, img0.rows), CV_32S, 1);
	cv::cvtColor(img, marker_mask, CV_BGR2GRAY);
	cv::cvtColor(marker_mask, img_gray, CV_GRAY2BGR);
	for (int i = 0; i < marker_mask.rows; i++)
		for (int j = 0; j < marker_mask.cols; j++)
			marker_mask.at<unsigned char>(i, j) = 0;
	for (int i = 0; i < g_markers.rows; i++)
		for (int j = 0; j < g_markers.cols; j++)
			g_markers.at<int>(i, j) = 0;
	cv::imshow("image", img);
	cv::imshow("watershed transform", wshed);
	cv::setMouseCallback("image", mouse_event, 0);
	for (;;)
	{
		int c = cv::waitKey(0);
		if ((char)c == 27)
			break;
		if ((char)c == 'r')
		{
			for (int i = 0; i < marker_mask.rows; i++)
				for (int j = 0; j < marker_mask.cols; j++)
					marker_mask.at<unsigned char>(i, j) = 0;
			img0.copyTo(img);
			cv::imshow("image", img);
		}
		if ((char)c == 'w' || (char)c == ' ')
		{
			vector<vector<cv::Point>> contours;
			CvMat* color_tab = 0;
			int comp_count = 0;
			cv::findContours(marker_mask, contours, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE, cv::Point(0, 0));
			for (int i = 0; i < g_markers.rows; i++)
				for (int j = 0; j < g_markers.cols; j++)
					g_markers.at<int>(i, j) = 0;
			vector<vector<cv::Point> >::iterator iter = contours.begin();
			for (int i = 0; i < (int)contours.size(); i++)
			{
				cv::drawContours(g_markers, contours, i, cv::Scalar::all(comp_count + 1),
					1, 8, vector<cv::Vec4i>());
				comp_count++;
			}

			if (comp_count == 0)
				continue;
			color_tab = cvCreateMat(1, comp_count, CV_8UC3);
			for (int 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);
			}
			cv::Mat temp = g_markers.clone();

			double t = (double)cvGetTickCount();
			//my_watershed(img0,g_markers,comp_count);
			cv::watershed(img0, g_markers);
			t = (double)cvGetTickCount() - t;
			std::cout << "exec time= " << t / (cvGetTickFrequency()*1000.) << std::endl;
			for (int i = 0; i < g_markers.rows; i++)
				for (int j = 0; j < g_markers.cols; j++)
				{
					int idx = g_markers.at<int>(i, j);
					uchar* dst = &wshed.at<uchar>(i, j * 3);
					if (idx == -1)
						dst[0] = dst[1] = dst[2] = (uchar)255;
					else if (idx <= 0 || idx > comp_count)
						dst[0] = dst[1] = dst[2] = (uchar)8;
					else {
						uchar* ptr = color_tab->data.ptr + (idx - 1) * 3;
						dst[0] = ptr[0]; dst[1] = ptr[1]; dst[2] = ptr[2];
					}
				}
			cv::addWeighted(wshed, 0.5, img_gray, 0.5, 0, wshed);
			cv::imshow("watershed transform", wshed);
			cvReleaseMat(&color_tab);
		}
	}
	return 0;
}

运行效果：

还有一种方法：

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

using namespace cv;
using namespace std;

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

Mat g_maskImage, g_srcImage;
Point prevPt(-1, -1);

static void ShowHelpText()
{
	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");
}

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 == EVENT_LBUTTONUP || !(flags & EVENT_FLAG_LBUTTON))
		prevPt = Point(-1, -1);
	else if (event == EVENT_LBUTTONDOWN)
		prevPt = Point(x, y);

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

int main(int argc, char** argv)
{
	//【0】显示帮助文字
	ShowHelpText();

	//【1】载入原图
	g_srcImage = imread("D:/0.jpg");

	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(WINDOW_NAME1, g_srcImage);
		}

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

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

			//轮廓为空时的处理
			if (contours.empty() || hierarchy.size() == 0)
				continue;

			int compCount = hierarchy.size();

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

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

			// 生成随机颜色
			vector<Vec3b> colorTab;
			for (int 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));
			}

			// 计算处理时间并输出到窗口中
			watershed(srcImage, maskImage);

			//双层循环，将分水岭图像遍历存入watershedImage中
			Mat watershedImage(maskImage.size(), CV_8UC3);
			for (int i = 0; i < maskImage.rows; i++)
				for (int j = 0; j < maskImage.cols; j++)
				{
					int index = maskImage.at<int>(i, j);
					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
						watershedImage.at<Vec3b>(i, j) = colorTab[index - 1];
				}

			//imshow("1", watershedImage);
			//imshow("2", grayImage);

			//混合灰度图和分水岭效果图并显示最终的窗口
			//watershedImage = watershedImage * 0.5 + grayImage * 0.5;
			addWeighted(watershedImage, 0.5, grayImage, 0.5, 0, watershedImage);
			imshow(WINDOW_NAME2, watershedImage);
		}
	}

	return 0;
}