SFCGen: A Framework for Efficient Generation of ... - ACM Digital Library

9 downloads 0 Views 393KB Size Report
GUOHUA JIN and JOHN MELLOR-CRUMMEY. Rice University. Because they are continuous and self-similar, space-filling curves have been widely used in ...
SFCGen: A Framework for Efficient Generation of Multi-Dimensional Space-Filling Curves by Recursion GUOHUA JIN and JOHN MELLOR-CRUMMEY Rice University

Because they are continuous and self-similar, space-filling curves have been widely used in mathematics to transform multi-dimensional problems into one-dimensional forms. For scientific applications, reordering computation by certain space-filling curves can significantly improve data reuse because of the locality properties of these curves. However, when space-filling curves are used in programs for reordering data, traversal or indexing of the curves must be efficient. To address this problem, we present the table-driven framework SFCGen to efficiently generate multi-dimensional space-filling curves on the fly. The framework is general and easy enough to be used in any application that can be partitioned recursively in multiple dimensions. We describe a movement specification table, a universal turtle algorithm to enumerate points along a space-filling curve, a table-based indexing algorithm to transform coordinates of a point into its position along the curve and an algorithm to pregenerate the table automatically. As examples, we show how high´ dimensional Hilbert, Morton, and Peano curves and a two-dimensional Sierpinski curve can be generated with our algorithms. We present performance results for Hilbert, Morton, and Peano curves and compare the efficiency of our curve generation algorithm with the most recent work on generating Hilbert curves. Our experimental results on three modern microprocessor-based platforms show that SFCGen performs up to 63% faster than the most recent recursive algorithm on 2D curve generation and up to a factor of 132 faster than two previous byte-oriented non-recursive implementations. On curve indexing, SFCGen performs as much as a factor of three faster than the byte-oriented implementation. Our results on 4D space-filling curves also show that SFCGen scales very well with curve level for higher dimensional spaces. Categories and Subject Descriptors: G.4 [Mathematical Software]—Algorithm design and analysis; efficiency; portability; I.3.3 [Computer Graphics]: Picture/Image Generation—Line and curve generation General Terms: Algorithms, Design, Performance Additional Key Words and Phrases: Space-filling curve

This work was supported in part by NCSA under National Science Foundation cooperative agreement ACI-9619019 and the Los Alamos National Laboratory Computer Science Institute (LACSI) through LANL contract number 03891-99-23 as part of the prime contract (W-7405-ENG-36) between the DOE and the Regents of the University of California. Authors’ address: Computer Science Department, Rice University, 6100 Main Street, Houston, TX 77005-1892; email: {jin;johnmc}@cs.rice.edu. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or direct commercial advantage and that copies show this notice on the first page or initial screen of a display along with the full citation. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, to redistribute to lists, or to use any component of this work in other works requires prior specific permission and/or a fee. Permissions may be requested from Publications Dept., ACM, Inc., 1515 Broadway, New York, NY 10036 USA, fax: +1 (212) 869-0481, or [email protected].  C 2005 ACM 0098-3500/05/0300-0120 $5.00 ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005, Pages 120–148.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



121

1. INTRODUCTION Space-filling curves map a compact interval to a multi-dimensional space by passing through every point of the space with positive Jordan content. Since Peano discovered the first space-filling curve over a century ago [Peano 1890], numerous types of space-filling curves have been created and much work has been done in studying mathematical properties and generation of different curves. Hilbert’s work is particularly important. Hilbert was the first to recognize a general geometrical generating procedure that enables the construction of an entire class of space-filling curves [Hilbert 1891]. Since then, spacefilling curves have been used in a variety of fields including mathematics [Butz 1968, 1969], algorithms [Platzman and Bartholdi 1989], geographical information systems [Abel and Mark 1990], image processing [Stevens et al. 1983; Moghaddam et al. 1991; Lamarque and Robert 1996; Velho and Gomes 1991; Zhang and Webber 1993], databases [B¨ohm et al. 2001; Lawder and King 2001; Mokbel and Aref 2001], circuit design [Alpert and Kahng 1994], cryptology [Matias and Shamir 1987; Bertilsson et al. 1989], and scientific computing [Mellor-Crummey et al. 1999; Chatterjee et al. 1999; Jin et al. 2001; Frens and Wise 1997; Kuo et al. 1999; Hu et al. 2000]. They also have been used for reordering both data and computation in regular and irregular problems [MellorCrummey et al. 1999; Chatterjee et al. 1999; Jin et al. 2001; Hu et al. 2000]. Space-filling curves have been studied in various forms. They have interesting locality properties that make them useful for partitioning or reordering data and computation. Their recursive structure can be exploited for generating the curves themselves. Bially [1969] presented an algorithm for converting back and forth between a point in an n-dimensional cube and a number representing a position along a space-filling curve. Bially described the algorithm in a diagrammatic form and gave a numerical procedure to construct the diagrams. Butz [1971] proposed an algorithm for generating a Hilbert curve in a byte-oriented manner. Butz’s algorithm maps a point from a onedimensional space into a point in an n-dimensional space by using circular shift and exclusive-or operations on bytes. Although one can generate a spacefilling curve or Hilbert curve in an arbitrarily high-dimensional space by following Bially’s or Butz’s description, their algorithms are complicated and inefficient. To improve efficiency, many researchers have studied recursive generation of space-filling curves in lower dimensional spaces [Sagan 1992, 1993, 1994; Ohno and Ohyama 1991; Prusinkiewicz et al. 1991; Goldschlager 1981; Cole 1981; Witten and Wyvill 1983; Breinholt and Schierz 1998]. Often, recursive generation of space-filling curves can be expressed in much simpler and more efficient forms than nonrecursive methods. The only limitation they have in common is that it is nontrivial to extend these algorithms to higher dimensional spaces. However, to apply space-filling curves to scientific programs, generation of curves in three or more dimensions is very important. To address this problem, in this article we present SFCGen—a table-driven framework to efficiently generate multi-dimensional self-similar space-filling curves on the fly. Examples of these curves include Hilbert, Peano, Morton ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

122



G. Jin and J. M.-Crummey

Fig. 1. Hilbert and Morton curves at three levels of recursion.

´ ´ [Morton 1966],1 and Sierpinski curves [Sierpinski 1912]. Although our framework is general enough to be used for any self-similar curves, we only describe table generation for three of them, namely, Hilbert, Morton, and Peano curves here, since they have been receiving increasing attention from scientific computing researchers. Recent studies [Abel and Mark 1990; Jagadish 1990; Kumar 1994] also show that Hilbert curve ordering yields better locality and is superior to other orderings for certain application domains. However, the lack of fast algorithms for evaluating a curve position for a cell or subquadrant from its coordinates and vice versa still remains a major weakness. Figure 1 shows 2D and 3D Hilbert and Morton curves.2 Our work differs from previous work in the following aspects. First, we use a universal table-driven turtle algorithm to enumerate points in a multidimensional space along a space-filling curve and a table-based indexing algorithm to transform coordinates of a point into its position along the curve. A movement specification table specifies the direction of movement at each step of the traversal and is pregenerated based on both geometric and arithmetic properties of the curve. Second, we describe how to generate the movement specification tables for high-dimensional Hilbert, Morton, and Peano curves. These precomputed tables enable us to generate space-filling curves in arbitrarily high-dimensional spaces much more efficiently on the fly than previous work. Third, we evaluated the overall performance for traversing and indexing these space-filling curves and compared it with the latest recursive algorithm available and two fast implementations of Butz’s byte-oriented non-recursive algorithm given recently by Moore [2000]. Experimental results on three modern microprocessor-based platforms show that our 2D curve generation algorithm performs up to 63% faster than the latest recursive curve generation algorithm and up to a factor of 132 faster than the two fast byte-oriented nonrecursive implementations. Our table-based indexing algorithm is as much as a factor of three faster than the previous algorithm on these platforms. Our 1 Morton

curves are also referred as Peano SFC by Mokbel et al. [2003]. The type of encoding using these curves is also known as Morton encoding, quad code, bit-interleaving, N-order, locational code or Z-order. 2 SFCGen outputs a list of points forming a space-filling curve. The curves were plotted from the lists of points with MATLABTM from MathWorks Inc. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



123

