胡建平，刘 凯,，郭新宇，吴 升，温维亮※

（1. 东北电力大学理学院，吉林 132012；2. 北京市农林科学院信息技术研究中心，北京 100097；3. 国家农业信息化工程技术研究中心，北京 100097；4. 数字植物北京市重点实验室，北京 100097）

0 引 言

植物表型研究通过集成自动化平台装备和信息化技术手段，获取多尺度、多生境、多源异构植物表型海量数据，解决了传统植物表型获取工作量大、效率低、误差大、适用性弱等问题，将极大的促进功能基因组学、作物分子育种与高效栽培的进程。随着激光雷达传感器和基于多视角三维重建技术的快速发展，利用三维点云进行3D表型研究成为高通量获取植物形态结构表型的一种有效手段，受到研究人员的广泛关注。

骨架是三维数据和模型的抽象表现形式，能够直观的表达物体的拓扑连接关系和几何结构，尤其适用于描述分枝形植物植株和长窄形叶片的形状结构。因此，研究人员围绕植物三维点云骨架提取开展了大量研究工作。

经典的三维点云骨架提取方法主要包括基于拉普拉斯矩阵和L1中值骨架提取算法。由于植物具有特定的形状结构特征，直接采用这些提取方法提取骨架往往难以得到理想的结果，需结合植物形状结构特征在原方法基础上进行改进和完善。植物点云骨架提取主要包括单株尺度和器官尺度。在单株尺度，宗泽等结合最小二乘法和遗传算法，提出了一种用于计算作物表型参数的骨架提取算法，能够通过提取的骨架得到较为准确的表型信息；对于整株植物，张卫正等通过骨架细化算法建立植株的骨架模型，利用二叉树分析骨架模型，最终确立各级节点及各个叶片实现茎叶分割；Xiang等通过分层聚类，实现高粱植株点云的骨架提取；Wu等通过改进拉普拉斯算法，对玉米植株点云提取骨架，可以得到精确的骨架结果。植株尺度点云骨架提取需要在确保植株拓扑结构准确的前提下，尽可能的提高骨架提取的精度，进一步用于器官分割和表型解析。在器官尺度，赵元棣等通过计算离散网格的平均曲率和网格曲面参数化方法，实现玉米叶片三维点云的主脉曲线提取；Xiao等通过检测点云首尾端点实现椭圆形叶片的叶脉骨架提取。器官尺度骨架提取对精度要求比植株尺度更高，因此对点云缺失更敏感。

叶片是植物的重要组成部分，其形状结构的差异直接影响植物的生长发育和最终产量。由于植物本身器官间存在大量遮挡、植物叶片表面材质属性复杂，导致所采集的植物三维点云数据通常存在缺失。此外，受品种和环境等因素的影响，植物叶片通常具有特定的形状结构特征。因此，直接应用计算机图形学领域一般化的骨架提取算法难以得到理想的植物叶片点云骨架提取结果。本研究针对玉米叶片缺失点云数据的骨架提取问题，考虑到长窄形叶片叶脉曲线两侧叶面结构的对称性，设计了一种新颖的自适应加权算子来计算玉米叶片点云的骨架约束点集，并引入主曲线对骨架约束点集进行优化，可实现缺失点云数据的玉米叶片骨架的鲁棒提取。

1 材料与方法

1.1 算法流程

本文提出的基于自适应加权算子的玉米叶片点云骨架提取流程可概括如下：首先通过K-means聚类将玉米叶片点云数据划分为若干个变化较为平缓的部分；然后对每个部分按设定步长计算多个与叶脉垂直的截面，根据距离提取每个截面附近的点作为当前截面的相关点集；随后通过设计的自适应加权算子提取每个截面相关点集对应的骨架约束点，将从叶片中提取的所有约束点形成骨架约束点集；最后通过计算骨架约束点集的主曲线得到最终玉米叶片点云骨架。算法流程如图1所示。

图1 基于自适应加权算子和主曲线的玉米叶片点云骨架提取流程图 Fig.1 Flow chart of skeleton extraction of maize leaf point cloud based on adaptive weighting operator and principal curve

1.2 数据获取与处理

