一种二阶几何连续的全凸闭合曲线的构建方法

技术领域

本发明涉及一种二阶几何连续的全凸闭合曲线的构建方法。

背景技术

地面三维激光扫描技术(Terrestrial Laser Scanning Technology，简称TLS)是一种新型的测绘技术，产生于20世纪90年代。TLS可快速准确的获得目标物体表面的点云数据信息，从而可以快速构建目标物体的三维模型。经过近20多年来的发展，三维激光扫描仪已经可以连续快速地对被观测物体进行非接触式测量，其通过获取物体表面至扫描仪的距离和发射强度获取大量物体表面的三维点云数据。使用TLS提取林业相关参数是近10年来的研究热点。

曲线插值是应用数学、计算机图形学及计算机辅助设计与制造领域中的经典问题。当给定的数据为凸的时候，构建的曲线也是凸的，此类曲线插值问题称为曲线保凸插值。现有的曲线保凸插值方法注重于构建一段非闭合的凸曲线，多通过在子区间增加节点、采用张力参数、曲线细分的方式实现。对于给定的凸包点，是否存在过凸包点的光滑的全凸闭合曲线，以及是否存在唯一的过凸包点的光滑的全凸闭合曲线，均有待于进一步研究。

发明内容

有鉴于此，本发明的目的是针对现有技术的不足，提供一种二阶几何连续的全凸闭合曲线的构建方法，用以通过在相邻凸包点间构建一条三次Bezier曲线，在曲线连接点处满足二阶几何连续性以构建连续的全凸闭合曲线。

为达到上述目的，本发明采用以下技术方案：

一种二阶几何连续的全凸闭合曲线的构建方法，包括如下步骤：

构建凸包点集P＝{p_i|i＝0,1,...,n-1}，且有当i≤0,或i>n时，i表示为i对n取模，如且p₀,p₁,...,p_n-1是按照逆时针次序排列的点集；

采用分段三次Bezier曲线构建全凸闭合曲线；

构建并求解方程组；

解的优化及全凸闭合曲线的计算，得到由分段三次Bezier曲线构成的二阶几何连续的全凸闭合曲线。

优选的，所述采用分段三次Bezier曲线构建全凸闭合曲线，包括：

在相邻凸包点p_i与p_i+1间构建一条三次Bezier曲线，以凸包点p_i与p_i+1作为所述三次Bezier曲线的第1个和第4个控制点，假设点q_2i与q_2i+1是所述三次Bezier曲线的第2个与第3个控制点；

点p_i与点q_2i构成的向量与与点p_i-1与点p_i+1构成的向量平行，V_i是向量的单位向量；点q_2i+1与点p_i+1构成的向量与与点p_i与点p_i+2构成的向量平行，V_i+1是向量的单位向量；点q_2i与点q_2i+1构成的向量与点p_i与点p_i+1构成的向量平行，L_i是向量的单位向量；

根据构建的几何关系，点q_2i+1与q_2i+2可表示为：

q_2i+1＝p_i+1-k_i+1V_i+1,k_i+1>0

q_2i+2＝f_i+1V_i+1+p_i+1,f_i+1>0；

依次在所有相邻的凸包点之间构建一条三次Bezier曲线。

优选的，所述构建并求解方程组，包括：

根据点p_i、q_2i、q_2i+1与点p_i+1之间的闭合关系、三次Bezier曲线的表达式、三次Bezier曲线的一阶二阶导数的函数表达式及分段曲线在连接点二阶几何连续的条件可得到如下方程组：

其中f_i、g_i与k_i均为大于0的实数；

方程组的解可简要描述如下：

任给μ>0，选择凸包点集P中第k(0≤k<n)个节点为初始节点，令实数α_i与β_i分别表示如下：

给凸包点集P中每一个点p_i设置一个大于0的参数值u_i，初始节点p_k的参数值为u_k，

令

令e_i是方程的正根，则有：

则方程组的解为

优选的，所述解的优化及全凸闭合曲线的计算，得到由分段三次Bezier曲线构成的二阶几何连续的全凸闭合曲线，包括：

选择比值β_k/α_k最小的节点μ_k作为初始节点，并设定μ_k＝μ；

(1)根据公式计算每一个凸包点对应的α_i、β_i与的值，并遍历得到的最小值

(2)对于给定的μ参数值，根据公式μ_k＝μ,计算每个凸包点对应的参数值μ_i；

(3)根据公式计算每个凸包点对应的参数值e_i；

