花老湿学习OpenCV:分水岭原理和实现watershed()

在很多实际应用中，我们需要分割图像，分割方法有多种经典的分割方法：

常见图像分割方法：

1、基于边缘检测的方法：

       此方法主要是通过检测区域的边缘进行分割，利用区域之间的特征的不一致性，首先检测图像中的边缘点，然后按照一定的方法把这些边缘点进行全部连接起来，从而构成分割区域。图像中的边缘通常是灰度，颜色或者纹理，其中基于灰度的方法很普遍，许多边缘检测算子利用灰度来检测图像的梯度，Roberts 算子、Laplace 算子、Prewitt 算子、Sobel 算子、Rosonfeld算子、Kirsch 算子以及Canny 。边缘检测算法比较适合边缘灰度值过渡比较显著且噪声较小的简单图像的分割。对于边缘比较复杂以及存在较强噪声的图像，则面临抗噪性和检测精度的矛盾。若提高检测精度，则噪声产生的伪边缘会导致不合理的轮廓：若提高抗噪性，则会产生轮廓漏检和位置偏差。

2.阈值分割方法：

       阈值分割是最古老的分割技术，也是最简单实用的。许多情况下，图像中目标区域与背景区域或者说不同区域之间其灰度值存在差异，此时可以将灰度的均一性作为依据进行分割。阈值分割即通过一个或几个阈值将图像分割成不同的区域。阈值分割方法的核心在于如何寻找适当的阈值。最常用的阈值方法是基于灰度直方图的方法，如最大类间方差法(OTSU)、最小误差法、最大熵法等。此类方法通常对整幅图像使用固定的全局阈值，如果图像中有阴影或亮度分布不均等现象，分割效果会受到影响。基于局部阈值的分割方法对图像中的不同区域采用不同的阈值，相对于全局阈值方法具有更好的分割效果，该方法又称为自适应阈值方法。

3.区域生长：

      区域生长方法[46]也是一种常用的区域分割技术，其基本思路是首先定义一个生长准则，然后在每个分割区域内寻找一个种子像素，通过对图像进行扫描，依次在种子点周围邻域内寻找满足生长准则的像素并将其合并到种子所在的区域，然后再检查该区域的全部相邻点，并把满足生长准则的点合并到该区域，不断重复该过程直到找不到满足条件的像素为止。该方法的关键在于种子点的位置、生长准则和生长顺序。

4.分水岭算法：

      是以数学形态学作为基础的一种区域分割方法。其基本思想是将梯度图像看成是假想的地形表面，每个像素的梯度值表示该点的海拔高度。原图中的平坦区域梯度较小，构成盆地，边界处梯度较大构成分割盆地的山脊。分水岭算法模拟水的渗入过程，假设水从最低洼的地方渗入，随着水位上升，较小的山脊被淹没，而在较高的山脊上筑起水坝，防止两区域合并。当水位达到最高山脊时，算法结束，每一个孤立的积水盆地构成一个分割区域。由于受到图像噪声和目标区域内部的细节信息等因素影响，使用分水岭算法通常会产生过分割现象，分水岭算法一般是作为一种预分割方法，与其它分割方法结合使用，以提高算法的效率或精度。
 

分水岭算法：

1.分水岭比较经典的计算方法是L．Vincent于1991年在PAMI上提出的[1]。传统的分水岭分割方法，是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆地，而集水盆地的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸人水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝如下图所示，即形成分水岭。

 

然而基于梯度图像的直接分水岭算法容易导致图像的过分割，产生这一现象的原因主要是由于输入的图像存在过多的极小区域而产生许多小的集水盆地，从而导致分割后的图像不能将图像中有意义的区域表示出来。所以必须对分割结果的相似区域进行合并。

[1]L.Vincent, P Soille. Watersheds in digital space: An efficientalgorithms based on immersion simulation[J]. IEEE Trans. on Pattern Analysisand Machine Intelligence, 1991, 13(6): 583-598.

2.改进的分水岭算法基本原理
因为传统分水岭算法存在过分割的不足，OpenCV提供了一种改进的分水岭算法，使用一系列预定义标记来引导图像分割的定义方式。使用OpenCV的分水岭算法cv::wathershed，需要输入一个标记图像，图像的像素值为32位有符号正数（CV_32S类型），每个非零像素代表一个标签。它的原理是对图像中部分像素做标记，表明它的所属区域是已知的。分水岭算法可以根据这个初始标签确定其他像素所属的区域。传统的基于梯度的分水岭算法和改进后基于标记的分水岭算法示意图如下图所示。

                                                   传统基于梯度的分水岭算法和基于标记的分水岭算法原理图

