【数字图像处理】7.10:灰度图像-图像分割 区域分割之分水岭算法

Abstract: 数字图像处理:第59天
Keywords: 分水岭算法

本文最初发表于csdn,于2018年2月17日迁移至此

灰度图像-图像分割 区域分割之分水岭算法

今天已经是第60篇博客了,这六十篇每一篇平均要两天左右,所以,在过去的四个月学到了这么多知识,想想挺开心,但学的越多就会发现自己不会的越多。从小学到大学,这么多年一直以学习为主要工作但学习又有很多阶段,对于通用知识,比如小学的语文数学此观点不适用,对于一些专业性较强的知识,感觉会有两个很主要的阶段,感觉自己目前处于入门阶段,由于数字图像涉及数学的知识较多,还有信号和信息论的知识,所以,感觉还是比较难学科,不然招聘公司也不会一年几十万的养着图像处理工程师。
下面的图是本人的一点点见解,只是自己总结的,没有实践,也没有科学依据,不喜勿喷:

我感觉图像处理能分成三个阶段,或者更多,第一阶段的人很多,听说这行前景好或者工资高的人,多半会学习点图像的知识,比如彩色空间啊,了解下OpenCV啊等等,还有一些属于纯属上课被逼无奈的,比如我们学校就对电子信息类专业和通信类开数字图像处理这门课,而且我当时考试还挂了。。。。这部分大家会接触一些名词,一些简单算法,有用心的同学可能会实现下代码,这阶段风景不错,而且好多名词可以拿来忽悠HR或者忽悠导师,得到一份不错的工作,或者做点导师的项目。
这个阶段多半使用现成的函数库,了解了基本算法或者听别人说一些算法,自己来跑结果,这个阶段Matlab的用户量较大。
第1阶段后半期也就是平台期,这个阶段是做了一段时间有一定算法使用基础和代码能力的工程师,多半在各企业负责最底层的图像算法编写,在他们面前是一座大山,和一个路口,继续做图像还是转管理。
如果继续选择图像,就会面临一个峭壁,具体是算法底层的数学,说的数学,总是困难的,爬这个峭壁的动力可以是挣更多的钱,或者爱好,因为一旦到达阶段2,就能成为首席图像处理工程师,到任何需要图像处理的公司,都能够独当一面,这些人已经到了不缺钱的地步,而且在行业内一定有一定名气。
第二阶段的一旦到达,可以说是事业的平稳期,或者巅峰,不缺钱,还能独自决定一些技术层面上的事,指挥手下阶段1的员工工作,这个阶段面临的也是一个选择,就是靠这个吃一辈子饭,绝对没问题,再有就是向更高的境界冲击。
第三阶段没有尽头,因为能促使进入阶段三的动力只有爱好,这部分挣的钱可能没有阶段2挣的多,而且难度更大,看不到尽头,所以这部分属于探索阶段,在这个阶段上看到的一些,都能推动未来学科的发展,所以这个阶段的人多半是在实验室和数学中忙碌一生,然后几十年后出现在各大论文和教材中。
以上属于个人猜想或者愚见,想法幼稚,仅供参考。

算法描述

今天废话太多了,哈哈,其实上面的可以新开一篇博客单独写出来,但觉得,首先自己年轻视野狭窄,第二水平太低,所以当做废话夹在本文中。
今天说分水岭算法,分水岭算法族是一个思想引出的一些列实现方法。
算法内容:首先将二维灰度图像抽象为三维地形,横纵坐标表示经纬度,灰度表示高度,这样就产生了一个地形图,地形中有高山有盆地,我们的任务就是找到盆地并且给盆地划定地盘。
算法过程,将最低点(灰度)作为起始点,开始从下注水,想泉水一样,水位逐渐上升,所有点都是上下透水的,但不会横向流动,每产生一个新的积水盆地时标记这个盆地使之与原有盆地区别开,当两个盆地内的水要汇集的时候,在此点建立一个水坝,将水分开,此水坝无限高,水永远不会溢出,依次建立所有水坝,最后的水坝就是所谓的分水岭或分水线,也就是得到的分割线。
算法步骤:

  1. 检测图像中所有的盆地,并依次标记,(此时最终的分割区域数已确定)
  2. 将盆地周围邻域的像素按照灰度值大小入队(队列,数据结构,先进先出)
  3. 按照灰度值顺序将像素出队,进行判断
    判断1:如果像素周围只存在一个盆地标记,将此像素划分到此盆地
    判断2:如果像素周围存在两个以及以上盆地,次像素为分水岭。
  4. 如果此像素是盆地像素,将该像素邻域像素入队。
  5. 重复步骤3,直到队列为空。
  6. 得出结果

