CN108180918B

CN108180918B - 一种点云测地路径正向跟踪生成方法及装置

Info

Publication number: CN108180918B
Application number: CN201711228825.9A
Authority: CN
Inventors: 李自胜; 肖晓萍
Original assignee: Southwest University of Science and Technology
Current assignee: Southwest University of Science and Technology
Priority date: 2017-11-29
Filing date: 2017-11-29
Publication date: 2021-04-30
Anticipated expiration: 2037-11-29
Also published as: CN108180918A

Abstract

本发明涉及点云数据处理领域。本发明选择测地路径起点p_s(x_s,y_s,z_s)和终点p_e(x_e,y_e,z_e)之间的主行进方向，以主行进方向来确定网格化区域；提取网格化区域内所有数据点坐标值，构成三个坐标分量数组，排序数组并剔除相同坐标分量值，找出网格化区域中所有点云数据点坐标的最大值与最小值；分别计算数据点在三个轴向的步长hx_i，hy_j和hz_k，非均匀网格化两点p_s和p_e之间的数据点，然后执行步骤3；根据步长hx_i，hy_j和hz_k，采用UNCDFMM计算波前网格各单元格的达到时间值；计算出各单元格的到达时间后，筛选出满足近似测地线性质的单元格，将这些单元格的顶点顺序连接构成测地路径。

Description

一种点云测地路径正向跟踪生成方法及装置

技术领域

本发明涉及点云数据处理领域，尤其是一种点云测地路径正向跟踪生成方法及装置。

背景技术

在包含待测两点内的地球面上测得的两点之间的最短线。测地线又称大底线或短程线，可以定义为空间中两点最短或最长路径。

航线：生活中飞机或轮船的航行都是测地线。

圆柱上的测地线：将圆柱面剪开铺平，给出了圆柱面到平面的一个等距变换，平面的测地线是直线。因此圆柱面上的测地线，将平面卷成圆柱面时，由平面上直线变过来的曲线，容易发现它们是直线(圆柱面的子母线)，平行圆或圆柱螺线。自然界的攀缘植物，沿螺线生长，就是测地线的一个有趣例子。

球面上的测地线：由于大圆的主法线就是球面的法线，由此推出大圆或其一部分就是球面上的测地线。在球面上，连接两点的大圆弧(测地线)有两条，这两条大圆弧一长一短，短的是最短线，也就是短程线。

物理学中的测地线：假设一个质点在曲面上自由运动，如无外力作用，则质点在曲面上的运动轨迹是一条测地线。再如一条无重量的弹性细线在一光滑曲面上自由移动，当这条细线受曲面上两点问的某张力作用处于平衡状态时，则它取测地线的形状。

测地线在工程技术中的应用：复杂曲面的外板展开，这种板金技术在飞机机身、汽车外壳、轮船船体、涡轮叶片、薄壳屋顶、机耕梨面等外形设计和制造。

测地线在制衣行业中的应用：帐篷脊线、鞋的足背线、衣服的腰线等。

1.2曲面测地线。

在微分几何中，如：

图1光滑曲面S的参数形式为r(u,v)，P₀＝P(u₀,v₀)为曲面S上与(u₀,v₀)对应的点，(u(s),v(s))是参数域D上过点(u₀,v₀)的一条曲线，且u₀＝u(0)，v₀＝v(0)，设s是曲线r(s)的弧长参数，则r(s)＝r(u(s),v(s))是过曲面S上过点P₀的一条曲线，且P₀对应参数s＝0。设T为曲线在P₀点的单位切向量，设n为曲面在P₀点的单位法向量，设β为曲线在P₀点的单位主法向量，设κ_β为曲线在P₀点的曲率向量。曲线r(s)在p₀点的单位切向量为：

r(s)在p₀点的曲率向量为：

曲线r(s)的弯曲由两部分产生：曲线随曲面的弯曲产生的法曲率k_n和曲线自身弯曲产生的测地曲率k_g。则κ_β可分解为向曲面法向量n的投影和过P₀切平面的投影(副法线方向)，也即是：

k_n＝k_β·n，k_g＝k_β·(n×T)；

测地线是曲面上测地曲率k_g＝0的曲线，即曲线的主法向量与曲面的法向量平行，同时也是曲面上局部距离最短的曲线，又称最短程线。

光滑曲面上可通过求解测地微分方程得到一条经过给定切方向的精确测地线，但对于网格模型，无法得到模型的参数表达式，因此只能采用近似算法来计算得到一组离散点，称为近似测地路径，测地路径上点的依次连线称为近似测地线。赵对离散网格模型上的测地线做了综述，指出离散网格模型上的近似测地线由于不能保证同时满足光滑参数曲面上精确测地线的条件，离散网格模型上的近似测地线就分为最直测地线和最短测地线。由于最短测地线满足三角不等式，因而是度量，因此离散模型的测地线研究中多研究最短测地线。

与离散网格模型不同，点云模型由散乱数据点构成，不具有模型表达式，也不具有网格结构，无法计算其最直测地线，因此对点云模型测地线的研究主要是利用测地线局部最短性质，近似计算其最短测地线。Memol定义球心在给定点的一组球来生成点云的偏置带(offset band)，通过fast marching method在偏置带构成的空间网格来计算近似测地线，其精度决定于点云采样密度、速度决定于网格数，对噪声的鲁棒性取决于球的半径，由于采用一组球面构成偏置带，无法区分角点、边或边界等特征。Klein通过在点云上构建几何近似图(Delaunay图和SIG(spheres-of-influence graph))来计算近似测地线，Hofer在点云上构建MSL曲面，将能量泛函作用于MSL曲面，并在曲面上施加约束使其最小化来计算测地线，该方法对噪声敏感，计算的测地线精度依赖于构建的MSL曲面。与Hofer类似，Ruggeri也运用能量泛函来计算测地线，能量函数由相连两点间的平方距离和L(P)、点贴近曲面片的惩罚项D_s(P)和两点间连线贴近曲面片的惩罚项U_s(P)构成，并在D_s(P)和U_s(P)上施加可调因子α和β来控制收敛速度和近似精度，初始路径利用Dijkstra算法计算。与Hofer不同的是后者无需将约束施加于构建的MSL曲面上，而是在能量泛函中加入约束，避免对点云模型曲面的依赖。肖采用MLS对点云进行采样和去噪处理，利用窄带Level Set生成一条虚拟路径，最后通过计算点云上点到虚拟路径最小的点来得到测地线。杜用Dijkstra算法找到两点间的初始测地路径，再用抛物面拟合点云邻近域，估算抛物面的法线，最小化平方距离度量函数来计算测地线。KEENAN采用热核方法计算点云测地线，但在热核离散化计算过程中利用了点云Laplace-Betrami算子矩阵，该矩阵计算需要利用点云的法线，Yu首先对点云进行空间单元格划分，然后利用连续Dijkstra算法计算一条近似测地路径，最后利用测地曲率流对测地线进行校正，该方法比较有效，但在校正过程中需要用到点云法线。

