CN103745459A - 一种非结构化点云特征点检测方法及其提取方法 - Google Patents

Abstract

本发明提供了一种非结构化点云特征点检测方法及其提取方法，（1）计算采样点在不同尺度空间中的Harris响应值；（2）选取最优尺度空间下的Harris响应值作为采样点的Harris响应值，得到特征点集合Q；（3）将Harris响应值在尺度空间邻域和几何邻域都具有极大性的点作为候选特征点，最后选择优化策略提取最终的特征点。在极坐标系下对得到的特征点所在切平面进行网格化分割，然后将特征点的邻域点投影到该切平面上，将每个网格中的投影点对应的投影长度向网格的四个顶点投票生成特征信息统计矩阵，然后分别对特征统计的行向量和列向量进行DCT变换和DFT变换，变换后的矩阵的左上角元素即特征描述向量。

Description

一种非结构化点云特征点检测方法及其提取方法

技术领域

本发明属于三维空间目标识别领域，涉及一种非结构化点云特征点检测方法及其提取方法。具体涉及一个面向非结构化三维点云数据的多尺度特征点检测算法和一个基于形状信息统计和空间变换思想的特征点描述算法。

背景技术

随着三维建模技术的大量普及，三维点云数据被广泛应用于文物保护、空间目标识别等诸多领域。面对信息量庞大的三维点云，如何提取符合实际应用要求的、有意义的信息是处理三维点云数据必须要解决的问题。

点云数据特征点提取作为三维点云数据处理中的关键技术是目前点云处理中的研究热点。现有的三维模型特征提取算法主要针对网格数据，无法满足对非结构化点云特征提取的要求。与网格化点云数据相比，非结构化点云数据量大、无点云拓扑信息、特征点检测容易受噪声影响，因此特征点检测和描述难度大。

发明内容

本发明提供了一种非结构化点云特征点检测方法及其提取方法，满足对非结构化点云特征的检测及提取要求。

为解决以上技术问题，本发明采用以下技术方案：

一种非结构化点云特征点检测方法，包括以下步骤：（1）利用采样点不同尺度空间的邻域信息计算采样点在不同尺度空间中的Harris响应值；（2）通过迭代算法选取最优尺度空间下的Harris响应值作为采样点的Harris响应值，得到特征点集合Q；（3）将Harris响应值在尺度空间邻域和几何邻域都具有极大性的点作为候选特征点，最后利用候选特征点选择优化策略提取最终的特征点。作为本发明的优选实施例，所述步骤（1）中，采样点在不同尺度空间中的Harris响应值根据以下方法计算：（1.1）设N_r(p_i)是p_i的邻域点集合，p_j∈N_r(p_i)，分别在p_j的[1,0,0]方向和[-1,0,0]方向确定两个点p_j+和p_j-，计算p_j+的邻域点集合N_r(p_j+)和p_j-的邻域点集合N_r(p_j+)之间的Hausdorff距离；其中，i、r、j均为正整数；（1.2）根据步骤（1.1）的Hausdorff距离计算p_j在X轴方向的偏导数f_x；（1.3）计算p_j在Y方向的偏导数f_y；（1.4）根据p_j在X轴方向和Y方向的偏导数f_x和f_y计算每个邻域点的高斯权重，即E矩阵；（1.5）根据E矩阵计算p_j的Harris响应值，计算公式为：Harris(x,y)=det(E)–k·(Trace(E))²，k为经验系数。

所述步骤（1.1）的Hausdorff距离根据以下公式计算：

$d_H(H_r\left(p_{j+}\right),N_r\left(p_{j-}\right))=\max \left\{\underset{a\∈N_r\left(p_{j+}\right)}{\sup }\underset{b\∈N_r\left(p_{j-}\right)}{\inf }\right|a-b|,\underset{b\∈N_r\left(p_{j+}\right)}{\sup }\underset{a\∈N_r\left(p_{j-}\right)}{\inf }|a-b\left|\right\}.$

所述（1.4）的高斯权重E矩阵根据以下公式计算：