(4)根据公式f_i＝||l_i-1+l_i||e_i与公式计算每个凸包点对应的参数值f_i与k_i；

(5)根据参数值f_i、k_i、公式q_2i+1＝p_i+1-k_i+1V_i+1与公式q_2i+2＝f_i+1V_i+1+p_i+1计算分别得到每段三次Bezier曲线的第2个与第3个控制点q_2i+1与q_2i+2；

(6)依次分别以点p_i、q_2i+1、q_2i+2与p_i+1为控制点构建一条三次Bezier曲线得到由分段三次Bezier曲线构成的二阶几何连续的全凸闭合曲线。

优选的，还包括：获取凸包点。

优选的，所述获取凸包点，包括：

从地面三维激光扫描仪获取的树木点云中提取树干点云；

计算树干在高度h处树干横断面的定位点；

计算树干点云中的最低点，并以所述最低点为基础计算树干的高度；

以最低点为基准，以向量为平面的法向量，在最低点的上方h处与h+预设距离处分别构建的两个平面；位于两个平面中间的树干点云构成了一个树干点云块，将树干点云块点云投影到树干高度h处的平面上得到一个平面点云集，计算平面点云集构成的凸包多边形的质心点，所述质心点为平面点云集的几何中心点与树干在高度h处树干横断面的定位点；

计算树干在高度h处的生长方向与树干横断面；

树干在高度h处的生长方向即是树干在高度h处的树干横断面的法向量；

采用迭代方法计算树干在高度h处的生长方向；第k次迭代计算为：构建过定位点且法向量为的一个平面，所述平面为基准平面；

计算树干在高度h处的生长方向其中，协方差矩阵的最大特征值对应的特征向量为树干在高度h处的生长方向平面法向量赋值为进行第k+1次迭代；初次迭代时法向量为迭代结束条件为：相邻两次迭代得到的两个树干生长方向与的夹角θ_k小于0.5度或者θ_k与θ_k+1的角度差小于0.5度；迭代结束时的树干生长方向即为树干在高度h处的生长方向；以树干在高度h处的生长方向为平面的法向量，通过定位点构建树干在高度h处的树干横断面；

获取用于计算树干直径的凸包点。

优选的，所述从地面三维激光扫描仪获取的树木点云中提取树干点云，包括：

使用地面三维激光扫描仪配备的点云处理软件对地面三维激光扫描仪在多站扫描的点云进行配准得到单木的点云数据，经移除地表点云、移除树枝点云、移除离散噪声点后得到树干点云。

优选的，遍历树干点云集得到z坐标值最小的一个点，并将得到的所述点设为树干点云的最低点。

优选的，所述计算树干在高度h处的生长方向包括：

以基准平面为基础，分别构建5个与此平行的平面，其中3个平面位于基准平面上方，2个平面位于基准平面下方，且相邻平面间的距离为0.5厘米；

得到5个树干点云块与5个几何中心点；

对5个几何中心点采用主成分分析方法计算树干在高度h处的生长方向

优选的，所述获取用于计算树干直径的凸包点，包括：

在树干横断面上方构建一个与树干横断面平行且距离为1厘米的平面，位于树干横断面与平面间的树干点云即为用于计算树干直径的树干点云，将此部分点云投影至树干横断面上得到一个三维空间上的平面点云集，采用几何变换中的旋转操作，将此平面点云集旋转为与水平面平行的一个二维空间上的平面点云集，根据凸包算法计算此平面点云集的凸包点集P。

本发明的有益效果是：

本发明提供一种二阶几何连续的全凸闭合曲线的构建方法，通过在相邻凸包点间构建一条三次Bezier曲线，在曲线连接点处满足二阶几何连续性以构建连续的全凸闭合曲线，该方法计算量小，无需求解方法组，简单实用，而且构建的曲线有物理学意义上的解释。

采用本发明的方法应用在在林业调查中，从树干点云中提取直径的准确性较高。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

图1本发明实施例1提供的方法流程图；

图2本发明构建的分段三次Bezier曲线构建的几何关系图；

图3相同凸包点不同μ值下构建的全凸闭合曲线的示例图；

图4模拟的围尺轨迹示例图；

图5本发明实施例2提供的方法流程图；

图6是定位树干横断面的示意图；

图7是获取用于计算树干直径的点云示例图。

具体实施方式

下面结合附图和实施例对本发明作进一步描述。

实施例1：