results also show our curve generation algorithm scales very well with curve level on the three space-filling curves. The rest of the article is organized as follows. In Section 2, we begin with a description of related work. In Section 3, we define a movement specification table and describe our universal turtle algorithm along with the table-based indexing algorithm used in SFCGen. In Section 4, we present an algorithm for generating movement specification tables of 2D and 3D Hilbert, Morton, and Peano curves. In Section 5, we extend our table generation algorithm to n-dimensional space. We report experimental results in Section 6 and conclude with Section 7. 2. RELATED WORK Recursive algorithms for generating space-filling curves have been studied previously by many researchers. Sagan [1992, 1993, 1994] discussed analyticarithmetic representations for various space-filling curves in two and three dimensions and showed how to generate these curves using either arithmetic or geometric constructions. SFCGen also generates curves by recursion. However, SFCGen is table driven. All geometric information for generating a spacefilling curve is included in a movement specification table. Ohno explained how to construct fractal curves through different types of segments [Ohno and Ohyama 1991]. Goldschlager [1981] presented short recursive algorithms that draw space-filling curves in a square region based on specifications of movement directions. Cole, Witten, and Wyvill described recursive procedures for drawing space-filling curves with graphics languages that have constructs representing geometrical transformations [Cole 1981; Witten and Wyvill 1983]. More recently, Breinholt and Schierz [1998] gave an algorithm for generating a 2D Hilbert space-filling curve using integer operations and recursion. In their work, the fundamental Hilbert shape, a line joining the four corners of a square, was represented by two variables with values of either 0 or 1 representing different orientations of the unit shape. The information was manipulated and passed to the next level when the recursive procedure was called. In contrast, SFCGen is a table-driven framework for efficient generation of space-filling curves in arbitrarily high-dimensional space. It uses a universal recursive turtle drawing scheme for elaborating recursive curves in conjunction with a table that defines the shape of primitive curve elements to be drawn. Bartholdi and Goldsman [2001] described a vertex-labeling method to generate algorithms for manipulating a Hilbert curve in 2D space and its application ´ to the generation of 2D Sierpinski and Peano curves. Their approach requires that each cell at a given level of the decomposition maintains a local frame of reference, which consists of vertex labels and codes that imply the spatial relationships of subcells. Labels and codes of subcells are determined by labeling rules. These rules are determined by inspection from a flow diagram, which shows how the curve orders cells at successive levels of the decomposition. It is clear from their experience in labeling a 3D Hilbert curve that their flow-diagram-based approach to generating the labeling rules is not applicable ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

124



G. Jin and J. M.-Crummey

to three or higher dimensional spaces. Rules for higher dimensional Hilbert curves are much more complicated than those for the 2D case. Drakopoulos et al. [1999] used linear fractal interpolation functions to construct space-filling curves in two dimensions. Starting with an initiator, which can be a straight line or a polygon, they replaced each side of the initiator by a connected set of straight lines based on an iterated function system coded in a table consisting of 6-tuples. Each 6-tuple specifies a transformation from the replaced line to a scaled straight line between interpolation points at the next level. The movement specification table we use records higher level information about each primitive curve and is precomputed for higher efficiency. Our table also includes useful information for efficient curve indexing given the coordinates of a point in an arbitrarily high-dimensional space. Prusinkiewicz et al. [1991] described how to use the formalism of L-systems proposed by the biologist Aristid Lindenmayer in 1968 [Lindenmayer 1968] for expressing and drawing FASS curves—curves that are space-Filling, selfAvoiding, Simple and self-Similar. Central to L-systems, is the notion of rewriting. One defines complex objects by successively replacing parts of a simple object using a set of rewriting rules. Prusinkiewicz et al. [1991] describe an algorithm for recursively constructing all connecting patterns suitable for the generation of FASS curves in a square grid using L-system formalisms. The generation of the rewriting rules is based on tile orientation and the connecting pattern between tiles. In addition to forward moves, rotation actions are embedded inside each rule. Their L-system for a 3D Hilbert curve is much more complicated than their 2D system and was constructed by solving the marked cubes puzzle by hand (using a set of appropriately marked wooden toy blocks) [Prusinkiewicz et al. 1991, page 362]. It is not clear how the manual approach they used to construct an L-system for a 3D Hilbert curve could be adapted for higher dimensional curves. Our movement specification table can be easily rewritten in an L-system form. The key difference between the Lsystem in Prusinkiewicz et al. [1991] and our movement specification table is that we specify a rule for each primitive curve; this is more readily extended to higher dimensions than their tile or cube orientation based approach. Since we explicitly enumerate each orientation of a primitive curve in our specification table, we have more table rows than rules needed in an L-system for the same curve. However, our table rows are simpler than L-system rules (and faster to elaborate) since they contain no rotations. Bially [1969] and Butz [1971] discussed generation of space-filling curves in a high-dimensional space. Bially used state diagrams to establish one-to-one mappings between points in spaces of different dimensionalities. Butz used a byte-oriented algorithm to map a binary number representing a point on a space-filling curve into a vector of binary numbers expressing coordinates in a multi-dimensional space. SFCGen is a table-driven framework. It generates curves recursively, based on a movement specification table that is pregenerated. Our first table row for a Hilbert curve is created based on Butz’s algorithm. The rest of the table is automatically generated afterwards. As a result, our curve generation algorithm performs significantly faster than two fast implementations of Butz’s algorithm. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



125

Ken Musgrave [1991] has implemented Butz’s algorithm for generating an n-dimensional Peano curve to m bits of precision. 3. A TABLE-DRIVEN FRAMEWORK In this section, we present SFCGen, a new framework for efficiently traversing and indexing multi-dimensional self-similar space-filling curves that can be constructed recursively. The framework includes a movement specification table for each type of space-filling curve, a universal turtle algorithm walking through the curve and a table-based indexing algorithm. The framework is general enough to be applied to many space-filling curves. We show four types of them in this article. 3.1 Movement Specification Table Geometrically, a level l (l = 2, 3, . . .) self-similar space-filling curve in an ndimensional space can be approximated by a level l − 1 curve, if we can evenly partition the n-dimensional space into (k l −1 )n subcubes and only keep the connection between the subcubes, where k is a refinement factor3 between adjacent levels of the space-filling curve. Similarly, the level l curve can also be approximated by a level 1 curve, if we can evenly partition the entire n-dimensional space into k n subcubes. We call level 1 the base level since it is the lowest level curves can be approximated; we call curves at the base level primitive curves. In our framework, a movement specification table T contains a list of table rows {T (i) | i = 1, . . . , t}. Each table row T (i) recursively defines a refinement j of a curve Cl (i) at the current level l (l = 1, 2, . . .) to a set of curves Cl +1 ( j = 1, . . . , r) of the same type at the1 next level connected with a set of moves M j ( j = M M2 M r−1 1 2 1, . . . , r − 1): Cl (i) : Cl +1 (i) −→ Cl +1 (i) −→ · · · −→ Clr+1 (i). r is the number of table columns for defining the refinement. In our movement specification table, each curve is represented by its approximation at the base level. The movement specification table is defined differently for each space-filling curve. 3.1.1 Hilbert Curves. The Hilbert curves start at one corner of a square and end at one of the two neighboring corners of the square. It is formed by  joining together square cups “ ” oriented in different directions as shown in Figure 2(a). Figure 2(b) shows a Hilbert curve marked with movement direction from the south-west corner to the south-east corner. The Hilbert curve can be approximated at any level greater than 1 by one or multiple square cups properly oriented and connected. There are four undirected primitive curves facing north, south, west, or east at base level as shown in Figure 2(a). Each cup can be traversed in two directions: clockwise and counterclockwise, which lead to a total of eight Hilbert directed primitive curves in 2D. Movement directions between these primitive curves include north, south, west, and east. Table I shows the full movement specification table for the Hilbert curves. The table contains eight table rows, each corresponding to one primitive curve, which is 3 The

refinement factor is a small integer constant that indicates how many instances of a curve at level l are instantiated along each dimension to construct a curve at level l + 1. For example, the refinement factor is 2 for the Hilbert and Morton curves and 3 for the Peano curves. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

126



G. Jin and J. M.-Crummey

Fig. 2. Hilbert curves. Table I. Movement Specification Table for Hilbert Curves Current Level ne sw wn es nw se ws en

Next Level en, ws, nw, se, wn, es, sw, ne,

n s w e n s w e

ne, sw, wn, es, nw, se, ws, en,

e w n s w e s n

ne, sw, wn, es, nw, se, ws, en,

s n e w s n e w

ws, en, se, nw, es, wn, ne, sw,

⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥

named by two letters representing its first two moves. The leftmost column of the table shows an approximation of a space-filling curve at the current level, while the rest of the columns show approximations of four curves at the next level connected by moves. “n”, “s”, “w” and “e” stand for north, south, west and east respectively. “⊥” represents no movement. For example, the Hilbert curve shown in Figure 2(b) is a level 2 curve approximated by an “ne” primitive curve. According to the movement specification in Table I, this curve can be traversed using the recipe given by the row of the table, labeled “ne” in the “Current Level” column. For each entry under the “Next Level” part of the row, we instantiate the primitive curve and connect it to the next primitive curve using the movement direction specified. The instantiation of a row in the table is performed left to right across the row. To instantiate a curve at level l , when traversing a row of the table, simply recursively replace the primitive curve shown in a column entry with the traversal of a row with that label up to a recursion depth of l . Defining a movement specification table for a 3D Hilbert curve is more complicated. First, there are two additional movement directions along the third dimension. Second, there are more primitive curves in 3D and each of them is more sophisticated than the “square cup” primitive curves for 2D. Figure 3(a) shows five of forty-eight different directed primitive curves embedded inside 2 × 2 × 2 cubes. Joining together these primitive curves we can form a 3D Hilbert curve at level 2 as shown in Figure 3(b). As a result, our full movement specification table for 3D Hilbert curves has six times as many table rows as the full table for 2D Hilbert curves. Table II shows only part of the whole table ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.



SFCGen: A Framework for Efficient Generation of Space-Filling Curves

127

