基于L1中值骨架提取的植物茎干补全研究

姜金岑，李联队，王美丽,3,4*

(1.西北农林科技大学信息工程学院，陕西 杨凌 712100；2.陕西省林业科学院，陕西 西安 710082；3.农业农村部农业物联网重点实验室，陕西 杨凌 712100；4.陕西省农业信息与智能服务重点实验室，陕西 杨凌 712100)

植物的数字化与可视化技术是实现数字林业的重要途径，对植物进行可视化建模能够对植物的生长过程和形态结构进行定量化研究和智能化处理。随着三维扫描设备的普及与技术的发展，基于植物点云数据的三维模型重建越来越受到研究者的重视，且具有精度高、真实感强以及细节特征显著等优点[1]。在对植物进行三维扫描过程中，由于外界环境的干扰、视点选择的不同，扫描得到的三维点云通常存在一定的噪声，噪声会使得模型失真，同时也会产生被测物体部分缺失的现象，直接对其进行重建将无法得到一个完整的模型。因此植物点云数据的处理及补全是植物三维建模中的重点和难点问题。

通过三维激光扫描技术对植物进行建模的研究一直都是计算机图形学的研究重点[2-7]。此外，点云分割是三维点云重建中的关键技术之一，也是三维重建的前提[8-13]。但重建过程中会发现扫描得到的三维点云模型常缺失部分细节信息，若需进行点云修补往往需要提取其骨架[14-19]。其中，Huang等[15]提出了基于 L1 中值的点云骨架提取方法和局部密度权重来解决点云密度不均匀的问题，但该方法需要调试大量参数来获取精确的骨架，实用性不高；Zhang等[17]使用概率图模型基于最小距离对骨架节点间边的权重进行计算，通过最小生成树判断每个节点间的拓扑连接情况。但通常这些方法仅适用于人体或其他简单物体点云数据，无法对植物点云数据的骨架进行拓扑连接。

近年来，Charles等[20]提出了可以直接处理原始三维点云数据的神经网络PointNet，能完成高层次的点云分类和语义分割任务。在此基础上许多研究者通过深度学习神经网络对点云进行修补[21-26]，但通过深度学习的方法往往是数据驱动的，需要对大量的数据进行训练，无法广泛应用到植物建模等领域中。针对上述研究现状，本研究拟解决的关键问题是：①构造植物点云的拓扑连接情况，判断缺失部位连接点；②确定待拟合点的点集，并保证断连部分衔接平滑且具有真实性；③对缺失部分进行拟合并调整点云密度，使补全结果与原始点云浑然一体。通过解决这些问题，以期能够在一定程度上弥补扫描时外界因素所带来的缺陷，使植物模型更为生动形象且有效地应用于教学、科研等方面，也更为直观地展现植物的真实形态，为数字林业提供依据，助力林业在信息化时代的发展。

1 材料与方法

1.1 植物点云数据选取

使用三维扫描仪对植物进行扫描，受植物形态结构、遮挡问题和扫描精度的影响，扫描时往往会丢失部分点云数据，其中植物的叶片与其茎干的连接处通常较细容易出现缺失情况，如图1A所示。为了进行定量指标验证，本研究使用一些完整的植物CAD模型，通过虚拟扫描技术获取其三维点云，去除茎干连接处部分点云后进行模拟，如图1B所示。

1.2 植物缺失部分定位

1.2.1 确定缺失部位连接点

不同植物经过扫描后其缺失的点云部分不尽相同，因此需要先定位出缺失的点云部分。对于输入的植物点云，首先对其坐标进行归一化，即将点云坐标映射到(-1,1)区间内。然后使用K近邻算法得到点云中每个点的最近邻点，计算两点间的距离，最终可由所有点与其最近邻点间距离的平均值，得到点云密度。点云密度计算公式：

(1)

其中：Dplant表示点云密度，Cplant为点云集合，‖·‖2表示任意两点间欧氏距离，NCplant表示点云集合中点的数量。得到点云密度后，可根据点云密度判断点云中任意两点间的距离并对点云进行聚类，聚类标准为：

(2)

δ=σ·Dplant。

(3)