现有计算点云测地线的方法，主要有两类。(1)依赖点云法线的方法：在测地线计算过程中需要用到点云法线的方法都存在先天缺陷，因为点云法线是不知道的，需要对目标点云的法线进行估算来得到，法线估算精度不高，法线方向不能全局一致，出现混乱，因而测地线计算精度差。(2)单元格划分的方法：在对实物模型进行采样过程中，各区域几何形状不同，采样密度不同。较平坦的区域通常用较大的间隔(较小的密度)，弯曲部分用较小的间隔(较大的密度)，由于采样的各向异性，现有网格均匀化方法对点云的进单元格划分，由大到小不断细分，直到单元格仅包含一个点云数据点，这种划分是均匀的(X,Y,Z三个方向的间距相等)，导致点云数据点往往位于单元格的内部，以单元格的位置来代替数据点的位置坐标，最后生成的测地线几乎不经过点云的任何数据点(点云数据不在测地线上)，未能忠实于点云数据，因此一种精度不高的计算方法。

发明内容

本发明所要解决的技术问题是：针对现有技术存在的问题，提供一种点云测地路径正向跟踪生成方法。根据点云的密度和间距实现非均匀网格划分，以使点云数据点全部位于网格的顶点上，保证了后续关于单元格的各项计算忠实于点云数据。

本发明采用的技术方案如下：

一种点云测地路径正向跟踪生成方法包括：

步骤2：提取网格化区域内所有数据点坐标值，构成三个坐标分量数组，排序数组并剔除相同坐标分量值，找出网格化区域中所有点云数据点坐标的最大值与最小值，记为：(x_max，y_max，z_max)，(x_min，y_min，z_min)；分别计算数据点在三个轴向的步长hx_i，hy_j和hz_k，非均匀化起点p_s(x_s,y_s,z_s)和终点p_e(x_e,y_e,z_e)之间的数据点，然后执行步骤3；

步骤3：根据步长hx_i，hy_j和hz_k，采用UNCDFMM(单向非均匀紧致差分快速行进法)计算波前网格各单元格的达到时间值；

步骤4：计算出各单元格的到达时间值后，筛选出满足近似测地线性质的单元格，将这些单元格的顶点顺序连接构成测地路径。

进一步的，所述步骤1具体过程是：

步骤12：取x轴为主行进方向，如果x_s<x_e，选择点云数据中x坐标满足x_s≤x≤x_e的点来当做网格化区域，x<x_s或x>x_e在测地路径计算中无意义，如果测地线路径过这类数据点，测地路径将会折回来，不满足测地路径最短性质；如果x_s>x_e，选择点云数据中x坐标满足x_s≥x≥x_e的点当做网格化区域；

步骤13：取y轴为主行进方向，如果y_s<y_e，选择点云数据中y坐标满足y_s≤y≤y_e的点来当做网格化区域；如果y_s>y_e，选择点云数据中y坐标满足y_s≥y≥y_e的点来当做网格化区域；

步骤14：取z轴为主行进方向，如果z_s<z_e，选择点云数据中z坐标满足z_s≤z≤z_e的点当做网格化区域；如果z_s>z_e，选择点云数据中z坐标满足z_s≥z≥z_e的点当做网格化区域。

进一步的，所述步骤2具体过程：

步骤21：找出所有点云数据点坐标的最大值与最小值，记为：(x_max，y_max，z_max)，(x_min，y_min，z_min)；

步骤22：假设网格区域内数据点数量为N，原数据点的坐标为p_m(x_m,y_m,z_m)，1≤m≤N，取出全部数据点的三个坐标分量，分别构成三个数组，排序数组并剔除同一数组内相同的坐标值，得到三个新的坐标分量数组

和

数组元素个数分别为I，J，K；由于在排序过程中坐标值相同的坐标分量只保留一个，因此坐标分量排序后，其坐标分量值数量小于等于N，则I≤N，J≤N，K≤N；新坐标分量数组在空间上构成三维非均匀化网格；

步骤23：分别计算新的坐标分量数组

和

形成的数据点在三个轴向的步长

令

此时2≤i≤I-1；2≤j≤J-1；2≤k≤K-1；

步骤24：在新的坐标分量数组

和

的最后一个分量末尾增加ε＝1.0×e^-3，如果x_I<x_max+ε，则

否则

如果y_J<y_max+ε，则

否则，

如果z_K<z_max+ε，则

否则，

以使步长

和

的数目与各轴向方向的单元格数目相同，方便UNCDFMM计算；x_max是x轴分量数组

的最大值；y_max是y轴分量数组

的最大值；z_max是z轴分量数组

的最大值。

进一步的，所述步骤3具体过程：

步骤31：采用紧致差分来计算新的坐标分量数组

和

构成的网格点数值，构造单向三点二阶紧致差分计算网格点数值；X轴、Y轴和Z轴方向上单向三点二阶紧致差分格式相同，以X轴向为例，令T_i为分量数组

中x_i处的数值，令

为分量数组

中x_i处的单向前向二阶紧致差分值，令

为分量数组

中x_i处的单向后向二阶紧致差分值，则表达式分别为：

和

分别为前向二阶紧致差分表达式(3)的系数，当i＝I-1时，前向二阶紧致差分不存在，用前向一阶差分

近似前向二阶紧致差分值，当i＝I时前向二阶紧致差分不存在。

和

分别为后向二阶紧致差分表达式(4)的系数，当i＝2时，后向二阶紧致差分不存在，用后向一阶差分

近似后向二阶紧致差分值，当i＝1时后向二阶紧致差分不存在。

由于X轴、Y轴和Z轴方向上单向二阶紧致差分格式相同，以X轴向为例，单向二阶紧致差分的系数计算如下：对公式(3)中T_i+1和T_i+2进行泰勒级数展开，得到公式(5)和(6)；对公式(4)中T_i-1和T_i-2进行泰勒级数展开，得到公式(7)和(8)。

其中T_i'、T_i”、T_i”'和

分别表示坐标分量数组

上x_i处数值T_i对应的一阶、二阶、三阶和四阶导数；

表示

同理可知

和

意义。

将式(5)和式(6)代入式(3)，合并相同项，式(3)中

即是T_i”，得到:

比较式(9)中T_i、T_i'和T_i”系数，构造系数方程组并求解，可得式(3)的系数

和

为：

同理，将式(7)和式(8)代入式(4)，合并相同项，式(4)中

