EN中文

地球信息科学学报 >

2014 , Vol. 16 >Issue 6: 846 - 851

DOI: https://doi.org/10.3724/SP.J.1047.2014.00846

基于状态转移矩阵的Hilbert码快速生成算法

  • 李绍俊 , 1, 2, 3 ,
  • 钟耳顺 1, 3 ,
  • 王少华 , 1, 3, * ,
  • 张珣 1, 2
展开
  • 1. 中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室,北京 100101
  • 2. 中国科学院大学,北京 100049
  • 3. 北京超图软件股份有限公司,北京 100015
*通讯作者：王少华(1983-),男,博士,研究方向为地理信息系统软件开发与应用。E-mail:

作者简介：李绍俊(1978-),男,博士生,研究方向为GIS软件技术。E-mail:

收稿日期: 2014-03-14

  要求修回日期: 2014-05-05

  网络出版日期: 2014-11-01

基金资助

国家科技支撑计划项目(2009AA12Z331、2011BAH06B03)

交通运输部科技项目(2012-364-X04-102)

收起

An Algorithm for Hilbert Ordering Code Based on State-Transition Matrix

  • LI Shaojun , 1, 2, 3 ,
  • ZHONG Ershun 1, 3 ,
  • WANG Shaohua , 1, 3, * ,
  • ZHANG Xun 1, 2
Expand
  • 1. State Key Laboratory of Resources and Environment Information System, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing 100101, China
  • 2. University of Chinese Academy of Sciences, Beijing 100039, China
  • 3. SuperMap Software Co., Ltd., Beijing 100015, China
*Corresponding author:WANG Shaohua, E-mail:

Received date: 2014-03-14

  Request revised date: 2014-05-05

  Online published: 2014-11-01

Copyright

《地球信息科学学报》编辑部 所有
Fold

摘要

空间填充曲线的空间排列码可实现多维空间到一维空间的线性映射,广泛应用于空间查询、空间索引、空间划分及影像编码等领域。Hilbert是一种优秀的空间填充曲线,具有非常好的空间聚集性。传统的Hilbert排列二进制循环位操作算法的算法复杂度为O(n2)。本文首先分析了Hilbert的分形自相似特性,推导并归纳出Hilbert状态转移矩阵,按位编码顺序定义了空间划分中的象限顺序,将Hilbert状态转移矩阵转换为C++中的数组运算,减少了Hilbert码计算过程中的嵌套循环及迭代处理,将算法复杂度降为O(n)。其次,采用位域共用体以数值计算替代了传统计算过程中的数值与字符串间类型转换,提高了Hilbert码生成算法的性能。最后,在C++环境下实现了Hilbert码快速生成算法的相关代码,并完成算法的正确性验证实验和性能对比实验。实验结果表明,本文提出的算法计算结果与二进制循环位算法的结果一致,在性能上本文算法与二进制循环位算法及空间层次分解算法相比有明显的优势。

关键词： 线性映射; 空间填充曲线; 状态转移矩阵; Hilbert排列码; QuickHilbertCode(QHC)

本文引用格式

李绍俊 , 钟耳顺 , 王少华 , 张珣 . 基于状态转移矩阵的Hilbert码快速生成算法[J]. 地球信息科学学报, 2014 , 16(6) : 846 -851 . DOI: 10.3724/SP.J.1047.2014.00846

Abstract

Spatial ordering based on space- filling curves is a linear mapping method from multi- dimensional space to one-dimensional space, which has been widely used in spatial querying, spatial indexing, space partition, image coding, and other related fields. The Hilbert curve is an excellent space-filling curve with remarkable spatial clustering properties. The traditional algorithm for Hilbert ordering code is based on binary circular bit manipulation, which has a complexity of O(n2). In this article, the Hilbert state- transition matrix is generated based on the fractal self- similarity, the spatial quadrants are redefined by the sequence of spatial partition, and the Hilbert state-transition matrix is translated into arrays in C++, all of which effectively reduce the nested loops and the iterations in the process of computing Hilbert code, thus decrease the complexity of algorithm to O(n). Also, the bit field union type is used to avoid the type conversion from number to string in the process of Hilbert code computation, which brings numerical calculation into full play and improves the performance of the Hilbert coding algorithm. Finally, the C++ code for Hilbert code generating algorithm is given, and various tests have been conducted to verify the correctness of the algorithm and to compare the performance with regard to the binary circular bit algorithm as well as the hierarchical space decomposition algorithm. Test results show that the algorithm demonstrates notable merits compared with the other two algorithms.