$E=\frac{1}{\sqrt{2\mathrm{\π\σ}}}\left[\begin{array}{cc}\underset{x_i,y_i,z_i}{\mathrm{\Σ}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\·{\mathrm{\σ}}^2}\·f_x^2& \underset{x_i,y_i,z_i}{\mathrm{\Σ}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\·{\mathrm{\σ}}^2}\·f_x\·f_y\\ \underset{x_i,y_i,z_i}{\mathrm{\Σ}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\·{\mathrm{\σ}}^2}\·f_x\·f_y& \underset{x_i,y_i,z_i}{\mathrm{\Σ}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\·{\mathrm{\σ}}^2}\·f_y^2\end{array}\right],$

其中，(x_i,y_i,z_i)∈N_r(p_i)，σ为尺度因子，第k个尺度因子根据以下公式计算：

${\mathrm{\σ}}_k=\frac{1}{\left|N_r^k\left(p_i\right)\right|}\underset{j=1}{\overset{\left|N_r^k\left(p_i\right)\right|}{\mathrm{\Σ}}}{\left|\right|p_i-p_i^j\left|\right|}_2.$

所述步骤（2）的方法为：首先选取一个K邻域点集合作为初始尺度空间大小，计算初始尺度因子σ₀；如果采样点p_i在当前尺度k下计算出的响应值h^k(p_i)不足以判定点p_i是否是特征点，那么尺度因子将更新到下一个尺度，继续计算下一个尺度p_i的响应值；随着尺度因子的增加，响应值变化最大的那个尺度称为最优的尺度；当迭代算法停止，将上一尺度及其响应值作为最优尺度和最优的响应值；具体方法为：（2.1）令k=0，特征点集合为空，其中，k为尺度因子索引；（2.2）判断σ₀<σ_max是否成立，其中，σ₀为初始尺度因子，σ_max为最大尺度因子，如果成立，则提取p_i的邻域点集合

其中，

其中，p_i为点云数据集P中的任意一点；（2.3）计算p_i的法向量

并建立局部坐标系；（2.4）将p_j的坐标平移、旋转变换转化为步骤（2.3）中建立的局部坐标系中的坐标，其中，

（2.5）计算p_j在X轴方向和Y轴方向的偏导数，更新E矩阵；（2.6）计算k尺度下点p_i的Harris特征值h^k(p_i)；（2.7）判断h^k(p_i)>τ·h^k-1(p_i)是否成立，如果成立，令h^k(p_i)=h^k-1(p_i)，将p_i，添加到Q中；（2.8）重复以上步骤，得到特征点集合Q。

所述步骤（3）中，候选特征点的优化策略方法为：（3.1）初始化后，判断特征点集合Q是否为空，如果不为空，则进行步骤（3.2）；（3.2）如果∈Q并且

则将q_i加入π中，即π=π∪q_i，其中，q_i∈Q，

为所有与q_i相连的边的最小权值，π为特征点选择优化后的特征点集合。

一种基于以上非结构化点云特征点检测方法对非结构化点云进行提取的方法，在极坐标系下对特征点所在切平面进行网格化分割，然后将特征点的邻域点投影到该切平面上，将每个网格中的投影点对应的投影长度向网格的四个顶点投票生成特征信息统计矩阵，然后分别对特征统计的行向量和列向量进行DCT变换和DFT变换，将由变换后得到的矩阵的左上角元素构成的向量作为该特征点的特征描述向量，所述左上角为A×A，其中，1<A<10。

利用双线性插值将投影长度向特征统计矩阵的对应位置进行投票的方法为：定义一个N×M的实数矩阵S，划分的网格在半径方向将半径划分为M等份，每等份长度用binR表示，则binR=r(k)/M，在每一个同心圆上，按逆时针顺序将圆周角划分为M等份，每等份的角度用binθ表示，则binθ=2π/N；记第j个特征点f_i ^j在f_i切平面上的投影点为

投票方法为：（A）计算f_i ^j^与坐标原点的连线从x方向按逆时针旋转到当前位置的旋转角度，记为thetaR；判断

