OpenCV入门（九）——图像分割技术之分水岭分割

目录

0x01 分水岭分割

0x02 分水岭分割合并

---

图像分割是利用图像特征灰度、颜色、纹理和形状等将图像中特定的具有独特性质的区域进行划分，进而实现感兴趣的目标的提取。根据分割成因可以分为连续分割和非连续分割。连续分割是指将具有同一灰度级或相同特征的像素划分为不同区域，常见的分割方法有区域生成、阈值分割及聚类分割等。非连续分割是利用像素值突变特性来呈现不同边界区域以实现图像分割，常见的分割方法有点线检测、边缘及能量等。

0x01 分水岭分割

（一）分水岭分割

分水岭分割是基于自然的启发算法来模拟水流通过地形起伏的现象从而研究总结出来的一种分割方法，其基本原理是将图像特征看作地理上的地貌特征，利用像素的**灰度值分布**特征，对每个符合特征的区域进行划分，形成边界以构成**分水岭**。分水岭分割的图像被认为是地形起伏，其中所述梯度大小被解释为**高度**相关信息。对于分水岭算法中的图像像素点，一般需要关注下面三种特征点：

（1）局部最值点：这类点反映的是图像中的局部最小值点或局部最大值点。
（2）交汇边缘点：这类的产生是由于灰度不均匀变化，对应于分水岭中不同地域形成的交接点。
（3）连接区域点：通过局部最小值的像素点向外慢慢扩展，集水区间内的像素点会承担连接作用，同时水域流向会通过连接区域点留到局部最小值的点。

分水岭算法是一种很好的分割相互接触的物体图像的方式，但在边缘分割时精度却存在一定问题，比如在实际应用过程中会出现由于噪点过多，或者其他信息的干扰，会导致局部极值点过多，那么图像分割就难以实现。对于图像某些小区域的背景与前景分割将意味着会有更少的交汇边缘点，这样容易造成过分隔现象，常用的解决方法是合并视觉上的小区域部分，将其视为独立部分进行分割操作。

（二）实现分水岭分割

在计算机视觉中我们常常关注的目标特征是**颜色**和**灰度**，其实也有很多其他的方法比如形状描述子、颜色特征、距离特征等。对于某种场景下的应用，具有独特的纹理对象可以使用一个很好的纹理描述符。针对颜色不同的区域中的单个对象相同的扩展，我们可以使用颜色特征来测量对象的不同部分的相似性。

• 实现使用图像距离变换来实现分水岭分割算法：
• 对源图像进行灰度化，使用OTSU进行二值化操作。
• 对二值化图像进行形态学开操作。
• 利用distanceTransform完成图像距离变换操作。
• 归一化距离变换的统计图像，并计算相应连通域分割块。

代码如下：

#include <iostream>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include "opencv2/core/core.hpp"
#include "opencv2/core/utility.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <opencv2/imgproc/imgproc_c.h>
#include <opencv2/opencv.hpp>
 
using namespace cv;
using namespace std;
 
cv::Mat displaySegResult(cv::Mat& segments, int numOfSegments) //, cv::Mat& image
{
	cv::Mat wshed(segments.size(), CV_8UC3); // 创建对于颜色分量 
	vector<Vec3b> colorTab;
 
	for (int i = 0; i < numOfSegments; 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));
	} //应用不同颜色对每个部分 
 
	for (int i = 0; i < segments.rows; i++)
	{
		for (int j = 0; j < segments.cols; j++)
		{
			int index = segments.at<int>(i, j);
			if (index == -1)
				wshed.at<Vec3b>(i, j) = Vec3b(255, 255, 255);
			else if (index <= 0 || index > numOfSegments)
				wshed.at<Vec3b>(i, j) = Vec3b(0, 0, 0);
			else
				wshed.at<Vec3b>(i, j) = colorTab[index - 1];
		}
	}
 
	//if (image.dims > 0)
		//wshed = wshed * 0.5 + image * 0.5;
	return wshed;
}
 