此算法描述只是分水岭算法之一,冈萨雷斯的书中使用的是形态学的描述方式,但与上述过程原理一致,感兴趣者可自行查阅。
解释下第一步,第一步是本算法的核心,此步骤首先要确定盆地还是非盆地,可以抽象成下面几种模型:

图像中的红色框内为非初始盆地,或叫做盆底,绿色为盆底,对于第二种盆地的确定需要使用图的深度或者广度优先搜索,判断所有候选点,来判断是第二种模型还是第三种。

代码

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
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
//
// Watershed
// tony.sheng.tan@gmail.com
// Created by 谭升 on 15/03/03.
// Copyright (c) 2015年 谭升. All rights reserved.
//
/*
typedef struct PriQueueNode_ PriQueueHead;
typedef struct NLevelPriQueueNode_ * NLevelPriQueue;
typedef struct NLevelPriQueueNode_ NLevelPriNode;
typedef struct ExPix_ Pix_Label;

* ___ ____________________
* | P |->| NLevelPriQueue |
* |---| |--------------------|
* | r |->| NLevelPriQueue |
* |---| |--------------------|
* | i |->| NLevelPriQueue |
* |---| |--------------------|
* | Q |->| NLevelPriQueue |
* |---| |--------------------|
* | u |->| NLevelPriQueue |
* |---| |--------------------|
* | e |->| NLevelPriQueue |
* |---| |--------------------|
* | u |->| NLevelPriQueue |
* |---| |--------------------|
* | e |->| NLevelPriQueue |
* |___| |____________________|

struct NLevelPriNode_{
int x;
int y;
NLevelPriQueue next;
};
struct PriQueueNode_{
int nodeNum;
NLevelPriQueue head;
NLevelPriQueue tail;
};
struct ExPix_{
int grayvalue;
int label;
};
*/
#include "watershed.h"
#include <stdio.h>

//将非极小值且与极小值相邻的元素入队
void inQueue(PriQueueHead* priQueue,int gray_level,int x,int y){
priQueue[gray_level].nodeNum++;
///malloc new node
NLevelPriNode* newNode=(NLevelPriNode *)malloc(sizeof(NLevelPriNode));
newNode->x=x;
newNode->y=y;
newNode->next=NULL;
if(priQueue[gray_level].head==NULL){
priQueue[gray_level].head=newNode;
priQueue[gray_level].tail=newNode;
}else{
priQueue[gray_level].tail->next=newNode;
priQueue[gray_level].tail=newNode;
}

}
//判断极小值是平底锅型还是台阶形状,使用图的深度优先搜索
int isPan(int *src,int width,int height,double value,int x,int y){
src[y*width+x]=-value;

for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++)
if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){
if(src[(j+y)*width+i+x]<value&&src[(j+y)*width+i+x]>0)
return 0;
}
for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++)
if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){
if(src[(j+y)*width+i+x]==value){
return(isPan(src,width, height, value, i+x, j+y));

}
}
return 1;

}
//由于判断平底锅时使部分数据损坏,现进行恢复
void repairPan(int *src ,int width,int height,double value,int x,int y){
src[y*width+x]=-value;
for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++)
if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){
if(src[(j+y)*width+i+x]==value){
repairPan(src, width, height, value, x+i, y+j);
}
}
}
//平底锅形状,标记极小值为255
void setMinimal(int *src,double *dst,int width,int height,int value,int x,int y){
dst[y*width+x]=255.0;
for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++)
if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){
if(src[(j+y)*width+x+i]==value&&dst[(j+y)*width+x+i]==0.0)
setMinimal(src,dst,width, height, value, x+i,y+j);
}
}
//台阶形状,标记非极小值127
void setUnMinimal(int *src,double *dst,int width,int height,int value,int x,int y){
dst[y*width+x]=127.0;
for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++)
if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){
if(src[(j+y)*width+x+i]==value&&dst[(j+y)*width+x+i]==0.0)
setUnMinimal(src,dst,width, height, value, x+i,y+j);
}
}

