CN108180918A - 一种点云测地路径正向跟踪生成方法及装置 - Google Patents
一种点云测地路径正向跟踪生成方法及装置 Download PDFInfo
- Publication number
- CN108180918A CN108180918A CN201711228825.9A CN201711228825A CN108180918A CN 108180918 A CN108180918 A CN 108180918A CN 201711228825 A CN201711228825 A CN 201711228825A CN 108180918 A CN108180918 A CN 108180918A
- Authority
- CN
- China
- Prior art keywords
- point
- value
- cell
- ijk
- difference
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/26—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for navigation in a road network
- G01C21/34—Route searching; Route guidance
- G01C21/3446—Details of route searching algorithms, e.g. Dijkstra, A*, arc-flags, using precalculated routes
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Automation & Control Theory (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及点云数据处理领域。本发明选择测地路径起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的主行进方向,以主行进方向来确定网格化区域;提取网格化区域内所有数据点坐标值,构成三个坐标分量数组,排序数组并剔除相同坐标分量值,找出网格化区域中所有点云数据点坐标的最大值与最小值;分别计算数据点在三个轴向的步长hxi,hyj和hzk,非均匀网格化两点ps和pe之间的数据点,然后执行步骤3;根据步长hxi,hyj和hzk,采用UNCDFMM计算波前网格各单元格的达到时间值;计算出各单元格的到达时间后,筛选出满足近似测地线性质的单元格,将这些单元格的顶点顺序连接构成测地路径。
Description
技术领域
本发明涉及点云数据处理领域,尤其是一种点云测地路径正向跟踪生成方法及装置。
背景技术
在包含待测两点内的地球面上测得的两点之间的最短线。测地线又称大底线或短程线, 可以定义为空间中两点最短或最长路径。
航线:生活中飞机或轮船的航行都是测地线。
圆柱上的测地线:将圆柱面剪开铺平,给出了圆柱面到平面的一个等距变换,平面的 测地线是直线.因此圆柱面上的测地线:将平面卷成圆柱面时,由平面上直线变过来的曲 线,容易发现它们是直线(圆柱面的子母线),平行圆或圆柱螺线.自然界的攀缘植物,沿螺线生长,就是测地线的一个有趣例子。
球面上的测地线:由于大圆的主法线就是球面的法线,由此推出大圆或其一部分就是 球面上的测地线。在球面上,连接两点的大圆弧(测地线)有两条,这两条大圆弧一长一短, 短的是最短线,也就是短程线。
物理学中的测地线:假设一个质点在曲面上自由运动,如无外力作用,则质点在曲面 上的运动轨迹是一条测地线。再如一条无重量的弹性细线在一光滑曲面上自由移动,当这 条细线受曲面上两点问的某张力作用处于平衡状态时,则它取测地线的形状。
测地线在工程技术中的应用:复杂曲面的外板展开,这种板金技术在飞机机身、汽车 外壳、轮船船体、涡轮叶片、薄壳屋顶、机耕梨面等外形设计和制造。
测地线在制衣行业中的应用:帐篷脊线、鞋的足背线、衣服的腰线等。
1.2曲面测地线。
在微分几何中,如:
图1光滑曲面S的参数形式为r(u,v),P0=P(u0,v0)为曲面S上与(u0,v0)对应的点,(u(s),v(s))是参数域D上过点(u0,v0)的一条曲线,且u0=u(0),v0=v(0),设s是曲线r(s)的 弧长参数,则r(s)=r(u(s),v(s))是过曲面S上过点P0的一条曲线,且P0对应参数s=0。设T 为曲线在P0点的单位切向量,设n为曲面在P0点的单位法向量,设β为曲线在P0点的单位主法向量,设κβ为曲线在P0点的曲率向量。曲线r(s)在p0点的单位切向量为:
r(s)在p0点的曲率向量为:
曲线r(s)的弯曲由两部分产生:曲线随曲面的弯曲产生的法曲率kn和曲线自身弯曲 产生的测地曲率kg。则κβ可分解为向曲面法向量n的投影和过P0切平面的投影(副法线方向),也即是:
kn=kβ·n,kg=kβ·(n×T);
测地线是曲面上测地曲率kg=0的曲线,即曲线的主法向量与曲面的法向量平行,同 时也是曲面上局部距离最短的曲线,又称最短程线。
光滑曲面上可通过求解测地微分方程得到一条经过给定切方向的精确测地线,但对于 网格模型,无法得到模型的参数表达式,因此只能采用近似算法来计算得到一组离散点, 称为近似测地路径,测地路径上点的依次连线称为近似测地线。赵i对离散网格模型上的测 地线做了综述,指出离散网格模型上的近似测地线由于不能保证同时满足光滑参数曲面上 精确测地线的条件,离散网格模型上的近似测地线就分为最直测地线和最短测地线。由于 最短测地线满足三角不等式,因而是度量,因此离散模型的测地线研究中多研究最短测地 线。
与离散网格模型不同,点云模型由散乱数据点构成,不具有模型表达式,也不具有网 格结构,无法计算其最直测地线,因此对点云模型测地线的研究主要是利用测地线局部最 短性质,近似计算其最短测地线。Memol定义球心在给定点的一组球来生成点云的偏置带 (offset band),通过fast marching method在偏置带构成的空间网格来计算近似测地线, 其精度决定于点云采样密度、速度决定于网格数,对噪声的鲁棒性取决于球的半径,由于 采用一组球面构成偏置带,无法区分角点、边或边界等特征。Kleinii通过在点云上构建几 何近似图(Delaunay图和SIG(spheres-of-influence graph))来计算近似测地线,Hofer在 点云上构建MSL曲面,将能量泛函作用于MSL曲面,并在曲面上施加约束使其最小化来计 算测地线,该方法对噪声敏感,计算的测地线精度依赖于构建的MSL曲面。与Hofer类似, Ruggeriiii也运用能量泛函来计算测地线,能量函数由相连两点间的平方距离和L(P)、点贴 近曲面片的惩罚项Ds(P)和两点间连线贴近曲面片的惩罚项Us(P)构成,并在Ds(P)和 Us(P)上施加可调因子α和β来控制收敛速度和近似精度,初始路径利用Dijkstra算法计 算。与Hofer不同的是后者无需将约束施加于构建的MSL曲面上,而是在能量泛函中加入 约束,避免对点云模型曲面的依赖。肖采用MLS对点云进行采样和去噪处理,利用窄带Level Set生成一条虚拟路径,最后通过计算点云上点到虚拟路径最小的点来得到测地线。杜用 Dijkstra算法找到两点间的初始测地路径,再用抛物面拟合点云邻近域,估算抛物面的法 线,最小化平方距离度量函数来计算测地线。KEENANiv采用热核方法计算点云测地线,但在 热核离散化计算过程中利用了点云Laplace-Betrami算子矩阵,该矩阵计算需要利用点云 的法线,Yuv首先对点云进行空间单元格划分,然后利用连续Dijkstra算法计算一条近似 测地路径,最后利用测地曲率流对测地线进行校正,该方法比较有效,但在校正过程中需 要用到点云法线。
现有计算点云测地线的方法,主要有两类。(1)依赖点云法线的方法:在测地线计算过 程中需要用到点云法线的方法都存在先天缺陷,因为点云法线是不知道的,需要对目标点 云的法线进行估算来得到,法线估算精度不高,法线方向不能全局一致,出现混乱,因而 测地线计算精度差(2)单元格划分的方法:在对实物模型进行采样过程中,各区域几何形 状不同,采样密度不同。较平坦的区域通常用较大的间隔(较小的密度),弯曲部分用较小 的间隔(较大的密度),由于采样的各向异性,现有网格均匀化方法对点云的进单元格划分, 由大到小不断细分,直到单元格仅包含一个点云数据点,这种划分是均匀的(X,Y,Z三个方 向的间距相等),导致点云数据点往往位于单元格的内部,以单元格的位置来代替数据点的 位置坐标,最后生成的测地线几乎不经过点云的任何数据点(点云数据不在测地线上),未 能忠实于点云数据,因此一种精度不高的计算方法。
发明内容
本发明所要解决的技术问题是:针对现有技术存在的问题,提供一种点云测地路径正 向跟踪生成方法。根据点云的密度和间距实现非均匀网格划分,以使点云数据点全部位于 网格的顶点上,保证了后续关于单元格的各项计算忠实于点云数据。
本发明采用的技术方案如下:
一种点云测地路径正向跟踪生成方法包括:
步骤1:根据dx=|xe-xs|,dy=|ye-ys|,dz=|ze-zs|确定测地路径起点ps(xs,ys,zs)和 终点pe(xe,ye,ze)之间的主行进方向,以主行进方向来确定网格化区域;
步骤2:提取网格化区域内所有数据点坐标值,构成三个坐标分量数组,排序数组并剔 除相同坐标分量值,找出网格化区域中所有点云数据点坐标的最大值与最小值,记为:(xmax, ymax,zmax),(xmin,ymin,zmin);分别计算数据点在三个轴向的步长hxi,hyj和hzk,非均 匀化起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的数据点,然后执行步骤3;
步骤3:根据步长hxi,hyj和hzk,采用UNCDFMM计算波前网格各单元格的达到时间值;
步骤4:计算出各单元格的到达时间后,筛选出满足近似测地线性质的单元格,将这些 单元格的顶点顺序连接构成测地路径。
进一步的,所述步骤1具体过程是:
步骤11:dx=|xe-xs|,dy=|ye-ys|,dz=|ze-zs|;:比较dx、dy及dz值大小,当dx 最大,执行步骤12;当dy最大时,执行步骤13;当dz最大时,执行步骤14;
步骤12:取x轴为主行进方向,如果xs<xe,选择点云数据中x坐标满足xs≤x≤xe的点来当做网格化区域,x<xs或x>xe在测地路径计算中无意义,如果测地线路径过这类数据点,测地路径将会折回来,不满足测地路径最短性质;如果xs>xe,选择点云数据中x坐 标满足xs≥x≥xe的点当做网格化区域;
步骤13:取y轴为主行进方向,如果ys<ye,选择点云数据中y坐标满足ys≤y≤ye的点来 当做网格化区域;如果ys>ye,选择点云数据中y坐标满足ys≥y≥ye的点来当做网格化区 域;
步骤14:取z轴为主行进方向,如果zs<ze,选择点云数据中z坐标满足zs≤z≤ze的点当做网格化区域;如果zs>ze,选择点云数据中z坐标满足zs≥z≥ze的点当做网格化区域。
进一步的,所述步骤2具体过程:
步骤21:找出所有点云数据点坐标的最大值与最小值,记为:(xmax,ymax,zmax),(xmin, ymin,zmin);
步骤22:假设网格区域内数据点数量为N,原数据点的坐标为pm(xm,ym,zm),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,如果xI<xmax+ε,则否则如果yJ<ymax+ε,则否 则,如果zK<zmax+ε,则否则,以使步长和的数目 与各轴向方向的单元格数目相同,方便UNCDFMM计算;xmax是x轴分量数组的最大值; ymax是y轴分量数组的最大值;zmax是z轴分量数组的最大值。
进一步的,所述步骤3具体过程:
步骤31:采用紧致差分来计算新的坐标分量数组和构成的网格点 数值,构造单向三点二阶紧致差分计算网格点数值;X轴、Y轴和Z轴方向上单向三点二阶 紧致差分格式相同,以X轴向为例,不妨令Ti为分量数组中xi处的数值,令为分量数组中xi处的单向前向二阶紧致差分值,令为分量数组中xi处的单向后向二阶紧致差分值,则表达式分别为:
和分别为前向二阶紧致差分表达式(1)的系数,当i=I-1时,前向二阶紧致 差分不存在,用前向一阶差分近似前向二阶紧致差分值,当i=I时前向 二阶紧致差分不存在。和分别为后向二阶紧致差分表达式(2)的系数,当i=2时, 后向二阶紧致差分不存在,用后向一阶差分近似后向二阶紧致差分值,当 i=1时后向二阶紧致差分不存在。
由于X轴、Y轴和Z轴方向上单向二阶紧致差分格式相同,以X轴向为例,单向二阶紧致差分的系数计算如下:对公式(1)中Ti+1和Ti+2进行泰勒级数展开,得到公式(3)和(4);对公式(2)中Ti-1和Ti-2进行泰勒级数展开,得到公式(5)和(6)。
其中T'i、T”i、T”'i和分别表示坐标分量数组上xi处数值Ti对应的一阶、二阶、 三阶和四阶导数;表示同理可知和意义。
将式(3)和式(4)代入式(1),合并相同项,式(1)中即是T”i,得到:
比较式(7)中Ti、T'i和T”i系数,构造系数方程组并求解,可得式(1)的系数和为:
同理,将式(5)和式(6)代入式(2),合并相同项,式(2)中即是T”i,得到:
比较式(9)中Ti、T'i和T”i系数,构造系数方程组并求解,可得式(2)的系数和为:
因此在X轴向上,(1)式和(2)式中单向前向二阶紧致差分和单向后向二阶紧致差分 分别为:
同理可计算Y轴向的单向二阶紧致差分值:
和Z轴向的单向二阶紧致差分值:
X轴、Y轴和Z轴坐标分量数组构成空间网格,X轴向第i分量对应数值Ti,Y轴向第j分 量对应数值Tj,Z轴向第k分量对应数值Tk,[i,j,k]在三维网格中对应数值Tijk的位置索引;
由Eikonal方程数字解公式:其中D-x(Tijk)、D+x(Tijk)、D-y(Tijk)、D+y(Tijk)、D-z(Tijk)和D+z(Tijk)分别表示Tijk在X轴向的后 向一阶差分、前向一阶差分、在Y轴向的后向一阶差分、前向一阶差分、在Z轴向的后向 一阶差分和前向一阶差分,Fijk为单元格[i,j,k]处的速度,max(D-x(Tijk),0)2表示取D-x(Tijk)与 0两者的最大值再平方,和min(-D+x(Tijk),0)2表示取-D+x(Tijk)与0两者的最小值再平方,同 理可知max(D-y(Tijk),0)2、min(-D+y(Tijk),0)2、max(D-z(Tijk),0)2、min(-D+z(Tijk),0)2表达的意义。
用二阶差分形式和代替Eikonal公式中一阶差分形式D-和D+,并分别赋予X、Y 和Z轴向和则得到在单元格[i,j,k]处方程:
(17)式是在X、Y和Z三个方向的方程;更一般地,在26领域方案中有26个方向,此时和分别表示以[i,j,k]为公共端点两共线射线方向上的单向前向二阶紧致 差分和单向后向二阶紧致差分。
步骤32:当前单元格邻近单元格分为波经过的单元格Frn和未经过的单元格O两类, 通过生成Frn的幂集Frnp,利用幂集的每个子集来计算当前点的数值,并取这些数值的最 小值作为当前单元格的达到时间值,最终得到波前网格各单元格的到达时间值;Frn代表 Frozen状态的单元格;O代表Narrow Band状态或Open状态;幂集Frnp包括至少一个单元格;Frozen状态表示波经过单元格,时间已经计算出,不会再改变,Narrow Band状态表 示波经过单元格,时间值已经计算出,需要改变更新,Open状态表示波还未经过,时间值 未知的状态。
进一步的,所述步骤32单元格的到达时间值计算过程具体包括:
步骤321:定义一个暂存单元格数值数组U,不妨令当前单元格为Cijk,依次从Frnp中 选择一个集合Set;
步骤322:如果集合Set的成员个数为1,即集合仅包含一个单元格Clmn(l,m和n为单元格空间位置索引),则将Tijk加入数组U中。其中|Cijk-Clmn|表示单元格Cijk和Clmn之间的距离,Fijk表示波在索引为[i,j,k]处单元格的传播速度;Tlmn表示索引为[l,m,n]单元格的时间值。
步骤323:如果集合Set的成员个数为2个,即集合包含两个单元格和如果与Cijk在同一方向上,则将和的值带入(17)式,构造一元二次方程求解单元格的到达时间值Tijk,加入数组U中;
三个单元格在同一个方向上判断规则:如果Cijk、和的下标依顺序或满足如下顺序关系(以顺序关系(1)为例,[i,j,k] 为Cijk的下标,[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={C1,C2,…CQ},分别将 Cq(1≤q≤Q)在Cijk处一阶渐进展开,可得:
其中sgn表示符号“+”或“-”,符号根据单元格位置三个索引之差的乘积符号来确定(例 如:[2,3,5]-[1,4,6]=[1,-1,-1],sgn=(+)×(-)×(-)=+);由可由下式计算:
其中
计算出就得到单元格Cijk邻近单元格集合Set的到达时间值Tijk,将其加入到数组U中; Mt表示矩阵M转置;(MtM)-1表示矩阵MtM的逆矩阵。
步骤324:幂集Frnp的每个集合Set都计算完毕,从数组U中选择数值最小且值大于Frn 中每个单元格数值作为当前单元格Cijk的达到时间值。
进一步的,所述步骤4具体包括:
步骤41:设S是E3曲面,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)测地曲率当kg=0时,曲线r为 测地线;将测地线起点ps与终点pe的连线近似为曲线的切线以ps为当前点,根据的正定向条件,可以计算出ps的一组方向和对应的为ps确定出了和就可以找到ps的一个或多个可能行进方向,得到一组测地线从ps可能经过的一组网格点{qj}j=1…m,然后依次将{qj}j=1…m中的点作为当前点,重复上述步骤,直到找到终点pe为止;上述测地线的经过的节点可用树状结构来描述,树的根节点为ps,树的叶节点均为pe,因此得到测地线从ps到pe的多条路径,深度遍历树并计算其长度,选择长度最小的路径作为最终的测地路径;其中测地曲率kg采用与当前单元格邻近单元格时间值的二阶差分近似, 选择测地曲率绝对值最小的点所在的方向作为尽管网格划分是非均匀的,但单元格之 间正交,对应的方向必在与正交方向上。
一种点云测地路径正向跟踪生成装置包括:
主行方向确认模块,用于根据dx=|xe-xs|,dy=|ye-ys|,dz=|ze-zs|确定测地路径起 点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的主行进方向,以主行进方向来确定网格化区域;
步长计算模块,用于提取网格化区域内所有数据点坐标值,构成三个坐标分量数组,排 序数组并剔除相同坐标分量值,找出网格化区域中所有点云数据点坐标的最大值与最小值, 记为:(xmax,ymax,zmax),(xmin,ymin,zmin);分别计算数据点在三个轴向的步长hxi,hyj和hzk,非均匀化起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的数据点;
达到时间值计算模块,用于根据所述步长hxi,hyj和hzk,采用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,其演化界面在二维情况下为一平面曲线,在三维情况下为一曲面。其离 散解在三维情况下为vi:
其中:
离散FMM在网格上采用迭代方法实现,每个单元格将被标记为三种类型之一:
Open:波还未传播到该单元格,时间T未知的;
Narrow Band:波界面下一步将会扩展的候选单元格,具有即将更新的时间T;
Frozen:波界面已经经过该网格,时间T是固定的。
对于OFMM(Ordinary Fast Marching Method),由(45)选择同维网格最小值T1、T2和T3, 将其带入(43)即可得到关于T的一元二次方程,方程系数分别为a,b,c,求解出方程解 不放假设T3>T2>T1(总可以通过排序使之成立),根据迎风条件,其解 为:
(1)当T>max(T1,T2,T3)时,
(2)当T2<T<T3时,选择T2和T1重新生成一元二次方程,方程系数分别为a1,b1,c1,此时
(3)当T1<T<T2时,
由于一阶差分近似计算结果误差较大,HAFMM(Higher Accuracy Fast MarchingMethod), 利用二阶差分来近似梯度计算网格点的数值,Rickett vii和Baerentzen viii利用 分别来近似网格点的后向差分和前向差分,该方法中网格等距且节点的系数分别定为3,-4和1。同样方式构造Y方向和Z方 向的前向和后向差分,由(46)式选择同维网格最小值T1、T2和T3,与OFMM方法类似,通过 一元二次方程的解,在满足迎风条件下得到网格点的值T。
在HAFMM中,计算一个单元格的值在一坐标轴方向上将涉及附近4个单元格,这些单元 格可能处于“Open”、“Narrow Band”或“Frozen”的状态之一,当处于前两种状态下时,二阶差分近似计算将得到一个无效值,因此HAFMM需要满足以下条件:
(1)与当前单元格距离为2个单位的单元格必须处于“Frozen”状态,如x方向上的Ti-2,j,k和Ti+2,j,k应处于“Frozen”状态;
(2)与当前单元格距离为2个单位的单元格的数值必须小于等于距离为1个单位的单元 格数值,如x方向上,Ti-2,j,k≤Ti-1,j,k与Ti+2,j,k≤Ti+1,j,k;
当HAFMM不满足以下条件时,用一阶差分近似梯度代替二阶差分近似。
1.2 Compact Difference(紧致差分)
紧致差分通常用Hermite公式表示
其实质就是f(x)邻近节点函数值、一阶导数值及二阶导数值组合而成的表达式。采用 相同网格构造出的紧致差分格式,可以达到比传统差分更高的精度,同时具有更高的尺度 分辨率以及更小的波相位误差。一般的四阶方法,在空间上要涉及到五个节点,而紧致差 分格式的四阶精度仅仅在空间上涉及到三个节点,大大简化了运算。Leleix探讨了均匀网格 下的七点差分格式,各项系数用Tailor级数展开得到。随后Chu讨论了在均匀网格和非均 网格条件下的三点六阶紧致差分格式。三点紧致差分中间节点如下:
αf”i-1+f”i+βf”i+1=Afi+1+Bfi+Cfi-1 (48)
系数:
边界节点i=1:
f”1+αf”2=Af1+Bf2+Cf3 (50)
其中系数:
边界节点i=n:
f”n+αf”n-1=Afn+Bfn-1+Cfn-2 (52)
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,那么如何来计算单元格长度来确保网格化后的单元格为立方体呢), 即使有的话,也将导致单元格太小使计算量大大增加,本专利将网格非均匀划分后统一用 网格左下角顶点坐标值表示网格坐标位值。
本发明技术实现过程:
本专利采用以下步骤来计算点云模型上任意两点p0和pn之间的测地线:
步骤101:根据dx=|xn-x0|,dy=|yn-y0|,dz=|zn-z0|,确定测地线起点p0(x0,y0,z0) 和终点pn(xn,yn,zn)之间的主方向,以主方向来确定网格化区域:以减少网格化时间、FMM 计算时间和路径跟踪时间;
步骤102:找出网格化区域中所有点云数据点坐标的最大值与最小值xmin、xmax、ymin、 ymax、zmin与zmax;分别计算数据点在三个轴向的步长hxi,hyj和hzk,非均匀化两点p0和pn之间的数据点,然后执行步骤103(点云数据点用数值坐标表示,如p点的坐标为 (0.5,-1.8,3.4),FMM计算是在单元格索引上进行,如网格化后p点的单元格在网格内的索引 可能为(3,5,12),测地路径也是根据单元格索引来跟踪的,路径跟踪完毕后,最后还得将单 元格索引换算为单元格的数值坐标位置,因此非均匀网格化可避免单元格索引与点云数据的坐标位置偏差,提高测地路径计算精度);
步骤103:根据步长hxi,hyj和hzk,采用UNCDFMM(Unilateral Nonuniform CompactDifference FMM,单向非均匀紧致差分)波前网格各单元格的达到最小时间T(FMM方法在单元格数值过程中,均需用到近似差分(前向差分、后向差分),均匀化网格由于单元格在各轴向的间距相等,不会出现精度问题。但非均匀网格由于单元格在各轴向的间距不等,再利用上述差分近似方法,将出现较大的精度误差,因此采用紧致差分);
步骤104:采用正向跟踪法生成测地路径:测地路径的反向跟踪法,采用梯度法(最速 下降法),由于数值计算误差,往往导致超越网格边界,怎么取下一个路径点,没有好的方 法进行解决。因此利用测地线的微分性质,可避免上述问题。
其中步骤101具体步骤是:
点云数据量通常比较大,估算任意两点间的测地线,将点云数据全部网格化没有必要: (1)因为测地线是一条短程线,这条曲线一定是某个区域上点与点之间的顺序连线,无须 对全部数据进行网格化;(2)利用FMM方法行进计算网格点数值时,意味着全部网格点都 要计算,这是一个非常耗时的过程;(3)FMM计算完成后,需要从网格上采用某种规则选择 相应的单元格构成测地路径,如果数据全部网格化,路径跟踪计算量也非常大。因此合理 选择网格化区域,减少单元格规模,可以减少网格的FMM计算时间和测地线跟踪时间,大大提升测地线生成速度。
本专利首先找出测地线起点p0和终点pn的主方向,以主方向来确定网格化区域。主方 向的确定方法如下:
步骤1011:dx=|xn-x0|,dy=|yn-y0|,dz=|zn-z0|;:判断dx、dy以及dz大小,当 dx最大,执行步骤1012;当dy最大时,执行步骤1013;当dz最大时,执行步骤1014;
步骤1012:则取x轴为主方向,如果x0<xn,选择点云数据中满足x0≤x≤xn的点来当 做网格化区域,x<x0或x>xn在测地线计算中无意义,如果测地线经过这类数据点,测地线将会折回来,不满足测地线最短性质;如果x0>xn,选择点云数据中满足x0≥x≥xn的点当做网格化区域;
步骤1013:dy最大,则取y轴为主方向,如果y0<yn,选择点云数据中满足y0≤y≤yn的点来当做网格化区域;如果y0>yn,选择点云数据中满足y0≥y≥yn的点来当做网格化区域;
步骤1014:dz最大,则取z轴为主方向,如果z0<zn,选择点云数据中满足z0≤z≤zn的点当做网格化区域;如果z0>zn,选择点云数据中满足z0≥z≥zn的点当做网格化区域。
其中,步骤102具体实现过程是:
网格的非均匀化过程:
步骤1021:找出所有点云数据点坐标的最大值与最小值xmin、xmax、ymin、ymax、zmin与zmax;
步骤1022:不妨假设区域内原数据点的坐标为(xm,ym,zm)(1≤m≤N),取出数据点的三 个坐标分量组成数组并进行排序并剔除相同的值,得到坐标分量数组,xi(1≤i≤I),yj(1≤j≤J)和zk(1≤k≤K)(由于在排序过程中坐标值相同的坐标分量(刻度)只保留一个,因此坐标分量排序后,其坐标刻度数量可能会小于N,I,J,K≤N)
步骤1023分别计算数据点在三个轴向的间隔(步长)hxi=xi+1-xi(2≤i≤I-1),hyj=yj+1-yj(2≤j≤J-1),hzk=zk+1-zk(2≤k≤K-1),令hx1=0,hy1=0,hz1=0。
步骤1024:在坐标刻度末尾补ε=1.0×e-3,如果xI<xmax+ε,则hxI=ε否则hxI=0,同理hyJ=ε否则hyJ=0,hzK=ε否则hzK=0,以使步长hxi,hyj和hzk的数目与各轴向方 向的单元格数目相同,方便FMM计算。
其中步骤103具体过程是:对点云区域进行网格化之后,FMM将根据波的传播原理在网 格上计算各单元格的到达时间。
OFMM和HAFMM是在正交均匀网格下(至少在同一坐标轴上的间距相等),通过选择当 前单元格每一维前后对称邻近单元格最小值构造一元二次方程,求解方程来得到网格点的 到达时间,与OFMM和HAFMM不同的是,MSFMx通过选择空间正交单元格对角线上邻近点来构 造方程,
本专利中网格是非均匀的,FMM需要在非均匀网格下完成计算,为提高计算精度:采用 如下步骤:
步骤1031:采用紧致差分来计算网格点的到达时间,因此采用单向三点二阶紧致差分 来计算网格点到达时间,前向紧致差分和后向紧致差分表达式分别为:
前向紧致差分当i=n-1时,用前向一阶差分近似,当i=n时前 向紧致差分不存在。后向紧致差分当i=2时,用后向一阶差分近似, 当i=1时前向紧致差分不存在。
前向紧致差分的系数计算如下:不妨假设区间[a,b]的非均匀划分为a=x0<x1…xn=b, 令hi=xi-xi-1(i=1,2,…,n),对Ti+1和Ti+2进行泰勒级数展开,得到:
比较式(58)的系数,可得到如下方程组:
求解方程组(59),可得到(1)式前向紧致差分系数:
后向紧致差分的系数计算如下:不妨假设区间[a,b]的非均匀划分为a=x0<x1…xn=b, 令hi=xi-xi-1(i=1,2,…,n),对Ti-1和Ti-2进行泰勒级数展开,得到:
比较式(63)的系数,可得到如下方程组:
求解方程组(64),可得到后向差分式(2)的系数
(1)式和(2)式的系数分别见(8)和(10)所示。因此与分别为:
将(66)式和(67)式带入(43)式(实际上是用(66)和(67)式的二阶差分代替(44)式的一 阶差分),得到在单元格[i,j,k]处方程:
(17)式是在X,Y,Z三个方向的方程,更一般地,在26领域方案中将有26个方向。和分别表示以[i,j,k]为公共端点两共线射线方向。
步骤1032:(OFMM和HAFMM在为波前单元格计算到达时间时,从波前单元格的6个或26个邻近单元格(窄带)选择到达时间值最小的单元格,以这些单元格的时间值来构成一元二次方程,求解该方程,取大于所有邻近单元格时间值且最小的解作为波前单元格的时间值。对于采用26领域方案中,由于邻近单元格与波前单元格非正交,同时多个可能的行进方向将使一元二次方程变得复杂,系数项数复杂,给方程求解带来麻烦)
本专利步骤将当前的邻近单元格分为波经过的单元格Frn(Frozen状态)和未经过的单 元格O(Narrow Band状态或Open状态)两类,邻近单元格划分示意图见图5所示。通过生 成Frn的幂集Frnp(假如Frn中的单元格分别为1,2,3,那么幂集Frnp为{{1},{2},{3},{1,2},{1,3},{2,3},{1,2,3}}),利用幂集的每个子集来计算当前点的到达时间,并取 这些时间的最小值作为当前单元格的到达时间。
单元格时间值计算过程:
定义一个暂存单元格时间值数组U,不妨令当前单元格为Ci,j,k,依次从Frnp中选择一 个集合S,
第1步:如果集合S的成员个数为1,即集合仅包含一个单元格Cl,m,n,则将Ti,j,k,加入数组U中。其中|α-β|表示单元格α和β之间的距 离,F表示波的传播速度;
第2步:如果集合S的成员个数为2个,即集合包含两个单元格和如果与Ci,j,k在同一方向上,则将和带入(17)式,构造一元二次方程求解网格时间值Ti,j,k,加入数组U中;
(在三维网格中Ci,j,k邻近单元格有26个,因此存在26个方向,如果把Ci,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]
)
如果Ci,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={C1,C2,…CQ},分别将Cq(1≤q≤Q) 在Ci,j,k处一阶渐进展开,可得:
其中sgn表示符号+或-,符号确定根据单元格位置三个索引之差的乘积符号来确定(例如:[2,3,5]-[1,4,6]=[1,-1,-1],那么sgn=(+)×(-)×(-)=+)。T'ijk由可由下式计算出:
其中
计算出T'ijk就得到单元格Ci,j,k邻近单元格集合S的时间值Tijk,将其加入到数组U中。
第4步.幂集Frnp的每个集合S都计算完毕,从数组U中选择时间值最小且值大于Frn 中每个单元格时间值的值作为当前单元格Ci,j,k的时间值。
步骤104具体包括:
设S是E3的曲面,r=r(u1,u2)是曲面S的参数表示。C:r(s)=r(u1(s),u2(s))是曲面上的 一条弧长参数曲线。沿曲线C取曲面的正交标架{e1,e2,e3},其中e3=n,且e1,e2, e3是正定向的,即是(e1,e2,e3)=<e1,e2∧e3>>0。曲面S上的弧长参数曲线r=r(s)的测地曲率 当kg=0时,曲线r为测地线。e1是曲线的切线,e2为副法线,e3为法线。将测地线起点p0与终点pn的连线近似为曲线的切线以p0为当前点,根据e1,e2, e3的正定向条件,可以计算出p0的一组e2方向和对应的e3。为p0确定出了e2和e3,就可以 找到p0的一个或多个可能行进方向,得到一组测地线从p0可能经过的一组网格点{qj}j=1…m, 然后依次将{qj}j=1…m中的点作为当前点,重复上述步骤,直到找到终点pn为止。上述正向 跟踪方法中,测地线的经过的节点可用树状结构来描述,树的根节点为p0,树的叶节点均 为pn,因此得到测地线从p0到pn的多条路径,深度遍历树并计算其长度,选择长度最小的 路径作为最终的测地路径。测地路径正向跟踪示意图见图6所示,实验结果见图7a到7e 所示。
与梯度下降反向跟踪不同,正向跟踪方法,测地曲率采用与当前单元格邻近单元格的二 阶差分近似,选择测地曲率绝对值最小的点所在的方向作为e2,尽管网格划分是非均匀的, 但单元格之间正交,对应的e3方向必在与e2正交方向上,放宽e1,e2,e3正交条件,即不 需要e1同时垂直于e2和e3,但满足混合积(e1,e2,e3)>0条件。
本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征 或任何新的组合,以及披露的任一新的方法或过程的步骤或任何新的组合。
Claims (7)
1.一种点云测地路径正向跟踪生成方法,其特征在于包括:
步骤1:根据dx=|xe-xs|,dy=|ye-ys|,dz=|ze-zs|确定测地路径起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的主行进方向,以主行进方向来确定网格化区域;
步骤2:提取网格化区域内所有数据点坐标值,构成三个坐标分量数组,排序数组并剔除相同坐标分量值,找出网格化区域中所有点云数据点坐标的最大值与最小值,记为:(xmax,ymax,zmax),(xmin,ymin,zmin);分别计算数据点在三个轴向的步长hxi,hyj和hzk,非均匀化起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的数据点,然后执行步骤3;
步骤3:根据步长hxi,hyj和hzk,采用UNCDFMM计算波前网格各单元格的达到时间值;
步骤4:计算出各单元格的到达时间后,筛选出满足近似测地线性质的单元格,将这些单元格的顶点顺序连接构成测地路径。
2.根据权利要求1所述的一种点云测地路径正向跟踪生成方法,其特征在于所述步骤1具体过程是:
步骤11:dx=|xe-xs|,dy=|ye-ys|,dz=|ze-zs|;:比较dx、dy及dz值大小,当dx最大,执行步骤12;当dy最大时,执行步骤13;当dz最大时,执行步骤14;
步骤12:取x轴为主行进方向,如果xs<xe,选择点云数据中x坐标满足xs≤x≤xe的点来当做网格化区域,x<xs或x>xe在测地路径计算中无意义,如果测地线路径过这类数据点,测地路径将会折回来,不满足测地路径最短性质;如果xs>xe,选择点云数据中x坐标满足xs≥x≥xe的点当做网格化区域;
步骤13:取y轴为主行进方向,如果ys<ye,选择点云数据中y坐标满足ys≤y≤ye的点来当做网格化区域;如果ys>ye,选择点云数据中y坐标满足ys≥y≥ye的点来当做网格化区域;
步骤14:取z轴为主行进方向,如果zs<ze,选择点云数据中z坐标满足zs≤z≤ze的点当做网格化区域;如果zs>ze,选择点云数据中z坐标满足zs≥z≥ze的点当做网格化区域。
3.根据权利要求1所述的一种点云测地路径正向跟踪生成方法,其特征在于所述步骤2具体过程:
步骤21:找出所有点云数据点坐标的最大值与最小值,记为:(xmax,ymax,zmax),(xmin,ymin,zmin);
步骤22:假设网格区域内数据点数量为N,原数据点的坐标为pm(xm,ym,zm),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,如果xI<xmax+ε,则否则如果yJ<ymax+ε,则否则,如果zK<zmax+ε,则否则,以使步长,和的数目与各轴向方向的单元格数目相同,方便UNCDFMM计算;xmax是x轴分量数组的最大值;ymax是y轴分量数组的最大值;zmax是z轴分量数组的最大值。
4.根据权利要求1所述的一种点云测地路径正向跟踪生成方法,其特征在于所述步骤3具体过程:
步骤31:采用紧致差分来计算新的坐标分量数组和构成的网格点数值,构造单向三点二阶紧致差分计算网格点数值;X轴、Y轴和Z轴方向上单向三点二阶紧致差分格式相同,不妨令Ti为分量数组中xi处的数值,令为分量数组中xi处的单向前向二阶紧致差分值,令为分量数组中xi处的单向后向二阶紧致差分值,则表达式分别为:
和分别为前向二阶紧致差分表达式(1)的系数,当i=I-1时,前向二阶紧致差分不存在,用前向一阶差分近似前向二阶紧致差分值,当i=I时前向二阶紧致差分不存在;和分别为后向二阶紧致差分表达式(2)的系数,当i=2时,后向二阶紧致差分不存在,用后向一阶差分近似后向二阶紧致差分值,当i=1时后向二阶紧致差分不存在;
由于X轴、Y轴和Z轴方向上单向二阶紧致差分格式相同,单向二阶紧致差分的系数计算如下:对公式(1)中Ti+1和Ti+2进行泰勒级数展开,得到公式(3)和(4);对公式(2)中Ti-1和Ti-2进行泰勒级数展开,得到公式(5)和(6);
其中Ti'、Ti”、Ti”'和分别表示坐标分量数组上xi处数值Ti对应的一阶、二阶、三阶和四阶导数;表示同理可知和意义;
将式(3)和式(4)代入式(1),合并相同项,式(1)中即是Ti”,得到:
比较式(7)中Ti、Ti'和Ti”系数,构造系数方程组并求解,可得式(1)的系数和为:
同理,将式(5)和式(6)代入式(2),合并相同项,式(2)中即是Ti”,得到:
比较式(9)中Ti、Ti'和Ti”系数,构造系数方程组并求解,可得式(2)的系数和为:
因此在X轴向上,(1)式和(2)式中单向前向二阶紧致差分和单向后向二阶紧致差分分别为:
同理可计算Y轴向的单向二阶紧致差分值:
和Z轴向的单向二阶紧致差分值:
X轴、Y轴和Z轴坐标分量数组构成空间网格,X轴向第i分量对应数值Ti,Y轴向第j分量对应数值Tj,Z轴向第k分量对应数值Tk,[i,j,k]在三维网格中对应数值Tijk的位置索引;
由Eikonal方程数字解公式:其中D-x(Tijk)、D+x(Tijk)、D-y(Tijk)、D+y(Tijk)、D-z(Tijk)和D+z(Tijk)分别表示Tijk在X轴向的后向一阶差分、前向一阶差分、在Y轴向的后向一阶差分、前向一阶差分、在Z轴向的后向一阶差分和前向一阶差分,Fijk为单元格[i,j,k]处的速度,max(D-x(Tijk),0)2表示取D-x(Tijk)与0两者的最大值再平方,和min(-D+x(Tijk),0)2表示取-D+x(Tijk)与0两者的最小值再平方,同理可知max(D-y(Tijk),0)2、min(-D+y(Tijk),0)2、max(D-z(Tijk),0)2、min(-D+z(Tijk),0)2表达的意义;
用二阶差分形式和代替Eikonal公式中一阶差分形式D-和D+,并分别赋予X、Y和Z轴向和则得到在单元格[i,j,k]处方程:
(17)式是在X、Y和Z三个方向的方程;更一般地,在26领域方案中有26个方向,此时和分别表示以[i,j,k]为公共端点两共线射线方向上的单向前向二阶紧致差分和单向后向二阶紧致差分;
步骤32:当前单元格邻近单元格分为波经过的单元格Frn和未经过的单元格O两类,通过生成Frn的幂集Frnp,利用幂集的每个子集来计算当前点的数值,并取这些数值的最小值作为当前单元格的达到时间值,最终得到波前网格各单元格的到达时间值;Frn代表Frozen状态的单元格;O代表Narrow Band状态或Open状态;幂集Frnp包括至少一个单元格;Frozen状态表示波经过单元格,时间已经计算出,不会再改变,Narrow Band状态表示波经过单元格,时间值已经计算出,需要改变更新,Open状态表示波还未经过,时间值未知的状态。
5.根据权利要求1所述的一种点云测地路径正向跟踪生成方法,其特征在于所述步骤32单元格的到达时间值计算过程具体包括:
步骤321:定义一个暂存单元格数值数组U,不妨令当前单元格为Cijk,依次从Frnp中选择一个集合Set;
步骤322:如果集合Set的成员个数为1,即集合仅包含一个单元格Clmn(l,m和n为单元格空间位置索引),则将Tijk加入数组U中;其中|Cijk-Clmn|表示单元格Cijk和Clmn之间的距离,Fijk表示波在索引为[i,j,k]处单元格的传播速度;Tlmn表示索引为[l,m,n]单元格的时间值;
步骤323:如果集合Set的成员个数为2个,即集合包含两个单元格和,如果与Cijk在同一方向上,则将和的值带入(17)式,构造一元二次方程求解单元格的到达时间值Tijk,加入数组U中;
三个单元格在同一个方向上判断规则:如果Cijk、和的下标依顺序或满足如下顺序关系,当是顺序关系(1)时,[i,j,k]为Cijk的下标,[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={C1,C2,…CQ},分别将Cq(1≤q≤Q)在Cijk处一阶渐进展开,可得:
其中sgn表示符号“+”或“-”,符号根据单元格位置三个索引之差的乘积符号来确定(例如:[2,3,5]-[1,4,6]=[1,-1,-1],sgn=(+)×(-)×(-)=+);由可由下式计算:
其中
计算出就得到单元格Cijk邻近单元格集合Set的到达时间值Tijk,将其加入到数组U中;Mt表示矩阵M转置;(MtM)-1表示矩阵MtM的逆矩阵;
步骤324:幂集Frnp的每个集合Set都计算完毕,从数组U中选择数值最小且值大于Frn中每个单元格数值作为当前单元格Cijk的达到时间值。
6.根据权利要求1所述的一种点云测地路径正向跟踪生成方法,其特征在于所述步骤4具体包括:
步骤41:设S是E3曲面,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)测地曲率当kg=0时,曲线r为测地线;将测地线起点ps与终点pe的连线近似为曲线的切线以ps为当前点,根据的正定向条件,可以计算出ps的一组方向和对应的为ps确定出了和就可以找到ps的一个或多个可能行进方向,得到一组测地线从ps可能经过的一组网格点{qj}j=1…m,然后依次将{qj}j=1…m中的点作为当前点,重复上述步骤,直到找到终点pe为止;上述测地线的经过的节点可用树状结构来描述,树的根节点为ps,树的叶节点均为pe,因此得到测地线从ps到pe的多条路径,深度遍历树并计算其长度,选择长度最小的路径作为最终的测地路径;其中测地曲率kg采用与当前单元格邻近单元格时间值的二阶差分近似,选择测地曲率绝对值最小的点所在的方向作为尽管网格划分是非均匀的,但单元格之间正交,对应的方向必在与正交方向上。
7.一种点云测地路径正向跟踪生成装置,其特征在于包括:
主行方向确认模块,用于根据dx=|xe-xs|,dy=|ye-ys|,dz=|ze-zs|确定测地路径起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的主行进方向,以主行进方向来确定网格化区域;
步长计算模块,用于提取网格化区域内所有数据点坐标值,构成三个坐标分量数组,排序数组并剔除相同坐标分量值,找出网格化区域中所有点云数据点坐标的最大值与最小值,记为:(xmax,ymax,zmax),(xmin,ymin,zmin);分别计算数据点在三个轴向的步长hxi,hyj和hzk,非均匀化起点ps(xs,ys,zs)和终点pe(xe,ye,ze)之间的数据点;
达到时间值计算模块,用于根据所述步长hxi,hyj和hzk,采用UNCDFMM计算波前网格各单元格的达到时间值;
测地路径形成模块,用于计算出各单元格的到达时间后,筛选出满足近似测地线性质的单元格,将这些单元格的顶点顺序连接构成测地路径。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711228825.9A CN108180918B (zh) | 2017-11-29 | 2017-11-29 | 一种点云测地路径正向跟踪生成方法及装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201711228825.9A CN108180918B (zh) | 2017-11-29 | 2017-11-29 | 一种点云测地路径正向跟踪生成方法及装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108180918A true CN108180918A (zh) | 2018-06-19 |
CN108180918B CN108180918B (zh) | 2021-04-30 |
Family
ID=62545544
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201711228825.9A Active CN108180918B (zh) | 2017-11-29 | 2017-11-29 | 一种点云测地路径正向跟踪生成方法及装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108180918B (zh) |
Cited By (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109241059A (zh) * | 2018-08-30 | 2019-01-18 | 百度在线网络技术(北京)有限公司 | 一种点云数据的构造方法、装置、电子设备及存储介质 |
CN109800814A (zh) * | 2019-01-25 | 2019-05-24 | 西南科技大学 | 曲线测量定位的不变特征量提取方法 |
CN110480075A (zh) * | 2019-08-26 | 2019-11-22 | 上海拓璞数控科技股份有限公司 | 基于点云数据的工件曲面轮廓补偿系统及方法及介质 |
CN111089592A (zh) * | 2019-12-13 | 2020-05-01 | 天津大学 | 一种离散曲面中测地线计算方法 |
CN111736602A (zh) * | 2020-06-16 | 2020-10-02 | 国网上海市电力公司 | 一种改进的轮式机器人路径规划方法 |
CN112799416A (zh) * | 2019-10-24 | 2021-05-14 | 广州极飞科技股份有限公司 | 航线生成方法、设备和系统、无人作业系统及存储介质 |
CN114270142A (zh) * | 2019-07-19 | 2022-04-01 | 祖克斯有限公司 | 非结构化的车辆路径规划器 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102306397A (zh) * | 2011-07-08 | 2012-01-04 | 中国科学院自动化研究所 | 点云数据网格化的方法 |
CN103810271A (zh) * | 2014-01-29 | 2014-05-21 | 辽宁师范大学 | 基于路径跟随的三维点云物体形状特征匹配方法 |
CN105631065A (zh) * | 2014-10-31 | 2016-06-01 | 北京临近空间飞行器系统工程研究所 | 一种基于背景网格的动网格方法 |
US20160154999A1 (en) * | 2014-12-02 | 2016-06-02 | Nokia Technologies Oy | Objection recognition in a 3d scene |
CN106600617A (zh) * | 2016-12-29 | 2017-04-26 | 中科宇图科技股份有限公司 | 基于曲率从Lidar点云数据提取建筑物轮廓线的方法 |
-
2017
- 2017-11-29 CN CN201711228825.9A patent/CN108180918B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102306397A (zh) * | 2011-07-08 | 2012-01-04 | 中国科学院自动化研究所 | 点云数据网格化的方法 |
CN103810271A (zh) * | 2014-01-29 | 2014-05-21 | 辽宁师范大学 | 基于路径跟随的三维点云物体形状特征匹配方法 |
CN105631065A (zh) * | 2014-10-31 | 2016-06-01 | 北京临近空间飞行器系统工程研究所 | 一种基于背景网格的动网格方法 |
US20160154999A1 (en) * | 2014-12-02 | 2016-06-02 | Nokia Technologies Oy | Objection recognition in a 3d scene |
CN106600617A (zh) * | 2016-12-29 | 2017-04-26 | 中科宇图科技股份有限公司 | 基于曲率从Lidar点云数据提取建筑物轮廓线的方法 |
Cited By (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109241059A (zh) * | 2018-08-30 | 2019-01-18 | 百度在线网络技术(北京)有限公司 | 一种点云数据的构造方法、装置、电子设备及存储介质 |
CN109800814A (zh) * | 2019-01-25 | 2019-05-24 | 西南科技大学 | 曲线测量定位的不变特征量提取方法 |
CN109800814B (zh) * | 2019-01-25 | 2022-08-09 | 西南科技大学 | 叶片曲线测量定位的不变特征量提取方法 |
CN114270142A (zh) * | 2019-07-19 | 2022-04-01 | 祖克斯有限公司 | 非结构化的车辆路径规划器 |
CN114270142B (zh) * | 2019-07-19 | 2024-01-23 | 祖克斯有限公司 | 非结构化的车辆路径规划器 |
CN110480075A (zh) * | 2019-08-26 | 2019-11-22 | 上海拓璞数控科技股份有限公司 | 基于点云数据的工件曲面轮廓补偿系统及方法及介质 |
CN112799416A (zh) * | 2019-10-24 | 2021-05-14 | 广州极飞科技股份有限公司 | 航线生成方法、设备和系统、无人作业系统及存储介质 |
CN112799416B (zh) * | 2019-10-24 | 2024-04-12 | 广州极飞科技股份有限公司 | 航线生成方法、设备和系统、无人作业系统及存储介质 |
CN111089592A (zh) * | 2019-12-13 | 2020-05-01 | 天津大学 | 一种离散曲面中测地线计算方法 |
CN111736602A (zh) * | 2020-06-16 | 2020-10-02 | 国网上海市电力公司 | 一种改进的轮式机器人路径规划方法 |
Also Published As
Publication number | Publication date |
---|---|
CN108180918B (zh) | 2021-04-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108180918A (zh) | 一种点云测地路径正向跟踪生成方法及装置 | |
CN103714577B (zh) | 一种适用于带纹理模型的三维模型简化方法 | |
CN102903145B (zh) | 植物群体形态结构三维重建方法 | |
CN105654483A (zh) | 三维点云全自动配准方法 | |
CN109325510B (zh) | 一种基于网格统计的图像特征点匹配方法 | |
Mangasarian et al. | Nonlinear knowledge in kernel approximation | |
CN105630905A (zh) | 一种基于散乱点云数据的分层式压缩方法及装置 | |
CN107729648A (zh) | 一种基于Shepard插值的曲线纤维复合结构设计瀑布型多级优化方法 | |
CN109948002A (zh) | 基于平衡kd树的非结构网格最近壁面距离求解方法 | |
CN105096159A (zh) | 一种区域售电量预测方法及装置 | |
Song et al. | Fuzzy C-means clustering analysis based on quantum particle swarm optimization algorithm for the grouping of rock discontinuity sets | |
CN109683552A (zh) | 一种基面曲线导向的复杂点云模型上的数控加工路径生成方法 | |
CN109285219B (zh) | 一种基于dem的网格型水文模型网格演算次序编码方法 | |
CN104317886A (zh) | 断层约束下网格节点插值时近邻条件数据点的搜索选取方法 | |
CN109685841A (zh) | 三维模型与点云的配准方法及系统 | |
CN108595631A (zh) | 基于图谱理论的三维cad模型双层检索方法 | |
CN112241676A (zh) | 一种地形杂物自动识别的方法 | |
CN105653881A (zh) | 基于多密度层次的流场可视化方法 | |
CN104731885A (zh) | 一种基于层次-语义的多尺度空间数据拓扑关系维护方法 | |
CN113987610A (zh) | 基于网格映射的不同分辨率服装曲面网格的匹配方法 | |
CN104036549B (zh) | 一种合轴分枝树木形态结构三维可视化模拟方法 | |
CN107492129A (zh) | 基于素描表示和结构化聚类的非凸压缩感知优化重构方法 | |
CN110322415A (zh) | 基于点云的高精度表面三维重构方法 | |
CN107463528A (zh) | 基于ks检验的高斯混合模型分裂与合并算法 | |
CN107063218A (zh) | 数据处理方法和装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |