CN111667529A - 一种植物点云叶片分割与表型特征测量方法 - Google Patents

Abstract

本发明提供了一种植物点云叶片分割与表型特征测量方法。首先利用多视角立体视觉算法，对盆栽植物进行了精确的点云三维重建。其次使用基于区域和颜色的滤波器和法向量差算法去除了植物点云中的非叶片部分。然后借助曲率信息和多特征区域生长算法得以实现对具有叶片重叠现象的冠层上的单叶片分割。对于每个单独叶片，利用PCA估计了叶片的3D包围盒，叶倾角可计算为包围盒高度方向与作物坐标系的Z轴的夹角，叶长与叶宽即为包围盒的长与宽。本发明能实现对叶片重叠成簇的盆栽植物自动叶片分割，并且对提取叶面积、叶片长、宽、和叶倾角等表型信息有较高的精度与实时性，适用于高通量的叶片表型分析。

Description

一种植物点云叶片分割与表型特征测量方法

技术领域

本发明涉及数字化农业、农业工程及其自动化、植物学等技术领域，特别是涉及一种基于多视角三维重建的植物点云自动叶片分割与表型测量方法。

背景技术

植物“表现型”指的是一组生物性状，这些性状可以将基因型与环境之间的相互作用联系起来，直观反映了植物的生长状态，也是所有育种过程的基础。植物表型是对植物生长发育、耐受性、抗性、结构、生理、产量等复杂性状的综合评价。通常植物表型的定量分析主要包含以下几个方面：叶片特征、茎杆特征、果实性状和根系形态。叶片作为表面积最大的植物器官，是植物光合作用和呼吸作用的主要场所。因此植株的叶面积、叶长、叶宽以及叶倾角都是十分重要的，生物特征随着对植物表型分析需求的日益增多，传统的手工测量方法操作繁冗而且时间较长，并不适用于高通量的表型参数采集。此外采集的精准度也有了更高的要求。

植物表型分析需要精确测量植物的几何特征，数据的采集过程叫做植物的数字化过程。目前植物的数字化方法主要分为两类：接触和非接触式的。前者依靠一些探针接触植物表面，通过记录接触点的空间位置，就可以得到植物的三维模型。接触式的方法非常依赖于手动操作与专业知识，是一个非自动化的枯燥过程。相较于接触式方式，非接触式方式具有更加明显的优势。目前激光雷达、基于结构光和飞行时间法(Time of Flight，简称ToF)的深度传感器以及双目立体视觉匹配算法都广泛应用于植物的三维重建。一方面这些设备较为昂贵，不适用于对成本敏感的农业领域。另一方面，由于成像角度限制，对于叶片严重遮挡的大型植物或叶密度较大的作物效果难以令人满意。此外现有叶片的自动分割方法对于具有高密度的植物冠层不太适用。业界也没有提出较好的一整套针对于植物叶片表型测量的方法。

发明内容

本发明的目的是：针对冠层密度较大的盆栽植物提供一种可靠的植物数字化、自动叶片分割和叶片表型测量的方法，该方法能够较为快速地测量植物叶片数量、叶面积、叶长、叶宽和叶倾角等植物表型参数，适用于高通量的植物表型分析平台。

为了达到上述目的，本发明的技术方案是提供了一种植物点云叶片分割与表型特征测量方法，其特征在于，包括以下步骤：

步骤1：采集n张植物的多视角图像I_n，其中n为正整数。图像序列中相邻两幅图像间应具有较窄的基线和较多相互重叠区域。根据植物大小与冠层密度，n的范围在[50，100]，能获得较好的重建效果。

步骤2：基于获取的图像序列，利用VisualSFM实现运动中三维重建(StructureFrom Motion，简称SFM)，得到稀疏的点云。

步骤2中得到稀疏的点云的具体步骤为：

步骤2.1：利用尺度不变特征变换匹配算法(Scale-invariant FeatureTransform，简称SIFT)对每幅图像进行特征点检测，并在多幅图像中进行特征点配对；

步骤2.2：选择特征点匹配对数最多的相邻两幅图像作为初始图像对，利用随机抽样一致算法(RANdom SAmple Consensus，简称RANSAC)计算其基础矩阵、相机参数和投影矩阵。对已匹配特征点进行集束调整；

步骤2.3：加入新的图像，根据其特征点与已重建三维点的对应关系，利用RANSAC算法求取新增图像的投影矩阵，从而计算新加入的三维点，并对已重建点进行集束调整；

步骤2.4：重复上述过程，直到所有的图像都加入现有的模型，并对所有的三维点进行集束调整，得到稀疏的三维点云。

步骤3：利用多视角立体匹配聚类算法(Clustering Views for Multi-viewStereo，简称CMVS)算法获得稠密的三维点云用集合χ表示，其中P_i是点集中的第i个点。

步骤4：对点云进行预处理，去除背景信息、花盆与离群噪声点。得到的点的集合用Γ表示，其中P_i是点集中的第i个点。

步骤4中对点云预处理的具体步骤为：

步骤4.1：空间区域滤波器，对于点集χ中的每个点P_i，如果同时满足条件：

