具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

在本发明的描述中，除非另有规定和限定，需要说明的是，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是机械连接或电连接，也可以是两个元件内部的连通，可以是直接相连，也可以通过中间媒介间接相连，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

本发明公开了一种基于分叉数量自动估计的点云骨架线提取方法，在一种优选实施方式中，如图1所示，包括：步骤S1，获取点云数据，基于KD－Tree结构构造点云数据的K近邻图；步骤S2，获取K近邻图的测地距离，根据K近邻图的测地距离将所有K近邻图中的点云划分为k个水平集，k为正整数；步骤S3，利用分叉数量自动估计算法自动将每个水平集内的点云划分为不同分支，获取每个分支的质心，将质心作为骨架点；步骤S4，连接步骤S3中获得的所有骨架点获得骨架线。

本发明公开的基于分叉数量自动估计的点云骨架线提取方法优选但不限于应用于人体骨骼重建、人体血管重建以及树木重建等场景中的骨架线提取。

在本实施方式中，由于扫描树木使用的是激光扫描仪，扫描的点数量一般很大，一般情况下，一棵树的点云数量会达到十万级，如果直接在树木原始的点云上进行试验，会由于树木采样点过多，计算的时间较长，同时原始点云存在大量噪声，噪声点也会影响点云骨架的提取，直接在不做任何处理的树木点云上进行树木的骨架提取，提取出来的骨架线效果并不好，提取出来的骨架线和树木真实的拓扑形状差距很大，此外，由于原始点云密度不均匀，提取出来的树木骨架一般会出现偏向点云密度较大的一侧的情况，为了减少原始点云数量过大且存在的噪声和密度不均匀对树木点云提取骨架线的影响，本方法首先对原始点云进行降采样。

在本实施方式中，优选的，如图1所示，在步骤S1中，获取点云数据的过程为对原始点云进行降采样的过程，具体过程包括：

步骤S11，获取原始点云，从原始点云中找到最大点位置坐标(x_max,y_max,z_max)和最小点位置坐标(x_min,y_min,z_min)；

