留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

面向Hilbert八叉树的邻近格元计算算法

吴宇豪 曹雪峰

吴宇豪, 曹雪峰. 面向Hilbert八叉树的邻近格元计算算法[J]. 武汉大学学报 ( 信息科学版), 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
引用本文: 吴宇豪, 曹雪峰. 面向Hilbert八叉树的邻近格元计算算法[J]. 武汉大学学报 ( 信息科学版), 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
WU Yuhao, CAO Xuefeng. A Hilbert-Octree Neighbor Grid Element Calculation Algorithm[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
Citation: WU Yuhao, CAO Xuefeng. A Hilbert-Octree Neighbor Grid Element Calculation Algorithm[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099

面向Hilbert八叉树的邻近格元计算算法

doi: 10.13203/j.whugis20200099
基金项目: 

国家自然科学基金 41401465

国家自然科学基金 41371384

详细信息
    作者简介:

    吴宇豪,硕士生,主要从事全球离散网格研究。514536575@qq.com

    通讯作者: 曹雪峰,博士,副教授。CAO_Xue_Feng@163.com
  • 中图分类号: P208

A Hilbert-Octree Neighbor Grid Element Calculation Algorithm

Funds: 

The National Natural Science Foundation of China 41401465

The National Natural Science Foundation of China 41371384

More Information
    Author Bio:

    WU Yuhao, postgraduate, specializes in the discrete global grid system. E-mail: 514536575@qq.com

    Corresponding author: CAO Xuefeng, PhD, associate professor. E‐mail: CAO_Xue_Feng@163.com
  • 摘要: 为提高线性八叉树邻近格元计算效率,利用Hilbert码标记格元,提出一种邻近格元Hilbert码快速计算方法。以Hilbert基元曲线为基础,引入状态向量的概念以记录Hilbert曲线对同属于一个父格元的所有子格元的填充顺序,从而建立状态向量的层级演进与退化函数,得到状态向量在m阶与m+1阶曲线中的层级映射关系,最终利用状态向量及其层级演进与退化函数实现邻近格元Hilbert码的计算。结果表明,所提算法计算结果正确;状态向量计算速度随层级提高而降低,在第20层级上1 ms内可完成4 201个格元的计算,对后续邻近格元计算影响较小;在指定层级上同等数量的邻近格元计算中,该算法的速度明显优于现有Morton码转换算法,在第15层级上百万级规模的邻近格元计算中,该算法的速度约为现有Morton码转换算法的2.1~2.4倍;在不同层级的百万级规模邻近格元计算中,该算法计算速度相比现有Morton码转换算法的提升倍数随层级提高而增大,在第20层级上该算法的效率提升达到2.6倍。
  • 图  1  一层级、二层级格元及其Hilbert码

    Figure  1.  First and Second Level Grid Elements and Their Hilbert Codes

    图  2  一次剖分得到的8个子格元

    Figure  2.  8 Sub-Grid Elements Obtained by One Split

    图  3  Hilbert基元曲线

    Figure  3.  Hilbert Primitive Curve

    图  4  状态向量演进

    Figure  4.  State Vector Evolution

    图  5  状态向量计算流程

    Figure  5.  Process of State Vector Calculation

    图  6  面邻近格元Hilbert码计算流程

    Figure  6.  Calculation Process of the Hilbert Code of the Neighboring Grid Element

    图  7  计算R方向上面邻近格元Hilbert码

    Figure  7.  Hilbert Code Calculation of Adjacent Grid Elements in the R Direction

    图  8  角邻近格元Hilbert码计算

    Figure  8.  Hilbert Code Calculation of Corner Adjacent Grid Element

    图  9  状态向量计算耗时与计算效率

    Figure  9.  Time-Consuming and Computational Efficiency of State Vector Calculation

    图  10  查询时间随格元个数变化

    Figure  10.  Query Time Varies with the Number of Grid Elements

    图  11  计算时间对比

    Figure  11.  Comparison of Calculation Time

    表  1  格元各方向面邻近格元位置

    Table  1.   Location of the Adjacent Grid Element in Each Direction of the Grid Element

    中心格元 邻近方向D
    T B L R I O
    LOB LOT LOT* ROB* ROB LIB LIB*
    LOT LOB* LOB ROT* ROT LIT LIT*
    ROT ROB* ROB LOT LOT* RIT RIT*
    ROB ROT ROT* LOB LOB* RIB RIB*
    RIB RIT RIT* LIB LIB* ROB* ROB
    RIT RIB* RIB LIT LIT* ROT* ROT
    LIT LIB* LIB RIT* RIT LOT* LOT
    LIB LIT LIT* RIB* RIB LOB* LOB
    下载: 导出CSV

    表  2  Hilbert码差值之和

    Table  2.   Sum of Hilbert Code Difference

    层级m 差值之和GVm 层级m 差值之和GVm
    1 0 6 0
    2 0 7 0
    3 0 8 0
    4 0 9 0
    5 0 10 0
    下载: 导出CSV

    表  3  不同格元数计算总时间统计

    Table  3.   Statistics of Total Time for Different Grid Elements

    格元数/106 计算总时间 提升倍数
    本文算法/μs 转换算法[13]/μs
    1 94.6 223.0 2.4
    2 183.5 425.7 2.3
    3 284.9 625.3 2.2
    4 367.9 845.6 2.3
    5 500.9 1 051.4 2.1
    6 536.5 1 251.4 2.3
    7 651.1 1 458.8 2.2
    8 753.3 1 660.8 2.2
    下载: 导出CSV

    表  4  算法效率提升倍数

    Table  4.   Improvement Multiple of Algorithm Efficiency

    层级 计算平均时间 提升倍数
    本文算法/μs 转换算法[13]/μs
    1 99.1 150.1 1.6
    2 99.0 152.4 1.7
    3 99.1 154.9 1.7
    4 102.4 157.2 1.6
    5 99.8 160.3 1.7
    6 97.1 168.4 1.9
    7 94.6 173.9 2.0
    8 98.1 180.5 2.0
    9 98.4 197.4 2.2
    10 97.9 197.2 2.2
    11 101.5 202.6 2.2
    12 107.6 205.0 2.1
    13 103.5 211.4 2.2
    14 102.5 219.1 2.3
    15 104.3 221.8 2.3
    16 101.7 230.0 2.4
    17 103.9 235.7 2.5
    18 105.8 256.7 2.6
    19 106.4 250.6 2.5
    20 104.9 258.6 2.6
    下载: 导出CSV
  • [1] 周果, 朱登明, 王兆其. 面向大规模场景的多片元效果高效绘制[J]. 计算机学报, 2017, 40(11): 2606-2618 doi:  10.11897/SP.J.1016.2017.02606

    Zhou Guo, Zhu Dengming, Wang Zhaoqi. Efficient Rendering of Multi-fragment Effects for Massive Scene[J]. Chinese Journal of Computers, 2017, 40 (11): 2606-2618 doi:  10.11897/SP.J.1016.2017.02606
    [2] Dutré J. Out-of-Core Construction of Sparse Voxel Octrees[J]. Computer Graphics Forum: Journal of the European Association for Computer Graphics, 2014, DOI: 10.1145/2492045.2492048
    [3] 吕梦雅, 刘丁, 唐勇, 等. 真实感海下光照效果实时绘制[J]. 小型微型计算机系统, 2018, 39(10): 2326-2329 doi:  10.3969/j.issn.1000-1220.2018.10.034

    Lü Mengya, Liu Ding, Tang Yong, et al. Improved Fast Rendering Algorithm for Large-Scale Underwater Illumination Effect[J]. Journal of Chinese Computer Systems, 2018, 39(10): 2326-2329 doi:  10.3969/j.issn.1000-1220.2018.10.034
    [4] Bi L, Zhao H, Jia M T. Database-Oriented Storage Based on LMDB and Linear Octree for Massive Block Model[J]. Transactions of Nonferrous Metals Society of China, 2016, 26(9): 2462-2468 doi:  10.1016/S1003-6326(16)64377-7
    [5] Moon B, Jagadish H V, Faloutsos C, et al. Analysis of the Clustering Properties of the Hilbert SpaceFilling Curve[J]. IEEE Trans Knowl Data Eng, 2001, 13(1): 124-141 doi:  10.1109/69.908985
    [6] 万刚, 曹雪峰. 地理空间信息网格的历史演变与思考[J]. 测绘学报, 2016, 45(S1): 15-22 doi:  10.11947/j.AGCS.2016.F002

    Wan Gang, Cao Xuefeng. The Historical Evolution and Reflection of Geospatial Information Grid[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45 (S1): 15-22 doi:  10.11947/j.AGCS.2016.F002
    [7] Guan X F, van Oosterom P, Cheng B. A Parallel N-Dimensional Space-Filling Curve Library and Its Application in Massive Point Cloud Management [J]. ISPRS International Journal of Geo-Information, 2018, 7(8): 327 doi:  10.3390/ijgi7080327
    [8] 曹布阳, 冯华森, 梁峻浩, 等. 利用Hilbert曲线与Cassandra技术实现时空大数据存储与索引[J]. 武汉大学学报·信息科学版, 2021, 46(5): 620-629 doi:  10.13203/j.whugis20200367

    Cao Buyang, Feng Huasen, Liang Junhao, et al. Hilbert Curve and Cassandra Based Indexing and Storing Approach for Large-Scale Spatiotemporal Data[J]. Geomatics and Information Science of Wuhan University, 2021, 46(5): 620-629 doi:  10.13203/j.whugis20200367
    [9] Mekuria R, Blom K, Cesar P. Design, Implementation, and Evaluation of a Point Cloud Codec for Tele-Immersive Video[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2017, 27(4): 828-842 doi:  10.1109/TCSVT.2016.2543039
    [10] Choi M G, Ju E, Chang J W, et al. Linkless Octree Using Multi-Level Perfect Hashing[J]. Computer Graphics Forum, 2009, 28(7): 1773-1780 doi:  10.1111/j.1467-8659.2009.01554.x
    [11] Meijers M, van Oosterom P. Clustering and Indexing Historic Vessel Movement Data with Space Filling Curves[J]. The International Archives of the Photogrammetry, Remote Sensing and Spatial Infor⁃ mation Sciences, 2018, XLII-4: 417-424
    [12] van Oosterom P, Martinez-Rubi O, Ivanova M, et al. Massive Point Cloud Data Management: Design, Implementation and Execution of a Point Cloud Benchmark [J]. Computers & Graphics, 2015, 49: 92-125 doi:  10.1016/j.cag.2015.01.007
    [13] Chen H L, Chang Y I. Neighbor-Finding Based on Space-Filling Curves [J]. Information Systems, 2005, 30(3): 205-226 doi:  10.1016/j.is.2003.12.002
    [14] 徐红波, 郝忠孝. 基于空间填充曲线网格划分的最近邻查询算法[J]. 计算机科学, 2010, 37 (1): 184-188 doi:  10.3969/j.issn.1002-137X.2010.01.043

    Xu Hongbo, Hao Zhongxiao. Nearest-Neighbor Query Algorithm Based on Grid Partition of SpaceFilling Curve[J]. Computer Science, 2010, 37(1): 184-188 doi:  10.3969/j.issn.1002-137X.2010.01.043
    [15] Sadeghi F, Kermani F, Rafsanjani M. Optimizing Image Steganography by Combining the GA and ICA [J]. ISC Int J Inf Secur, 2015, 7: 47-58 doi:  10.1158/1078-0432.CCR-04-2378
    [16] Huang W, Zhang W, Zhang D Y, et al. Elastic Spatial Query Processing in OpenStack Cloud Computing Environment for Time-Constraint Data Analysis [J]. ISPRS International Journal of Geo-Information, 2017, 6(3): 84 doi:  10.3390/ijgi6030084
    [17] 刘辉, 冷伟, 崔涛. 高维Hilbert曲线的编码与解码算法设计[J]. 数值计算与计算机应用, 2015, 36(1): 42-58 https://www.cnki.com.cn/Article/CJFDTOTAL-SZJS201501005.htm

    Liu Hui, Leng Wei, Cui Tao. Development of Encoding and Decoding Algorithms for High Dimensional Hilbert Curves[J]. Journal on Numerical Methods and Computer Applications, 2015, 36(1): 42-58 https://www.cnki.com.cn/Article/CJFDTOTAL-SZJS201501005.htm
    [18] Chen H L, Chang Y I. All-Nearest-Neighbors Finding Based on the Hilbert Curve[J]. Expert Systems with Applications, 2011, 38(6): 7462-7475 doi:  10.1016/j.eswa.2010.12.077
    [19] 戴晶, 吴明光, 郑培蓓, 等. 基于Hilbert曲线的STR索引改进算法[J]. 武汉大学学报·信息科学版, 2014, 39(7): 777-781 http://ch.whu.edu.cn/article/id/3019

    Dai Jing, Wu Mingguang, Zheng Peibei, et al. An Improved STR-Tree Spatial Index Algorithm Based on Hilbert-Curve[J]. Geomatics and Information Science of Wuhan University, 2014, 39(7): 777-781 http://ch.whu.edu.cn/article/id/3019
    [20] Shen J H, Lu C T, Jian M S. Neighbor-Index Method for Continuous Window Queries over Wireless Data Broadcast[J]. Applied Mechanics and Materials, 2013, 284: 3295-3299 doi:  10.4028/www.scientific.net/AMM.284-287.3295
    [21] Stommel M, Edelkamp S, Wiedemeyer T, et al. Fractal Approximate Nearest Neighbour Search in Log-Log Time[C]//The British Machine Vision Conference 2013, Bristol, UK, 2013
    [22] Liu X, Schrack G F. An Algorithm for Encoding and Decoding the 3-D Hilbert Order[J]. IEEE Transactions on Image Processing, 1997, 6(9): 1333-1337 doi:  10.1109/83.623197
    [23] Campbell P M, Flaherty J E, Gervasio L G. Dynamic Octree Load Balancing Using Space-Filling Curves[J]. Theory of Computing Systems, 2010, DOI: 10.1007/s00224-010-9306-3
    [24] Butz A R. Alternative Algorithm for Hilbert's Space-Filling Curve[J]. IEEE Transactions on Computers, 1971, C-20(4): 424-426 doi:  10.1109/T-C.1971.223258
  • [1] 徐光晶, 周坚鑫, 舒晴.  四叉树分解在海空重力测网交叉点搜索中的应用 . 武汉大学学报 ( 信息科学版), 2022, 47(11): 1847-1853. doi: 10.13203/j.whugis20200404
    [2] 吴宇豪, 曹雪峰.  虚拟战场环境时空数据的Hilbert码索引方法 . 武汉大学学报 ( 信息科学版), 2020, 45(9): 1403-1411. doi: 10.13203/j.whugis20190394
    [3] 牛磊, 宋宜全, 张宏敏, 侯绍洋.  一种针对室内疏散的集成Hilbert曲线的R*树空间索引 . 武汉大学学报 ( 信息科学版), 2018, 43(9): 1416-1421. doi: 10.13203/j.whugis20160352
    [4] 赵龙飞, 赵学胜, 朱思坤, 付瑞全.  一种球面退化四叉树格网的多层次邻近搜索算法 . 武汉大学学报 ( 信息科学版), 2018, 43(4): 529-535. doi: 10.13203/j.whugis20150611
    [5] 王金鑫, 禄丰年, 郭同德, 陈 杰.  球体大圆弧QTM八叉树剖分 . 武汉大学学报 ( 信息科学版), 2013, 38(3): 344-348.
    [6] 球体大圆弧QTM八叉树剖分 . 武汉大学学报 ( 信息科学版), 2013, 38(3): 344-.
    [7] 罗广祥, 刘苗, 樊鸿宇, 杨芳.  全球等面积四叉树离散格网建模与编码体系研究 . 武汉大学学报 ( 信息科学版), 2012, 37(10): 1252-1255.
    [8] 汪柱, 李树军, 张立华, 李宁.  基于航路二叉树的航线自动生成方法 . 武汉大学学报 ( 信息科学版), 2010, 35(4): 407-410.
    [9] 赵学胜, 崔马军, 李昂, 张美娟.  球面退化四叉树格网单元的邻近搜索算法 . 武汉大学学报 ( 信息科学版), 2009, 34(4): 479-482.
    [10] 喻永平, 陈晓勇, 刘经南, 都洁.  状态扩展元胞自动机模型在时空数据挖掘中的应用 . 武汉大学学报 ( 信息科学版), 2008, 33(6): 592-595.
    [11] 吴芳, 芮国胜.  基于四叉树和纠错编码的数字图像水印算法 . 武汉大学学报 ( 信息科学版), 2007, 32(3): 208-211.
    [12] 王永杰, 孟令奎, 赵春宇.  基于Hilbert空间排列码的海量空间数据划分算法研究 . 武汉大学学报 ( 信息科学版), 2007, 32(7): 650-653.
    [13] 吴波, 张良培, 李平湘.  基于支撑向量机概率输出的高光谱影像混合像元分解 . 武汉大学学报 ( 信息科学版), 2006, 31(1): 51-54.
    [14] 乔朝飞, 赵仁亮, 陈军, 陈云浩.  基于Voronoi内邻近的等高线树生成法 . 武汉大学学报 ( 信息科学版), 2005, 30(9): 801-804.
    [15] 白建军, 赵学胜, 陈军.  基于线性四叉树的全球离散格网索引 . 武汉大学学报 ( 信息科学版), 2005, 30(9): 805-808.
    [16] 徐毓, 李群, 周焰, 杨瑞娟.  按分量加权的探测目标状态线性融合 . 武汉大学学报 ( 信息科学版), 2002, 27(4): 420-423.
    [17] 盛业华, 唐宏, 杜培军.  线性四叉树快速动态编码及其实现 . 武汉大学学报 ( 信息科学版), 2000, 25(4): 324-328.
    [18] 李清泉, 李德仁.  八叉树的三维行程编码 . 武汉大学学报 ( 信息科学版), 1997, 22(2): 102-106.
    [19] 谈国新, 林宗坚.  二值图像的紧凑二叉树表示及其编码方法 . 武汉大学学报 ( 信息科学版), 1995, 20(3): 219-223,281.
    [20] 邓朝晖, 贾华.  直接表达区域的四叉树链式编码 . 武汉大学学报 ( 信息科学版), 1995, 20(3): 224-227.
  • 加载中
图(11) / 表(4)
计量
  • 文章访问数:  282
  • HTML全文浏览量:  84
  • PDF下载量:  32
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-03-17
  • 刊出日期:  2022-04-05

面向Hilbert八叉树的邻近格元计算算法

doi: 10.13203/j.whugis20200099
    基金项目:

    国家自然科学基金 41401465

    国家自然科学基金 41371384

    作者简介:

    吴宇豪,硕士生,主要从事全球离散网格研究。514536575@qq.com

    通讯作者: 曹雪峰,博士,副教授。CAO_Xue_Feng@163.com
  • 中图分类号: P208

摘要: 为提高线性八叉树邻近格元计算效率,利用Hilbert码标记格元,提出一种邻近格元Hilbert码快速计算方法。以Hilbert基元曲线为基础,引入状态向量的概念以记录Hilbert曲线对同属于一个父格元的所有子格元的填充顺序,从而建立状态向量的层级演进与退化函数,得到状态向量在m阶与m+1阶曲线中的层级映射关系,最终利用状态向量及其层级演进与退化函数实现邻近格元Hilbert码的计算。结果表明,所提算法计算结果正确;状态向量计算速度随层级提高而降低,在第20层级上1 ms内可完成4 201个格元的计算,对后续邻近格元计算影响较小;在指定层级上同等数量的邻近格元计算中,该算法的速度明显优于现有Morton码转换算法,在第15层级上百万级规模的邻近格元计算中,该算法的速度约为现有Morton码转换算法的2.1~2.4倍;在不同层级的百万级规模邻近格元计算中,该算法计算速度相比现有Morton码转换算法的提升倍数随层级提高而增大,在第20层级上该算法的效率提升达到2.6倍。

English Abstract

吴宇豪, 曹雪峰. 面向Hilbert八叉树的邻近格元计算算法[J]. 武汉大学学报 ( 信息科学版), 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
引用本文: 吴宇豪, 曹雪峰. 面向Hilbert八叉树的邻近格元计算算法[J]. 武汉大学学报 ( 信息科学版), 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
WU Yuhao, CAO Xuefeng. A Hilbert-Octree Neighbor Grid Element Calculation Algorithm[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
Citation: WU Yuhao, CAO Xuefeng. A Hilbert-Octree Neighbor Grid Element Calculation Algorithm[J]. Geomatics and Information Science of Wuhan University, 2022, 47(4): 613-622. doi: 10.13203/j.whugis20200099
  • 线性八叉树只记录叶结点的编码,不记录中间结点的编码及层次关系,相比普通四叉树、八叉树大大节省了存储空间,且包含层次信息,因此在实际工作中得到了广泛的应用[1-4]。在众多线性八叉树叶节点编码方法中,三维Hilbert曲线映射得到的Hilbert码聚簇效果最优[5],使得空间数据对应的格元分布更加紧凑,因而受到许多研究人员关注[6-9]

    三维Hilbert曲线将连续三维空间划分为规则的格元,且不遗漏、不交叉地穿行所有格元,线性八叉树叶节点Hilbert码即为叶节点所占据格元在Hilbert曲线上的排列次序。八叉树邻域计算是空间聚类、索引、查询等空间操作的基础[10]。Hilbert曲线编码的线性八叉树邻域计算过程一般分为两步[11-12]:(1)给定中心格元的Hilbert码,计算其邻近格元Hilbert码;(2)遍历八叉树,寻找邻近格元Hilbert码对应的叶节点。可见,邻近格元Hilbert码的计算效率是实现高效八叉树邻域计算的关键。

    典型的邻近格元Hilbert码计算是基于Morton码与Hilbert码之间的空间变换关系实现的[13-14],首先通过空间变换关系将Hilbert码转换为Morton码;然后使用二进制位操作进行邻近Morton码计算过程;最后将邻近格元Morton码转换为Hilbert码。文献[13]的算法已被广泛应用于图像数据隐写[15]、时空数据查询[16]等领域。上述算法通过将Hilbert码转换为Morton码实现邻近Hilbert码计算,但是这种转换算法中存在以下两个缺点:(1)构造三维Hilbert码的过程较为复杂,Hilbert码层级拆解转换时需要进行多次分支判断和嵌套循环运算,计算复杂度为O3。(2)当中心格元层级为m时,两种转换算法均需要进行两次步长为m的循环转换过程,因此转换算法的总时间复杂度为O3×2m=O6m[17]。随着格元层级的增高,转换算法计算指定邻近方向上邻近Hilbert码所需的时间将会线性增长,效率随之线性降低。

    为了提升八叉树邻域计算效率,不需要将Hilbert码转换为其他形式的坐标或空间排列码,仅使用Hilbert曲线自身性质实现邻近格元Hilbert码计算,文献[18]在研究两个二维点数据集之间的全近邻问题时,通过记录中心格元在各层级中的相对位置以及Hilbert曲线的填充顺序实现了邻近格元Hilbert码计算,相比转换算法在效率上有较大提升;文献[18]提出的二维Hilbert码计算高效,已被应用于空间数据索引[19]、无线广播数据窗口查询[20]、分形空间邻近搜索[21]等领域。但是二维与三维Hilbert曲线的构造方法全然不同,层级之间的演进关系也因此不同,且各个叶节点之间并不同属于一个层级,无法直接将文献[18]所提算法套用至本文所研究的八叉树邻近计算中。

    为提高线性八叉树中邻近格元Hilbert码的计算效率,本文在研究三维Hilbert基元曲线的基础上,引入状态向量记录曲线在格元间的填充顺序,建立层级演进与退化函数确定状态向量在m阶与m+1阶曲线中的映射关系,实现了一种Hilbert曲线编码的线性八叉树邻近快速计算方法,并设计了相关实验对算法的正确性以及效率进行分析。

    • Hilbert曲线先将n维空间规则划分为无缝无叠的2mn个超立方体网格单元(称为格元),再利用一维曲线不遗漏、不交叉地穿行所有格元,实现了n维空间Rn与一维空间R1的一一映射。三维Hilbert曲线实现了三维空间至一维空间的一一映射,其对边长为T的立方体空间作嵌套八叉剖分,第一次剖分得到8个子格元,经过m次剖分后得到8m个子格元,每个子格元边长为T/2mm层级格元对应于m阶曲线。层级m的取值范围为0,+),当m=0时,格元为根格元。

      将曲线在各层级格元中的填充顺序ε相互衔接作为m层级格元的Hilbert码,表示为h=ε1ε2εm,如图 1所示。因此,Hilbert码涵盖了两个重要层级剖分性质:(1)Hilbert码各个位上的数值为曲线在各层级格元中的填充顺序;(2)两个同属于一个父格元的子格元,其Hilbert码除最后一位外前缀完全相同。

      图  1  一层级、二层级格元及其Hilbert码

      Figure 1.  First and Second Level Grid Elements and Their Hilbert Codes

      将层级为m、Hilbert码为h的格元记作vhm。根据上述曲线与格元性质,作以下计算、概念约定:(1)与格元vhm层级差为k的父格元记作Pkv=vhdiv8m-k,其中hdiv8=h/8k,表示除去Hilbert码h的最后k位数。本文中的父格元指与格元vhm层级差为1的父格元,简记为Pv。(2)格元vhm进行一次剖分得到的8个子格元中,填充顺序为i的子格元Cv,i=vhim+1Cv,i称作格元v的第i个子格元。(3)格元vhmPv剖分得到的8个子格元中的填充顺序Iv=hmod8,其中mod为取余运算,Iv即为Hilbert码h的最后1位数。(4)为区分不同的邻近关系,本文将同一层级格元间的邻近关系分为兄弟邻近以及堂兄弟邻近两种关系。由同一格元经过一次剖分得到的称为兄弟邻近格元,其余称为堂兄弟邻近格元。例如格元v302v312v322v362为兄弟邻近关系,与v472为堂兄弟邻近关系。兄弟邻近格元Hilbert码除最后一位外前缀完全相同。

    • 为区分不同的空间方向D,将空间方向共分为上(top,T)、下(bottom,B)、左(left,L)、右(right,R)、内(in,I)、外(out,O)6个方向。格元进行一次剖分得到8个子格元,将其按照方向位置分别称为LOB、LOT、ROT、ROB、RIB、RIT、LIT、LIB,如图 2所示。线性八叉树邻近格元分为面相邻的面邻近、边相邻的边邻近、角相邻的角邻近3种,例如LOT为ROT的面邻近,LOB为ROT的边邻近,LOT为RIB的角邻近。在同层级下,线性八叉树任一格元的邻近格元集合是由6个面邻近、12个边邻近、8个角邻近等共26个邻近格元组成。本文主要研究的邻域计算问题就是快速计算任一格元的26个邻近格元Hilbert码。

      图  2  一次剖分得到的8个子格元

      Figure 2.  8 Sub-Grid Elements Obtained by One Split

    • 本文的研究目的是仅使用Hilbert曲线自身性质实现邻近格元Hilbert码计算。因此,了解掌握Hilbert曲线的填充规律是进行邻近Hilbert码计算的前提。首先对Hilbert曲线的填充规律进行分析,提出状态矩阵以及层级演进退化函数概念,然后将其应用于邻近格元Hilbert码的计算中。

    • Hilbert曲线的自相似与自复制特征决定了Hilbert曲线中只有24种形态的基元曲线,不管多复杂的Hilbert曲线都由这24种基元曲线递归派生而来[22],分别为φ1,φ2φ24,如图 3所示,每种形态基元曲线的区分标准是Hilbert曲线在由一个格元剖分得到的LOB、LOT、ROT、ROB、RIB、RIT、LIT、LIB中的填充状态。

      图  3  Hilbert基元曲线

      Figure 3.  Hilbert Primitive Curve

      为描述这24种形态基元曲线,引入状态矩阵S记录24种基元曲线,按照LOB、LOT、ROT、ROB、RIB、RIT、LIT、LIB的顺序,记录第k种基元曲线在填充至当前子格元时共经过了多少个子格元,记作状态向量skk1, 224,然后将24个行向量合成矩阵S,矩阵S的第u行对应状态向量su+1,对应基元曲线为φu+1

      针对一个格元,Hilbert曲线在其子格元间的填充状态为φk,Hilbert曲线在各个子格元再次剖分后得到的格元间的填充状态由φk决定,这种映射关系不随层级变化而发生改变[23]

      图 4所示,二阶曲线的填充状态为φ2,对应状态向量为s2;剖分后Hilbert曲线在各个子格元间的填充状态分别为φ3φ1φ2φ10φ3φ2φ5φ10,状态向量分别对应s3s1s2s10s3s2s5s10,则存在状态向量层级映射关系T(s2)s3,s1, s2,s10,s3,s2,s5,s10。这种状态向量映射关系并不会随层级变化而发生改变,即无论格元层级的高低,若其状态向量为s2,其子格元的状态向量均为s3s1s2s10s3s2s5s10

      图  4  状态向量演进

      Figure 4.  State Vector Evolution

      引入状态演进矩阵E记录24种状态向量与下一层级格元状态向量的映射关系。Euu0, 1,223对应映射关系Tsu+1,即:

      Tsu+1sE[u][0],sE[u][1],sE[u][2],sE[u][3],sE[u][4],sE[u][5],sE[u][6],sE[u][7],sE[u][8]
      S=016745230123456707432561034765214523016747032165076125436107432567452301612543074567012323016745216547036521034725430761432561070321654765470321430761252103476547652103670123452345670125610743
      E=721871118312103251024392369173418174191815115161552162019621206321117722172122222081228181102393109223913102910142245111924111191922122419122324211013202113242018241417181410175181523515202323211652316175471714417161414151841418841261911121941168201362015131316216132176812227822971692315162312151114241211241312

      根据E可知,若已知格元vhm的状态向量sk,可求得其填充顺序为i的子格元的状态向量sCv,i为:

      sCv,i=sTcsk,i=sE[k-1][i]

      式中,Tcsk,i为状态演进函数。使用状态演进函数对格元vhm进行一次演进,可得到格元的第i个子格元的Hilbert码以及状态向量。借助函数Tcsk,i,以格元vhm的Hilbert码h为输入,计算格元vhm状态向量的算法流程如图 5所示,其中hlength为Hilbert码的长度,hbiti为Hilbert码h的第i位。

      图  5  状态向量计算流程

      Figure 5.  Process of State Vector Calculation

      对函数Tcsk,i求逆,可得其逆函数为:

      sv=sTpsCv,i,i

      式中,sCv,i为格元vhm的第i个子格元的状态向量;TpsCv,i,i为状态退化函数。使用状态退化函数对子格元Cv,i进行一次退化,可得到格元vhm的Hilbert码以及状态向量。

    • §1.1中已经按照空间方向对一次剖分得到的8个子格元的位置进行区分,由此推导得到邻近方向D上各个子格元的同层级面邻近子格元位置,结果如表 1所示

      表 1  格元各方向面邻近格元位置

      Table 1.  Location of the Adjacent Grid Element in Each Direction of the Grid Element

      中心格元 邻近方向D
      T B L R I O
      LOB LOT LOT* ROB* ROB LIB LIB*
      LOT LOB* LOB ROT* ROT LIT LIT*
      ROT ROB* ROB LOT LOT* RIT RIT*
      ROB ROT ROT* LOB LOB* RIB RIB*
      RIB RIT RIT* LIB LIB* ROB* ROB
      RIT RIB* RIB LIT LIT* ROT* ROT
      LIT LIB* LIB RIT* RIT LOT* LOT
      LIB LIT LIT* RIB* RIB LOB* LOB

      因此,根据Hilbert曲线所涵盖的层级剖分性质,可在已知中心格元状态向量的前提下,计算面邻近格元Hilbert码,步骤如下:(1)以中心格元为起点,层级退化搜索共同父格元。(2)以共同父格元为起点,层级演进计算得到面邻近格元的Hilbert码。层级退化演进过程中,利用状态向量判断当前格元所处的位置,而层级退化演进函数则用于计算不同层级格元的状态向量。每一次退化都将中心格元在D方向上的邻近格元位置存入数组Direction[ ]中。演进时,每一次选择的演进格元位置只需从Direction[ ]的末尾开始往前取值即可。此外,还需在各次层级退化演进中检查八叉树中是否存在当前格元。面邻近格元Hilbert码计算算法具体流程如图 6所示。

      图  6  面邻近格元Hilbert码计算流程

      Figure 6.  Calculation Process of the Hilbert Code of the Neighboring Grid Element

      以计算中心格元v0423在R方向上面邻近格元Hilbert码为例。详细阐述本文面邻近格元Hilbert码计算过程,步骤如下:(1)以格元v0423为起点退化一层级,计算得到其父格元v042的状态向量为s2,并推得v0423属于ROT。查询表 1可知,ROT在R方向上的面邻近格元为LOT*,LOT*与ROT不属于同一个父格元,LOT*与ROT之间属于堂兄弟面邻近关系,无法使用格元v042演进一层级的方法计算LOT*的Hilbert码。(2)继续以格元v042为起点退化一层级,计算得到其父格元v01的状态向量为s3。在状态向量s3中,填充顺序为404̲2的子格元为ROT,v042在R方向上的面邻近格元仍然与v042不同属于一个父单元。(3)继续以格元v01为起点退化1层级,计算得到根格元v0的状态向量为s2。在状态向量s2中,填充顺序为00̲42的子格元为LOB,LOB在R方向上的子格元是ROB,LOB与ROB同属于一个父格元,LOB与ROB之间属于兄弟面邻近关系,可使用演进1层级的方法计算LOB的Hilbert码。(4)在状态向量s2中,ROB的填充顺序为3。以根格元v0为起点,向下演进1层级,得到根格元v0填充顺序为3的子格元Cv0,3=v31v31的状态向量sv31=sC(v0,3)=sTcs2,3=sE[1][3]=s10。(5)在状态向量s10中,LOT的填充顺序为1。以格元v31为起点,向下演进1层级,计算格元v31填充顺序为1的子格元Cv31,1=v312v312的状态向量sv312=sC(v31,1)=s13。(6)在状态向量s13中,LOT的填充顺序为1。以格元v312为起点,向下演进一层级,计算格元v312填充顺序为1的子格元Cv312,1=v3113v3113的状态向量sv3113=sC(v312,1)=s10。检测八叉树中不存在格元v3113,则停止演进,返回Hilbert码31,最终计算得到格元v0423在R方向上面邻近格元为v312。整体计算过程如图 7所示,黄色格元为中心格元v0423

      图  7  计算R方向上面邻近格元Hilbert码

      Figure 7.  Hilbert Code Calculation of Adjacent Grid Elements in the R Direction

      分析上述面邻近格元Hilbert码计算算法流程可知,面邻近格元Hilbert码计算算法的演进与退化过程都最多为k次循环(k为中心格元vhm与共同父格元之间的层级差),因此计算同层级邻近格元Hilbert码时,面邻近格元Hilbert码计算算法的时间复杂度为O2k

      相比现有的转换算法[13],本文算法的优势如下:(1)定义了状态向量与层级演进退化函数,算法中不需要将Hilbert码转换为其他形式的坐标或空间排列码,避免了复杂Hilbert码层级拆解转换运算,不存在多种分支判断计算或嵌套循环计算,算法计算复杂度较小。(2)算法只具有单层循环,共需要进行两次步长为k的循环计算,算法的时间复杂度为O2k。相比于转换算法的时间复杂度O6m,层级差k必定不大于中心格元层级m,故本文算法的循环计算次数必定不大于转换算法,且本文算法的计算复杂度更小,整体算法效率更高。

    • 角邻近格元Hilbert码计算过程可以由3次面邻近格元Hilbert码计算组成,如图 8所示,首先以黄色格元为中心格元,计算B方向上的面邻近无色格元Hilbert码,然后再以无色格元为中心格元计算L方向上的面邻近灰色格元Hilbert码,最后再以灰色格元为中心格元计算I方向上的面邻近绿色格元Hilbert码,得到黄色格元角邻近格元中的绿色格元Hilbert码。

      图  8  角邻近格元Hilbert码计算

      Figure 8.  Hilbert Code Calculation of Corner Adjacent Grid Element

    • 为了验证算法正确性、性能与效率,设计了邻近格元计算正确性、状态向量计算效率、Hilbert码计算效率3组实验进行测试。实验硬件环境为CPU Intel Core i7-7700K(双核4.2 GHz),内存为64 GB,硬盘速度为7 200 r/min;软件环境为Visual Studio 2015,Release版本,X64,C++。

    • Hilbert曲线具有空间填充曲线的分形特征,具有自相似与自复制特性,高层级曲线由低层级曲线平移、求反与交换推导得出。对比不同算法的计算结果,若低层级结果一致,由分形特征可知,其高层级结果也一致[24]。因此,使用本文算法与文献[13]算法分别计算不同层级中心格元邻近格元的Hilbert码并进行对比验证。

      将层级为m的所有网格单元集合记为Vm,将文献[13]算法求得的Hilbert码记作I(v),将本文算法生成的Hilbert码记作H(v)。对于同一邻近格元,定义两种算法计算Hilbert码差值ΔCode(v)为:

      ΔCodev=absIv-Hv

      定义m层级单元的Hilbert码差值之和GVm为:

      GVm=vVmΔCodev

      计算了1~10层的Hilbert码差值之和,结果如表 2所示。由表 2可知,各层级GVm均为0,证明本文算法与文献[13]算法生成的Hilbert码结果一致。10层级之内本文算法结果正确,由分形特征推知更高层级结果仍正确。

      表 2  Hilbert码差值之和

      Table 2.  Sum of Hilbert Code Difference

      层级m 差值之和GVm 层级m 差值之和GVm
      1 0 6 0
      2 0 7 0
      3 0 8 0
      4 0 9 0
      5 0 10 0
    • 相比转换算法,本文算法执行之前需要已知中心格元的状态向量。为探究中心格元状态向量计算步骤的执行效率是否会对后续邻近格元Hilbert码计算产生影响,设计实验进行测试。在1~20层级的所有格元中,随机挑选1×106个格元(存在重复格元)进行状态向量计算,记录其计算总耗时,并计算各个层级内的计算效率,即1 ms时间内可进行计算的次数。图 9为1~20层级状态向量计算耗时与计算效率,展现了计算总时间与计算效率随层级变化的趋势。

      图  9  状态向量计算耗时与计算效率

      Figure 9.  Time-Consuming and Computational Efficiency of State Vector Calculation

      图 9可知,计算平均时间对应的曲线随着层级的增加而逐步上升,表明进行一次状态向量计算消耗的时间与层级呈正相关,其中第20层级进行一次计算的时间为238 μs。计算效率对应曲线随层级的增加而逐步下降,表明状态向量计算效率与层级呈负相关,其中第20层级1 ms内仍可对4 201个格元进行计算。

    • 为分析本文算法计算邻近格元Hilbert码的效率性能,将本文算法与现有的Morton码转换算法[13]进行对比实验。实验1:对比相同层级下,两种算法计算不同次数所需的总耗时;实验2:对比相同计算次数下,两种算法在不同层级下计算所需的总耗时。

    • 在层级15情况下,随机选取当前层级的n个格元,计算其状态序列,并采用两种算法计算26个邻近格元Hilbert码。实验过程中,n的取值分别为1×106,2×1068×106,统计各数量格元邻近格元Hilbert码的总耗时,结果如表 3所示,查询时间随格元个数变化如图 10所示。由图 10可知,两种算法对应曲线随着格元个数的增加均呈现近似线性上升,表明两种算法的计算所需时间均与格元个数呈正相关;对比图 10中相同格元数下两种算法所需的计算时间,本文算法的计算时间均要少于转换算法,计算速度为转换算法的2.1~2.4倍。

      表 3  不同格元数计算总时间统计

      Table 3.  Statistics of Total Time for Different Grid Elements

      格元数/106 计算总时间 提升倍数
      本文算法/μs 转换算法[13]/μs
      1 94.6 223.0 2.4
      2 183.5 425.7 2.3
      3 284.9 625.3 2.2
      4 367.9 845.6 2.3
      5 500.9 1 051.4 2.1
      6 536.5 1 251.4 2.3
      7 651.1 1 458.8 2.2
      8 753.3 1 660.8 2.2

      图  10  查询时间随格元个数变化

      Figure 10.  Query Time Varies with the Number of Grid Elements

    • 在层级1~20下,随机选取当前层级中的1×106个格元作为样本集,计算样本集内所有格元状态序列,并分别采用两种算法计算其26个邻近格元的Hilbert码,统计两种算法计算Hilbert码总耗时后,对一次邻域计算的平均时间进行计算,结果如图 11所示。将Morton码转换算法的计算时间除以本文算法的计算时间,作为本文算法效率提升的倍数,各层级结果如表 4所示。

      图  11  计算时间对比

      Figure 11.  Comparison of Calculation Time

      表 4  算法效率提升倍数

      Table 4.  Improvement Multiple of Algorithm Efficiency

      层级 计算平均时间 提升倍数
      本文算法/μs 转换算法[13]/μs
      1 99.1 150.1 1.6
      2 99.0 152.4 1.7
      3 99.1 154.9 1.7
      4 102.4 157.2 1.6
      5 99.8 160.3 1.7
      6 97.1 168.4 1.9
      7 94.6 173.9 2.0
      8 98.1 180.5 2.0
      9 98.4 197.4 2.2
      10 97.9 197.2 2.2
      11 101.5 202.6 2.2
      12 107.6 205.0 2.1
      13 103.5 211.4 2.2
      14 102.5 219.1 2.3
      15 104.3 221.8 2.3
      16 101.7 230.0 2.4
      17 103.9 235.7 2.5
      18 105.8 256.7 2.6
      19 106.4 250.6 2.5
      20 104.9 258.6 2.6

      图 11可知,转换算法对应曲线随层级增加呈现较为均匀的增长趋势,表明转换算法的邻近格元Hilbert码计算效率随层级的增加而线性降低;本文算法对应曲线较为平缓,随着层级的增加存在较为微小的增长,但基本处于100 μs以下,表明本文算法受层级影响较小,更加满足线性八叉树邻域计算需求。对比图 11中相同层级下两种算法所需的计算时间,本文算法所需时间均小于转换算法,且本文算法较转换算法的效率提升倍数随层级增加有扩大趋势,由表 4可知,在第20层级,本文算法效率相比于转换算法提升了2.6倍。

      由§3.2中实验结果可知,计算一个格元的状态向量所需时间的数量级为10-4 μs,进而可推得26个格元的状态向量所需时间的数量级为10-3 μs,而进行一次邻近Hilbert计算所需平均时间数量级为102 μs,对比两者可知,计算格元状态向量所需时间可在计算邻近Hilbert码时忽略不计。

      综上所述,在计算邻近格元Hilbert码的效率上,本文算法的效率要优于转换算法。

    • 为提高基于Hilbert曲线编码的线性八叉树邻近格元计算效率,本文基于Hilbert曲线的填充规律设计了一种邻近格元Hilbert码的快速计算算法,并详细阐述了该算法的流程。两组实验对比结果表明,本文算法在邻近格元Hilbert码计算的效率上优于现有的Morton码转换算法,可作为高效数据检索、三维拓扑关系确定与连通性判断等空间操作的基础。

参考文献 (24)

目录

    /

    返回文章
    返回