Key words： linear mapping; space- filling curves; Hilbert ordering code; state- transition matrix; Quick Hilbert Code (QHC)

1 引言

空间填充曲线的线性映射可将多维数据序列化为一维数据,便于多维数据在磁盘上的存储与索引[1],这种映射模式在影像压缩、行程编码形式的栅格数据表达、空间划分[2-6]、空间查询[7]、空间索引[8]、计算几何中的启发式搜索等领域应用广泛。空间排列码(Spatial Ordering Code)以连续整数与空间填充曲线实体集元素建立起一一对应的可逆对应关系,奠定了空间填充曲线相关算法的基础。目前,应用最广泛的空间排列方法包括Morton码、Gray码及Hilbert码。经大量应用验证,Hilbert空间排列码具有最好的空间聚集及空间连续性[7-10]。因此,本文提出一种基于状态转移矩阵的快速Hilbert码计算方法,减少了算法过程中的递归和迭代次数,提升了Hilbert码的计算效率,有利于Hilbert码的推广及深入应用。

2 Hilbert空间填充曲线及其空间排列码算法

2.1 Hilbert空间填充曲线

Hilbert空间填充曲线属于Peano曲线族,具有FASS(Space Filling, Self-Avoiding, Simple and Self-Similar)特征,Hilbert曲线可以不间断、不自交地通过栅格所有格元,且每一格元仅被通过一次[11-12]。Hilbert曲线对完整空间作逐级双向等分分割,第1次分割可得到4个区域,将每个区域称为格元,经过k次分割后可得到4k个子格元,每个子格元的边长为原完整空间边长的2-k。将这些格元记为 A i k ( i = 1,2 , … , 4 k ) ,则这些格元满足如下条件：
（1） k 级格元 A i k 完整覆盖k+1级的4个子格元;
A i k = A 4 i k + 1 ⋃ A 4 i + 1 k + 1 ⋃ A 4 i + 2 k + 1 ⋃ A 4 i + 3 k + 1 （1）
（2）同一级中下标为连续的格元 A i k 与 A i + 1 k 必相邻（有公共边）;
（3） k 级4个子格元的形状可以由k-1级格元的形状进行旋转及对调0和3子格元得到。
图1给出了k=1,2,3时的Hilbert曲线。

2.2 Hilbert空间排列码算法

二进制位操作的Hilbert空间排列码生成算法,由Faloutsos和Roseman在1989年提出并通过计算机实现[1],它基于格网行列坐标进行二进制位操作,需进行迭代循环,其算法时间复杂度为O(n2)。陆锋等在2001年提出基于空间层次分解的生成算法[13],王笋等在2006年提出基于曲线扫描矩阵的生成算法[14]。鉴此,本文提出了一种基于状态转移矩阵的Hilbert码快速生成算法,避免了算法中的嵌套循环,将Hilbert码生成算法的时间复杂度降为稳定的O(n);同时,通过引入位域结构体,避免了计算过程中的数值与字符串间的类型转换,进而提高了Hilbert码生成算法的性能。
2.2.1 二进制位操作的Hilbert码生成算法[1]
算法描述如下（层级：k;列坐标：x;行坐标：y）：
（1）将列坐标（x）和行坐标（y）表达为二进制字符串,并进行位交叉,转化为对应的Morton码;
（2）Morton码的二进制字符串以2位为单位分组,将分组结果中的编码′10′和编码′11′互换;
（3）自右侧第k个分组 G k 开始,逐个向右进行迭代,并作如下处理：若编码 G k = ′00′,将后续各组中的编码′11′和′01′互换;若 G k = ′11′,将后续各组中的编码′10′和′00′互换;
（4）将结果数组连起来,并转化为整型,即为k级（x, y）格元的Hilbert空间排列码。
2.2.2 空间层次分解的Hilbert码生成算法
陆锋等通过对Hilbert空间填充曲线特征的分析得到曲线共有4种子象限形态,并决定了状态转移向量的4种取值形式,设计实现了基于空间层次分解的Hilbert码生成算法,具体计算时采用栅格空间层次分解与格元状态转移向量相结合的方法,其中状态转移向量表明空间层次分解中子曲线的旋转与反射[13]。具体图示、算法描述及伪码请参阅文献[13]。