//数据类型从double到int
void Double2Int(double *src,int* dst,int width,int height){
for(int i=0;i<width*height;i++)
dst[i]=(int)src[i]+1;
}
//寻找极小值,包括单点极小值和平底锅
void findMinimal(double *src,double *dst,int width,int height){

int *temp=(int *)malloc(sizeof(int)*width*height);
double *dsttemp=(double *)malloc(sizeof(double)*width*height);
Zero(dsttemp, width, height);
Double2Int(src, temp, width, height);
int lessthan=0;
int equ=0;
double min=findMatrixMin(src, width, height);
for(int i=0;i<width*height;i++)
if(src[i]==min)
dsttemp[i]=255.0;
for(int j=0;j<height;j++){
for(int i=0;i<width;i++){
lessthan=0;
equ=0;
int pix=temp[j*width+i];
if(dsttemp[j*width+i]==0.0){
for(int m=-1;m<2;m++)
for(int n=-1;n<2;n++)
if(j+m>=0&&i+n>=0&&j+m<height&&i+n<width){
if(m!=0||n!=0){
if(temp[(j+m)*width+i+n]<pix)
lessthan=1;
if(temp[(j+m)*width+i+n]==pix)
equ=1;
}
}

if(equ==1&&lessthan==0){
if(isPan(temp, width, height,pix, i, j)){
//repairPan(temp,width, height, -pix, i,j);
setMinimal(temp,dsttemp, width, height, -pix, i, j);
}else {
repairPan(temp,width, height, -pix, i,j);
setUnMinimal(temp,dsttemp, width, height, pix, i, j);
}
}
if(lessthan==1)
dsttemp[j*width+i]=127.0;
if(0==lessthan&&0==equ)
dsttemp[j*width+i]=255.0;

}
}
}
matrixCopy(dsttemp, dst, width, height);
free(dsttemp);
free(temp);
}
//标记极小值的label,从-1开始向下增长 -2 -3 -4 -5 -6 -7.....
void LabelMinimal(double * src,Pix_Label* dst,int width,int height,int x,int y,int label){
dst[y*width+x].label=label;
for(int i=-1;i<2;i++){
for(int j=-1;j<2;j++)
if(x+i>=0&&x+i<width&&y+j>=0&&y+j<height&&(i!=0||j!=0))
if(src[(y+j)*width+x+i]==255.0&&dst[(y+j)*width+x+i].label==0){
LabelMinimal(src, dst, width, height, x+i, y+j, label);
}
}

}
//初始化label数组,此数组与图像数组多加了label
void InitLabelMat(double *src,Pix_Label* dst,double *mask,int width,int height){
for(int i=0;i<width*height;i++){
dst[i].grayvalue=src[i];
dst[i].label=0;
}
int label_minimal=-1;
for(int j=0;j<height;j++)
for(int i=0;i<width;i++){
if(mask[j*width+i]==255.0&&dst[j*width+i].label==0){
LabelMinimal(mask, dst, width,height,i,j, label_minimal);
label_minimal--;
}
}
}
//初始化队列头数组
void InitPriQueue(PriQueueHead* priQueue,Pix_Label *srclabel,int width,int height){
for(int i=0;i<GRAY_LEVEL;i++){
priQueue[i].head=NULL;
priQueue[i].nodeNum=0;
priQueue[i].tail=NULL;
}
int inqueue=0;
for(int j=0;j<height;j++)
for(int i=0;i<width;i++){
inqueue=0;
if(srclabel[j*width+i].label==0){
for(int m=-1;m<2;m++)
for(int n=-1;n<2;n++){
if(m+j>=0&&m+j<height&&n+i>=0&&n+i<width&&(m!=0||n!=0))
if(srclabel[(m+j)*width+n+i].label<0)
inqueue=1;
}
if(inqueue){
inQueue(priQueue,srclabel[j*width+i].grayvalue,i,j);
srclabel[j*width+i].label=INQUEUE;
}
}
}
//for(int i=0;i<GRAY_LEVEL;i++)
// printf("g:%d NumofNode:%d\n",i,priQueue[i].nodeNum);

}