是否成立，若成立，令

否则，令 $\mathrm{thetaR}=2\mathrm{\&pi;}-{\cos }^{-1}\left(\frac{{f_i}^{j^}\&CenterDot;u}{{\left|\right|{f_i}^{j^}\left|\right|}_2\&CenterDot;{\left|\right|u\left|\right|}_2}\right);$ （B）计算行索引h和列索引l，其中，

（C）计算行方向的投票权重a和列方向的投票权重b，其中， $a=\frac{\mathrm{thetaR}-h\&CenterDot;\mathrm{bin\&theta;}}{\mathrm{bin\&theta;}},b=\frac{{\left|\right|{f_i}^{j^}\left|\right|}_2-j\&CenterDot;\mathrm{binR}}{\mathrm{binR}};$ （D）利用双线性差值和当前投影点周围的四个网格节点的投票权重计算当前投影点的位置；（E）根据步骤（D）计算投影点周围的四个网格节点形成的区域内的投影点个数；（F）当所有投影点投票结束后，对S矩阵中的元素进行平均，得到特征统计矩阵。

步骤（D）根据以下公式计算： $S(h,l)=S(h,l)+(1-a)(1-b)\&CenterDot;{\mathrm{proj}}_i^j;$ $S(h,l+1)=S(h,l+1)+(1-a)b\&CenterDot;{\mathrm{proj}}_i^j;$ $S(h+1,l)=S(h+1,l)+a(1-b)\&CenterDot;{\mathrm{proj}}_i^j;$ $S(h+1,l+1)=S(h+1,l+1)+\mathrm{ab}\&CenterDot;{\mathrm{proj}}_i^j;$ 式中：

为点f_i ^j的法向量

在点f_i所在切平面的投影。

所述DCT变换和DFT变换的方法为：首先对特征统计矩阵S的行向量进行DCT变换，用DCT变换后的行向量替换原始S矩阵中的行向量，得到新的S矩阵，然后对S矩阵的列向量进行DFT变换，DFT变换后的结果是复数，将变换后复数的模作为新的列向量元素。

本发明所提供的非结构化点云特征点检测方法及其提取方法至少具有以下优点：（1）能够直接应用于非结构化点云数据，完成非结构化点云特征的检测和描述；（2）检测效率高，本发明所提供的方法不受点云大小的限制；（3）算法鲁棒性强，本发明所提供的方法受噪声干扰的影响较小，结果准确、稳定。

附图说明

图1是本发明非结构化点云特征点的检测方法流程图；

图2是本发明非结构化点云特征点的描述方法流程图；

图3是本发明特征点的投影和投票示意图；

图4是本发明投票权重计算示意图。

具体实施方式

本发明提供了一种非结构化点云特征点提取方法，设计并实现了面向非结构化点云的特征检测和描述算法。具体内容如下：

一是设计并实现了面向非结构化三维点云数据的多尺度3D Harris特征点检测算法。通过利用采样点不同尺度空间的邻域信息计算采样点在不同尺度空间中的Harris响应值，然后通过迭代算法选取最优尺度空间下的Harris响应值作为采样点的Harris响应值，得到特征点集合Q，然后将Harris响应值在尺度空间邻域和几何邻域都具有极大性的点作为候选特征点，最后利用特征点选择优化策略提取最终的特征点。根据实际应用需求的不同，特征点选择优化策略可以采用最小生成树算法或者聚类算法。

二是设计并实现了ST-RIFT(Spatial Transform-Rotation Invariant FeatureTransform)特征点描述算法。首先，在极坐标系中将特征点所在切平面进行网格化划分，然后将特征点的邻域点投影到特征点所在切平面，然后利用双线性插值将领域点在切平面的投影长度向投影点所在网格的角点进行投票，统计采样点所在局部区域投影分布信息，得到特征统计矩阵。然后对特征统计矩阵的每一行进行DCT(Discrete Cosine Transform)变换，再对特征统计矩阵的每一列进行DFT(Discrete Fourier Transform)变换。最后，取变换后的特征统计矩阵左上角区域元素作为特征点的特征描述向量。