即是T_i”，得到:

比较式(11)中T_i、T_i'和T_i”系数，构造系数方程组并求解，可得式(4)的系数

和

为：

因此在X轴向上，(3)式和(4)式中单向前向二阶紧致差分

和单向后向二阶紧致差分

分别为：

同理可计算Y轴向的单向二阶紧致差分值：

和Z轴向的单向二阶紧致差分值：

X轴、Y轴和Z轴坐标分量数组构成空间网格，X轴向第i分量对应数值T_i，Y轴向第j分量对应数值T_j，Z轴向第k分量对应数值T_k，[i,j,k]在三维网格中对应数值T_ijk的位置索引；

由Eikonal方程数字解公式：

其中D^-x(T_ijk)、D^+x(T_ijk)、D^-y(T_ijk)、D^+y(T_ijk)、D^-z(T_ijk)和D^+z(T_ijk)分别表示T_ijk在X轴向的后向一阶差分、前向一阶差分、在Y轴向的后向一阶差分、前向一阶差分、在Z轴向的后向一阶差分和前向一阶差分，F_ijk为单元格[i,j,k]处的速度，max(D^-x(T_ijk),0)²表示取D^-x(T_ijk)与0两者的最大值再平方，和min(-D^+x(T_ijk),0)²表示取-D^+x(T_ijk)与0两者的最小值再平方，同理可知max(D^-y(T_ijk),0)²、min(-D^+y(T_ijk),0)²、max(D^-z(T_ijk),0)²、min(-D^+z(T_ijk),0)²表达的意义。

用二阶差分形式

和

代替Eikonal公式中一阶差分形式D^-和D⁺，并分别赋予X、Y和Z轴向

和

则得到在单元格[i,j,k]处方程：

(19)式是在X、Y和Z三个方向的方程；更一般地，在26领域方案中有26个方向，此时

和

分别表示以[i,j,k]为公共端点两共线射线方向上的单向前向二阶紧致差分和单向后向二阶紧致差分。

步骤32：当前单元格邻近单元格分为波经过的单元格Frn和未经过的单元格O两类，通过生成Frn的幂集Frn^p，利用幂集的每个子集来计算当前单元格的数值，并取这些数值的最小值作为当前单元格的达到时间值，将每个未计算出到达时间值的单元格作为当前单元格，采用该步骤进行计算，最终得到整个网格各单元格的到达时间值；Frn代表Frozen状态的单元格；O代表Narrow Band状态或Open状态；幂集Frn^p包括至少一个单元格；Frozen状态表示波经过单元格，时间已经计算出，不会再改变，Narrow Band状态表示波经过单元格，时间值已经计算出，需要改变更新，Open状态表示波还未经过，时间值未知的状态。

进一步的，所述步骤32单元格的到达时间值计算过程具体包括：

步骤321：定义一个暂存单元格数值数组U，令当前单元格为C_ijk，依次从Frn^p中选择一个集合Set；

步骤322：如果集合Set的成员个数为1，即集合仅包含一个单元格C_lmn(l，m和n为单元格空间位置索引)，则

将T_ijk加入数组U中。其中|C_ijk-C_lmn|表示单元格C_ijk和C_lmn之间的距离，F_ijk表示波在索引为[i,j,k]处单元格的传播速度；T_lmn表示索引为[l,m,n]单元格的时间值。

步骤323：如果集合Set的成员个数为2个，即集合包含两个单元格

和

如果

与C_ijk在同一方向上，则将

和

的值带入(19)式，构造一元二次方程求解单元格的到达时间值T_ijk，加入数组U中；

三个单元格在同一个方向上判断规则：如果C_ijk、

和

的下标依顺序

或

满足如下顺序关系(以顺序关系(1)为例，[i,j,k]为C_ijk的下标，[i,j,k+1]为

下标，[i,j,k+2]为

下标；或者[i,j,k+1]为

下标，[i,j,k+2]为

下标)：

{

(1)[i,j,k]→[i,j,k+1]→[i,j,k+2]，

(2)[i,j,k]→[i,j,k-1]→[i,j,k-2]，

(3)[i,j,k]→[i,j+1,k]→[i,j+2,k]，

(4)[i,j,k]→[i,j+1,k+1]→[i,j+2,k+2]，

(5)[i,j,k]→[i,j+1,k-1]→[i,j+2,k-2]，

(6)[i,j,k]→[i,j-1,k]→[i,j-2,k]，

(7)[i,j,k]→[i,j-1,k+1]→[i,j-2,k+2]，

(8)[i,j,k]→[i,j-1,k-1]→[i,j-2,k-2]，

(9)[i,j,k]→[i+1,j,k]→[i+2,j,k]，

(10)[i,j,k]→[i+1,j,k+1]→[i+2,j,k+2]，

(11)[i,j,k]→[i+1,j,k-1]→[i+2,j,k-2]，

(12)[i,j,k]→[i+1,j+1,k]→[i+2,j+2,k]，

(13)[i,j,k]→[i+1,j+1,k+1]→[i+2,j+2,k+2]，

(14)[i,j,k]→[i+1,j+1,k-1]→[i+2,j+2,k-2]

(15)[i,j,k]→[i+1,j-1,k]→[i+2,j-2,k]，

(16)[i,j,k]→[i+1,j-1,k+1]→[i+2,j-2,k+2]，

(17)[i,j,k]→[i+1,j-1,k-1]→[i+2,j-2,k-2]，

(18)[i,j,k]→[i-1,j,k]→[i-2,j,k]，

(19)[i,j,k]→[i-1,j,k+1]→[i-2,j,k+2]，

(20)[i,j,k]→[i-1,j,k-1]→[i-2,j,k-2]，

(21)[i,j,k]→[i-1,j+1,k]→[i-2,j+2,k]，

(22)[i,j,k]→[i-1,j+1,k+1]→[i-2,j+2,k+2]，

(23)[i,j,k]→[i-1,j+1,k-1]→[i-2,j+2,k-2]，

(24)[i,j,k]→[i-1,j-1,k]→[i-2,j-2,k]，

(25)[i,j,k]→[i-1,j-1,k+1]→[i-2,j-2,k+2]，

(26)[i,j,k]→[i-1,j-1,k-1]→[i-2,j-2,k-2]

}

步骤324：如果集合Set的成员个数大于2个，令Set＝{C₁,C₂,…C_Q}，分别将C_q(1≤q≤Q)在C_ijk处一阶渐进展开，可得：

其中sgn表示正负符号，“+”表示正号，“-”表示负号，正负符号根据单元格位置三个索引之差的乘积的正负符号来确定，例如：[2,3,5]-[1,4,6]＝[1,-1,-1]，sgn＝(+)×(-)×(-)＝,+