3 状态转移矩阵的Hilbert码快速生成算法

状态转移矩阵的Hilbert码快速生成算法(Quick Hilbert Coding,QHC)更适合于计算机编码实现,并可获得更高的运行效率。

3.1 QHC中的象限定义

为便于计算机算法实现中的位运算,QHC算法中对每一级空间的分割不采用笛卡尔平面直角坐标系的象限划分,而采用图2所示的四分象限规则,水平方向为记为 x , ( x ∈ { 0,1 } ) ,垂直方向记为 y , ( y ∈ { 0,1 } ) 。
Fig. 2 Quadrant of QHC partition

图2 QHC中的象限划分

定义1：QHC算法中4个象限编码顺序依次为：
I : x = 0 y = 0 ， II : x = 0 y = 1 ， III : x = 1 y = 0 ， IV : x = 1 y = 1 （2）
四分象限对于任意限定层级（k）的格元行列号（GridX, GridY）,自高位（右侧第k位）逐个取得GridX和GridY的相应位值,即可逐层计算格元在各层级中所处的象限。取位的函数（取Value的自右第index位的值）可记为Bit（Value, index）。
以3层填充空间为例,取坐标为（6,5）的格元,即(GridX=6, GridY=5),其格元坐标的二进制表达为： GridX = 0110 , GridY = 0101 。自高位（k=3）逐层取位结果为：[1,1],[1,0],[0,1],可知格元在第1层位于第IV象限,第2层位于第III象限,第3层位于第II象限。图1可验证此计算结果的正确性。
Fig. 1 Hilbert space-filling curve (k=1,2,3)

图1 Hilbert空间填充曲线(k=1,2,3)

3.2 位域共同体

QHC算法中需要频繁调用取位函数,传统算法为实现取位的结果需将整型数据转换为二进制的字符串数据,取得相应位的字符,最后再将字符转换为数值型来实现。对计算机而言,类型转换是强CPU耗用和长耗时的操作。
为提升QHC算法中取位运算的效率,本文定义了位域共同体[15]来减少数值型与字符串型之间的转换,进而减少CPU占用并提升算法效率,所定义的位域共用体及取位的宏如表1所示。
Tab. 1 Definition of Bit-Field structure and macros

表1 位域共用体及宏定义

名称 伪代码
位域共同体宏定义 union unGridPos
{
unsigned short s;
struct
{
unsigned a1:1;
unsigned a2:1;
…
unsigned a16:1;
} a;
};
//定义取位运算的宏
#define GetBit(x,y) switch(x){\
case 1: y=X.a.a##1;break; \
case 2: y=X.a.a##2;break; \
…
case 16: y=X.a.a##16;break; \
}

3.3 演化规则与复制特性

Hilbert曲线具有显著的自相似及自复制特征[16]。若将 k 取不同值时的Hilbert曲线记为 F k ,可用 F k + 1 表示下一层的Hilbert曲线,用 F k + 表示将 F k 顺时针旋转90°,用 F k - 表示将 F k 逆时针旋转90°。Hilbert曲线的自相似与自复制特征体现为如下规律：
（1） F k + 1 按象限等分为4个部分;
（2）其中,第II和第IV象限为 F k 的平移;
（3）第I象限是 F k 顺时针旋转90°的结果;
（4）第III象限是 F k 逆时针旋转90°的结果;
当k≥0时,Hilbert曲线的每一层都符合上述规律,所以,Hilbert曲线的每一层都可由上一层的曲线推导出来（图3）。
Fig. 3 Self-similarity characteristics of Hilbert