则保留该点。其中P_i.x、P_i.y和P_i.z分别是第i个点的XYZ坐标值。X_min、X_max、X_min、Y_max和Z_min、Z_max分别是XYZ三个方向上的限定阈值。其取值根据植株在点云中的实际位置确定。

步骤4.2：颜色滤波器，由于植物叶片与茎干边缘易与白色背景形成明显的纹理信息。因此重建生成的点云会在植物边缘形成一圈白色冗余点。如果满足RGB三通道值都大于190则删除该点。

步骤4.3：半径滤波器，对于步骤4.2处理后点集中的每一个点，以该点为圆心，r₁为半径的球形域中寻找邻近点，邻近点的个数用k₁表示。如果k₁≤N₁，就认为该点为离群点。其中N₁为常数，是设置的点数阈值。通常r₁设置为点云中点的平均距离的5～10倍，N₁取值在[5，20]之间。得到的三维点集合用Γ表示，其中P_i是点集中的第i个点。

步骤5：计算点云尺度，与点云缩放。由多视角三维重建方法得到的植物点云与世界坐标系尺度不能唯一确定。需要将点云缩放到与世界坐标系下相同的比例，计算出来的表型参数才有意义。

步骤5中植物点云尺度变换的具体步骤为：

步骤5.1：以花盆为参照物，测量现实中花盆的直径，计作l₁；

步骤5.2：测量植物点云中花盆的直径，计作l₂；

步骤5.3：计算世界坐标系与点云坐标系尺度比例s＝l₂/l₁；

步骤5.4：对点集Γ中的所有点P_i的坐标除以s，将点云坐标系的尺度缩放到与世界坐标系尺度相同。

步骤6：利用法线差分法(Difference of Normals，简称DoN)去除植物的茎干。DoN是一种多尺度滤波操作，通过计算同一点在不同尺度下的法向差，可以判断该点是否在平面上。对于大多数植物来说，叶片相对扁平，同一点在合理尺度范围下的法线方向变化不大。但由于植物茎干近似于圆柱形，不同的尺度法向量会产生较大差异。因此利用法线差可以提取出植物叶片部分。

步骤6中利用DoN去除植物茎干的具体步骤为：

步骤6.1：分别在两个半径尺度上计算点集Γ中的所有点P_i的法向量方向。以P_i为圆心，分别以r₂和r₃为半径，搜索在半径范围内的所有邻近点。以所有邻近点作为支撑点利用PCA进行平面拟合，所得平面法向量就代表该点的法向量。

步骤6.2：计算所有的点P_i的法线差Δ_n，公式如下：

其中，r₂和r₃是计算P_i法向量的两个半径。满足

且r₂＜r₃。n(P_i，r)是点P_i在支撑半径r下的法向量。Δn(P_i，r₂，r₃)是同一点在两种半径尺度下的法线差。实验表明，只有当小半径r₂为点平均距离的2倍以上，同时大半径r₃也不应该超过茎秆直径的时候，所估计出的两个尺度的法向量才有意义。同时，还应保证两半径比值要在可能范围内尽可能地大以加强特征区分度。经过参数调谐，我们发现当r₂为点云中点平均间距的两倍，且r₃等于两倍的r₂时，DoN去除植物茎干效果达到最好。

步骤6.3：由于两个尺度计算出的法向量均为单位向量，而且其夹角位于0至90度之间变化，因此DoN值在0与1.44之间变化。为了便于阈值的选取，我们将所有的DoN计算值归一化到[0，1]区间。取阈值为μ，当|Δ_n(x_i，r₂，r₃)|＜μ时，该点位于叶片部位。反之，当法线差大于阈值μ时，就认为该点属于茎干。一般的，阈值μ的区间在[0.03，0.1]之间。叶片越大茎干越粗的植株选择越大的阈值。

步骤6.4：在DoN处理之后，一些孤立的茎部，特别是叶片附近的端部，可能仍然留在点云中。我们使用欧式聚类删除点云中剩余的茎部分，在此之后得到的点云用集合Ψ表示，其中P_i代表点集中的第i个点。

步骤7：受成像精度的影响，生成的原始植物点云叶片部位也会存在点的波动。部分叶片上的点也被错误的认为是茎干而被移除。为了减少漏报率，对点集Ψ中的所有点，以P_i为球心，r₄为半径在点集Γ中搜寻邻近点，并将搜索到的不在Ψ中的点加入到Ψ点集中。得到完整的没有茎干的植物冠层点云，用集合Φ表示。其中r₄对结果影响不大，一般可以取平均点距的4～10倍。

步骤8：对于步骤7得到的冠层点云需要分割成单个叶片。通过先去除曲率过大的点，随后综合考虑曲率、法向量夹角、点间距离进行区域生长实现了单叶片点云分割。

步骤8中单叶片分割的具体步骤为：