Fig. 3. 3D Hilbert directed curves. Table II. Part of the Movement Specification Table for 3D Hilbert Curves Current Level neb ebn bne swb fnw efs

Next Level ebn, bne, neb, wbs, nwf, fse,

n e b s f e

bne, neb, ebn, bsw, wfn, sef,

e b n w n f

bne, neb, ebn, bsw, wfn, sef,

s w f n b w

swb, wfn, fse, neb, bsw, wbs,

b n e b w s

swb, wfn, fse, neb, bsw, wbs,

n e b s f e

fnw, sef, wbs, fse, efs, neb,

w f s e s b

fnw, sef, wbs, fse, efs, neb,

s w f n b w

efs, bsw, nwf, wfn, neb, fnw,

⊥ ⊥ ⊥ ⊥ ⊥ ⊥

describing the formation of the level 2 curve in Figure 3(b). Each primitive curve listed in Table II is named by three letters. The first letter represents the first move. The second letter represents the second move, which is the move between the two parallel lines in the first plane. The last letter represents the fourth move that is the move between the two planes. “e”, “w”, “b”, “f”, “n”, and “s” stand for east, west, back, front, north, and south respectively. 3.1.2 Morton Curves. Movement specification tables for Morton curves are simple for any number of dimensions. A 2D Morton curve starts from one corner of a square and ends at its diagonally opposed corner. Any 2D Morton curve at a level greater than 1 can be approximated by a 2D Morton primitive curve, an “N” or flipped “Z” curve, if it starts from the south-west corner, or multiple copies of the primitive curve diagonally connected based on approximations at lower levels. Figure 4(a) shows a level 2 Morton curve formed by four copies of a primitive curve. Similarly, a 3D Morton curve may start at the south-frontwest corner and end at the north-back-east corner and can be approximated by one or multiple copies of a 3D primitive curve as shown in Figure 4(b). One copy of the primitive curves in both cases is displayed in bold. Table III shows movement specification tables for the 2D and 3D Morton curves.4 Each table contains one table row, since there is only one primitive curve in each case. 4 There

are eight Morton directed primitive curves in two dimensions and forty-eight in three dimensions. Since the set of curves they can approximate are disjoint, we show in Table III only the primitive curves that we used to create the curves in Figure 4. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

128



G. Jin and J. M.-Crummey

Fig. 4. Morton curves at level 2. Table III. Movement Specification Tables for 2D and 3D Morton Curves Current Level ne

Next Level ne, n

Current Level nbe

ne, se

ne, n

ne, ⊥

Next Level nbe, n

nbe, sb

nbe, n

nbe, sfe

nbe, n

nbe, sb

nbe, n

nbe, ⊥

The movement direction “n”, “se”, “sb”, or “sfe” shown in the tables stands for going north, south-east, south-back, or south-front-east after finishing the current subspace (subsquare for 2D or subcube for 3D). Within each primitive curve, the north moves in both 2D and 3D are axial moves. All moves between primitive curves are diagonal (The north move in 2D becomes a north-west move and the north move in 3D becomes a north-front-west move.). In our implementation, we represent each diagonal move with an axial move along the axis of the lowest traversal priority. For example, the “se”, “sb”, and “sfe” moves shown above are represented with “e”, “b”, and “e” respectively. Using this notation enables us to simplify and improve the underlying data structures and the traversal algorithm in the next subsection. 3.1.3 Peano Curves. Peano curves in 2D and 3D spaces are similar to Hilbert and Morton curves. Like a Hilbert curve, all moves are axial in a Peano curve. On the other hand, it starts at one corner and ends at the farthest corner from its origin, a common feature of Morton curves. Peano curves also move along axes in a designated order. However, Peano curves differ from Hilbert and Morton curves in that they have a refinement factor of three. As a result, a 2D or 3D Peano primitive curve connects nine points in a 3 × 3 square and twenty-seven points in a 3 × 3 × 3 cube. Figure 5(a) shows one of the eight 2D Peano directed primitive curves: an ne curve, inside a 3 × 3 square. An ne curve moves towards the north-east corner from the south-west corner. It moves from west to east by columns and from south to north or from north to south interleavingly within each column. A 2D Peano curve at any level of refinement greater than 1 can be approximated by a 2D Peano primitive curve, or multiple primitive curves properly connected. There are forty-eight ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.



SFCGen: A Framework for Efficient Generation of Space-Filling Curves

129

Fig. 5. Peano directed primitive curves. Table IV. Movement Specification Tables for 2D and 3D Peano Curves Current Level ne nw se sw

Next Level ne, n nw, n se, s sw, s

nw, n ne, n sw, s se, s

ne, e nw, w se, e sw, w

se, s sw, s ne, n nw, n

Current Level neb

nwf

sef

swb

sw, s se, s nw, n ne, n

se, e sw, w ne, e nw, w

ne, n nw, n se, s sw, s

nw, n ne, n sw, s se, s

ne, ⊥ nw, ⊥ se, ⊥ sw, ⊥

Next Level neb, n swb, s neb, n nwf, n sef, s nwf, n sef, s nwf, n sef, s swb, s neb, n swb, s

nwf, n sef, s nwf, n neb, n swb, s neb, n swb, s neb, n swb, s sef, s nwf, n sef, s

neb, e swb, w neb, e nwf, w sef, e nwf, w sef, e nwf, w sef, e swb, w neb, e swb, w

sef, s nwf, n sef, s swb, s neb, n swb, s neb, n swb, s neb, n nwf, n sef, s nwf, n

swb, s neb, n swb, s sef, s nwf, n sef, s nwf, n sef, s nwf, n neb, n swb, s neb, n

sef, e nwf, w sef, e swb, w neb, e swb, w neb, e swb, w neb, e nwf, w sef, e nwf, w

neb, n swb, s neb, n nwf, n sef, s nwf, n sef, s nwf, n sef, s swb, s neb, n swb, s

nwf, n sef, s nwf, n neb, n swb, s neb, n swb, s neb, n swb, s sef, s nwf, n sef, s

neb, b swb, b neb, ⊥ nwf, f sef, f nwf, ⊥ sef, f nwf, f sef, ⊥ swb, b neb, b swb, ⊥

3D Peano directed primitive curves moving in different axial directions or orders. Figure 5(b) shows one of them, an neb curve, inside a 3 × 3 × 3 cube. An neb curve moves towards the east-back-north corner from the west-front-south corner in a 3D space. As it moves from the front plane to the back plane, it moves through the first plane from west to east, the middle plane east to west, and the back plane west to east again. Like a 2D curve, a 3D Peano curve at any level of refinement greater than 1 can be formed or approximated by one or multiple primitive curves. Movement specification tables for generating an ne curve in 2D space and neb curve in 3D space by recursion are described in Table IV. The tables only show the table rows that are required to generate the two specific curves. Each table row of the 2D specification table includes nine table entries, while each table row of the 3D specification table contains twenty-seven entries. (For the 3D table, a stack of three “Next Level” rows should be read left to right and top to bottom to form a sequence of twenty-seven entries.) ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

130



G. Jin and J. M.-Crummey

´ Fig. 6. 2D Sierpinski curves. Table V. Movement Specification Table for 2D ´ Sierpinski Curves Current Level w e s n

Next Level w, e, s, n,

sw ne se nw

s, n, e, w,

w e s n

n, s, w, e,

nw se sw ne

w, e, s, n,

⊥ ⊥ ⊥ ⊥

´ 3.1.4 Sierpinski ´ Curves. A Sierpinski curve differs from the other three types of curves in several ways. First, it is a closed curve—it walks through every point of a square and returns to its origin. Second, it walks through one half (isosceles right triangles) of each square first and then the other half. Third, moves can be axial or diagonal with different distances. To generate a ´ Sierpinski curve, we divide a square into four isosceles right triangles. Each of the isosceles right triangles is further divided into four smaller isosceles right triangles. We call the curve segment in each of these triangles Sierpinski ´ quarter curve and generate it based on a movement specification table. ´ Figure 6(a) shows four undirected Sierpinski primitive quarter curves inside ´ a 2 × 2 square and Figure 6(b) shows a level 2 Sierpinski curve formed by four quarter curves connected with dotted lines. “w”, “e”, “s”, and “n” represent quarter curves moving towards west, east, south, and north. “sw”, “ne”, “se”, and “nw” are diagonal moves towards south-west, north-east, south-east, and north-west. 3.1.5 Multi-level Specification Table. Our framework can also generate variants of different space-filling curves using a multi-level table strategy. Similar to the one-level tables, each level of our multi-level movement specification table defines a curve refinement. However, with a multi-level table, we can specify different refinements at different levels. Figure 7(a) shows a 2D Morton curve at level 4, with the three outer levels instantiated as “N” curves and the innermost level instantiated as flipped “Z” curves. Figure 7(b) shows two different types of Peano curves. The first Peano curve is created with ne curves at outer levels and en curves at the innermost level. The second Peano curve is created by using Table IV for the outer levels and Table VI for the innermost level. Each move in Table IV is further instantiated as a curve defined interleavingly as the first or second half of the corresponding table row in Table VI. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



131

Fig. 7. Space-filling curves constructed with multi-level movement specification tables.

Table VI. A Level Two Movement Specification Table for Peano Curves Current Level w e s n

