分水岭算法及其实现

1 - 算法描述
1.1 分水岭算法的原理
  分水岭的概念是以三维方式来形象化一幅图像为基础的:两个空间坐标再加上强度。在这种“地形学”解释中,考虑三种类型的点:(a)局部最小值点,该点对应一个盆地的最低点,当我们在盆地里滴一滴水的时候,由于重力作用,水最终会汇聚到该点。注意:可能存在一个最小值面,该平面内的都是最小值点。(b)盆地的其它位置点,该位置滴的水滴会汇聚到局部最小点。(c)盆地的边缘点,是该盆地和其它盆地交接点,在该点滴一滴水,会等概率的流向任何一个盆地。对于一个特定的区域最小值,满足条件(b)的点的集合称为该最小值的汇水盆地或分水岭。满足条件(c)的点形成地面的锋线,称为分割线或分水线。


图1 分水岭算法中三种类型的点

  基于这些概念的分割算法的主要目标是找出分水线。假设在每个区域的最小值上打一个洞,并且让水通过洞以均匀的速率上升,从低到高淹没整个地形。当不同汇水盆地中上升的水聚集时,修建一个水坝来阻止这种聚合。水将达到在水线上只能见到各个水坝的顶部的程度。这些大坝的边界对应于分水岭的分割线。因此,他们是由分水岭算法提取出来的边界。

1.2 水坝构建
  在实际操作中,水位的上升常利用灰度0-255的分层处理来实现,当处理到第N层时,灰度小于N的像素就会被淹没。对于大多数情况,聚水盆的生长速度是不定的,也许在第N-1次处理时,两个聚水盆之间还有相当的距离,而在第N次处理中他们就已经相互连通了,这为堤坝的建立带来了一定的难度。


图2 水坝构造示意图

  图2中a显示了两个相邻的聚水盆,其中颜色越深的区域海拔越低,当水位逐渐上升到第N-1阶段时,黑色区域首先被淹没,如图2中b所示,这时图中存在两个彼此分开的水域,他们之间还有相当一段距离。然而当水面继续上升至第N阶段时,图中的水域已经连通,如图2中c所示。
  为了在水域之间适当的位置构建堤坝防止水域的连通,可以借助图像形态学处理中的膨胀方法。由于已知两个水域会在第N-1阶段至第N阶段间连通,回到第N-1阶段的图像如图2中b所示,尝试对图中不同水域分别进行膨胀操作,并在膨胀中寻找两个区域的交点。如图2中d所示为两个水域分别进行一次膨胀操作的结果,从结果中可以看出两个区域仍然没有交点,因此还不能形成堤坝。继续对两个水域进行膨胀操作,在对右边区域的膨胀操作中,两个区域出现了交点,标记这些交点并在交点处修建堤坝,如图2中e所示,其中堤坝位置用交叉图案表示。继续进行膨胀操作,并不断标记新的交点,最终就能够得到完整的水坝轮廓。值得注意的是,在算法实现中,应对每次膨胀结果与第N阶段的淹没区域进行比较,保证堤坝位置在第N阶段的淹没区域之内,算法一直进行直到经过膨胀的水域完全覆盖第N阶段中淹没区域为止。通过以上步骤就能够获得较为准确的堤坝位置,如图2中f所示。
  由于实际中处理的图像往往比较复杂,图中灰度极小值点较多,这就导致聚水盆数目过多,从而造成图像的过度分割,比如下面的图3,这样的分割效果是毫无用处的。解决图像的过度分割问题最直接的方法就是进行区域的合并,如果在分水岭算法中定义聚水盆的最小深度,就能比较好地解决过度分割的问题。


图3 过度分割的结果

  聚水盆的深度可以定义为盆地中海拔最低点与当前水面高度的差值,在计算时可以用当前已淹没的灰度值减去区域的灰度最小值来得到。当水位上升造成两个不同的水域连接时,判断两个水域所在聚水盆的深度是否都大于定义的最小深度,若是则在两个区域间构建堤坝,否则对两个水域进行合并。运用这样的方法,就可以通过改变最小深度的值来控制图像分割的程度,从而获得图像最有效的分割结果。

