基于四叉树-八叉树联合索引的法向量计算与聚类优化

建立四叉树-八叉树联合索引后，由于单站扫描数据中的扫描线所定义的点和点的邻接关系在单视角情况下是物理准确的，因此可以直接用这些邻接关系来对点云进行法向量计算，不需要借助KNN等第三方库进行高计算成本的点云邻接关系计算，可以提高效率。前面已经提到，在本索引条件下，点云中的每个点都包含如下信息（u，v，x，y，z），u和v表示该点在单站点云矩阵中的平面坐标，（x，y，z）表示该点的三维坐标，具体的计算方法如下：

（1）根据计算准确度需求和计算成本来平衡选择点的法向量估算的邻域范围，在点云矩阵中，可选择8邻域或者24邻域来进行计算。

（2）由于邻域内的点和点的拓扑关系均非常清楚，本书采用构建局部三角网的方式来计算法向量，计算过程中遵循同一方向进行计算，避免出现法向量计算方向不统一的情况。

计算过程如图4-9所示。

图4-9　单站点云快速法向量估算

计算法向量时，以8邻域计算为例，图中P（i，j）为经过四叉树-八叉树联合索引组织后点云矩阵中的第i列j行的点，根据统一方向，均顺时针或者逆时针，选取6个邻近点构建三角形，分别计算P（i，j），P（i，j+1），P（i+1，j）与P（i，j-1），P（i-1，j）等以P（i，j）为共同顶点的6个三角形的法向量并计算其平均值，由此算得P（i，j）的法向量。对于其他点，采用相同的方式计算其他点的法向量，边界上的点采用最邻近点赋值的策略来计算其法向量。

从上面的算法描述可以看出估算法向量时不需要像最邻近搜索算法那样计算点和点的邻接关系，而是直接利用前文中描述的四叉树-八叉树联合索引结构进行计算，由于点和点的物理邻接关系已非常明确，计算结果准确性和效率都随之提高。法向量估算出来以后，需要根据法向量对所有点进行统计聚类。基于前文所述的数据组织索引结构，本书提出了一种二次聚类策略来提高效率和聚类的可靠性，即单站法向量经过计算之后，依然在单站点云矩阵内部进行点云聚类计算，主要通过法向量统计以及点云的邻接关系来进行邻域增长聚类，同时去除局部噪声点，单站内的聚类完成之后，再通过配准变换矩阵，将法向量转换到全局坐标系。对多站内的数据再进行二次聚类，其示意图如图4-10所示。

图4-10　快速法向量聚类流程图

在四叉树-八叉树联合索引下的聚类过程如图4-11所示，具体描述如下：

（1）以A、B两站点云数据为例，对于A站扫描点云与B站扫描点云，从*.PTX数据格式中读取出扫描线后，按照四叉树进行索引编码，其中A1、A2、B1、B2分别为两站中的两个四叉树叶节点，存储了顺序排列且含有相邻拓扑关系的点云数据。

（2）第一次聚类的对象为A1、A2、B1、B2这一层级的所有节点，本次聚类采用基于法向量统计的区域增长方法，由于点云为有序排列，该聚类过程的计算比较快，可以迅速地将节点内的点云按照设定的法向量夹角阈值聚类成多个点云块（子集）。如图4-11中A1被分为（1）、（2）、（3）三块，B1被分为（4）、（5）、（6）、（7）四块，A2被分为（8）（9）两块，B2被分为（10）、（11）、（12）、（13）四块。每一块点云都有其对应的平均法向量和轮廓范围。

（3）再进入多站配准参数，所有的聚类块将通过配准变换矩阵变换到全局坐标系下，图4-11中A1、A2、B1、B2正好属于同一个八叉树叶节点。由于第一次聚类是在单站点云内完成的，配准后需要将多站点云的数据进行融合，因此在八叉树叶节点内需要对之前的多个点云块进行二次聚类。

（4）二次聚类利用每个点云块的平均法向量和轮廓范围来决定多块点云是否合并在一起。示意图中最后的聚类结果为：一次聚类的点云块（1）、（4）、（8）、（12）由于法向量夹角小于一定阈值且轮廓范围相邻而被合并在一起，同理，点云块（2）、（5）、（9）、（13）也由于法向量夹角小于一定阈值且轮廓范围相邻而被合并在一起。此外，（3）、（6）、（10）点云块以及（7）、（11）点云块也被分别合并在一起。至此，完成了二次聚类的过程，得到了全局坐标系下的聚类结果。

在通过8邻域或者24邻域计算出法向量NV之后，可以计算出每个点的曲率C，本算法先选曲率最小的点（C=min）作为种子点P，选取该种子点的所有邻域点与该点进行法向量夹角比较，当夹角小于一定阈值θ时，将该点加入本子集，继续执行该算法，最后得到所有的分割后的点子集。经过该步骤后，可以通过最小二乘等方法对岩体结构面进行初始的拟合计算，从而计算分组平面的相关参数，用以结构面的最终拟合。具体过程如下：

由前文可知，岩体结构面的方程用式（4-1）表示为：

Ax+By+Cz+D=0(www.xing528.com)

其中，A，B，C是结构面的法向量，D是初始点与平面的垂直距离。

（1）对于每一个点组子集，应用最小二乘法计算结构面方程的向量参数。

PL=（A，B，C，D）=（nx，ny，nz，d）；PL∈｛｝

图4-11　点云聚类示意图

（2）将每个子集的法向量NV赋给该点组里的第一个初始点op。

（3）op与拟合平面PL之间的垂直距离PD可以通过公式（4-4）求得。

（4）定义相邻子集Pi的平面模型参数，包含有坐标、向量和垂直距离这几个属性：

Pi=（X，Y，Z，nx，ny，nz，PD）；Pi∈｛P1，P2，P3，…，Pn｝

根据计算的PD值将整个点云数据分割成主要的子集并且包含适当的点数。其次，根据式（4-5）计算和比较每个子集中PD值最小的点的法向量NV与其他每个点的法向量NV之间夹角θth。

通过适当的设置夹角值θth，所有的点可以法向量NV聚类分成相互平行的平面子集：

PSi=（）；PSi∈｛PS1，PS2，PS2，…，PSn｝

分组后，相互平行的理想平面子集如图4-12所示。

经过以上处理，每个点组子集的迭代次数会减小很多，从而减少总的计算时间，提高计算效率。

图4-12　分组后相互平行的子集