OpenCV提取最大联通面积--铁轨轨道特征(二)

挺突然的。

讲道理我这个算法写的并不好 比官方的findcounters差远了。。。

关联项目:
http://www.h13studio.com/基于Linux嵌入式开发板的轨道检测机器人-大创项目-持续更新中
http://www.h13studio.com/OpenCV在Arm-Linux上的安装
http://www.h13studio.com/OpenCV提取铁轨轨道特征-一

基本原理/目的

找出一张二值图的联通部分以及求出对应面积

也就是说找到下图中若干条 曲线 所构成的白色联通区域 并分别求出他们的面积。

当然 图片中只有两条 但是实际中却不是。
不过,我做这个的目的只是为了提取线段,假设线段宽度恒定,那么面积即可代表线段长度。

也许用”面积”来求”长度”不严谨
但是因为OpenCV有 Canny 等边缘提取算法的存在,我们可以很方便的把一条有宽度的线段计算生成为一条宽度为1或者宽度相近的线段。
为了美观,后面的图片均用黑色代表有效数据。

逻辑思考

先认真的思考再开始编程 比先编程再思考要好得多 – 沃·滋基索德

要达到这个目的,我们就肯定需要遍历这幅图中所有的像素点。

方便起见,这里先按照先上后下 然后先左后右 也就是先rows再cols的遍历扫描
也就是先红后蓝

代码上就直接两个for解决。

1
2
for (int temprows = 0; temprows < input.rows; temprows++)
for (int tempcols = 0; tempcols < input.cols; tempcols++)

不过方向无所谓,大不了提前旋转镜像处理一下。

大体思路无非就两种
1.先计算联通面积 再提取线条,分成两步计算,这样效率高,但是用起来有一点点门槛
2.计算联通面积的同时提取线条,然后将线条中的点存入vector中,然后线条再组合成一个vector。这个和OpenCV官方的 findContours 有点像,但是我写出来之后时间消耗120ms,性能很差,所以我

累加部分

思路

计算联通面积其实很简单,也就是累加
比如这个线头

在二值化的 Mat 对象中可能是这样

然后我们就可以直接按照上文中的

1
方便起见,这里先按照先上后下 然后先左后右 也就是先rows再cols的遍历扫描

顺序扫描累加
很简单地就成了这样

这样一来,在线段宽度基本一致的情况下,就可以直接获取各个线段的起始点和终点,以及线段面积。
只要线段长短差距足够大,就可以很简单的筛选和提取线段
比如下图中的黄线和蓝线就很容易区分

那么 原理确定了,开始写程序吧。

这里的程序比较要求严谨性,很容易诱发各种bug,所以要兼顾下列所有情况。

第一种情况

也就是最简单的情况:
当我们扫描到红色点时,直接读取上一行的这个点的位置是否有有效像素即可判断是否为新线段起点。
然后,还需要遍历上一行的最大值来进行累加。

第二种情况

因为扫描顺序,在我们判断这个点是不是新点时,需要遍历本行临近有效像素点的上一行是否有接触点。

第三种情况

这个点是新点,需要从零赋值

第四种情况

两条线相交,或因扫描顺序造成的汇合。

这种情况其实我没有写,因为我用不到 →_→

第五种情况

一条线分叉为两条
直接单独继续累加就OK。

第六种情况

一条线分叉为两条后又汇合

….

第七种情况

因为图片质量或图像预处理造成的断线。
因为有 膨胀 算法等的存在,外加阈值不好判定,故不写。

第****种情况…

QNMD,自己写去…

代码实现

Windows版本:
(因为用到了微秒级计时,所以把计时部分替换或者删除就可以在Linux上跑)

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
#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui_c.h>
#include <vector>

#include <windows.h>

