Hierarchical Hexagonal Clustering and Indexing

: Space-ﬁlling curves (SFCs) represent an efﬁcient and straightforward method for sparse-space indexing to transform an n -dimensional space into a one-dimensional representation. This is often applied for multidimensional point indexing which brings a better perspective for data analysis, visualization and queries. SFCs are involved in many areas such as big data analysis and visualization, image decomposition, computer graphics and geographic information systems (GISs). The indexing methods subdivide the space into logic clusters of close points and they differ in various parameters including the cluster order, the distance metrics, and the pattern shape. Beside the simple and highly preferred triangular and square uniform grids, the hexagonal uniform grids have gained high interest especially in areas such as GISs, image processing and data visualization for the uniform distance between cells and high effectiveness of circle coverage. While the linearization of hexagons is an obvious approach for memory representation, it seems there is no hexagonal SFC indexing method generally used in practice. The main limitation of hexagons lies in lacking inﬁnite decomposition into sub-hexagons and similarity of tiles on different levels of hierarchy. Our research aims at deﬁning a fast and robust hexagonal SFC method. The Gosper fractal is utilized to preserve the beneﬁts of hexagonal grids and to efﬁciently and hierarchically linearize points in a hexagonal grid while solving the non-convex shape and recursive transformation issues of the fractal. A comparison to other SFCs and grids is conducted to verify the robustness and effectiveness of our hexagonal method.