本发明面向非结构化点云数据提出了一种有效的三维特征点检测算法和一种特征点描述算法。具体技术包括：

1）利用Hausdorff距离对采样点邻域点云集合相似性测量，代替采样点的Harris响应值计算中偏导数的计算。

2）计算采样点在多个尺度空间中的Harris响应值，采用迭代算法提取采样点最优尺度下的Harris响应值。

3）根据不同的应用需求，候选特征点优化策略可以采用基于MST算法的优化策略或基于欧氏距离的聚类策略。前者检测到的特征点主要分布在三维点云的边缘等表面变化明显的区域；后者检测到的特征点均匀分布在三维点云的表面上。

4）采样点邻域信息统计策略。在极坐标系中对特征点所在切平面进行网格化划分，将特征点邻域点投影到特征点所在切平面，统计所有领域点的投影信息，得到特征统计矩阵。

5）对特征统计矩阵的行向量进行DCT变换，再对特征统计矩阵的列向量进行DFT变换。

按照上述技术方案，本发明给出了以下的实施例。

图1给出了多尺度3D Harris特征点检测算法流程。

设计思想是利用不同尺度下采样点的邻域信息计算多个尺度下采样点的Harris响应值，然后利用迭代算法的思想提取最优尺度下的Harris响应值作为采样点的Harris响应值。将Harris响应值在几何邻域和尺度邻域中具有局部极大性的采样点作为候选特征点。在候选特征点选择优化策略中，可以由用户根据实际应用的不同选择MST算法或聚类算法对候选特征点进行选择优化，提取最终的特征点。

图2给出了ST-RIFT特征描述算法流程。设计思想是首先在极坐标系下对特征点所在切平面进行网格化分割，然后将特征点的邻域点投影到该切平面上，将每个网格中的投影点对应的投影长度向网格的四个顶点投票生成特征信息统计矩阵，然后分别对特征统计的行向量和列向量进行DCT变换和DFT变换，将由变换后得到的矩阵的左上角元素构成的向量作为该特征点的特征描述向量。

多尺度3D Harris特征点检测算法设计，具体包括：

多尺度3D Harris响应值计算。即在采样点的多个尺度空间中分别计算Harris响应值。设N_r(p_i)是p_i的邻域点集合，p_j∈N_r(p_i)，分别在p_j的[1,0,0]方向和[-1,0,0]方向确定两个点p_j+和p_j-，p_j+的邻域点集合N_r(p_j+)和p_j-的邻域点集合N_r(p_j+)之间的Hausdorff距离为： $d_H(H_r\left(p_{j+}\right),N_r\left(p_{j-}\right))=\max \left\{\underset{a\&Element;N_r\left(p_{j+}\right)}{\sup }\underset{b\&Element;N_r\left(p_{j-}\right)}{\inf }\right|a-b|,\underset{b\&Element;N_r\left(p_{j+}\right)}{\sup }\underset{a\&Element;N_r\left(p_{j-}\right)}{\inf }|a-b\left|\right\},$ 其中，i、r、j均为正整数，因此，p_j在X轴方向的偏导数为：f_x=d_H(N_r(p_j+),N_r(p_j-))，p_j在Y方向的偏导数计算方法与X方向相似，这里不再赘述。

为了消除噪声对计算偏导数的影响，为每个邻域点分配一个高斯权重。因此E矩阵可以写成如下形式：

$E=\frac{1}{\sqrt{2\mathrm{\&pi;\&sigma;}}}\left[\begin{array}{cc}\underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_x^2& \underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_x\&CenterDot;f_y\\ \underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_x\&CenterDot;f_y& \underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_y^2\end{array}\right]$ 其中，

(x_i,y_i,z_i)∈N_r(p_i)，σ的选择与点云噪声情况以及扫描的质量有关，其计算公式为：

因此可得p_i的Harris响应值Harris(x,y)=det(E)–k·(Trace(E))²。

