二值图像连通域标记算法优化
文章概要
非常感谢☆ronny丶博主在其博文中对二值图像连通域的介绍和算法阐述，让我这个毫无数据结构算法底子的小白能够理解和复现代码。本文的目的是基于我自己的理解，对该博文中two-pass算法的一些优化和补充，同时也希望帮助更多像我一样的人较快地掌握连通域标记。
连通域标记是图像分割计数的一个重要环节，在工业上应用非常地多。例如像硬币的计件，在二值化处理后，为了能够感知数量，就得对硬币区域进行标记（当然标记前可能还要经过一系列的形态学处理）。另外，还有一个我想到的，更有趣、也更具有挑战性的例子——二维码连通域标记，这用来检验算法的性能是再合适不过了。言归正题——本文介绍了两大流行算法，一个是利用dfs的seed-filling算法，另一个是two-pass算法。后者因为处理等价对的方法不同，又细分为dfs two-pass（使用dfs处理等价对）和union-find two-pass（使用并查集处理等价对）。如果硬要给这三种算法排序的话，大概是union-find two-pass > seed-filling > dfs two-pass，反正我写的程序是这样的速度排序。
seed-filling算法
这个算法其实实质就是dfs，笔者曾经有幸做过一个“水洼连通”的算法题，当时就是用dfs或者bfs来做的，显然，“水洼连通”也是属于连通域标记问题的。dfs在这个问题上的思路是：优先地寻找一个完整连通域，在找的同时把他们都标记一下，找完一个完整连通域， 再去找下一个连通域。按照这个想法，程序无非就是维护一个堆栈或者队列罢了，写起来相对简洁易懂。要说缺点的话，就是频繁的堆栈操作可能会拉低程序的性能。
简要地说明一下这部分代码含义，故事就定义成小明踩水坑吧，虽然小明对我表示自己很文静，只喜欢做数学题。首先，定义了一个二维矩阵labels，大小跟二值图一样。一开始labels都是标签0，这是一个无效标签，可以理解为是充满迷雾的未知区域或者是已确定的非水坑区域。小明每到达一个新的单位域（也就是一个新像素），首先要先看看这个域是不是未曾踩过的水坑（未曾踩过的水坑其标签为0且灰度值为255），如果是的话，那么小明就原地开心地踩水坑了，踩过以后还不忘给它画上一个大于0的标记（以标签1为例）。接下来，小明回顾四周，又发现了接壤的另一个水坑， 于是又在该水坑上留下了标记1······这样看似单调的循环，在小明眼里却是一次次奇妙的冒险。愉快的时光很短暂，小明不一会儿就发现身边已经没有“新鲜”的水坑了，伤心的同时回到最初的那个水坑，继续朝远方走去。渐渐地，眼前依稀出现了陌生又熟悉的水坑，重现微笑的小明决定要开启新的旅途，因此标记1.0进化至2.0。
故事的结束，要额外补充一点，程序里要不停地将新的单位域加入队列， 因此队列遍历其上限是动态的。
1 vector<vector<int>> seedfilling(mat src)
2 {
3
4 // 标签容器，初始化为标记0
5 vector<vector<int>> labels(src.rows, vector<int>(src.cols, 0));
6 // 当前的种子标签
7 int curlabel = 1;
8 // 四连通位置偏移
9 pair<int, int> offset[4] = {make_pair(0, 1), make_pair(1, 0), make_pair(-1, 0), make_pair(0, -1)};
10 // 当前连通域中的单位域队列
11 vector<pair<int, int>> templist;
12
13 for (int i = 0; i < src.rows; i++)
14 {
15 for (int j = 0; j < src.cols; j++)
16 {
17 // 当前单位域已被标记或者属于背景区域, 则跳过
18 if (labels[i][j] != 0 || src.at<uchar>(i, j) == 0)
19 {
20 continue;
21 }
22 // 当前单位域未标记并且属于前景区域, 用种子为其标记
23 labels[i][j] = curlabel;
24 // 加入单位域队列
25 templist.push_back(make_pair(i, j));
26
27 // 遍历单位域队列
28 for (int k = 0; k < templist.size(); k++)
29 {
30 // 四连通范围内检查未标记的前景单位域
31 for (int m = 0; m < 4; m++)
32 {
33 int row = offset[m].first + templist[k].first;
34 int col = offset[m].second + templist[k].second;
35 // 防止坐标溢出图像边界
36 row = (row < 0) ? 0: ((row >= src.rows) ? (src.rows - 1): row);
37 col = (col < 0) ? 0: ((col >= src.cols) ? (src.cols - 1): col);
38
39 // 邻近单位域未标记并且属于前景区域, 标记并加入队列
40 if (labels[row][col] == 0 && src.at<uchar>(row, col) == 255)
41 {
42 labels[row][col] = curlabel;
43 templist.push_back(make_pair(row, col));
44 }
45 }
46 }
47 // 一个完整连通域查找完毕，标签更新
48 curlabel++;
49 // 清空队列
50 templist.clear();
51 }
52 }
53
54 return labels;
55 }
two-pass算法
等价对生成
关于two-pass的算法原理可以参考上面提到的博文，原文还是很详细的，唯一的遗憾就是后面程序的注释有点少，看起来会吃力些，说白了就是自己菜。要找一张二维图像中的连通域，很容易想到可以一行一行先把子区域找出来，然后再拼合成一个完整的连通域，因为从每一行找连通域是一件很简单的事。这个过程中需要记录每一个子区域，为了满足定位要求，并且节省内存，我们需要记录子区域所在的行号、区域开始的位置、结束的位置，当然还有一个表征子区域总数的变量。需要注意的就是子区域开始位置和结束位置在行首和行末的情况要单独拿出来考虑。
1 // 查找每一行的子区域
2 // numberofarea：子区域总数 starea：子区域开始位置 enarea：子区域结束位置 rowarea：子区域所在行号
3 void searcharea(const mat src, int &numberofarea, vector<int> &starea, vector<int> &enarea, vector<int> &rowarea)
4 {
5 for (int row = 0; row < src.rows; row++)
6 {
7 // 行指针
8 const uchar *rowdata = src.ptr<uchar>(row);
9
10 // 判断行首是否是子区域的开始点
11 if (rowdata[0] == 255)
12 {
13 numberofarea++;
14 starea.push_back(0);
15 }
16
17 for (int col = 1; col < src.cols; col++)
18 {
19 // 子区域开始位置的判断：前像素为背景，当前像素是前景
20 if (rowdata[col - 1] == 0 && rowdata[col] == 255)
21 {
22 // 在开始位置更新区域总数、开始位置vector
23 numberofarea++;
24 starea.push_back(col);
25 // 子区域结束位置的判断：前像素是前景，当前像素是背景
26 }else if (rowdata[col - 1] == 255 && rowdata[col] == 0)
27 {
28 // 更新结束位置vector、行号vector
29 enarea.push_back(col - 1);
30 rowarea.push_back(row);
31 }
32 }
33 // 结束位置在行末
34 if (rowdata[src.cols - 1] == 255)
35 {
36 enarea.push_back(src.cols - 1);
37 rowarea.push_back(row);
38 }
39 }
40 }
另外一个比较棘手的问题，如何给这些子区域标号，使得同一个连通域有相同的标签值。我们给单独每一行的子区域标号区分是很容易的事， 关键是处理相邻行间的子区域关系（怎么判别两个子区域是连通的）。
主要思路：以四连通为例，在上图我们可以看出be是属于同一个连通域，判断的依据是e的开始位置小于b的结束位置，并且e的结束地址大于b的开始地址；同理可以判断出ec属于同一个连通域，cf属于同一个连通域，因此可以推知becf都属于同一个连通域。
迭代策略：寻找e的相连区域时，对前一行的abcd进行迭代，找到相连的有b和c，而d的开始地址已经大于了e的结束地址，此时就可以提前break掉，避免不必要的迭代操作；接下来迭代f的时候，由于有e留下来的基础，因此对上一行的迭代可以直接从c开始。另外，当前行之前的一行如果不存在子区域的话，那么当前行的所有子区域都可以直接赋新的标签，而不需要迭代上一行。
标签策略：以上图为例，遍历第一行，a、b、c、d会分别得到标签1、2、3、4。到了第二行，检测到e与b相连，之前e的标签还是初始值0，因此给e赋上b的标签2；之后再次检测到c和e相连，由于e已经有了标签2，而c的标签为3，则保持e和c标签不变，将（2,3）作为等价对进行保存。同理，检测到f和c相连，且f标签还是初始值0，则为f标上3。如此对所有的子区域进行标号，最终可以得到一个等价对的列表。
下面的代码实现了上述的过程。子区域用一维vector保存，没办法直接定位到某一行号的子区域，因此需要用currow来记录当前的行，用firstareaprev记录前一行的第一个子区域在vector中的位置，用lastareaprev记录前一行的最后一个子区域在vector中的位置。在换行的时候，就去更新刚刚说的3个变量，其中firstareaprev的更新依赖于当前行的第一个子区域位置，所以还得用firstareacur记录当前行的第一个子区域。
1 // 初步标签,获取等价对
2 // labelofarea:子区域标签值， equallabels：等价标签对 offset:0为四连通,1为8连通
3 void markarea(int numberofarea, vector<int> starea, vector<int> enarea, vector<int> rowarea, vector<int> &labelofarea, vector<pair<int, int>> &equallabels, int offset)
4 {
5 int label = 1;
6 // 当前所在行
7 int currow = 0;
8 // 当前行的第一个子区域位置索引
9 int firstareacur = 0;
10 // 前一行的第一个子区域位置索引
11 int firstareaprev = 0;
12 // 前一行的最后一个子区域位置索引
13 int lastareaprev = 0;
14
15 // 初始化标签都为0
16 labelofarea.assign(numberofarea, 0);
17
18 // 遍历所有子区域并标记
19 for (int i = 0; i < numberofarea; i++)
20 {
21 // 行切换时更新状态变量
22 if (currow != rowarea[i])
23 {
24 currow = rowarea[i];
25 firstareaprev = firstareacur;
26 lastareaprev = i - 1;
27 firstareacur = i;
28 }
29
30 // 相邻行不存在子区域
31 if (currow != rowarea[firstareaprev] + 1)
32 {
33 labelofarea[i] = label++;
34 continue;
35 }
36 // 对前一行进行迭代
37 for (int j = firstareaprev; j <= lastareaprev; j++)
38 {
39 // 判断是否相连
40 if (starea[i] <= enarea[j] + offset && enarea[i] >= starea[j] - offset)
41 {
42 if (labelofarea[i] == 0)
43 // 之前没有标记过
44 labelofarea[i] = labelofarea[j];
45 else if (labelofarea[i] != labelofarea[j])
46 // 之前已经被标记，保存等价对
47 equallabels.push_back(make_pair(labelofarea[i], labelofarea[j]));
48 }else if (enarea[i] < starea[j] - offset)
49 {
50 // 为当前行下一个子区域缩小上一行的迭代范围
51 firstareaprev = max(firstareaprev, j - 1);
52 break;
53 }
54 }
55 // 与上一行不存在相连
56 if (labelofarea[i] == 0)
57 {
58 labelofarea[i] = label++;
59 }
60 }
61 }
dfs two-pass算法
通过上面的努力，标记任务并没有做完，最核心的部分正是如何处理等价对。这里简单贴上原博主说的dsf方法，又是深搜啊。相比于直接dfs标记连通域，先找等价对再深搜减少了大量的堆栈操作，因为前者深度取决于连通域的大小，而后者是连通域数量，显然这两个数量级的差距挺大的。
原博主的想法是建立一个bool型等价对矩阵，用作深搜环境。具体做法是先获取最大的标签值maxlabel，然后生成一个$maxlabel*maxlabel$大小的二维矩阵，初始值为false；对于例如（1，3）这样的等价对，在矩阵的（0，2）和（2，0）处赋值true——要注意索引和标签值是相差1的。就这样把所有等价对都反映到矩阵上。
深搜的目的在于建立一个标签的重映射。例如4、5、8是等价的标签，都重映射到标签2。最后重映射的效果就是标签最小为1，且依次递增，没有缺失和等价。深搜在这里就是优先地寻找一列等价的标签，例如一口气把4、5、8都找出来，然后给他们映射到标签2。程序也维护了一个队列，当标签在矩阵上值为true，而且没有被映射过，就加入到队列。
当然不一定要建立一个二维等价矩阵，一般情况，等价对数量要比maxlabel来的小，所以也可以直接对等价对列表进行深搜，但无论采用怎样的深搜，其等价对处理的性能都不可能提高很多。
1 // 等价对处理，标签重映射
2 void replaceequalmark(vector<int> &labelofarea, vector<pair<int, int>> equallabels)
3 {
4 int maxlabel = *max_element(labelofarea.begin(), labelofarea.end());
5 // 等价标签矩阵，值为true表示这两个标签等价
6 vector<vector<bool>> eqtab(maxlabel, vector<bool>(maxlabel, false));
7 // 将等价对信息转移到矩阵上
8 vector<pair<int, int>>::iterator labpair;
9 for (labpair = equallabels.begin(); labpair != equallabels.end(); labpair++)
10 {
11 eqtab[labpair->first -1][labpair->second -1] = true;
12 eqtab[labpair->second -1][labpair->first -1] = true;
13 }
14 // 标签映射
15 vector<int> labelmap(maxlabel + 1, 0);
16 // 等价标签队列
17 vector<int> templist;
18 // 当前使用的标签
19 int curlabel = 1;
20
21 for (int i = 1; i <= maxlabel; i++)
22 {
23 // 如果该标签已被映射，直接跳过
24 if (labelmap[i] != 0)
25 {
26 continue;
27 }
28
29
30 labelmap[i] = curlabel;
31 templist.push_back(i);
32
33 for (int j = 0; j < templist.size(); j++)
34 {
35 // 在所有标签中寻找与当前标签等价的标签
36 for (int k = 1; k <= maxlabel; k++)
37 {
38 // 等价且未访问
39 if (eqtab[templist[j] - 1][k - 1] && labelmap[k] == 0)
40 {
41 labelmap[k] = curlabel;
42 templist.push_back(k);
43 }
44 }
45 }
46
47 curlabel++;
48 templist.clear();
49 }
50 // 根据映射修改标签
51 vector<int>::iterator label;
52 for (label = labelofarea.begin(); label != labelofarea.end(); label++)
53 {
54 *label = labelmap[*label];
55 }
56
57 }
union-find two-pass算法
如果读者看到了这里，真的要感谢一下您的耐心。two-pass算法的代码要比直接深搜来得多，用不好甚至性能还远不如深搜。原博主在文中提及了可以用稀疏矩阵来处理等价对，奈何我较为愚钝，读者可以自研之。
讲到等价对，实质是一种关系分类，因而联想到并查集。并查集方法在这个问题上显得非常合适，首先将等价对进行综合就是合并操作，标签重映射就是查询操作（并查集可以看做一种多对一映射）。并查集具体算法我就不唠叨了，毕竟不怎么打程序设计竞赛。不过，采用并查集的话，函数定义估计就满天飞了，这里我包装了一下，做成了类——实在是有点强迫症，其中等价对生成的函数方法跟上面的是一样的。
网上有一些代码也实现了这个算法，但是有的牺牲了秩优化，合并时让树指向较小的根，个人认为这样做太不值了。所以为了解决这个，我在并查集映射后，又用labelremap来进行第二次的映射，主要的步骤跟前面的差不多。
然后，自己跑了一下这代码，不算画图标记的时间，效率要比上面的快四五倍左右，实时性上肯定是绰绰有余了。
1 class areamark
2 {
3 public:
4 areamark(const mat src,int offset);
5 int getmarkedarea(vector<vector<int>> &area);
6 void getmarkedimage(mat &dst);
7
8 private:
9 mat src;
10 int offset;
11 int numberofarea=0;
12 vector<int> labelmap;
13 vector<int> labelrank;
14 vector<int> starea;
15 vector<int> enarea;
16 vector<int> rowarea;
17 vector<int> labelofarea;
18 vector<pair<int, int>> equallabels;
19
20 void markarea();
21 void searcharea();
22 void setinit(int n);
23 int findroot(int label);
24 void unite(int labela, int labelb);
25 void replaceequalmark();
26 };
27
28 // 构造函数
29 // imageinput:输入待标记二值图像 offsetinput:0为四连通，1为八连通
30 areamark::areamark(mat imageinput,int offsetinput)
31 {
32 src = imageinput;
33 offset = offsetinput;
34 }
35
36 // 获取颜色标记图片
37 void areamark::getmarkedimage(mat &dst)
38 {
39 mat img(src.rows, src.cols, cv_8uc3, cv_rgb(0, 0, 0));
40 cvtcolor(img, dst, cv_rgb2hsv);
41
42 int maxlabel = *max_element(labelofarea.begin(), labelofarea.end());
43 vector<uchar> hue;
44 for (int i = 1; i<= maxlabel; i++)
45 {
46 // 使用hsv模型生成可区分颜色
47 hue.push_back(uchar(180.0 * (i - 1) / (maxlabel + 1)));
48 }
49
50 for (int i = 0; i < numberofarea; i++)
51 {
52 for (int j = starea[i]; j <= enarea[i]; j++)
53 {
54 dst.at<vec3b>(rowarea[i], j)[0] = hue[labelofarea[i]];
55 dst.at<vec3b>(rowarea[i], j)[1] = 255;
56 dst.at<vec3b>(rowarea[i], j)[2] = 255;
57 }
58 }
59
60 cvtcolor(dst, dst, cv_hsv2bgr);
61 }
62
63 // 获取标记过的各行子区域
64 int areamark::getmarkedarea(vector<vector<int>> &area)
65 {
66 searcharea();
67 markarea();
68 replaceequalmark();
69 area.push_back(rowarea);
70 area.push_back(starea);
71 area.push_back(enarea);
72 area.push_back(labelofarea);
73 return numberofarea;
74 }
75
76 void areamark::searcharea()
77 {
78 for (int row = 0; row < src.rows; row++)
79 {
80 // 行指针
81 const uchar *rowdata = src.ptr<uchar>(row);
82
83 // 判断行首是否是子区域的开始点
84 if (rowdata[0] == 255)
85 {
86 numberofarea++;
87 starea.push_back(0);
88 }
89
90 for (int col = 1; col < src.cols; col++)
91 {
92 // 子区域开始位置的判断：前像素为背景，当前像素是前景
93 if (rowdata[col - 1] == 0 && rowdata[col] == 255)
94 {
95 // 在开始位置更新区域总数、开始位置vector
96 numberofarea++;
97 starea.push_back(col);
98 // 子区域结束位置的判断：前像素是前景，当前像素是背景
99 }else if (rowdata[col - 1] == 255 && rowdata[col] == 0)
100 {
101 // 更新结束位置vector、行号vector
102 enarea.push_back(col - 1);
103 rowarea.push_back(row);
104 }
105 }
106 // 结束位置在行末
107 if (rowdata[src.cols - 1] == 255)
108 {
109 enarea.push_back(src.cols - 1);
110 rowarea.push_back(row);
111 }
112 }
113 }
114
115
116
117 void areamark::markarea()
118 {
119 int label = 1;
120 // 当前所在行
121 int currow = 0;
122 // 当前行的第一个子区域位置索引
123 int firstareacur = 0;
124 // 前一行的第一个子区域位置索引
125 int firstareaprev = 0;
126 // 前一行的最后一个子区域位置索引
127 int lastareaprev = 0;
128
129 // 初始化标签都为0
130 labelofarea.assign(numberofarea, 0);
131
132 // 遍历所有子区域并标记
133 for (int i = 0; i < numberofarea; i++)
134 {
135 // 行切换时更新状态变量
136 if (currow != rowarea[i])
137 {
138 currow = rowarea[i];
139 firstareaprev = firstareacur;
140 lastareaprev = i - 1;
141 firstareacur = i;
142 }
143
144 // 相邻行不存在子区域
145 if (currow != rowarea[firstareaprev] + 1)
146 {
147 labelofarea[i] = label++;
148 continue;
149 }
150 // 对前一行进行迭代
151 for (int j = firstareaprev; j <= lastareaprev; j++)
152 {
153 // 判断是否相连
154 if (starea[i] <= enarea[j] + offset && enarea[i] >= starea[j] - offset)
155 {
156 if (labelofarea[i] == 0)
157 // 之前没有标记过
158 labelofarea[i] = labelofarea[j];
159 else if (labelofarea[i] != labelofarea[j])
160 // 之前已经被标记，保存等价对
161 equallabels.push_back(make_pair(labelofarea[i], labelofarea[j]));
162 }else if (enarea[i] < starea[j] - offset)
163 {
164 // 为当前行下一个子区域缩小上一行的迭代范围
165 firstareaprev = max(firstareaprev, j - 1);
166 break;
167 }
168 }
169 // 与上一行不存在相连
170 if (labelofarea[i] == 0)
171 {
172 labelofarea[i] = label++;
173 }
174 }
175 }
176
177
178 // 并查集初始化
179 void areamark::setinit(int n)
180 {
181 for (int i = 0; i <= n; i++)
182 {
183 labelmap.push_back(i);
184 labelrank.push_back(0);
185 }
186 }
187
188 // 查根
189 int areamark::findroot(int label)
190 {
191 if (labelmap[label] == label)
192 {
193 return label;
194 }
195 else
196 {
197 //路径压缩优化
198 return labelmap[label] = findroot(labelmap[label]);
199 }
200 }
201
202 // 合并
203 void areamark::unite(int labela, int labelb)
204 {
205 labela = findroot(labela);
206 labelb = findroot(labelb);
207
208 if (labela == labelb)
209 {
210 return;
211 }
212 // 秩优化，秩大的树合并秩小的树
213 if (labelrank[labela] < labelrank[labelb])
214 {
215 labelmap[labela] = labelb;
216 }
217 else
218 {
219 labelmap[labelb] = labela;
220 if (labelrank[labela] == labelrank[labelb])
221 {
222 labelrank[labela]++;
223 }
224 }
225
226 }
227
228 // 等价对处理，标签重映射
229 void areamark::replaceequalmark()
230 {
231 int maxlabel = *max_element(labelofarea.begin(), labelofarea.end());
232
233 setinit(maxlabel);
234
235 // 合并等价对，标签初映射
236 vector<pair<int, int>>::iterator labpair;
237 for (labpair = equallabels.begin(); labpair != equallabels.end(); labpair++)
238 {
239 unite(labpair->first, labpair->second);
240 }
241
242 // 标签重映射，填补缺失标签
243 int newlabel=0;
244 vector<int> labelremap(maxlabel + 1, 0);
245 vector<int>::iterator old;
246 for (old = labelmap.begin(); old != labelmap.end(); old++)
247 {
248 if (labelremap[findroot(*old)] == 0)
249 {
250 labelremap[findroot(*old)] = newlabel++;
251 }
252 }
253 // 根据重映射结果修改标签
254 vector<int>::iterator label;
255 for (label = labelofarea.begin(); label != labelofarea.end(); label++)
256 {
257 *label = labelremap[findroot(*label)];
258 }
259
260 }
最后的最后，这些代码都没有经历过“岁月的历练”，如果存在不合理之处，请读者指正！