Mat watershedSegment(Mat& srcImage, int& noOfSegments)
{
	Mat grayMat;
	Mat otsuMat;
	cvtColor(srcImage, grayMat, CV_BGR2GRAY);
	imshow("grayMat", grayMat);
	// 阈值操作
	threshold(grayMat, otsuMat, 0, 255,
		CV_THRESH_BINARY_INV + CV_THRESH_OTSU);
	imshow("otsuMat", otsuMat);
	// 形态学开操作
	morphologyEx(otsuMat, otsuMat, MORPH_OPEN,
		Mat::ones(9, 9, CV_8SC1), Point(4, 4), 2);
	imshow("Mor-openMat", otsuMat);
	// 距离变换
	Mat disTranMat(otsuMat.rows, otsuMat.cols, CV_32FC1);
	distanceTransform(otsuMat, disTranMat, CV_DIST_L2, 3);
	// 归一化
	normalize(disTranMat, disTranMat, 0.0, 1, NORM_MINMAX);
	imshow("DisTranMat", disTranMat);
	// 阈值化分割图像
	threshold(disTranMat, disTranMat, 0.1, 1, CV_THRESH_BINARY);
	//归一化统计图像到0-255
	normalize(disTranMat, disTranMat, 0.0, 255.0, NORM_MINMAX);
	disTranMat.convertTo(disTranMat, CV_8UC1);
	imshow("TDisTranMat", disTranMat);
	//计算标记的分割块
	int i, j, compCount = 0;
	vector<vector<Point> > contours;
	vector<Vec4i> hierarchy;
	findContours(disTranMat, contours, hierarchy,
		CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);
	if (contours.empty())
		return Mat();
	Mat markers(disTranMat.size(), CV_32S);
	markers = Scalar::all(0);
	int idx = 0;
	// 绘制区域块
	for (; idx >= 0; idx = hierarchy[idx][0], compCount++)
		drawContours(markers, contours, idx,
			Scalar::all(compCount + 1), -1, 8,
			hierarchy, INT_MAX);
	if (compCount == 0)
		return Mat();
	//计算算法的时间复杂度
	double t = (double)getTickCount();
	watershed(srcImage, markers);
	t = (double)getTickCount() - t;
	printf("execution time = %gms\n", t * 1000. /
		getTickFrequency());
	Mat wshed = displaySegResult(markers,compCount);
	imshow("watershed transform", wshed);
	noOfSegments = compCount;
	return markers;
}
 
 
 
int main()
{
	Mat srcImage = imread("./image/beauty2.png");
	imshow("src", srcImage);
	int resultall = 0;
	Mat result = watershedSegment(srcImage, resultall);
 
	waitKey(0);
 
}

以这副图像为例：

 第一步为灰度化操作，第一步后的图像是这样的：

 第二步为阈值化操作，第二步的图像是这样的：

 第三步为形态学开操作，其图像为：

这一步的作用为填充操作，为了其更好的分割。

第三步为距离变换：

距离变换是针对二值化图像的一种变换，是计算标识空间点（对目标点）距离的过程，它最终把二值化图像变换为灰度图像（其中每个栅格的灰度值等于它到最近目标点的距离）。在二维空间中，一幅二值化图像可以认为仅仅包含目标和背景两种像素，目标的像素值为1，背景的像素值为0；距离变换的结果不是另一幅二值化图像，而是一幅灰度级图像，即距离图像，图像中的每个像素的灰度值为该像素与距离其最近的背景像素间的距离。

现有的距离变换算法主要采用两类距离测度：非欧式距离和欧式距离。前者常用的有城市街区、棋盘、倒角等距离，算法采用串行扫描实现距离变换，在扫描过程中传递最短距离信息。这些算法简单快速，易于实现，但是得到的仅仅是欧式距离变换（EDT）的一种近似值，在很多应用中不能满足精度要求。

void distanceTransform( InputArray src, 	//8通道二值化图像
						OutputArray dst,	//输出距离图像
                        int distanceType, 	//距离类型
                        int maskSize, 		//距离变换掩膜大小
                        int dstType=CV_32F);

