[发明专利]面向激光点云数据的阔叶树真实叶片建模与形变方法

摘要

本发明公开了一种面向激光点云数据的阔叶树真实叶片建模与形变方法，首先从扫描获取的叶子点云中采用多项式拟合的方法获取精确的树叶边界，并结合计算机图形学的方法计算叶面的主叶脉；其次采用双三次广义张量积Bezier曲面对叶面点云数据进行拟合，进而去除由于风吹抖动造成的扫描误差并消除噪声点云；再次，根据固体力学受力模型，定义叶脉与叶肉不同的力学属性，并构造非线性有限元的受力形变方程，计算施加不同载荷力后的叶子形态，进而模拟树叶在自然环境中真实的形变。本发明利用计算机自动分析获取林学指标，能准确描述不同林分条件下的活立木动态生长变化的叶面积指数，从实验结果看，此方法是切实可行并且高效的，并为林学参数的精确求取打下基础。

主权项

1.一种面向激光点云数据的阔叶树真实叶片建模与形变方法，其特征在于：首先运用高精度的激光扫描仪获取不同阔叶树种的三维叶面点云数据；其次从扫描获取的叶子点云中采用多项式拟合的方法获取精确的树叶边界，并结合计算机图形学的方法计算叶面的主叶脉；接着采用双三次广义张量积Bezier曲面对叶面点云数据进行拟合，进而去除由于风吹抖动造成的扫描误差并消除噪声点云点；再次，根据固体力学受力模型，定义叶脉与叶肉不同的力学属性，并构造非线性有限元的受力形变方程，计算施加不同载荷力后的叶子形态，进而模拟树叶在自然环境中真实的形变；该方法包括如下步骤：(1)数据获取(11)扫描获取整株树木，提取叶面点云的两个端点，p_e＝(x_e，y_e，z_e)^T，p_s＝(x_s，y_s，z_s)^T，分别认为p_e是叶尾点，p_s是叶脉的顶端，通过计算p_e和p_s之间的连线L₁来确定主叶脉上的点，L₁：p＝p_e+t×(p_s‑p_e)，t为连线的自变量，同时确定与矢量相垂直的法向量将主叶脉L₁等间隔n等份后取线上面的n+1个等分点，这些等分点与法向量构成了叶子宽度的n+1条扫描线L_2，i，其中i＝1，2，3...，n+1，定位求取L_2，i左右两端的端点p_2，li和p_2，ri，即可以获取叶面在不同部分的边缘点；(12)对步骤(11)的多条L2，i扫描线，取每条扫描线的左右两端的端点p2，li和p2，ri，分别记为p2，li：(xli，yli，zli)和p2，ri：(xri，yri，zri)，其中i＝1，2，3...，n+1；除去两侧的两条扫描线，若某条扫描线的左端点相对于相邻两条扫描线的左端点向主叶脉内凹，且右端点相对于相邻两条扫描线的右端点也向主叶脉内凹，则认为该扫描线存在遮挡或由误差造成，剔除该扫描线；(13)基于步骤(12)，剔除完所有存在遮挡或由误差造成的扫描线后，对所有扫描线按序重新编号为p2，lj和p2，rj，左右两端的端点分别记为p2，lj：(xlj，ylj，zlj)和p2，rj：(xrj，yrj，zrj)；最终得到：左端点的集合为p2l＝{(xl1，yl1，zl1)，(xl2，yl2，zl2)，...，(xlj，ylj，zlj)，...}；右端点的集合为p2r＝{(xr1，yr1，zr1)，(xr2，yr2，zr2)，...，(xrj，yrj，zrj)，...}左端点的集合和右端点的集合构成了叶面的初始扫描轮廓；对于一片叶子的左/右边部分轮廓，采用不同投影面拟合与求交的方法定位真实叶子的边缘；具体步骤如下：(131)对于左/右边部分轮廓上的边缘点p2，j＝(x，y，z)，将(xi，yi，zi)i＝1，2，3...，n作为叶子几何边缘的参数，运用多项式曲线拟合的方法，把边缘点y值作为输入参数，计算拟合系数来求取相应的x′和z′，公式为：整体算式表示为：x≈x′＝vx(y)＝vx1yn+vx2yn‑1+vx3yn‑2+...+vxn′‑1y+vxn′z≈z′＝vz(y)＝vz1yn+vz2yn‑1+vz3yn‑2+...+vzn′‑1y+vzn′其中，vx1vx2vx3...vxn′‑1vxn′为计算得到的多项式拟合系数；(132)得到新的叶面边缘点为P′edge＝{x′l，yl，z′l；x′r，yr，z′r}，从而得到平滑和无偏差的叶子边缘点；在这一部分中，也可以理解为运用了多项式拟合分别把叶子边缘点投影到X‑Y平面和Y‑Z平面上，接着对两个投影面求交后，从而定位真实的拟合后的叶子边缘；(2)叶面拟合与叶面重建在传统的张量积Bezier曲面的基础上，提出了双三次广义张量积Bezier曲面拟合方法；广义张量积Bezier曲面公式为：其中，和分别代表按照参数m、n等间隔分割：取m＝n＝3，得到U，V{0≤U≤1；0≤V≤1}是采样值，采用精度根据m、n的大小设置，D_i，_j和E_i，j分别为处的x方向的目标偏导矢和y方向的目标偏导矢；张量积得Bezier曲面，在单只叶面以4*4的16个控制点，从每条扫描线中依次选取4个点，分别为p0，j，p1，j，p2，j，p3，j，再选择相邻的四条线段，共取16个点构成网格的16个点，把这16个点代入广义张量积Bezier曲面公式pi，j中，进而采用广义双三次Bezier曲面对其拟合，最后，重复上述过程直到所有扫描线上的扫描点都遍历一遍；(3)叶面受力形变(31)用Ω来表示叶面形变前的状态，p∈Ω是叶面上的某个点p：(p_x，p_y，p_z)；用Γ表示树叶形变后的状态且p′∈Γ是形变后树叶上的某个点p′(p′_x，p′_y，p′_z)，与p对应；根据Ω，用位移u来描述Γ，即u：Ω→Γ，p′＝p+u(p)，其中：设p₁，p₂∈Ω，p′₁，p′₂∈Γ，形变前p₁，p₂之间的空间位移d表示为：d＝p₂‑p₁，同理，形变后p′₁，p′₂之间的空间距离位移d′：d′＝p′₂‑p′₁，由上式可得：其中：I为单位矩阵，为形变梯度，用F来表示，即通过d，d′之间的变化来表示其形变：|d′|2‑|d|2＝d′Td′‑dTd＝dT(FTF‑I)d根据弹性力学知识，可由上式得出格林应变E∈R3×3为：从上式可看出E是对称矩阵，对于以四面体为基本单元组成的叶面模型而言，采用线性插值函数来表示每个节点的位移向量u＝(u，v，w)；对于一个由4个节点组成的四面体，每个节点有3个自由度，分别沿节点坐标的x，y，z方向，设四面体四个顶点分别表示ti(xi，yi，zi)，tj(xj，yj，zj)，tk(xk，yk，zk)，tl(xl，yl，zl)，通过采用四面体的四个顶点来表现位移形式的模式称为单元位移模式或位移函数，一般用一次多项式表示：由上式，得到4个单元位移场函数为：S_i为单元的插值函数或形函数，V为四面体单元的体积，而a_i，b_i，c_i，d_i，…d_l为：至于a_j，b_j，c_j，d_j，…和d_l各项按右手定则轮换下标i，j，k，l即可求得，例如：由上面式子可以得到：其中：是四面体每个顶点的位移向量，I为单位矩阵，它们组成的列向量称为四面体单元每个顶点的位移列阵，在三维空间中，Piola‑Kirchhoff应力σ被计算出，其表示为3×3的对称矩阵，如下所示：(32)树叶变形(321)模型表示将形变模型表示为(V，G，P)，其中V＝(1，2，…，n)表示顶点集，G＝{(i，j)|i，j∈V，i≠j}是边的集合，P＝{Pi∈R3|1≤i≤n}是模型中每个顶点的坐标集合；对于每个顶点i，用Ni＝{j|(i，j)∈G}表示i的领域的集合；同时，为了方便计算，设定形变模型中所有四面体的质量都被平均分配到四个顶点上，则mi表示i点的质量，有：其中，ρ_j为四面体密度，为常量；表示与i相邻的第j个四面体的体积；使用r向量表示顶点的初始位置，p向量表示顶点形变后的位置，则r_i和p_i分别表示顶点i的初始位置和形变后的位置；(322)定义形变模型中的势能按下式计算形变模型的弹性势能：V(P)＝Vr(P)+Vv(P)其中，Vr(P)表示使形变后的模型恢复到平衡或静止状态的势能，即恢复势能；Vv(P)表示保持模型形变所需要的能量，即体积势能；(323)恢复势能的定义按下式计算恢复势能：其中，Li表示i节点的微分算子，di表示i节点在平衡或静止状态时的微分坐标；Ri(P)是i节点的旋转矩阵，λ为形变模型的杨氏模量；对于树叶这种固体模型，采用拉普拉斯算子和拉普拉斯坐标作为微分算子和微分坐标，计算方法如下：其中，为边长的权值，l_k为边e_ij的对边e_pq的边长，θ_k为对边e_pq的二面角，是中心顶点V_i的归一化权值，V_dual是中心顶点V_i的对偶的Voronoi体域的体积，为标量场函数，Φ_i(·)是分段线性基函数，为中心顶点V_i处的基函数，v为四面体体积，S_i为相对的底三角形的面积，为底三角形的指向四面体外的法线方向；将表示成矩阵形式为：(δ^(x)，δ^(y)，δ^(z))＝L(P^(x)，P^(y)，P^(z))＝M^‑1L_s(P^(x)，P^(y)，P^(z))，其中，(P^(x)，P^(y)，P^(z))和(δ^(x)，δ^(y)，δ^(z))分别表示初始网格顶点坐标和拉普拉斯坐标的三个分量的向量表示形式，M是体积权值的对角线矩阵L_s是网格权值的对称稀疏线性矩阵：则拉普拉斯算子矩阵可以表示为：L＝M‑1Ls；根据网面拉普拉斯坐标性质可知，体网格的拉普拉斯坐标也具有平移不变性，不具有旋转不变性和缩放不变性；为了计算恢复势能，还要计算每个顶点i的旋转矩阵Ri(P)，计算方式如下：rij＝mj(rj‑ri)pij＝mj(pj‑pi)ri和pi分别表示节点i在变形前后的坐标，则节点i的线性变换矩阵Ai通过求解下式的极小值获得：当获得上式的极小值后，可得：其中，由于为对称矩阵，但此矩阵不包含任何的旋转形变，所以旋转矩阵应通过矩阵进行极分解获得；同时设定模型的旋转矩阵为T，则R_i(TP)＝TR_i(P)，即：L_iTP‑R_i(TP)d_i＝T(L_iP‑R_id_i)由前述可知Vr(TP)＝Vr(P)，即Vr(P)是旋转不变的；(4)模型形变模型形变模拟由下面的欧拉‑拉普拉斯运动方程控制：其中，m是系统的质量矩阵，D是系统的阻尼矩阵，V(P)是系统的势能，f_ext是系统所受的外力，是系统的内力，同时系统的内力的雅克比矩阵为即系统的刚度矩阵K，对于系统的阻尼矩阵D，采用Rayleigh阻尼：D＝αM+βK其中，K∈R3n，3n是系统的刚度矩阵，M是系统的质量矩阵，α，β是常数；求解过程中设定β＝0，采用隐式的欧拉方法或者采用隐式的Newmark方法求解上式；最后将体模型剖分成四面体模型，并对该模型施加外力，通过求解方程获取受力变形后的叶片。