分水岭分割方法（Watershed Segmentation），是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸入水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。
分水岭的计算过程是一个迭代标注过程。分水岭比较经典的计算方法是L. Vincent提出的。在该算法中，分水岭计算分两个步骤，一个是排序过程，一个是淹没过程。首先对每个像素的灰度级进行从低到高排序，然后在从低到高实现淹没过程中，对每一个局部极小值在h阶高度的影响域采用先进先出（FIFO）结构进行判断及标注。
分水岭变换得到的是输入图像的集水盆图像，集水盆之间的边界点，即为分水岭。显然，分水岭表示的是输入图像极大值点。因此，为得到图像的边缘信息，通常把梯度图像作为输入图像，即：
grad(f(x,y))=((f(x-1,y)-f(x+1,y))^2 + (f(x,y-1)-f(x,y+1))^2)^0.5
式中，f(x,y)表示原始图像，grad(.)表示梯度运算。
分水岭算法对微弱边缘具有良好的响应，图像中的噪声、物体表面细微的灰度变化，都会产生过度分割的现象。但同时应当看出，分水岭算法对微弱边缘具有良好的响应，是得到封闭连续边缘的保证的。另外，分水岭算法所得到的封闭的集水盆，为分析图像的区域特征提供了可能。
为消除分水岭算法产生的过度分割，通常可以采用两种处理方法，一是利用先验知识去除无关边缘信息。二是修改梯度函数使得集水盆只响应想要探测的目标。
为降低分水岭算法产生的过度分割，通常要对梯度函数进行修改，一个简单的方法是对梯度图像进行阈值处理，以消除灰度的微小变化产生的过度分割。即：
g(x,y)=max(grad(f(x,y)),gθ)
式中，gθ表示阈值。
程序可采用方法：用阈值限制梯度图像以达到消除灰度值的微小变化产生的过度分割，获得适量的区域，再对这些区域的边缘点的灰度级进行从低到高排序，然后在从低到高实现淹没的过程，梯度图像用Sobel算子计算获得。对梯度图像进行阈值处理时，选取合适的阈值对最终分割的图像有很大影响，因此阈值的选取是图像分割效果好坏的一个关键。缺点：实际图像中可能含有微弱的边缘，灰度变化的数值差别不是特别明显，选取阈值过大可能会消去这些微弱边缘。
下面用C++实现分水岭算法：
|
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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
|
#define _USE_MATH_DEFINES #include <cstddef> #include <cstdlib> #include <cstring> #include <climits> #include <cfloat> #include <ctime> #include <cmath> #include <cassert> #include <vector> #include <stack> #include <queue> using namespace std; typedef void GVVoid; typedef bool GVBoolean; typedef char GVChar; typedef unsigned char GVByte; typedef short GVInt16; typedef unsigned short GVUInt16; typedef int GVInt32; typedef unsigned int GVUInt32; typedef long long GVInt64; typedef unsigned long long GVUInt64; typedef float GVFloat32; typedef double GVFloat64; const GVBoolean GV_TRUE = true; const GVBoolean GV_FALSE = false; const GVByte GV_BYTE_MAX = UCHAR_MAX; const GVInt32 GV_INT32_MAX = INT_MAX; const GVInt32 GV_INT32_MIX = INT_MIN; const GVInt64 GV_INT64_MAX = LLONG_MAX; const GVInt64 GV_INT64_MIN = LLONG_MIN; const GVFloat32 GV_FLOAT32_MAX = FLT_MAX; const GVFloat32 GV_FLOAT32_MIN = FLT_MIN; const GVFloat64 GV_FLOAT64_MAX = DBL_MAX; const GVFloat64 GV_FLOAT64_MIN = DBL_MIN; class GVPoint; class GVPoint { public: GVInt32 x; GVInt32 y; public: GVPoint() : x(0), y(0) { } GVPoint(const GVPoint &obj) : x(obj.x), y(obj.y) { } GVPoint(GVInt32 x, GVInt32 y) : x(x), y(y) { } public: GVBoolean operator ==(const GVPoint &right) const { return ((x == right.x) && (y == right.y)); } GVBoolean operator !=(const GVPoint &right) const { return (!(x == right.x) || !(y == right.y)); } }; /* * <Parameter> * <image> image data; * <width> image width; * <height> image height; * <label out> image of labeled watersheds. */GVVoid gvWatershed( const GVByte *image, GVInt32 width, GVInt32 height, GVInt32 *label) { // Local constants const GVInt32 WSHD = 0; const GVInt32 INIT = -1; const GVInt32 MASK = -2; const GVPoint FICT_PIXEL = GVPoint(~0, ~0); // Image statistics and sorting GVInt32 size = width * height; GVInt32 *image_stat = new GVInt32[GV_BYTE_MAX + 1]; GVInt32 *image_space = new GVInt32[GV_BYTE_MAX + 1]; GVPoint *image_sort = new GVPoint[size]; ::memset(image_stat, 0, sizeof (GVInt32) * (GV_BYTE_MAX + 1)); ::memset(image_space, 0, sizeof (GVInt32) * (GV_BYTE_MAX + 1)); ::memset(image_sort, 0, sizeof (GVPoint) * size); for (GVInt32 i = 0; !(i == size); ++i) { image_stat[image[i]]++; } for (GVInt32 i = 0; !(i == GV_BYTE_MAX); ++i) { image_stat[i + 1] += image_stat[i]; } for (GVInt32 i = 0; !(i == height); ++i) { for (GVInt32 j = 0; !(j == width); ++j) { GVByte space = image[i * width + j]; GVInt32 index = image_stat[space] - (++image_space[space]); image_sort[index].x = j; image_sort[index].y = i; } } for (GVInt32 i = GV_BYTE_MAX; !(i == 0); --i) { image_stat[i] -= image_stat[i - 1]; } // Watershed algorithm GVPoint *head = image_sort; GVInt32 space = 0; GVInt32 *dist = new GVInt32[size]; GVInt32 dist_cnt = 0; GVInt32 label_cnt = 0; std::queue<GVPoint> queue; ::memset(dist, 0, sizeof (GVInt32) * size); ::memset(label, ~0, sizeof (GVInt32) * size); for (GVInt32 h = 0; !(h == (GV_BYTE_MAX + 1)); ++h) { head += space; space = image_stat[h]; for (GVInt32 i = 0; !(i == space); ++i) { GVInt32 index = head[i].y * width + head[i].x; GVInt32 index_l = ((head[i].x - 1) < 0) ? -1 : ((head[i].x - 1) + head[i].y * width); GVInt32 index_r = !((head[i].x + 1) > width) ? -1 : ((head[i].x + 1) + head[i].y * width); GVInt32 index_t = ((head[i].y - 1) < 0) ? -1 : (head[i].x + (head[i].y - 1) * width); GVInt32 index_b = !((head[i].y + 1) > height) ? -1 : (head[i].x + (head[i].y + 1) * width); label[index] = MASK; if ( (!(index_l < 0) && !(label[index_l] < WSHD)) || (!(index_r < 0) && !(label[index_r] < WSHD)) || (!(index_t < 0) && !(label[index_t] < WSHD)) || (!(index_b < 0) && !(label[index_b] < WSHD))) { dist[index] = 1; queue.push(head[i]); } } dist_cnt = 1; queue.push(FICT_PIXEL); while (GV_TRUE) { GVPoint top = queue.front(); GVInt32 index = top.y * width + top.x; GVInt32 index_l = ((top.x - 1) < 0) ? -1 : ((top.x - 1) + top.y * width); GVInt32 index_r = !((top.x + 1) > width) ? -1 : ((top.x + 1) + top.y * width); GVInt32 index_t = ((top.y - 1) < 0) ? -1 : (top.x + (top.y - 1) * width); GVInt32 index_b = !((top.y + 1) > height) ? -1 : (top.x + (top.y + 1) * width); queue.pop(); if (top == FICT_PIXEL) { if (queue.empty()) break; else { ++dist_cnt; queue.push(FICT_PIXEL); top = queue.front(); queue.pop(); } } if (!(index_l < 0)) { if ((dist[index_l] < dist_cnt) && !(label[index_l] < WSHD)) { if (label[index_l] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_l]; else if (!(label[index] == label[index_l])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_l] == MASK) && (dist[index_l] == 0)) { dist[index_l] = dist_cnt + 1; queue.push(GVPoint(top.x - 1, top.y)); } } if (!(index_r < 0)) { if ((dist[index_r] < dist_cnt) && !(label[index_r] < WSHD)) { if (label[index_r] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_r]; else if (!(label[index] == label[index_r])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_r] == MASK) && (dist[index_r] == 0)) { dist[index_r] = dist_cnt + 1; queue.push(GVPoint(top.x + 1, top.y)); } } if (!(index_t < 0)) { if ((dist[index_t] < dist_cnt) && !(label[index_t] < WSHD)) { if (label[index_t] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_t]; else if (!(label[index] == label[index_t])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_t] == MASK) && (dist[index_t] == 0)) { dist[index_t] = dist_cnt + 1; queue.push(GVPoint(top.x, top.y - 1)); } } if (!(index_b < 0)) { if ((dist[index_b] < dist_cnt) && !(label[index_b] < WSHD)) { if (label[index_b] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_b]; else if (!(label[index] == label[index_b])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_b] == MASK) && (dist[index_b] == 0)) { dist[index_b] = dist_cnt + 1; queue.push(GVPoint(top.x, top.y + 1)); } } } for (GVInt32 i = 0; !(i == space); ++i) { GVInt32 index = head[i].y * width + head[i].x; dist[index] = 0; if (label[index] == MASK) { label_cnt++; label[index] = label_cnt; queue.push(head[i]); while (!queue.empty()) { GVPoint top = queue.front(); GVInt32 index_l = ((top.x - 1) < 0) ? -1 : ((top.x - 1) + top.y * width); GVInt32 index_r = !((top.x + 1) > width) ? -1 : ((top.x + 1) + top.y * width); GVInt32 index_t = ((top.y - 1) < 0) ? -1 : (top.x + (top.y - 1) * width); GVInt32 index_b = !((top.y + 1) > height) ? -1 : (top.x + (top.y + 1) * width); queue.pop(); if (!(index_l < 0) && (label[index_l] == MASK)) { queue.push(GVPoint(top.x - 1, top.y)); label[index_l] = label_cnt; } if (!(index_r < 0) && (label[index_r] == MASK)) { queue.push(GVPoint(top.x + 1, top.y)); label[index_r] = label_cnt; } if (!(index_t < 0) && (label[index_t] == MASK)) { queue.push(GVPoint(top.x, top.y - 1)); label[index_t] = label_cnt; } if (!(index_b < 0) && (label[index_b] == MASK)) { queue.push(GVPoint(top.x, top.y + 1)); label[index_b] = label_cnt; } } } } } // Release resources delete[] image_stat; delete[] image_space; delete[] image_sort; delete[] dist; } |
以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持服务器之家。
原文链接：http://blog.csdn.net/nagao_kagetora/article/details/5995925