距离类型：**CV_DIST_L1 /CV_DIST_L2  /CV_DIST_C**

maskSize：是 3 或 5，对 CV_DIST_L1 或 CV_DIST_C 的情况，参数值被强制设定为 3, 因为 3×3 mask（由两个数(水平/垂直位量，对角线位移量）组成） 给出 5×5 mask（由三个数组成(水平/垂直位移量，对角位移和 国际象棋里的马步(马走日))） 一样的结果，而且速度还更快。

> 总结来说：distanceTransform方法用于计算图像中每一个非零点距离离自己最近的零点的距离，distanceTransform的第二个Mat矩阵参数dst保存了每一个点与最近零的距离信息，图像上越亮的点，代表了离零点的距离越远。

之后将距离变换后归一化：

normalize函数：该函数归一化输入数组使得它的范数或者数值范围在一定的范围内。简单来说，就是把0-255的数值，到0-1范围内，方便后面的处理。

最后显示的图像如下：

 第四步是，继续阈值化分割图像，使用普通的二值化，再次进行归一：

是距离图像能够更好地表达。

第五步是开始计算标记的分割块，也就是找连通域，最后使用不同的颜色进行区分，效果如下：

 使用源图像来看可能更加直观：

看到这里大概对分水岭算法有些大概的了解，那么就继续深入详述吧：

在上面我们也看到了分水岭处理的效果，那么它的实现的具体依据在哪：分水岭算法中会用到一个很重要的概念——测地线距离(Geodesic Distance)。

测地线距离就是地球表面两点之间的最短路径（可执行路径）的距离，在图论中，Geodesic Distance就是图中两节点之间最短路径的距离，这与平时在几何空间通常用到的Euclidean Distance（欧氏距离），即两点之间最短距离有所区别。

可以看看下图，两个黑点Euclidean Distance是用虚线所表示的距离长度d15，而Geodesic Distance作为实际路径的最短距离，其距离应为沿途实线段距离之和的最小值，即d12 + d23 + d34 + d45。

 

在三维曲面空间中两点间的测地距离就是两点间沿着三维曲面的表面走的最短路径。

那么再回来看看这个算法：

图像的灰度空间很像地球表面的整个地理结构，每个像素的灰度值代表高度。其中**灰度值较大的像素连成的线可以视为分水岭**。其中的水就是二值化阈值后所得到的值，我们将它定为水平面，如果比水平面低的区域，那么就淹没，之后用水填充每一个孤立的山谷（也就是局部最小值）。那么当水平面上升到一定高度的时候，水就会溢出当前的山谷，那么我们就可以在分水岭上修大坝，从而避免两个山谷的水汇集，那么这样的图像就被分成两个像素集，那么一个就是被水淹没的山谷像素集，一个是分水岭线像素集。最终这些大坝的线就对整个图像进行了分区，就实现了对图像的分割。

在该算法中，空间上相邻并且灰度值相近的像素被划分为一个区域。

那么分水岭算法的整个过程：

• 把梯度图像中的所有像素按照灰度值进行分类，并设定一个阈值。（其实也就是我们二值化所得到的阈值），在上面进行了两次分割。
• 找到灰度值最小的像素点（默认标记为灰度值最低点），让threashold从小值开始增长，这些点为起始点。（其实就是每次调用threashold函数的时候，往上抬高一些，把水位再往上涨）
• 水平面在增长的过程中，会碰到周围的邻域像素，测量这些像素到起始点（灰度值最低点）的测地距离，如果小于设定阈值，则将这些像素淹没，否则就在这些像素上设置大坝，就是画上了白点，这样就对这些邻域值进行了分类。
• 随着水平面越来越高，会设置更多更高的大坝，直到灰度值的最大值，所有区域都在分水岭线上分出来了，也就是分区成功。

用上面的算法对图像进行分水岭运算，由于噪声点或其他因素的干扰，可能会得到密密麻麻的小区域，就是把图像分的太细了（过度分割），这因为图像中有非常多的局部极小值点，每个点都会自成一个小区域。

其中的解决方法：

1. 对图像进行高斯平滑操作，抹除很多小的最小值，这些小分区就可以合并。
2. 不从最小值开始增长，可以将相对较高的灰度值的像素作为起始点，那就是我们需要手动设定一个阈值，从标记处开始进行淹没，很多小区域都会被合并为一个区域，这被称为基于图像标记（mark）的分水岭算法。

0x02 分水岭分割合并

上面的方法是使用色度饱和度直方图来提取图像每个区域，并测量了其巴特查里亚距离用于其距离相似性测量。巴特查里亚距离测量的是两个分离或连续区域部分概率分布的相似性，当测量距离很小时，两个区域将被合并成单个部分，当图像的过分隔现象明显时，这种合并可重复进行。分水岭分割合并首先计算分割部分的像素归属，利用直方图信息统计相关特征，对每个分割部分统计直方图对比的相似性，然后根据**相似性**判断分割的两个部分是否需要合并成一个区域。

看看代码：

void segMerge(Mat& image, Mat& segments, int& numSeg)
{
	// 对一个分割部分进行像素统计
	vector<Mat> samples;
	// 统计数据更新
	int newNumSeg = numSeg;
	// 初始化分割部分
	for (int i = 0; i <= numSeg; i++)
	{
		Mat sampleImage;
		samples.push_back(sampleImage);
	}
	// 统计每一个部分
	for (int i = 0; i < segments.rows; i++)
	{
		for (int j = 0; j < segments.cols; j++)
		{
			// 检查每个像素的归属
			int index = segments.at<int>(i, j);
			if (index >= 0 && index < numSeg)
				samples[index].push_back(image(Rect(j, i, 1, 1)));
		}
	}
	// 创建直方图
	vector<MatND> hist_bases;
	Mat hsv_base;
	// 直方图参数设置
	int h_bins = 35;
	int s_bins = 30;
	int histSize[] = { h_bins, s_bins };
	// hue 变换范围 0 to 256, saturation 变换范围0 to 180
	float h_ranges[] = { 0, 256 };
	float s_ranges[] = { 0, 180 };
	const float* ranges[] = { h_ranges, s_ranges };
	// 使用第0与1通道
	int channels[] = { 0, 1 };
	// 直方图生成
	MatND hist_base;
	for (int c = 1; c < numSeg; c++)
	{
		if (samples[c].dims > 0) {
			// 将区域部分转换成hsv
			cvtColor(samples[c], hsv_base, CV_BGR2HSV);
			// 直方图统计
			calcHist(&hsv_base, 1, channels, Mat(),
				hist_base, 2, histSize, ranges, true, false);
			// 直方图归一化
			normalize(hist_base, hist_base, 0, 1,
				NORM_MINMAX, -1, Mat());
			// 添加到统计集
			hist_bases.push_back(hist_base);
		}
		else
		{
			hist_bases.push_back(MatND());
		}
		hist_base.release();
	}
	double similarity = 0;
	vector<bool> mearged;
	for (int k = 0; k < hist_bases.size(); k++)
	{
		mearged.push_back(false);
	}
	// 统计每一个部分的直方图相似
	for (int c = 0; c < hist_bases.size(); c++)
	{
		for (int q = c + 1; q < hist_bases.size(); q++)
		{
			if (!mearged[q])
			{
				// 判断直方图的维度
				if (hist_bases[c].dims > 0 && hist_bases[q].dims > 0)
				{
					// 直方图对比
					similarity = compareHist(hist_bases[c],
						hist_bases[q], CV_COMP_BHATTACHARYYA);
					if (similarity > 0.8)
					{
						mearged[q] = true;
						if (q != c)
						{
							// 区域部分减少
							newNumSeg--;
							for (int i = 0; i < segments.rows; i++)
							{
								for (int j = 0; j < segments.cols; j++)
								{
									int index = segments.at<int>(i, j);
									// 合并
									if (index == q)
									{
										segments.at<int>(i, j) = c;
									}
								}
							}
						}
					}
				}
			}
		}
	}
	numSeg = newNumSeg;
}