其中：dP,Q为P、Q两点间距离，δ表示用于聚类分类的距离标准，σ表示阈值(默认为3.0)。若dP,Q≤δ，则P、Q两点属于同一簇；若dP,Q>δ，则P、Q两点不属于同一簇。

使用概率图模型确定每个点云簇间的拓扑连接情况，同时确定每个点云簇的连接点。将每个点云簇都视为一个图结构中的节点，每两个节点之间边的权重定义为：

(4)

其中:C1、C2是两个点云簇，WC1,C2为两个节点间边的权重，P1是集合C1中的点，P2是集合C2中的点，NC1、NC2分别是集合C1、C2中点的数量。

如果两节点间边的权重越小，那么这两个节点越有可能相连接；反之，如果两节点间边的权重越大，那么这两个节点越不可能相连接。通常而言，每两个节点间边权重的定义应当仅与这两个节点对应的点云簇中点的最近距离有关[17]，即：

(5)

然而植物的叶片、花朵等部位不一定与茎干距离最近，存在叶片与叶片距离最近或叶片与花朵距离最近等情况，故本研究采取最近距离与点的数量相结合的权重定义方式。在概率图模型中使用Prim算法构建最小生成树，使得生成的树形结构边的权重和最小，即得到概率图中最大概率的拓扑连接情况。

构建最小生成树后，可依据最小生成树中节点的度判定茎干点云。Huang等[15]对点云使用基于L1中值的局部迭代方法，可得到邻域内点的局部中心，并逐渐扩大邻域范围，最终获得点云大致的骨架。然而该方法中需要调节的参数较多，整体流程上不够鲁棒，故本研究仅对茎干点云使用L1中值局部迭代方法提取茎干部分的骨架。对于提取出的茎干骨架点集，整体连线近似于茎干的轴线，点集中的每个点则近似于茎干的中心点，为后续确定茎干上的连接点及使用Bezier曲线对缺失部分进行拟合等步骤奠定基础。使用该方法提取得到的茎干点云的骨架点集如图2所示。

由于点云具有无序性，对后续步骤的实施造成影响，因此需要对茎干骨架点集进行排序。本研究提出一种基于最近邻点距离的搜索算法，首先迭代最近邻点判断输入点集的端点，进而从端点出发寻找当前点的最近邻点即可对茎干骨架点集进行排序，且能满足旋转不变性。该算法具有较高的鲁棒性，能够应对不同顺序、不同稀疏程度的骨架点集。通过该算法寻找点集端点并对点集进行排序的流程图如图3所示。

对任意输入点集随机选取其中一点寻找点集端点的过程如图4B所示，图中点颜色的渐变(从红色到白色)为搜索算法寻找过程，最终对点集排序的结果如图4C所示，图中颜色的渐变为点排序顺序的结果。

本研究需要对叶片点云及茎干点云间缺失的部分进行补全，找到叶片点云及茎干点云中的连接点。茎干的拟合及连接应当从茎干中心出发，即茎干部分的连接点应当在茎干骨架点集上。本研究利用叶片点云重心与叶片连接点的连线同茎干轴线的交点作为连接点。对于植物叶片而言，让连接茎干指向叶片重心在现实中能起到对叶片的支撑作用，因此按该方法定义的茎干骨架连接点在物理层面更具有实际意义。对叶片点云重心的求取是计算叶片点云中所有点坐标的平均值。

但实际上三维空间中两条直线不一定存在交点，因此对于叶片重心P与叶片连接点Q的连线PQ与茎干轴线不一定存在交点。故在实际算法执行中要遍历茎干骨架点集，寻找点集中距离连线PQ最近的点R，则点R为茎干骨架点集连接点。本研究使用海伦公式计算三维空间中点到直线的距离(d)：

(6)

(7)

(8)

式中，p为半周长。对茎干骨架点集中的每个点R，计算其到连线PQ的距离d，在茎干骨架点集中到连线距离最小的点即为茎干骨架点集连接点。

1.2.2 确定待拟合点的点集

