CN105654543B - 面向激光点云数据的阔叶树真实叶片建模与形变方法 - Google Patents

Abstract

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

Description

面向激光点云数据的阔叶树真实叶片建模与形变方法

技术领域

本发明涉及一种面向激光点云数据的阔叶树真实叶片建模与形变方法，尤其涉及一种基于LTS(Transport Layer Security，激光扫描)点云数据的阔叶树叶面建模和形变方法。

背景技术

阔叶树一般指双子叶植物类的树木，具有扁平、较宽阔叶片，叶脉成网状，叶常绿或落叶，一般叶面宽阔，叶形随树种不同而有多种形状的多年生木本植物。有的常绿、落叶类大多在秋冬季节叶从枝上脱落。阔叶树的经济价值很大，不少为重要用材树种，而且其中有些为名贵木材、景观植物、各种水果等，还有一些阔叶树用作行道树或庭园绿化树种。本发明研究是阔叶树的景观植物含笑和樱花树作为叶面重建和形变的研究对象。

数据获取借助地面激光扫描仪，地面激光扫描仪不会对被测物造成任何损伤，且能以点云的形式精确还原出目标体的三维数据，而且三维激光扫描仪在计测学中具有无可比拟的优势，因此国外许多林业科研工作者就地面三维激光扫描技术在林业中的应用进行了深入研究和探讨。但目前为止尚未发现利用离散点云叶片数据进行建模，这主要是因为树木外形特征无规律且形态复杂，并且外界环境对树木的状态产生着持续的影响，获取没有外界环境影响的点云数据技术很难；因此考虑使用计算机自动对阔叶树叶片进行建模和变形外，需要解决树木受到外部环境如风吹抖动及遮挡的影响；树木枝繁叶茂，叶子的形态及方位角度不固定；如何去除偏差并得到真实的形变的叶面数据是需要考虑的问题；上述外在因素都是使用计算机研究林木的阻力，因此如何从离散的激光点云中自动获取精确林学指标是亟待解决的问题。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种面向激光点云数据的阔叶树真实叶片建模与形变方法，旨在借助激光扫描仪技术，搭建一个精确可行的三维活立木数据采集分析平台，融合图形图像学的最新方法，通过计算机自动分析获取精确林学指标，从而准确描述不同林分条件下的活立木动态生长变化的叶面积指数。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种面向激光点云数据的阔叶树真实叶片建模与形变方法，首先运用高精度的激光扫描仪获取不同阔叶树种的三维叶面点云；其次从扫描获取的叶子点云中采用多项式拟合的方法获取精确的树叶边界，并结合计算机图形学的方法计算叶面的主叶脉；接着采用双三次广义张量积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条扫描线，其中i＝1，2，3...，n+1，定位求取L_2，i左右两端的端点p_2，li和p_2，ri，即可以获取叶面在不同部分的边缘点；

(12)对步骤(11)的多条L_2，i扫描线，取每条扫描线的左右两端的端点p_2，li和p_2，ri，分别记为p_2，li：(x_li，y_li，z_li)和p_2，ri：(x_ri，y_ri，z_ri)，其中i＝1，2，3...，n；除去两侧的两条扫描线，若某条扫描线的左端点相对于相邻两条扫描线的左端点向主叶脉内凹，且右端点相对于相邻两条扫描线的右端点也向主叶脉内凹，则认为该扫描线存在遮挡或由误差造成，剔除该扫描线；

(13)基于步骤(12)，剔除完所有存在遮挡或由误差造成的扫描线后，对所有扫描线按序重新编号为p_2，lj和p_2，rj，左右两端的端点分别记为p_2，lj：(x_lj，y_lj，z_lj)和p_2，rj：(x_rj，y_rj，z_rj)；最终得到：

左端点的集合为p_2l＝{(x_l1，y_l1，z_l1)，(x_l2，y_l2，z_l2)，...，(x_lj，y_lj，z_lj)，...}；