步骤8.1：植物叶片分割的主要难点在于将相互粘连的叶片分割开。针对现实中叶片的重叠情况，将叶片的重叠方式定义为三类。共面相接，两片叶片基本处于同一个平面，叶片边缘部分相互接触；交叉相接，两叶平面法向量差异很大，只有交叉部位少量点粘连在一起；平行相接，两叶片平面近似平行，上下重叠，底部叶片边缘接触于顶部叶片中部。无论哪种粘连方式，叶片的粘连部位点往往具有较大的曲率且集中在0.4～0.6之间。这一步骤中，首先计算冠层中每个点的曲率值c_i，其中i代表第i个点。将c_i＞c_max的点移除，其中c_max为曲率阈值，其取值范围在[0.1，0.2]之间。目的是将粘连的叶片分隔开，得到的点云用集合Θ表示。

步骤8.2：对于点集Θ，利用区域增长算法，将植物冠层分割成单个叶片。算法主要需要考虑两个参数：法向量夹角Δ_n，搜寻半径r₅。其中法向量夹角Δ_n取值一般在20度左右，搜寻半径r₅为平均点距的2倍。分割得到的单叶片点云用Leaf_i，其中i代表第i个叶片。

步骤9：对于每一个单独的叶片点云Leaf_i，先对其进行面片化过分割，将叶片分割成为若干个小面片。其次利用罗德里格旋转公式将小面片分别投影到二维平面上。接着对二维平面上的小面片进行德劳内三角化，通过计算每个小三角面积和得到小面片面积。最后累加所有小面片面积，就可以得到叶面积s_i。为了消除叶边缘部位噪点产生的冗余三角形对面积计算精度的影响，该步骤中还引入了三角形边长阈值，作为判断冗余三角形的条件。

步骤9叶面积计算的具体步骤为：

步骤9.1：面片过分割，首先利用迭代主成分分析法(Iterative PrincipleComponent Analysis，简称IPCA)计算点云Leaf_i中每个点的空间特征(包括法向量和平滑度)。选择一系列平滑度较大区域的点作为种子点开始面片的初步生长，从种子点开始的生长包含三个条件：1)候选点的距离与种子点的距离必须小于一定的距离r₆；2)候选点与种子点的法向量夹角小于阈值θ；3)候选点到种子点所在的平面之间的距离小于某个阈值σ₁。同时符合上述三个条件的候选点将被纳入当前面片，如果候选点不被归属于任何种子点，则其变为一个新的种子点重新参与计算。最后对已经初步划分的所有面片使用局部K均值聚类(Local K-means Clustering)以进一步改善面片的边缘结构。面片的粒度主要与参数K与r₇相关。K代表用于计算点的空间特征的邻域点数量，r₇代表了允许的面片最小半径。其中r₆取值为平均点距的10～20倍；θ取值为23左右；σ₁取值为平均点距的10～20倍；K取10～20；r₇取r₆的2倍大小。得到的小面片用facet_i表示，其中i代表该叶片上分割出来的第i个小面片。

步骤9.2：面片旋转投影到二维平面。对于每一个小面片facet_i，都需要将其投影到二维平面以计算其面积。我们选取投影后的法向量为q＝[0，0，1]^T，投影前的面片法向量p由对面片上所有的点做PCA估计得到。根据罗德里格旋转公式(Rodrigues RotationFormula)，可以计算出一个3×3的旋转矩阵R(θ)：

R(θ)＝I+Ωsin+Ω²(1-cosθ)

其中，I为3-by-3的单位矩阵。ω为单位长度的旋转轴向量：

Ω和ω的关系如下：

为旋转角，即向量p和q之间的夹角，可由向量点乘得到。将旋转矩阵R(θ)与当前面片所有点的XYZ坐标进行相乘，得到的新面片几乎平行于XOY平面。由于面片并不严格等于一个平面，因此旋转后的所有点的z-axis坐标会在0左右微小波动，直接令所有z坐标为0，就实现了将面片旋转投影到二维平面XOY上。

步骤9.3：德劳内三角化(Delaunay triangulation)，对步骤9.2中得到的二维小面片进行德劳内三角化。即将面片上所有的点为顶点，构成三角形网格面。通过计算所有三角形面积和，就能得到面片面积。

所述步骤9.3中，需要去除由叶边缘噪声点产生的冗余三角形，这些三角形往往具有细长的形态特征，而且面积较大，对叶面积计算精度影响较大。需要去除该部分三角形的影响。在每个面片被Delaunay三角化算法剖分后，会将形成的所有三角形的边长进行排序。此时用r₈表示边长的中值，并定义一个阈值倍数为λ。如果三角形任意边长都小于λr₈，则该三角形正常参与面积计算；否则从面片中删除该三角形。该方法能够较好地适应于不同品种的植物叶片和点云，同时能够避免离群点的影响。通过对不同植物不同叶片定量分析，λ＝3时叶面积计算精度最高。

步骤9.4：计算同一叶片上所有面片facet_i的面积和，得到的面积即为该叶片面积s_i。

步骤10：计算所有单叶片Leaf_i的叶长I_l，叶宽L_w和叶倾角α。考虑到PCA能用于计算数据聚类的正交主成分方向，而叶长方向、叶宽方向、叶法向量方向正好两两正交，所以PCA非常适合用于这三个向量方向的估计。叶片一般具有细长型片状结构，所以针对单叶上所有点进行PCA计算出的具有最大特征值的主成分可以理解为代表了叶长方向，而第二主成分则代表了叶宽方向，最小的主成分方向为叶法向量方向。通过将叶片点云质心、PCA坐标系分别与坐标原点和世界坐标系对齐，寻找叶片在XYZ方向上的最大最小值，构建3D包围盒。包围盒的长代表叶长，包围盒的宽代表叶宽。叶法线方向与z轴的夹角代表叶倾角。