确定缺失部位连接点后可得到茎干骨架连接点，但仍需确定茎干骨架点集上参与拟合的待拟合点。根据植物的生长特征，茎干是由底部向上生长，因此需在茎干连接点下方的点中选取待拟合点。对于已经过排序的茎干骨架点集，仍无法确定其顺序的方向，在保证旋转不变性的前提下，仍需对茎干骨架的排序顺序进行调整，确保其在茎干轴线上从下到上排序的顺序。通常对于植物的形态而言，植物的上半部分存在叶片、花朵等部位较为丰富而下半部分仅存在茎干较为稀疏的现象，因此可通过植物点云中点的密集程度来判断植物在形态意义上的上下位置关系。

为了近似将植物点云从中线划分成上下两个部分，取茎干骨架点集的中点，及当前点集顺序下一点，以这两点连线的方向向量(指向茎干骨架点集的中点)表示茎干轴线方向向量，过茎干骨架点集的中点作茎干的横截面，以该截面为界分为两个部分的点集，如图5所示。通过两部分点集内点的数量来确定上半部分点集与下半部分点集，若上半部分点集内的点与茎干的中点连线的方向向量(指向茎干骨架点集的中点)与茎干轴线方向向量相同，则表明当前茎干轴线方向与植物形态方向相同；反之若两向量方向相反，则表明当前茎干轴线方向与植物形态方向相反，则需将茎干骨架点集顺序反转，从而将茎干骨架点集的方向与植物形态意义上的方向调整一致。

本研究使用球邻域搜索选取待拟合点，保证拟合时茎干骨架上参与拟合的点与缺失部分需要拟合生成的点所占比例的权重近似相等。以茎干骨架点集连接点为球心，以该点与叶片点云连接点间的距离为半径，该球范围内在茎干骨架点集中的点都将参与拟合。通过这种方法对缺失部分拟合得到的茎干半径与茎干轴线将更为平滑，具有较强的真实感。

接着，需要确定缺失部分待拟合点，即叶片点云与茎干骨架间的待拟合点。为保证拟合效果更为平滑，叶片点云连接点到茎干骨架点集连接点的距离与茎干骨架点集连接点到上述选择的待拟合点点集中最下方点的距离应相当。故依据茎干骨架点集上的待拟合点的数量，在叶片点云连接点与茎干骨架点集连接点的连线段上均匀取点，这些点将作为缺失部分的待拟合点参与茎干半径的拟合。

1.3 植物缺失部位拟合及补全

1.3.1 拟合缺失茎干半径

缺失茎干往往在现实植物中较细，在扫描结果中是需要补全的部分。对于该部分茎干，靠近叶片处茎的半径应当较小、较贴合叶片的厚度，靠近主轴茎干部分茎半径应当较大并且较贴合主轴茎干的粗细。首先定义叶片点云连接处茎干半径，从植物主轴茎干连向叶片部分茎干的粗细应当是逐渐减小的，对于该部分半径的拟合应由主轴茎干上的半径(即茎干骨架点集对应横截面拟合圆的半径)与叶片点云连接处半径共同决定。

对于主轴茎干上的半径，要得到茎干横截面上近似的轮廓线，而轮廓线则是由茎干点云在茎干横截面上投影得到的一系列点。在植物原始点云中，以茎干骨架点集中每个点为中心，使用K近邻(默认K为30，即周围30个点参与拟合)算法计算得到与之相距最近的K个点。以茎干骨架点集中的点P0(x0,y0,z0)及在该点处的方向向量n0(a0,b0,c0)为横截面的法向量，由三维空间平面点法式得到该横截面S的方程为：

a0(x-x0)+b0(y-y0)+c0(z-z0)=0。

(9)

由K近邻算法得到的每个点Pk(xk,yk,zk)，过该点且垂直于横截面直线的方向向量即为横截面法向量n0，故该直线L的方程为：

(10)

由式(9)与式(10)可得K近邻点集中每个点在茎干横截面上的投影点P′k的坐标为：

(11)

由式(11)可以看出，对于投影点P′k(x′k,y′k,z′k)的坐标公式具有相似性，且该公式对于三维空间内任意一点同样适用，因此对于点在茎干横截面上投影的求取算法具有旋转不变性的特征。

根据植物的特征，植物的茎干在局部范围内大多呈圆柱体状，因此以茎干轴线作茎的横截面，截得植物茎干点云形成的轮廓线近似可看作圆，故需要在横截面内以上述K近邻点的投影点拟合圆。在三维空间内对圆的拟合，需满足待拟合点点集中所有点在同一平面内，先拟合球体，再联合球与平面的方程得到三维空间内圆的方程。使用最小二乘法对点拟合三维空间内球体：