多尺度空间中最优尺度下的候选特征点检测策略。为了自适应选择最优尺度下的特征点，算法首先选取一个K邻域点集合作为初始尺度空间大小，计算初始尺度因子σ₀。如果采样点p_i在当前尺度k下计算出的响应值h^k(p_i)不足以判定点p_i是否是特征点，那么尺度因子将更新到下一个尺度，继续计算下一个尺度p_i的响应值。随着尺度因子的增加，响应值变化最大的那个尺度称为最优的尺度。当迭代算法停止，将上一尺度及其响应值作为最优尺度和最优的响应值。算法的实现过程如下：

（1）令K=0，令特征点集合为空，其中，k为尺度因子索引；

（2）判断σ₀<σ_max，其中，σ₀为初始尺度因子，σ_max为最大尺度因子，如果是，进行进入第（3）步骤；

（3）对点云数据中任意一个点p_i∈P，提取pi的邻域点集合

$N_{{\mathrm{\&sigma;}}_k}^k\left(p_i\right)=\left\{p_j\right|{\left|\right|p_j-p_i\left|\right|}_2\&le;{\mathrm{\&sigma;}}_k\};$

（4）利用PCA计算pi的法向量

并建立局部坐标系；

（5）对任意的一个点

将pj的坐标平移、旋转变换转化为步骤（4）中建立的局部坐标系中的坐标；

（6）对任意的一个点

计算p_j在X轴方向和Y轴方向的偏导数，更新E矩阵；

（7）计算k尺度下点p_i的Harris特征值h^k(p_i)；

（8）判断h^k(p_i)>τ·h^k-1(p_i)是否成立，如果成立，令h^k(p_i)=h^k-1(p_i)，从集合P中删除p_i，将p_i添加到Q中；

（9）重复以上步骤，得到特征点集合Q。

基于MST算法的候选特征点优化策略。首先，设Q表示将候选特征点按照Harris响应值按降序排列得到的特征点集合，设q_i∈Q，

是q_i的K邻域N_k(q_i)中的点。定义

表示q_i和

之间边的权值，那么

的计算公式为： $\cos t(q_1,q_i^k)=\frac{1}{{(\mathrm{harris}\left(q_i\right)\&CenterDot;\mathrm{harris}\left(q_i^k\right))}^2}+\mathrm{\&gamma;}\&CenterDot;{\left|\right|q_i-q_i^k\left|\right|}_2,$ 式中：γ为常系数；设q表示候选特征点集合Q中Harris响应值最大的候选特征点；key[q_i]表示所有与q_i相连的边的最小权值，约定如果不存在这样的边，key[q_i]=∞；设π[q_i]表示q_i的父节点，那么特征点选择优化算法如下：

（1）初始化；

（2）判断Q是否为空，，如果不为空，则根据以下步骤寻找Q中harris(q_i)最小的候选特征点q_i；

（3）如果

∈Q并且

则将qi加入π中，即π=π∪qi，其中，qi∈Q，为所有与qi相连的边的最小权值，π为特征点选择优化后的特征点集合。

ST-RIFT特征描述算法设计，主要涉及的技术有：

1、在极坐标系中，对特征点所在切平面进行网格化分割。附图1显示了邻域点在特征点切平面的投影和投票过程。其中红色点代表N×M个网格节点，黑色点表示邻域点的投影点。

2、统计特征点邻域点在特征点切平面的投影长度信息，然后利用双线性插值将投影长度向特征统计矩阵的对应位置进行投票。定义一个N×M的实数矩阵S，划分的网格在半径方向将半径划分为M等份，每等份长度用binR表示，则binR=r(k)/M，在每一个同心圆上，按逆时针顺序将圆周角划分为M等份，每等份的角度用binθ表示，则binθ=2π/N。记f_i ^j在f_i切平面上的投影点为f_i ^j^

投票过程原理如下：

i.计算f_i ^j^与坐标原点的连线从x方向按逆时针旋转到当前位置的旋转角度，记为thetaR。判断是否成立，若成立，令

