CN105405162A - 基于局部结构和方向感知的树点云三维重建方法 - Google Patents

Abstract

本发明实现了<b>基于局部结构和方向感知的树点云三维重建方法</b>，属于空间信息技术领域。考虑到局部结构内树枝方向和点云密度具有连续性的特点，该方法推演出了数据缺失处树枝的连接关系，设计了点云修复的优化方程和迭代过程，利用骨架点云驱使点云数据中的点沿着骨架线方向边收缩边扩散，使得优化后的点云不会产生额外的噪声，不依赖于优化后的点云而采用收缩的骨架继承相邻迭代之间的空间关联，充分保留了原始点云的结构信息。采用植物生长模型计算树骨架上每一节点的半径，根据每个节点的半径，将骨架膨胀成三维树模型。与已有树木三维重建方法相比，该方法更能精确的修复缺失区域的点云，获得更好的建模结果。

Description

基于局部结构和方向感知的树点云三维重建方法

一、技术领域

涉及基于局部结构和方向感知的树点云重建方法，属于空间信息技术领域。

二、背景技术

近年来激光雷达技术得到了快速发展，它通过高速激光扫描测量，高精度的快速获取大面积被测对象表面的三维坐标数据，成为快速获取精细的、大规模三维城市场景数据的重要手段。地面激光雷达是一种非破坏性的高分辨率三维测量手段，可弥补现有观测手段的不足，实现植被等几何结构参数的自动获取，为重建树木的三维真实结构提供了可能。

由于测量条件的限制导致测量数据经常存在较高噪声，同时周围场景地物的复杂性也决定了测量过程中的遮挡不可避免，所以测量得到的数据通常是不完整的，存在部分点云缺失，再加上，树木一般具有复杂的三维结构、拓扑关系和丰富的细节，这就决定了基于单方向扫描的地面激光雷达树木点云的重建具有极大的挑战性，有许多关键技术问题需要解决。

在数据缺失区域的点云修复和枝干结构连接方式上，预先无法确定树枝相互间的连接关系是否正确，传统的树木建模方法通常使云中的每点向着树枝两端运动，进行数据修复，造成树的骨架几何拓扑关联产生较大畸变。事实上，树枝之间在局部结构上存在一定的规律性，在树点云三维重建中，可以借助树的局部结构表达树枝的延伸和逻辑关系，基于此，申请人给出了一种能描述树枝的角度以及每段树枝点密度的树枝局部结构的表达方法，让每段树枝局部结构点云中的点沿着主方向向该树枝两端移动，由缺失处两端点云共同负责其修复，并且点的移动方式是依据局部结构中点云的分布，不受周围其它树枝点云是否缺失的影响，所以减少了周围点云噪声、点密度不均一或点云缺失对该树枝提取的影响，极大提高了树三维重建的精度。提出的方法在植被遥感、林业调查、数字城市等方面具有广阔的应用前景。

三、发明内容

1、目的：树木作为自然界中最普遍存在的景观之一，构造其三维模型已成为现实世界三维数字化工作中不可或缺的一部分。由于树木普遍具有复杂的几何形状和拓扑形态，构建树木三维模型一直以来具有很大的挑战性，尤其是重建真实的三维树木模型。申请人基于局部结构内树枝方向和点云密度具有连续性的特点，推演出缺失处树枝的连接关系，设计了点云修复的优化方程和迭代过程，充分保留了原始点云的结构信息，获得了精确的三维建模结果，满足树木三维建模在植被遥感、林业调查、数字城市等方面的需求。

2、技术方案：

基于局部结构和方向感知的树点云重建方法，其特征在于，包括如下步骤：

步骤一：初步提取树的骨架

用K最邻近结点算法将树点云中每个点与其周围的k个点连接起来形成一个权重图G(V，E)，V表示顶点集合，E表示边的集合。G中各个联通分量之间最近的两个点依次连接起来，以此提取初始树骨架、每一点的主方向以及每一段树枝上包含的点。

步骤二：树点云收缩形成点云骨架，获取骨架上点的主方向