Introduction
A space-filling curve (SFC) represents a bijective transformation of multidimensional points onto 1D representation by computing an index (hash code) for each point defining the linear order of points in the space [1].A SFC hashing algorithm usually follows a pattern of some tree or pyramid structure.SFCs are utilized, e.g., for space indexing and range querying [1][2][3] and also for building spatial hierarchies especially for parallel platforms [4,5].While the SFCs based on orthogonal and triangular uniform grids are well investigated and preferred in the applications, we propose an effective and robust hexagonal Node-Gosper SFC based on the Gosper fractal [6] (Figure 1).This paper analyzes the problems of hierarchical hexagonal indexing and proposes solutions of issues related to transformations of multi-resolution hexagonal grids and fractal shape of Gosper tiles.Our novel method overcomes the lacks of naive Gosper-based approaches [7,8] and competes well with other non-hexagonal SFCs.A point cloud (PC) is a collection of points of the same type (2-3D real points) [9] representing a discrete form of real objects.Beside traditional PCs produced mostly by 3D scanners [9], a set of points can also represent general spatial information [10] such as location, GPS coordinates on the map [11], swarm particles [12], nodes of graph generated by visualization algorithms [13], etc.We focus on 2D point clouds as they are highly investigated in the areas such as data mining [2,3,10], computer graphics [14,15], image analysis [16][17][18], geographic information systems (GISs) [19][20][21][22][23], sensor networks [24,25], triangulation [26][27][28][29] and data visualization [12].In the mentioned application fields, the point indexing/clustering algorithms are crucial, e.g., for efficient organization of data in the memory, analyzing their properties, visualization of clusters, computation of forces in particle systems and searching near neighbors.Point indexing algorithms decompose the space into subspaces or tiles [1,10] containing a subset and create structures for effective addressing of the points.Usually, the structures utilize the uniform regular grids that are represented by list of clusters [30], or the clusters are organized in the hierarchical manner to reduce the costs of point queries [2,10] by omitting irrelevant clusters such as empty or distant ones (Figure 1a).The hierarchies offer good query times, esp.for sparse datasets, but they are often memory consuming and unbalanced [5,31].Linear methods save memory usage and provide constant query time for uniformly distributed datasets [30].However, they are inefficient for sparse point clouds.
The space-filling curves combine both approaches as they define the linear order of clusters and preserve the hierarchical information about points in the hash codes.To compute an index uniquely addressing a point in a hierarchy covering the whole space, the multi-resolution grids are often used.Most widely used SFCs including Peano, Z-order, Hilbert and Sierpi ński SFC [1,[32][33][34] are for orthogonal or triangular grids.The multi-resolution hexagonal grids were utilized for image processing [18,[35][36][37].Burt [38] proposed a structure for coding binary images with hexagonal grids and Gibson and Lucas [39] extended this structure to the more sophisticated Generalized Balanced Ternary (GBT).Those methods serve to decompose and represent images with pixels having positive coordinates, but they do not offer any general method for real points indexing and querying.Several papers discussed the problems of different grid systems and their mutual conversion to combine their best properties [21][22][23], but they all miss the correct hierarchical 1-to-7 refinement in the hexagonal grid (Figure 2) for a correct hexagonal SFC.During the last years, the hexagons have gained increased interest in scientific areas including geographic information systems [19][20][21][22][23]40,41], data clustering and querying [36,38,39,42,43], data visualization [44,45], computer graphics [46] and computer or sensor networks [24,25].Recently, the transportation network company Uber Technologies Inc. introduced a discrete global grid system [11] based on multi-resolution hexagonal grids.They use it to localize cars and customers on the map, to efficiently optimize ride pricing and dispatch and to visualize and explore spatial data.The GPS coordinates of locations on the map represent the points that are indexed by hexagonal grids.The hexagons also form the basis for analysis of the Uber marketplace, they can easily approximate zones defined by edges of required area on the map and they minimize the quantization error [11].Unlike squares and triangles, hexagons have uniform distance between neighboring cells and they have the optimal perimeter/area ratio which leads to good approximation of circles [47].
While the hexagonal grids have a wide range of applications and a convenient linear representation could simplify performing point indexing, the hexagonal SFCs are not utilized very often.The main obstacle is that a hexagon is not a rep-tile.A rep-tile is a tile which is infinitely decomposable into similar tiles.A hexagon can be decomposed only into triangles, thus the hierarchy has to be defined by bottom-up composition of hexagons creating a system of multi-resolution grids [10].After exploring different hexagonal patterns, we decided to follow the concept of the Gosper fractal (or Flowsnake) which is a self-similar recursive fractal discovered by William Gosper in 1973 and introduced by Martin Gardner in 1976 [6,48] (Figure 2).The fractal connects the vertices of hexagonal lattice and creates a pattern of 7 hexagons called the Gosper island.These islands can be composed again into 7-pattern to create an island of greater level (Figure 2).We define a procedure extending the Gosper fractal to its node form (Node-Gosper curve) for point indexing in the manner of SFCs.The key tasks are: point localization in the hierarchy and inverse point mapping algorithm.The first published method [7] simply decomposes the points into 7 disjoint areas in the top-down manner to approximate the hexagonal 7-pattern.The paper [7] shows that the top-down decomposition is not reliable and it cannot produce a correct hexagonal hierarchy.The Gosper islands are non-convex fractal tiles (Figure 2) with jagged edges that are fitted into surrounding ones.This leads to point mislocalization in the hierarchy (esp.at the edges) which produces many non-hexagonal clusters in the lower levels of decomposition.Thus, the hexagonal properties are not preserved.Later, the inverse bottom-up mapping procedure fixing those problems was theoretically proposed [8], but its mapping algorithm is inaccurate and inefficient.To the best of our knowledge, there is no straightforward hashing algorithm mapping the points onto the Gosper curve.This paper summarizes a background about hexagonal grids and Gosper fractal to understand the advantages and lacks of hexagons and hexagonal multi-resolution grids.The main contribution of this work lies in the derivation of the hexagonal space-filling curve solving the lacks by using the pattern of the Gosper fractal to make a hexagonal linearization available to wide range of applications.We focus on correct point localization in the Gosper islands, stable point mapping onto the Gosper curve, fast determination of related hexagonal pattern, its rotation and indexation and efficient design of these algorithms.Our Gosper SFC method is supported by several theoretical metrics, e.g., Honeycomb conjecture [47] proving good coverage of a circle by hexagonal grid, distance metrics [49] declaring good real distance between points neighboring along the curve and construction times comparison.The metrics are verified by practical experiments comparing different grids and SFCs which show that our Node-Gosper SFC is robust and it produces the correct hexagonal hierarchy without mislocalization much faster than previous approaches [7,8].
The paper is organized as follows.Section 2 describes the related theoretical background we need for hierarchical hexagonal clustering.Section 3 gives the details of our hierarchical hexagonal composition and the corresponding Node-Gosper SFC construction.Finally, Section 4 evaluates our SFC by comparing against other grids and space-filling curves.

Background of Hierarchical Hexagonal Clustering
This section summarizes the theoretical and mathematical background needed to define and understand our Node-Gosper SFC.At first, we state why we prefer hexagonal grids over orthogonal and triangular ones (Section 2.1) and briefly survey selected properties of hexagonal grids.Section 2.2 overviews the existing SFCs and defines the criteria that we generally expect from a good SFC.Afterwards, we include a brief analysis of hexagonal multi-resolution grids and explanation why the Gosper fractal is selected for point indexing (Section 2.3).

Hexagonal Grids
In mathematics, a tiling means decomposition of the space into collection of geometric shapes representing some logical subspaces covering the whole space.There are only three types of regular tilings [10]: Triangular, square and hexagonal (Figure 3).Unlike other tilings, hexagons have uniform distances between tiles, a hexagon is the regular shape closest to the circle and one lattice vertex is shared by only three hexagons.The outstanding circular coverage by hexagonal grids is supported by several theoretical metrics including isoperimetric inequality [50], Honeycomb conjecture [47] and circle packing problem solution [51] which show that the best regular shape having the maximal area with the least perimeter is the hexagon.The chart in Figure 4 confirms experimentally the predominance of the hexagonal grid over the orthogonal grid.These properties are useful for disc queries and nearest neighbors search as the number of tested tiles should be minimized [42].On the other hand, hexagons are not infinitely decomposable.We solve this by following the principles of Gosper fractal that gives hexagons a hierarchical organization.As the points in a PC are real points, the discrete hexagonal coordinate system has to be defined to localize them in the hexagonal grid and give them integer coordinates of the corresponding tile.The axial (q, l) ∈ Z 2 and cube (x, y, z) ∈ Z 3 coordinates are used to address the hexagons.The cube coordinates extend the axial ones by one redundant coordinate to get the six-directional system simplifying some calculations.If the condition x + y + z = 0 is met, the axial coordinates are calculated simply by neglecting the redundant y coordinate q = x, l = z.We refer to papers [43,52] for more details and visualizations.

Space-Filling Curves
Space-filling curves represent a family of point hashing algorithms which define the basic standards of space linearization [1].These standards are carefully taken into account when designing our Node-Gosper SFC (see, e.g., Figure 1b).Space-filling curves transform a general n-dimensional point into a one-dimensional (1D) index by localizing it in multi-resolution grids (Figure 5).For each point, a set of hierarchical indices addressing the sub-clusters containing the hashed point is computed.The hash code is usually a bit string merged from the indices in the top-down manner so that the root index occupies the most significant bits [1].After sorting the points by codes an SFC is constructed and the points from the same cluster (with the same hash) are aligned in the memory.An SFC order of clusters represents the depth first passage of the corresponding tree (Figure 5).The nearby points are stored close to each other in the memory.However, there can be some points lying very close on the curve, but in reality, they are far away (see, e.g., Figure 5a).The shape of clusters and their indexation are varied, and thus, there are many different curves.The popular ones are Z-order, Hilbert or Sierpinski space-filling curves [1,[32][33][34].SFCs are very straightforward and efficient methods for sparse-space clustering.Such a structure design is also much more convenient for massive parallel processing of clusters [4,5,12,13,53].
The linear point memory representation is used, e.g., for space indexing and range querying [1-3] or sophisticated triangulation algorithms [26,27].SFCs are often utilized for building spatial hierarchies especially for parallel platforms [4,5].The linearly stored points can be grouped together in parallel to build a space tree (e.q.Octree) instead of a sequential particle insertion [53].Papers [12,13] discuss the utilization of SFCs for speeding up the nearest neighbors (NNs) algorithm and real time graph layout visualization.To bring a solid hexagonal method, we define the general criteria for good SFCs [1] in Section 2.2.

Space-Filling Curves Criteria
The SFCs differ in their properties, so that there is no SFC that works efficiently for all kinds of applications.However, several criteria should be met to define a high-quality SFC [1].The conditions 1-3 are mandatory followed by quality criteria: 1.
Mapping between a point and SFC code is bijective-each point is assigned a unique hash code defining the cluster along the SFC.This has to be possible for greater depths of hierarchy.

2.
The curve is space-filling-every point of space lies on the curve.

3.
Total ordering-each cluster is visited by an SFC only once.4.
Mapping between points and codes should be easy to compute. 5.
Two points close along the curve should be also close in the multidimensional space.6.
Sequentialisation of a contiguous data set should also lead to a contiguous sequence of code numbers.This is handy for efficient memory addressing.

Hexagonal Curve Selection
To best fulfill the requirements for good SFCs, many hexagonal and combined fractal tiling patterns have been studied (see, e.g., [17,37,38,42]).The fractal patterns usually have very complex shapes, and thus, they are not appropriate for point clustering.There are three reasonable variants of multi-resolution hexagonal grids often used in the Discrete Global Grid Systems [21][22][23]41,54].Generally, the grids of aperture 3 and 4 are preferred (see Figure 6) to avoid the problems with fractal shapes.While these grids can be used for localization, they do not preserve the basic properties of the regular hexagonal grids.A grid of level i ∈ N\{0} divides the hexagons of level i + 1 into non-hexagonal polygons, which have to be joined together on the level i + 1.These hierarchies are incongruent, and therefore, the classical tree-based algorithms cannot be utilized with them.Thus, only the hexagonal grid of aperture 7 (Figure 6C) is applicable for an SFC construction (Figure 2) because the composites are completely disjoint (not overlapping).As a result, the points of each island can be aligned and linearly stored in the memory in the manner of an SFC.As mentioned before, the Gosper curve is a curve based on such a hierarchy.The space filled by the i-th recursion of this curve is bounded by the level-i Gosper island (Figure 2).The Gosper curve and its islands have many notable properties deriving from the essence of the hexagonal grid and its distance characteristics (see Figure 3) [10,37,42].The Gosper curve satisfies the conditions of the good SFCs such as the continuity, uniformity, scale invariance, etc. [49,55] (see Section 2.2).Each point in the space is visited just once, and the distance between adjacent clusters along the curve is uniform.Any hexagon (or Gosper island) shares edges with six other identical shapes so that the Gosper islands can be simply addressed by the same coordinate system as the regular hexagons.All the points of the dataset lie on the curve.The whole point cloud just has to be fitted into the circle inscribed into the Gosper island.The mapping is bijective for sufficiently deep hierarchies.The specific indexation of clusters and mapping complexity depend on the final algorithm.The Gosper curve is also held up by several theoretical metrics.The worst-case locality value (W L) [49] (smaller is better) examines the rate between the Euclidean distance of two points in the real 2D space and the distance along an SFC.The Gosper curve (W L = 6.35) is comparable with sophisticated non-hexagonal SFCs, e.g., Hilbert (W L = 6) or Sierpi ński (W L = 4) curve.The Arrwwid number [42] (smaller is better) represents the maximum tiles hit by reasonably small radius and judges the suitability of a tiling (or an SFC) for a circular neighborhood querying.The Gosper curve has optimal Arrwwid = 3 while other named SFCs have at least Arrwwid = 4.

Node-Gosper Space-Filling Curve
We mentioned the main issues related to hexagonal hierarchy such as impossibility of recursive decomposition and fractal shape of Gosper islands that make it difficult to design a correct point mapping algorithm.The solutions proposed by previous work are insufficient.The approximative top-down decomposition [7] following the fractal drawing algorithm produces deformed polygons and obviously incorrect point localization.The naive bottom-up approach [8] solving the issue with point localization and producing the correct Node-Gosper SFC proposes only a theoretical pipeline with bad accuracy and performance.
This section represents the main contribution of the paper addressing all the named issues related to hierarchical hexagonal point indexing.We improve the existing simple concepts from papers [7,8] and define the complete hexagon-based system in detail with mathematical definitions and pseudocodes.
At first, the mathematical definition and the recursive construction of the Gosper fractal is reminded in Section 3.1.Following the transformations defined for the Gosper fractal, we introduce our Node-Gosper indexation model in Section 3.2 which gives the right order to the hashed points.Next, we define the procedure computing a hexagon centroid according to its hash code based on the indexation model and top-down drawing algorithm of the Gosper fractal.Section 3.4 defines our novel mapping algorithm representing the hashing function which correctly maps a point onto the Node-Gosper curve and returns the corresponding hash code.

Gosper Fractal
Our Node-Gosper SFC is inspired by the Gosper fractal drawing algorithm which is reminded here.The Gosper fractal is based on the pattern of seven hexagons (Figure 7a).The fractal is recursively constructed by replacing the oriented red arrows with the properly rotated and scaled pattern in the top-down manner (Figure 7b).The arrows connect the vertices of the hexagonal grid, and thus, this iterative method secures the geometrical continuity of the patterns and draws the correct image of the fractal.Alternatively, the same procedure can be applied to the node-based pattern (Figure 7c) connecting the hexagon centers which was probably first mentioned by Ventrella [56] and its design is better for point indexing than the vertex one.The derivation of exact mathematical properties of the Gosper fractal is crucial for the Node-Gosper SFC definition and we briefly remind them here.See papers [7,8,57,58] for details. Figure 8 describes the transformation of the pattern replacing the base hexagon inscribed in the unit circle of size = 1 (black dashes) by seven smaller hexagons (level-1 Gosper island).A hexagon itself practically represents a level-0 Gosper island.The figure also shows two features of the Gosper islands: The Gosper island preserves the area of the parent hexagon, but its perimeter limit is infinity.The scale factor r = 1 √ 7 can be computed according to the blue triangle (Figure 8) using the Law of Cosine and it also represents the size of child hexagons.The angle α = arcsin can be computed using the Law of Sine and the scale factor r. Figure 7a shows that the patterns replacing differently oriented arrows have to be additionally rotated.The angles of orientation are defined by function f (y) in (1).Going through the basic Gosper pattern (Figure 7a) from the first vertex (0, 0) to the last vertex (1, 0), the hexagons 0 and 3 are rotated by −120 • and the hexagon 5 by 120 • .The transformations can be easily computed using the well-known matrices for rotation (2), scale and translation.

Node-Gosper Indexation
While the recursive method described in the previous section serves only to draw an image of the Gosper fractal with geometrical continuity, the indexation model must precisely state the order of the clusters along the SFC.Each point p of the point cloud is localized in the hexagonal grid and the corresponding hexagon is assigned a bit hash code representing its serial number.This section explains our indexation model of the hexagons ordered by the Node-Gosper SFC.
The basic patterns defining rotation and indexation are shown in Figure 9.The main point is that the subordinate patterns have different order of indices depending on the indexation of the parent pattern.For the indexation model, the geometrical continuity is insufficient.Thus, two structural patterns are described in Figure 9b,c.They define both forward (F) and backward (B) indexing variants.The F-pattern is initial and the curve indexation is calculated in the top-down manner.Each hexagon is recursively replaced by a corresponding subordinate pattern depending on the designation of its parent.Notice that the F/B designation in the backward indexation is just the opposite direction to the forward indexation.This secures the smooth passage of all hexagons along the curve without any large leaps and crossings.The final hierarchical indexation based on the patterns from Figure 9 is displayed in Figure 10.A hash code consists of n ∈ N\{0} indices k = (k 1 , k 2 , . . . ,k n ), where k i ∈ N, addressing the Gosper islands containing a query point on each level i ∈ {1, 2, . . ., n} of the hexagonal hierarchy.The indices of the upper hierarchical levels are represented by more significant trinities and the lower ones by less significant trinities.The hash code is computed Inversely, an index k i can be obtained from code by simple bit operations.An index k i is a number from 0 to 6, but the three bits can represent maximally the number 7. Thus, there is unfortunately a numerical gap breaking the requirement of a contiguous sequence of indices (codes) stated in Section 2.2.However, it does not have to be a crucial problem.While having a 64-bit unsigned integer, a 21-level hexagonal indexation can be defined.

Localization with Node-Gosper Curves
As emphasized before, the hexagons are not hierarchically decomposable.However, the Gosper fractal (Section 3.1) and Node-Gosper indexation (Section 3.2) are defined recursively in the top-down manner.At first, we define a top-down algorithm computing a centroid of hexagon hierarchically addressed by indices k.In the next section, we use this method to define the inverse algorithm mapping the points onto the Node-Gosper SFC.
Let us consider the case of the base hexagon with size = 1 (Figure 11), which is recursively replaced by the Node-Gosper patterns (Figure 9) up to the level n.Given a vector of indices k = (k 1 , k 2 , . . . ,k n−1 , k n ) addressing a hexagon in a level-n Gosper island, a center point of such a hexagon is computed.This algorithm represents the backward mapping which verifies the bijectivity of the Gosper mapping.The subordinate patterns can be indexed in the forward (F) or backward (B) manner (Figure 9), which affects the addressing order of the hexagons by k.A pattern passing order F/B of the i-th recursion can be determined depending on k as follows: Each k i index has to be recalculated to index k i ∈ N in the F order (Figure 9) to unify addressing of the corresponding Gosper islands in the top-down manner.The calculation of k = (k 1 , k 2 , . . ., k n−1 , k n ) is defined as follows: For example, a point p hashed onto indices k p = (5, 1) according to the level-2 Node-Gosper curve (Figure 10) belongs to the hexagon 5 (Order(2) = B) on the first hexagonal level (Order(1) = F), whose subordinate island thus has reverse indexation of hexagons.Therefore, the corresponding index is 1 on the second level of recursion and not 5, as it would be according to the F order, which means k p = (5, 5).A pattern rotation is considered in the following transformations.
k is used to compute the center point of the addressed hexagon.The algorithm recursively transforms two vectors: Position vector c i and directional vector d i (see Figure 11), where i is the current level of recursion.A vector c i represents the initial center position of the current Gosper island.The vector d i represents the rotation and scale of the grid on the i-th level of recursion.The initial directional vector is pointing to the neighboring hexagon and it is defined as d 0 = (w, 0), where w is the width of the base hexagon.The d 0 defines the basic orientation and distance between two neighboring hexagon centers before the first transformation.Depending on k , d i is defined by the recurrent formula: where 1 √ 7 • R α represents the scale by r = 1 (Section 3.1), and the R f (k i−1 ) represents the additional rotation defined by (1).The grid transformation can be simply defined by the following matrix: Next, a vector d k i i pointing to the specific neighboring hexagon addressed by k i is defined as: where γ(k i ) represents the corresponding angle by which d i has to be rotated, so that the d k i i points to the hexagon addressed by the k i index.γ(y) is defined as: γ(y) = γ y , where y ∈ {0, 1, . . ., 6} and γ = (0 where γ(y) equals to a multiple of 60 • .The hexagon with y = 0 is initial, and therefore, d 0 i = d i .For example, d i is rotated by 120 • to obtain the second hexagon center d 2 i (see Figure 11).The center point c i is computed by a translation of c i−1 by the properly rotated directional vector The center point c n represents the final hexagon center addressed by indices k.

Node-Gosper Point Mapping
In the previous section, we described the algorithm computing a hexagon center according to k.In this section, we define the inverse mapping of a point onto the Node-Gosper curve, so that k is computed and the corresponding hash is generated.
The mapping onto the Node-Gosper curve can be done only in the bottom-up manner, because hexagons cannot be decomposed into sub-hexagons (Figure 12).As it was explained in the previous sections, the hexagonal 7-pattern (Gosper island) is rotated and scaled to replace a parent hexagon which defines the required hierarchy of hexagons.The constructed Gosper islands are non-convex fractal tiles, so that they have to be fitted into neighboring islands to tile the plane.The top-down point localization is incorrect, because the jagged edges of islands overstep the borders of parent hexagons, so that a point previously localized in the hexagon (especially close to the edges) may belong to different branch of hierarchy in the next iteration (Figure 12).An important fact about the hexagonal hierarchy is that seven hexagon centers of each Node-Gosper pattern are always located inside their parent hexagon.It means that the upward hexagon centers localization defines numerically stable inverse procedure to the top-down algorithm and it secures the correct hierarchical assignment regardless of the jagged edges.There are still several problems that have to be solved: 1.
The relative position of hexagons of two subsequent grids at level i and i − 1 is not obvious, because information of the complete parental hierarchy is unknown at level i (Figure 12). 2.
The grid rotation varies in the different branches of hierarchy.The specific grid transformation has recursive character and is unknown at level i (Section 3.2).

3.
For the same reason, the indexation and passage order are unknown at level i too.To achieve both the correctness of hierarchical localization and recursive principles of the Gosper fractal, we separate the hashing algorithm into two phases.At first, we localize a point p in the multi-resolution hexagonal grid in the bottom-up manner using the inverse formulas for the Gosper fractal construction (Figure 12).The hierarchical localization gives us just a vector of indices defining the relative location of a point p in the multi-resolution grid.This step secures the correct localization inside the fractal shape of Gosper islands and we describe the procedure in detail in Section 3.4.1.As the Gosper fractal is defined recursively, the specific pattern transformation and indexation for the different recursive branches can be computed only in the top-down manner.Thus, the indices returned by the first part of the algorithm have to be recursively recalculated in the top-down manner (see Section 3.4.2).The second part of the algorithm returns the correct indices that are used to compute the final hash code defining the point order along the Node-Gosper SFC.

Gosper Hierarchical Localization
The task is to determine the relative position of the current hexagon in the context of the parent hexagon.We ignore the difficulties with rotation and passage order in different branches of recursion for start.First, a point p is localized in the level-n hexagonal grid (details in [8,43,52]) using the matrix (12) which returns the axial coordinates a n ∈ Z 2 of the hexagon containing the point p.The next procedure computes only the relative transformation between axial coordinates which is independent on the specific size of hexagons.As mentioned before, the center points of hexagons can be reliably localized inside the parent hexagon.Thus, our transformation between levels i and i − 1 consists of three steps: 1.
Calculation of hexagon center c i according to axial coordinates a i .

2.
Transformation of the c i according to the orientation of the parent grid.

3.
Localization of the transformed center point in the parent grid and calculation of the new axial coordinates a i−1 .
The superior grid of level-(i − 1) should be scaled by √ 7 and rotated by −α.However, it is easier to preserve the same transformation matrix for calculation of all indices and alternatively transform the center point (scaled by 1 , rotated by +α), so the matrix G ( 6) is used instead.The transformation between levels i and i − 1 can be written as follows: where C is the matrix (11) computing the center of a hexagon according to the axial coordinates, C −1 is the inverse matrix ( 12) computing the axial coordinates from the center point and G is the forward transformation matrix (6) (read [8,43,52] for details).
The final constant matrix T is calculated as: a i and a i−1 are the integer axial coordinates of hierarchically corresponding hexagons at level i and i − 1.As the transformation matrix T produces real unrounded coordinates, we denote the parent real axial coordinates as a R i−1 ∈ R 2 that have to be rounded to the closest integers a i−1 ∈ Z 2 .Their hierarchical relationship with child axial coordinates a i ∈ Z 2 can be derived as: The Equation ( 14) computes the transformed axial coordinates, but we also need to compute the hexagon index corresponding to the Node-Gosper pattern described in Section 3.2.Let the parent axial coordinates denote as a R i−1 = (q R , l R ).They are recalculated to the cube coordinates as The decimal part of coordinates represents the relative position between the child and parent hexagon.The relative decimal coordinates As each hexagon center is close to one vertex and edge of the parent grid (see Figure 12), the dominant axis dominant = max{abs(x d ), abs(y d ), abs(z d )} and its sign determine the index of the child hexagon in the view of the parent hexagon.The advantage is that this numeric transfer takes into account the negative coordinates of points, so that the dataset can be positioned anywhere in the space.The relative decimal coordinates v d = (x d , y d , z d ) can be recalculated to the index b i = B( v d ) of the default pattern (Figure 13 P1), where B( v d ) is defined as: The threshold constant Thd → 0 represents some small decimal number (Thd = 10 −5 works well).
As the coordinate transformation works on the level of a unitary grid, the method is numerically stable.The optimized algorithm is summarized in Algorithm 1.
(0,0) (1,0) returns the indices k = (k 0 , . . ., k n ) of a point p {coordinate transformation from cube v to axial a} 1: a.q = v.x, a.l = v.z {for all levels} 2: for i ← 1 to n do {transformation of axial coordinates a} k n−i = B(v d ) 10: end for 3.4.2.Index Conversion Section 3.4.1 described a general method for hierarchical localization of a point p in the level-n Gosper island.This section shows how to recalculate the default indexation model to other models.Especially, it defines a method solving the difficulties like pattern rotation and passage order to construct the correct Node-Gosper SFC.
Let b = (b 1 , . . ., b n ) be a vector of indices returned by Algorithm 1, where n is the number of hierarchical levels.Each b i ∈ 0, 6 can be further transfered to final index k i = E(t, b i ), where E(t, b i ) is the function defined in Table 1.The different indexation patterns are shown in Figure 13 and their conversion is very simple.The final vector of indices k = (k 1 , . . . ,k n ) defines the hierarchical address of a hexagon in the Gosper island to define the final order of hexagons.Each k i can be represented by three bits, so all the indices are masked to the corresponding bits of hash according to the index significance (read Section 3.2).
The procedure computing the correct Node-Gosper curve is more complicated.There are two tasks: 1.
The indices have to be hierarchically recalculated according to the rotation specified by patterns at different levels of recursion (see Figure 9a).

2.
The order of hexagons has to be hierarchically changed according to the forward and backward indexation (see Figure 9b,c).
These properties have recursive nature and they cannot be solved by bottom-up localization.The b has to be recalculated to obtain the indexation model described in Section 3.2.The whole method is described in Algorithm 2. Here, we summarize it in two steps: 1.
As the patterns are rotated only by ±120 • (1), the rotation can be done by index shifting.
The default indexation has the index 0 on the central hexagon and the indices 1-6 form a ring of hexagons.Let the rotDir i ∈ {−1, 0, 1} be a number symbolizing the pattern rotation on the i-th level by {−120 • , 0 • , +120 • }.The shifted index is computed as explained in Algorithm 2 (lines 2-11).

2.
The passage order is reversed depending on the current pattern as explained in Section 3.3.The algorithm is described in Algorithm 2 (lines 12-16).if b i = 0 and rotDir = 0 then 3: end if {transfer to Node-Gosper basic pattern} 8: end if 17: end for

Experiments and Discussion
This section evaluates our novel Node-Gosper SFC based on the hexagonal grid defined in Section 3. Subsection 4.1 compares the complexity and distance metrics of our algorithm with older methods and SFCs.Subsection 4.2 evaluates the pros and cons of the presented curve.

Comparison
Our method is compared with the older approximative top-down Gosper SFC [7], the naive bottom-up approach [8] and several widely used triangular and orthogonal SFCs.After extensive research, we have not found any other Gosper-based SFC method for comparison.The section is divided into three subsections comparing correctness of point localization, execution time and distance metrics.
All experiments run on the following hardware: Intel Core i7-7700HQ @ 2.8 GHz, 16 GB RAM, Windows 10 64-bit 4.1.1.Jagged Edges Incorrect hierarchical point localization during the construction of the Node-Gosper SFC induces wrong point mapping along the curve [7].The general Gosper islands are non-convex fractal tiles, so that the jagged edges of islands overstep the borders of parent hexagons.The previous algorithms [7,8] could not efficiently deal with jagged edges (Figure 14).The experiment presented in Figure 14 uses a point cloud representing the nodes of the level-5 Node-Gosper curve that are hashed according to the top-down and bottom-up methods.The figure confirms that the correct Node-Gosper SFC cannot be constructed by top-down decomposition of hexagons [7], because it produces different types of polygonal clusters that do not preserve the features of the hexagonal grids.The misclassified points are accumulated close to the solid hexagonal edges between the Gosper islands (Figure 14A).Our method uses the bottom-up approach (Figure 14B) composing the Gosper islands from subordinate hexagons, which secures right point localization in the hexagonal hierarchy and produces the correct Node-Gosper SFC passing through the nodes of a regular hexagonal grid.The Gosper islands are undamaged and respecting the jagged edges (Figure 14B).

Efficiency
The efficiency of the curve construction is also a crucial problem, because an SFC with expensive hashing function is not suitable for practical usage regardless of the outstanding theoretical metrics.Next experiment is conducted to compare the execution times of our new optimized Node-Gosper method with other algorithms including the naive implementation of Node-Gosper curve [8] and the Node-Gosper simple curve, which uses our new mapping algorithm (Algorithm 1) without additional index conversion (second recursive part) required for the correct Node-Gosper SFC.The simple variant is faster and can be used for general hierarchical hexagonal model without quality requirements on the final SFC.As the complexity of a hashing function is considered to be constant, the construction times are measured using a randomly generated dataset of 10 6 uniformly distributed points.The chart in Figure 15 shows that our new algorithm is four times faster than the naive old one.In rough comparison with non-hexagonal SFCs (Z-order, Hilbert, Sierpi ńsky), the Node-Gosper hashing algorithm is approximately two times slower.However, this is a logical consequence of applied tiling.Unlike with quads (4-pattern) or triangles (2-pattern), seven sub-clusters have to be considered at each level, thus the number of clusters generated by Node-Gosper curve is much greater for the same hierarchical level than the number of quads, e.g., for n = 8 a complete SFC contains 2 8  The graph shows a comparison of construction times of our Node-Gosper (optimized) method with selected space-filling curves including the previous naive Node-Gosper algorithm [8], top-down Node-Gosper algorithm [7] and simple Node-Gosper (Figure 13 P2) algorithm grouping the hexagons into Gosper islands without additional pattern rotation and clusters reordering.Each measurement was performed on a random dataset with 10 6 points for different levels of clustering.
Definition 1.Let time be the total construction time of the SFC, let level be the hierarchical level until which the SFC is generated and let clusters be the number of sub-clusters of the tiling pattern, the time per cluster is computed as: The metric in Definition 1 takes into account the number of clusters that have to be checked during the evaluation of hashing function which is applicable for any SFC based on the different types of grids.This metric is computed and compared for all SFCs in the chart in Figure 16.The chart shows that our optimized Node-Gosper SFC is actually faster than other good SFCs such as Hilbert and Sierpi ńsky curve and way faster than older variants of the Node-Gosper SFC.Moreover, the tiling with Gosper islands has Arrwwid number α = 3 [42], which means that a sufficiently small disc hits at most 3 tiles at each level, while α = 4 for orthogonal and α = 6 for triangular grids.This feature can be successfully utilized to define an efficient query structure based on the Node-Gosper curve in the future.

Distance Metrics
The last test is focused on the distance metrics of SFCs.The basic requirement for a good SFC is that two points along a curve should be close in the multidimensional space as well.The principal metric here is the sum of distances between the query point and 64 closest points along the curve (32 antecedent, 32 following).The metric average of all the points in the dataset is displayed in the charts in Figure 17.The charts compare four variants of the Node-Gosper curve (patterns Figure 13 P1-P4) and other 3 commonly used SFCs tested with six different datasets summarized in Table 2.The first three ones are randomly generated datasets with uniform or Gaussian distribution.The Gaussian islands is a dataset containing several clusters with Gaussian distribution.The datasets S1, S2 and S3 were downloaded from the web page [59].Experiments show that the Node-Gosper curve is comparable with Hilbert or Sierpi ński curves, its metrics are stable according to the standard deviation and its much better than other simple variants of SFCs.The top-down variant seems to have comparable metrics, but points are not localized correctly along the curve, so that it is not a good option.

Properties
Finally, we summarize the main properties of our proposed curve which fulfill the requirements for a good SFC.

1.
Bijective mapping-Sections 3.2-3.4defined how to map a point onto the Node-Gosper curve and how to compute a cluster point according to the given code.

2.
The curve is space-filling-every point of space lies on the curve.The dataset has to be fitted into the circle inscribed into the Gosper island.

3.
Total ordering-there are no pattern crossings, thus each cluster is visited by SFC only once.4.
Simplicity and efficiency-the designed algorithms are easy and efficient to compute according to Figures 15 and 16.

5.
Localization-points close on the curve are also spatially close (see Figure 17).6.
Sequentialisation-the Node-Gosper SFC does not lead to continuous sequence of code numbers.
The only trouble here is the sequentialisation, because our SFC is based on the 7-pattern, so there is a code gap between the cluster levels.However, it does not have to be a real problem for general clustering purposes.Sparse datasets lead to an incomplete SFC anyway, because the empty areas make an indexation gap along the curve, so most of the clustering methods do not rely on it.Overall, the Node-Gosper curve is a very good hexagonal SFC.

Conclusions
The main contribution of this paper is the new hexagonal space-filling curve (SFC) algorithm called the Node-Gosper SFC based on the recursive Gosper fractal.To the best of our knowledge, there is no paper defining a Gosper-based point indexing method except our previous naive proposals.As the multi-resolution hexagonal grids are highly required in various scientific fields such as computer graphics, image processing, data visualization and geographic information systems (GISs), our hexagonal Node-Gosper SFC represents an effective point indexing method applicable in those areas.This paper summarizes the background about the hierarchical hexagonal point clustering and the corresponding theory.It shows that the multi-resolution hexagonal grids have several outstanding properties including uniform distance between cells and the best coverage of circular queries.The paper surveys the existing SFC methods and highlights the current problems of the hexagonal clustering including the issues with hierarchical decomposition of hexagons, fractal shape of hexagonal composites and their reasonable indexation.Our method solves all those issues.In comparison with our older approaches, the novel method can deal with jagged edges and it is four times faster.The tests showed that the curve has similar distance properties as widely preferred non-hexagonal SFCs (Hilbert, Sierpinski, etc.) and its construction is even faster.Thus, our Node-Gosper SFC represents a scalable, rapid and stable method for hexagonal hierarchical indexing with good quality parameters.
The future work aims at hierarchical hexagonal structures construction and their application for the practical tasks including data visualization or location indexing in GISs.The discussed metrics and measurements already bring us promising results in areas like nearest neighbors search or range query.

Figure 1 .
Figure 1.Hexagonal clustering: (a) Circular disc query.(b) Points ordered along the Node-Gosper space-filling curve according to the corresponding cluster codes.

Figure 2 .
Figure 2. Row 1: Three iterations of the Gosper curve.Row 2: Tiling with Gosper islands.Seven joined hexagons represent the basic level-1 island.

Figure 7 .
Figure 7.The basic patterns of the Gosper curve: (a) 7-pattern of vertex representation, (b) the second iteration of vertex representation, (c) 7-pattern of node representation, (d) the second iteration of node representation.

Figure 9 .
Figure 9. Indexation and passage order of the Node-Gosper pattern.Figure (a) shows the basic rotation of subordinate patterns.The other two figures represent the patterns with (b) forward (F) and (c) backward (B) indexation.

Figure 11 .
Figure 11.The transformation of directional d i and position c i vectors.

Figure 12 .
Figure 12.The localization of p in level-n Gosper island and the search of the corresponding centers.

Figure 13 .Algorithm 1
Figure 13.Four different hexagonal indexation patterns.The pattern (P1) is the default one.

Figure 14 .
Figure 14.The comparison of (a) top-down and(b) bottom-up level-5 Node-Gosper SFC.The top-down method obviously cannot solve the problem with jagged edges, because it produces non-hexagonal sub-clusters.
Figure 15.The graph shows a comparison of construction times of our Node-Gosper (optimized) method with selected space-filling curves including the previous naive Node-Gosper algorithm[8], top-down Node-Gosper algorithm[7] and simple Node-Gosper (Figure13 P2) algorithm grouping the hexagons into Gosper islands without additional pattern rotation and clusters reordering.Each measurement was performed on a random dataset with 10 6 points for different levels of clustering.

Figure
Figure Cont.

d 1 d 2 d 2 d 1 S≈0.41πr 2 S≈0.64πr 2 S≈0.83πr 2 d 1 ≠d 2 d 1 ≠d 2 d 1 =d 2
[7,8]raph shows total construction times per cluster.Each method has a different number of sub-clusters per pattern (triangular Sierpi ńsky: 2, orthogonal Z-order and Hilbert: 4, hexagonal Node-Gosper: 7), and thus, produces different number of clusters on the same hierarchical level.This metric is computed as (total time/(sub-clusters × curve level)) and it helps to compare the real time complexity of different point indexing methods.The graph compares of our novel Node-Gosper (optimized) method with the older ones[7,8]and other SFCs.Each measurement was performed on a random dataset with 10 6 points for different levels of clustering.

Table 2 .
[59]list of tested datasets with corresponding point number and the length of the bounding box diagonal.The datasets were obtained from website[59].