$\mathrm{thetaR}={\cos }^{-1}\left(\frac{{f_i}^{j^}\&CenterDot;u}{{\left|\right|{f_i}^{j^}\left|\right|}_2\&CenterDot;{\left|\right|u\left|\right|}_2}\right),$ 否则， $\mathrm{thetaR}=2\mathrm{\&pi;}-{\cos }^{-1}\left(\frac{{f_i}^{j^}\&CenterDot;u}{{\left|\right|{f_i}^{j^}\left|\right|}_2\&CenterDot;{\left|\right|u\left|\right|}_2}\right).$

ii.计算行索引h和列索引l。

iii.计算行方向的投票权重a和列方向的投票权重b。

$a=\frac{\mathrm{thetaR}-h\&CenterDot;\mathrm{bin\&theta;}}{\mathrm{bin\&theta;}},b=\frac{{\left|\right|{f_i}^{j^}\left|\right|}_2-j\&CenterDot;\mathrm{binR}}{\mathrm{binR}}$

iv.利用双线性差值和当前投影点周围的四个网格节点的投票权重计算当前投影点的位置。图2表示了投票权重的计算过程。如图2所示，可以得出：

$S(h,l)=S(h,l)+(1-a)(1-b)\&CenterDot;{\mathrm{proj}}_i^j;$

$S(h,l+1)=S(h,l+1)+(1-a)b\&CenterDot;{\mathrm{proj}}_i^j;$

$S(h+1,l)=S(h+1,l)+a(1-b)\&CenterDot;{\mathrm{proj}}_i^j;$

$S(h+1,l+1)=S(h+1,l+1)+\mathrm{ab}\&CenterDot;{\mathrm{proj}}_i^j$

式中：

——点f_i ^j的法向量

在点f_i所在切平面的投影。

v.再定义一个N×M维的实数矩阵Count，Count中的每一个元素Count(h,l)存储向S(h,l)投票的投影点的个数。同样的，为了准确计算投影点个数，利用双线性插值对Count矩阵赋值（即计算投影点周围的四个网格节点形成的区域内的投影点个数）。计算公式为：

Count(h,l)=Count(h,l)+(1-a)(1-b)·1；

Count(h,l+1)=Count(h,l+1)+(1-a)b·1；

Count(h+1,l)=Count(h+1,l)+a(1-b)·1；

Count(h+1,l+1)=Count(h+1,l+1)+ab·1。

vi.当所有投影点投票结束后，对S矩阵中的元素进行平均，即：

S(h,l)=S(h,l)/Count(h,l)

对特征统计矩阵的行和列进行DCT变换和DFT变换。首先对特征统计矩阵S的行向量进行DCT变换，用DCT变换后的行向量替换原始S矩阵中的行向量，得到新的S矩阵，然后对S矩阵的列向量进行DFT变换，DFT变换后的结果是复数，将变换后复数的模作为新的列向量元素。

以上所述仅为本发明的一种实施方式，不是全部或唯一的实施方式，本领域普通技术人员通过阅读本发明说明书而对本发明技术方案采取的任何等效的变换，均为本发明的权利要求所涵盖。

Claims (10)

1.一种非结构化点云特征点检测方法，其特征在于：包括以下步骤：

（1）利用采样点不同尺度空间的邻域信息计算采样点在不同尺度空间中的Harris响应值；

（2）通过迭代算法选取最优尺度空间下的Harris响应值作为采样点的Harris响应值，得到特征点集合Q；

（3）将Harris响应值在尺度空间邻域和几何邻域都具有极大性的点作为候选特征点，最后利用候选特征点选择优化策略提取最终的特征点。

2.根据权利要求1所述的一种非结构化点云特征点检测方法，其特征在于：所述步骤（1）中，采样点在不同尺度空间中的Harris响应值根据以下方法计算：（1.1）设N_r(p_i)是p_i的邻域点集合，p_j∈N_r(p_i)，分别在p_j的[1,0,0]方向和[-1,0,0]方向确定两个点p_j+和p_j-，计算p_j+的邻域点集合N_r(p_j+)和p_j-的邻域点集合N_r(p_j+)之间的Hausdorff距离；其中，i、r、j均为正整数；