2 - 分水岭算法的OpenCv实现
  函数watershed实现的分水岭算法是基于标记的分割算法中的一种。在把图像传给函数之前,需要大致勾画标记出图像中的期望进行分割的区域,他们被标记为正指数。所以,每一个区域都会被标记为像素值1、2、3等,表示成为一个或者多个连接组件。这些标记的值可以使用findContours()函数和drawContours()函数由二进制的掩码检测出来。不难理解,这些标记就是即将绘制出来的分割区域的“种子”,而没有标记清楚的区域,被置为0。在函数输出中,每一个标记中的像素被设置为“种子”的值,而区域间的值被设置为-1。
  C++: void watershed(InputArray image, InputOutputArray markers)
  第一个参数,InputArray类型的src,输入图像,即原图像,填Mat类的对象即可,且需为8位三通道的彩色图像。
  第二个参数,InputOutputArray类型的markers,函数调用后的运算结果存在这里,输入/输出32位单通道图像的标记结果。即这个参数用于存放函数调用后的输出结果,需和原图像有一样的尺寸和类型。


图4 原始图

  编写程序,对图4进行基于分水岭算法的图像分割。


图5 期望分割的区域

  直接使用图像处理工具箱的提供的现成函数进行分割往往很难达到预期效果,而预先勾画标记出图像中期望进行分割的区域,在应用分水岭算法就会取得较好的分割效果。如图5所示,其中左图显示了原始图中期望进行分割的各个部分,右图更清楚的显示了期望分割部分的空间位置。


图6 分割结果

  图6显示了分割的结果。分割的各个部分使用随机的颜色进行表示;若存在没有标记清楚的区域,则将此区域用纯黑色表示;用于划分区域的线条,即被设置为“种子”的像素值,用纯白色表示。左图显示了原始图像的分割结果,右图显示了分割结果和原始图按等比例进行混合的效果。

3 - OpenCv代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;
using namespace std;
//-----------------------------------【宏定义部分】--------------------------------------------
// 描述:定义一些辅助宏
//------------------------------------------------------------------------------------------------
#define WINDOW_NAME1 "【程序窗口1】" //为窗口标题定义的宏
#define WINDOW_NAME2 "【分水岭算法效果图】" //为窗口标题定义的宏
//-----------------------------------【全局函变量声明部分】--------------------------------------
// 描述:全局变量的声明
//-----------------------------------------------------------------------------------------------
Mat g_maskImage, g_srcImage;
Point prevPt(-1, -1);
static void on_Mouse(int event, int x, int y, int flags, void*);
int main(int argc, char** argv)
{
//【0】改变console字体颜色
system("color 6F");
//【1】载入原图并显示,初始化掩膜和灰度图
g_srcImage = imread("1.jpg", 1);
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, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);
//轮廓为空时的处理
if (contours.empty())
continue;
//拷贝掩膜
Mat maskImage(g_maskImage.size(), CV_32S);
maskImage = Scalar::all(0);
//循环绘制出轮廓 //感觉这里绘制轮廓不重要,重要的是几个 因为后面又覆盖了maskImage watershed( srcImage, maskImage );
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);
for (i = 0; i < maskImage.rows; i++)
for (j = 0; j < maskImage.cols; j++)
{
int index = maskImage.at<int>(i, j); //不同的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
watershedImage.at<Vec3b>(i, j) = colorTab[index - 1];
}
//混合灰度图和分水岭效果图并显示最终的窗口
watershedImage = watershedImage * 0.5 + grayImage * 0.5;
imshow(WINDOW_NAME2, watershedImage);
}
}
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), 5, 8, 0);
line(g_srcImage, prevPt, pt, Scalar::all(0), 5, 8, 0);
prevPt = pt;
imshow(WINDOW_NAME1, g_srcImage);
}
}