图3 Hilbert曲线的自相似特征

Hilbert曲线的自相似与自复制特征决定了Hilbert曲线中只有4种形态的基元曲线,不管多复杂的Hilbert曲线皆由这4种基元曲线递归派生而来[7,17]。这4种基元曲线形态如图3所示,本文将曲线形态分别记为type0、type1、type2、type3（图4）。
Fig. 4 Four state vector types in Hilbert curve

图4 Hilbert曲线中4种状态向量

若type类型确定,便可确定其4个子象限中的key值,4种type对应的4个子象限key值状态向量分别为：{0,1,3,2}、{0,3,1,2}、{2,3,1,0}、{2,1,3,0};进一步分析研究发现,当type类型确定,其4个子象限的type值亦就确定下来,4种type对应的4个子象限type值状态向量分别为：{1,0,3,0}、{0,2,1,1}、{2,1,2,3}、{3,3,0,2}。
基于上述分析,可形成Hilbert曲线状态转移向量矩阵,如式（3）所示。
M = 0 0 0 0 1 0 0 1 1 0 0 1 0 3 3 0 1 1 2 0 1 0 0 0 0 1 0 1 3 2 1 1 0 1 1 1 1 1 2 1 2 0 0 2 2 2 0 1 3 1 2 1 0 1 2 2 1 1 0 3 3 0 0 2 3 3 0 1 1 3 3 1 0 3 0 3 1 1 0 2 （3）

3.4 QHC算法流程及伪码

状态转移向量矩阵可逐层分解计算Hilbert编码,其算法流程描述如下：
（1）设定需计算编码的格元坐标为GridX, GridY,分别为在栅格中的列号和行号,格元所在层级为k,k>0,有 GridX < 2 k , GridY < 2 k 。
（2）最顶层(k=1)的 type 为0,记为：
type(k=1)=0（4）
（3）四分象限中各象限的 key 可以由 ( type , x , y ) 唯一决定,可通过查询状态转移向量矩阵获得,记为：
key = K type , x , y （5）
（4）四分象限中各象限的 type 也可以由(type,x,y)唯一决定,可通过查询状态转移向量矩阵获得,记为：
nextType = T type , x , y （6）
（5）若曲线只有一层(k=1),则通过式（7）计算所得到的结果即为Hilbert编码。
key = K 0 , x , y （7）
（6）若曲线多于一层( k > 1 ),算法伪码如表2所示。
Tab. 2 Pseudo-code of QHC algorithm

表2 QuickHilbertCode算法伪码

算法 伪代码
QHC Input:
GridX, GridY:格元所在的列号和行号
k:填充空间的层级
Initial:
arKey[4][2][2] = {0,1,3,2,0,3,1,2,2,3,1,0,2,1,3,0}：key向
量矩阵
arType[4][2][2] = {1,0,3,0,0,2,1,1,2,1,2,3,3,3,0,2}：type
向量矩阵
Output:
resKey：格元的Hilbert码
nType=0; bitX=0; bitY=0; resKey=0;
for level in k..1 do
bitX=GetBit(GridX, level); bitY=GetBit(GridY, level);
resKey = (resKey<<2)|arKey[nType][bitX][bitY];
nType = arType[nType][bitX][bitY];
end for
return resKey;

3.5 QHC算法分析

从上述算法流程描述及伪码可看出,算法内只有一层for循环,不需要迭代过程,故算法时间复杂度为O(n),n为填充空间的层级。对于给定层级的任意格元,算法具有稳定一致的时间耗费。与之前的研究成果相比,本文提出的算法具有如下优势：
（1）算法无递归操作,具有单层for循环,算法时间复杂度为O(n),且稳定;
（2）使用位域结构体及宏实现取位运算,避免了类型转换,减少了代码量,降低了内存占用,且提高了计算效率;
（3）定义了Hilbert的key及type的转移向量矩阵,减少了代码中的分支及判断,提高了效率;
（4）算法由单个函数实现,独立性强,耦合度低,不依赖于其他函数或库,更易于实现并行计算及GPU计算。