int PirQueueisEmpty(PriQueueHead* priqueue){
int sum=0;
for(int i=0;i<GRAY_LEVEL;i++)
sum+=priqueue[i].nodeNum;
return !sum;
}


NLevelPriNode* outQueue(PriQueueHead* priqueue){
NLevelPriNode* node=NULL;
if(!PirQueueisEmpty(priqueue))
for(int i=0;i<GRAY_LEVEL;i++)
if(priqueue[i].nodeNum!=0){
node=priqueue[i].head;
priqueue[i].head=node->next;
priqueue[i].nodeNum--;
break;
}
return node;
}


void findWaterShed(Pix_Label * srclabel,PriQueueHead* priqueue,int width,int height){
NLevelPriNode* node=outQueue(priqueue);
while(node!=NULL){
int y=node->y;
int x=node->x;
//printf("x:%d y:%d \n",x,y);
int label=0;
int isWatershed=0;
for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++){
if(j+y>=0&&j+y<height&&i+x>=0&&i+x<width&&(j!=0||i!=0)){
if(srclabel[(j+y)*width+x+i].label<0){
if(label==0)
label=srclabel[(j+y)*width+x+i].label;
else if(label!=srclabel[(j+y)*width+x+i].label){
isWatershed=1;
}
}
}
}
if(isWatershed)
srclabel[y*width+x].label=WATERSHED;
else if(label<0){
srclabel[y*width+x].label=label;
for(int j=-1;j<2;j++)
for(int i=-1;i<2;i++){
if(j+y>=0&&j+y<height&&i+x>=0&&i+x<width&&(j!=0||i!=0)){
if(srclabel[(j+y)*width+x+i].label==0){
inQueue(priqueue, srclabel[(j+y)*width+i+x].grayvalue, i+x, j+y);
srclabel[(j+y)*width+i+x].label=INQUEUE;
}
}
}
}
else if(label==0&&isWatershed==0){
srclabel[y*width+x].label=0;

}
free(node);
node=outQueue(priqueue);

}

}
//meyer分水岭方法
void MeyerWatershed(double *src,double *dst,int width,int height){
double *dst_temp=(double *)malloc(sizeof(double)*width*height);
Zero(dst_temp, width, height);
Pix_Label * srclabel=(Pix_Label*)malloc(sizeof(Pix_Label)*width*height);
PriQueueHead priqueue[GRAY_LEVEL];
double *minimal=(double *)malloc(sizeof(double)*width*height);
Zero(minimal, width, height);
findMinimal(src, minimal, width, height);
InitLabelMat(src, srclabel, minimal, width, height);
//for(int j=0;j<height;j++){
// for(int i=0;i<width;i++)
// printf(" l:%3d|",srclabel[j*width+i].label);
// printf("\n");
//}
InitPriQueue(priqueue, srclabel, width, height);
/*for(int i=0;i<GRAY_LEVEL;i++){
NLevelPriNode *node=priqueue[i].head;
printf("%d:",i);
while (node!=NULL) {
printf("x:%d,y:%d ",node->x,node->y);
node=node->next;
}
printf("\n");

}*/
findWaterShed(srclabel, priqueue, width, height);
for(int i=0;i<width*height;i++)
if(srclabel[i].label==WATERSHED)
dst_temp[i]=255.0;
matrixCopy(dst_temp, dst, width, height);
free(dst_temp);
free(srclabel);
free(minimal);

}

实现效果

总结

以上实现目前还有点缺陷,就是盆底面积很大的时候,递归计算会crash,初步设想的解决办法是使用迭代替换递归,分割结果相对来讲对于前面的算法来说相对较好,但运算速度较慢,对噪声敏感。
灰度图像的内容简单的介绍至此,下篇开始介绍彩色基础知识和算法,欢迎收看。
待续。。。

原文地址1:https://www.face2ai.com/DIP-7-10-灰度图像-图像分割-区域分割之分水岭算法转载请标明出处

0%