由可由下式计算：

其中

M⁺＝(M^tM)^-1M^t

计算出

就得到单元格C_ijk邻近单元格集合Set的到达时间值T_ijk，将其加入到数组U中；M^t表示矩阵M转置；(M^tM)^-1表示矩阵M^tM的逆矩阵。

步骤325：幂集Frn^p的每个集合Set都计算完毕，从数组U中选择数值最小且值大于Frn中每个单元格数值作为当前单元格C_ijk的达到时间值。

进一步的，所述步骤4具体包括：

步骤41：设S是E³曲面，r＝r(u,v)是S关于空间参数u和v的曲面表示，假设弧长参数为s；将空间参数u和v分别表示为弧长参数为u(s)和v(s)，参数曲线Γ:r(s)＝r(u(s),v(s))，简记为r＝r(s)，是曲面S上的一条弧长参数曲线；沿曲线Γ取曲面的正交标架

其中切向

法向

且

正定向，即向量混合积

是曲线的切向，

为测地曲率矢量，

为法向；

步骤42：曲面S上的弧长参数曲线r＝r(s)测地曲率

当k_g＝0时，曲线r为测地线；将测地线起点p_s与终点p_e的连线近似为曲线的切线

以p_s为当前点，根据

的正定向条件，可以计算出p_s的一组

方向和对应的

为p_s确定出了

和

就可以找到p_s的一个或多个可能行进方向，得到一组测地线从p_s可能经过的一组网格点{q_j}_j＝1…m，然后依次将{q_j}_j＝1…m中的点作为当前点，重复上述步骤，直到找到终点p_e为止；上述测地线的经过的节点可用树状结构来描述，树的根节点为p_s，树的叶节点均为p_e，因此得到测地线从p_s到p_e的多条路径，深度遍历树并计算其长度，选择长度最小的路径作为最终的测地路径；其中测地曲率k_g采用与当前单元格邻近单元格时间值的二阶差分近似，选择测地曲率绝对值最小的点所在的方向作为

尽管网格划分是非均匀的，但单元格之间正交，对应的

方向必在与

正交方向上。

一种点云测地路径正向跟踪生成装置包括：

步长计算模块，用于提取网格化区域内所有数据点坐标值，构成三个坐标分量数组，排序数组并剔除相同坐标分量值，找出网格化区域中所有点云数据点坐标的最大值与最小值，记为：(x_max，y_max，z_max)，(x_min，y_min，z_min)；分别计算数据点在三个轴向的步长hx_i，hy_j和hz_k，非均匀化起点p_s(x_s,y_s,z_s)和终点p_e(x_e,y_e,z_e)之间的数据点；

达到时间值计算模块，用于根据所述步长hx_i，hy_j和hz_k，采用UNCDFMM计算波前网格各单元格的达到时间值；

测地路径形成模块，用于计算出各单元格的到达时间值后，筛选出满足近似测地线性质的单元格，将这些单元格的顶点顺序连接构成测地路径。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

(1)点云非均匀网格化：本专利根据点云的密度和间距实现非均匀网格划分，以使点云数据点全部位于网格的顶点上，保证了后续关于单元格的各项计算忠实于点云数据。

(2)UNCDFMM(Unilateral Nonuniform Compact Difference FMM，单向非均匀紧致差分快速行进法)：OFMM适用于在均匀划分网格上行进，用的一阶差分(一阶前向，后向或中心差分)来为当前单元格的每个方向选择最小值的单元格构造一元二次方程，HAFMM用于均匀划分网格，采用二阶差分或二阶致紧差分，紧致差分用的固定系数3，-4和1，来为当前单元格的每个方向选择最小值的单元格构造一元二次方程。本专利针对非均匀化网格，采用前向(后向)三点非均匀紧致差分快速行进法(统称：UNCDFMM)来为当前单元格的每个方向选择最小值的单元格构造一元二次方程，计算精度提高。

(3)邻近单元格分类与单元格求解：在计算单元格的到达时间值过程中，OFMM和HAFMM从当前单元格的邻近三个坐标轴(X,Y,Z)方向6个单元格(6-邻近单元格)上选择各一个时间最小的数值来构造一元二次方程，求解出当前单元格的时间值。上述行进方法对6-邻近单元格模式容易实现每个方向最小值选择，对于26-单元格模式，将有13个方向，按照6-邻近单元格模式很难准确选出13个最小值来构造一元二次方程，邻近单元格中存在“Frozen”,“Narrow Band”或“Open”状态的单元格，邻近最小值单元格筛选困难，求解过程容易出现混乱。本专利将当前单元格的邻近单元格(6-邻近单元格或26-邻近单元格模式)分为两类，“Frozen”状态和非“Frozen”状态，将“Frozen”状态的单元格构成集合，计算其幂集，依次在幂集的成员上计算当前单元格的时间到达值，选择这些值中最小者作为当前单元格的值。在每个幂集成员上计算单元格值过程中，将当前单元格用集合成员Tailor级数展开，构造线性方程组来求解。

(4)正向跟踪生成测地路径：现有基于FMM的方法生成测地路径都是采用反向跟踪法(梯度下降法)，基本思想是：令测地线起点值为0，从测地线终点开始，计算当前点在各个方向的梯度值，选择梯度值最大的方向作为行进方向，计算下一个点的位置坐标，取整该坐标值作为下一个单元格的索引位置，如此循环下去，直到遇到值为0的单元格停止跟踪，将所有点连接起来构成测地线。这种反向跟踪方法，理论上可行，但实际上网格点的数值通过FMM计算得到，FMM是一种数值近似计算，单元格的数值没有内在规律，采用梯度值计算下一个点的位置坐标为浮点数，取整后往往超出网格边界，或者网格点之间来回震荡，无法确定曲线的下一步走向，导致反向跟踪无法继续进行。本专利采用正向跟踪方法，避免了无法确定测地线走向的确定。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是曲线曲率、法曲率与测地曲率示意图。

图2是测地线主方向示意图。

图3a是均匀网格划分示意图。

图3b是非均匀网格划分示意图。

图4点云非均匀网格化示意图。

图5是邻近单元格划分为Frn和O两类示意图。

图6a是测地路径正向路径跟踪示意图一。

图6b是是测地路径正向路径跟踪示意图二。

图7a是球面模型(401点)测地路径。

图7b是圆环面模型(2500点)测地路径。

图7c是飞机操作杆手柄模型(5006点)测地路径。

图7d是吹风机出风口模型(9458点)测地路径。

图7e是钣金结构件模型(10086点)测地路径。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

1、本发明基础背景：

1.1Fast Marching Method(快速行进法)