如图1所示，本实施例提供的一种二阶几何连续的全凸闭合曲线的构建方法，包括如下步骤：

步骤S101，对于给定的凸包点集P＝{p_i|i＝0,1,...,n-1}，且有当i≤0,或i>n时，i表示为i对n取模，如且p₀,p₁,...,p_n-1是按照逆时针次序排列的点集；

步骤S102，构建全凸闭合曲线的几何关系图，采用分段三次Bezier曲线构建全凸闭合曲线；

在相邻凸包点p_i与p_i+1间构建一条三次Bezier曲线，以凸包点p_i与p_i+1作为此条三次Bezier曲线的第1个和第4个控制点，假设点q_2i与q_2i+1是此条三次Bezier曲线的第2个与第3个控制点，其中，点p_i与点q_2i构成的向量与与点p_i-1与点p_i+1构成的向量平行(其中V_i是向量的单位向量)，点q_2i+1与点p_i+1构成的向量与与点p_i与点p_i+2构成的向量平行(其中V_i+1是向量的单位向量)，点q_2i与点q_2i+1构成的向量与点p_i与点p_i+1构成的向量平行(其中L_i是向量的单位向量)，由此，点q_2i+1与q_2i+2可表示如下：

q_2i+1＝p_i+1-k_i+1V_i+1,k_i+1>0

q_2i+2＝f_i+1V_i+1+p_i+1,f_i+1>0

依次在所有相邻的凸包点之间构建一条三次Bezier曲线，曲线的构建示意图如图2所示，点p_i、q_2i、q_2i+1与p_i+1依次是一条三次Bezier曲线的四个控制点，其中点p_i与p_i+1是相邻的两个凸包点，点q_2i与q_2i+1是需要计算得到的控制点；点p_i与点q_2i构成的向量与与点p_i-1与点p_i+1构成的向量平行，V_i是向量的单位向量；点q_2i+1与点p_i+1构成的向量与与点p_i与点p_i+2构成的向量平行，V_i+1是向量的单位向量；点q_2i与点q_2i+1构成的向量与点p_i与点p_i+1构成的向量平行，L_i是向量的单位向量；

步骤S103，构建并求解方程组；

根据点p_i、q_2i、q_2i+1与点p_i+1之间的闭合关系、三次Bezier曲线的表达式、三次Bezier曲线的一阶二阶导数的函数表达式及分段曲线在连接点二阶几何连续的条件可得到如下方程组：

其中f_i,g_i,k_i>0；

方程组的解可简要描述如下：

任给μ>0，选择凸包点集P中第k(0≤k<n)个节点为初始节点，令：

给凸包点集P中每一个点p_i设置一个大于0的参数值u_i，初始节点p_k的参数值为u_k，

令

令e_i是方程的正根，则有：

则方程组的解为

步骤S104，解的优化及全凸闭合曲线的计算，得到由分段三次Bezier曲线构成的二阶几何连续的全凸闭合曲线。

根据上述方程组解的描述可以看出，在μ值确定的情况下，选择不同的节点p_k(0≤k<n)作为初始节点，参数值u_i的值不同，方程组的解不同，即在μ值确定的情况下，曲线的形状与初始节点的选择相关。为避免曲线形状与初始节点选择的相关性，本发明选择比值β_k/α_k最小的节点μ_k作为初始节点，并设定μ_k＝μ，使得曲线的形状只与参数μ值相关。实验表明，对于给定的凸包点，μ值越大，相邻凸包点间的曲线长度越小，从而闭合曲线的总长度也越小，这相当于曲线绷的越紧，如图3所示，图(a)、(b)、(c)分别表示相同凸包点不同μ值下构建的全凸闭合曲线的示例图，图(a)中，μ＝1；图(b)中，μ＝100；图(c)中，μ＝10000。

全凸闭合曲线的计算方法为：

(1)根据公式计算每一个凸包点点对应的α_i，β_i与的值，并遍历得到的最小值

(2)对于给定的μ参数值，根据公式μ_k＝μ,计算每个凸包点对应的参数值μ_i；

(3)根据公式计算每个凸包点对应的参数值e_i；

(4)根据公式f_i＝||l_i-1+l_i||e_i与公式计算每个凸包点对应的参数值f_i与k_i；

(5)根据参数值f_i、k_i、公式q_2i+1＝p_i+1-k_i+1V_i+1与公式q_2i+2＝f_i+1V_i+1+p_i+1计算分别得到每段三次Bezier曲线的第2个与第3个控制点q_2i+1与q_2i+2；