Next Level w, e, s, n,

nw ne se ne

w, e, s, n,

w e s n

w, e, s, n,

sw se sw nw

w, e, s, n,

sw se sw nw

w, e, s, n,

w e s n

w, e, s, n,

nw ne se ne

Fig. 8. Universal turtle algorithm.

3.2 A Turtle Traversal Algorithm A key part of our table-driven framework is the universal turtle traversal algorithm, which traverses a self-similar space-filling curve based on a movement specification table. The algorithm is general enough to be used for generating a variety of curves and efficient enough to be integrated into applications so that computation is performed in a specific curve ordering. Figure 8 shows a brief description of the algorithm. To traverse a curve of n-dimensions, it starts  and ex representing the integer extent along from the n-dimensional point 0 each dimension that the curve will cover. The initial value of the table index ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

132



G. Jin and J. M.-Crummey

variable row should point to a table row matching the primitive curve approximation of the curve. The algorithm first checks if it has reached the last recursion level by calling IsBaseCase. If true, it calls BaseAction to print the primitive curve or perform certain calculations by calling the routine prog, otherwise, it divides the current space into subspaces and calls Uturtle for each subspace recursively. Move connects or moves from an ending point of a curve in one subspace to a starting point of another curve in the next subspace. Our algorithm is more general than previous algorithms in several aspects. First, the algorithm itself is curve- and dimensionality-independent. Second, each dimension may have a different extent. Third, the number of recursion levels can be controlled by adjusting the subspace size at the last recursion level. To improve efficiency of the algorithm, each move is described by a 2-tuple (dim, distance) instead of a distance vector, where d im represents the axial direction of the move or the axial direction of the lowest traversal priority of a diagonal move, and distance represents a movement distance that can be a unit forward or backward along the axis for an axial move or axes for a diagonal move. This axial movement specification reduces cost from unnecessary element-wise checking and calculation for each move, particularly at the last recursion level. The representation naturally describes moves in many space-filling curves including the Hilbert and the Peano curves. It substantially reduces the low-level ´ operations for other curves including the Morton and the Sierpinski curves. In Section 6, we present performance measurements for traversing a Hilbert, Peano, or Morton curve. 3.3 A Curve Indexing Algorithm In the previous subsection, we presented a universal turtle algorithm for traversing a space-filling curve based on a movement specification table. By recording the relative position of each part of the curve at a particular recursion level, the movement specification table can be used for other purposes as well. In this section, we describe how to use the table to find the position of a point along a curve given the coordinates of the point. Consider the level 2 Hilbert curve shown in Figure 2(b) as an example. Assume the coordinates of the starting point are (0, 0). To find the position along the curve of the point at the north-east corner with coordinates (3, 3), we start from the table row for “ne” and initialize the position variable pos as 0. From the coordinates (3, 3) and the table row, we know that the point is located in the third subsquare (the third point in the level 1 approximation). We adjust the position by 8 since each of the subsquares has four points. Next, we check the curve at the next level (level 2). The point at (3, 3) is the third point within the north-east subsquare. We update the position by adding 3 to it and find the position of the point: 11. The complete algorithm for computing the position along a space-filling curve from a point’s coordinates using a movement specification table is shown in Figure 9. Inputs to the algorithm include level and row in addition to the coordinates and the table. For each recursion level, it computes a position offset and finds a new table row for the next level based on information from the current table row and the coordinates at the current ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



133

Fig. 9. Indexing algorithm.

level. Preliminary results from evaluating the algorithm are presented in Section 6. 4. GENERATION OF MOVEMENT SPECIFICATION TABLES In this section, we describe how to automatically generate movement specification tables for Hilbert, Morton, and Peano curves. For Morton curves, the movement specification table has only one table row, which defines the primitive curve with designated axial movement direction and traversal priority order of the axes. The primitive curve approximation and moves within and between multiple copies of the primitive curve are determined by the traversal priority order of the axes according to the following rule: always move along the axes with higher traversal priorities before the axes with lower traversal priorities. Moving along the axis with the highest traversal priority stands for moving along the axis in its designated direction within a primitive curve. Moving along any of the other axes requires moving diagonally across the axes, or more precisely, moving along the axis in its designated direction and along all the axes with lower traversal priorities in directions opposite to their designated directions. Moving between primitive curves requires moving along the axis with the highest traversal priority in its designated direction and moving along all other axes in the directions opposite to their designated directions. Generating the movement specification table for Hilbert and Peano curves is more complicated and is discussed in the rest of this section. Our table generation process for Hilbert and Peano curves includes two steps: creating an initial table row and generating the rest of the table rows by transformations. 4.1 Defining Geometric Transformations In this subsection, we define two basic geometric transformations, namely rotation and reflection. We use these transformations to create a new curve or movement specification from an existing one. 4.1.1 Rotation. Rotation in a multi-dimensional space is a plane rotation and can be defined by a rotation matrix R. Below is a rotation matrix ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

134



G. Jin and J. M.-Crummey

representing the xy-plane rotation in 2D.   cos θ − sin θ x,y R2 (θ) = sin θ cos θ There are three 3D rotation matrices representing the xy-, yz-, and zx-plane rotations as shown below. In each of the matrices, θ represents a rotation degree according to the following right hand rule: point your right thumb into the positive direction of the axis perpendicular to the plane of rotation and curl your fingers into the forward direction of the rotation.   cos θ − sin θ 0   x,y cos θ 0  R3 (θ) =  sin θ 0 0 1 

   1 0 0 cos θ 0 sin θ     y,z 0 1 0 . R3 (θ ) =  0 cos θ − sin θ  , R3z,x (θ) =  0 sin θ cos θ − sin θ 0 cos θ 4.1.2 Reflection. We are interested in two types of reflections in 2D and 3D: reflection through an axis and reflection through a plane. Reflection through an axis in 2D or 3D is defined by one of the following reflection matrices:     1 0 −1 0 x y R2 = , R2 = 0 −1 0 1 

     1 0 0 −1 0 0 −1 0 0 x y z       R 3 =  0 −1 0  , R 3 =  0 1 0  , R 3 =  0 −1 0  . 0 0 −1 0 0 −1 0 0 1 Similarly, we define a reflection through xy-, yz- or zx-plane in 3D by one of the following matrices.       1 0 0 −1 0 0 1 0 0 x,y y,z z,x       R 3 =  0 1 0  , R 3 =  0 1 0  , R 3 =  0 −1 0  . 0 0 −1 0 0 1 0 0 1 Reflection through any axis can also be accomplished by two reflections through the planes that the axis lies in. We use reflections for generating Hilbert and Peano curves and rotations for Hilbert curves. We only use 90◦ rotations in our table generation for Hilbert curves. 4.2 Changing a Movement by Transformations To determine the movement direction from point a to point b, we compute the  = cb − ca . The relationship difference between their coordinates ca and cb,  between movement directions and the coordinate differences for 2D and 3D ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



135

cases are east or west if x > 0 or x < 0,  y = 0; north or south if  y > 0 or  y < 0, x = 0; and east or west if x > 0 or x < 0,  y = 0, and z = 0; north or south if z > 0 or z < 0, x = 0, and  y = 0; back or front if  y > 0 or  y < 0, x = 0, and z = 0. Transformations change movement directions by changing the difference t =   t . For example, in two dimensions, an east between coordinates R  move is changed to a north move after a 90◦ xy-plane rotation since x,y x,y (x , y )t = R2 (90◦ )(x ,  y )t = R2 (90◦ )(1, 0)t = (0, 1)t . In 3D space, an east move is not changed after a 90◦ yz-plane rotation since (x , y , z )t = y,z y,z R3 (90◦ )(x ,  y , z )t = R3 (90◦ )(1, 0, 0)t = (1, 0, 0)t . After applying reflection through the y axis, an original east move in two dimensions is changed into y a west move since (x , y )t = R 2 (1, 0)t = (−1, 0)t . Similarly, we can compute how other movement directions are changed by a rotation or reflection transformation. 4.3 Creating an Initial Table Row Creating an initial table row in our movement specification table is important for constructing a Hilbert or Peano curve since the rest of the table is generated based on the initial table row by transformations. 4.3.1 Hilbert Curve. An initial table row for Hilbert curves is created using the byte-oriented algorithm by Butz [1971]. The algorithm maps a point r from the one-dimensional space R1N into a point a in the n-dimensional space Rnm in the Hilbert curve ordering, where N and m are the numbers of binary number bits needed to represent points in the one and n-dimensional space, respectively, and N = mn. We assume m = 2, since each table row involves only two recursion levels. We discuss the creation of an initial table row in threedimensional space (n = 3). According to Butz [1971], the point r is represented by a binary number: r = .ρ 1 ρ 2 = .ρ11 ρ21 ρ31 ρ12 ρ22 ρ32 , while the point a is a vector of three elements a1 , a2 , a3 , with each of them expressed as a 2-bit binary number a j = .α 1j α 2j ( j = 1, 2, 3). Let α i = α1i α2i α3i , where (i = 1, 2). We have



