具体实施方式
下面结合附图和具体实施方式对本发明的技术方案作进一步具体说明,由此,本发明的优点和特点将会随着描述而更为清楚。
本领域的技术人员能够理解,尽管以下的说明涉及到有关本发明的实施例的很多技术细节,但这仅为用来说明本发明的原理的示例、而不意味着任何限制。本发明能够适用于不同于以下例举的技术细节之外的场合,只要它们不背离本发明的原理和精神即可。
另外,为了避免使本说明书的描述限于冗繁,在本说明书中的描述中,可能对可在现有技术资料中获得的部分技术细节进行了省略、简化、变通等处理,这对于本领域的技术人员来说是可以理解的,并且这不会影响本说明书的公开充分性。
图1为根据本发明的实施例的裂缝检测方法的流程示意图。如图1所示,本方法主要分为三个阶段,其中第一阶段为裂缝增强(线特征增强),主要涉及点云数据滤波、最小二乘拟合等算法;第二阶段为裂缝边缘提取,主要应用梯度方向直方图来提取边缘;第三阶段为裂缝目标的提取,主要运用方向分水岭的算法来提取目标。
下面,依次说明上述三个阶段的实现方法。
(1)裂缝增强
为了有利于后续的裂缝识别处理,可对路面深度图像线特征进行增强处理。这里,可采用传统的中值滤波和最小二乘拟合相结合的方法。首先,从局部特征出发,用中值滤波对线特征进行增强,消除系统和外界噪声;然后,从全局特征出发,运用最小二乘拟合的方法进行线特征的增强,消除数据采集过程中的无效值。
根据线结构光的三角测量原理,路面深度图像描述高精密三维物体表面空间信息分布。其测量基本原理如图3所示,由线结构光光源在空间投射出一个光平面,当光平面与被测物体表面相交时在物体表面产生一个亮光条,利用3D相机采集物体表面的光条图像;若被测物体表面的几何形状变化,投射的光条发生形变,根据三角测量原理,从变形的光条图像信息中获取被测物体表面的三维轮廓信息。并将深度数据转换灰度数据,从而形成被测表面的深度图像,即,每个像素的灰度大小代表该点的深度大小。
中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术,中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替,从而消除孤立的噪声点。这里使用二维滑动模板,将模板内深度值按照大小进行排序,生成单调增大(或减小)的为二维数据序列。二维中值滤波输出为g(x,y)=med{f(x-k,y-l),(k,l∈W)},其中,f(x,y)、g(x,y)分别表示原始深度图像中心点邻近的深度值和处理后中心点处的深度值,med{}表示中值滤波操作,W为二维滤波模板,通常为3*3,5*5区域,也可以是不同的形状,如线状、圆形、十字形等。
由于路面局部区域的强反射、强吸收和裂缝两侧遮挡,路面裂缝深度图像存在许多值为零值点。这里采用最小二乘法对扫描线上的数据(深度数据)进行曲线拟合。从整体上考虑近似函数p(x)同所给数据点(xi,yi)(i=0,1,…,m)误差ri=p(xi)-yi(i=0,1,…,m)的大小,常用的方法有以下三种:一是误差ri=p(xi)-yi(i=0,1,…,m)绝对值的最大值即误差向量r=(r0,r1,…rm)T的∞—范数;二是误差绝对值的和即误差向量r的1—范数;三是误差平方和的算术平方根,即误差向量r的2—范数;前两种方法简单、自然,但不便于微分运算,后一种方法相当于考虑2—范数的平方,因此在曲线拟合中常采用误差平方和来度量误差ri(i=0,1,…,m)的整体大小。
数据拟合的具体作法是:对给定路面深度数据(xi,yi)(i=0,1,…,m),在取定的函数类Φ中,求p(x)∈Φ,使误差ri=p(xi)-yi(i=0,1,…,m)的平方和最小,即,使最小,
并输出值p(xi)为拟合值,其逼近于真实路面值。
(2)裂缝边缘提取
通过对路面深度图像处理的实验和分析,运用梯度方向直方图算法(参见参考文献[1])提取裂缝边缘。在不同方向上计算梯度直方图,得到影像的梯度图,然后对梯度影像进行滤波,对边缘信号进行平滑,最后计算方向非极大值抑制,即将不同方向上的非极大值设置为零值,从而实现对裂缝边缘的提取。其技术路线图如图2所示。
梯度方向直方图的边缘提取主要包括以下几个步骤:
(2-1)计算直方图梯度
在影像的多个尺度空间中,计算八个方向下对所有扫描点的直方图梯度。首先,计算方向为θ的每个像素左右两个半圆盘区域内的直方图,获得相应邻域内的深度值统计信息。利用高斯函数对直方图的直方柱计算一维滤波卷积,形成平滑的函数曲线。该高斯函数的参数由宽度因子σ和直方图的直方柱(bin)个数共同确定。接着,需要归一化直方图的直方柱。归一化在每一个圆盘中进行,一般采用的归一化函数有以下四种:
a)L1-norm:v→v/(||v||1+ε);
b)L2-sqrt:
c)L2-norm:
d)L2-Hys:方法同上,限制v的最大值到0.2,并再次归一化处理;
其中,v为直方柱的数值,||v||1是v的1范,是v的2范,ε可为任意值。
归一化直方图的深度值统计信息,让裂缝深度特征对边缘变化具有鲁棒性。本发明可使用L1-norm方法。
最后,由下面的公式(2)计算方向为θ的每个像素左右两个半圆盘区域内的直方图梯度。
其中,gθ和hθ是邻域半径为r、方向为θ的左右两个半圆盘的直方图。K是直方图直方柱的个数。θ表示特定方向上的弧度,且m∈[0,n),表示m某个特定的分片,n表示对π等分的个数。
(2-2)梯度影像滤波
由于路面深度图像成像时自然光照明条件影响、地面地形影响、裂缝多样性等原因引入各种噪声,它会产生多个检测峰,让边缘曲线变得不平滑。由于噪声的影响会带入梯度影像中,需要检测并予以去除。
常用的曲线平滑法中,均值滤波器平滑法和加权平均法都没有考虑曲线本身的趋势,而最小二乘平滑法则假设曲线具有某种数学特性,采用多项式拟合曲线趋势。在曲线平滑方法中,Savitzky-Golay滤波算法是一种经典的最小二乘平滑方法,它使用简化的最小二乘拟合卷积方法对曲线进行平滑处理并可计算平滑后曲线各阶导数。
通常,可假设曲线为p次多项式,即
Yi=a0+a1i+a2i2+…+apip (3)
式中,Yi代表第i点平滑后的数值,用上述多项式拟合曲线的误差为:
式中,yj表示平滑前的数值,平滑窗口大小k=2m+1。为使误差S为最小,对S作偏微分,并使S的各项偏微分等于零。
Savitzky和Golay对整个求解过程进行了推导,给出了平滑窗中心点平滑后数值的最后公式及公式中系数的计算方法。虽然高次多项式能够捕获高峰和窄峰,但是在平滑宽峰上平滑性能不足。本发明采用二阶Savitzky-Golay滤波器对每一个边缘扫描点在8个方向分别进行拟合。在每个扫描点处,使用抛物曲面拟合椭圆区域。椭圆区域的长轴为圆盘的半径(尺度),短轴为长轴的四分之一,拟合方向为π/2、3π/8、π/4、π/8、0、7π/8、3π/4、5π/8。
实验结果表明使用二阶Savitzky-Golay滤波对增强边缘曲线的局部极值,去除噪声,平滑去掉多个检测峰,同时最大程度地保留边缘的细节。
(2-3)方向非极大值抑制并生成边缘影像
使用线性插值,在3×3邻域内对给定方向的垂直方向上执行非极大值抑制。
沿着该方向的垂直方向,局部最大值应远远大于它的邻接元素的插值。为每一个扫描点指定一个方向,默认该方向小于等于π/2作为方向矢量。沿着局部方向进行矢量投影,投影的长度用来确定局部极值。在实际投影过程中,实际角度等于其方向角度减去投影线的角度。对于非极大值设置为0,剩余的像素作为边缘。注意,影像的梯度矩阵应保证为非负矩阵。
(3)裂缝目标提取
分水岭变换(参见参考文献[2])是一种基于数学形态学的影像分割方法,本质上是模拟水淹没地表的过程。传统分水岭变换具有简单、速度快、能检测出弱边缘对象及能获得对象完整边界等诸多优点,因其一般在梯度影像进行变换,受噪声等因素的影响较大,存在大量伪局部极小值区域,会出现过分割现象。由于路面状况复杂,路面深度影像在成像过程容易产生噪声,特别是裂缝空间分布的不确定性和复杂性,利用传统分水岭方法处理深度影像时容易出现大量破碎区域,同时,在强边界附近易于造成交叉,使分割效果变差。
为改善分割效果,本发明应用边缘方向改进分水岭变换算法,其主要思想:由传统分水岭变换对上一节得到的边缘影像进行区域提取,通过边缘方向的拟合并修正边缘的强弱,获取区域分割结果。
利用前一步得到的裂缝边缘方向,下面采用分水岭算法获得裂缝的连通域,得到裂缝区域的轮廓边界。
下面说明本发明采用的离分水岭算法的原理。
首先,假设路面深度影像f是完全下降的,也就是说非局部极小值的点云位置点必然拥有一个小于其深度值的邻域位置点。当影像不服从这一假设时,可以利用完全下降变换对原始影像进行一次变换。该变换fLC的公式为:
公式中d(p)的定义为:
上式中l(π)为路线π的测地距离,表达了所有起点为p,终止于D某点q(f(p)>f(q))的下降路线集合。另外,参数
假设完全下降的影像的局部极小值集为{mi|i∈I},那么任意一个局部极小值点对应的积水盆地区域可以写作:
分水线定义为:
利用分水岭变换算法获得裂缝的连通域,将连通域边界进行标识,表达裂缝边缘点信息,得到裂缝区域的轮廓边界,从而提取裂缝目标。
由于数字影像并非连续的表面,本发明运用离散域分水岭变换算法,提取路面深度影像的裂缝目标区域,涉及到的概念与定义可参考文献[3]。
综上所述,本发明提出的基于方向分水岭变换算法提取裂缝目标方法主要包括以下几个步骤:
1)根据上一节HOG(方向梯度直方图)的处理结果可知,提取路面梯度影像所有位置点八个方向下的全局梯度的最大值,将其作为边缘,得到裂缝边缘影像。
2)将深度影像中深度极小值作为种子位置,采用分水岭算法得到分割区域。
3)提取裂缝边界和邻接区域信息,作为裂缝目标。
下面具体说明步骤3)中裂缝目标提取的实现方法。
3-1)计算离散边界
a、在图像整数矩阵中对图像中的连接成分进行标识。对于整数矩阵,相邻元素具有相同值的才被认为是相连接的。零元素的位置是空白区域,里面没有任何地物类型。每一个像素赋给一个整数值作为标识,用来表明它属于图像的哪一个成分。整数值为零的标识的像素不属于图像的任何成分的元素。标识的整数值的取值范围为[0,n],n是图像中成分的个数。
利用掩膜矩阵来确定中心像元与哪些邻域内的像元连接情况。这里,使用模板为大小3×3的矩阵。注意,保证模板元素非空,且关于中心像元对称。
b、从已标识的连接成分中提取边界,并将每一个连接成分拆分为一组细小边界。当且仅当该小边界在节点相交处才发生拆分。
c、部分拆分边界
按照下述准则,利用递归方法将小边界细分更小的边界,直到这些边界近似为直线。准则1:如果边界两个端点连线与一个端点到边界上其余点连线之间的夹角大于给定角度,就将该边界断开更小的边界。准则2:边界上某点到连接边界两端点线段的距离超出该线段乘以某个给定系数,同时大于给定的距离阈值。
按照准则1或准则2,一条边界总是在边界点与边界两端点连线距离最大的该点处进行分段。
d、整体拆分边界
按照当边界两端点的连线仅仅在边界节点相交时整体分开的准则,利用递归方法将边界分为更小的边界,直到它们的直线段近似反映真实的边界拓扑关系。
e、利用约束Delaunay三角剖分方法确定边缘
计算轮廓边缘近似为直线的约束Delaunay三角剖分(Constrained DelaunayTriangulation,CDT)。将上述拆分的所有边界进行CDT计算来修正边缘。Delaunay三角网具有空外接圆,以及最小角最大的性质,可最大限度的保证网中三角形满足近似等边(角)性,避免了过于狭长和尖锐的三角形的出现,是公认的最优三角网。
以近似直线段的轮廓边缘的两个端点,作为二维实数域上点集V。以相应的轮廓边缘进行约束边的嵌入,对于约束边,其两端点必定在剖分结果中。计算CDT,得到Delaunay三角剖分DT=(V,E),点集V是相应的边界顶点,从而修正边缘。这样,确立约束的边缘构建过程。在整个计算过程中,为保证CDT计算结果的存在,让直线段能够近似反映真实的边界拓扑关系。针对这一点,在步骤(4)的过程中已经增强这个条件。
在不改变边界点集的情况下实现约束边的嵌入,这符合数据采集及数模生成的规则,可以方便地处理图像边缘上的间断线、陡峭线等,真实地反映边缘情况。
3-2)提取边界两侧的区域
提取每一条给定边界的左右两侧的区域。左右区域的宽度由相应边缘的长度的比例因子确定,短边缘确定宽度较小的区域。通过这种方式,由边界的不同长度来确定区域的大小。
3-3)创建顶点/边缘映射图
检测某个扫描点是否为顶点;检测扫描点是否位于边缘上;并求出顶点扫描点的标识,确定扫描点的标识。
3-4)提取顶点坐标这里,顶点坐标以图像的行列号进行表示。
3-5)提取边缘的端点
3-6)提取边缘上各个扫描点的坐标
3-7)提取封闭边缘的标识
(4)最后,进行边界拟合(边缘拟合)。
从每个弧段的局部几何形状出发,估计弧上每个扫描点的方向。用线段去拟合每条弧段。具体方法:弧段端点连线与弧段上任意一点的距离大于某个阈值就表示需要使用线段拟合该弧段,通过迭代,当距离不大于设定阈值时,停止拟合过程。这样,通过逼近方法,将多条线段替代弧段,弧段就表示称为尺度不变的分段折线。这样,得到相应折线段上的像元(x,y)的方向o(x,y)∈[0,π)。
实施边界拟合的具体步骤:
1)计算每条小边缘的方向取近似直线段的小边缘任意两点的坐标,计算两点的夹角,并转化为相应的弧度,范围在[0,π)之间。
2)按照弧度的大小分别将这两点之间的方向设置为图4中0-7的对应区域方向,即将小边缘的方向记为ID,ID∈[1,8]。
3)取对应区域边界方向的边缘强度值,作为小边缘的边缘强度值。
参考文献列表
[1]Dalal N.,Triggs B..2005.Histograms of Oriented Gradients for HumanDetection[C]//IEEE.Computer Vision and Pattern Recognition.San Diego:IEEE,2005:886-893.
[2]Vincent L.,Soille P..Watersheds in Digital Spaces:An EfficientAlgorithm based on Immersion Simulations[J].IEEE Transactions on PatternAnalysis and Machine Intelligence,1991,13(6):583-598.
[3]Kanda F.,Kubo M.,Muramoto K.,2004.Watershed segmentation andclassification of tree species using high resolution forest imagery[C]//IEEEinternational Proceedings Geoscience and Remote sensing Symposium,2004:3822-3825.
最后所应说明的是,以上具体实施方式仅用以说明本发明的技术方案而非限制,尽管参照较佳实施例对本发明进行了详细说明,本领域的普通技术人员应当理解,可以对本发明的技术方案进行修改或者等同替换,而不脱离本发明技术方案的精神和范围,其均应涵盖在本发明的权利要求范围当中。例如,在裂缝增强部分的技术方案中,中值滤波方法可以由方向滤波替代,最小二乘拟合方法也可以由其他多项式拟合方法替代。在裂缝边缘提取的技术方案中,除梯度方向直方图能够检测边缘外,还有边缘方向直方图方法、SIFT描述子和上下文形状方法,同样可以实现所提的裂缝边缘识别方法。