由于树木点云中部分区域往往存在着数据缺失，所以会造成权重图G构建的树骨架在这些区域存在着错误的连接关系，但是在没有数据缺失的区域里，树枝连接方式还是相似于真实的树枝形状，利用这些正确的树枝连接方式对数据缺失区域的枝干进行修复。基于已有的骨架以及原始点云中每一点的主方向，采用拉普拉斯算法对点云收缩，形成骨架点云，减弱原始点对树枝方向的干扰。为了确保收缩后的点云真实反映原始骨架，点云应向着已知的骨架进行收缩，并且收缩到骨架线上。采用的拉普拉斯算子如方程(1)所示。

$L_{ij}=\left\{\begin{array}{cc}{\mathrm{\&omega;}}_{ij}=e^{-d(i,j)}\&times;tg\mathrm{\&theta;}& if(i,j)\&Element;E\\ {\&Sigma;}_{(i,k)\&Element;E}^k-{\mathrm{\&omega;}}_{ik}& ifi=j\\ 0& otherwise\end{array}\right.---\left(1\right)$

在(1)式中，L_ij是拉普拉斯算子L的一个元素，ω_ij是一个变量，点i、j、k分别是树上的三个点，θ是i点与其在骨架线段上投影的连线和边d(i，j)的夹角，d(i，j)是点i、j之间的欧式距离，E是边的集合。

W_ij中加入tgθ后，点云中的点会朝着骨架线的方向进行收缩。求解拉普拉斯方程(2)，获得收缩后的骨架点云。

$\left[\begin{array}{c}{W}_{L}L\\ W_H\end{array}\right]V^{\&prime;}=\left[\begin{array}{c}0\\ W_{H}V\end{array}\right]---\left(2\right)$

在(2)式中，W_L和W_H是对角矩阵，W_L控制平滑的程度，W_H控制和原始点相似的程度，W_L中第i个对角元素的的值定义成W_L，i，W_H中第i个对角元素的值定义成W_H，i。V’是收缩后的点云，V是原始的点云。

点云收缩以后，尖端的树枝结构也会收缩消失，所以需要增加尖端的权重。因为尖端的点突出于周围的点，使得尖端点和周围点的夹角比较小，而非尖端的点和周围点的夹角会比尖端点大。基于这个原因，W_L的值由下式(3)定义。

$W_{L,i}={\mathrm{\&Sigma;}}_{j\&Element;{\mathrm{\&Omega;}}_i}{\mathrm{sin\&theta;}}_2\&times;e^{-d(i,j)}/\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{e}^{-d(i,j)}---\left(3\right)$

在(3)式中，θ₂是边d(i，j)和i点主方向的夹角，Ω_i是与E中点i相邻点的集合。

经过拉普拉斯方程收缩后，骨架点云和骨架的连接方式非常相似，并且数据缺失区域没有点，保证了骨架点云不会添加入额外的噪声。

步骤三：基于树局部结构的点密度的修复

虽然骨架点云提供了良好的方向特性，但是通过骨架上每一段包含的点数量来确定点密度并不准确。这一步设计局部结构，利用树本身存在的重复性结构对点密度进行估算，获得更加精确的点密度。初始点密度由初始骨架获得，表示一段骨架之间包含的点云数量。首先定义和描述树木的局部结构，使得该结构尽可能真实的表达不同形态局部树枝空间几何拓扑关系；接下去，采用稀疏编码学习和重建该局部结构，实现点密度的修复，使得点密度更加精确，以此进行精确点云修复。

通过有向树结构描述一棵树的骨架，树形结构上每一条边的方向表示枝条方向。有向树一条边的两端点中，靠近根节点的端点称之为前向点，远离根节点的端点作为后向点，该边的方向是从前向点指向后向点，骨架上每一条边与其首尾相连的边作为该边的局部结构。定义了边的局部结构后，就需要表达每一个局部结构的特征，过程如下：

(1)树枝之间的角度具有一定的规律性，这些角度用来作为局部结构的特征。

计算每个局部结构中衔接边与所有承接边之间的夹角。因为夹角在0-180°之间，所以每10°作为一个采样间隔，每个间隔记录夹角在该区域的边，通过这种方法构造这个局部结构的特征向量。为了统计每个树枝局部结构前向和后向的树枝分布，需要区分前向和后向的树枝，分开进行编码记录。这个特征向量是一个36维的向量，前18维存储了前向的承接边和衔接边的关系，后18维存储了后向的承接边和衔接边的关系。

(2)一个局部结构的向量除了记录角度关系外，也记录与树枝半径相关的信息，而在点云中相关量为局部点密度，每一维中记录这个方向间隔中所有枝的平均点数量。

这样一个局部结构就被编码成一个特征向量，通过这个特征向量来描述该特征。此特征向量跟该结构衔接边树枝的角度及其旋转是无关的，可以描述树上不同方向以及不同旋转过后的局部结构。

由于获得了树上每一条边对应的局部结构，通过这些局部结构训练计算过冗余字典。为了让字典结果稳定，连接叶节点以及根节点的边，即前18个向量或者后18个向量均为0的特征向量，不进行字典训练，采用k-SVD方法从这些向量中获得自适应过冗余字典。

获得过冗余字典以后，任意的一个局部结构的特征向量α1，可以被过冗余字典重构获得一个新的特征向量α2，用α2中对应的值代替α1中的值，作为重建后的点密度。为了让点密度具有连续性和减少初始点归属的影响，在进行稀疏字典训练和重建前，采用拉普拉斯算法平滑点密度。

步骤四：修复缺失的树点云

把上一步收缩得到的点云和原始的点云合并起来，形成一个新的点云P。为了防止点云修复过程中造成点的扩散，每一点在沿着主方向运动的同时，沿着轴向进行收缩。最小化下式(4)和(5)，对缺失点云进行修复。

如果点i不在骨架上，

$\left|\right|\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{O}_i^T\&times;\frac{P_i-P_j}{|P_i-P_j|}\&times;O_i\&times;e^{-|P_i-P_j|}\&times;{\mathrm{\&rho;}}_j/\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{e}^{-|P_i-P_j|}+\mathrm{\&lambda;}\&times;{\mathrm{\&rho;}}_{max}((U_i-P_i)+\mathrm{\&phi;}\&times;(P_{si}-P_i))|{|}^2---\left(4\right)$

如果该点在骨架上，

$\left|\right|\underset{i\&Element;\mathrm{\&Omega;}}{\&Sigma;}{Q}_i^T\&times;\frac{P_i-P_j}{|P_i-P_j|}\&times;O_i\&times;e^{-|P_i-P_j|}/\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{e}^{-|P_i-P_j|}+\mathrm{\&lambda;}\&times;(U_i-P_i)|{|}^2---\left(5\right)$

O_i是点i的主方向，U_i是原始点云中i点的位置，P_i、P_j分别表示优化后点i、j的坐标，P_si表示收缩后的点在优化点云中的坐标，p_j表示点j的点密度，λ表示向着轴向和原始点收缩的系数，Φ表示向着轴向收缩和向着原始点收缩的比例。力是从点j指向点i的，所以力的大小和j处的点密度有关，ρ_max表示点i周围点最大的那个点密度。采用ρ_max后会压制点密度小处的点运动，防止了末端点的扩散，并且鼓励点密度高处的点运动。在优化中加入了骨架点云，使每一个点在优化过程中沿着骨架点云进行运动，可以达到轴向收缩的效果。在优化过程中，骨架上的点较为密集，这些点在优化后延伸较大，通过这些点再拉动与其对应的原始点云，会让整个点云沿着各自的骨架线方向进行扩散，完成点云的修复过程。

为了防止点云沿着骨架方向收缩过大或者沿着不正确的骨架方向移动，在优化过程中控制点云移动的步长，每次只移动一小段，迭代进行修复和重建过程，直到获得良好的树模型，所以点云缺失的修复需要迭代进行。由于骨架点云是从点云沿着骨架线的收缩得到的，用以描述树木主要结构的骨架点云保留了大量的原始点云位置和结构信息，充分利用骨架点云的这个特点设计了迭代过程。

(1)当点云经过一次修复后，把原始输入点云、上一步的骨架点云叠加到优化后的点云上，形成新的树木点云。

(2)利用前一步获得的树骨架，对优化后的点云进行收缩，得到新的骨架点云。

(3)把第(2)步产生的骨架点云、原始点云叠加起来，获得这些点的连接关系，进行点云修复。修复后的点云包括了两部分：一部分是骨架点云运动后的点集，另一部分是原始点云运动后的点集。

(4)依据第(3)步获得的连接关系，对原始点云运动后的点进行收缩，得到新的骨架点云。

由第(3)步得到的连接关系可获得第(2)步骨架点云中存在的边e1，即e1的两端点连接的是第(2)步骨架点的边，计算e1的两个端点和这一步骨架点以及原始点中任两点的距离之和，搜索距离最邻近的两个点，这两个点形成边e2，将边e2的权重W_e2变成下式(6)。

W_e2＝W_e2(1-cosθ₃/(d_e1/d_e2+d_e2/d_e1))(6)

在(6)式中，d_e1是边e1在边e2上的投影长度，d_e1是边e2的长度，W_e2是边e2的权重，θ₃是边e1和e2之间的夹角。

(5)当修改完所有相关边的权重后，提取新骨架，然后转到第(2)步，比较相邻两次迭代的骨架点云末端点的移动情况，当所有末端点的运动小于阈值δ时，整个迭代结束。

步骤五：重建树的三维模型

当获得树骨架以后，采用植物生长模型计算每一个节点的半径，根据每个节点的半径，将骨架膨胀成三维树模型。

3、优点及功效：提出了树局部结构平滑点云密度的方法，基于点云密度和点的主方向，设计了数据缺失区域点云修复的优化方程和迭代过程，利用骨架点云驱使点云中的点沿着骨架线方向边收缩边扩散，使得优化后的点云不会产生额外的噪声，不依赖于优化后的点云而采用收缩的骨架继承相邻迭代之间的几何拓扑关联关系，充分保留了原始点云的结构信息。采用植物生长模型计算树骨架上每一节点的半径，根据每个节点的半径，将骨架膨胀成三维树模型。与已有树三维重建方法相比，提出的方法能更精确的修复缺失区域的点云，获得精确的三维建模结果。

四、附图说明

图1方法的流程示意图

图2存在自遮挡的树点云建模

图2(a)原始点云

图2(b)数据修复后的点云

图2(c)重建的三维树模型。

图3存在较大数据缺失的树点云建模

图3(a)原始点云

图3(b)数据修复后的点云

图3(c)重建的三维树模型。

图4树冠处点较为密集的多棵树的三维重建

图4(a)原始点云

图4(b)数据修复后的点云

图4(c)重建的三维树模型

图5手工删去部分点后的树建模和原始点云建模结果的比较

图5(a)原始树点云

图5(b)原始树点云的三维建模结果

图5(c)手工删除一部分点后的点云数据

图5(d)数据修复后的树点云

图5(e)数据缺失修复后树点云的三维建模结果。

五、具体实施方式

涉及一种基于地面激光雷达点云数据自适应构建三维树木骨架的方法，该方法具体步骤如下(如图1)：

步骤一：初步建立树木三维几何结构

对于包含多棵树的地面激光点云，需要先识别树的株数及其对应的底部。获取一定高度层的点，然后聚类提取单棵树的点云集合。计算每棵树点云集合的中点，作为该树的底点。对于复杂的区域需要交互式指定底部来确定树的棵数和底点。

在树木点云中，利用k-NearestNeighbor(KNN)算法计算每个点最邻近的n个点，一般的，n＝50；设其中两点分别为i点和j点，KNN算法计算出的最邻近点i的最相邻点包含j，而j的最邻近不包括i，所以让i的最相邻点包含j，且点的相邻点包含i，点与点连线形成边的长度作为该边的权重，以此构建一个无向图，通过该无向图构建最小生成树(MST)。用Diikstra算法计算每个底点到MST各点的最短距离，把每个点归属到离所对应底点最近的树上去。这样就获得了一株树的点云以及其连接方式。

步骤二：构建距离最小生成树

从树木点云中用Dijkstra算法获得树的三维结构缺乏相邻点的信息，而MST也不能很好描述树的伸展形状。树枝通常是向上生长的，小范围内树干的分布应该是相对光滑的，树上的点被认为是以离根节点最近的方式连接，并且越相临的点在同一树枝上的可能性越大。采用距离最小生成树方法(DMst)来构建子树。DMst的方法流程如下：

(1)对单株树的点云构建无向图，其中顶点集合为V。设每个顶点到根节点的距离D_dij，每个顶点到当前DMst(此时，DMst为空)的最小距离为D_mst，每个顶点指向为P，每个顶点的MST指向为P_mst，每个顶点Dijkstra树的指向为P_dij；

(2)把底点做为DMst的根节点，初始化V，除底点以外的所有顶点放入V中，计算所有点的D_dij和D_mst，把根节点相邻点的P、P_dij和P_mst指向根节点。

(3)在V中，选取D_mst+D_dij之和最小的点，把该点从V中移出。

(4)更新P_mst和D_mst，设点i、j之间的距离为d_ij，如果D_mst[i]＞d_ij，那么P_mst指向j，D_mst[i]＝di_j，更新P_dij和D_dij；如果D_dij[i]＞D_dij[j]+d_ij，那么P_mst指向j，D_dij[i]＝D_dij[j]+d_ij；当D_mst[i]+(D_mst[i]+D_dij[j])＜＝D_dij[j]-D_dij[k]＋D_dij[i]时，P[i]指向距树最小的位置j，否则，P[i]指向距根节点距离最小的位置。

(5)重复步骤(3)和(4)，直到V中没有点为止，通过P获取所有点的组织结构，完成DMst构建。

构建DMst时，用D_mst[i]+w×D_dij[i]代替D_mst+D_dij控制入选的点。更新时，因为在DMst中每次获取i点和j点的公共前向，计算复杂度大，而用D_mst[i]+(D_mst[i]+D_dij[j]-D_dij[k]＜＝w×(D_dij[i]-D_dij[k]+D_dij[i]-D_dij[k])作为判定条件，可以获得较好效果。

DMst获得的树形结构中，大量的树节点点只有少量的子节点，去除那些只有少量子节点的树节点，剩下的点构成一个近似的树骨架。设置每个点的权值c_i，对没有子节点的点来说，c_i＝0.1，对有子节点的点，

$c_i=\underset{j\&Element;\mathrm{\&Omega;}}{\&Sigma;}(c_j+d_{ij})---\left(1\right)$

式(1)中，Ω是i点相邻子节点的集合，j点是i点相邻的子节点，d_ij是i点和j点之间的距离。

设阈值为δ，剔除c_i＜δ的点而保留c_i＞δ的点。保留的点构成了树的近似骨架，并且骨架点是树的表面点，骨架节点之间的方向很好反映了树枝的伸展方向。DMst构成树的连接方式具有方向性，因为考虑了每个点到根部的距离，形成树枝的方向主要是向上的，符合树的生长特征，这也为缺失点云区域的修补提供了帮助。

步骤三：优化树骨架结构

一般的，树的表面是连续的，树枝的变化是平滑的，不会出现粗枝瞬间变细枝或者细枝瞬间变成粗枝的情况，所以用如下假设来描述树的表面：地面激光雷达扫描仪扫描树木表面形成点云的过程可看作连续表面离散化的过程，点云类似于粒子，它们之间有相互的斥力，所在位置对其本身有一定的束缚作用，这些粒子所构成的表面处于受力平衡的状态，因此有点云缺失的树木区域，本应该存在的表面点没有了，整个表面的受力平衡状态就破坏了，粒子会产生运动，最后会获得一个新的受力平衡状态，而这个平衡状态就是修复了缺失的状态。

(1)提取枝干主方向

因为点密集分布在树的表面，相邻点之间的方向并不能代表树枝的伸展方向，所以点与点之间的方向不能代表表面上每一个点受力方向，用主方向来反映树枝的伸展方向。

为了获得主方向，对去除c_i＜δ点的点云进行重采样。将c_i＞δ且所有后向点的c_i＜δ的点作为末端点，如果一个点有两个或两个以上后向点c_i＞δ，这个点是树的分叉点。在所有采样点构成的树骨架上，每隔距离d再采一个点，d是点与点之间树干上的距离而不是直线距离，把这些点按DMst的结构连接起来，d＝0.2。设相邻两点为i，j，离根节点近的点设为i，另一点设为j，由i指向j的方向作为这两点之间的方向，该方向作为i点的所有子节点的主方向；取末枝点的前向节点到末枝点的方向作为末枝点的所有子节点的主方向。

(2)优化点云缺失的树木

对于完整点云集合来说，其中每个点处于受力平衡状态，所以一个点所连接的所有点对该点的合力在主方向上的投影F_r与其所在位置对该点的束缚力F_s之和F＝F_r+F_s，作为该点所受的外力，显然F＝0时，该点处于平衡状态。类似于同性电荷之间的力，点之间的斥力与其距离成反比，通过公式(2)计算每个点的电荷量q_i，

$q_i=\frac{1}{r_i}---\left(2\right)$

式(2)中，q_i是对应点i的电荷量，r_i是对应点i到n个最临近点的平均距离，相连i与点j之间的斥力是，

$F_r(i,j)=\frac{q_i{q}_j}{{\left|\right|P_i-P_j\left|\right|}^2}(P_i-P_j)---\left(3\right)$

式(3)中P_i，P_j是优化后对应点i和j的坐标，q_i，q_j是对应点的电荷量。

i点所在的原始位置对其的束缚力应该和i点偏离原始位置的距离成正比。通过公式(3)求当前点i的阻尼系数K_i。

$K_i=\frac{(\underset{j\&Element;\mathrm{\&Omega;}}{\&Sigma;}{q}_j)q_i}{{r_m}^2}---\left(4\right)$

式(4)中，r_m是树上两点之间的平均距离，Ω是与i点相连接的所有点集合，j是i相连的点，i所在位置对i的束缚力为，

F_S(i)＝K_i(Q_i-P_i)(5)

式(5)中，P_i是优化后对应点i的坐标，Q_i是为优化前点i的坐标。最后的优化函数：

式(6)，Ω是i点所有相连接点的集合，j是i相连的点，O_i是i点的主方向，w是权重，w的改变可以控制点的偏离程度。

对DMst进行重采样，并且对重采样的每一个末枝点进行延伸，寻找末枝点的后向子节点中ci最大的点，作为待处理点，再搜素待处理点的后向子节点中ci最大的点，如果点到末枝点在树干上的距离超过d，那么把这个点作为新的末枝点加到骨架中，直到找到一个没有后向点的点为止，并把该点加入到骨架。这样就得到最后的树骨架。

步骤四：重建树木三维骨架

用拉普拉斯方程(7)对上面获取的树骨架进行平滑。

$L_{ij}=\left\{\begin{array}{cc}{\mathrm{\&omega;}}_{ij}=1& if(i,j)\&Element;E\\ {\&Sigma;}_{(i,k)\&Element;E}^k-{\mathrm{\&omega;}}_{ik}& ifi=j\\ 0& otherwise\end{array}\right.---\left(7\right)$

L是n×n的拉普拉斯算子矩阵，E是i相邻点的集合，j和k是i的相邻点。

$\left[\begin{array}{c}{W}_{L}L\\ W_H\end{array}\right]V^{\&prime;}=\left[\begin{array}{c}0\\ W_{H}V\end{array}\right]---\left(8\right)$

公式(8)中，W_L和W_H是对角矩阵，用来平衡收缩力和吸引力；V’是优化后的点，V是初始点；W_L对角线上的值设为1，根据需要设置平滑的大小。

因为点云数据包含噪声或存在缺失，难以确定树枝干的宽度。一般的，树底部附近的点比较密集，并且周围干扰点较少，可以估算出树底部的半径，根据树的异速生长模型计算所有树枝的半径。如果父节点只有一个子节点，那么，

$r_c=r_p{\left(\frac{l_c}{l_p}\right)}^{\frac{3}{2}}---\left(9\right)$

r_c是子节点的半径，r_p是父节点的半径，l_c是子节点的c，l_p是父节点的c。

如果一个父节点i有多个子节点，那么，

$r_{ci}=r_p{\left(\frac{l_{ci}}{\underset{j}{\&Sigma;}{l}_{cj}}\right)}^{\frac{3}{2}}---\left(10\right)$

j是子节点的个数。

至此计算了每个节点的半径。对于有叶子的树，认为树干半径小但周围点云多的位置长叶子的可能较大，所以，在优化后点云中，树骨架上的每点的q_i除以该点处树干半径r_i，对这些值进行降序排序，根据树叶的疏密，添加叶子。

实施例1：

采用RIEGLLMS-Z620地面激光雷达扫描仪单站式扫描获得树点云，扫描时的垂直分辨率和水平分辨率都是0.057°，场景树木距离扫描仪的距离远近不同，远处的树被前面的地物遮挡，造成数据缺失和点密度不均一。在一台InterCore(TM)i7-47903.6GHZ、内存4GB的计算机上对上述树点云进行三维建模。

图2(a)是距离扫描仪较近的一棵树的点云数据，树的点密度较高，但是树冠部分树枝较多且细长，自遮挡造成了局部区域点密度的跳变，并且背向扫描仪的一面点云极为稀疏。图2(b)是对原始点云修复后的新点云数据，从该数据中重建了该树的三维模型(图2c)。从图2(c)可以看出，重建后的树干形状得到较好的保持，没有受到自遮挡的影响。

如图3(a)给出了数据缺失严重的树点云。红框内树枝点云缺失了，这棵树距离扫描仪较远，整棵树的点密度不高。因为迭代的进行修复，收缩得到骨架点云，很好的保持了原始点云中较明显的主干结构，点云修复后(图3b)可以实现较好的建模结果(图3c)。

图4(a)是相互间距离较近的三棵树点云，由于每颗树冠覆盖面积较大，所以在树冠部分处点云数据相互重叠，并且树冠处点较为密集，因为扫描仪前面地物以及自身相互遮挡，造成树冠具有条带状的缺失区域。采用提出的方法，数据缺失区域的点云能够较好的恢复(图4b)，重建后的三维树模型的树枝形状保持的十分完好(图4c)。

图5(b)是基于图5(a)所示的完整树点云的三维重建结果，在图5(a)的红框中树枝延伸和彼此间的连接关系较为复杂，为了验证申请人的方法，申请人删除了红框内的点(如图5c)，通过点云修复最终得到了图5(d)所示的点云，红框区域处的数据得到了较好的修复，图5(e)是利用该点云数据重建的三维模型，两种数据的建模结果具有很好的相似性。

Claims (1)

1.基于局部结构和方向感知的树点云重建方法，其特征在于，包括如下步骤：

步骤一：初步提取树的骨架

用K最邻近结点算法将树点云中每个点与其周围的k个点连接起来形成一个权重图G(V，E)，V表示顶点集合，E表示边的集合，G中各个联通分量之间最近的两个点依次连接起来，以此提取初始树骨架、每一点的主方向以及每一段树枝上包含的点；

步骤二：树点云收缩形成点云骨架，获取骨架上点的主方向

基于已有的骨架以及原始点云中每一点的主方向，采用拉普拉斯算法对点云收缩，形成骨架点云；点云向着已知的骨架进行收缩，并且收缩到骨架线上，采用的拉普拉斯算子如方程(1)所示：

$L_{ij}=\left\{\begin{array}{cc}{\mathrm{\&omega;}}_{ij}=e^{-d(i,j)}\&times;tg\mathrm{\&theta;}& if(i,j)\&Element;E\\ {\mathrm{\&Sigma;}}_{(i,k)\&Element;E}^k-{\mathrm{\&omega;}}_{ik}& ifi=j\\ 0& otherwise\end{array}\right.---\left(1\right)$

在方程(1)中，L_ij是拉普拉斯算子L的一个元素，ω_ij是一个变量，点i、j、k分别是树上的三个点，θ是i点与其在骨架线段上投影的连线和边d(i，j)的夹角，d(i，j)是点i、j之间的欧式距离，E是边的集合；

W_ij中加入tgθ后，点云中的点会朝着骨架线的方向进行收缩；求解拉普拉斯方程(2)，获得收缩后的骨架点云；

$\left[\begin{array}{c}{W}_{L}L\\ W_H\end{array}\right]V^{\&prime;}=\left[\begin{array}{c}0\\ W_{H}V\end{array}\right]---\left(2\right)$

在方程(2)中，W_L和W_H是对角矩阵，W_L控制平滑的程度，W_H控制和原始点相似的程度；W_L中第i个对角元素的的值定义成W_L，i；W_H中第i个对角元素的值定义成W_H，i；V’是收缩后的点云，V是原始的点云；

W_L的值由下式(3)定义；

$W_{L,i}={\mathrm{\&Sigma;}}_{j\&Element;{\mathrm{\&Omega;}}_i}{\mathrm{sin\&theta;}}_2\&times;e^{-d(i,j)}/\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{e}^{-d(i,j)}---\left(3\right)$

在方程(3)中，θ₂是边d(i，j)和i点主方向的夹角，Ω_i是与E中点i相邻点的集合；

步骤三：基于树局部结构的点密度的修复

通过有向树结构描述一棵树的骨架，树形结构上每一条边的方向表示枝条方向；有向树一条边的两端点中，靠近根节点的端点称之为前向点，远离根节点的端点作为后向点，该边的方向是从前向点指向后向点；骨架上每一条边与其首尾相连的边作为该边的局部结构，表达每一个局部结构的特征，过程如下：

(1)树枝之间的角度具有一定的规律性，这些角度用来作为局部结构的特征；

计算每个局部结构中衔接边与所有承接边之间的夹角，每10°作为一个采样间隔，每个间隔记录夹角在该区域的边，通过这种方法构造这个局部结构的特征向量；为了统计每个树枝局部结构前向和后向的树枝分布，需要区分前向和后向的树枝，分开进行编码记录；这个特征向量是一个36维的向量，前18维存储了前向的承接边和衔接边的关系，后18维存储了后向的承接边和衔接边的关系；

(2)一个局部结构的向量除了记录角度关系外，也记录与树枝半径相关的信息，而在点云中相关量为局部点密度，每一维中记录这个方向间隔中所有枝的平均点数量；

这样一个局部结构就被编码成一个特征向量，通过这个特征向量来描述该特征，采用k-SVD方法从这些向量中获得自适应过冗余字典；获得过冗余字典以后，任意的一个局部结构的特征向量α1，可以被过冗余字典重构获得一个新的特征向量α2，用α2中对应的值代替α1中的值，在进行稀疏字典训练和重建前，采用拉普拉斯算法平滑点密度；

步骤四：修复缺失的树点云

把步骤三收缩得到的点云和原始的点云合并起来，形成一个新的点云P；为了防止点云修复过程中造成点的扩散，每一点在沿着主方向运动的同时，沿着轴向进行收缩；最小化下式(4)和(5)，对缺失点云进行修复；

如果点i不在骨架上，

$\left|\right|\underset{j\&Element;Q_i}{\&Sigma;}{O_i}^T\&times;\frac{P_i-P_j}{|P_i-P_j|}\&times;O_i\&times;e^{-|P_i-P_j|}\&times;{\mathrm{\&rho;}}_j/\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{e}^{-|P_i-P_j|}+\mathrm{\&lambda;}\&times;{\mathrm{\&rho;}}_{max}\left(\right(U_i-P_i)+\mathrm{\&phi;}\&times;(P_{si}-P_i\left)\right)|{|}^2---\left(4\right)$

如果该点在骨架上，

$\left|\right|\underset{i\&Element;Q}{\&Sigma;}{O_i}^T\&times;\frac{P_i-P_j}{|P_i-P_j|}\&times;O_i\&times;e^{-|P_i-P_j|}/\underset{j\&Element;{\mathrm{\&Omega;}}_i}{\&Sigma;}{e}^{-|P_i-P_j|}+\mathrm{\&lambda;}\&times;\left(U_i-P_i\right)|{|}^2---\left(5\right)$

O_i是点i的主方向，U_i是原始点云中i点的位置，P_i、P_j分别表示优化后点i、j的坐标，P_si表示收缩后的点在优化点云中的坐标，ρ_j表示点j的点密度，ρ_max表示点i周围点最大的那个点密度；λ表示向着轴向和原始点收缩的系数，φ表示向着轴向收缩和向着原始点收缩的比例；

在优化过程中控制点云移动的步长，每次只移动一小段，迭代进行修复和重建过程，直到获得良好的树模型，充分利用骨架点云的这个特点设计了迭代过程：

(1)当点云经过一次修复后，把原始输入点云、上一步的骨架点云叠加到优化后的点云上，形成新的树木点云；

(2)利用前一步获得的树骨架，对优化后的点云进行收缩，得到新的骨架点云；

(3)把第(2)步产生的骨架点云、原始点云叠加起来，获得这些点的连接关系；

(4)依据第(3)步获得的连接关系，对原始点云运动后的点进行收缩，得到新的骨架点云；

从第(3)步得到的连接关系中，可获得第(2)步骨架点云中存在的边e1，计算e1的两个端点和这一步骨架点以及原始点中任两点的距离之和，搜索距离最邻近的两个点，这两个点形成边e2，将边e2的权重w_e2变成下式(6)；

w_e2＝w_e2(1-cosθ₃/(d_e1/d_e2+d_e2/d_e1))(6)

在方程(6)中，d_e1是边e1在边e2上的投影长度，d_e1是边e2的长度，w_e2是边e2的权重，θ₃是边e1和e2之间的夹角；

(5)当修改完所有相关边的权重后，提取新骨架，然后转到第(2)步；当所有末端点的运动小于阈值δ时，整个迭代结束；

步骤五：重建树的三维模型

当获得树骨架以后，采用植物生长模型计算每一个节点的半径，根据每个节点的半径，将骨架膨胀成三维树模型。