α 1 = ρ11 ρ21 ⊕ ρ11 ρ31 ⊕ ρ21 ; “⊕” represents bitwise exclusive-or operation. The value of α 2 is computed as follows based on the principal position5 J of ρ 1 , which is the last position j in ρ 1 such that ρ 1j = ρ31 ( j = 1, 2). 

1

1

1 1 2 1 2 2 1 2 2    ρ1 ⊕ ρ3 ⊕ ρ1 ρ2 ⊕ ρ1 ⊕ ρ2 ⊕ ρ1 ρ3 ⊕ ρ2 ⊕ ρ3 ⊕ ρ2 if J = 1 α2 = ρ11 ⊕ ρ32 ⊕ ρ22 ρ21 ⊕ ρ11 ⊕ ρ31 ⊕ ρ12 ρ31 ⊕ ρ21 ⊕ ρ22 ⊕ ρ12 if J = 2  





 1 ρ1 ⊕ ρ22 ⊕ ρ12 ρ21 ⊕ ρ11 ⊕ ρ32 ⊕ ρ22 ρ31 ⊕ ρ21 ⊕ ρ31 ⊕ ρ12 if J = 3. 5 Since

m = 2, we use J to represent J1 in Butz’s notation [Butz 1971]. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

136



G. Jin and J. M.-Crummey

J is 3 if all bits of ρ 1 are equal. An overbar stands for bitwise negation. Due to the depth-first recursive nature of the Hilbert curve construction, points within each subcube at level l will be referenced consecutively after its predecessor subcubes at level l , but before its successor subcubes at level l . This implies that points within the same subcube have the same value of ρ 1 . Values of ρ 2 determine the order in which points within a subcube are visited. For example, ρ 2 = 000 represents the first point accessed in a subcube while ρ 2 = 100 starts the second half points in the subcube. To determine the primitive curve for each of the eight subcubes, we need to find movement direction from the first point to the second for each of the following three pairs of points ρ11 ρ21 ρ31 0 0 0 and ρ11 ρ21 ρ31 0 0 1 ρ11 ρ21 ρ31 0 0 1 and ρ11 ρ21 ρ31 0 1 0

(∗)

ρ11 ρ21 ρ31 0 1 1 and ρ11 ρ21 ρ31 1 0 0. As an example we here show how to determine the primitive curve for the first subcube. We have ρ 1 = 000 for points in the first subcube. The principal position J for ρ 1 is 3. So, the previous formulas for α 1 and α 2 can be simplified. Then, a1 , a2 , and a3 can be computed based on α 1 and α 2 .



α 1 = 0 0 0, α 2 = ρ22 ⊕ ρ12 ρ32 ⊕ ρ22 ρ12



a1 = 0 ρ22 ⊕ ρ12 , a2 = 0 ρ32 ⊕ ρ22 , a3 = 0 ρ12 . In the first pair of points in (∗), only ρ32 is flipped. ρ12 and ρ22 are not changed. Therefore, a1 and a3 are not changed. a2 is increased by 1 from 00 to 01. This is a back move within the first yz-plane. For the second pair of points, ρ22 and ρ32 are flipped, but ρ12 is not changed. So, a2 and a3 are not changed, while a1 is increased by 1 from 00 to 01. From this change, we learn that it is a north move in the west yz-plane. Finally, in the third pair of points, ρ12 , ρ22 , and ρ32 are all flipped. Therefore, a1 and a2 are not changed and a3 is increased by 1 from 00 to 01. This change of a implies a move from the west plane to the east plane in the subcube. As a result, we recognize it as a “bne” primitive curve in the first subcube. Similarly, we can find the primitive curves in the other seven subcubes and the movement directions between the subcubes. 4.3.2 Peano Curve. Creating an initial table row for a Peano curve is relatively easy. The primitive curve and its definition are determined by the designated axial movement directions and traversal priority order of the axes. The primitive curve in the first subspace at the next level is the same as the primitive curve at the current level. Any of the primitive curves in the rest of the subspaces is generated from its immediate predecessor by a reflection transformation through the axis along which a move flows from the predecessor curve. Moves between primitive curves are determined as follows: always move along the axes with higher traversal priorities before the axes with lower traversal priorities, move along each of the axes (except the axis with the lowest traversal priority) first in its designated axial direction and then the opposite direction, in an interleaving fashion. For example, to create the table entries under the ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



137

“Next Level” part of the first row “ne” in Table IV, we generate the same primitive curve in the first table entry, then an “nw” curve in the second table entry by a reflection transformation through y axis, an “ne” curve in the third table entry from the “nw” curve in the second table entry by the same reflection, and then similarly the curves in the rest of the table entries by reflections. Between the primitive curves, we move north twice, then east, then south twice, followed by an east move and two north moves. 4.4 Generating Table Rows by Transformations After we create an initial table row for one primitive curve, we generate the rest of the table by transforming the initial primitive curve and other curves generated previously during the process. The impact of rotation and reflection transformations on a primitive curve or a primitive curve approximation can be inferred by changes to movement directions and relative coordinate positions of the curve. For example, a 2D sw Hilbert curve consists of a south and a west move followed by a north move. Rotating the xy-plane by 90◦ produces the following coordinate changes   1x 2x 3x = 1 y 2 y 3 y  x,y R2 (90◦ )

1x 2x 3x 1 y 2 y 3 y



 =

x,y R2 (90◦ )

0 −1 0 −1 0 1



 =

1 0 −1 0 −1 0

 ,

where each column of the coordinate matrix describes a move along the primitive curve. The curve has been changed to an east, south, west curve, which is an es curve based on the specification shown in Table I. Similarly, applying a 3D rotation to a 3D Hilbert primitive curve may generate a new 3D primitive curve. We apply reflection transformations through axes and planes to both Hilbert and Peano curves. For example, by applying a reflection through the z axis, we can transform the first table row that defines the “neb” Peano curve in Table IV into the second table row for “nwf” curve. To generate the movement specification table from an initial table row, we used the demand-driven table generation algorithm shown in Figure 10. To build the table more efficiently, the algorithm only creates curve specifications that are actually needed during the recursive curve generation. For curves of three or more dimensions, many primitive curves are not needed because they will never be used if only one orientation for the outermost level curve is always used. TableGen starts from an initial table row created in Section 4.3 as the initial set of existing table rows Sin . For each table row in Sin , it checks the definition of the curve called Cthis in each table entry. If no curve specification has been created for a curve Cnext in the definition, it computes the transformations that are needed to transform the known curve Cthis to Cnext . Applying these transformations to the definition of Cthis will produce the definition of Cnext . Finally, it creates a new table row by inserting Cnext and its definition to the table and adds Cnext into a newly created table row set Sout . It switches the input and ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

138



G. Jin and J. M.-Crummey

Fig. 10. Demand-driven table generation algorithm.

output table row sets and resets the output set after checking all the existing table rows in the input set Sin . 5. HIGHER DIMENSIONAL CURVES So far, we have presented the movement specification table and table generation for 2D and 3D space-filling curves in Sections 3 and 4, respectively, because their construction for these dimensional spaces is intuitive. Our table-driven framework can also be applied to create tables for space-filling curves of higher dimensionality. In this section, we extend the representation of movement specification and table generation to describe Morton, Hilbert, and Peano curves6 and their formations in n-dimensional space more precisely though less intuitively. 5.1 Extension of Movement Specification Table Assume d 1 , d 2 , . . . , d n are the n axes in n-dimensional space. A move M along axis d i is represented by d i+ or d i− in the movement specification table if M increases or decreases the value of d i .7 Primitive curves are represented by the 6 We

´ are unaware of a specification for n-dimensional Sierpinski curves. simplicity, we use a uni-axial direction to represent each move. This presentation matches well with curves moving along only one axis in each step such as a Hilbert or Peano curve. For ´ others with diagonal moves such as a Morton curve, a Sierpinski curve, or the variant type of the Peano curve shown in Figure 7(b), each move is represented by the movement direction along the axis with the lowest traversal priority. Movement distance is not specified for each move here. It is ´ assumed as a unit size along each of the axes involved in the move. For a Sierpinski curve, or the second Peano curve in Figure 7(b), we use two different unit sizes. 7 For

ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



139

minimum number of moves sufficient to distinguish each individual primitive curve. An n-dimensional Morton, Hilbert, or Peano primitive curve is represented by n moves d i∗1 d i∗2 . . . d i∗n , where d i∗j describes the move from the k j −1 -th node to its immediate successor node and k is the refinement factor of the curve. “∗” can be “+” or “−”. With this notation, the 3D Hilbert primitive curve “bne” in the previous example will be represented by “ y + z + x + ” or “d 2+ d 3+ d 1+ ”—meaning first move along dimension 2, then along dimension 3 and then along the first dimension. 5.2 Extension of Table Generation Table generation for Morton and Peano curves as discussed in the previous section can be directly applied to n-dimensional space. Here, we only describe how to extend the work presented in Section 4.3 and 4.4 to generate a movement specification table for an n-dimensional Hilbert curve. As discussed in Section 4.3, our initial table row for the Hilbert curve is created based on Butz’s byte-oriented algorithm [Butz 1971]. Again, we only need to consider m = 2 because of the recursive property of the curves. Assume binary representations of points r ∈ R1N and a ∈ Rnm are: r = .ρ 1 ρ 2 = .ρ11 ρ21 . . . ρn1 ρ12 ρ22 . . . ρn2 and a = (a1 , a2 , . . . , an ), where a j = .α 1j α 2j (1 ≤ j ≤ n). According to Butz [1971], we have the following formulas