步骤10中计算单叶片叶长，叶宽和叶倾角的具体步骤为：

步骤10.1：首先计算单叶片Leaf_i的质心点c＝[x_c，y_c，z_c]^T；

步骤10.2：对叶片点集进行PCA计算可得到大小呈递减趋势的三个特征值λ₁，λ₂，λ₃，其对应的特征向量分别表示为α₁，α₂，α₃，将特征向量按列排列可得到3-by-3的特征向量矩阵A＝[α₁，α₂，α₃]。此时将叶片坐标系变换到世界坐标系的齐次旋转平移矩阵V可表示为：

V是一个4-by-4的矩阵，可以实现齐次的点云坐标值变换。点云中任意某一点在变换前的齐次坐标可定义为x_i＝[x_i，y_i，z_i，1]^T，被旋转平移矩阵左乘后可得到变换后的齐次点坐标

有如下的变换对实现叶片上的点在两个坐标系下的相互转换：

步骤10.3：经过旋转平移变换，遍历所有叶片点并分别记录点集在XYZ坐标轴上达到的最大最小值即可建立唯一的3D包围盒。其中叶长L_l为包围盒的长，叶宽L_w为包围盒的宽。

步骤10.4：PCA计算得到的第三分量的特征向量为α₁，可以近似认为是叶片腹面方向的法向量。计算α₁与Z轴方向的夹角，即为叶倾角α。

由于采用了上述的技术方法，本发明与现有技术相比，具有以下优点：

1)在基于MVS方法得到的高精度植物三维点云上，通过一种改进的DoN算法实现了作物冠层点云中的叶片与非叶部分(比如茎干)的分离，提取了只包含叶片部分的纯净冠层点云。方法能保留冠层叶片完整的结构信息，能消除单视角叶片遮挡导致的叶片结构残缺问题，适用于高精度植物叶片重建。

2)设计了一种全自动的单叶片(individual leaf segmentation)分割算法，首先使用点云曲率特征去除叶片相互重叠的边界，然后基于多特征区域生长方法实现了高精确度的单叶片分割。该分割算法有能力处理叶片重叠成簇的情况。

3)在叶片分割结果上，提出了利用面片过分割和德劳内三角化(Delaunaytriangulation)精确计算植株单叶面积的方法。方法同样适用于不同结构且具有较大弯曲度和不平整度的叶片且达到了较高精度。

4)设计了一种利用PCA产生叶片空间最小包围盒的方法提取叶片长、宽、和叶倾角等表型信息。方法具有良好的精度与实时性，适用于高通量的叶片表型分析。

附图说明

图1为基于多视角三维重建的数字化植物自动叶片分割与叶片表型分析方法流程图。

图2为实例植物部分多角度图像与利用VisualSFM生成的稀疏三维点云，其中，(a)为多角度图片的相机位置分布和生成的稀疏点云；(b)～(e)是抽取的4帧用于重建的实例植物图像。

图3为CMVS生成的稠密点云，其中，(a)是含有背景信息的植物点云，(b)、(c)、(d)是去除背景信息与花盆的植物点云，其中，(b)是前视图，(c)是俯视图，(d)是侧视图。

图4为实例植物叶片茎干分离结果，其中，(a)是经过DoN之后保留的点云；(b)是经过整体茎干去除处理之后的冠层点云结果侧视图；(c)是去除的茎干部分的侧视图；(d)是(b)的俯视图。

图5为在三种叶片重叠情况下的叶分割原理及结果展示。第1列到第3列分别显示了三种不同的重叠情况。(a1)为两个叶片在边界处重叠。在(a2)中，我们进一步展示了重叠区域，分别使用点和网格来显示该区域中两片叶子的三维方向。网格平面代表叶平面。在(a3)中，我们从两片叶子中计算出每一点的归一化曲率，并用不同的颜色表示。与叶表面的点相比，重叠边界周围的点具有较高的曲率值。去除高曲率点后，利用区域生长技术对两片叶片进行分割，最终得到(a1)的单叶片分割结果(a4)，不同的叶片用不同的颜色表示。(b1)显示了三个叶片交叉重叠的情况，(b2)是其交叉情况的模拟图，(b3)用不同颜色展示曲率值，(b4)是(b1)的分割结果。(c1)显示了两个叶片，从两个平行的层相互接触。(c2)是其交叉情况的模拟图，(c3)用不同颜色展示曲率值，(c4)是分割结果。

图6为单叶片分割的结果图，其中，(a)为分割结果俯视图，(b)为分割结果侧视图，其中不同的颜色代表分割出来的不同叶片。

