3D Hough Transform在激光点云特征提取中的应用

李明磊，李广云，王 力，李海波，范哲瑞

(1.信息工程大学导航与空天目标工程学院，河南郑州450052;2.68029部队，甘肃兰州730000;3.66240部队，北京 100042)

一、引 言

近年来，激光扫描测量技术得到了飞速发展，测量速度从每秒几千个点提升到百万量级，测量精度在一定距离内可达到毫米级甚至亚毫米级;其应用范围几乎涵盖所有传统测绘的应用领域，包括地形测绘［1］、矿山测量［2］、变形监测［3］、文物保护［4］、虚拟现实［5］、数字城市［6］、逆向工程［7］等。激光扫描测量技术大幅减少了外业数据采集的时间和工作量，内业数据处理成为项目的工作重心，也是难点所在。对激光点云数据进行处理、挖掘点云数据信息是当前激光扫描测量技术的研究热点之一。特征提取方法作为二维图像处理的重点研究对象，在三维激光点云数据处理中也存在类似的研究价值，其不仅可以用来进行多站点云的快速自动拼接，还可以用于点云的简化及后续的建模等任务。Paul Hough［8］于 1962 年提出 Hough transform(HT)，用于二维图像中的直线检测，目前已发展出多种形式，包括直接由二维HT导出的standard Hough transform(SHT)，以及改进的 fast Hough transform(FHT)［9］、probabilistic Hough transform(PHT)［10］和randomized Hough transform(RHT)［11］等。二维 HT 可以容易地扩展到三维空间，即三维Hough transform(3D-HT)。本文在分析三维激光点云特点的基础上，以点云中平面特征的提取为例，研究将3D-HT用于点云特征提取所涉及的关键技术。

二、Hough transform原理

1.二维HT直线探测原理

在二维图像空间，每一个像素点在笛卡尔坐标系中都有一个明确的二维坐标(x，y)，如果图像中某一特征能用函数F进行数学表示，则此特征上的所有像素点都应满足由函数F变换得到的等式。以二维直线特征为例，函数F用式(1)表示。式中，θ取任意角度图像中每一个像素点(x，y)都能在二维Hough空间找到对应的Hough变换后的坐标(θ，ρ)。其中θ和ρ的几何意义是:从原点O向过点P的直线F作垂线L，L与x轴之间的夹角即为θ，ρ为原点O到直线F的垂距，如图1所示，图中FON(foot of normal)表示L与F的交点。

之后利用二维累加器(θ方向按角度划分为若干段，ρ方向按距离数值划分为若干段)将二维Hough空间的点坐标θ和ρ分别按照角度和距离进行累加统计，即对于累加器中的任一角度θ1，利用式(1)求出对应的 ρ1，累加器中(θ1，ρ1)处进行一次累加，这个过程也被称作投票。如果二维图像中存在一条直线，其参数为(θ0，ρ0)，则累加器累加结果在(θ0，ρ0)处存在一个峰值。每个峰值代表一条直线，因此下一步任务就是对累加器的累加结果进行局部峰值的探测。峰值探测最常用的方式是通过划定一个阈值来确定各个峰值的分布，提取出各个票数高于阈值的峰值。

图1 二维空间直线的Hough变换

2.三维SHT平面识别原理

对于三维激光点云的特征识别，对二维HT增加一个参数φ就可以比较容易地扩展到三维HT。如果三维空间中某一特征可以用函数F进行数学表示，如平面、球、圆柱等，3D-HT就可以以某种适当的方式对此特征进行识别。以平面特征为例，在三维笛卡尔坐标系中，F一般表示为

为确定某一平面，必须利用这4个参数a、b、c、d，由于4个参数的取值范围是不确定的，直接以这4个参数进行投票时累加器无法确定范围和投票间隔。因此，3D-HT采用式(3)表示平面特征，将平面参数转换为角度信息。

式中，参数θ表示平面的法向n在xoy平面中的投影和x轴正向之间的夹角;φ表示n和xoy平面的夹角;ρ表示原点O距离平面的距离，如图2所示。

以nφ表示在φ方向将三维Hough空间划分的段数，以nθ表示在θ方向划分的段数，以nρ表示在ρ方向划分的段数，则三维Hough空间共被三维划分为nθ×nφ×nρ块，用作累加器的分区，以进行下一步投票。事实上三维空间中过点P的平面有无穷多个，而由于对三维Hough空间的离散化划分，平面个数被限制到 nθ×nφ个(θ和φ 确定后代入式(3)，ρ也可确定，故个数为 nθ×nφ)，这也就意味着对于每一个点P，需要找到这nθ×nφ个满足式(3)的平面(按划分间隔逐个寻找)，即每个点P在三维累加器中需要投票nθ×nφ次。假设点云中共有N个点，则累加器整个过程就需要N×nθ×nφ次。最终在累加器累加结果中，局部峰值点表示可能存在的平面。

图2 三维空间平面的Hough变换