（1.2）根据步骤（1.1）的Hausdorff距离计算p_j在X轴方向的偏导数f_x；

（1.3）计算p_j在Y方向的偏导数f_y；

（1.4）根据p_j在X轴方向和Y方向的偏导数f_x和f_y计算每个邻域点的高斯权重，即E矩阵；

（1.5）根据E矩阵计算p_j的Harris响应值，计算公式为：Harris(x,y)=det(E)–k·(Trace(E))²，k为经验系数。

3.根据权利要求2所述的一种非结构化点云特征点检测方法，其特征在于：所述步骤（1.1）的Hausdorff距离根据以下公式计算：

$d_H(H_r\left(p_{j+}\right),N_r\left(p_{j-}\right))=\max \left\{\underset{a\&Element;N_r\left(p_{j+}\right)}{\sup }\underset{b\&Element;N_r\left(p_{j-}\right)}{\inf }\right|a-b|,\underset{b\&Element;N_r\left(p_{j+}\right)}{\sup }\underset{a\&Element;N_r\left(p_{j-}\right)}{\inf }|a-b\left|\right\}.$

4.根据权利要求2所述的一种非结构化点云特征点检测方法，其特征在于：所述（1.4）的高斯权重E矩阵根据以下公式计算：

$E=\frac{1}{\sqrt{2\mathrm{\&pi;\&sigma;}}}\left[\begin{array}{cc}\underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_x^2& \underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_x\&CenterDot;f_y\\ \underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_x\&CenterDot;f_y& \underset{x_i,y_i,z_i}{\mathrm{\&Sigma;}}{e}^{-(x_i^2+y_i^2+z_i^2)/2\&CenterDot;{\mathrm{\&sigma;}}^2}\&CenterDot;f_y^2\end{array}\right],$

其中，(x_i,y_i,z_i)∈N_r(p_i)，σ为尺度因子，第k个尺度因子根据以下公式计算：

${\mathrm{\&sigma;}}_k=\frac{1}{\left|N_r^k\left(p_i\right)\right|}\underset{j=1}{\overset{\left|N_r^k\left(p_i\right)\right|}{\mathrm{\&Sigma;}}}{\left|\right|p_i-p_i^j\left|\right|}_2.$

5.根据权利要求1所述的一种非结构化点云特征点检测方法，其特征在于：所述步骤（2）的方法为：首先选取一个K邻域点集合作为初始尺度空间大小，计算初始尺度因子σ₀；如果采样点p_i在当前尺度k下计算出的响应值h^k(p_i)不足以判定点p_i是否是特征点，那么尺度因子将更新到下一个尺度，继续计算下一个尺度p_i的响应值；随着尺度因子的增加，响应值变化最大的那个尺度称为最优的尺度；当迭代算法停止，将上一尺度及其响应值作为最优尺度和最优的响应值；具体方法为：

（2.1）令k=0，特征点集合为空，其中，k为尺度因子索引；

（2.2）判断σ₀<σ_max是否成立，其中，σ₀为初始尺度因子，σ_max为最大尺度因子，如果成立，则提取p_i的邻域点集合

其中， $N_{{\mathrm{\&sigma;}}_k}^k\left(p_i\right)=\left\{p_j\right|{\left|\right|p_j-p_i\left|\right|}_2\&le;{\mathrm{\&sigma;}}_k\},$ 其中，p_i为点云数据集P中的任意一点；

（2.3）计算p_i的法向量

并建立局部坐标系；

（2.4）将p_j的坐标平移、旋转变换转化为步骤（2.3）中建立的局部坐标系中的坐标，其中， $p_j\&Element;N_{{\mathrm{\&sigma;}}_k}^k\left(p_i\right);$

（2.5）计算p_j在X轴方向和Y轴方向的偏导数，更新E矩阵；

（2.6）计算k尺度下点p_i的Harris特征值h^k(p_i)；

（2.7）判断h^k(p_i)>τ·h^k-1(p_i)是否成立，如果成立，令h^k(p_i)=h^k-1(p_i)，将p_i，添加到Q中；