采用作物单株高通量表型平台MVS-Pheno，于2020年获取吐丝期玉米植株点云数据，并利用玉米植株点云标注软件Label3DMaize分割得到玉米叶片三维点云。由于MVS-Pheno平台是利用多视角成像原理得到的植株三维点云，受获取光环境和叶片遮挡等因素影响，所得到的植株叶片点云存在少量缺失。为验证算法的有效性，从多个玉米植株上选取30个具有典型形状结构特征的叶片，各叶片点云形状结构特征描述及相关信息见表1。

表1 玉米叶片形状基本信息介绍 Table 1 Introduction to the basic information of maize leaf shape

对获取得到的玉米叶片点云数据，首先通过去噪和降采样的方式去除噪声点，并在不降低点云数据质量的情况下减少点的数量，以期提高后续计算效率；然后，采用K近邻（K-Nearest Neighbor）算法确定每个点的邻域，并使用主成分分析（Principal Component Analysis，PCA）算法估计点云中各点的法向；最后，通过最小生成树算法来调整点云中每个点的法向，满足法向在方向上的一致性，最终得到点云数据中所有点的法向信息。对于预处理后的玉米叶片点云数据，再通过K均值聚类的方式将单个叶片分为若干部分（一般叶片分为3～4个部分即可），使每部分的形状变化较为平缓，如图2所示。

图2 叶片聚类结果示例图 Fig.2 Example diagram of leave clustering results

1.3 截面相关点集计算

三维物体的骨架可通过计算垂直于骨架方向的截面与三维物体相交所得曲线的中心点计算（如图3所示）。然而，在欧式空间中垂直于叶片伸长方向（骨架方向）的截面与三维物体的相交曲线不易显式计算，因此，本文采用垂直于叶片伸长方向相关点集替代截面与三维物体相交所得曲线，并用其进行后续骨架约束点的计算。

图3 提取截面相关点及骨架的示意图 Fig.3 Schematic diagram of extracting the related points of cross sections and the corresponding skeleton

本文采用改进的距离场方法计算各截面对应的相关点集，该距离场通过计算点云中的点与其包围盒的距离构建出该点云的距离场，这种方法与传统方法相比可以极大的降低计算成本。

对于通过K均值聚类得到的原始玉米叶片变化比较平缓的某个子叶片点集= {p|= 1,2,...,}，首先通过PCA算法计算该点云数据的三个特征向量1、、，并将其进行标准正交化构成正交矩阵=[′′]，然后以标准正交基1′、′、′构建新的笛卡尔坐标系，构建方向包围盒（Oriented Bounding Box）。计算该部分点云数据中每个点到与叶片伸长方向垂直的正交面的距离并进行归一化，得到该部分点云数据的距离场，然后通过距离场从叶片点云中抽提出每个截面对应的相关点集。距离场计算公式如下：

其中D表示点p在距离场中的归一化距离，表示与叶片伸长方向垂直的正交面（由于整个叶片点云被分为若干个较为平缓的部分，所以该正交面为当前点所在部分包围盒的一端，如图4所示，(p,)表示p到正交面的距离。

图4 子叶片的距离场 Fig.4 Distance field of each sub-leaf

该距离场可将整个叶片点云划分为多个子集，每个子集为一个相关点集，可以近似为一个垂直于叶片伸长方向的截面点集合，具体形式如下：

其中为设定的最大截面数，试验中设置为100，t表示第个截面对应的起始距离。

1.4 基于自适应加权算子的骨架约束点提取

为了克服玉米叶片点云数据采样不均匀以及数据缺失等影响，本文设计了一个自适应加权算子来计算每个截面相关点集S的骨架约束点c，其实质是根据叶片点云数据的空间分布特征自适应的计算每个截面相关点集S的中心，并根据点云数据的缺失情况对计算出来的中心点进行适当的移动，使得骨架约束点c接近于玉米叶片中间的叶脉点，其计算方法如下：

空间距离权重主要解决点云数据的采样不均匀导致骨架约束点提取偏移的问题，其计算公式为：

其中p为截面相关点集S中的点，q为p点的K近邻集合(p)中的点，表示点p在空间距离上的影响参数，在本文的试验中设置为整个玉米叶片点云数据中点的平均距离。

法向差异权重通过当前点法向与叶片朝向的差异获得，目的是为了使获得骨架约束点位于叶片点云数据的内部，其计算公式为：

