探地雷达(GPR)检测地下钢筋，并可视化钢筋

先说一下，这篇博文是关于 地下钢筋串线 的，就是用 探地雷达（GPR）检测到地下的钢筋 ，成像是多张 Bscan图像，每张Bscan图像中都有许多钢筋的形状。将多张Bscan中属于同一条钢筋的串起来，并做一个可视化的效果。

---

前言

首先来看一下我们用探地雷达（GPR）采集到的地下病害数据。假设 多条钢筋在地下的朝向 与 探地雷达行进的方向一致 ，那么在不同位置，雷达就会呈现出不同的Bscan图像。如下图就是多条钢筋，不同部位的图像。

（注意：这里，我们需要通过网络模型去检测到每张Bscan的钢筋，并且得到钢筋的框坐标，然后才能进行后续工作。目标检测的模型很多，比如Faster-rcnn，YOLO等，本文只讲如何串线）

我现在要做的就是将多张Bscan图像中属于同一条的钢筋串起来，并且可视化。

---

一、问题分析，钢筋串线

我们分析一下上面的三幅Bscan图像，可以得知，属于同一条钢筋的话，那么它在各自Bscan图像上的位置大致是一样的，至少差别不大。所以，我们可以通过位置信息来判断它属于哪一条钢筋。可以在X和Y方向上都设定一个阈值，只要位置(坐标)之差在阈值内则属于同一条钢筋，否则不属于同一条钢筋。

（先不考虑钢筋条斜着的情况）。

二、钢筋串线步骤

1.找到每张Bscan中钢筋的顶点

为什么要找到每张Bscan中钢筋的顶点？第一，因为能更加方便的串线，也比较简单。第二，实现也比较容易，可以通过传统的图像处理方式实现。找到顶点之后，记录它在整张Bscan图像上的图像坐标（像素坐标）。

代码如下（示例）：

def get_top(src, w_th=3, h_th=3):
	#src代表钢筋图像
    ret, thresh1 = cv2.threshold(src, 180, 255, cv2.THRESH_BINARY)
    height, weight = np.shape(thresh1)
    h = np.zeros(height, dtype=np.int)
    w = np.zeros(weight, dtype=np.int)
    for i in range(height):
        for j in range(weight):
            if src[i][j] == 255:
                h[i] += 1
                w[j] += 1
    max = 0
    y = 0
    for i in range(height):
        if h[i] > h_th:
            y = i
            break
    for j in range(weight):
        if w[j]>max:
            max = w[j]
    x1 = x2 = 0
    for j in range(weight):
        if w[j] < max-w_th and x2 == 0:
            x1 += 1
        elif w[j] >= max-w_th:
            x2 += 1
    return [y, int(x1+x2/2)]  #返回顶点坐标

结果如下图所示：

2.根据顶点坐标，通过X和Y方向的阈值判断是否属于统一条钢筋

这一步比较简单，就是 比较前后Bscan图相应位置上顶点的X和Y方向的差值 ，差值小于阈值则认为是同一条钢筋，否则不属于同一条钢筋。这里，我把所有的顶点按照钢筋进行了区分，并且把区分的结果存在一个字典中，型如：

testpipline={
“0”:[ #第0条钢筋有哪些顶点
{‘class’: ‘00’, ‘id’: 0, ‘pointy’: 85, ‘pointx’: 606, ‘d’: 1, ‘z’: 0}, # 'z’代表属于哪一张Bscan图
{‘class’: ‘00’, ‘id’: 0, ‘pointy’: 83, ‘pointx’: 605, ‘d’: 1, ‘z’: 1}
]}

代码如下（示例）：

def getpipline(result_list_2d_all):
	# result_list_2d_all是一个字典，里面包含所有的顶点信息
    x_threhold = 10  #x方向阈值
    y_threhold = 20  #y方向阈值
    # 先把所有的PIP按照属于哪一张Bscan进行分类
    split_list_2d = splitlsit_2d(result_list_2d_all)
    list_len = len(split_list_2d)
    # 把第一个bscan的所有顶点赋给pipline，认为初始的钢筋条数等于第一张Bscan的Pip的总数
    init_pipline(split_list_2d)
    index = 1
    if list_len == 0:
        print("result_list_2d_all is empty!")
    if list_len == 1:
        return split_list_2d[0]
    else:
        while index < list_len:
            #从第二张Bscan开始，依次和前面所有PIP顶点进行匹配
            findpoints(split_list_2d[index], x_threhold, y_threhold)
            index += 1

    w = open('result_piplines2.txt', 'w')
    for kk in pipline:
        w.write('{}:\n'.format(kk))
        for jj in pipline[kk]:
            w.write('{},{},{}\n'.format(jj["pointx"],jj["pointy"],jj["z"]))
    w.close()

    return pipline

3.将像素坐标转化为实际地理坐标（重点）

这里说明一下，探地雷达(GPR)在行进的过程中，我们生成了一个txt文件，记录了采集时的实际地理位置（x,y,z=0）,z=0是因为GPR始终在地表面。如下图，探地雷达向地下发射电磁波，它在地表面有真实的（x，y，z=0）坐标，生成的Bscan图的钢筋顶点有（pointx，pointy）像素坐标。
因此，我们可以根据钢筋顶点的 pointx，在 txt 文件的第 pointx 行找到该点真实坐标（x,y）。至于该点的深度，通过公式：deep= pointy / 450 * 1.5 ，这个公式是根据GPR设置的。
（这里其实不难，只要明白GPR采集数据的过程就能懂了。）

4.可视化串线结果

通过上面三步，我们已经能得到所有Bscan图中，每个钢筋顶点的实际坐标，并且将这些顶点按照是否属于一条钢筋进行了分类。

#获取转换为实际坐标后的顶点，返回的一个字典类型
res=changePipCor()

figure = plt.figure()
axes = Axes3D(figure)  #使用了Axes3D进行可视化
for k in res:
    x=[]
    y=[]
    z=[]
    print("{}:".format(k))
    for l in res[k]:
        print("l:",l)
        x.append(float(l["pip_factx"]))
        y.append(float(l["pip_facty"]))
        z.append(float(l["pip_factz"]))
    x=np.array(x)
    y = np.array(y)
    z = np.array(z)
    cc=randomcolor()
    axes.scatter(x, y, z, color=cc) #可视化

plt.show()

结果如下图所示：

三、总结

总的来说：
1.我先通过深度学习模型检测出每一张Bscan图像里面的钢筋，并且返回目标框的坐标。
2.通过检测的目标框，找到每个钢筋顶点的像素坐标。
3.通过像素坐标，设定X和Y方向阈值，判断哪些顶点属于同一条钢筋，进行分类。
4.将分好类的像素坐标转化为实际坐标，（这里要了解GPR采集数据的过程）
5.根据实际坐标可视化

（注意：当然，我这种方法还没能解决钢筋在地下斜放着的情况，这可能涉及到聚类，优化等方法，还需要考虑到更多信息，如果您有更好的方法，欢迎交流。希望本文对您有帮助！）

---