int main()
{
using namespace cv;
Mat input;
//input = imread("C:/Users/h13/Desktop/test2.png");
input = imread("C:/Users/h13/Desktop/pathway.jpg");
//input = imread("C:/Users/h13/Desktop/lightrail.jpg");
//imshow("input", input);

//开始计算时间
clock_t tstart, tsobel0, tthreshold0, tcvtColor0, tline0, tcontours0, tdilate0, tend0;//用于分别计算各个环节的时间消耗。


//求导图像
Mat sobel0,sobel1;
Sobel(input, sobel0, CV_8UC1, 1, 0); //Opencv的Sobel函数
Sobel(input, sobel1, CV_8UC1, 0, 1); //Opencv的Sobel函数
tsobel0 = clock();
sobel0 += sobel1;
//imshow("sobel0", sobel0); //性能测试时请注释本行代码。

//滤掉低亮度部分
Mat thershold0;
threshold(sobel0, thershold0, 200, 255, THRESH_TOZERO);
tthreshold0 = clock();
//imshow("thershold0", thershold0);//性能测试时请注释本行代码。

//转灰度图
Mat gray0;
cvtColor(thershold0, gray0, CV_BGR2GRAY);
tcvtColor0 = clock();
//imshow("gray0", gray0); //性能测试时请注释本行代码。

Mat thershold1;
threshold(gray0, thershold1, 200, 255, CV_THRESH_BINARY);

////膨胀(加粗线条)。如非必要,就别开启
//Mat dilate0, structure_element3;
//structure_element3 = Mat::zeros(3, 3, CV_8UC1);

//structure_element3.at<uchar>(1, 0) = 1;
//structure_element3.at<uchar>(1, 1) = 1;
//structure_element3.at<uchar>(1, 2) = 1;

//dilate(thershold1, dilate0, structure_element3);

//std::cout << "start" << std::endl;

////tthreshold0 = clock();
input = thershold1;
//imshow("thershold1", input);//性能测试时请注释本行代码。


LARGE_INTEGER freq_;
QueryPerformanceFrequency(&freq_);
LARGE_INTEGER begin_time,end_time;
QueryPerformanceCounter(&begin_time);
//findlines测试,只能输入单通道图片,目前只支持二值线图。
Mat temp = Mat::zeros(input.rows,input.cols, CV_32SC1);
Mat output = Mat::zeros(input.rows, input.cols, CV_8UC1);
int32_t maxvalue = 0;
int maxcols = 0, maxrows = 0;
typedef struct maxpoint_t{
Point maxpoint;
int value;
};
std::vector<maxpoint_t> maxpoints;
/*****************************************/
//扫描方式:
// 先行后列
// 1>---------------
// 2>---------------
// 3>--->
/*****************************************/
for (int temprows = 0; temprows < input.rows; temprows++)
{
uchar* data = input.ptr<uchar>(temprows);
int32_t* tempdata = temp.ptr<int32_t>(temprows);
if (temprows)
{
//当不是第一行的时候,开始遍历此行:

//tempdata指向本行,用于生成累加结果,lasttempdata指向上一行的累加结果,用于计算累加值。
int32_t* lasttempdata = temp.ptr<int32_t>(temprows - 1);

//existnextpixel用于表示本行是否为线段端点,true表示本行为线段端点,需要记录端点坐标,false则什么都不做。
bool existnextpixel = false;

//nexttempdata指向输入数据的下一行,用于计算本行是否含有端点。
uchar* nexttempdata = input.ptr<uchar>(temprows + 1);

//当是最后一行时,直接判断其为端点
if (temprows == input.rows - 1)
existnextpixel = true;

for (int tempcols = 0; tempcols < input.cols; tempcols++)
{
if (data[tempcols])
{
/*****************************************/
//情况1:
// > ...???
// > ...
//其中 '.' 代表有效像素 ' ' 代表无效像素
//
//思路:
// 直接算即可。
/*****************************************/
if (lasttempdata[tempcols] > 0.1)
{
float maxpixel = 0;
int templastcols = tempcols;
int32_t temp0 = lasttempdata[templastcols];

//查找上一行累加结果最大值
while ((temp0) && (templastcols > 0))
{
if (temp0 > maxpixel)
maxpixel = temp0;
temp0 = lasttempdata[templastcols--];
}

while ((temp0) && (templastcols > 0))
{
if (temp0 > maxpixel)
maxpixel = temp0;
temp0 = lasttempdata[templastcols++];
}

//开始计算累加值,并判断是否为端点
tempdata[tempcols] = maxpixel + 1;
if ((!existnextpixel) && (nexttempdata[tempcols]))
existnextpixel = true;

while ((data[tempcols]) && (tempcols < input.cols - 1))
{
tempcols++;
tempdata[tempcols] = tempdata[tempcols - 1] + 1;
if ((!existnextpixel) && (nexttempdata[tempcols]))
existnextpixel = true;
}

//只有当下一行没有相连的pixel的时候才记录最大值
if (!existnextpixel)
{
maxpoint_t temppoint;
temppoint.maxpoint = Point(temprows,tempcols);
maxpoints.push_back(temppoint);

if (tempdata[tempcols] > maxvalue)
{
maxvalue = tempdata[tempcols];
maxcols = tempcols;
maxrows = temprows;
}
}
}


/*****************************************/
//情况2:
// > ???
// > ..????
//其中 '.' 代表有效像素 ' ' 代表无效像素
//
//思路:
// 算到第二行时,先找到最右像素再开始向左算。
/*****************************************/
else
{
int startcol = tempcols;
float maxpixel = 0;
//如果本行为最右行
if (tempcols == input.cols - 1)
{
tempdata[tempcols] = 1;
}
else {
//切换到最右边的有效像素点
while ((data[tempcols] > 0) && (tempcols < input.cols - 1))
{
tempcols++;
//std::cout << maxpixel << std::endl;
//先检测相接的最大值
if (lasttempdata[tempcols] > maxpixel)
maxpixel = lasttempdata[tempcols];
}
//再算不重合的
int temp0 = tempcols;
while ((lasttempdata[temp0]) && (temp0 < input.cols - 1))
{
temp0++;
if (lasttempdata[temp0] > maxpixel)
maxpixel = lasttempdata[temp0];
}
}

/*****************************************/
//情况2.1
// > ...
// > ....
//
//情况2.2
// >
// > ..
//
//思路:
// 既然 maxpixel 已经初始化为0,
// 那么下列写法对于这两种情况是等价的。
/*****************************************/

//开始进行累加运算
tempdata[tempcols] = maxpixel + 1;
if ((!existnextpixel) && (nexttempdata[tempcols]))
existnextpixel = true;

int temp0 = tempcols;
while (temp0 > startcol)
{
temp0--;
if ((!existnextpixel) && (nexttempdata[temp0]))
existnextpixel = true;
tempdata[temp0] = tempdata[temp0 + 1] + 1;
}

//只有当下一行没有相连的pixel的时候才记录最大值
if (!existnextpixel)
{
maxpoint_t temppoint;
temppoint.maxpoint = Point(temprows, tempcols);
maxpoints.push_back(temppoint);

if (tempdata[tempcols] > maxvalue)
{
maxvalue = tempdata[tempcols];
maxcols = tempcols;
maxrows = temprows;
}
}

/*if (tempdata[startcol] > maxvalue)
{
maxvalue = tempdata[startcol];
maxcols = startcol;
maxrows = temprows;
}*/
}

}
}
}
else {
for (int tempcols = 0; tempcols < input.cols; tempcols++)
{
if (tempcols)
{
//当他不是第一行第一个时:
if (data[tempcols])//当这个点是有效点是:
tempdata[tempcols] = tempdata[tempcols - 1] + 1;
}
else
{
//当他是第一行第一个时:
if (data[tempcols])//当这个点是有效点是:
tempdata[tempcols] = 1;
}

if (tempdata[tempcols] > maxvalue)
{
maxvalue = tempdata[tempcols];
maxcols = tempcols;
maxrows = temprows;
}
}
}
}
QueryPerformanceCounter(&end_time);
double timecost = (end_time.QuadPart - begin_time.QuadPart) * 1000000.0 / freq_.QuadPart;

system("pause");

temp = temp*255 / maxvalue;
temp.convertTo(output, CV_8UC1);
imshow("findlines0", output);
//imshow("findlines0", temp);
std::cout << "findlines: " << std::endl;
std::cout << "maxvalue: " << maxvalue << std::endl;
std::cout << "maxpoint: " << maxcols << " " << maxrows << std::endl;
std::cout << "timecost: " << timecost << "ns" << std::endl;
std::cout << std::endl << "maxpoint: " << std::endl << " number: " << maxpoints.size() << std::endl;


//collectlines开始收集直线
//创建一个按照面积大小排序的引索
//应该有更快的算法 但是代价就是程序复杂度,所以还是直接反向算回去吧


waitKey(0);
}