(6)依次分别以点p_i、q_2i+1、q_2i+2与p_i+1为控制点构建一条三次Bezier曲线得到由分段三次Bezier曲线构成的二阶几何连续的全凸闭合曲线。

解的正确性：

对于给定的4个控制点p₀、p₁、p₂与p₃，三次Bezier曲线C(u)可定义如下：

C(u)＝(1-u)³p₀+3u(1-u)²p₁+3u²(1-u)p₂+u³p₃,0u>

三次Bezier曲线C(u)在参数u处的一阶参数导数与二阶参数导数分别为：

C′(u)＝3(p₁-p₀)(1-u)²+6(p₂-p₁)(1-u)u+3(p₃-p₂)u²

C(u)＝6(p₂-2p₁+p₀)(1-u)+6(p₃-2p₂+p₁)u

根据本发明中的构建关系，对于由p_i、q_2i+1、q_2i+2与p_i+1为控制点构建一条三次Bezier曲线C_i，曲线在节点p_i与p_i+1处的一阶与二阶导数可分别表示为：

C_i+1′(0)＝3(q_2i+2-p_i+1)，C_i′(1)＝3(p_i+1-q_2i+1)

C_i+1(0)＝6(q_2i+3-2q_2i+2+p_i+1)，C_i(1)＝6(p_i+1-2q_2i+1+q_2i)

证明：

闭合性：曲线C_n-1过点p_n，曲线C₀过点p₀，因，从而曲线C_n-1与C₀在点p₀点连接。这也证明了两段曲线C_i与C_i+1的0阶几何连续性。

一阶几何连续性：根据相邻三次Bezier曲线在连接点处的一阶导数有：

即本发明构造的分段Bezier曲线在连接点处已满足一阶几何连续性。

二阶几何连续性：

根据相邻三次Bezier曲线在连接点处的二阶导数有：

C_i+1(0)＝6(g_i+1L_i+1-f_i+1V_i+1)，C_i(1)＝6(k_i+1V_i+1-g_i+1L_i+1)

根据二阶几何连续性的条件，在满足1阶几何连续的基础上，若要求两条曲线在连接点处达到二阶几何连续，则要求两条曲线在结合处点有共同的曲率向量：

上式可化简为：

列式计算，上式可改写为：

即有

因根据行列式初等变换有：

等价于存在λ_i+1，使得：；

根据与有：

等价于

因l_i与l_i+1线性无关，因此有

全凸性：

参数曲线C(u)，若满足当u递增时，曲线是逆时针方向，则相对曲率可表示为：

则当若κ≥0，区域是凸区域，若κ<0，区域是非凸区域。

因此，对于本发明构建的曲线，只需要证明任意点处的κ≥0。

根据相对曲率的计算公式，因其分母恒大于0，因此对于一段曲线C_i，若曲线C_i满足公式则曲线C_i是凸曲线。

有

进一步化简有：

根据公式有

上述有|l_i,l_i+1|>0，α_i/β_i>0，e_i>0且有μ>0与μ_i>0，因此，上述3个式子均大于0，因此有|C_i′(u),C_i(u)|≥0，即曲线C_i是凸曲线，从而闭合曲线是全凸曲线。

实施例2：

在林业调查中，经常使用钢制围尺测量树干的直径，围尺在人为拉力的作用下环绕树干横断面一周，从而读出树干的直径值。钢制围尺具有弹性，在拉力的作用下，围尺轨迹建立在树干横断面的凸出部分上，其是一个全凸的且连续的闭合轨迹，此时，围尺轨迹可简化为一条光滑的全凸曲线。将本发明的全凸闭合曲线的构建方法应用于从地面三维激光扫描仪获取的树干点云中提取树干直径，以外业调查时获取的31条树干直径数据作为真实值，将不同μ值下构建的围尺轨迹提取的直径与实测直径相比，得出本次外业调查最优的μ值为0.3。再分别使用拟合圆、凸包折线与本发明中的模拟围尺测径的方法提取树干直径，与外业调查的实测直径相比，3种方法提取树干直径的精度比较如下表所示。

表1不同方法实验结果对照比

表中评测指标的含义及计算公式如下：平均绝对百分比误差MAPE(mean absolutepercentage error)、均方根误差RMSE(root mean square error)和决定系数R²；对应计算公式如下：