上述过程即为利用三维标准Hough变换(3DSHT)对平面特征进行识别。从其过程的实现来看，此方法耗费的时间和计算资源较大，因此许多改进算法被用于改善标准Hough变换效率低、耗空间等问题。

三、3D-HT平面提取新方法

仅用三维点的笛卡尔坐标进行3D-SHT平面特征识别还可能会存在一些问题。如图3所示，内层平面和外层平面(就平面本身来说是三维无限的)虽然相交，但事实上却没有相交的点，但3D-SHT识别的结果是将外层平面上与内层平面理论相交的点也识别为内层平面上的点。为解决这个问题就需要对3D-SHT额外增加相应的限制条件，这样的限制条件并不容易找到，而且适用性不一定广泛。从3D-SHT原理来看，此方法仅利用各个点的三维坐标信息，并且将各个点单独考虑，没有利用各点之间的相互关系，而且每个点对360°范围内任意倾角φ和θ的平面都要投票(求出任意倾角φ和θ相应的距离常数ρ)，存在变换效率低、耗空间等缺点。事实上对激光点云进行处理时至关重要的一步就是确定各点之间的关系，以此为基础才能对点云进行下一步分析处理，本文尝试以此种思路出发提出3D-HT的新方法。

图3 3D-SHT平面识别结果

首先利用数据结构快速恢复点云之间的邻近关系，进而估计每个点的法矢(法向计算与重定向)。激光点云中每个点的法向是指该点及其临近的若干点确定的三维表面在此点处的法向。激光点云的法向信息是激光点云数据处理中的一个极其重要的信息，对激光点云进行法向计算也是其数据处理中很重要的一步，之后利用恢复出的每个点的法向信息进行HT的投票。每个点只投票一次，降低投票的次数;利用法向信息进行投票，利用到了点与点之间的邻接关系，每个点与其邻近点的法向比较近似，一般不会发生突变，图3所示的错误识别问题也可解决。

具体过程如下:

1)利用数据结构对点云进行分块，从而可以迅速确定点之间的邻接关系，加快对每个点邻近点的搜索速度，本文中数据结构采用的是八叉树［12］。利用八叉树对某模型的点云数据进行分块的效果如图4所示(视图的方向逆着y轴方向，故只能看到二维效果)。

图4 八叉树分割效果

2)利用数据结构寻找点P0的k个邻近点，用这些点进行最小二乘平面拟合(按式(2)进行)［13］，并依次以此方式计算各个点的法向。算得的法向方向具有二向性(可能与实际情况相反)，通过从某点出发进行邻近传递可对其进行调整。而对于激光点云，每个点的法向n=(nx，ny，nz)T必定逆向激光的入射方向r=(x，y，z)T，因此只需保证式(4)中 d为负，即可调整法向到正确指向。

3)将各个点的法向转换为式(5)所示的以角度和距离参数表示的形式，以利于投票过程的执行。

4)确定累加器的划分间隔，利用各点变换得到的(θ，φ，ρ)进行投票。

5)局部峰值探测，从而找出点云中可能存在的各个平面。

6)利用归类到各个平面上的点集重新按照式(2)进行抗差最小二乘平面拟合［13］，以使提取的平面具有更好的平面度。

从执行过程可以看出，该方法比较充分地利用了激光点云预处理得到的成果(邻近关系的确立、法向的计算与重定向等)，在步骤4)的投票过程中假设点云共有N个点，此方法也就仅需要投票N次，远远少于3D-SHT的N×nθ×nφ次，大大降低了投票过程的复杂度。

由于法向计算方式的影响，在墙角等特征突变点附近的法向计算结果一般是不正确的，因此墙角等处的点集也就不能通过HT进行有效识别。本文方法在步骤6)重新进行了抗差最小二乘平面拟合，得到的不仅有平面法向参数(a、b、c)，还有平面到原点的距离参数d，即激光三维点云笛卡尔坐标系下完整的平面方程式，在平面上的点一定满足相应的平面方程。据此可对法向计算结果错误的特征突变点及其附近点等进一步归类，分离出墙角点等特征突变点:

1)由于墙角点满足步骤6)结果中两个以上的平面方程(平面的交点)，依此可以有效地将其分离出。

2)将墙角点附近点的坐标分别代入步骤6)平面方程式，可以将其重新归类到其应在的平面，并且其法向计算结果也可以得到修正。

四、试验验证

首先是累加器划分间隔的具体确定方法。本文利用Rigel公司的VZ-400扫描仪对某实验室的普通墙面(一般建筑内部的常用墙面具有较好的代表性)进行精扫，如图5所示。图5(a)是对墙面拍摄的照片;图5(b)中上图为在扫描结果中选取墙面比较平整的点云块，下图为点云法向计算结果添加光照显示的效果;图5(c)为将对每个点算得的单位化后的法向(包含3个方向的分量)作为三维坐标显示(分布在半径为1的球面上)的效果，由于墙面不平整、激光扫描仪误差、法向计算误差等因素的影响，对墙面点云算得的法向并不完全相同(若完全相同，图5(c)所有点都重合为一个点)。