图7为对叶片面片化，并对有代表性的小面片进行Delaunay三角剖分后的结果示意图，其中，(a)是叶片点云进行面片过分割后的示意图，其中原始点云在下，面片结果在上，我们选取了分别位于边缘和中心的两个小面片进行分析。(b1)与(b2)是对(a)中选取的两个小面片经过旋转投影后进行德劳内三角化剖分的结果图。(c1)与(c2)分别对应(b1)与(b2)中通过边长阈值去除了离群三角形后的结果图，其中深黑色部分是去除的离群三角形，浅灰色部分是保留的面片部分，计算叶面积时将只计算浅灰色部分。

图8演示了提出的叶片表型特征提取方法，其中，(a)表示植物冠层在世界坐标下的原始点云，其中z轴沿重力方向；(b)中选取了冠层中的一片叶子，中部白色圆点是叶片点云的质心。(c)显示了建立的叶片坐标系，其中，利用主成分分析法计算叶片点云三个正交主成分。在(d)中，我们使用旋转平移矩阵将叶坐标系与世界坐标系对齐。(e)为三个坐标中求最大值和最小值构造的三维边界框，叶长和宽度分别为边界框的两个水平边长。(f)表示通过逆旋转平移变换将叶片与包围盒旋转到原始位置，叶倾角为叶片法线方向与世界坐标系中z轴的夹角。

图9为用基于PCA的3D包围盒求取植物叶片表型特征效果图，其中，(a)(b)(c)分别是带有叶片包围盒与法向量的竹芋冠层点云的俯视图，前视图和侧视图。

具体实施方式

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

如图1所示，本实施例公开了一种植物点云叶片分割与表型特征测量方法，其包括以下步骤：

步骤1：对盆栽植物竹芋(Maranta arundinacea)采集了102张窄基线和较多相互重叠区域且分辨率为4000×3000的图像。其部分视角图像如图2(b)～(e)所示。

步骤2：基于获取的图像序列，利用VisualSFM软件进行了特征点选取和特征匹配后得到了含有23458个特征点的稀疏点云，如图2(a)。

步骤3：利用CMVS实现稠密点云的构建，得到了含有1400333个点的稠密点云集合，如图3(a)所示。

步骤4：对点云进行预处理，去除背景信息、花盆与离群噪声点。在空间滤波器部分，选取的XYZ三个方向上的限定阈值分别为：X_min＝-11.5、X_max＝1.7，Y_min＝14.6、Y_max＝21.4和Z_min＝-3.0、Z_max＝10.2。其取值根据植株在点云中的实际位置确定。得到的去除背景信息点云的点数为701253。由于植物叶片与茎干边缘易与白色背景形成明显的纹理信息。因此重建生成的点云会在植物边缘形成一圈白色冗余点。如果满足RGB三通道值都大于190则删除该点。得到去除白边的点云点数为689014。通过计算，植物点云的平均点距为0.03m，最后利用半径滤波器，设置搜寻半径r₁＝0.2，邻近点的个数k₁＝15，去除离群点。最终得到含有688374个点的去除背景信息的植物点云。得到的三维点集合用Γ表示，设P_i是三维点集合Γ中的第i个点。如图3(a)(b)(c)所示分别为前视图，俯视图和侧视图。

步骤5：计算点云尺度，与点云缩放。由多视角三维重建方法得到的植物点云与世界坐标系尺度不能唯一确定。需要将点云缩放到与世界坐标系相同比例。以花盆为参照物，测量现实中花盆的直径l₁＝0.2m；测量植物点云中花盆的直径l₂＝3.03m。故世界坐标系与点云坐标系尺度比例s＝l₂/l₁＝15.15。对三维点集合Γ中的所有点P_i的坐标除以s，将点云坐标系的尺度缩放到与世界坐标系尺度相同。此时植物平均点间距为0.002m。

步骤6：利用DoN去除植物茎干部分。分别在两个尺度上计算三维点集合Γ中的所有点P_i的法向量方向。以P_i为圆心，分别以r₂＝0.004和r₃＝0.008为半径，搜索在半径范围内的所有邻近点。以所有邻近点作为支撑点利用PCA进行平面拟合，所得平面法向量就代表该点的法向量。法线差的模长归一化之后，取阈值μ＝0.031，删除所有法线差大于阈值的点。最终得到的剩余点云点数为126940，如图4(a)所示。在DoN处理之后，一些孤立的茎部，特别是叶片附近的端部，可能仍然留在点云中。我们使用欧式聚类去除包含少于200个点的孤立聚类，这步的目的是去除点云中仍残留的小茎段和聚成团状的噪声区域，在此之后得到的点云用集合Ψ表示，设P_i代表集合Ψ中的第i个点，共有125870个点。

步骤7：虽然DoN处理后对植物叶片部分的精度接近100％，但召回率仅为19.4％。因此，有必要在DoN处理后对点云补充叶片结构，在保持精确度的情况下提升召回率。对集合Ψ中的所有点，以P_i为球心，r₄＝0.01为半径在集合Ψ中搜寻邻近点，并将搜索到的不在集合Ψ中的点加入到集合Ψ中。得到完整的没有茎干的植物冠层点云，用集合Φ表示，共有637353个点，召回率提高到97.4％。如图4(b)(c)(d)所示，其中(b)为只含冠层叶片点云的侧视图，(d)为俯视图，(c)是去除的茎干部分点云。