1 α 1 = ρ11 ρ21 ⊕ ρ11 . . . ρn1 ⊕ ρn−1       1 2 1 2  ρ11 ⊕ ρn1 ⊕ ρ12 . . . ρi1 ⊕ ρi−1 ⊕ ρi2 ⊕ ρi−1 ⊕ ρn2 ⊕ ρn−1 . . . ρn1 ⊕ ρn−1 if J = 1           1 2 1  ⊕ ρi+1 ⊕ ρi2 . . . ρn1 ⊕ ρn−1 ⊕ ρn1 ⊕ ρ12 if J = n  ρ11 ⊕ ρ22 ⊕ ρ12 . . . ρi1 ⊕ ρi−1      2 2 1 1 2 2 α2 = ρ11 ⊕ ρn−J +2 ⊕ ρn−J +1 . . . ρi  ⊕ ρi  −1 ⊕ ρn−J +i  +1 ⊕ ρn−J +i        if 1 < J < n . . . ρ 1J ⊕ ρ 1J −1 ⊕ ρn1 ⊕ ρ12 . . .         1  1 2 2  ⊕ ρn−J ρi ⊕ ρi1 −1 ⊕ ρi2 −J +1 ⊕ ρi2 −J . . . ρn1 ⊕ ρn−1 +1 ⊕ ρn−J ,

where 2 ≤ i ≤ n − 1, 2 ≤ i  ≤ J − 1, J + 1 ≤ i  ≤ n − 1, and J is the principal position of ρ 1 . Each hypercube at the current level is partitioned into 2n subcubes at the next level represented by ρ 1 . Positions within the same subcube at the next level have the same values of ρ 1 . Primitive curve approximation in each of the subcubes is determined by examining the movement direction from the first point to the second for each of the following n pairs of points: ρ11 . . . ρn1 0 0 . . . 0 and ρ11 . . . ρn1 0 0 . . . 1 ρ11 . . . ρn1 0 0 . . . 1 and ρ11 . . . ρn1 0 . . . 1 0 ...... 1 1 ρ1 . . . ρn 0 1 . . . 1 and ρ11 . . . ρn1 1 0 . . . 0. Movement directions between the primitive curve approximations at the next level can be found from the moves out of points ρ 1 11 . . . 1 with different values of ρ 1 . To generate the rest of the table from the initial table row, we define plane rotation and reflection in n-dimensional space by the following rotation and ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

140



G. Jin and J. M.-Crummey

reflection transformation matrices R and R    ri,i = 1   r = cos θ, r = cos θ  b,b  a,a  a,b Rn (θ ) = ri, j   ra,b = − sin θ, rb,a = sin θ   r =0 i, j  Q  R n = r i, j

  r i,i = 1    r i,i = −1   r i, j = 0

i = a, b

elsewhere

i∈Q i ∈ Q elsewhere,

where 1 ≤ i, j, a, b ≤ n. Rna,b represents an ab-plane rotation and θ is a rotation Q degree. R n is a reflection through a k-dimensional space Q, where 1 ≤ k ≤ n − 1. By applying the table generation algorithm presented in Section 4.4 with this extended matrix representation, we can generate a movement specification table for a Hilbert or Peano curve in an arbitrarily high-dimensional space. 5.3 4D Space-Filling Curves By using the table generation and the universal turtle algorithms, we can generate higher dimensional space-filling curves as easily and efficiently as their lower dimensional counterparts. However, presenting even a 4D curve is much more challenging, especially rendering the 4th dimension in an intuitively understandable way. In this section, we present Hilbert, Morton, and Peano curves in a 4D space. As with the other curves shown in this article, we used MATLAB for drawing them. Our curves position along the 4th dimension by contracting the three other dimensions towards a point at the center of the curve. Figure 11(a) shows the 4D Hilbert, Morton, and Peano primitive curves. The interior portions of the Hilbert and Morton curves are presented with dotted lines, while in the Peano curve, the two levels of interior portions are presented in dashed lines and dotted lines. In MATLAB, we plotted the coordinates in the first three dimensions by contracting them 40% towards the center of the 3D space curve for each unit increment along the fourth dimension. Similarly, in Figure 11(b), we present the 4D curves at level 2 with a contraction rate of 25% inward. Although it is impossible to follow the flow of the curves at level 2, the purpose of the figures is to provide a sense of these 4D curves. 6. EXPERIMENTAL RESULTS We evaluated the performance of our universal curve generation algorithm for Hilbert, Morton, and Peano curves, and the performance of our indexing algorithm on Hilbert curves. We conducted experiments on three modern microprocessor-based platforms: an SGI Origin 2000, an Intel Pentium 4 PC, and an Intel Itanium 2 workstation. Table VII summarizes the characteristics of each platform. The SGI origin 2000 we used contains 16 MIPS R12000 processors. Each R12000 processor can issue four instructions per clock cycle. There are two floating-point units and two integer ALUs(Arithmetic/Logic Units) for ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



141

Fig. 11. 4D Hilbert (leftmost), Morton (center), Peano (rightmost) curves.

Table VII. Machine Characteristics Machine SGI Origin 2000 Intel Pentium 4 Intel Itanium 2

Clock Rate

Cache Sizes

Cache Associativities

Miss Latencies (Cycles)

300 MHz 2.4 GHz 900 MHz

32KB/8MB 8KB/512KB 16KB/256KB/1.5MB

2/2 4/8 4/8/12

8/100 18/92 5/12/120+

address computation and integer operations. The Pentium 4 processor decodes instructions at the maximum rate of one instruction per cycle. There are one floating-point unit and three ALUs for float-pointing and integer operations. The Itanium 2 processor has six integer and four floating-point units. The processor can issue up to six instructions per cycle. For the Hilbert curve, we compared the performance of our universal curve generation algorithm to the latest recursive algorithm and Butz’s non-recursive algorithm. We also compared the performance of our indexing algorithm with a recent implementation of Butz’s byte-oriented algorithm. Our algorithms were implemented in C++, whereas the other two were implemented in C. All programs are compiled with the vendors’ compilers. Table VIII lists the C/C++ compilers and optimization options we used in the experiments. For each executable, we collected the number of graduated instructions and the execution time using hardware performance counters. We used ssrun on the SGI origin 2000 and PAPI [Browne et al. 2000] on the Pentium 4 and Itanium 2 to access the hardware counters. Each experiment was repeated multiple times. Variations between different executions were small and average measurements are presented. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

142



G. Jin and J. M.-Crummey Table VIII. Compilers and Optimizations Machine SGI Origin 2000 Intel Pentium 4 Intel Itanium 2

Compiler

Compiler Optimizations

MIPSpro cc/CC 7.3.1.3m Gnu C++ V2.96 Intel C++ V7.1

−O3 −O3 −O3

Table IX. Execution Metrics for Traversing a 2D Hilbert Curve at Different Levels in Millions of Cycles and Graduated Instructions Breinholt & Schierz Origin 2000 Cycles GInst

Intel Pentium 4 Cycles GInst

Itanium 2 Cycles GInst

6 8 10 12 14 Curve Level

0.20 2.59 41.0 655 10487

0.23 3.51 56.0 894 14311

0.19 0.17 3.00 2.73 46.8 43.7 749 699 11929 11185 Butz & Moore

0.09 1.49 23.8 380 6086

0.19 3.02 48.2 772 12348

6 8 10 12 14 Curve Level

1.43 25.4 458 8147 143K

1.90 35.5 661 11957 213K

6.06 4.08 112 79.3 2140 1512 38836 27778 695K 502K Butz & Moore (incr)

1.50 27.4 500 8883 156K

2.77 52.1 967 17461 311K

6 8 10 12 14 Curve Level

1.31 25.9 506 9584 177K

2.10 43.5 851 16149 298K

9.53 5.85 200 123 3975 2422 75671 46150 1403K 857K SFCGen

1.86 38.6 759 14423 267K

3.12 64.7 1273 24155 447K

6 8 10 12 14

0.16 1.87 29.0 464 7402

0.16 2.65 42.3 676 10823

0.16 2.63 41.0 658 10461

0.06 0.91 14.6 233 3727

0.14 2.16 34.6 554 8858

Curve Level

0.16 2.64 42.3 677 10827

6.1 Performance of Traversal Table IX shows execution cost of traversing 2D Hilbert curves at different levels using various algorithms. Results shown on the top are from Breinholt and Schierz’s recursive algorithm [Breinholt and Schierz 1998]. Results shown at the bottom are from our universal curve generation algorithm. Results shown in the middle of the table are from two implementations of Butz’s byte-oriented algorithm [Butz 1971] by Doug Moore at Rice University [Moore 2000]: one computes coordinates of a point from an index number representing the order of the point along a curve and the other, denoted by “incr”, computes coordinates of each point incrementally from its predecessor’s coordinates. Columns marked as “Cycles” and “GInst” represent millions of execution cycles and graduated instructions respectively. As the table shows, our algorithm uses many fewer graduated instructions than Breinholt and Schierz’s algorithm on Origin and Itanium processors. On Pentium 4, we see a very small difference of ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.



