C++实现分水岭算法(Watershed Algorithm)

2025-05-27 0 62

分水岭分割方法(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;

}

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持快网idc。

收藏 (0) 打赏

感谢您的支持,我会继续努力的!

打开微信/支付宝扫一扫,即可进行扫码打赏哦,分享从这里开始,精彩与您同在
点赞 (0)

声明:本站所有文章,如无特殊说明或标注,均为本站原创发布。任何个人或组织,在未征得本站同意时,禁止复制、盗用、采集、发布本站内容到任何网站、书籍等各类媒体平台。如若本站内容侵犯了原著者的合法权益,可联系我们进行处理。

快网idc优惠网 建站教程 C++实现分水岭算法(Watershed Algorithm) https://www.kuaiidc.com/72850.html

相关文章

发表评论
暂无评论