点云完整性权重主要解决玉米叶片点云数据存在缺失对骨架约束点提取干扰的问题。对于完整的叶片点云，该权重为0；随着缺失部分的增大，点云完整性权重也会相应的增加。对于缺失的叶片点云，仅利用空间距离权重和法向差异权重提取的骨架约束点在缺失部分产生较大偏差，因此，点云完整性权重作为一个补偿，可通过该权重让骨架约束点朝缺失方向移动来进行校准。具体的，本文将截面相关点集投影到与叶片伸长方向垂直的正交面上并求其凹壳。凹壳作为一个多边形，存在大量缺失部分对应的边长是较长的，可计算凹壳最长的边与总边长的比值来判断缺失情况，具体如下：

其中l表示凹壳所形成的多边形的第条边的长度，为判断点云数据缺失的阈值，在本文的试验中设置为0.1。

使用公式（4）对玉米叶片点云每个截面的相关点集进行运算，能够得到一组玉米叶片的骨架约束点集。

1.5 基于主曲线的骨架优化

主曲线是线性主成分分析的非线性推广，主要通过寻找到一条曲线，使得数据集中的点到该曲线的距离和最小，实质上是寻找通过数据集中间位置的一条光滑曲线。由于前面提取的骨架约束点集是通过对每个子叶片分别计算的，不同子叶片之间的约束点集可能会存在一些不连续的情况，因此可以通过计算上述所有玉米子叶片的骨架约束点集的主曲线，得到更加接近叶脉位置的骨架曲线。

本文使用Verbeek等提出的迭代插入新线段的方式来计算骨架约束点集的主曲线，以优化调整玉米叶片骨架。首先通过PCA 算法从样本集中提取主成分线段，然后依据条件重新划分样本集，并调整对应的Voronoi区域得到新的主成分线段集合，直到满足指定线段数目为止，最后通过贪心算法构建Hamiltonian路径得到最终的主曲线。主要步骤如下：

1）初始化。令=1，将所有子叶片的骨架约束点集作为样本集C，并用其构成初始Voronoi区域V（包含所有骨架约束点），通过对其进行主成分分析，以第一主成分为方向，并从质心开始向其正负方向分别取长度为1.5倍V标准差长度的线段作为初始的主成分线段L，并计算中每个骨架约束点到主成分线段的投影及距离。距离计算公式为：

其中表示点云中的一个点，(（)∊）表示直线段所在直线的参数方程。

2）插入新的主成分线段。首先从骨架约束点集中找到离主成分线段集合 {,...,L}最远的点，同时该点周围有3个以上点（包括3个），用其构建新的样本集C；然后在新的样本集C中计算其第一主成分线段L加入到主成分线段集合中得到 {,...,L}，并重新计算新的Voronoi区域 {,...,V}，直到达到设定的最大线段数目或者使如下的能量函数最小为止：

其中为骨架约束点集的方差，是所有主成分线段的长度之和，表示骨架约束点集中点的数目，c表示 Voronoi区域中的点，L是第条主成分线段。

3）生成主曲线。通过贪心算法对所计算的主成分线段集合构建Hamiltonian路径，并使用2-opt优化算法对得到的路径进行优化，得到最终的主曲线。

1.6 骨架后处理

上述计算骨架约束点集的主曲线生成的骨架可能会在叶片基部和叶尖部产生收缩，故在得到骨架后，需要通过后处理对骨架两端的位置进行校准。将骨架约束点集首尾两点作为骨架的起始点和结束点，搜索骨架两端延长线方向到首尾骨架约束点即可实现该后处理。首尾两端点校准后，可提高骨架提取精度，以及后续表型（如叶长）计算的精度。

2 结果与分析

2.1 玉米叶片点云骨架提取结果

利用本文所提出的基于自适应加权算子和主曲线的玉米叶片点云骨架提取方法，对所选取的玉米叶片点云数据进行骨架提取，结果如图5所示。

对于输入的玉米叶片点云数据，其理想的玉米叶片骨架应该是玉米叶片的叶脉。从图5可以看出，对于常见玉米叶形的点云数据，均能提取出较好的结果；对于存在点云缺失的数据，如叶片6~8，本文方法也能较好地提取出叶片骨架。

图5 玉米叶片骨架提取可视化结果 Fig.5 Visual results of maize leaf skeleton extraction