(12)

其中：Pb(xb,yb,zb)为拟合球体的球心，Q为点P在K近邻范围内点集中的点，rb为拟合球体的半径。由式(11)可得点Pb在茎干横截面上的投影点P′(x′0,y′0,z′0)，该点将作为迭代更新优化后新的茎干骨架点集中的点，由此可计算拟合球体与茎干横截面所在平面相交圆，即茎干在点P处拟合圆的半径rP：

(13)

式中，h为拟合球心到茎干横截面的距离。

根据最小二乘法拟合得到的圆如图6所示。

对于叶片点云连接处茎干的直径，应当近似为叶片点云的厚度。由于点云的无序性且叶片存在弯折情况，因此无法根据叶片点云的边界来判定叶片的厚度。考虑到叶片点云与叶片点云连接处周围茎干部分点云的密度应当近似，即两部分在连接处应当分布均匀且较为平滑，故求取叶片点云中所有点K近邻(经过实验分析设定K为3)范围内点间距离平均值的平均值，此距离长度则近似为叶片点云连接处茎干的直径。叶片点云连接处茎干的半径(rleaf)为：

(14)

其中：Cleaf为叶片点云，P为叶片点云中的点，CKNN(P,K)为点P在叶片点云中K近邻范围内的点集，Q为该点集中的点，NCleaf为叶片点云中点的数量。

由于茎干缺失部分的半径应当是逐渐减小的，对于缺失部分的起始部位，即靠近主轴茎干部分的茎干半径应当贴合于主轴茎干半径；对于缺失部分的终止部位，即靠近叶片点云连接处部分的茎干半径应当贴合于叶片厚度。故对于待拟合半径集合：

R={rt|rt∈Rstem,t∈[0,t0]}∩{rt|rt=rleaf,t=1};

(15)

(16)

其中:Rstem为茎干骨架点集中待拟合点对应的茎干半径集合，即Rstem={rP|P∈Cstem_fit};t为Bezier曲线的参数，t0为茎干骨架点集中最后一个待拟合点对应的Bezier曲线参数，NCstem_fit为茎干骨架上待拟合点点集的数量，NCleaf_fit为缺失部分待拟合点点集的数量。

进而对待拟合半径集合R与Bezier曲线参数t使用最小二乘法拟合二次曲线R′(t)，根据待生成点集中点的数量计算出待生成的Bezier曲线参数t′，从而计算出待补全部分的茎干半径集合R′。

1.3.2 拟合密度均匀的缺失茎干

参与Bezier曲线拟合的控制点数越多，即Bezier曲线的次数越高，拟合生成的曲线段越贴近控制点，且在空间内弯折较多。若将待拟合点点集中所有点均作为控制点参与拟合，得到的缺失部分茎干轴线将过于弯曲，真实感欠佳，如图7A所示。因此，将拟合二次Bezier曲线作为缺失部分茎干的轴线，该曲线为抛物线，得到的拟合轴线将更为平滑，真实感较强。二次Bezier曲线的3个控制点分别为叶片点云连接点、茎干骨架点集连接点和茎干骨架点集中最下方点。由Bezier曲线锁定端点控制点性质，Bezier曲线拟合生成的曲线段必过叶片点云连接点，故茎干在缺失部分与叶片部分的衔接较完整且平滑。最终使用二次Bezier曲线拟合生成的缺失部分茎干轴线段如图7B所示。

由二次Bezier曲线拟合生成的曲线段作为缺失部分茎干的轴线，以曲线上的点作关于该曲线的切线为法向量，在三维空间中做一系列圆，在圆周上均匀取点，所有点构成的点集即为补全缺失部分茎干所生成的点集。在三维空间直角坐标系中圆的参数方程为：

(17)

其中：圆心为C(x0,y0,z0)，圆的半径为r，向量n1=(a,b1,c1)与向量n2=(a2,b2,c2)为圆面上的两个单位向量，θ为参数，指圆上某点所对应的弧度。

二次Bezier曲线方程为：

B(t)=(1-t)2P0+2t(1-t)P1+t2P2。

(18)