性能:
平台:
i5-9300h
测试:

1
2
3
4
5
6
7
8
请按任意键继续. . .
findlines:
maxvalue: 72
maxpoint: 7 591
timecost: 1972.5ns

maxpoint:
number: 655

虽然才写一半 但是单线程不到2ms的成绩非常棒!!!

效果:

线条处理部分

然后就需要把累加的结果再返还给线条

TODO: 自定义扫描方向
暂时解决方案: 先旋转/镜像图片再输入
思路: 用1和-1表示方向 1表示正方向 //好像有点问题 先不慌

接下来,关于如何输出数据,有两种方案。
方案1:
先累加运算,用Mat对象保存累加值,然后再逆向

方案2:

1
2
3
4
findlines:
maxvalue: 2957
maxpoint: 639 0
timecost: 119355ns

目前思路: 从上到下从左到右扫描,扫描时就把线条数据存到vector内。

随便想了想 大概有这么多种可能以及对应解决方案。

分叉(一定要做):
永久分叉: -> 全都单独拷贝保存
保留全部/保留长线/保留短线/保留左终止点/保留右终止点/单独保存

分叉后汇合:
单独保存

多线分叉:
单独保存

汇合:
单线汇合:
汇合后终止 -> 合并为一条线/汇合后继续

回合后分离(也就是交叉):
TODO

断线:
断线连接:
距离判定:
1.平断 最左最右相减 断线垂直距离相减 阈值判断 how?
2.斜断 TODO

断线不连接:
简单的一批。

最后 考虑到实用性 我决定这样处理

1
2
3
4
5
6
7
8
9
10
11
12
13
从上到下从左到右扫描,扫描时就把线条数据存到vector内。

分叉(一定要做):
永久分叉: -> 全都单独拷贝保存

汇合:
单线汇合:
汇合后终止 -> 合并为一条线 -> maxvalue如何处理?

断线:
断线连接:
距离判定:
1.平断 最左最右相减 断线垂直距离相减 阈值判断 how?

逆向筛选部分

思路

代码实现

我打算分为两个版本,一个只累加,一个累加并筛出线段

后半段还没写完,等我两天 (咕咕咕)

中间断线是因为 Soble 之后 二值化没有处理好,然后传入函数的图片本来就断线了,想要修复这个问题很简单,等我把程序后半段写好就OK了。
至于横向的干扰,再简单不过,直接侵蚀算法干他就完事。