从上图可以看出，传统基于梯度的分水岭算法由于局部最小值过多造成分割后的分水岭较多。而基于标记的分水岭算法，水淹过程从预先定义好的标记图像（像素）开始，较好的克服了过度分割的不足。本质上讲，基于标记点的改进算法是利用先验知识来帮助分割的一种方法。因此，改进算法的关键在于如何获得准确的标记图像，即如何将前景物体与背景准确的标记出来。
 

API：

算法会根据markers传入的轮廓作为种子（也就是所谓的注水点），对图像上其他的像素点根据分水岭算法规则进行判断，并对每个像素点的区域归属进行划定，直到处理完图像上所有像素点。而区域与区域之间的分界处的值被置为“-1”，以做区分。 

 

代码示例：

#include "pch.h"
#include <iostream>
#include "opencv2/opencv.hpp"

using namespace std;
using namespace cv;


int main()
{
	Mat src;
	src = imread("F:\\visual studio\\Image\\coin3.jpg");
	if (!src.data)
	{
		printf("could not load image...\n");
		return -1;
	}

	namedWindow("src", CV_WINDOW_AUTOSIZE);
	imshow("src", src);
	
	//转化为灰度图
	Mat graysrc;
	cvtColor(src, graysrc, COLOR_BGR2GRAY);
	//imshow("graysrc", graysrc);

	//转化为二值图像
	Mat bin;
	threshold(graysrc, bin, 0, 255, CV_THRESH_BINARY_INV|CV_THRESH_OTSU);
	imshow("bin1", bin);
	
	//去除噪声 开运算：先腐蚀后膨胀
	Mat element = getStructuringElement(MORPH_RECT, Size(3, 3), Point(-1, -1));
	morphologyEx(bin, bin, MORPH_OPEN, element, Point(-1, -1));
	imshow("bin2", bin);

	//膨胀得到确定为背景的区域
	Mat imgbg;
	dilate(bin, imgbg, element, Point(-1, -1), 3);
	imshow("imgbg", imgbg);

	//距离变换加腐蚀得到确定为前景的区域
	Mat imgthin = Mat::zeros(src.size(), CV_32FC1);
	distanceTransform(bin, imgthin, DIST_L1, 3);
	normalize(imgthin, imgthin, 0, 1, NORM_MINMAX); //浮点类型需要转化为0-1之间
	imshow("imgthin", imgthin);
	
	Mat imgfg;
	threshold(imgthin, imgfg, 0.4, 1, CV_MINMAX);
	erode(imgfg, imgfg, element, Point(-1, -1));
	imshow("imgfg", imgfg);

	//相减得到未知区域，这部分可能是前景可能是背景
	Mat unkown(src.size(),CV_8UC1);
	imgfg.convertTo(imgfg, CV_8UC1);
	normalize(imgfg,imgfg, 0, 255, CV_MINMAX);
	unkown = imgbg - imgfg;
	imshow("unkown", unkown);

	//创建标记图像
	vector<vector<Point>> Contours;
	vector<Vec4i> hiearchy;
	
	//首先标记确定的前景图像 编号分别为2,3,4...
	findContours(imgfg, Contours, hiearchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE, Point(-1, -1));

	Mat markers = Mat::zeros(src.size(), CV_32SC1);
	for (int i = 0; i < Contours.size(); i++)
	{
		drawContours(markers, Contours, i, Scalar::all(i + 2), FILLED, LINE_AA, hiearchy);
	}
	
	//然后标记确定的背景图像 编号为1
	bitwise_not(imgbg, imgbg);
	imgbg.convertTo(imgbg, CV_32SC1);
	normalize(imgbg, imgbg, 0, 1, CV_MINMAX);
	//完成标记 此时只有不确定区域的值为0
	markers = markers + imgbg;
	imshow("markers", markers * 10000);

	//分水岭运算  此时不确定区域被填充为1,2,3...
	watershed(src, markers);
	
	RNG rng(getTickCount());
	Mat dst = Mat::zeros(markers.size(), CV_8UC3);
	
	// 生成随机颜色
	vector<Vec3b> colors;
	for (size_t i = 0; i < Contours.size(); i++) {
			int r = rng.uniform(0, 255);
		    int g = rng.uniform(0, 255);
		    int b = rng.uniform(0, 255);
		    colors.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));		
	}


	for (int row = 0; row < markers.rows; row++) {
		for (int col = 0; col < markers.cols; col++) {
			int index = markers.at<int>(row, col);		
			if (index > 1 && index <= (Contours.size() + 1)) {
				//此时为我们的前景区域
				dst.at<Vec3b>(row, col) = colors[index-2];
			}
			else {
				//此时为我们的背景区域
				dst.at<Vec3b>(row, col) = Vec3b(0, 0, 0);
			}
		}
	}
	imshow("dst", dst);
	waitKey(0);
}

效果展示：