其导数为：

B′(t)=-2(1-t)P0+2(1-2t)P1+2tP2。

(19)

拟合缺失部分茎干的生成点集的密度应当与植物点云密度近似，因此需根据点云密度确定生成点集中点的密度。生成点集的密度由两部分组成：

1)Bezier曲线拟合生成曲线段中点的数量。计算生成曲线段中点的间隔(d=1/N)，其将影响生成点集中点在茎干轴线上的横向距离。由植物点云密度可得植物点云中最近邻点的平均距离，而此距离应当与曲线段中点的间隔相同。由式(19)可得B′(t)：

t∈[0,1]。

(20)

其中曲线段的长度为：

t∈[0,1]。

(21)

因此，曲线段中点的间隔为：

(22)

由此即可确定Bezier曲线拟合时参数t′的集合，进而确定生成点集中茎干轴线上点的间隔。

2)圆参数方程得到圆周上点的数量。圆周上点的数量将影响生成点集中在茎干轴线上的纵向间隔，根据点云密度确定圆周上点间距离：

(23)

基于上述点间距离可确定生成点集密度，使生成的缺失部分茎干点云能够与原始植物点云相匹配。

本研究采用的方法仅需补全缺失部分点集，而使用Bezier曲线拟合得到的曲线段中既包含缺失部分曲线段也包含茎干内部曲线段，因此使用上述方法生成的部分圆周上的点会与原始植物点云中茎干部分的点重合，甚至存在于茎干点云内部，故需删除茎干点云上或其内部的生成点。

通过判断生成点集中每个点到茎干的距离，若该距离小于该处茎干半径则将该点删除；反之，则保留。对于求取某个生成点P到茎干的距离，需先得到茎干骨架点集中到点P距离最近的点Q，及茎干骨架点集中下一点R。根据海伦公式求得点P到茎干线段QR的距离d。

遍历所有生成点集，若d≤Rstem(P)则将当前点从生成点集中删去。

2 结果与分析

2.1 缺失部位连接点确定结果

在输入植物点云并经过预处理后，对点云进行聚类，确定点云在植物中所属的不同部位。对不同点云的分类结果如图9所示，其中σ表示阈值，用来调节点云聚类的精密程度。σ越小说明两点间的距离较小才能属于同一簇；反之，σ越大说明两点间距离较大也能属于同一簇。本研究设定σ的默认值为3.0。

随后使用最小生成树判断点云间的拓扑连接情况，仅基于最近距离和将最近距离与点的数量相结合的结果如图10所示。

将叶片点云与茎干点云进行连接，需要先找到相连的叶片点云及茎干点云中的连接点。依据植物的生长特性，叶片往往是以距离茎干最近的部位与茎干相连接。因此，对于叶片点云连接点而言，该连接点应当是当前集合中与相连集合中所有点的最小距离中最小的点。

在得到叶片部位的连接点后，根据本文所述方法确定的茎干骨架点集连接点及其对应的最终补全后的茎干效果如图11所示。

2.2 待拟合点集确定结果

待拟合点点集由茎干骨架上的待拟合点及缺失部分的待拟合点构成。两部分待拟合点如图12所示，其中：红点表示缺失部分待拟合点与茎干骨架上待拟合点，这些点都将参与缺失部分茎干半径的拟合；黄点表示叶片点云连接点、茎干骨架点集连接点和茎干骨架上待拟合点中最下方点，这些点将作为Bezier曲线拟合的控制点参与缺失部分茎干轴线的拟合。

2.3 拟合结果密度及细节调整

获得待拟合点点集后，通过Bezier曲线对该点集进行拟合，并基于曲线方程用在其切线方向上作圆拟合缺失部分茎干点云。根据本研究所述方法，对拟合的生成点集的密度进行调整，结果如图13所示。

最后，删去原始茎干点云上或其内部的生成点，结果如图14所示。

2.4 植物茎干缺失部分补全结果

2.4.1 补全效果分析

为了能对实验结果进行可视化比较和定量指标验证，使用一些完整的植物CAD模型，去除其茎干连接处部分点云后进行模拟，得到的实验结果如图15所示。其中图15A为输入的茎干与叶片连接处缺失的待补全点云，图15B为通过本文方法拟合的结果，图15C为用于参照的真实值点云。