SFCGen: A Framework for Efficient Generation of Space-Filling Curves

143

Table X. Execution Time for Traversing 3D Hilbert Curves at Different Levels Butz & Moore Origin 2000 Cycles GInst

Intel Pentium 4 Cycles GInst

Itanium 2 Cycles GInst

4 5 6 7 8 9 Curve Level

1.30 11.1 97.9 844 7151 61.6K

1.70 15.6 138 1213 10363 90.8K

5.42 3.62 49.8 33.7 437 304 3841 2698 33402 23740 298K 210K SFCGen

1.30 11.7 103 892 7648 66.4K

2.53 22.9 202 1767 15351 134K

4 5 6 7 8 9

0.18 0.92 6.89 54.6 436 3488

0.16 1.28 10.1 80.3 642 5137

0.15 1.07 8.52 67.9 538 4280

0.05 0.57 2.63 36.0 168 2300

0.1 1.44 6.35 92.1 407 5893

Curve Level

0.14 1.15 9.18 73.4 587 4698

instruction counts between Breinholt and Schierz’s algorithm and ours. However, we achieve a 14% to 19% performance improvement over their algorithm mainly because our table-driven algorithm uses simpler operations. Overall, our algorithm is up to 63% faster than Breinholt and Schierz’s recursive algorithm. In comparison with Butz’s non-recursive byte-oriented algorithm, our algorithm uses only 1% to 8% of the number of instructions their algorithm uses. Computing principal position and calculating exclusive-or operations tend to be very expensive. Overall, our algorithm improves Moore’s fast implementations of Butz’s byte-oriented algorithm from a factor of 20 to a factor of 132. Both can generate Hilbert curves in arbitrarily high-dimensional spaces. Note that the implementation that computes coordinates incrementally performs worse than the other version, which computes coordinates from an indexing number, although the former seems to be a more reasonable way to traverse a curve. Table X shows performance metrics for traversing a 3D Hilbert curve. We present results only from the faster implementation of Butz’s algorithm and ours. Breinholt and Schierz’s algorithm can only generate a 2D Hilbert curve. The results show that our algorithm performs a factor of 9 to 65 times faster than Butz’s algorithm for traversing a 3D Hilbert curve. In addition to high-dimensional Hilbert curves, our algorithm can also generate other self-similar space-filling curves. In Table XI, we report results of SFCGen for traversing 3D Morton curves at different levels. Overall, generating a Morton curve is more expensive than generating a Hilbert curve at the same level with our algorithm because the diagonal moves in the Morton curve need more calculations than the orthogonal unit-size moves along the axes in the Hilbert curve. 6.2 Performance of Indexing We also measured performance of our table-based indexing algorithm. Table XII shows the execution metrics for translating tuples of coordinates into positions along a 2D Hilbert curve of different levels. We compared performance of ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.



144

G. Jin and J. M.-Crummey

Table XI. Execution Time of SFCGen for Traversing 3D Morton Curves at Different Levels Level 4 5 6 7 8 9

Origin 2000 Cycles GInst 0.23 1.36 10.6 84.5 674 5395

0.23 1.67 13.2 105 844 6749

G I nst C yc

Cycles

Pentium 4 GInst

G I nst C yc

Cycles

Itanium 2 GInst

G I nst C yc

1.0 1.2 1.2 1.2 1.3 1.3

0.24 1.86 13.8 98.5 810 6238

0.20 1.59 12.7 101 812 6493

0.8 0.9 0.9 1.0 1.0 1.0

0.09 0.87 5.54 55.1 398 3529

0.22 2.28 13.9 146 892 9331

2.4 2.6 2.5 2.6 2.2 2.6

Table XII. Execution Time for Indexing a 2D Hilbert Curve at Different Levels Butz & Moore Origin 2000 Cycles GInst

Intel Pentium 4 Cycles GInst

Itanium 2 Cycles GInst

6 8 10 12 14 Curve Level

1.44 26.9 491 8.59K 149K

1.86 36.7 693 12.0K 213K

6.90 4.08 129 78.6 2224 1444 37.7K 24.8K 674K 440K SFCGen

1.79 35.2 635 11.2K 193K

2.93 58.5 1059 18.5K 326K

6 8 10 12 14

0.54 8.73 170 2.99K 55.0K

1.02 18.1 344 6.62K 111K

3.54 73.7 1284 21.7K 377K

1.11 21.8 370 6.14K 104K

2.41 47.0 818 13.9K 239K

Curve Level

2.49 49.1 881 14.8K 261K

SFCGen with an implementation of coordinates-to-index converter hilbert c2i based on Butz’s byte-oriented algorithm. The results show that SFCGen performs from 61% to over a factor of 3 faster than the byte-oriented algorithm across all three platforms. This is mainly because our indexing is much simpler (60% to a factor of 2 fewer graduated instructions) and more efficient with the precomputed table information. The less significant improvement observed on the Intel Pentium 4 and Itanium 2 is caused by the higher cost of bit transpose that both algorithms use to retrieve coordinate information at each level from the original coordinates. Table XIII shows the execution time and the number of graduated instructions for indexing 3D Hilbert curves. Overall, SFCGen outperforms the byte-oriented algorithm by from 42% up to a factor of 3 on the three machines. 6.3 Performance of Higher Dimensional Traversal Our general table-based framework for curve generation and the algorithm that we used to generate movement specification tables enable us to generate different space-filling curves in higher dimensional spaces. Table XIV shows the execution time and the number of graduated instructions of SFCGen for traversing 4D Hilbert, Morton, and Peano curves at different refinement levels. We only show four levels for Peano curves since it has a bigger refinement factor. The results show that the performance of SFCGen scales well with curve level in higher dimensional spaces for all three space-filling curves. On average, ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



145

Table XIII. Execution Time for Indexing 3D Hilbert Curves at Different Levels Butz & Moore Origin 2000 Cycles GInst

Intel Pentium 4 Cycles GInst

Itanium 2 Cycles GInst

4 5 6 7 8 9 Curve Level

0.94 11.2 98.9 937 7.43K 63.8K

1.80 16.1 140 1244 10.6K 92.7K

6.85 3.94 61.9 35.7 503 304 4294 2702 37.6K 23.5K 322K 205K SFCGen

1.63 14.3 125 1083 10.0K 84.6K

2.80 24.4 212 1857 16.9K 145K

4 5 6 7 8 9

0.44 4.04 33.8 338 2.48K 24.4K

0.75 8.65 71.5 628 4.91K 55.0K

3.46 30.7 258 2273 21.2K 186K

1.15 10.6 78.3 722 6.10K 54.9K

2.55 22.6 183 1607 14.1K 121K

Curve Level

2.61 23.9 197 1772 15.5K 137K

Table XIV. Execution Time of SFCGen for Traversing 4D Curves at Different Levels Hilbert curve Origin 2000 Cycles GInst

Pentium 4 Cycles GInst

Itanium 2 Cycles GInst

3 4 5 6 7 Curve Level

0.16 1.61 24.8 396 6330

0.16 2.26 35.8 572 9159

0.17 0.14 2.46 2.18 39.4 34.9 607 558 10.3K 8930 Morton curve

0.07 0.64 17.4 164 4449

0.18 1.82 46.6 465 11.9K

3 4 5 6 7 Curve Level

0.23 2.86 45.7 720 11.7K

0.26 3.77 60.2 963 15.4K

0.25 0.21 4.02 3.30 64.2 52.8 1027 845 16.4K 13.5K Peano curve

0.17 1.92 42.0 491 10.8K

0.47 4.53 121 1158 30.9K

2 3 4 5

0.18 11.4 912 73.8K

0.23 17.9 1447 117K

0.20 13.4 1118 88.1K

0.07 8.12 401 53.3K

0.18 22.6 1196 148K

Curve Level

0.19 15.1 1224 99.1K

the traversal algorithm performs most efficiently on Peano curves based on the ratio of visited points per CPU cycle. Traversing Morton curves is least efficient among the three space-filling curves we studied because of expensive calculations for diagonal moves. 7. CONCLUSIONS Because of their mathematical elegance and their many interesting properties, the subject of space-filling curves has been studied by mathematicians for over a century. Various space-filling curves have been created and many of them have been used in a variety of application areas. For this reason, efficient algorithms ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

146



G. Jin and J. M.-Crummey

for traversal and indexing of space-filling curves are of interest. Motivated by the use of SFC in scientific computation, we proposed a table-driven framework for enumerating and indexing points of a multi-dimensional space along spacefilling curves. Our framework is general enough for use with any self-similar space-filling curve and can be used in any application in which computation can be recursively partitioned along multiple dimensions. We described movement specification tables for four space-filling curves and an algorithm for automatically generating the tables for Hilbert, Morton, and Peano curves. We also described the table generation algorithm for arbitrarily high-dimensional Hilbert, Morton, and Peano curves. We evaluated performance of our framework on three modern microprocessor-based platforms and compared our results with the latest algorithms for traversing and indexing space-filling curves. Our performance results show up to a 63% improvement over the most recent recursive algorithm on 2D curve generation and up to a factor of 132 improvement over two fast implementations of a non-recursive byte-oriented algorithm for curve generation and as much as a factor of three improvement over the byte-oriented algorithm on curve indexing. The results also show our curve generation algorithm scales very well with curve level in high-dimensional spaces. ACKNOWLEDGMENTS