右端点的集合为p_2r＝{(x_r1，y_r1，z_r1)，(x_r2，y_r2，z_r2)，...，(x_rj，y_rj，z_rj)，...}

左端点的集合和右端点的集合构成了叶面的初始扫描轮廓；对于一片叶子的左/右边部分轮廓，采用不同投影面拟合与求交的方法定位真实叶子的边缘；具体步骤如下：

(131)对于左/右边部分轮廓上的边缘点p_2，j＝(x，y，z)，将(x_i，y_i，z_i)i＝1，2，3...，n作为叶子几何边缘的参数，运用多项式曲线拟合的方法，把边缘点y值作为输入参数，计算拟合系数来求取相应的x′和z′，公式为：

整体算式表示为：

x≈x′＝v_x(y)＝v_x1yⁿ+v_x2y^n-1+v_x3y^n-2+...+v_xn′-1y+v_xn′

z≈z′＝v_z(y)＝v_z1yⁿ+v_z2y^n-1+v_z3y^n-2+...+v_zn′-1y+v_zn′

其中，v_x1 v_x2 v_x3...v_xn′-1 v_xn′为计算得到的多项式拟合系数；

(132)得到新的叶面边缘点为P′_edge＝{x′_l，y_l，z′_l；x′_r，y_r，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个点，分别为p_0，j，p_1，j，p_2，j，p_3，j，再选择相邻的四条线段，共取16个点构成网格的16个点，把这16个点代入广义张量积Bezier曲面公式p_i，j中，进而采用广义双三次Bezier曲面对其拟合，最后，重复上述过程直到所有扫描线上的扫描点都遍历一遍；

(3)叶面受力形变

根据权利要求2所述，从扫描获取的叶子点云中采用多项式拟合的方法获取精确的树叶边界，在精确的树叶里进行叶面受力形变。

(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′|²-|d|²＝d′^Td′-d^Td＝d^T(F^TF-I)d

根据弹性力学知识，可由上式得出格林应变E∈R^3×3为：

从上式可看出E是对称矩阵；

对于以四面体为基本单元组成的叶面模型而言，采用线性插值函数来表示每个节点的位移向量u＝(u，v，w)；对于一个由4个节点组成的四面体，每个节点有3个自由度，分别沿节点坐标的x，y，z方向，设四面体四个顶点分别表示t_i(x_i，y_i，z_i)，t_j(x_j，y_j，z_j)，t_k(x_k，y_k，z_k)，t_l(x_l，y_l，z_l)，通过采用四面体的四个顶点来表现位移形式的模式称为单元位移模式或位移函数，一般用一次多项式表示：

由上式，得到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＝{P_i∈R³|1≤i≤n}是模型中每个顶点的坐标集合；对于每个顶点i，用N_i＝{j|(i，j)∈G}表示i的领域的集合；同时，为了方便计算，设定形变模型中所有四面体的质量都被平均分配到四个顶点上，则m_i表示i点的质量，有：

其中，ρ_j为四面体密度，为常量；表示与i相邻的第j个四面体的体积；使用r向量表示顶点的初始位置，p向量表示顶点形变后的位置，则r_i和p_i分别表示顶点i的初始位置和形变后的位置；

(322)定义形变模型中的势能

按下式计算形变模型的弹性势能：

V(P)＝V_r(P)+V_v(P)

其中，V_r(P)表示使形变后的模型恢复到平衡或静止状态的势能，即恢复势能；V_v(P)表示保持模型形变所需要的能量，即体积势能；

(323)恢复势能的定义

按下式计算恢复势能：

其中，L_i表示i节点的微分算子，d_i表示i节点在平衡或静止状态时的微分坐标；R_i(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^-1L_s；

根据网面拉普拉斯坐标性质可知，体网格的拉普拉斯坐标也具有平移不变性，不具有旋转不变性和缩放不变性；

为了计算恢复势能，还要计算每个顶点i的旋转矩阵R_i(P)，计算方式如下：

r_ij＝m_j(r_j-r_i)

p_ij＝m_j(p_j-p_i)

r_i和p_i分别表示节点i在变形前后的坐标，则节点i的线性变换矩阵A_i通过求解下式的极小值获得：

当获得上式的极小值后，可得：

其中，

由于为对称矩阵，但此矩阵不包含任何的旋转形变，所以旋转矩阵应通过矩阵进行极分解获得；同时设定模型的旋转矩阵为T，则R_i(TP)＝TR_i(P)，即：

L_iTP-R_i(TP)d_i＝T(L_iP-R_id_i)

由前述可知V_r(TP)＝V_r(P)，即V_r(P)是旋转不变的；

(4)模型形变

根据权利要求3所述，固体力学受力模型，定义叶脉与叶肉不同的力学属性，并构造非线性有限元的受力形变方程，模拟树叶形变。

模型形变模拟由下面的欧拉-拉普拉斯运动方程控制：

其中，M是系统的质量矩阵，D是系统的阻尼矩阵，V(P)是系统的势能，f_ext是系统所受的外力，是系统的内力，同时系统的内力的雅克比矩阵为即系统的刚度矩阵K，对于系统的阻尼矩阵D，采用Rayleigh阻尼：

D＝αM+βK

其中，K∈R^3n，3n是系统的刚度矩阵，M是系统的质量矩阵，α，β是常数；求解过程中设定β＝0，采用隐式的欧拉方法或者采用隐式的Newmark方法求解上式；

最后将体模型剖分成四面体模型，并对该模型施加外力，通过求解方程获取受力变形后的叶片。

具体实施过程，根据叶面积测量仪和各种算法得到叶面积指数比对，得到本发明实验数据有一定的可行性。

有益效果：本发明提供的面向激光点云数据的阔叶树真实叶片建模与形变方法，借助激光扫描仪技术，搭建了一个精确可行的三维活立木数据采集分析平台，融合图形图像学的最新方法，通过计算机自动分析获取精确林学指标，能够准确描述不同林分条件下的活立木动态生长变化的叶面积指数。

附图说明

图1根据边缘点进行扫描线的选取；

图2广义Bezier曲面算法实施图；

图3叶面的边缘定位与重建，由步骤(3a)向步骤(3f)依次进行；

图4为叶边缘确定的过程图，由步骤(4a)向步骤(4b)依次进行；

图5含笑树叶面三角体网格化与受力形变图，由步骤(5a)向步骤(5f)依次进行；

图6樱花树叶面三角体网格化与受力形变图，由步骤(6a)向步骤(6f)依次进行；

具体实施方式

下面结合附图对本发明作更进一步的说明。

一种面向激光点云数据的阔叶树真实叶片建模与形变方法，首先运用高精度的激光扫描仪获取不同阔叶树种的三维叶面点云；其次从扫描获取的叶子点云中采用多项式拟合的方法获取精确的树叶边界，并结合计算机图形学的方法计算叶面的主叶脉；接着采用双三次广义张量积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条扫描线，其中i＝1，2，3...，n+1，定位求取L_2，i左右两端的端点p_2，li和p_2，ri，即可以获取叶面在不同部分的边缘点；

(12)对步骤(11)的多条L_2，i扫描线，取每条扫描线的左右两端的端点p_2，li和p_2，ri，分别记为p_2，li：(x_li，y_li，z_li)和p_2，ri：(x_ri，y_ri，z_ri)，其中i＝1，2，3...，n；除去两侧的两条扫描线，若某条扫描线的左端点相对于相邻两条扫描线的左端点向主叶脉内凹，且右端点相对于相邻两条扫描线的右端点也向主叶脉内凹，则认为该扫描线存在遮挡或由误差造成，剔除该扫描线；

(13)基于步骤(12)，剔除完所有存在遮挡或由误差造成的扫描线后，对所有扫描线按序重新编号为p_2，lj和p_2，rj，左右两端的端点分别记为p_2，li：(x_lj，y_lj，z_lj)和p_2，rj：(x_rj，y_rj，z_rj)；最终得到：

左端点的集合为p_2l＝{(x_l1，y_l1，z_l1)，(x_l2，y_l2，z_l2)，...，(x_lj，y_lj，z_lj)，...}；

右端点的集合为p_2r＝{(x_r1，y_r1，z_r1)，(x_r2，y_r2，z_r2)，...，(x_rj，y_rj，z_rj)，...}

左端点的集合和右端点的集合构成了叶面的初始扫描轮廓；对于一片叶子的左/右边部分轮廓，采用不同投影面拟合与求交的方法定位真实叶子的边缘；具体步骤如下：

(131)对于左/右边部分轮廓上的边缘点p_2，j＝(x，y，z)，将(x_i，y_i，z_i)i＝1，2，3...，n作为叶子几何边缘的参数，运用多项式曲线拟合的方法，把边缘点y值作为输入参数，计算拟合系数来求取相应的x′和z′，公式为：

整体算式表示为：

x≈x′＝v_x(y)＝v_x1yⁿ+v_x2y^n-1+v_x3y^n-2+...+v_xn′-1y+v_xn′

z≈z′＝v_z(y)＝v_z1yⁿ+v_z2y^n-1+v_z3y^n-2+...+v_zn′-1y+v_zn′

其中，v_x1 v_x2 v_x3...v_xn′-1 v_xn′为计算得到的多项式拟合系数；

(132)得到新的叶面边缘点为P′_edge＝{x′_l，y_l，z′_l；x′_r，y_r，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个点，分别为p_0，j，p_1，j，p_2，j，p_3，j，再选择相邻的四条线段，共取16个点构成网格的16个点，把这16个点代入广义张量积Bezier曲面公式p_i，j中，进而采用广义双三次Bezier曲面对其拟合，最后，重复上述过程直到所有扫描线上的扫描点都遍历一遍；

(3)叶面受力形变

根据权利要求2所述，从扫描获取的叶子点云中采用多项式拟合的方法获取精确的树叶边界，在精确的树叶里进行叶面受力形变。

(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′|²-|d|²＝d′^Td′-d^Td＝d^T(F^TF-I)d

根据弹性力学知识，可由上式得出格林应变E∈R^3×3为：

从上式可看出E是对称矩阵；

对于以四面体为基本单元组成的叶面模型而言，采用线性插值函数来表示每个节点的位移向量u＝(u，v，w)；对于一个由4个节点组成的四面体，每个节点有3个自由度，分别沿节点坐标的x，y，z方向，设四面体四个顶点分别表示t_i(x_i，y_i，z_i)，t_j(x_j，y_j，z_j)，t_k(x_k，y_k，z_k)，t_l(x_l，y_l，z_l)，通过采用四面体的四个顶点来表现位移形式的模式称为单元位移模式或位移函数，一般用一次多项式表示：

由上式，得到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＝{P_i∈R³|1≤i≤n}是模型中每个顶点的坐标集合；对于每个顶点i，用N_i＝{j|(i，j)∈G}表示i的领域的集合；同时，为了方便计算，设定形变模型中所有四面体的质量都被平均分配到四个顶点上，则m_i表示i点的质量，有：

其中，ρ_j为四面体密度，为常量；表示与i相邻的第j个四面体的体积；使用r向量表示顶点的初始位置，p向量表示顶点形变后的位置，则r_i和p_i分别表示顶点i的初始位置和形变后的位置；

(322)定义形变模型中的势能

按下式计算形变模型的弹性势能：

V(P)＝V_r(P)+V_v(P)

其中，V_r(P)表示使形变后的模型恢复到平衡或静止状态的势能，即恢复势能；V_v(P)表示保持模型形变所需要的能量，即体积势能；

(323)恢复势能的定义

按下式计算恢复势能：

其中，L_i表示i节点的微分算子，d_i表示i节点在平衡或静止状态时的微分坐标；R_i(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^-1L_s；

根据网面拉普拉斯坐标性质可知，体网格的拉普拉斯坐标也具有平移不变性，不具有旋转不变性和缩放不变性；

为了计算恢复势能，还要计算每个顶点i的旋转矩阵R_i(P)，计算方式如下：

r_ij＝m_j(r_j-r_i)

p_ij＝m_j(p_j-p_i)

r_i和p_i分别表示节点i在变形前后的坐标，则节点i的线性变换矩阵A_i通过求解下式的极小值获得：

当获得上式的极小值后，可得：

其中，

由于为对称矩阵，但此矩阵不包含任何的旋转形变，所以旋转矩阵应通过矩阵进行极分解获得；同时设定模型的旋转矩阵为T，则R_i(TP)＝TR_i(P)，即：

L_iTP-R_i(TP)d_i＝T(L_iP-R_id_i)

由前述可知V_r(TP)＝V_r(P)，即V_r(P)是旋转不变的；

(4)模型形变

根据权利要求3所述，固体力学受力模型，定义叶脉与叶肉不同的力学属性，并构造非线性有限元的受力形变方程，模拟树叶形变。

模型形变模拟由下面的欧拉-拉普拉斯运动方程控制：

其中，M是系统的质量矩阵，D是系统的阻尼矩阵，V(P)是系统的势能，f_ext是系统所受的外力，是系统的内力，同时系统的内力的雅克比矩阵为即系统的刚度矩阵K，对于系统的阻尼矩阵D，采用Rayleigh阻尼：

D＝αM+βK

其中，K∈R^3n，3n是系统的刚度矩阵，M是系统的质量矩阵，α，β是常数；求解过程中设定β＝0，采用隐式的欧拉方法或者采用隐式的Newmark方法求解上式；

最后将体模型剖分成四面体模型，并对该模型施加外力，通过求解方程获取受力变形后的叶片。

具体实施过程，根据叶面积测量仪和各种算法得到叶面积指数比对，得到本发明实验数据有一定的可行性。

表.1 AM-300测量出的叶面积与重建后的叶面积

表.2点云预处理中点云数目和重建各个阶段所需时间

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

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)的多条L_2，i扫描线，取每条扫描线的左右两端的端点p_2，li和p_2，ri，分别记为p_2，li：(x_li，y_li，z_li)和p_2，ri：(x_ri，y_ri，z_ri)，其中i＝1，2，3...，n+1；除去两侧的两条扫描线，若某条扫描线的左端点相对于相邻两条扫描线的左端点向主叶脉内凹，且右端点相对于相邻两条扫描线的右端点也向主叶脉内凹，则认为该扫描线存在遮挡或由误差造成，剔除该扫描线；

(13)基于步骤(12)，剔除完所有存在遮挡或由误差造成的扫描线后，对所有扫描线按序重新编号为p_2，lj和p_2，rj，左右两端的端点分别记为p_2，lj：(x_lj，y_lj，z_lj)和p_2，rj：(x_rj，y_rj，z_rj)；最终得到：

左端点的集合为p_2l＝{(x_l1，y_l1，z_l1)，(x_l2，y_l2，z_l2)，...，(x_lj，y_lj，z_lj)，...}；

右端点的集合为p_2r＝{(x_r1，y_r1，z_r1)，(x_r2，y_r2，z_r2)，...，(x_rj，y_rj，z_rj)，...}左端点的集合和右端点的集合构成了叶面的初始扫描轮廓；对于一片叶子的左/右边部分轮廓，采用不同投影面拟合与求交的方法定位真实叶子的边缘；具体步骤如下：

(131)对于左/右边部分轮廓上的边缘点p_2，j＝(x，y，z)，将(x_i，y_i，z_i)i＝1，2，3...，n作为叶子几何边缘的参数，运用多项式曲线拟合的方法，把边缘点y值作为输入参数，计算拟合系数来求取相应的x′和z′，公式为：

整体算式表示为：

x≈x′＝v_x(y)＝v_x1yⁿ+v_x2y^n-1+v_x3y^n-2+...+v_xn′-1y+v_xn′

z≈z′＝v_z(y)＝v_z1yⁿ+v_z2y^n-1+v_z3y^n-2+...+v_zn′-1y+v_zn′

其中，v_x1v_x2v_x3...v_xn′-1v_xn′为计算得到的多项式拟合系数；

(132)得到新的叶面边缘点为P′_edge＝{x′_l，y_l，z′_l；x′_r，y_r，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个点，分别为p_0，j，p_1，j，p_2，j，p_3，j，再选择相邻的四条线段，共取16个点构成网格的16个点，把这16个点代入广义张量积Bezier曲面公式p_i，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′|²-|d|²＝d′^Td′-d^Td＝d^T(F^TF-I)d

根据弹性力学知识，可由上式得出格林应变E∈R^3×3为：

从上式可看出E是对称矩阵，对于以四面体为基本单元组成的叶面模型而言，采用线性插值函数来表示每个节点的位移向量u＝(u，v，w)；对于一个由4个节点组成的四面体，每个节点有3个自由度，分别沿节点坐标的x，y，z方向，设四面体四个顶点分别表示t_i(x_i，y_i，z_i)，t_j(x_j，y_j，z_j)，t_k(x_k，y_k，z_k)，t_l(x_l，y_l，z_l)，通过采用四面体的四个顶点来表现位移形式的模式称为单元位移模式或位移函数，一般用一次多项式表示：

由上式，得到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＝{P_i∈R³|1≤i≤n}是模型中每个顶点的坐标集合；对于每个顶点i，用N_i＝{j|(i，j)∈G}表示i的领域的集合；同时，为了方便计算，设定形变模型中所有四面体的质量都被平均分配到四个顶点上，则m_i表示i点的质量，有：

其中，ρ_j为四面体密度，为常量；表示与i相邻的第j个四面体的体积；使用r向量表示顶点的初始位置，p向量表示顶点形变后的位置，则r_i和p_i分别表示顶点i的初始位置和形变后的位置；

(322)定义形变模型中的势能

按下式计算形变模型的弹性势能：

V(P)＝V_r(P)+V_v(P)

其中，V_r(P)表示使形变后的模型恢复到平衡或静止状态的势能，即恢复势能；V_v(P)表示保持模型形变所需要的能量，即体积势能；

(323)恢复势能的定义

按下式计算恢复势能：

其中，L_i表示i节点的微分算子，d_i表示i节点在平衡或静止状态时的微分坐标；

R_i(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^-1L_s；

根据网面拉普拉斯坐标性质可知，体网格的拉普拉斯坐标也具有平移不变性，不具有旋转不变性和缩放不变性；

为了计算恢复势能，还要计算每个顶点i的旋转矩阵R_i(P)，计算方式如下：

r_ij＝m_j(r_j-r_i)

p_ij＝m_j(p_j-p_i)

r_i和p_i分别表示节点i在变形前后的坐标，则节点i的线性变换矩阵A_i通过求解下式的极小值获得：

当获得上式的极小值后，可得：

其中，

由于为对称矩阵，但此矩阵不包含任何的旋转形变，所以旋转矩阵应通过矩阵进行极分解获得；同时设定模型的旋转矩阵为T，则R_i(TP)＝TR_i(P)，即：L_iTP-R_i(TP)d_i＝T(L_iP-R_id_i)

由前述可知V_r(TP)＝V_r(P)，即V_r(P)是旋转不变的；

(4)模型形变

模型形变模拟由下面的欧拉-拉普拉斯运动方程控制：

其中，m是系统的质量矩阵，D是系统的阻尼矩阵，V(P)是系统的势能，f_ext是系统所受的外力，是系统的内力，同时系统的内力的雅克比矩阵为即系统的刚度矩阵K，对于系统的阻尼矩阵D，采用Rayleigh阻尼：

D＝αM+βK

其中，K∈R^3n，3n是系统的刚度矩阵，M是系统的质量矩阵，α，β是常数；求解过程中设定β＝0，采用隐式的欧拉方法或者采用隐式的Newmark方法求解上式；最后将体模型剖分成四面体模型，并对该模型施加外力，通过求解方程获取受力变形后的叶片。