对于补全后的缺失部分茎干，半径由茎干至叶片衔接平滑，补全的点云与原始点云密度分布均匀，且补全的茎干整体上弯曲自然，具有一定的物理意义与美观效用。缺失部分补全的茎干与原始点云浑然一体，具有较强的真实感。

2.4.2 定量指标分析

采用倒角距离和均方误差两种定量指标对实验结果进行分析。

倒角距离(chamfer distance，CD)用来评估两个点集S1和S2之间距离(其中S1表示点云真实值数据，S2表示本实验的结果)。eCD表示倒角距离误差，具体评价公式如下：

(24)

式中：|S1|、|S2|分别表示S1、S2点集中元素的数量。采用均方误差(mean square error，MSE，式中记为eMSE)来评估实验结果，公式如下：

(25)

式中，CNN(x)表示S1中的点x在S2中的最近邻点点集，在此使用K=10的最近邻搜索，即|CNN(x)|=10。

在两种评价指标中，本实验的结果都取得了较好的结果，如表1所示(数值越低表示误差越小)。

表1 实验结果误差分析表

2.4.3 鲁棒性分析

本研究方法适用于大多数存在茎干与叶片断连情况的植物点云，提出的算法具有旋转不变性，具有较高的鲁棒性。

为了验证该方法的鲁棒性，对输入点云进行随机角度的旋转，实验结果如图16所示，发现对于旋转过后的点云进行补全仍能达到预期效果。

3 讨 论

本研究以茎叶体植物为研究对象，通过三维扫描设备得到植物点云数据。首先基于点云密度使用并查集对点云进行聚类，通过概率图模型计算点云簇间边的权重，使用最小生成树确定点云的拓扑连接情况；然后对于已聚类的点云簇，针对茎干部分点云使用基于L1中值的局部迭代方法提取其骨架，提出一种基于最近邻点距离的搜索算法对茎干骨架点集进行排序，进而根据球邻域搜索确定待拟合点集；最后使用最小二乘法拟合茎干半径，基于Bezier缺陷拟合缺失部分茎干轴线，以茎干切线方向使用三维参数方程生成补全点集，对点云进行补全，为后续的三维建模提供数据基础和技术支撑。基于L1中值骨架提取的对植物茎干缺失部分进行拟合补全的方法，具有旋转不变性及较高的鲁棒性。实验结果在保证茎干骨架点集提取的真实性及茎干半径计算的准确性的前提下，完成了对缺失部位的补全，并基于点云密度对生成点的点集进行调整，使得补全茎干与植物整体连接恰当，浑然一体，可为后续植物的三维重建与植物的可视化奠定基础。

在补全植物点云的过程中，提取植物的骨架是关键步骤之一。由于植物的复杂性和点云扫描的不完整性，直接使用L1中值局部迭代方法提取的植物骨架会出现骨架不完整、偏移过大等问题。笔者先基于点云密度对点云聚类，得到植物叶片、植物茎干等点云簇，然后仅提取植物茎干点云骨架。但对植物完整骨架的提取仍是一个亟待研究的问题，提取植物点云的骨架将会是对植物点云重建、植物点云分类分割等工作的基础，也是植物可视化领域研究的重难点之一。

深度学习相关方法具有较高的鲁棒性及泛化能力，但多数方法是数据驱动的，需要大量的相关数据进行训练。通过本研究提出的方法可简单、直接地获得完整且效果良好的植物点云模型，进而构建植物模型数据集。进一步可考虑在植物模型数据集中，使用深度学习相关方法进行训练，对点云的分类及点云间拓扑连接情况的判断。如可使用PointNet等深度学习网络，直接对点云进行语义分割，将点云进行分类的同时判断该部分点云的种类，同时基于分割的点云种类及点云的空间集合结构可推导出点云间的拓扑链接情况，此类方法将具有更高的鲁棒性。

此外，由于植物本身结构较为复杂，且现实生活中的植物种类繁多，笔者仅针对茎叶体植物完成对其叶片和茎干连接处的补全，在今后将考虑实现对更多类型植物点云的补全工作，增加本算法的泛化能力，实现对多样性植物点云的处理及建模，从而促进林业信息化的发展。