步骤8：进行单叶片分割，首先计算冠层点云每个点的曲率值c_i，其中，i代表冠层点云中的第i个点。将c_i＞c_max的点移除，其中，c_max＝0.01。目的是将粘连的叶片分隔开，得到的点云用集合Θ表示。然后利用区域增长算法，将植物冠层分割成单个叶片。算法主要的两个参数：法向量夹角Δ_n＝20、搜寻半径r₅＝0.05。分割得到的单叶片点云用Leaf_i，其中，i代表第i个叶片。植株原始冠层中共有23个叶片，经过算法分割后最终得到21片可进行表型分析的叶片，正确率为91.3％。如图6(a)(b)所示，为分割的结果图。

步骤9：对步骤8中的每个叶片分别进行面片化过分割，设置的主要参数为：候选点的距离与种子点的距离必须小于一定的距离r₆＝0.03m；候选点与种子点的法向量夹角小于阈值θ＝23；候选点到种子点所在的平面之间的距离小于某个阈值σ₁＝0.025m；用于计算点的空间特征的邻域点数量K＝15；允许的面片最小半径r₇＝0.1。最终每个叶片平均被分割为100片左右的小面片。得到的小面片用facet_i表示，其中，i代表该叶片上分割出来的第i个小面片。如图7(a)展示的面片化效果。

得到的每个小面片都需要将其投影到二维平面以计算其面积。我们选取投影后的法向量为q＝[0，0，1]^T，投影前的面片法向量p由对面片上所有的点做PCA估计得到。根据罗德里格旋转公式实现将面片旋转投影到二维平面XOY上。

最后对上述二维小面片进行德劳内三角化。即将面片上所有的点为顶点，构成三角形网格面。如图7(b1)(b2)所示，通过计算所有三角形面积和，就能得到面片面积。该步骤中取λ＝3去除由叶边缘噪声点产生的冗余三角形，如图7(c1)(c2)所示。将同一片叶片上的所有面片面积相加，就得到了该叶片面积。

步骤10：计算单叶片叶长，叶宽和叶倾角。首先计算单叶片Leaf_i的质心点c＝[x_c，y_c，z_c]^T，如图8(b)所示；对叶片点集进行PCA计算可得到大小呈递减趋势的三个特征值λ₁，λ₂，λ₃，其对应的特征向量分别表示为α₁，α₂，α₃。将特征向量按列排列可得到3-by-3的特征向量矩阵A＝[α₁，α₂，α₃]。此时将叶片坐标系变换到世界坐标系的齐次旋转平移矩阵V可表示为：

V是一个4-by-4的矩阵，可以实现齐次的点云坐标值变换；c表示单叶片Leaf_i的质心点。点云中任意某一点在变换前的齐次坐标可定义为x_i＝[x_i，y_i，z_i，1]^T，被旋转平移矩阵左乘后可得到变换后的齐次点坐标

如图8(c)(d)所示过程，有如下的变换实现叶片上的点在两个坐标系下的相互转换：

经过旋转平移变换，遍历所有叶片点并分别记录点集在XYZ坐标轴上达到的最大最小值即可唯一地建立三维包围盒，其中叶长L_l为三维包围盒的长，叶宽L_w为三维包围盒的宽，如图8(e)所示。

PCA计算得到的第三分量的特征向量为α₁，可以近似认为是叶片腹面方向的法向量。计算α₁与Z轴方向的夹角，即为叶倾角α，如图8(f)所示。图9是对实例植物所有叶片计算3D包围盒的效果图，其中(a)(b)(c)分别是俯视图，前视图和侧视图。

从定量结果上看，本发明可对盆栽植物实现多视角三维重建和数字化植物自动叶片分割与叶片表型分析。并且达到了较高的计算精度，叶片分割成功率达到91.3％；叶面积平均计算精度达到98.48％；叶长与叶宽的平均测量误差在4％以下，叶倾角最大误差不超过3.9度，适用于高通量植物表型分析。

Claims (10)

1.一种植物点云叶片分割与表型特征测量方法，其特征在于，包括以下步骤：

步骤1：采集n张植物的多视角图像I_n组成图像序列，其中，n为正整数，图像序列中相邻两幅多视角图像间具有较窄的基线和较多相互重叠区域；

步骤2：基于获取的图像序列，利用VisualSFM实现运动中三维重建得到稀疏的点云；

步骤3：利用多视角立体匹配聚类算法获得稠密的三维点云χ，设P_i是三维点云χ中的第i个点；

步骤4：对三维点云χ进行预处理，用空间滤波器、颜色滤波器和半径滤波器，去除背景信息、花盆与离群噪声点，得到的三维点云Γ，设P_i是三维点云Γ中的第嘻个点；

步骤5：计算三维点云Γ缩放尺度，将三维点云Γ缩放到与世界坐标系相同比例；

步骤6：利用法线差分法去除植物茎干部分获得点集Ψ，法线差分法是一种多尺度滤波操作，通过计算同一点在不同尺度下的法向差，判断该点是否在平面上，利用法线差分法提取出植物叶片部分；

步骤7：叶片点云回填