Fast Marching Method(简称：FMM)通过求解Eikonal Equation边值问题,计算界面在网格法向扩展到达各网格点时间的数值方法。FMM广泛用于计算机图形学、图像处理等研究领域。Eikonal Equation如下：

其中F(x)≥0是波界面在网格位置x的速度函数，T(x)是界面到达网格位置x的时间函数。到达时间函数的梯度幅值反比于速度，即是：

波在起点处T＝0，由于F(x)≥0，波在沿法向方向演化过程中，除原点外任一点处的到达时间值均大于0，其演化界面在二维情况下为一平面曲线，在三维情况下为一曲面。其离散解在三维情况下为：

其中：

离散FMM在网格上采用迭代方法实现，每个单元格将被标记为三种类型之一：

Open:波还未传播到该单元格，时间T未知的；

Narrow Band:波界面下一步将会扩展的候选单元格，具有即将更新的时间T；

Frozen：波界面已经经过该网格，时间T是固定的。

对于OFMM(Ordinary Fast Marching Method)，由(26)选择同维网格最小值T₁、T₂和T₃，将其带入(24)即可得到关于T的一元二次方程，方程系数分别为a，b，c，求解出方程解

不放假设T₃>T₂>T₁(总可以通过排序使之成立)，根据迎风条件，其解为：

(1)当T>max(T₁，，T₂)T₃时，

(2)当T₂<T<T₃时，选择T₂和T₁重新生成一元二次方程，方程系数分别为a₁，b₁，c₁，此时

(3)当T₁<T<T₂时，

由于一阶差分近似计算结果误差较大，HAFMM(Higher Accuracy Fast MarchingMethod)，利用二阶差分来近似梯度计算网格点的数值，Rickett和

利用

分别来近似网格点的后向差分和前向差分，该方法中网格等距且节点的系数分别定为3，-4和1。同样方式构造Y方向和Z方向的前向和后向差分，由(27)式选择同维网格最小值T₁、T₂和T₃，与OFMM方法类似，通过一元二次方程的解，在满足迎风条件下得到网格点的值T。

在HAFMM中，计算一个单元格的值在一坐标轴方向上将涉及附近4个单元格，这些单元格可能处于“Open”、“Narrow Band”或“Frozen”的状态之一，当处于前两种状态下时，二阶差分近似计算将得到一个无效值，因此HAFMM需要满足以下条件：

(1)与当前单元格距离为2个单位的单元格必须处于“Frozen”状态，如x方向上的T_i-2,j,k和T_i+2,j,k应处于“Frozen”状态；

(2)与当前单元格距离为2个单位的单元格的数值必须小于等于距离为1个单位的单元格数值，如x方向上，T_i-2,j,k≤T_i-1,j,k与T_i+2,j,k≤T_i+1,j,k；

当HAFMM不满足以下条件时，用一阶差分近似梯度代替二阶差分近似。

1.2 Compact Difference(紧致差分)

紧致差分通常用Hermite公式表示

其实质就是f(x)邻近节点函数值、一阶导数值及二阶导数值组合而成的表达式。采用相同网格构造出的紧致差分格式，可以达到比传统差分更高的精度，同时具有更高的尺度分辨率以及更小的波相位误差。一般的四阶方法，在空间上要涉及到五个节点，而紧致差分格式的四阶精度仅仅在空间上涉及到三个节点，大大简化了运算。Lele探讨了均匀网格下的七点差分格式，各项系数用Tailor级数展开得到。随后Chu讨论了在均匀网格和非均网格条件下的三点六阶紧致差分格式。三点紧致差分中间节点如下：

αf″_i-1+f″_i+βf″_i+1＝Af_i+1+Bf_i+Cf_i-1 (29)

系数：

边界节点i＝1:

f″₁+αf″₂＝Af₁+Bf₂+Cf₃ (31)

其中系数:

边界节点i＝n:

1.3均匀化网格：Fast Marching Method需要在网格上完成各点的到达时间值计算，而不是在散乱数据点上进行计算，因此在FMM计算之前需要对数据进行网格化，现有网格大多采用均匀网格划分，每个单元格成为一个正方形(2维)或正方体(3维)，将单元格顶点坐标索引(i,j,k)作为网格坐标，由于点云采样的非均匀性或各向异性，如果对点云进行均匀网格划分，一方面将导致真实数据点偏离网格顶点，因而测地路径计算出来的点可能不是点云数据的真实坐标，导致测地路径计算出现偏差。如图3a所示，在均匀网格化下数据点(0.8,1.1)将被网格点(0.5,1.0)所替代，图3b的非均匀网格下点(0.8,1.1)的坐标值不变；另一方面数据点在各坐标轴之间的间隔不能保证各坐标轴向间隔有公共的长度因子(比如点云数据p(x,y,z)的坐标范围x∈[1,3]，y∈[2,6.8]，z∈[-2.3,7.4]，x、y、z轴的坐标长度分别为2、4.8和9.7，那么如何来计算单元格长度来确保网格化后的单元格为立方体呢)，即使有的话，也将导致单元格太小使计算量大大增加，本专利将网格非均匀划分后统一用网格左下角顶点坐标值表示网格坐标位值。

本发明技术实现过程：

本专利采用以下步骤来计算点云模型上任意两点p₀和p_n之间的测地线：

步骤101：根据dx＝|x_n-x₀|，dy＝|y_n-y₀|，dz＝|z_n-z₀|，确定测地线起点p₀(x₀,y₀,z₀)和终点p_n(x_n,y_n,z_n)之间的主方向，以主方向来确定网格化区域：以减少网格化时间、FMM计算时间和路径跟踪时间；

步骤102：找出网格化区域中所有点云数据点坐标的最大值与最小值x_min、x_max、y_min、y_max、z_min与z_max；分别计算数据点在三个轴向的步长hx_i，hy_j和hz_k，非均匀化两点p₀和p_n之间的数据点，然后执行步骤103(点云数据点用数值坐标表示，如p点的坐标为(0.5,-1.8,3.4)，FMM计算是在单元格索引上进行，如网格化后p点的单元格在网格内的索引可能为(3,5,12)，测地路径也是根据单元格索引来跟踪的，路径跟踪完毕后，最后还得将单元格索引换算为单元格的数值坐标位置，因此非均匀网格化可避免单元格索引与点云数据的坐标位置偏差，提高测地路径计算精度)；

步骤103：根据步长hx_i，hy_j和hz_k，采用UNCDFMM计算波前网格各单元格的达到最小时间T(FMM方法在单元格数值过程中，均需用到近似差分(前向差分、后向差分)，均匀化网格由于单元格在各轴向的间距相等，不会出现精度问题。但非均匀网格由于单元格在各轴向的间距不等，再利用上述差分近似方法，将出现较大的精度误差，因此采用紧致差分)；