n为样本单元，D是直径的实测值，其平均值用表示，表示从树干点云提取的直径值。上述评价指标中MAPE与RMSE越小越好(R²则越大越好)，这说明从树干点云中提取直径的准确性越高。

在围尺测量直径时，围尺受力越大，围尺环绕树干横断面越紧，从而树干横断面周长越小，树干直径也最小。根据上述，μ值越大，曲线的长度越小，构建曲线的μ参数恰好反映了曲线在构建时所受的拉力情况。由此，μ值可认为是闭合曲线的拉力系数，能反映曲线在插值点受力与曲线的紧绷程度。

如图5所示，本实施例中一种二阶几何连续的全凸闭合曲线的构建方法，包括如下步骤：

步骤S201，获取凸包点；

步骤S202，构建凸包点集P＝{p_i|i＝0,1,...,n-1}，且有当i≤0,或i>n时，i表示为i对n取模，如且p₀,p₁,...,p_n-1是按照逆时针次序排列的点集；

步骤S203，采用分段三次Bezier曲线构建全凸闭合曲线；

步骤S204，构建并求解方程组；

步骤S205，解的优化及全凸闭合曲线的计算，得到由分段三次Bezier曲线构成的二阶几何连续的全凸闭合曲线。

在一个实施例中，步骤S201，获取凸包点，包括以下步骤：

步骤1，从地面三维激光扫描仪获取的树木点云中提取树干点云；

使用地面三维激光扫描仪配备的点云处理软件对地面三维激光扫描仪在多站扫描的点云进行配准得到单木的点云数据，经移除地表点云、移除树枝点云、移除离散噪声点后得到树干点云。

步骤2，计算树干在高度h处树干横断面的定位点。

根据z轴的坐标值，计算树干点云中的最低点，并以此点为基础计算树干的高度；以最低点为基准，以向量为平面的法向量，在最低点的上方h处与h+0.5厘米处分别构建的两个平面，树干高度h处的平面如图6所示。图6中，虚线2为计算定位点时的下平面，点1为定位点，虚线3为树干在高度h处的的横断面，箭头所指方向为树干横断面的法向量。图7是获取用于计算树干直径的点云示例图。图7中，虚线2是树干的横断面，虚线1是与虚线2平行且到虚线2所在平面距离为1厘米的一个平面，位于虚线1与虚线2中间的树干点云为即为用于计算树干直径的点云。位于两个平面中的树干点云构成了一个树干点云块，将此部分点云投影到树干高度h处的平面上得到一个平面点云集，计算此平面点云集构成的凸包多边形的质心点，质心点即为此平面点云集的几何中心点，同时也是树干在高度h处树干横断面的定位点。

步骤3，计算树干在高度h处的生长方向与树干横断面。

树干在高度h处的生长方向即是树干在高度h处的树干横断面的法向量。采用迭代方法计算树干在高度h处的生长方向，第k次迭代计算为：构建过定位点且法向量为的一个平面，称此平面为基准平面。

以基准平面为基础，分别构建5个与此平行的平面，其中3个平面位于基准平面上方，2个平面位于基准平面下方，且相邻平面间的距离为0.5厘米；由此得到5个树干点云块与5个几何中心点，对5个几何中心点采用主成分分析方法计算树干在高度h处的生长方向其中，协方差矩阵的最大特征值对应的特征向量为树干在高度h处的生长方向平面法向量赋值为进行第k+1次迭代；初次迭代时法向量为迭代结束条件为：相邻两次迭代得到的两个树干生长方向与的夹角θ_k小于0.5度或者θ_k与θ_k+1的角度差小于0.5度；迭代结束时的树干生长方向即为树干在高度h处的生长方向；以树干在高度h处的生长方向为平面的法向量，过定位点构建树干在高度h处的树干横断面。

步骤4，获取用于计算树干直径的凸包点。

在得到树干在高度h处的树干横断面的基础上，在此树干横断面上方构建一个与树干横断面平行且距离为1厘米(围尺的宽度)的平面，位于树干横断面与平面间的树干点云即为用于计算树干直径的树干点云，将此部分点云投影至树干横断面上得到一个三维空间上的平面点云集，采用几何变换中的旋转操作，将此平面点云集旋转为与水平面平行的一个二维空间上的平面点云集，根据凸包算法计算此平面点云集的凸包点集P。

需要说明的是，以上实施例仅用以说明本发明的技术方案而非限制，本领域普通技术人员对本发明的技术方案所做的其他修改或者等同替换，只要不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。