为了减少漏报率，对点集Ψ中的所有点，以P_i为球心、r₄为半径在点集Ψ中搜寻邻近点，并将搜索到的不在点集Ψ中的点加入到点集Ψ中，得到完整的没有茎干的植物冠层点云Φ；

步骤8：将步骤7得到的植物冠层点云Φ分割成单个叶片，通过先去除曲率过大的点，随后综合考虑曲率、法向量夹角、点间距离进行区域生长实现单叶片点云分割，获得单独的叶片点云；

步骤9：对于步骤8得到的每一个单独的叶片点云Leaf_i，先对其进行面片化过分割，将叶片分割成为若干个小面片，其次利用罗德里格旋转公式将小面片分别投影到二维平面上，接着对二维平面上的小面片进行德劳内三角化，通过计算每个小三角面积和得到小面片面积，最后累加所有小面片面积，得到叶面积s_i，为了消除叶边缘部位噪点产生的冗余三角形对面积计算精度的影响，在每个面片被德劳内三角化算法剖分后，将形成的所有三角形的边长进行排序，用r₈表示边长的中值，并定义一个阈值倍数为λ，如果三角形任意边长都小于λr₈，则该三角形正常参与面积计算；否则从面片中删除该三角形；

步骤10：计算所有单独的叶片点云Leaf_i的叶长L_l、叶宽L_w和叶倾角α，针对单叶片上所有点进行PCA计算出的具有最大特征值的主成分理解为代表了叶长方向，而第二主成分则代表了叶宽方向，最小的主成分方向为叶法向量方向，通过将叶片点云质心、PCA坐标系分别与坐标原点和世界坐标系对齐，寻找叶片在XYZ方向上的最大最小值，构建3D包围盒，包围盒的长代表叶长，包围盒的宽代表叶宽，叶法线方向与z轴的夹角代表叶倾角。

2.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤1中，根据植物大小与冠层密度，采集50～100张图像，能获得较好的重建效果。

3.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤2中，利用VisualSFM实现运动中三维重建得到稀疏的点云具体包括以下步骤：

步骤201、利用尺度不变特征变换匹配算法对每幅图像进行特征点检测，并在多幅图像中进行特征点配对；

步骤202、选择特征点匹配对数最多的相邻两幅图像作为初始图像对，利用随机抽样一致算法计算其基础矩阵、相机参数和投影矩阵，对已匹配特征点进行集束调整；

步骤203、加入新的图像，根据其特征点与已重建三维点的对应关系，利用RANSAC算法求取新增图像的投影矩阵，从而计算新加入的三维点，并对已重建点进行集束调整；

步骤204、重复步骤201至步骤203，直到所有的图像都加入现有的模型，并对所有的三维点进行集束调整，得到稀疏的三维点云。

4.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤4中，空间滤波器的具体使用方法为：对于三维点云χ中的每个点P_i，如果同时满足条件：

则保留该点，其中P_i.x、P_i.y和P_i.z分别是第i个点的XYZ坐标值；X_min、X_max，Y_min、Y_max和Z_min、Z_max分别是XYZ三个方向上的限定阈值，其取值根据植株在点云中的实际位置确定；

颜色滤波器具体使用方法为：由于植物叶片与茎干边缘易与白色背景形成明显的纹理信息，因此重建生成的三维点云χ在植物边缘形成一圈白色冗余点，如果满足RGB三通道值都大于190则删除该点；

半径滤波器具体使用方法为：对于颜色滤波器处理后点集中的每一个点，以该点为圆心，r₁为半径的球形域中寻找邻近点，邻近点的个数用k₁表示，如果k₁≤N₁，就认为该点为离群点，其中，N₁是设置的点数阈值，得到的三维点云Γ。

5.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤5中，计算计算三维点云Γ缩放尺度的具体步骤为：

步骤501：以花盆为参照物，测量现实中花盆的直径，计作l₁；

步骤502：测量植物点云中花盆的直径，计作l₂；

步骤503：计算世界坐标系与点云坐标系尺度比例s＝l₂/l₁；

步骤504：对三维点云Γ中的所有点P_i的坐标除以s，将点云坐标系的尺度缩放到与世界坐标系尺度相同。

6.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤6中利用法线差分法去除植物茎干部分获得点集Ψ具体包括以下步骤：

步骤601、分别在两个尺度上计算三维点云Γ中的所有点P_i的法向量方向，以P_i为圆心，分别以r₂和r₃为半径，搜索在半径范围内的所有邻近点，以所有邻近点作为支撑点利用PCA进行平面拟合，所得平面法向量就代表该点的法向量；

步骤602、计算所有的点P_i的法线差Δ_n，如下式(1)所示：

式(1)中，r₂和r₃是计算P_i法向量的两个半径，满足

r₂＜r₃；n(P_i，r)是点P_i在支撑半径r下的法向量；Δn(P_i，r₂，r₃)是同一点在两种半径尺度下的法线差；只有当小半径r₂为点平均距离的2倍以上，同时大半径r₃也不超过茎秆直径的时候，所估计出的两个尺度的法向量才有意义；同时，还应保证两半径比值在可能范围内尽可能地大以加强特征区分度；