步骤104：采用正向跟踪法生成测地路径：测地路径的反向跟踪法，采用梯度法(最速下降法)，由于数值计算误差，往往导致超越网格边界，怎么取下一个路径点，没有好的方法进行解决。因此利用测地线的微分性质，可避免上述问题。

其中步骤101具体步骤是：

点云数据量通常比较大，估算任意两点间的测地线，将点云数据全部网格化没有必要：(1)因为测地线是一条短程线，这条曲线一定是某个区域上点与点之间的顺序连线，无须对全部数据进行网格化；(2)利用FMM方法行进计算网格点数值时，意味着全部网格点都要计算，这是一个非常耗时的过程；(3)FMM计算完成后，需要从网格上采用某种规则选择相应的单元格构成测地路径，如果数据全部网格化，路径跟踪计算量也非常大。因此合理选择网格化区域，减少单元格规模，可以减少网格的FMM计算时间和测地线跟踪时间，大大提升测地线生成速度。

本专利首先找出测地线起点p₀和终点p_n的主方向，以主方向来确定网格化区域。主方向的确定方法如下：

步骤1011：dx＝|x_n-x₀|，dy＝|y_n-y₀|，dz＝|z_n-z₀|；:判断dx、dy以及dz大小，当dx最大，执行步骤1012；当dy最大时，执行步骤1013；当dz最大时，执行步骤1014；

步骤1012：则取x轴为主方向，如果x₀<x_n，选择点云数据中满足x₀≤x≤x_n的点来当做网格化区域，x<x₀或x>x_n在测地线计算中无意义，如果测地线经过这类数据点，测地线将会折回来，不满足测地线最短性质；如果x₀>x_n，选择点云数据中满足x₀≥x≥x_n的点当做网格化区域；

步骤1013：dy最大，则取y轴为主方向，如果y₀<y_n，选择点云数据中满足y₀≤y≤y_n的点来当做网格化区域；如果y₀>y_n，选择点云数据中满足y₀≥y≥y_n的点来当做网格化区域；

步骤1014：dz最大，则取z轴为主方向，如果z₀<z_n，选择点云数据中满足z₀≤z≤z_n的点当做网格化区域；如果z₀>z_n，选择点云数据中满足z₀≥z≥z_n的点当做网格化区域。

其中，步骤102具体实现过程是：

网格的非均匀化过程：

步骤1021：找出所有点云数据点坐标的最大值与最小值x_min、x_max、y_min、y_max、z_min与z_max；

步骤1022：假设区域内原数据点的坐标为(x_m,y_m,z_m)(1≤m≤N)，取出数据点的三个坐标分量组成数组并进行排序并剔除相同的值，得到坐标分量数组，x_i(1≤i≤I)，y_j(1≤j≤J)和z_k(1≤k≤K)(由于在排序过程中坐标值相同的坐标分量(刻度)只保留一个，因此坐标分量排序后，其坐标刻度数量可能会小于N，I,J,K≤N)

步骤1023：分别计算数据点在三个轴向的间隔(步长)hx_i＝x_i+1-x_i(2≤i≤I-1)，hy_j＝y_j+1-y_j(2≤j≤J-1)，hz_k＝z_k+1-z_k(2≤k≤K-1)，令hx₁＝0，hy₁＝0，hz₁＝0。

步骤1024：在坐标刻度末尾补ε＝1.0×e^-3，如果x_I<x_max+ε，则hx_I＝ε否则hx_I＝0，同理hy_J＝ε否则hy_J＝0，hz_K＝ε否则hz_K＝0，以使步长hx_i，hy_j和hz_k的数目与各轴向方向的单元格数目相同，方便FMM计算。

其中步骤103具体过程是：对点云区域进行网格化之后，FMM将根据波的传播原理在网格上计算各单元格的到达时间值。

OFMM和HAFMM是在正交均匀网格下(至少在同一坐标轴上的间距相等)，通过选择当前单元格每一维前后对称邻近单元格最小值构造一元二次方程，求解方程来得到网格点的到达时间值，与OFMM和HAFMM不同的是，MSFM通过选择空间正交单元格对角线上邻近点来构造方程，

本专利中网格是非均匀的，FMM需要在非均匀网格下完成计算，为提高计算精度：采用如下步骤：

步骤1031：采用紧致差分来计算网格点的到达时间值，因此采用单向三点二阶紧致差分来计算网格点到达时间值，前向紧致差分和后向紧致差分表达式分别为：

前向紧致差分当i＝n-1时，用前向一阶差分

近似，当i＝n时前向紧致差分不存在。后向紧致差分当i＝2时，用后向一阶差分

近似，当i＝1时前向紧致差分不存在。

前向紧致差分的系数计算如下：假设区间[a,b]的非均匀划分为a＝x₀<x₁…x_n＝b，令h_i＝x_i-x_i-1(i＝1，2,…,n)，对T_i+1和T_i+2进行泰勒级数展开，得到：

比较式(39)的系数，可得到如下方程组：

求解方程组(40)，可得到(3)式前向紧致差分系数：

后向紧致差分的系数计算如下：假设区间[a,b]的非均匀划分为a＝x₀<x₁…x_n＝b，令h_i＝x_i-x_i-1(i＝1，2,…,n)，对T_i-1和T_i-2进行泰勒级数展开，得到：

比较式(44)的系数，可得到如下方程组：

求解方程组(45)，可得到后向差分式(4)的系数

(3)式和(4)式的系数分别见(41)和(46)所示。因此

与

分别为：

将(47)式和(48)式带入(24)式(实际上是用(47)和(48)式的二阶差分代替(25)式的一阶差分)，得到在单元格[i,j,k]处方程：

(19)式是在X,Y,Z三个方向的方程，更一般地，在26领域方案中将有26个方向。

和

分别表示以[i,j,k]为公共端点两共线射线方向。

步骤1032：(OFMM和HAFMM在为波前单元格计算到达时间值时，从波前单元格的6个或26个邻近单元格(窄带)选择到达时间值最小的单元格，以这些单元格的时间值来构成一元二次方程，求解该方程，取大于所有邻近单元格时间值且最小的解作为波前单元格的时间值。对于采用26领域方案中，由于邻近单元格与波前单元格非正交，同时多个可能的行进方向将使一元二次方程变得复杂，系数项数复杂，给方程求解带来麻烦)

本专利步骤将当前的邻近单元格分为波经过的单元格Frn(Frozen状态)和未经过的单元格O(Narrow Band状态或Open状态)两类，邻近单元格划分示意图见图5所示。通过生成Frn的幂集Frn^p(假如Frn中的单元格分别为1,2,3，那么幂集Frn^p为{{1}，{2}，{3}，{1,2}，{1,3}，{2,3}，{1,2,3}})，利用幂集的每个子集来计算当前点的到达时间值，并取这些时间的最小值作为当前单元格的到达时间值。