4 算法验证与结果分析

算法实验环境为1台CPU Intel(R) Core(TM) i7-4770 CPU @3.40 GHz,16 GB内存的台式计算机。
（1）算法正确性验证
Hilbert曲线具有分形特征,体现为自相似与自复制特性,第k+1层由第k层平移并旋转推导得出。不同算法间,若前5层结果一致,由算法的递归特性可知,其后续结果也是一致的。本文对QHC算法与二进制位操作算法的结果进行一一比对来验证QHC算法的正确性。
定义2：二进制位操作算法(Faloutsos&Roseman)计算所得Hilbert码为 FR k , x , y 。
定义3：本文算法(QHC)计算所得Hilbert码为 QHC ( k , x , y )
定义4：格元码差值 ΔH ( x , y ) 。对同一格元,两种算法所得结果差值的绝对值：
ΔH k , x , y = abs ( FR k , x , y - QHC ( k , x , y ) ) （8）
定义5： k 级填充空间的格元码差值之和 G k 为：
G k = ∑ x = 0 2 k ∑ y = 0 2 k ΔH k , x , y （9）
本文实现了定义2至定义5的C++代码[18],通过实验验证了1-16级的计算结果,实验结果显示当k值取1-16时,所得到的结果 G k 均为0,证明本算法计算所得到的Hilbert码与二进制位操作算法计算所得的Hilbert码一致,证明16级之内本文算法结果正确。本文算法的C++代码可参阅文献[18]。
（2）算法性能验证
陆锋等提出了基于空间层次分解的Hilbert码生成算法,并将之与二进制位操作的Hilbert码生成算法作了性能对比实验[13],实验结果如表3所示。
Tab. 3 Comparison of performance between spatial hierarchical decomposition and FR (s)

表3 Hilbert空间排列码算法效率比较(s)

格网数 512×512 768×768 1024×1024 2048×2048
位操作 3 7 12 53
层次分解 1 4 7 31
本文采用类似的实验思路,自1-16各层级以2种算法计算填充覆盖完整空间的Hilbert码,计算20次取平均时间,验证算法伪码如表4所示,具体代码可参阅文献[18]。
Tab. 4 Pseudo-code of verification algorithm

表4 验证算法伪码

名称 伪代码
验证算法 Input:
k: 填充空间的层级
nTimes: 重复计算的次数,重复计算后取平均值以提高实验结果精确度
Initial:
dQHCTime = 0.0: QHC算法耗费的时间
dFRTime = 0.0: 二进制位运算算法耗费的时间
Output:
for level in 1..k do
for t in 1..nTimes do
dQHCTime+=QHC(level); 以QHC算法计算第level层所有格元耗费的时间
dFRTime+=FR(level); 以二进制位运算算法计算第level层所有格元耗费的时间
end for
end for
returndQHCTime, dFRTime;
从实验结果可知,与二进制位操作的Hilbert码生成算法相比,本算法的效率倍数平均提升为2036.4倍,而且随层级逐渐增大,填充空间中的格元数（Nk）快速增加,本算法的加速效果稳定提升,如图5所示。实验结果摘录（取第7至第14级）如表5所示。本文算法正确性验证实验及性能对比实验的完整C++代码可参阅文献[18]。
Tab. 5 Result of performance test (s)

表5 性能对比实验结果

层级 格元数 QHC算法 FR算法 性能提升
7 128×128 0.000010 0.016548 1591
8 256×256 0.000038 0.068585 1792
9 512×512 0.000152 0.285086 1842
10 1024×1024 0.000604 1.189585 1971
11 2048×2048 0.002365 4.976408 2104
12 4096×4096 0.009638 21.498139 2231
13 8192×8192 0.038625 89.662764 2321
14 16384×16384 0.151542 369.571435 2439
Fig. 5 Speed-up ratio of QHC

图5 算法性能提升

5 结语