步骤603、由于两个尺度下计算出的法向量均为单位向量，而且其夹角位于0至90度之间变化，因此法线差分法值在0与1.44之间变化，为了便于阈值的选取，将所有的法线差分法计算值归一化到[0，1]取阈值为μ，当|Δ_n(x_i，r₂，r₃)|＜μ时，该点位于叶片部位，反之，当法线差大于阈值μ时，就认为该点属于茎干，叶片越大茎干越粗的植株选择越大的阈值；

步骤604、在法线差分法处理之后，使用欧式聚类删除点云中剩余的茎部分，在此之后得到的点云用点集Ψ表示，设P_i代表集合Ψ中的第i个点。

7.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤8中，单叶片分割具体包括以下步骤：

步骤801、将叶片的重叠方式定义为三类：共面相接，两片叶片基本处于同一个平面，叶片边缘部分相互接触；交叉相接，两叶平面法向量差异很大，只有交叉部位少量点粘连在一起；平行相接，两叶片平面近似平行，上下重叠，底部叶片边缘接触于顶部叶片中部；首先计算植物冠层点云Φ中每个点的曲率值c_i，其中i代表第i个点，将c_i＞c_max的点移除，c_max为曲率阈值，得到的点云用集合Θ表示；

步骤802、对于点集Θ，利用区域增长算法，将植物冠层分割成单个叶片。算法主要需要考虑两个参数：法向量夹角Δ_n、搜寻半径r₅，分割得到的单叶片点云用Leaf_i表示，其中，i代表第i个叶片。

8.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤9中，面片过分割具体包括以下步骤：首先利用迭代主成分分析法计算点云Leaf_i中每个点的空间特征，选择一系列平滑度较大区域的点作为种子点开始面片的初步生长，从种子点开始的生长包含三个条件：

条件1)候选点的距离与种子点的距离必须小于一定的距离r₆；

条件2)候选点与种子点的法向量夹角小于阈值θ；

条件3)候选点到种子点所在的平面之间的距离小于某个阈值σ₁；

同时符合上述三个条件的候选点被纳入当前面片，如果候选点不被归属于任何种子点，则其变为一个新的种子点重新参与计算；

最后对已经初步划分的所有面片使用局部K均值聚类以进一步改善面片的边缘结构；

面片的粒度与参数K与r₇相关，K代表用于计算点的空间特征的邻域点数量，r₇代表了允许的面片最小半径，其中，r₆取值为平均点距的10～20倍；σ₁取值为平均点距的10～20倍；K取10～20；r₇取r₆的2倍大小，得到的小面片用facet_i表示，其中，i代表该叶片上分割出来的第i个小面片。

9.如权利要求1所述的一种植物点云叶片分割与表型特征测量方法，其特征在于：所述步骤9中，将小面片旋转投影到二维平面包括以下步骤：

对于每一个小面片facet_i，将其投影到二维平面以计算其面积，选取投影后的法向量为q＝[0，0，1]^T，投影前的面片法向量p由对面片上所有的点做PCA估计得到，根据罗德里格旋转公式计算出一个3×3的旋转矩阵R(θ)：

R(θ)＝I+Ωsinθ+Ω²(1-cosθ) (2)

式(2)中，I为3-by-3的单位矩阵；θ表示旋转角，即向量p和q之间的夹角，由向量点乘得到；ω为单位长度的旋转轴向量：

Ω和ω的关系如下式(3)所示：

将旋转矩阵R(θ)与当前面片所有点的XYZ坐标进行相乘，得到的新面片几乎平行于XOY平面，由于面片并不严格等于一个平面，因此旋转后的所有点的z-axis坐标会在0左右微小波动，直接令所有z坐标为0，实现了将面片旋转投影到二维平面XOY上。

10.如权利要求1所述的一种基于多视角三维重建的植物点云自动叶片分割与表型测量方法，其特征在于：所述步骤9中，计算单叶片叶长、叶宽和叶倾角的具体步骤为：

步骤901、首先计算单叶片Leaf_i的质心点c＝[x_c，y_c，z_c]^T；

步骤902、对叶片点集进行PCA计算得到大小呈递减趋势的三个特征值λ₁，λ₂，λ₃，其对应的特征向量分别表示为α₁，α₂，α₃，将特征向量按列排列得到3-by-3的特征向量矩阵A＝[α₁，α₂，α₃]，此时将叶片坐标系变换到世界坐标系的齐次旋转平移矩阵V如下式(4)所示：

V是一个4-by-4的矩阵，可以实现齐次的点云坐标值变换，，点云中任意某一点在变换前的齐次坐标定义为x_i＝[x_i，y_i，z_i，1]^T，被旋转平移矩阵左乘后得到变换后的齐次点坐标

有下式(5)所示的变换对实现叶片上的点在两个坐标系下的相互转换：

步骤903、经过旋转平移变换，遍历所有叶片点并分别记录点集在XYZ坐标轴上达到的最大最小值即可唯一地建立三维包围盒，其中叶长I_l为三维包围盒的长，叶宽L_w为三维包围盒的宽；

步骤904、PCA计算得到的第三分量的特征向量为α₁，可以近似认为是叶片腹面方向的法向量，计算α₁与Z轴方向的夹角，即为叶倾角α。

Images