单元格时间值计算过程：

定义一个暂存单元格时间值数组U，令当前单元格为C_i,j,k，依次从Frn^p中选择一个集合S；

第1步：如果集合S的成员个数为1，即集合仅包含一个单元格C_l,m,n，则

将T_i,j,k，加入数组U中。其中|α-β|表示单元格α和β之间的距离，F表示波的传播速度；

第2步：如果集合S的成员个数为2个，即集合包含两个单元格

和

如果

与

在同一方向上，则将

和

带入(49)式，构造一元二次方程求解网格时间值T_i,j,k，加入数组U中；

(在三维网格中C_i,j,k邻近单元格有26个，因此存在26个方向，如果把C_i,j,k简记为[i,j,k],则26个方向为

[i,j,k+1]、[i,j,k-1]、[i,j+1,k]、[i,j+1,k+1]、[i,j+1,k-1]、

[i,j-1,k]、[i,j-1,k+1]、[i,j-1,k-1]、[i+1,j,k]、[i+1,j,k+1]、

[i+1,j,k-1]、[i+1,j+1,k]、[i+1,j+1,k+1]、[i+1,j+1,k-1]、[i+1,j-1,k]、

[i+1,j-1,k+1]、[i+1,j-1,k-1]、[i-1,j,k]、[i-1,j,k+1]、[i-1,j,k-1]、

[i-1,j+1,k]、[i-1,j+1,k+1]、[i-1,j+1,k-1]、[i-1,j-1,k]、[i-1,j-1,k+1]、[i-1,j-1,k-1]

)

如果C_i,j,k与

或

的下标能够满足如下顺序，则三个三元格在同一个方向上：

(

(1)[i,j,k]、[i,j,k+1]、[i,j,k+2]

(2)[i,j,k]、[i,j,k-1]、[i,j,k-2]

(3)[i,j,k]、[i,j+1,k]、[i,j+2,k]

(4)[i,j,k]、[i,j+1,k+1]、[i,j+2,k+2]

(5)[i,j,k]、[i,j+1,k-1]、[i,j+2,k-2]

(6)[i,j,k]、[i,j-1,k]、[i,j-2,k]

(7)[i,j,k]、[i,j-1,k+1]、[i,j-2,k+2]

(8)[i,j,k]、[i,j-1,k-1]、[i,j-2,k-2]

(9)[i,j,k]、[i+1,j,k]、[i+2,j,k]

(10)[i,j,k]、[i+1,j,k+1]、[i+2,j,k+2]

(11)[i,j,k]、[i+1,j,k-1]、[i+2,j,k-2]

(12)[i,j,k]、[i+1,j+1,k]、[i+2,j+2,k]

(13)[i,j,k]、[i+1,j+1,k+1]、[i+2,j+2,k+2]

(14)[i,j,k]、[i+1,j+1,k-1]、[i+2,j+2,k-2]

(15)[i,j,k]、[i+1,j-1,k]、[i+2,j-2,k]

(16)[i,j,k]、[i+1,j-1,k+1]、[i+2,j-2,k+2]

(17)[i,j,k]、[i+1,j-1,k-1]、[i+2,j-2,k-2]

(18)[i,j,k]、[i-1,j,k]、[i-2,j,k]

(19)[i,j,k]、[i-1,j,k+1]、[i-2,j,k+2]

(20)[i,j,k]、[i-1,j,k-1]、[i-2,j,k-2]

(21)[i,j,k]、[i-1,j+1,k]、[i-2,j+2,k]

(22)[i,j,k]、[i-1,j+1,k+1]、[i-2,j+2,k+2]

(23)[i,j,k]、[i-1,j+1,k-1]、[i-2,j+2,k-2]

(24)[i,j,k]、[i-1,j-1,k]、[i-2,j-2,k]

(25)[i,j,k]、[i-1,j-1,k+1]、[i-2,j-2,k+2]

(26)[i,j,k]、[i-1,j-1,k-1]、[i-2,j-2,k-2]

)

第3步：如果集合S的成员个数大于2个，令S＝{C₁,C₂,…C_Q}，分别将C_q(1≤q≤Q)在C_i,j,k处一阶渐进展开，可得：

其中sgn表示正负符号，“+”表示正号，“-”表示负号，正负符号确定根据单元格位置三个索引之差的乘积的正负符号来确定，例如：[2,3,5]-[1,4,6]＝[1,-1,-1]，sgn＝(+)×(-)×(-)＝,+T′_ijk由可由下式计算出：

其中

M⁺＝(M^tM)^-1M^t

计算出T′_ijk就得到单元格C_i,j,k邻近单元格集合S的时间值T_ijk，将其加入到数组U中。

第4步：幂集Frn^p的每个集合S都计算完毕，从数组U中选择时间值最小且值大于Frn中每个单元格时间值的值作为当前单元格C_i,j,k的时间值。

步骤104具体包括：

设S是E³的曲面，r＝r(u¹,u²)是曲面S的参数表示。C:r(s)＝r(u¹(s),u²(s))是曲面上的一条弧长参数曲线。沿曲线C取曲面的正交标架{e₁,e₂,e₃}，其中

e₃＝n，且e₁，e₂，e₃是正定向的，即是(e₁,e₂,e₃)＝<e₁,e₂∧e₃>>0。曲面S上的弧长参数曲线r＝r(s)的测地曲率

当k_g＝0时，曲线r为测地线。e₁是曲线的切线，e₂为副法线，e₃为法线。将测地线起点p₀与终点p_n的连线近似为曲线的切线

以p₀为当前点，根据e₁，e₂，e₃的正定向条件，可以计算出p₀的一组e₂方向和对应的e₃。为p₀确定出了e₂和e₃，就可以找到p₀的一个或多个可能行进方向，得到一组测地线从p₀可能经过的一组网格点{q_j}_j＝1…m，然后依次将{q_j}_j＝1…m中的点作为当前点，重复上述步骤，直到找到终点p_n为止。上述正向跟踪方法中，测地线的经过的节点可用树状结构来描述，树的根节点为p₀，树的叶节点均为p_n，因此得到测地线从p₀到p_n的多条路径，深度遍历树并计算其长度，选择长度最小的路径作为最终的测地路径。测地路径正向跟踪示意图见图6所示，实验结果见图7a到7e所示。