（2.8）重复以上步骤，得到特征点集合Q。

6.根据权利要求1所述的一种非结构化点云特征点检测方法，其特征在于：所述步骤（3）中，候选特征点的优化策略方法为：

（3.1）初始化后，判断特征点集合Q是否为空，如果不为空，则进行步骤（3.2）；（3.2）如果

∈Q并且

则将q_i加入π中，即π=π∪q_i，其中，q_i∈Q，

为所有与q_i相连的边的最小权值，π为特征点选择优化后的特征点集合。

7.一种基于权利要求1所述的非结构化点云特征点检测方法对非结构化点云进行提取的方法，其特征在于：在极坐标系下对权利要求1得到的特征点所在切平面进行网格化分割，然后将特征点的邻域点投影到该切平面上，将每个网格中的投影点对应的投影长度向网格的四个顶点投票生成特征信息统计矩阵，然后分别对特征统计的行向量和列向量进行DCT变换和DFT变换，将由变换后得到的矩阵的左上角元素构成的向量作为该特征点的特征描述向量，所述左上角为A×A，其中，1<A<10。

8.根据权利要求7所述的方法，其特征在于：利用双线性插值将投影长度向特征统计矩阵的对应位置进行投票的方法为：定义一个N×M的实数矩阵S，划分的网格在半径方向将半径划分为M等份，每等份长度用binR表示，则binR=r(k)/M，在每一个同心圆上，按逆时针顺序将圆周角划分为M等份，每等份的角度用binθ表示，则binθ=2π/N；记第j个特征点f_i ^j在f_i切平面上的投影点为f_i ^j^, ${f_i}^{j^}=(x_i^{j^},y_i^{j^}),$ 投票方法为：

（A）计算f_i ^j^与坐标原点的连线从x方向按逆时针旋转到当前位置的旋转角度，记为thetaR；判断

是否成立，若成立，令

否则，令 $\mathrm{thetaR}=2\mathrm{\&pi;}-{\cos }^{-1}\left(\frac{{f_i}^{j^}\&CenterDot;u}{{\left|\right|{f_i}^{j^}\left|\right|}_2\&CenterDot;{\left|\right|u\left|\right|}_2}\right);$

（B）计算行索引h和列索引l，其中，

（C）计算行方向的投票权重a和列方向的投票权重b，其中， $a=\frac{\mathrm{thetaR}-h\&CenterDot;\mathrm{bin\&theta;}}{\mathrm{bin\&theta;}},b=\frac{{\left|\right|{f_i}^{j^}\left|\right|}_2-j\&CenterDot;\mathrm{binR}}{\mathrm{binR}};$

（D）利用双线性差值和当前投影点周围的四个网格节点的投票权重计算当前投影点的位置；

（E）根据步骤（D）计算投影点周围的四个网格节点形成的区域内的投影点个数；

（F）当所有投影点投票结束后，对S矩阵中的元素进行平均，得到特征统计矩阵。

9.根据权利要求8所述的方法，其特征在于：步骤（D）根据以下公式计算：

$S(h,l)=S(h,l)+(1-a)(1-b)\&CenterDot;{\mathrm{proj}}_i^j;$

$S(h,l+1)=S(h,l+1)+(1-a)b\&CenterDot;{\mathrm{proj}}_i^j;$

$S(h+1,l)=S(h+1,l)+a(1-b)\&CenterDot;{\mathrm{proj}}_i^j;$

$S(h+1,l+1)=S(h+1,l+1)+\mathrm{ab}\&CenterDot;{\mathrm{proj}}_i^j;$

式中：

为点f_i ^j的法向量在点f_i所在切平面的投影。

10.根据权利要求7所述的方法，其特征在于：所述DCT变换和DFT变换的方法为：首先对特征统计矩阵S的行向量进行DCT变换，用DCT变换后的行向量替换原始S矩阵中的行向量，得到新的S矩阵，然后对S矩阵的列向量进行DFT变换，DFT变换后的结果是复数，将变换后复数的模作为新的列向量元素。

Images