为了定量评价对于缺失点云数据骨架提取的准确性，对完整叶片点云手动去除部分点进行试验，并引入Hausdroff距离对比去除前后点云所提取骨架差异。Hausdorff距离是一种度量空间中任意两点云之间的距离，其值越小表示两点云的相似程度越高，如果为0则说明两点云完全相同。图6给出了叶片1和叶片2去除部分点前后所提取骨架的可视化差异，它们在完整和缺失情况下所提取骨架的Hausdorff距离和叶长的比值分别为1.3%和0.3%。这些结果表明，本文方法对于存在缺失数据的点云数据仍能提取出较好的骨架，并且与完整的叶片点云提取结果差异不大。

图6 完整点云数据与缺失点云数据提取骨架结果对比 Fig.6 Comparison of skeletons extracted from complete point cloud data and missing point cloud data

为测试该算法对于整株玉米的叶片骨架提取的效果，选取2株玉米植株点云数据进行测试。测试时首先对植株进行器官点云分割，再对植株上的所有叶片进行骨架提取，结果如图7所示。可视化结果表明，利用本文方法所提取的植株上所有叶片的骨架均能准确反映其三维形状结构特征。此外，通过骨架计算的叶长与实测叶长的单株平均绝对百分比误差（MAPE）为1.79%，标准化均方根误差（NRMSE）为1.99%。

图7 两株玉米植株叶片骨架提取结果 Fig.7 Skeleton extraction results of the leaves in two maize plants

2.2 完整性权重和主曲线对骨架提取的影响

本文通过自适应加权算子和主曲线实现了玉米叶片点云数据的骨架提取，下面针对自适应加权算子中完整性权重和主曲线对骨架提取的贡献进行分析。

对于图6中手动删除部分点后所得的有较大的数据缺失的叶片1，图8给出了对其在自适应加权算子中是否使用点云完整性权重所提取的骨架约束点集的对比结果。从图8中可以看出，包含完整性权重的自适应加权算子所提取的骨架约束点（蓝色）在叶片1缺失处更接近叶片中间，而不包含完整性权重的自适应加权算子所提取的骨架约束点（红色）在叶片1缺失处明显偏离叶片中间位置。

图8 有无完整性权重提取的骨架约束点结果 Fig.8 Skeleton constraint point results extracted with completeness weight or not

图9给出了对于图8中提取的骨架约束点集所生成的主曲线的对比结果。从中可以看出，通过骨架约束点集生成的主曲线能够在一定程度上对骨架进行校正，并且使得提取的骨架更加均匀平滑，方便后续应用。然而，对于不包含完整性权重的自适应加权算子所提取的骨架约束点（蓝色），由于骨架约束点在叶片1缺失处明显偏离叶片中间位置，所生成的主曲线（红色）仍然存在较大的偏离。相反的，对于包含完整性权重的自适应加权算子所提取的骨架约束点，所生成的主曲线（黑色）则接近于玉米叶片的中间位置。

图9 有无完整性权重提取的骨架约束点集生成的主曲线 Fig.9 Results of principle curves extracted by skeleton constraint point sets with completeness weight or not

2.3 与基于Laplace算子的骨架提取方法对比结果

为评估本文提出的玉米叶片点云骨架提取方法先进性，用本方法和Cao等提出的利用Laplace算子收缩进行点云提取骨架方法进行对比，如图10所示。由图10可知，本文方法提取的骨架更加呈现曲线形式，在形态上较Laplace方法以点集的形式表示的骨架更为精细；而且对于存在缺失的叶片点云数据，Laplace方法在缺失处的提取的骨架明显存在偏离中心的情况，并且经常出现分叉结构，如图10所示，而本文方法不会出现分叉情况，说明本文方法可以提取更加准确反映玉米叶片三维形状特征的骨架。

图10 两种骨架提取方法的结果对比 Fig.10 Comparison of results of two skeleton extraction methods

2.4 基于提取骨架的叶片表型计算结果

玉米叶片的骨架是计算诸多叶片表型的有效途径，如叶长、叶倾角和叶方位角等。其中，骨架提取精度对叶长的计算影响最大。为此，对比了通过本方法提取骨架计算得到的叶长与实测叶长的差异，如表2所示。本文使用了几种典型的误差度量来测量上述30个叶片的叶长计算误差，包括平均绝对百分比误差（Mean Absolute Percent Error，MAPE）、均方根误差（Root Mean Square Error，RMSE）和标准化均方根误差（Normalized Root Mean Square Error，NRMSE）。本方法计算叶长的MAPE为2.10%，RMSE为2.21cm，NRMSE为2.89%，而对比的Laplace方法计算叶长的MAPE为6.16%，RMSE为6.31cm，NRMSE为8.26%。上述结果表明本方法提取叶片点云骨架可以得到更好的结果，利用本方法提取骨架所计算的叶长具有较高的精度。