与梯度下降反向跟踪不同，正向跟踪方法，测地曲率采用与当前单元格邻近单元格的二阶差分近似，选择测地曲率绝对值最小的点所在的方向作为e₂，尽管网格划分是非均匀的，但单元格之间正交，对应的e₃方向必在与e₂正交方向上，放宽e₁，e₂，e₃正交条件，即不需要e₁同时垂直于e₂和e₃，但满足混合积(e₁,e₂,e₃)>0条件。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims

1.一种点云测地路径正向跟踪生成方法，其特征在于包括：

步骤3：根据步长hx_i，hy_j和hz_k，采用UNCDFMM计算波前网格各单元格的达到时间值；UNCDFMM指的是非均匀单向紧致差分快速行进法；

2.根据权利要求1所述的一种点云测地路径正向跟踪生成方法，其特征在于所述步骤1具体过程是：

步骤12：取x轴为主行进方向，如果x_s＜x_e，选择点云数据中x坐标满足x_s≤x≤x_e的点来当做网格化区域，x＜x_s或x＞x_e在测地路径计算中无意义，如果测地线路径过这类数据点，测地路径将会折回来，不满足测地路径最短性质；如果x_s＞x_e，选择点云数据中x坐标满足x_s≥x≥x_e的点当做网格化区域；

步骤13：取y轴为主行进方向，如果y_s＜y_e，选择点云数据中y坐标满足y_s≤y≤y_e的点来当做网格化区域；如果y_s＞y_e，选择点云数据中y坐标满足y_s≥y≥y_e的点来当做网格化区域；

步骤14：取z轴为主行进方向，如果z_s＜z_e，选择点云数据中z坐标满足z_s≤z≤z_e的点当做网格化区域；如果z_s＞z_e，选择点云数据中z坐标满足z_s≥z≥z_e的点当做网格化区域。

3.根据权利要求1所述的一种点云测地路径正向跟踪生成方法，其特征在于所述步骤2具体过程：

和

步骤23：分别计算新的坐标分量数组

和

形成的数据点在三个轴向的步长

令

此时2≤i≤I-1；2≤j≤J-1；2≤k≤K-1；

步骤24：在新的坐标分量数组

和

的最后一个分量末尾增加ε＝1.0×e^-3，如果x_I＜x_max+ε，则

否则

如果y_J＜y_max+ε，则

否则，

如果z_K＜z_max+ε，则

否则，

以使步长

和

的最大值；y_max是y轴分量数组

的最大值；z_max是z轴分量数组

的最大值。

4.根据权利要求1所述的一种点云测地路径正向跟踪生成方法，其特征在于所述步骤3具体过程：

步骤31：采用紧致差分来计算新的坐标分量数组

和

构成的网格点数值，构造单向三点二阶紧致差分计算网格点数值；X轴、Y轴和Z轴方向上单向三点二阶紧致差分格式相同，令T_i为分量数组

中x_i处的数值，令

为分量数组

中x_i处的单向前向二阶紧致差分值，令

为分量数组

中x_i处的单向后向二阶紧致差分值，则表达式分别为：

和

分别为前向二阶紧致差分表达式(1)的系数，当i＝I-1时，前向二阶紧致差分不存在，用前向一阶差分

近似前向二阶紧致差分值，当i＝I时前向二阶紧致差分不存在；

和

分别为后向二阶紧致差分表达式(2)的系数，当i＝2时，后向二阶紧致差分不存在，用后向一阶差分

近似后向二阶紧致差分值，当i＝1时后向二阶紧致差分不存在；

由于X轴、Y轴和Z轴方向上单向二阶紧致差分格式相同，单向二阶紧致差分的系数计算如下：对公式(1)中T_i+1和T_i+2进行泰勒级数展开，得到公式(3)和(4)；对公式(2)中T_i-1和T_i-2进行泰勒级数展开，得到公式(5)和(6)；

其中T_i′、T_i″、T_i″′和T_i ⁽⁴⁾分别表示坐标分量数组

上x_i处数值T_i对应的一阶、二阶、三阶和四阶导数；

表示

同理可知

和

意义；

将式(3)和式(4)代入式(1)，合并相同项，式(1)中

即是T_i″，得到:

比较式(7)中T_i、T_i′和T_i″系数，构造系数方程组并求解，可得式(1)的系数

和

为：

同理，将式(5)和式(6)代入式(2)，合并相同项，式(2)中

即是T_i″，得到:

比较式(9)中T_i、T_i′和T_i″系数，构造系数方程组并求解，可得式(2)的系数

和

为：

因此在X轴向上，(1)式和(2)式中单向前向二阶紧致差分

和单向后向二阶紧致差分

分别为：

同理可计算Y轴向的单向二阶紧致差分值：

和Z轴向的单向二阶紧致差分值：

由Eikonal方程数字解公式：

其中D^-x(T_ijk)、D^+x(T_ijk)、D^-y(T_ijk)、D^+y(T_ijk)、D^-z(T_ijk)和D^+z(T_ijk)分别表示T_ijk在X轴向的后向一阶差分、前向一阶差分、在Y轴向的后向一阶差分、前向一阶差分、在Z轴向的后向一阶差分和前向一阶差分，F_ijk为单元格[i,j,k]处的速度，max(D^-x(T_ijk),0)²表示取D^-x(T_ijk)与0两者的最大值再平方，和min(-D^+x(T_ijk),0)²表示取-D^+x(T_ijk)与0两者的最小值再平方，同理可知max(D^-y(T_ijk),0)²、min(-D^+y(T_ijk),0)²、max(D^-z(T_ijk),0)²、min(-D^+z(T_ijk),0)²表达的意义；

用二阶差分形式

和

和

则得到在单元格[i,j,k]处方程：

(17)式是在X、Y和Z三个方向的方程；更一般地，在26领域方案中有26个方向，此时

和

分别表示以[i,j,k]为公共端点两共线射线方向上的单向前向二阶紧致差分和单向后向二阶紧致差分；

步骤32：当前单元格邻近单元格分为波经过的单元格Frn和未经过的单元格O两类，通过生成Frn的幂集Frn^p，利用幂集的每个子集来计算当前单元格的数值，并取这些数值的最小值作为当前单元格的达到时间值。将每个未计算出到达时间值的单元格作为当前单元格，采用该步骤进行计算，最终得到波前网格各单元格的到达时间值；Frn代表Frozen状态的单元格；O代表Narrow Band状态或Open状态；幂集Frn^p包括至少一个单元格；Frozen状态表示波经过单元格，时间已经计算出，不会再改变，Narrow Band状态表示波经过单元格，时间值已经计算出，需要改变更新，Open状态表示波还未经过，时间值未知的状态。

5.根据权利要求4所述的一种点云测地路径正向跟踪生成方法，其特征在于所述步骤32中当前单元格C_ijk的到达时间值计算过程具体包括：