步骤S12，设置第一包围盒l′l′l，l＝max(L_x,L_y,L_z)，L_x＝|x_max-x_min|，L_y＝|y_max-y_min|，L_z＝|z_max-z_min|，将第一包围盒细分出多个第二包围盒l'′l'′l'，l'＝min(L_x,L_y,L_z)，如n'个，所有第二包围盒的体积总和等于第一包围盒体积，满足其中，i'∈[1,n']，V_i'代表第i'个第二包围盒的体积，每一个第二包围盒中都有一定数量的点云，对于包含n”个原始点云点的第i'个第二包围盒。

步骤S13，求取所有第二包围盒的中心点，将中心点作为步骤S1中获得的点云数据中点云点。使用公式计算出第i'个第二包围盒的中心点Q_i'，从而实现了降采样，第i'个第二包围盒中包含有n”个原始点云点，i”∈[1,n”]。

在本实施方式中，由于树木点云的点是没有规则的，同时点与点之间是没有联系的，为了解决树木点云之间的无关联性，首先根据KD－Tree建立出点云数据的K近邻图，使用K近邻图(K优选但不限于为4)来解决树木点云之间的无关联性，本申请使用了KD－Tree结构来加快构建树木点云的K近邻图。相对于普通的K近邻算法，由于需要多次遍历树木点云，时间复杂度较高，而本申请使用基于KD－Tree来构造树木点云K近邻图，时间复杂度和原始构造方法的遍历相比，时间复杂度低得多，本文构建领域图方法如下，给定一个需要构建K近邻图的树木点云的原始数据集，本文定义的KD－Tree结构如下表所示。

表一 KD－Tree结构

构建过程的伪代码为：

KD－Tree建立完成以后，本文可以通过维度信息来选择每个点云最近邻的4个点，建立出来的K近邻图如图2所示。

在本实施方式中，建立好K近邻图以后，需要把K近邻图上距离相近的点云找出来，方便后面的骨架提取，本文采用的是水平集的方法来寻找K近邻图上距离相近的点，水平集算法是一种隐式的表现曲线的方法，就是把低维目标用比他高一维的水平集函数的零水平来表示。水平集可以表示成一个具有n个变量的实值函数f的水平集，一般可以用公式其中，c是水平集常数，当n＝2时称为水平曲线，当n＝3时，称为水平曲面，当n＞3时，水平集被叫做水平超曲面。并且集合称为f的子水平集。对于本文，本文根据测地距离的不同划分出不同的水平集，按照测地距离，这里的测地距离表示为从树木的根节点(人体骨骼骨架线或血管重建等情况时，根节点可作为高度最低点，即Z轴坐标最小的点)到每一个点云点的k近邻图上的距离的不同，将点云数据划分为k个不同的水平集，如图3所示，用公式可以表示为L_(f,c,p)＝{p|c≤f_(p)≤c+ρ}，f_(p)表示水平集L_(f,c,p)中点云点p的测地距离，f_(p)∈[c,c+ρ]；c表示水平集L_(f,c,p)的距离常数，可预先设置，每个水平集对应一个距离常数；ρ表示距离增量阈值，可预先设置；设置邻域半径ε，L_(f,c,p)代表的是一个水平集。如图4所示为最终获得的水平集。

在一种优选实施方式中，每个水平集执行分叉数量自动估计算法的过程包括：

步骤S31，设水平集表示为：L_(f,c,p)＝{p|c≤f_(p)≤c+ρ}，f_(p)表示水平集L_(f,c,p)中点云点p的测地距离，f_(p)∈[c,c+ρ]；设置邻域半径ε；

步骤S32，找出水平集中所有的核心对象点云点，核心对象点云点的寻找过程为：

对于水平集中的点云点p，将与点云点p的距离Dist_(p,q)不大于邻域半径ε的所有点云点构成点云点p的邻域集合n_ε(p)＝{q∈L_(f,c,p)|Dist_(p,q)≤ε}，若点云点p的邻域集合n_ε(p)中点云点数量|n_ε(p)|大于等于邻域数量阈值MinPTs，即|n_ε(p)|≥MinPTs，则点云点p是核心对象点云点，否则，点云点p不是核心对象点云点；

步骤S33，若点云点q满足q∈n_ε(p)，则认为点云点q与点云点p密度可达；如果点云点o∈L_(f,c,p)，并且点云点o与点云点p密度可达，点云点o与点云点q密度可达，则称点云点p和点云点q密度相联；

将水平集中密度可达和密度相联的所有点云点聚集为一个点云点集合，将点云点集合作为一个分支，一个水平集可获得至少一个点云点集合，每个点云点集合作为一个分支。

在本实施方式中，可以以任一核心对象点云点出发，找出所有与该核心对象点云点密度可达和密度相联的点云点作为一个分支，判断所有点云点是否均加入了分支，若没有，以没加入分支的任一核心对象点云点出发重复之前步骤聚类为分支，直到所有点云点加入对应分支，实现了自动分支。

在本实施方式中，分叉数量自动估计算法是根据密度划分分支的算法，而且该算法是对噪声具有良好的鲁棒性算法，能够把密度达到一定程度(包括密度可达和密度相联)的区域划分为对应的一个分支，并且具有发现任何分支的能力。分叉数量自动估计算法的一个最大的优点就是在进行树枝点云聚类的时候，并不需要预先知道需要聚类的树枝的数量，选择合适的聚类参数(包括邻域半径ε和邻域数量阈值MinPTs)，算法就能自动聚类出所有分支。可以定义密度阈值公式为：式中Den是指密度阈值，调整Den可以聚类出不同的分支。算法描述如下：

在一种优选实施方式中，在步骤S4中，连接步骤S3中获得的所有骨架点获得骨架线的过程包括：

步骤S41，对于获取的骨架点，将根节点作为初始搜索点，从根节点开始使用KD－Tree进行半径为r的近邻搜索，将搜索到的骨架点标记为主支骨架点并加入骨架点链表，将搜索到的骨架点作为新的搜索点进行半径为r的近邻搜索直到遇到分叉点，如图5所示，当在搜索点半径为r的邻域内Z轴坐标大于搜索点的Z轴坐标的骨架点的数量大于1时，将每个Z轴坐标大于搜索点的Z轴坐标的骨架点作为一个分叉点。步骤S42，设搜索点为p'，将搜索点p'的所有分叉点存入待搜索集合S_(p')，将搜索点p'分别与待搜索集合S_(p')中的每个分叉点相连获得多个分叉向量，将所有分叉向量之和作为骨架生长方向向量如图6所示，将上一历史搜索点与当前搜索点连接后得到向量作为第一向量，求取骨架生长方向向量与第一向量的锐角夹角θ。

步骤S43，设分叉点为q'，分叉点q'对应的分叉向量为若将分叉点q'作为主支骨架点，若将分叉点q'作为分支骨架点，将主支骨架点和分支骨架点加入骨架点链表；表示求取的方向向量，表示求取的方向向量。

步骤S44，以主支骨架点作为新的搜索点继续进行主支的近邻搜索，以分支骨架点作为所在分支的根节点开始分支的近邻搜索，直到完成所有骨架点搜索，依次连接骨架点链表中的主支骨架点和各分支的分支骨架点得到骨架线。骨架点链表中的主支骨架点可按照搜索到的顺序存储，作为交叉点的分支骨架点可与在与其相邻的Z轴坐标较小的主支骨架点关联存储，每个分支的分支骨架点可按照搜索到的顺序存储。

在一种优选实施方式中，还包括骨架线优化步骤；骨架线优化步骤包括：

步骤A，将步骤S4中获取的骨架线作为初始骨架线，沿初始骨架线生长方向将点云数据划分为多个局部点云集合，所述局部点云集合中仅包含主支的点云或一个分支的点云；

步骤B，将每个局部点云集合内的骨架线近似为三维空间直线，所述三维空间直线表示为：

其中，A1表示第一x轴参数，B1表示第一y轴参数，C1表示第一z轴参数，D1第一平面常数，A2表示第二x轴参数，B2表示第二y轴参数，C2表示第二z轴参数，D2表示第二平面常数；

根据局部点云集合内的初始骨架线分别获得A1、B1、C1、D1、A2、B2、C2、D2的初始值；

分别设置A1、B1、C1、A2、B2、C2的取值区间，每个参数的取值区间包含所述参数的初始值，分别在A1、B1、C1、A2、B2、C2的取值区间进行取值，并根据A1、B1、C1、A2、B2、C2的取值获得D1和D2的取值，当A1、B1、C1、D1、A2、B2、C2、D2的取值对应的三维空间直线满足目标函数时，将此时A1、B1、C1、D1、A2、B2、C2、D2的取值对应的三维空间直线作为该局部点云集合的优化骨架线；

目标函数为：其中，arg min(·)表示使表达式取值最小的函数；d_j表示局部点云集合中的第j个点云点到局部点云集合内三维空间直线的距离，设局部点云集合中包含m个点云点，j∈[1,m]；d_mean表示局部点云集合中点云点到局部点云集合内三维空间直线的距离平均值，

步骤C，按照骨架线生长方向连接所有局部点云集合的优化骨架线，连接所有局部点云集合的优化骨架线得到整个的优化骨架线，如图7所示。

在本实施方式中，当A1、B1、C1、A2、B2、C2取值确定后，根据三维空间直线方程和局部点云集合内的点云点坐标可以求解得到D1和D2。

在本实施方式中，优选的，A1、B1、C1、A2、B2、C2的取值区间分别以各自的初始值为中心点，在初始值上下波动，具体波动范围可根据实际情况设置，如初始值的正负20％。优选的，可设置步进值，各参数依次或同时按照步进值取值，每取值一次，将获得一个三维空间曲线，同时判断是否满足目标函数，若满足停止步进，若不满足继续步进取值，直到满足目标函数为止。

在本实施方式中，优选的，在步骤B中通过如下步骤获得每个局部点云集合的优化骨架线：

步骤B1，设置n个粒子，所述n为正整数，以6维变量(A1,B1,C1,A2,B2,C2)作为粒子的位置变量，粒子的速度变量为6维的，初始化每个粒子的位置和速度，可将每个粒子的速度初始化为0，将每个粒子的位置分别初始化为在A1、B1、C1、A2、B2、C2的取值区间内随机选取值；将历史最优位置和当前最优位置初始化为(A1₀,B1₀,C1₀,A2₀,B2₀,C2₀)，A1₀、B1₀、C1₀、D1₀、A2₀、B2₀、C2₀、D2₀分别表示根据局部点云集合内的初始骨架线分别获得A1、B1、C1、D1、A2、B2、C2、D2的初始值；

步骤B2，在A1、B1、C1、A2、B2、C2的取值区间内按照如下公式更新每个粒子的位置和速度：

v_i(t+1)＝v_i(t)+c₁′rand()′(pbest-x_i(t))+c₂′rand()′(gbest-x_i(t))；

x_i(t+1)＝x_i(t)+v_i(t+1)；

其中，v_i(t)、v_i(t+1)分别表示第t次迭代、第t+1次迭代中第i个粒子的速度，x_i(t)、x_i(t+1)分别表示第t次迭代、第t+1次迭代中第i个粒子的位置，i∈[1,n]，pbest为当前最优位置，pbest为历史最优位置，rand()为随机设置值，t为迭代次数；c₁表示第一学习因子，c₂表示第二学习因子；粒子的位置必须在A1、B1、C1、A2、B2、C2的取值区间内进行更新；

步骤B3，计算每个粒子的适应度值，适应度值计算公式为其中，将当前最优位置pbest更新为适应度值最小的粒子位置，若当前最小的适应度值小于历史最优位置对应的适应度值，将历史最优位置gbest更新为适应度值最小的粒子位置；

步骤B4，若当前迭代次数t+1没有达到预设的最大迭代次数T，令t＝t+1，返回执行步骤B2、步骤B3和步骤B4，若当前迭代次数t+1达到预设的最大迭代次数T，t∈[1,T-1]，停止迭代，将根据历史最优位置获得三维空间直线作为局部点云集合的优化骨架线。

在本实施方式的一种应用场景中，由于在自然界中，树木生长或骨骼或血管是有一定的生成规律的，大多可近似为圆柱形，又由于我们知道圆柱柱面上的点到轴线的距离是相等的，由此本文采用了PSO算法来进行优化初始连接的骨架线，只有当骨架线为圆柱的中轴线时，圆柱面上的各个点到骨架线的方差才能达到最小，由于开始我们已经提取出来了一个初始的骨架线，我们可以计算出初始骨架线的方程，那么我们可以对初始骨架线进行优化来得到居中性良好的骨架线。为此可以采用空间点云粒子群算法不断来优化骨架线参数，使得距离的方差尽可能最小，即：通过降低loss，可以不断逼近真正的骨架线，粒子群算法过程描述如下：首先我们可以把树木点云在三维空间的位置和速度分别表示为矢量，同时，每一个树木点云都有一个目标函数决定的适应值，而且每一个树木点云点知道自己目前发现的最好的位置和现在自己的位置x_i，同时，每一个树木点云也知道目前一个点云簇点云发现的最好的位置，即同伴的飞行经验，每一个点云点x_i就是通过自己和同伴的经验来决定下一步的运动情况。点云粒子群算法首先初始化一群随机的粒子，然后通过不断的迭代优化来得到最优的解。在每一次的迭代过程中，每一个点云点根据自己目前发现的最好位置和所有粒子发现的最好位置来更新自己的位置和速度。用公式可以描述如下：

v_i(t+1)＝v_i(t)+c₁′rand()′(pbest-x_i(t))+c₂′rand()′(gbest-x_i(t))；

x_i(t+1)＝x_i(t)+v_i(t+1)；

其中，v_i(t)指的是记忆项，表示上次的速度大小和方向；c₁′rand()′(pbest-x_i(t))是自身认知项，是目前点云点自己凭借经验找到的最好的位置和速度；c₂′rand()′(gbest-x_i(t))是群体认知项，指的是树木点云簇中最好的位置和速度。

在一种优选实施方式中，还包括：步骤S5，采用三次B样条插值算法顺滑骨架线。B样条插值定义为在已知n₁+1个控制点w_o(o＝0,1,2,3,...,n₁)，则可以定义k'次B样条曲线的表达式为：其中N_o,w(μ)为每一个控制点对于的一个基函数，μ是自变量，o为第o个节点，w_o为节点的次数。

在本发明的一种应用场景中，为了验证本文算法的正确性和有效性，本文选用了几种不同的树木点云来进行算法测试，这几种点云分别是点云数据有噪声点的点云、点云数据缺失的点云、密度分布不均的点云。其中，实验平台使用ubuntu16.04，CPU使用i7－6700HQ，内存为12G，实现本文中的算法使用的是PCL1.8与C++，实验使用的点云数据来自激光扫描仪获取的点云，对树木点云骨架提取出来的结果采用的是视觉评估的方法来评价提取出来骨架结构的准确性。

在本应用场景中，为了试验骨架提取算法的正确性和有效性，试验在有数据缺失的树木点云、有噪声的树木点云、数据分布不均匀的树木点云上进行。将本申请方法和提到的算法和其它骨架提取算法进行了对比：图8为数据有缺失的树木点云通过本申请方法提取的骨架线与现有其它骨架提取算法提取骨架线的对比图，其中，第一列为原始的数据有缺失的树木点云，第二列为现有其它骨架提取算法提取骨架线，第三列为本申请方法提取的骨架线；图9为密度不均的树木点云通过本申请方法提取的骨架线与现有其它骨架提取算法提取骨架线的对比图，其中，第一列为原始的密度不均的树木点云，第二列为现有其它骨架提取算法提取骨架线，第三列为本申请方法提取的骨架线；图10为有噪声的树木点云通过本申请方法提取的骨架线与现有其它骨架提取算法提取骨架线的对比图，其中，第一列为原始的有噪声的树木点云，第二列为现有其它骨架提取算法提取骨架线，第三列为本申请方法提取的骨架线。从图8到图10中可以看出，在大多数情况下，本申请方法提取出来的树骨架具有良好的居中性，效果更好。

在本应用场景中，本申请方法先使用水平集和分叉数量自动估计的点云骨架提取算法来得到初始的骨架，然后使用PSO算法来优化得到的初始骨架，然后得到了居中性良好的树木点云骨架线。本申请方法有以下的优点：

1、在树木点云存在噪声时、数据分布不均匀或者数据存在少量缺失时，可以很好的提取得到居中性良好的骨架线。

2、使用了基于KD－Tree的点云结构来加速寻找点云点，明显提高了分叉数量自动估计的点云骨架提取算法的速度，加快了树木点云骨架提取的速度。

3、使用了PSO算法来优化初始的骨架线，保证了树木骨架线的居中性。本发明还公开了一种基于分叉数量自动估计的点云骨架线提取系统，包括点云数据获取模块和处理模块，点云数据获取模块获取点云数据并传输至处理模块，处理模块执行上述基于分叉数量自动估计的点云骨架线提取方法的步骤获得骨架线。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。