表2 本文方法和Laplace方法的叶长计算结果对比 Table 2 Comparison of leaf length calculation results generated by the method in this paper and Laplace method

3 讨 论

由于植物叶片表面光学属性和环境因素的影响，以及部分数据获取手段受获取角度和植物器官自身遮挡的限制，数据缺失在植物叶片点云中普遍存在。如何利用缺失点云数据尽可能高精度地实现叶片骨架提取对于进一步表型信息的提取尤为重要。直接基于点云位置收缩的方法通过计算邻域中心点计算骨架特征点，如分层聚类、Laplace骨架提取、L1中值和平面参数化等方法，其难以避免因点云缺失产生骨架点偏移和叶脉曲线扭曲问题。尽管Wu等针对玉米叶片形态特征对提取的骨架进行矫正，提升了骨架与叶片的重合度和叶尖点位置，但其仍无法解决部分点云缺失产生的骨架偏移问题。叶脉骨架的偏移和扭曲将直接导致提取的叶长参数偏大、骨架驱动的叶片三维建模失真等问题。基于缺失点云的骨架提取是计算机图形学领域的经典问题，但已有研究多侧重于对具有圆柱形对称性的人体、动物和文物开展，或用于解决圆柱形植物枝干和器官的骨架提取，直接应用其解决薄面形叶片的骨架提取时难以达到理想的效果。Jin等在玉米植株点云分割的基础上，通过寻找叶片在不同投影方向上的最小损失，可得到整体光顺的玉米叶片骨架。该方法对于基于缺失点云的表型解析是一个较为理想的解决方案，但所提取骨架与实际叶脉多是在某个投影方向上拟合效果较好，在其他方向仍有偏离。本方法在骨架约束点提取过程中，针对性地设计了点云完整性权重，在一定程度上解决了玉米叶片点云数据缺失的骨架提取问题；同时，通过引入主曲线方法，利用叶片完整部分点云所提取的骨架约束点对缺失部分提取的骨架进行再次矫正，因此本文方法对于点云缺失具有更好的鲁棒性。

长窄形植物叶片的叶脉具有一定的光滑性，为提升骨架提取与叶面的重合度，研究人员常利用样条曲线方法对所提取的骨架特征点进行优化处理。主曲线方法与样条曲线方法具有较高的相似性。样条曲线是通过控制点构建样条曲线，但其形态对控制点尤为敏感，控制点分布的变化会形成不同的样条曲线，导致利用其提取骨架产生较大偏差。而主曲线方法不存在该问题，可以更好的保证提取骨架的准确性。同时本方法提取的骨架可以呈现曲线形式，较传统的点集的形式更为精细。

本文提出的骨架提取方法可实现玉米叶片点云骨架的自动提取，无需后期手动调整，对于处理表型大数据非常重要，尤其是在可实现玉米单株3D数据高通量获取的背景下，自动化的表型解析方法对于推动植物表型组学研究具有重要作用。

4 结 论

1）本文提出了一种基于自适应加权算子和主曲线的玉米叶片点云骨架提取方法，可实现玉米叶片骨架的精确提取。该方法使用的点云完整性权重可以提升点云缺失叶片的骨架提取效果。对比利用完整叶片点云和手动去除部分点形成缺失的叶片点云所提取骨架，二者的Hausdorff距离和叶长的比值低于1.3%。方法中主曲线的引入可以提升骨架的光滑性和骨架点分布的均匀性。

2）利用本方法针对典型的玉米叶片所提取骨架计算叶长与实测叶长的平均绝对百分比误差为2.10%，均方根误差为2.21 cm，标准化均方根误差为2.89%，与同类方法相比可提升叶片表型解析精度。本方法也可很好地应用在整株玉米的表型解析中，通过用其对整株玉米叶片进行测试，所提取骨架计算的叶长与实测叶长的单株平均绝对百分比误差和标准化均方根误差均在2%以内。本方法提取骨架过程中无需手动交互，可实现自动化的骨架提取和表型解析。