图5 实验室墙面照片与点云

对墙面的法向计算结果进行统计，结果见表1。由于y轴与墙面的夹角最大，接近垂直，故其方向误差较小。由于 sin(1°)=0.017 5，对于两个单位向量，坐标差别0.017 5相当于夹角差别接近1°。根据法向计算的中误差将累加器中角度的划分间隔设置为1°(法向夹角小于1°的两个点当作处于两个平行平面上的点，根据距离参数确定两平面是否重合)，根据实际平面特征分布，将距离划分间隔设为5 cm(距离小于5 cm的两平行平面当作同一平面，间隔太小无意义且增加空间消耗，太大则不能有效识别平面特征)。

表1 墙面法向中误差

为测试本文方法的有效性，首先利用两层立方体模拟数据(共12个平面，303 612个点)进行试验。图6即为平面识别的效果，不同的平面之间用不同的颜色来区分。图6(a)所示，由于一般的邻近点平面拟合不能有效计算得到正确的法向，故在平面相交处的点(边角点)并未被有效识别。而通过本文特征突变点及其附近点重新归类的思路(修正方案)，可以有效地对其重新归类，将平面和墙角分离，效果如图6(b)所示。表2为对模拟数据进行平面提取的输出结果(转换成式(2)的4个参数，并记录平面上识别到的点数)，与图6数据模拟时的设定值基本符合(内层平面完全符合，外层平面由于内层平面与其交线的影响有个别点仍未有效归类)。

图6 两层立方体数据平面提取效果

表2 两层立方体数据平面提取结果

为测试本文方法的实用性，利用某实验室内部激光扫描实测数据进行试验。图7即为平面识别的效果，图7(a)为仅依据法向投票识别的效果，图7(b)为利用修正方案识别的效果，图7(c)为将分离出的墙角显示的效果。最终平面提取结果为6个，与实际符合，而且墙面墙角特征可以有效分离。

图7 激光实测数据平面提取效果

五、结束语

本文根据激光点云数据处理的特点，结合激光点云数据预处理中点与点之间邻近关系的确立及法向的计算与重定向，提出了直接根据每个点法向进行Hough变换的方法，极大地降低了3D HT的复杂度，并可以有效解决3D-HT存在的错误识别问题。对特征突变点及其附近点的重新归类方法可以进一步识别大部分没有正确识别的特征突变点，并可以对其法向计算结果进行修正。本文提出的方法仅受法向计算结果制约，具有良好的适用性。试验中首先利用普通常见墙面的激光扫描仪精扫数据确定了三维累加器的划分间隔，之后分别利用模拟数据和实测数据证明了该方法的有效性和实用性。

本文的工作主要集中于激光三维点云平面特征及相应墙角轮廓的提取，下一步可以以此为基础研究更多种类特征的提取方法。

［1］许映林.3维激光扫描技术在温泉水电站大比例尺地形图测量中的应用［J］.测绘通报，2007(6):9-13.

［2］杨青，张亮，王振.基于KD树和CDT的露天矿三维建模［J］.计算机技术与发展，2011，21(12):224-226.

［3］谭国铨.3D激光扫描仪在电厂冷却塔变形监测中的应用探讨［J］.勘测设计，2009(4):30-32.

［4］周俊召，郑书民，胡松，等.地面三维激光扫描在石窟石刻文物保护测绘中的应用［J］.测绘通报，2008(11):68-69.

［5］李晖，吴禄慎.三维激光扫描技术在虚拟现实中的应用［J］.南昌大学学报:工科版，2007，29(3):239-242.

［6］赵海莹，张正鹏.三维激光点云数据在城市建模中的应用［J］.城市勘测，2009(1):69-72.

［7］马有良，李占峰.三维激光扫描仪在曲面逆向工程中的应用研究［J］.机械设计与制造，2009(10):85-87.

［8］HOUGH PV C.Method and Means for Recognizing Complex Patterns:US，Patent 3069654［P］.1962-12-18.

［9］OGUNDANA OO，COGGRAVE CR，BURGUETE R L，et al.Automated Detection of Planes in 3-D Point Clouds Using Fast Hough Transforms［J］.Optical Engineering，2011，50(5):053609.

［10］KIRYATI N，ELDAR Y，BRUCKSTEIN A M.A Probabilistic Hough Transform ［J］.Pattern Recognition，1991，24(4):303-316.

［11］XU L，OJA E，KULTANEN P.A New Curve Detection Method:Randomized Hough Transform(RHT)［J］.Pattern Recognition Letters，1990(11):331-338.

［12］邢坤.海量激光扫描测量数据的处理［D］.郑州:信息工程大学，2009.

［13］李明磊.点云平面拟合新方法［J］.测绘通报，2012(S1):84-87.