本文基于状态转移向量矩阵设计了一种快速计算Hilbert码的算法,并结合位域共用体技术,简化了算法代码,提高了算法执行效率。本文对1-16级的完整填充空间进行了验证测试,实验结果表明,本文算法的计算结果正确,与二进制位操作的Hilbert码生成算法相比,在各层级均有很高的加速效果,且随着待计算数据量增大算法加速效果稳定提升,与空间层次分解的Hilbert码生成算法相比其性能显著提升。同时本文提供了算法及相关实验的完整C++代码[18]。
本文算法面向二维平面的空间填充曲线编码计算,后续研究工作中,还需要对三维[19]以及多维空间的Hilbert编码[20]实现相应的优化改进和实验。

The authors have declared that no competing interests exist.

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
Christos Faloutsos and Shari Roseman. Fractals for secondary key retrieval[C]. Proceedings of the 8th ACM SIGACT-SIGMOD-SIGART Symposium on Principles of Database Systems, Philadelphia, 1989,1:247-252.

[2]
齐琳,沈婕,郭立帅,等.面向D-TIN并行构建的动态条带数据划分方法与实验分析[J].地球信息科学学报,2012,14(1):55-61.

[3]
许虎,聂云峰,舒坚.基于中间件的瓦片地图服务设计与实现[J].地球信息科学学报,2010,12(4):562-567.

[4]
周艳,朱庆,张叶廷.基于Hilbert曲线层次分解的空间数据划分方法[J].地理与地理信息科学,2007,23(4):13-17.

[5]
Meng L K, Huang C Q, Zhao C Y, et al.An improved Hilbert curve for parallel spatial data partitioning[J]. Geo-Spatial Information Science, 2007,10(4):282-286.

[6]
Castro J, Georgiopoulos M, Demara R, et al.Data-partitioning using the Hilbert space filling curves: Effect on the speed of convergence of Fuzzy ARTMAP for large database problems[J]. Neural Networks, 2005,18(7):967-984.

[7]
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.

[8]
陆锋,周成虎.一种基于Hilbert排列码的GIS空间索引方法[J].计算机辅助设计与图形学学报, 2001,13(5):424-429.

[9]
Abel D J and Mark D M. A comparative analysis of some two-dimensional orderings[J]. International Journal of Geographical Information System, 1990,4(1):21-31.

[10]
王永杰,孟令奎,赵春宇.基于Hilbert空间排列码的海量空间数据划分算法研究[J].武汉大学学报(信息科学版),2007,32(7):650-653.

[11]
Chen H L and Yein Chang. Neighbor-finding based on space-filling curves[J]. Information Systems, 2005,30(3):205-226.

[12]
陈宁涛,王能超,陈莹. Hilbert曲线的快速生成算法设计与实现[J].小型微型计算机系统,2005(10):1754-1757.

[13]
陆锋,周成虎.一种基于空间层次分解的Hilbert码生成算法[J].中国图象图形学报,2001,6(5):465-469.

[14]
王笋,徐小双. Hilbert曲线扫描矩阵的生成算法及其MATLAB程序代码[J].中国图象图形学报,2006,11(1):119-122.

[15]
管晓春,陈孝敬,吴桂初.基于外设位域结构体的DSP的C语言程序设计的教学研究[C]. 2010 Third International Conference on Education Technology and Training (ETT)论文集,中国湖北武汉,2010,6:208-211.

[16]
汪璇,曹万强. Hilbert变换及其基本性质分析[J]. 湖北大学学报(自然科学版),2008(1):53-55.

[17]
赵春宇,孟令奎,林志勇.一种面向并行空间数据库的数据划分算法研究[J].武汉大学学报(信息科学版),2006,31(11): 962-965.

[18]
Li S J.An Algorithm for Hilbert Ordering Code Based on State-Transition Matrix[EB/OL].2013.

[19]
Zhang J and Kamata S I. A generalized 3-D Hilbert scan using look-up tables[J]. Journal Of Visual Communication And Image Representation, 2012,23(3):418-425.

[20]
李晨阳,段雄文,冯玉才. N维Hilbert曲线生成算法[J].中国图象图形学报,2006,11(8):1068-1075.

Options
文章导航

/

〈 〉