We thank Richard Hanson for his comments on improving the presentation of an early draft of this article and Ronald Goldman for his careful reading and comments on the relationship between movement specification tables and L-systems. We also thank the anonymous referees and the editors for their comments. One referee tested SFCGen and provided four pages of comments, which greatly helped us to improve the software and presentation of this article. REFERENCES ABEL, D. J. AND MARK, D. M. 1990. A comparative analysis of some two-dimensional orderings. Int. J. Geog. Info. Syst. 4, 1, 21–31. ALPERT, C. J. AND KAHNG, A. B. 1994. Multi-way partitioning via spacefilling curves and dynamic programming. In Proceedings of the 31st Annual Conference on Design Automation Conference. 652–657. BARTHOLDI III, J. AND GOLDSMAN, P. 2001. Vertex-labeling algorithms for the Hilbert spacefilling curve. Software–Practice and Experience 31, 395–408. BERTILSSON, M., BRICKELL, E., AND INGEMARSSON, I. 1989. Cryptanalysis of video encryption based on space-filling curves. In Advances in Cryptology—EUROCRYPT’89. 403–411. BIALLY, T. 1969. Space-filling curves: Their generation and their application to bandwidth reduction. IEEE Trans. Info. Theory IT-15, 6 (Nov.), 658–664. ¨ , C., BERCHTOLD, S., AND KEIM, D. A. 2001. Seaching in high-dimensional spaces: Index strucBOHM tures for improving the performance of multimedia databases. ACM Comput. Surv. 33, 322–373. BREINHOLT, G. AND SCHIERZ, C. 1998. Algorithm 781: Generating Hilbert’s space-filling curve by recursion. ACM Trans. Math. Soft. 24, 2, 184–189. BROWNE, S., DONGARRA, J., GARNER, N., LONDON, K., AND MUCCI, P. 2000. A scalable cross-platform infrastructure for application performance tuning using hardware counters. In Proceedings of 2000 ACM/IEEE Conference on Supercomputing. Dallas, TX. BUTZ, A. R. 1968. Space filling curves and mathematical programming. Information and Control 12, 314–330. BUTZ, A. R. 1969. Convergence with Hilbert’s space filling curve. J. Comput. Syst. Sci. 3, 2, 128– 146. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

SFCGen: A Framework for Efficient Generation of Space-Filling Curves



147

BUTZ, A. R. 1971. Alternative algorithm for Hilbert’s space-filling curve. IEEE Trans. Comput. C-20, 4 (Apr.), 424–426. CHATTERJEE, S., LEBECK, A. R., PATNALA, P. K., AND THOTTETHODI, M. 1999. Recursive array layouts and fast parallel matrix multiplication. In ACM Symposium on Parallel Algorithms and Architectures. 222–231. COLE, A. J. 1981. A note on space filling curves. Software–Practice and Experience 13, 12 (Dec.), 1181–1189. ¨ , A., AND DALLA, L. 1999. Fractal interpolation techniques for DRAKOPOULOS, V., TZIOVARAS, A., BOHM the generation of space-filling curves. In Hellenic European Research on Computer Mathematics and Its Applications, E. A. Lipitakis, Ed. LEA, 843–850. FRENS, J. AND WISE, D. 1997. Auto-blocking matrix multiplication or tracking BLAS3 performance from source code. In Proceedings of the Sixth ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming. Las Vegas, NV, 206–216. GOLDSCHLAGER, L. M. 1981. Short algorithms for space-filling curves. Software–Practice and Experience 11, 1 (Jan.), 99–100. ¨ ¨ die stetige abbildung einer linie auf ein flachenstuck. Math. Ann. 38, HILBERT, D. 1891. Uber 459–460. HU, Y. C., COX, A., AND ZWAENEPOEL, W. 2000. Improving fine-grained irregular shared-memory benchmarks by data reordering. In Proceedings of the 2000 ACM/IEEE Conference on Supercomputing. Dallas, TX. JAGADISH, H. V. 1990. Linear clustering of objects with multiple attributes. In Proceedings of the 1990 ACM SIGMOD Conference. 332–342. JIN, G., MELLOR-CRUMMEY, J., AND FOWLER, R. 2001. Increasing temporal locality with skewing and recursive blocking. In Proceedings of the 2001 ACM/IEEE Conference on Supercomputing. Denver, CO. KUMAR, A. 1994. A study of spatial clustering techniques. In Proceedings of the 5th International Conference on Database and Expert System Applications. 57–71. KUO, S., WINSLETT, M., CHO, Y., LEE, J., AND CHEN, Y. 1999. Efficient input and output for scientific simulations. In Proceedings of the 6th Workshop on I/O in Parallel and Distributed Systems. 33–44. LAMARQUE, C. H. AND ROBERT, F. 1996. Image analysis using space-filling curves and 1D wavelet bases. Pattern Recognition 29, 8, 1309–1322. LAWDER, J. K. AND KING, P. J. H. 2001. Querying multi-dimensional data indexed using the Hilbert space-filling curve. ACM SIGMOD Record 30, 19–24. LINDENMAYER, A. 1968. Mathematical models for cellular interaction in development, parts I and II. J. Theore. Biol. 18, 280–315. MATIAS, Y. AND SHAMIR, A. 1987. A video scrambling technique based on space filling curve. In Advances in Cryptology—CRYPTO’89. 398–416. MELLOR-CRUMMEY, J., WHALLEY, D., AND KENNEDY, K. 1999. Improving memory hierarchy performance for irregular applications. In Proceedings of the 13th ACM International Conference on Supercomputing. Rhodes, Greece, 425–433. MOGHADDAM, B., HINTZ, K. J., AND STEWART, C. V. 1991. Space-filling curves for image compression. In Proceedings of the SPIE. 414–421. MOKBEL, M. F. AND AREF, W. G. 2001. Multimedia information processing: Irregularity in multi-dimensional space-filling curves with applications in multimedia database. In Proceedings of the Tenth International Conference on Information and Knowledge Management. 512– 519. MOKBEL, M. F., AREF, W. G., AND KAMEL, I. 2003. Analysis of multi-dimensional space-filling curves. GeoInformatica 7, 179–209. MOORE, D. 2000. http://www.caam.rice.edu/~dougm/twiddle/Hilbert. MORTON, G. M. 1966. A computer oriented geodetic data base and a new technique in file sequencing. Technique Report, IBM, Canada. MUSGRAVE, K. 1991. A Peano curve generation algorithm. In Graphics Gems II, J. Arvo, Ed. Academic Press, San Diego, CA, 25. OHNO, Y. AND OHYAMA, K. 1991. A catalog of symmetric self-similar space-filling curves. J. Recreational Math. 23, 161–173. ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.

148



G. Jin and J. M.-Crummey

PEANO, G. 1890. Sur une courbe qui remplit toute une aire plane. Mathematishe Annalen 36, 157–160. PLATZMAN, L. K. AND J. J. BARTHOLDI, III. 1989. Spacefilling curves and the planar travelling salesman problem. J. ACM 36, 719–737. PRUSINKIEWICZ, P., LINDENMAYER, A., AND FRACCHIA, F. D. 1991. Synthesis of space-filling curves on the square grid. In Fractals in the Fundamental and Applied Sciences, H. O. Peitgen, J. M. Henriques, and L. F. Penedo, Eds. Elsevier Science Publisher BV, Amsterdam, The Netherlands, 341–366. SAGAN, H. 1992. On the geometrization of the Peano curve and the arithmetization of the Hilbert curve. Int. J. Math. Educ. Sci. Tech. 23, 403–411. SAGAN, H. 1993. A three-dimensional Hilbert curve. Int. J. Math. Educ. Sci. Tech. 24, 541–545. SAGAN, H. 1994. Space-Filling Curves. Springer-Verlag, New York, NY. ´ , W. 1912. Sur une nouvelle courbe countinue qui remplit toute une aire plane. Bull. SIERPINSKI Acad. Sci. de Cracovie (Sci. math. et nat., S´erie A), 462–478. STEVENS, R. J., LEHAR, A. F., AND PRESTON, F. H. 1983. Manipulation and presentation of multidimensional image data using the Peano scan. IEEE Trans. Pattern Anal. Mach. Intell. 5, 520–526. VELHO, L. AND GOMES, J. M. 1991. Digital halftoning with space filling curves. In Proceedings of the 18th Annual Conference on Computer Graphics and Interactive Techniques. 81–90. WITTEN, I. H. AND WYVILL, B. 1983. On the generation and use of space-filling curves. Software– Practice and Experience 13, 519–525. ZHANG, Y. AND WEBBER, R. E. 1993. Space diffusion: An improved parallel halftoning technique using space-filling curves. In Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques. 305–312. Received November 2003; revised July 2004 and October 2004; accepted October 2004

ACM Transactions on Mathematical Software, Vol. 31, No. 1, March 2005.