Cyclic Interlaced Quadtree Algorithms for Quincunx Multiresolution

JOURNAL OF ALGORITHMS ARTICLE NO.

27, 97]128 Ž1998.

AL980925

Cyclic Interlaced Quadtree Algorithms for Quincunx MultiresolutionU D. J. Hebert Department of Mathematics, Uni¨ ersity of Pittsburgh, Pittsburgh, Pennsyl¨ ania 15260 Received January 9, 1997; revised December 30, 1997

Recent advances in wavelet theory and in finite element computations draw attention to a well-known, simple, and computationally efficient triangulation method. We take a new look at this triangulation, which is obtained by repeated symmetric bisection, starting with a half square. The cells form the leaves of a binary tree and the nodes of a directed graph consisting of a single simple cycle. Computational speed is facilitated by the binary and quad-digit expression of triangle vertices, which reduce all vertex calculations to simple integer and logical operations. The leaf cycle interlaces a pair of quadtrees whose orientations differ by pr4. Detailed analysis leads to algorithms which exploit the structure and computational efficiencies in calculations such as pyramid algorithms for image processing with non-separable wavelets. Q 1998 Academic Press

1. INTRODUCTION Consider a simple triangulation of the unit square obtained by repeated bisection of triangular cells starting with the triangle which divides the square along the diagonal containing the origin. Such triangulations appeared in standard tesselations of the plane Žcf. w1x., and have been used for efficient computations of approximations of functions and solutions of partial differential equations by the finite element method. For example, it has long been recognized that the grids may be generated by reflections, leading to directed paths through the cells which facilitate program organization and computing speed Žcf. w2x.. Such paths appeared earlier as approximations to a space filling curve now known as the Sierpinski curve Žcf. w3x.. For finite element calculations the grids offer advantages such as the data structures of binary trees, easy neighbor finding, and avoidance of * ICMA-96-200, version 2.0. E-mail: [email protected]. 97 0196-6774r98 $25.00 Copyright Q 1998 by Academic Press All rights of reproduction in any form reserved.

98

D. J. HEBERT

global communications Žcf. w4, 5x.. Maubach w6x has shown that these triangulations form a special case of a general n-dimensional structure. Further, as noted in w7x, a natural symbolic labeling of the tetrahedral cells allows refinement, pruning, vertices, and neighbor-finding to be computed with simple integer and logical operations. This labeling is closely related to the labeling of quadtree elements by location codes which has been well developed as an efficient computational tool for applications to computer graphics and image processing Žsee w8x for detailed studies and an extensive bibliography.. Recent developments and applications may also be found in w9x, w10x, and w11x. Quadtrees are also especially suited for some types of parallel computing Žcf. w12x and its references .. The same triangulations have appeared as a natural setting for a special class of box splines including the quadratic Zwart spline w13x generated by four direction vectors and supported by a triangulated octahedron Žfor a full study of box splines see the book by de Boor et al. w14x and its references .. Chiu, Jetter, and Stockler, have shown that the four-direc¨ tional box splines may be used to generate a multiresolution analysis Žcf. w15x. and a wavelet which generates a frame for L2 . To generate a multiresolution analysis, the wavelet function is scaled by powers of the matrix Ms

1 1

1 y1

and translated to cartesian integer lattice points along with the M-transforms of these points which comprise the ‘‘quincunx lattice’’ Žcf. w16x.. Very recently, Ron and Shen w17x have shown how to construct directional box spline wavelets based on the same scaling. The wavelets of Ron and Shen have small support and generate a tight frame with frame bound 1. Practical wavelet calculations, in terms of discrete wavelet filters for image processing applications, are usually restricted to the separable tensor products of standard one-dimensional wavelets. However, for application to edge-detection, orientation analysis, and motion analysis Žcf. w18x., the separable wavelet filters do not perform as well as specialized non-separable directional wavelet filters such as those based on the quincunx lattice Žcf. w19x.. Computations with non-separable wavelets are generally much more complex than those of the separable case. Even in the highly structured quincunx lattice setting the organization of efficient computations has not been apparent. The purpose of the present article is to show that four-directional calculations for wavelet transforms and other applications can be obtained by a further detailed study of binary triangulation trees. The result is an algorithm which serially traverses a uniform triangulation, interlacing the

INTERLACED QUADTREE ALGORITHMS

99

leaves of two quadtrees, using memory stacks and only a small number of simple logical and small-table operations to perform calculations similar to those of the most efficient wavelet image decomposition, the pyramid algorithm of the ‘‘fast wavelet transform.’’ In contrast to the usual line-byline scan calculation of the pyramid algorithm in which the row and column calculations can be performed separately, a triangular scanning with temporary stack storage makes possible the serial calculation of two-dimensional non-separable filter masks such as the four-directional box splines w20x. The serial processing of a single scan path also suggests applications to image analysis and compression with implementations via DSP chips, array processors, or specialized hardware. In addition to wavelet analysis, there are natural applications of interlaced quadtree algorithms in a wide range of computational studies. For example, typical linear calculations of solutions of algebraic and differential equations on rectangular and triangular grids are expressed in terms of matrices whose rows and columns are indexed by cells of the grid or by their vertices and special points. Discrete space-filling curves such as the Sierpinski scan paths provide a convenient ordering for rows and columns of matrices. In this ordering the matrices can be segmented into blocks for performing local calculations in parallel. In many cases of multiscale analysis there are large collections of contiguous cells for which no fine-scale calculations are necessary. In these cases simple binary tree pruning as in w21x provides a natural method of load balancing for parallel calculations. Section 2 introduces the Cartesian and quincunx quadtree lattice origin Žlo. codes and reviews the following properties noted in w21x: The cells of a binary triangulation form the leaves of a binary tree and the nodes of a directed graph consisting of a single simple cycle called the leaf cycle. The triangles are uniquely located by a label consisting of an lo code along with a triangle number and the parity of the level. From this label we may determine all the properties of the triangle by small dyadic integer calculations. The tree position of a triangular cell is determined by its level and a level leaf cycle index. Section 3 presents some examples from wavelet-based image processing: a simple triangular Haar wavelet transform, a compass edge-detecting square Haar wavelet transform, a discrete edge-detecting box spline wavelet related to the wavelets of Ron and Shen, and a four-directional box spline pyramid algorithm. The remainder of the paper is devoted to a detailed analysis of the computational structures in the binary triangulation tree. Section 4 gives an algorithm for incrementing lo codes for redundant paths through the lattice origins called lo paths. The removal of redundancy by examining turn flags gives an algorithm for traversing a leaf cycle in the triangulation

100

D. J. HEBERT

tree. Finally the triangle vertices and midpoints are coded in terms of lo codes eliminating the need for numerical vertex calculations. In this odd-level triangles split a cartesian quadtree square block into four parts sharing a vertex at the center of the block while even-level triangles split a diagonally oriented square into four parts sharing a vertex which is the center of a quincunx lattice quadtree block. Thus two related quadtrees are interlaced by the leaf cycle paths of the odd- and even-level uniform triangular grids. A walk through the leaf cycle at a given level forms a path, called a Sierpinski path, which fills the square, a path which is similar in nature to the well-known Peano]Hilbert scan path recently used for improvement of image compression Žcf. w22, 23x.. Section 5 is devoted to the proof that the stack-based memory operations are correct. Each triangle forms half of a quadtree block, and the path through the triangles visit the two halves in such a manner that data from one half may be pushed to a stack to be popped on the second visit for completion of a calculation with data from the whole block. Section 6 classifies the points in an lo path by a sequence of four neighboring left and right turns called a turn list. Section 7 establishes the patterns necessary to step through adjacent leaf cycle levels for the interlaced quadtree algorithm of Section 8 which uses only logical decisions, digit increments, and small table references to calculate triangle vertices as it steps through adjacent Cartesian and quincunx quadtrees. The final Section 9 shows how the interlaced quadtree algorithms are implemented in the examples of Section 2.

2. BINARY TRIANGULATION TREES AND LATTICE ORIGIN QUADTREES 2.1. A Binary Triangulation Tree Consider a binary tree constructed as follows: The root is the unit square I 2 s w0, 1x2 . The two children of the root, called level zero triangles, split the unit square and share the diagonal ŽŽ0, 0., Ž1, 1... To obtain the children of the level zero triangles, we divide each of these into two by bisecting the diagonal. We continue at each level l to bisect each triangle at the longest edge to obtain two subtriangles which are the children at level l q 1 in a binary tree. The triangles of levels 0]3 are shown in Fig. 1. Note that for l ) 0 each triangle of level l q 2 is a ‘‘half-scale’’ version of some triangle of level l. Each triangle obtained by this procedure is one of the forms shown in Fig. 1 Žb, c. located in a square with edge of length 1r2 m , where l s 2 m or l s 2 m q 1. At level l all triangles have the same similarity class and are said to be of type l mod 2. Thus, at each level l the

INTERLACED QUADTREE ALGORITHMS

101

FIG. 1. Triangles of levels 0]3.

triangles form a uniform grid which subdivides the unit square. Indeed, the leaves of any finite binary subtree form a triangulation of I 2 . With the additional edge compatibility condition that no triangle shares an edge with more than one triangle, such a triangulation will be called a binary triangulation. The special vertex at the center of the square with edge length 1r2 m , yky1 called the Žcartesian . lattice origin, has the form Ý m s k , where s k ks 0 2 is an ordered pair of signs Ži.e., pairs of elements of the set  1, y14.. A binary triangulation cell is thus identified by a list of digits Ž l, t, j1 , . . . , jm . which we shall call its label, where l is the level Žmod 2., t is the triangle number as in Fig. 1Žb, c., and ji is a 4-ary digit representing the pair si according to the following assignment of sign pair to Boolean pair to 4-ary digit: Žy1, y1. ª Ž0, 0. ª 0, Žy1, 1. ª Ž0, 1. ª 1, Ž1, y1. ª Ž1, 0. ª 2, Ž1, 1. ª Ž1, 1. ª 3. The number m is called the depth of the lattice origin. 2.2. Le¨ el Leaf Cycles The set of triangles of level l form the leaves of a finite full binary tree whose root is the unit square. The triangles of level l also form the nodes of an undirected graph whose links connect triangles which share a common edge. At each level there is a natural ordering of the two children of a given triangle given geometrically by counterclockwise motion around the cartesian lattice origin. Making use of this ordering of children we may establish a natural linear ordering of the leaves of the full simplex tree of any depth. PROPOSITION 1. For each l ) 0 the triangles of le¨ el l form the nodes of a directed graph consisting of a single simple cycle in which linked triangles share an edge. The element s of this le¨ el leaf cycle are indexed in integer order by the positi¨ e integers whose l binary digits locate the triangles as elements of a binary tree. Proof. Figure 1Žb, c. shows an ordering and numbering at levels 1 and 2. If a level leaf cycle is given at level l with triangles numbered 0, . . . , 2 lq1 y 1 then for triangle i with binary digits bi : i s Ý kjs0 bj 2 j, we assign

102

D. J. HEBERT

to child 0 and child 1 Žordered counterclockwise about the lattice origin. the numbers 2 i and 2 i q 1 to establish an ordering at level l q 1. The set of triangles at a given level form the nodes of a graph whose links represent shared edges. The number of triangles gives a direction to a set of links in the graph to form a directed subgraph consisting of a single closed path. Thus, the linear ordering of the leaves on each level gives a direction to a set of links in the leaf graph so that the leaves form a directed subgraph consisting of a single closed path. We shall call such a path a level leaf cycle. The level leaf cycle at level 5 with cycle index numbering is illustrated in Fig. 2. The curve which connects the centers of the triangles will be called a Sierpinski path w3x. 2.3. Cartesian and Quincunx Lattice Origin Quadtrees If Z = Z is the Cartesian lattice of integer points, then the set of quincunx lattice points is  Mx: x g Z = Z 4 , where Ms

1 1

1 y1

is the quincunx matrix Žcf. w15x.. We consider now certain lattice points in the unit square which are dyadic scalings of the cartesian integer lattice points and the quincunx lattice points. We define a cartesian lattice origin Žclo. point of depth m to be a point in w0, 1x 2 of the form m

Ý ks0

1 2

kq1

sk ,

FIG. 2. Cycle index of triangles at level 5.

INTERLACED QUADTREE ALGORITHMS

103

where s k is an ordered pair of signs, i.e., of elements of  y1, 14 , and s 0 s Ž1, 1.. A quincunx lattice origin Žqlo. point of depth m is a point in w0, 1x2 of the form 1 1 , q 2 2

ž /

m

Ý ks1

1 2k

fk ,

where f k is an ordered pair of the form Ž1r2. Ms kT , s k is a pair of signs, and M is the quincunx matrix. The clo points form nodes of a quadtree whose root is Ž1r2, 1r2. such that the children of a clo point of depth m are the four elements of depth m q 1 whose digits numbered 0 through m coincide with those of x. Similarly the qlo points form the nodes of a quadtree with root Ž1r2, 1r2. Žsome of whose nodes lie outside the unit square, hence outside the set of qlo points.. Figure 3 shows the quadtree nodes as dots whose sizes decrease as level numbers increase. For simplicity of terminology let us agree now that a qlo point is one which is within the unit square. We represent the clo points and the qlo points by lists of 4-ary digits as follows: First define s Ž0. s Žy1, y1., s Ž1. s Žy1, 1., s Ž2. s Ž1, y1., s Ž3. s Ž1, 1.. Next define f Ž j . s Ž1r2. Ms Ž j .T for j s 0, 1, 2, 3, so that f Ž0. s Žy1, 0., f Ž1. s Ž0, y1., f Ž2. s Ž0, 1., and f Ž3. s Ž1, 0.. Then for a list j s Ž j0 , . . . , jm . of 4-ary digits we assign the clo point m

x Ž j. s

Ý ks0

1 2

kq1

s Ž jk . .

The clo point x is thus represented by a list j Ž x . s Ž j0 Ž x . , j1 Ž x . , . . . , jm Ž x . . , where j0 Ž x . s 3. For calculations it is convenient to also represent x by the 2 = m binary matrix JŽ x . whose kth column consists of the two binary

FIG. 3. qlo and clo quadtrees.

104

D. J. HEBERT

digits of jk . Similarly, we assign to the list q s Ž q0 , . . . , qm . the qlo point y Ž q. s

1 1 , q 2 2

ž /

m

Ý ks1

1 2k

f Ž qk .

and write q Ž y . s Ž q0 Ž y . , . . . , qm Ž y . . , where q0 Ž y . s 2 as a convenient way to distinguish a qlo code from a clo code by its 0th digit. The lists jŽ x . and qŽ y . are called the clo code of x and the qlo code of y and are equivalent to the quadtree locational codes described, for example, in w8x.

3. EXAMPLES FROM WAVELET-BASED IMAGE PROCESSING Let us view a two-dimensional gray scale image as a matrix of size 2 m = 2 m , whose entries are pixel intensities. We imagine that the pixel intensities are located at the centers of cartesian quadtree blocks, the clo points of depth m. Now let us associate with each clo point x the triangle of level 2 m y 1 which contains x on an edge and which follows x in a counterclockwise circling of the parent of x. Figure 4 shows an arrow from each clo point of depth 3 into the triangle of level 5 which it represents. In this way we associate pixel intensities with triangles. The ordering of the pixels according to the leaf cycle ordering at level 2 m y 1 is called a Sierpinski scan path. In the examples below, we use a Sierpinski scan path to reduce the computation of two-dimensional non-separable wavelet transforms of images to one-dimensional processing algorithms.

FIG. 4. Arrow from the clo point into the triangle it represents.

INTERLACED QUADTREE ALGORITHMS

105

3.1. Triangular Haar Wa¨ elets The Žone-dimensional. fast wavelet transform Žcf. w24x. is a fast and simple method for image encoding and decoding. The transform is determined by two ‘‘filter masks’’ in the form of sequences of finite support  g i 4 , the wavelet Žor high-pass. filter mask, and  h i 4 , the smoothing function Žor low-pass. filter mask. An input ‘‘signal,’’ a sequence real numbers  a im 4 of length 2 m is transformed as follows: For l s m, m y 1, . . . , 0 the sequence  a ly1 4 is computed as the convolution of the signal sequence  a il 4 at level l i with a filter sequence: a ly1 s i

Ý h j a2l iyj , d ily1 s Ý g j a2l iyj .

Ž 1. Ž 2.

The output is the structured sequence

 a00 , d00 , d01 , d11 , d02 , . . . , d32 , . . . , d0my 1 , . . . , d nmy 1 4 , where n s 2 my 1. The simplest example is the ‘‘Haar wavelet filter’’ for which g 0 s y1r2, g 1 s 1r2, and h 0 s h1 s 1r2, g i s h i s 0 otherwise. The transform is given by a ly1 s i

Ž a2l iq1 q a2l i . , d ily1 s 21 Ž a2l iq1 y a2l i . . 1 2

Ž 3. Ž 4.

The detail coefficients d il along with the final result a0i s const are stored for each l as a transform Žan encoding. of the signal. The set of coefficients is the same size as the index set of the signal values. In the case of Haar wavelet filters the inverse transform, which reconstructs the signal from the stored coefficients, is simply implemented by the formulae a2l iq1 s a ily1 q d ily1 ,

Ž 5.

a2l i s a ly1 y d ily1 . i

Ž 6.

A simple extension of the one-dimensional Haar wavelet pyramid filter for image processing is given in w21x. The idea is to apply the one-dimensional fast wavelet transform for the Haar wavelet directly to the one-dimensional set of intensities of a Sierpinski scan of the pixels, viewing each pixel as a small triangle. This is equivalent to performing a two-dimensional wavelet transform in which the wavelets are the step functions supported by the triangles of a binary triangulation with values y1 on the first child and q1 on the second child in the binary tree of triangles.

106

D. J. HEBERT

3.2. Wa¨ elet Edge and Feature Detection Next we consider a standard method of edge detection in two-dimensional images based on the application of discrete directional or gradient filters Žsee w25x.. A filter mask is defined in terms of a small matrix D s  d i j 4 which is skew-symmetric about a central horizontal, vertical, or diagonal line. For a pixel image F s  f i j 4 the magnitude of the sum Ý i j f pqi, qqj d i j indicates the prominence of an ‘‘edge,’’ i.e., of the variation of image intensity across the direction of the line of skew-symmetry. Other local features such as ‘‘texture’’ may be detected by other symmetries and asymmetries. We are interested in discrete wavelet filters, i.e., filters which may be obtained by scalings of linear transformations of a simple mask. For computational work we restrict attention to wavelet filters of finite support. 3.2.1. Example: Haar Wa¨ elets Supported by Cartesian and Quincunx Squares Let a2l i be a number Žintensity . associated with triangle t 2 i of cycle index 2 i in a level leaf cycle path through the triangles of level l. The triangles t 2 i and t 2 iq1 share a lattice origin l as a vertex and, for some index 2 iU , the triangles t 2 iU and t 2 iU q1 also share the vertex l. We associate three Haar wavelet coefficients with the quadtree block centered at l: wll , 1 s

Ž a2l i y a2l iq1 q a2l i y a2l i q1 . s Ž yd ily1 y d ily1 . , wll , 2 s 12 Ž a2l i q a2l iq1 y a2l i y a2l i q1 . s Ž a ly1 y a ly1 ., i i wll , 3 s 21 Ž a2l i y a2l iq1 y a2l i q a2l i q1 . s Ž yd ily1 q d ily1 . , 1 2

U

U

U

U

U

U

U

U

U

where a ly1 s 12 Ž a2l i q a2l iq1 ., d ily1 s 12 Ž a2l iq1 y a2l i ., etc. The wavelet coi efficients wll, 2 and wll, 3 represent the application of a two-valued filter mask supported by a quadtree block detecting an edge which Žroughly. splits the block along a diagonal line. If slly1 s

Ž a2l i q a2l iq1 q a2l i q a2l i q1 . s aily1 q aily1 , w T s Ž slly2 , wll , 1 , wll , 2 , wll , 3 . , and aT s Ž a2l i , a2l iq1 , a2l i , a2l i q1 . 1 2

U

U

U

U

then w s Ba,

U

INTERLACED QUADTREE ALGORITHMS

107

where Bs

1 2



1 1 1 1

1 y1 1 y1

1 1 y1 y1

1 y1 y1 1

0

is the symmetric, orthogonal Hadamard matrix so that B s B T s By1 is non-singular. The invertibility of the matrix B leads to the invertibility of a wavelet transform in which the stored coefficients are the wavelet coefficients wlj, 1 and wlj, 3 for j F l. As in the triangular Haar wavelet transform we traverse the leaf cycles at each level. As we reach a2j i in the traversal we compute d ijy1 and a ijy1 from a2j i and a2j iq1 , saying d ijy1, a2j i , a2j iq1 , and a ijy1 in temporary storage. When we reach the square again at index iU we use U U d ijy1 along with d ijy1 to compute Ž wlj, 1 , wlj, 3 .. We also compute a ijy1 to put in temporary storage. This procedure is continued for each level j until level 0 is reached, and we must save a00 and a10 . The triangle intensities a ij may be reconstructed from the coefficients a00 , a10 and wlj, 1 , wlj, 3 for j G 0 by a similar procedure. The special structure of the binary triangulation tree makes it possible to replace the temporary storage mentioned above with more efficient stack storage, as we shall see from the detailed analysis of later sections. For edge detection applications of this and other examples, see w26x. For other recent developments concerning Haar wavelets based on triangulations see w27x and w28x. 3.2.2. Example: Box Spline Wa¨ elets Experiments show that the Haar wavelet edge detection mask is far from ideal. It is known that the ideal edge detecting filters have the shape of directional derivatives of the two-dimensional Gaussian density function. The recent work of Ron and Shen w17x shows the construction of a computationally simple approximation to the compass directional derivatives of the Gaussian in the form of box spline wavelets of small support. These box spline wavelets may be computed by iterated subdivision algorithms applied to skew symmetric arrangements of q and y signs on the filter mask 1 16



0 1 1 0

1 2 2 1

1 2 2 1

0 1 . 1 0

0

The above mask without signs is the discrete four-directional box spline. A subdivision algorithm applied to this mask leads to the well-known

108

D. J. HEBERT

Zwart]Powell elements studied, for example, in w14x and w29x. The interlaced quadtree algorithms to be developed in later sections provide a structure for simply computing, for each level l, a set of coefficients sil which represent the coefficient of a discrete box spline centered at the kth lo point at level l. The discrete box spline masks are naturally interpreted as discrete low-pass filters for a possible pyramid algorithm. Since the integer translates of Zwart splines are linearly dependent Žcf. w14x. we cannot expect to construct a fast wavelet transform based on an orthonormal wavelet basis, nor even an efficient biorthogonal wavelet transform Žcf. w24x.. However, since the half-integer translates form a frame Žcf. w29x., we may hope to construct some sort of redundant set of wavelets for the representation of images. To this end let a ly1 s Ž1r2.Ž s2l iq1 q s2l i . and d ily1 s i l l Ž1r2.Ž s2 iq1 y s2 i .. Consider the Haar-like coefficient a ily1 as a prediction of the box spline coefficient sily1, and let ¨ ily1 s sily1 y a ly1 . A pyramid i algorithm is obtained as follows: Given intensities sil at level l, we calculate sily1 by a quadtree interlace algorithm and save the coefficients d ily1 and ¨ ily1 as wavelet coefficients. To reconstruct we start with sily1, d ily1, and ¨ ily1 and compute a ly1 s sily1 y ¨ ily1. Then we have s2l iq1 s i ly1 ly1 l ly1 ly1 a i q d i and s2 i s a i y d i . In terms of the coefficients a ly1 and d ily1 we can define wavelets wll, j as i in Section 3.2.1. These wavelets can be interpreted as differences of discrete box splines with geometric meaning as detectors of shape and texture. They also turn out to be the masks which generate the edge detecting box splines of Ron and Shen w17x.

4. LATTICE ORIGIN PATHS AND TRIANGULAR SCAN PATHS 4.1. Cartesian Lattice Origin Paths We now consider certain paths through lattice origins defined in terms of the locational codes. We first define a redundant path through the cartesian lattice origins in which some points occur twice. Then, by removing the redundancies we define simple closed paths which are similar to the Sierpinski paths. Finally, we associate each lattice origin of depth m with a binary triangle of level 2 m y 1 and find that the nonredundant lattice origin path lists the triangles in the leaf cycle ordering. The triangles of level 2 m are then associated in a similar manner with the quincunx lattice origins. A 4-ary digit j is naturally associated with the quadrant containing s Ž j .. We refer to the cyclic ordering 0, 2, 3, 1 of the digits as the Žcounterclockwise. quadrant ordering and abbreviate the successor and predecessor

INTERLACED QUADTREE ALGORITHMS

109

functions by writing 0qs 2, 2qs 3, 3qs 1, 1qs 0, and 2ys 0, 0ys 1, 1ys 3, 3ys 2. Each quadrant has its complement or negative, and we write 0 c s 3, 3 c s 0, 1c s 2, 2 c s 1. The paths of interest move counterclockwise around the quadtree parent until reaching a turn into another block or an edge of the unit square. Here are the essential definitions A turning point is a list of the form Ž j0 , . . . , jiy1 , ji , . . . , jm ., where ji s jiq1 s ??? s jm and jic s jiy1. The following is a definition of a redundant path through the clo points Žcalled a clo path. which was introduced in w21x DEFINITION 1 Žclo path.. Ž1. The clo path through the cartesian lattice origins of depth m starts at Ž3, 0, . . . , 0.. Ž2. If j s Ž j0 , . . . , jm . is in the path then the next element jq in the path is determined as follows: if j is a turning point of the form Ž j0 , . . . , jiy1 , ji , . . . , jm ., where ji s ??? s jm and the precedessor of j is not a turning point, q q. then jqs Ž j0 , . . . , jq iy1 , j i , . . . , j m , . else jqs Ž j0 , . . . , jmy1 , jq m . An example of a clo path is given in Fig. 5. The clo path starts at the lower left corner and follows the arrows turning left at each turning point and right at other points. We may note that each turning point is visited twice. The number of steps in a clo path may be calculated by an induction argument.

FIG. 5. clo path at depth 4.

110

D. J. HEBERT

L cy1 PROPOSITION 2. The number of clo points at le¨ el m is 2 2 m . If  x km 4ks0 is the list of clo points in the clo path at le¨ el m, then

L c s 43 Ž 2 2 m y 1 . . 4.2. Quincunx Lattice Origin Paths To keep the paths through the quincunx lattice origins within the unit square we need to identify points adjacent to corners: A qlo point y with qlo code q s Ž q0 , . . . , qn . is a corner point if q is not a turning point, if q 1 s qq 2 and for each k G 2, q 2 s qk . DEFINITION 2 Žqlo path.. Ž1. The qlo path through the quincunx lattice origins of depth m starts at Ž2, 1, 0, . . . , 0.. Ž2. If q s Ž q0 , . . . , qm . is in the path then the next element qq in the path is determined as follows: if q is a turning point of the form Ž q0 , . . . , qiy1 , qi , . . . , qm ., where qi s ??? s qm , and the precedessor of q is not a turning point, q q. then qqs Ž q0 , . . . , qq iq1 , qi , . . . , qm , . else qqs Ž q0 , . . . , qmy1 , qq m . q q q. Ž3. If q is a corner point then qq s Ž qy 0 , q1 , q 2 , . . . , q m . An example qlo path is given in Fig. 6. The qlo path also starts at the lower left corner and follows arrows, but turns left at each turning point and right at other points. The number of steps in a qlo path may be calculated by a recursion formula. PROPOSITION 3. The number of distinct qlo points at le¨ el m is L q Ž m.y1 2 m Ž2 my 1 q 1.. If  y k 4ks0 is the list of qlo points at le¨ el m then the number L q Ž m. of steps in the qlo path is gi¨ en by a recursion formula, Lq Ž m. s

4,

if m s 1,

4 L q Ž m y 1. y 1. ,

for m ) 1,

½Ž

and by the explicit formula L q Ž m . s 43 Ž 2 2 my1 q 1 . . 4.3. Non-redundant Listing of Lattice Origins If the unit square w0, 1x 2 is subdivided into 2 2 m squares of side 1r2 m by the lines x s ir2 m and y s jr2 m , for i s 0, . . . , 2 m and j s 0, . . . , 2 m , then there is a cartesian lattice origin of depth m located at the center of

INTERLACED QUADTREE ALGORITHMS

111

FIG. 6. qlo path at depth 4.

each subdividing square. Indeed, the clo points are the points of the form ŽŽ2 i y 1.r2 mq 1 , Ž2 j y 1.r2 mq 1 . for i s 1, . . . , 2 m and j s 1, . . . , 2 m . We may make a non-redundant listing of these clo points of depth m by following the clo path and by eliminating the first turn point of each consecutive pair. N cy1 ALGORITHM 1 ŽNon-redundant Listing of clo Points.. Let  x k 4ks0 be a N y1 c listing of the points in the clo path of depth m, and let t k 4ks0 be the corresponding list of turn flags Žt k g  0, 14 , t k s 1 iff x k is a turn point.. Define y 0 s x 0 , p s 1. For 1 F k F Nc y 1 if t k s 0 or t ky1 s 0 then y p s x k , p ª p q 1 endif 2 y1 Then  y p 4ps0 is a non-redundant listing. The qlo points may be listed non-redundantly in almost the same manner. 2m

N cy1 ALGORITHM 2 ŽNon-redundant listing of qlo points.. Let  x k 4ks0 be a N y1 c listing of the points in the qlo path of depth m, and let t k 4ks0 be the corresponding list of turn flags, where corner points are included as turn points. Define y 0 s x 0 , p s 1. For 1 F k F Nc y 1 if t k s 0 or t ky1 s 0 then y p s x k , p ª p q 1 endif

112

D. J. HEBERT

4.4. Listing of Triangles of Le¨ els 2 m y 1 and 2 m The triangles of levels 2 m q 1 and 2 m q 2 all share, as a common vertex, a clo point of depth m. For each clo point x of depth m there are four quadtree children x j s x q Ž1r2 mq 2 . s Ž j ., where s Ž j . is the sign pair corresponding to the 4-ary digit j. We may associate with each 4-ary digit j the level 2 m q 1 triangle which contains x j and which follows x j in a counterclockwise circling of x. A non-redundant listing of the clo points of depth m q 1 then gives a simultaneous listing of the triangles of level 2 m q 1. Note that the triangles of level 2Ž m q 1. may be listed by substituting, for each level 2 m q 1 triangle, the two binary simplex tree children in counterclockwise ordering, as the level m q 1 clo path is traversed non-redundantly. The following proposition, which may be proved by induction, says that the listings of triangles above are precisely the listings of the level leaf cycles. PROPOSITION 4. The non-redundant listing of the clo points of depth m q 1 uniquely determines a listing of the triangles of le¨ el 2 m q 1 which coincides with the listing of the le¨ el leaf cycle at le¨ el 2 m q 1 in the binary simplex tree. A similar proposition makes a connection between level 2 m triangles and the non-redundant listing of qlo points. PROPOSITION 5. The non-redundant listing of the qlo points of depth m q 1 uniquely determines a listing of the triangles of le¨ el 2 m which coincides with the listing of the le¨ el leaf cycle at le¨ el 2 m in the binary simplex tree. Figures 7 and 8 show the paths through the clo points and the qlo points of depth 4.

FIG. 7. Non-redundant path through the clo points.

INTERLACED QUADTREE ALGORITHMS

113

FIG. 8. Non-redundant path through the qlo points.

4.5. Triangle Vertex Codes LEMMA 1. The ¨ ertices of a triangle t of le¨ el 2 m q 1 are v

yky1 a clo point ¨ 0 s Ý m s k , the lattice origin of t, ks0 2

v

a point of the form ¨ 1 s ¨ 0 q 2ym y1smq1 ,

v

q a point of the form ¨ 2 s ¨ 0 q 2ym y1smq1 .

and

1 yky1 The point x Ž t . s Ý mq s k , is a clo point of depth m q 1 located at ks 0 2 the center of the edge of t which joins the vertices ¨ 0 and ¨ 1. The clo code of this point is lŽ t . s Ž j0 , . . . , jmq1 ., where jk s jŽ s k .. We shall call lŽ t . the clo code of the triangle t. Geometrically, if we draw a small arrow from x Ž t . pointing into the triangle t, then the arrow points in a counterclockwise direction around the lattice origin ¨ 0 . Figure 4 shows the arrows attached to clo points associated with the triangles of level 5. For later reference we now picture some important points for a binary simplex in the plane. The triangle vertices, vectors, and points of a typical type 1 triangle are illustrated in Fig. 9. Each binary simplex in the plane is the image of this reference triangle t under a suitable affine transformation consisting of compositions of translations, rotations, dyadic scalings, and quincunx transformation. The corresponding points for triangles of other types and levels are obtained by identifying the corresponding points in the image of t under the appropriate transformation. The vertices are ¨ 0 , ¨ 1 , and ¨ 2 ; the vectors are e1 , e2 , e3 , and s; and the midpoints are m1 ,

114

D. J. HEBERT

FIG. 9. Vertices, midpoints, and vectors of a reference triangle.

m 2 , m 3 , c0 , c1 , and c 2 . The calculations are as follows: mq1

m1 s x Ž t . s

Ý ks0

m

¨0 s

Ý ks0

e1 s

1 2

mq2

s Ž jmq1 . ,

¨ 1 s ¨ 0 q 2 e1 ,

e3 s 12 Ž e1 q e2 . , m1 s ¨ 0 q e 1 , c 0 s ¨ 1 q s,

1 2

1 2

kq1

kq1

sk ,

sk ,

e2 s

1 2

mq2

s Ž jq mq1 . ,

¨ 2 s ¨ 0 y 2 e2 ,

s s 12 Ž e2 y e1 . ,

m 2 s ¨ 0 q e2 , c1 s ¨ 0 q e 3 ,

m 3 s ¨ 0 q 2 e3 , c 2 s m 3 q s.

The point m1 , called the rep point of the triangle, is the clo point which represents the triangle in the correspondence between the non-redundant clo path and the level leaf cycle. The point m 3 , called the split point, is a qlo point which bisects the longest edge and is the quincunx lattice origin

INTERLACED QUADTREE ALGORITHMS

115

for the children of the triangle t. The points c1 and c 2 are the rep points of the children of t in the full binary simplex tree. The lattice origin vertex ¨ 0 is the quadtree parent of m1; so the clo code of ¨ 0 may be obtained from that of m1 by dropping the last digit. The vector s, called the step vector, points in the direction of the lattice origin path link which begins at m1 and crosses the triangle. Thus s points from m1 into the triangle t in counterclockwise direction about the lattice origin ¨ 0 .

5. STACK THEOREM Figure 10 shows the depth 4 clo and qlo paths on the same picture. This figure suggests that we may create an algorithm to step through the two quadtree leaf sets simultaneously, thus visiting the level 2 m and level 2 m q 1 triangles in one sweep. As we step through the paths at a given level we serially visit triangular halves of quadtree blocks. In this section we shall see that as we pass through one triangular half-block we may push local data onto a stack so that when we reach the second half-block the data may be popped from the stack to do a calculation for the whole block. Consider the following procedure: Start with an empty stack. Step through the clo path at depth m. At each turn point x, if lŽ x . Žthe clo code. is not on the top of the stack then push it onto the stack, else pop it off the stack. This section is devoted to the proof of the following theorem

FIG. 10. Interlaced paths at depth 4.

116

D. J. HEBERT

which will allow the construction of algorithms for serially processing data indexed by leaf cycles which, in turn, interlace the clo and qlo paths. THEOREM 1 ŽStack Theorem.. At the second ¨ isit to a turn point x in the clo path, the clo code is on top of the stack. The maximum size of the stack while tra¨ ersing the clo path of le¨ el m is 2 m y 1. The proof of the theorem is based on the following lemma and its corollaries: LEMMA 2 ŽHalf-loop Sublemma.. Ž1. If j k s Ž j0k , j1k , . . . , jnk , a, a, . . . , a., then for some p ) 0, j kq p s Ž j0k , j1k , . . . , jnk , a, b, . . . , b . is a first turn and, for q - p, kq q j kq q s Ž j0k , j1k , . . . , jnk , a, jnq2 , . . . , jmkq q . .

Ž2. p ) 0,

If j k s Ž j0k , j1k , . . . , jnk , a, b, . . . , b . is a second turn, then for some j kq p s Ž j0k , j1k , . . . , jnk , a, a, . . . , a . ,

and, for q - p, kq 1 , . . . , jmkq q . . j kq q s Ž j0k , j1k , . . . , jnk , a, jnq2

Proof. For n s m y 2, statements Ž1. and Ž2. follow from the fact that the following pattern for the pair Ž a, b . s Ž a, a c . which occurs when m s 2: Ž a, a. ª Ž a, aq. ª Ž a, aqq . s Ž a, b ., a first turn, and if Ž a, b . is a second turn, then Ž a, b . ª Ž a, bq. ª Ž a, bqq . s Ž a, a.. Assume Ž1. and Ž2. above as an induction hypothesis, and suppose that j k s Ž j0k , j1k , . . . , k . Ž . jny 1 , a, a, a, . . . , a . Then by 1 there is a p 1 such that k j kq p 1 s Ž j0k , j1k , . . . , jny1 , a, a, b, . . . , b .

is a first turn and, for q - p1 , k kq q j kq q s Ž j0k , j1k , . . . , jny1 , a, a, jnq2 , . . . , jmkq q . .

It follows from the definition of clo path that k , a, aq, bq, . . . , bq . j kq p 1q1 s Ž j0k , j1k , . . . , jny1

is a second turn. From Ž2. it follows that, for some p 2 , k j kq p 1q1 qp 2 s Ž j0k , j1k , . . . , jny1 , a, aq, aq, . . . , aq . .

INTERLACED QUADTREE ALGORITHMS

117

From Ž1., then, there is p 3 such that k j kq p 1q1 qp 2qp 3 s Ž j0k , j1k , . . . , jny1 , a, aq, bq, . . . , bq . ,

a first turn, and k j kq p 1q1 qp 2qp 3q1 s Ž j0k , j1k , . . . , jny1 , a, b, a, . . . , a .

is a second turn. By Ž2. again, there is p4 such that k j kq p 1q1 qp 2qp 3q1 qp 4 s Ž j0k , j1k , . . . , jny1 , a, b, b, . . . , b . ,

and Ž1. holds for n y 1. To show that Ž2. also holds for n y 1 assume that j k s Ž j0k , j1k , . . . , k jny 1 , a, b, b, . . . , b . is a second turn, then by Ž1. for some p1 ) 0, k j kq p 1 s Ž j0k , j1k , . . . , jny1 , a, b, a, . . . , a . ,

a first turn. Then k j kq p 1q1 s Ž j0k , j1k , . . . , jny1 , a, bq, aq, . . . , aq .

is a second turn and by Ž2. there is p 2 such that k j kq p 1q1 qp 2 s Ž j0k , j1k , . . . , jny1 , a, bq, bq, . . . , bq . .

Once more, by Ž1. we have p 3 such that k j kq p 1q1 qp 2qp 3 s Ž j0k , j1k , . . . , jny1 , a, bq, aq, . . . , aq . ,

a first turn, so that k j kq p 1q1 qp 2qp 3q1 s Ž j0k , j1k , . . . , jny1 , a, a, b, b, . . . , b .

is a second turn. Finally, from Ž2. again, there is p4 such that k j kq p 1q1 qp 2qp 3qp 4q1 s Ž j0k , j1k , . . . , jny1 , a, a, a, . . . , a . .

A subpath of the clo path starting and ending as in Ž2. is called a first triangle path of level m y n, and a subpath starting and ending as in Ž1. is called a second triangle path of level m y n. The combination of a second triangle path followed by a first triangle path is called a quadloop Žof level m y n.. The loop sublemma has the following corollaries: COROLLARY 1 ŽA Second Turn Begins a Loop.. If j k s Ž j0k , j1k , . . . , jnk , a, b, . . . , b . is a second turn, then for some p ) 0, j kq p s j k and for q - p, kqq j kq q s Ž j0k , j1k , . . . , jnk , a, jnq2 , . . . , jmkqq ..

118

D. J. HEBERT

COROLLARY 2 ŽTriangle: Second, Loop, First.. A first or second triangle path of le¨ el k begins with a second triangle path of le¨ el k y 1 followed by a quadloop of le¨ el k y 1 followed by a first triangle path of le¨ el k y 1 at its end. COROLLARY 3 ŽLoop: First, Loop, Loop, Loop, Second.. A quadloop of le¨ el k begins with a first triangle path of le¨ el k y 1 followed by three quadloops of le¨ el k y 1 followed by a second triangle path of le¨ el k y 1 at its end. The following two corollaries refer to the setup for the stack theorem at the beginning of this section. COROLLARY 4 ŽStack Growth.. In the tra¨ ersal of a first triangle path of le¨ el k the stack grows by an amount 2 k if the last point of the path is a turn point and by the amount 2 k y 1 otherwise. COROLLARY 5 ŽLoop-Stack.. In tra¨ ersing a quadloop at any le¨ el the stack is the same size at the beginning and at the end. If we push the clo code of the beginning point of a quadloop we will pop the same clo code as we exit. These corollaries either follow from the lemma or can be proved by calculations similar to those in the proof of the lemma. The theorem follows immediately. An analogous lemma and its analogous corollaries can be stated and proved for qlo paths.

6. TURN LISTS 6.1. Properties of the Turn list at a Gi¨ en Le¨ el From the definition of clo path we see that if l is a turning point of the form l s Ž j0 , . . . , jiy1 , ji , . . . , jm . , where ji s ??? s jm and the predecessor of l is not a turning point, then the successor lq is also a turn point and lqq is not a turn point. Hence, in the list of clo codes in the clo path, the turn points come in pairs preceded and followed by non-turn points. The following proposition says that there are only two possible local patterns in the list of turn flags. Indeed, a pair of turn points is followed either by one or three non-turn points before the next pair of turn points. PROPOSITION 6. If L s Ž l k , . . . , l kq5 . is a list of clo codes in the clo path ordering, and if T s Žt k , . . . , t kq5 . is the corresponding list of turnflags with t k s t kq1 s 1, then there are only two possibilities: T s Ž1, 1, 0, 0, 0, 1. or T s Ž1, 1, 0, 1, 1, 0..

INTERLACED QUADTREE ALGORITHMS

119

Proof. Assume that l k s Ž j0 , . . . , jm . where jk s a and jkq1 s ??? s jkq n s jm s b s a c. First we consider the case n s 1. In this case the list of the last three digits of l k is Ž x, a, b .. If x / aq and x / bq then the lists of the last three digits of the entries of L Žin order. are Ž x, a, b ., Ž x, aq, bq ., Ž x, aq, a., Ž x, aq, aq ., Ž x, aq, b ., Ž x, aq, bq ., so that T s Ž1, 1, 0, 0, 0, 1.. If x s aq and t kq 3 s 0 we get the same result. If x s aq and t kq 3 s 1 we get T s Ž1, 1, 0, 1, 1, 0. because turns come in pairs separated by non-turns. If x s bq we have Ž bq, a, b ., Ž bq, aq, bq ., Ž bq, aq, a., Ž bq, aq, aq ., Ž a, b, b ., Ž a, b, bq. so that T s Ž1, 1, 0, 1, 1, 0.. Next we consider the case n ) 1. In this case we have lists of the last n q 1 digits l k to be Ž a, b, . . . , b ., Ž aq, bq, . . . , bq ., Ž aq, bq, . . . , a., Ž aq, bq, . . . , bq, aq ., Ž aq, bq, . . . , a, b ., Ž aq, bq, . . . , a, bq ., so that T s Ž1, 1, 0, 1, 1, 0.. For calculations dependent on local properties of a lattice origin path we will use a list of four local turn flags. For this purpose we define a 4-turn list to be a list of the form Tk s Žt ky1 , t k , t kq1 , t kq2 .. If Tk s Ž0, 1, 1, 0. then l ky 1 is called a preturn point, l k is called a first turn, l kq1 is called a second turn, and l kq 2 is called a post turn. PROPOSITION 7. There are only se¨ en possible 4-turn lists:

Ž 1, 1, 0, 0 . , Ž 1, 0, 0, 0. , Ž 0, 0, 0, 1. , Ž 0, 0, 1, 1. , Ž 0, 1, 1, 0 . , Ž 1, 1, 0, 1 . , Ž 1, 0, 1, 1 . . Proof. By Proposition 6 there are two possible 4-turn lists which begin with a pair of 1s, namely, Tk s Ž1, 1, 0, 0. and Tk s Ž1, 1, 0, 1.. If Tk s Ž1, 1, 0, 0. then the last two zeroes begin a sequence of three so that Tkq 1 , Tkq2 , Tkq3 , Tkq4 are necessarily Ž1, 0, 0, 0., Ž0, 0, 0, 1., Ž0, 0, 1, 1., Ž0, 1, 1, 0., while Tkq 5 is either Ž1, 1, 0, 0. or Ž1, 1, 0, 1.. If Tk s Ž1, 1, 0, 1. then necessarily Tkq 1 s Ž1, 0, 1, 1. and Tkq2 s Ž0, 1, 1, 0., and Tkq3 is either Ž1, 1, 0, 0. or Ž1, 1, 0, 1.. This takes care of all the possibilities. COROLLARY 6. In a lattice origin path the 4-turn lists follow in the following two possible orderings:

Ž 0, 1, 1, 0 . ª Ž 1, 1, 0, 0 . ª Ž 1, 0, 0, 0 . ª Ž 0, 0, 0, 1 . ª Ž 0, 0, 1, 1 . ª Ž 0, 1, 1, 0 . or

Ž 0, 1, 1, 0 . ª Ž 1, 1, 0, 1 . ª Ž 1, 0, 1, 1 . ª Ž 0, 1, 1, 0 . .

120

D. J. HEBERT

The following table numbers and names the types of points for the 4-turn lists: Number

4-turn list

Name

1 2 3 4 5 6 7

Ž0, 1, 1, 0. Ž1, 1, 0, 0. Ž1, 0, 0, 0. Ž0, 0, 0, 1. Ž0, 0, 1, 1. Ž1, 1, 0, 1. Ž1, 0, 1, 1.

First turn point Preloop second turn point First 0 loop point Second 0 loop point Third 0 preturn point Single 0 second turn point Single 0 preturn point

6.2. Triangles and 4-turn Lists We may associate the 4-turn lists with triangles in the nonredundant scanning order by writing TiŽ k . s Tk when l k is not a first turn point. Since we eliminate first turns, by Proposition 7 there are six possibilities for TiŽ k . and by Corollary 6 there are two possible orderings:

Ž 1, 1, 0, 0 . ª Ž 1, 0, 0, 0 . ª Ž 0, 0, 0, 1 . ª Ž 0, 0, 1, 1 . and Ž 1, 1, 0, 1 . ª Ž 1, 0, 1, 1 . . The final list in each of these orderings corresponds to a preturn point and the first list in each ordering is a second turn point. If l k is a preturn point then we have TiŽ k .q1 s TiŽ kq2. , which may be either Ž1, 1, 0, 1. or Ž1, 1, 0, 0. depending on the value of t Ž kq4. . To compute t Ž kq4. , however, we do not have to leave the triangle t iŽ k . , as we shall see in Proposition 8, Section 8. 7. TURN SEQUENCES FOR INTERLACED PATHS In this section we list some observations of importance in constructing algorithms which simultaneously step through the clo and qlo paths. ŽFor example, vertex listing algorithms, spline and wavelet construction algorithms, refining and pruning algorithms.. Since the observations hold both for type 1 triangles associated with clo points and paths Žwhich we shall call type 1 lo points and paths. and for type 0 triangles associated with qlo points and paths Žwhich we shall call type 0 lo points and paths. we henceforth refer to type A and type B triangles, points, and paths where A, B g  0, 14 and A / B. The fact that we may simultaneously step through type A paths of level l and type B paths of level l y 1 follows from the fact that if m1 is the rep point of a type A triangle t of level l, then m 3 is a type B lo point. If m 3 is

INTERLACED QUADTREE ALGORITHMS

121

not a first turn point in its type B lo path then it is the rep point of a type B triangle of level l y 1. Stepping through the type A lo path implies stepping through a set of type B lo points. Fortunately, by Proposition 5, as we step through the rep points m1Ž k . of the type A triangles t k the split points m 3 Ž k . visit all the type B lo points in the lo path ordering. 7.1. Some Notation Ly 1 To make more precise statements we need some notation. Let  l k 4ks0 Ly 1 be the list of type A lo codes for a type A lo path. Let t k 4ks0 be the corresponding list of turnflags. Let iŽ k . be the cycle index of the triangle whose rep point is l k , so that if t k s t kq1 s 1 Žwhich means that l k is eliminated as redundant in the list of triangle rep points. we have iŽ k . s iŽ k y 1.. If m1Ž k . s m1Ž t iŽ k . . is the rep point for the triangle t iŽ k . of cycle index iŽ k ., and l k is not a first turn point, we have m1Ž k . s l k . Let m 3 Ž k . be the split point of the triangle t iŽ k . . Since m 3 Ž k . is a type B lo point let kX Ž k . denote its index in its type B lo path at level l y 1. We then have lXk X s mX1Ž kX . s m 3 Ž k ., where kX s kX Ž k ., lXk X is the lo code of the level l y 1 triangle of index kX and mX1Ž kX . denotes the rep point of the appropriate triangle of level l y 1. We also use t kX X to denote the turnflag ˆˆk be the level l q 1 type B lo code of index ˆk whose of lXk X . Finally, let l turnflag is tˆˆk Žcorner flags are included as turnflags in the list of qlo points.. If m1Ž k . s m 3 Ž ˆ k . then we write k s k Ž ˆ k .. The following lemma lists some properties of the turnlists at levels l y 1, l, and l q 1:

LEMMA 3. If t k s 0, where k s k Ž ˆ k ., then ˆ t k s 0, tˆkq1 s ˆ t kq2 s 1, X ˆ ˆ ˆ Ž . Ž . Ž . and ˆ t kq 3 s 0, and k k s k k q 1 s k k q 2 . If t k s 1, where k s k Ž ˆ k ., then ˆ t k s 0, k Ž ˆ k . s kŽ ˆ k y 1. q 1, and k Ž ˆ kq ˆ . Ž . 1 s k k q 1. If t k s t kq1 s t kq2 s 0, then t kX X s t kX Xq1 s 1, t kX Xq2 s 0, kX s kX Ž k ., X k q 1 s kX Ž k q 1., and kX q 2 s kX Ž k q 2.. If t k s t kq1 s 1, then m 3 Ž k y 1. s m 3 Ž k . s m 3 Ž k q 1. s mX1Ž kX Ž k y 1.. and kX Ž k y 1. s kX Ž k . s kX Ž k q 1.. Table 1 shows a typical pattern of turnflags at levels l y 1, l, and l q 1: l q 1 00011 00011 00011 011 011 00011 00011 00011 011 011 00011 011 0 l 110 110 110 0 0 110 110 110 0 0 110 0 0 l y 1 y0 y0 y1 1 0 y0 y0 y1 1 0 y1 1 0

122

D. J. HEBERT

7.2. Obser¨ ations about m1Ž k . and m 3 Ž k . If m1Ž k . is a preturn point then m 3 Ž k . s m 3 Ž k q 2. is Ža. the common split point of the triangles t iŽ k . s t iŽ kq1. and t iŽ kq2. , Žb. the rep point of the parent of t iŽ kq2. which is its first child, Žc. not a first turn point Žsince it is a rep point.. v

mŽ k . is a first loop point iff m 3 Ž k . is a first turn point.

v

m1Ž k . is a second loop point iff m 3 Ž k . is a second turn point.

If m1Ž k . is a preturn then it is the second child of its parent Žtherefore, has parity 1.. v

v

A first loop point has parity 1 and an second loop point has parity 0.

v

A first turn point has parity 1 and a second turn point has parity 0.

Finally, we note that a triangle of level 2 m has a vertex ¨ 0 a qlo point at depth m and as a rep point a qlo point at depth m q 1. A triangle of level 2 m y 1 has as vertex ¨ 0 a clo point at depth m and as a rep point m1 a clo point at depth m q 1.

8. AN INTERLACED QUADTREE ALGORITHM 8.1. Setting The general setting for a cyclic quadtree interlace algorithm is as follows: We associate with the triangle t il of cyclic index i at level l a list of data l h i , and with the quadtree block l k of lo path index k at level l a list of data w kl . We also assume that there are functions H and W such that lq1 . l Ž l l . where t 2lq1 h li s H Ž h lq1 and t 2lq1 2 i , h 2 iq1 and w k s W h i , h iU i iq1 are the l l l children of t i , and where t i and t iU are the triangular halves of the quadtree block l k centered at m 3 Ž t il .. An interlaced quadtree algorithm starts with coefficients h li and steps through the cycle index values i, computing h ly1 and pushing both h ly1 and h li to stacks for later usage in i n ly1 ly1 l l calculating w k s w kU and wn s wnU as the indices kU and nU of the second visit to the quadtree blocks of indices k and n are reached in the stepping procedure. In our examples the function H represents the computation of half wavelets from local triangle data, and the function W computes the wavelets from the half wavelet data on a stack and the local data.

INTERLACED QUADTREE ALGORITHMS

123

8.2. Algorithm to Compute m1Ž k ., m 3 Ž k ., and Tk from m1Ž k y 1., m 3 Ž k y 1., and Tky1 An important part of a quadtree interlace algorithm is the simultaneous stepping through the clo and qlo paths with guidance from the 4-turn lists. We shall find that as we step through the lattice origin path at level l a small computation obtains, for each lo path index k, the triangle rep point m1Ž k ., the split point m 3 Ž k ., and the 4-turn list Tk from their values at the index k y 1. Indeed, except in the preturn cases, the vertices m1Ž k . and m 3 Ž k . are obtained by simply incrementing the last digit of the lo code, and Tk is obtained from Tky1 by referring to the list sequences in Corollary 6. The preturn cases are handled by the following proposition: PROPOSITION 8. If t iŽ ky1. is a preturn point then Tkq1 s Ž1, 1, 0, d 4 ., where d 4 s t Ž m 3 Ž k .q. c. Proof. By Corollary 6 the incrementation of Tky 1 is automatic except in the case of a preturn, i.e., Tky 1 s Ž1, 0, 1, 1. or Tky1 s Ž0, 0, 1, 1.. In these cases the triangles t iŽ ky1. and t iŽ kq1. share an edge, m 3 Ž k y 1. s Ž m 3 Ž k q 1. and mq 1 is not the next triangle rep point for example, refer to . Ž . Fig. 11 . Thus Tkq 1 is of the form 1, 1, 0, d k , where d k is the turnflag of m1Ž k y 1.qqqq. Note that m1Ž k y 1.qqq is either a first 0 loop point or a single 0 preturn point. By Table 1 Žin Section 7.1. we see that the turnflag of its successor is the complement of the turnflag of m 3 Ž k y 1.q. ALGORITHM. If Tny 1 s Ž0, 0, 1, 1. or Tny1 s Ž1, 0, 1, 1. Žpreturn case. then m 3 Ž n. s m 3 Ž n y 1.; m1Ž n. s incturn ŽinclŽ m1Ž n y 1...; d 4 s t ŽinclŽ m 3 Ž n... c ; Tn s Ž1, 1, 0, d 4 .;

FIG. 11. The cases where t ky 1 is a preturn.

124

D. J. HEBERT

else m1Ž n. s inclŽ m1Ž n y 1.. If Tny 1 s Ž1, 0, 0, 0. then Tn s Ž0, 0, 0, 1.; m 3 Ž n. s incturn Ž m 3 Ž n y 1.. else m 3 Ž n. s inclŽ m 3 Ž n y 1.. If Tny 1 s Ž1, 1, 0, 0. then Tn s Ž1, 0, 0, 0.; If Tny 1 s Ž0, 0, 0, 1. then Tn s Ž0, 0, 1, 1.; If Tny 1 s Ž1, 1, 0, 1. then Tn s Ž1, 0, 1, 1.. Notice that this algorithm requires only a few computer operations, at most five logical decisions, three digit increments, and an lo code turn increment. 8.3. Stack Algorithm We consider now the stack usage by focusing on the step in which we calculate Tk from Tky1. Here we have i s iŽ k . s iŽ k y 1. q 1. We use a stack called stack0 to store h i and m 3 Ž t i . for later computation of wm 3 s w km Žwhere 2 j s i y 1. and ¨ 0 Ž t i . for computing w¨ . and stack1 to store h ly1 j 0 To indicate a push or a pop of the pair Ž m, h. for stack0 we write push0Ž m,h. or pop0Ž m, h.. To indicate a non-destructive reading of the top of the stack we write top0 s Ž m, h. or top1 s Ž m, h.. For each k we perform the following stack test: First we check for a possible calculation of wm 3 by examining the top of stack0 to see whether m 3 is present. If top0 s Ž m 3 , h. then we use pop0 to remove the top and compute wm 3 s W Ž h, h li .. Otherwise we use push0 to store Ž m 3 Ž k ., h li . on stack 0. In the case that m1Ž k y 1. is a preturn point we find top0 s Ž m 3 , h liy1 . so that the stack usage could be skipped. l The values h ly1 s H Ž jly1 , h li ., where 2 j s i y 1, are calculated when i j Ž is odd so that t iy1 and t i have the same parent.. The possible values of Tky 1 in these cases are Ž1, 1, 0, 0., Ž0, 0, 0, 1., and Ž1, 1, 0, 1.. We then do a stack test on stack1. If top1 s Ž ¨ 0 Ž t i ., h. then we pop1 and compute ., Again, in the two cases Ž1, 1, 0, x . we may skip the stack w¨ 0 s W Ž h, h ly1 j ly1 . usage and compute w¨ 0 s W Ž h ly1 directly. iy1 , h j The stack-usage part of the algorithm, without the skipping of stack usage in special cases may be described briefly as follows: ALGORITHM. For each k, top0 s Ž m, h. If m s m 3 Ž k . then pop0Ž m, h.; wm 3 s W Ž h, h li .. else push0Ž m 3 Ž k ., h li .; If i mod 2 s 1 then j s Ž i y 1.r2, h ly1 s H Ž h lly1 , h li ., where 2 j s i y j 1;

INTERLACED QUADTREE ALGORITHMS

125

top1 s Ž ¨ , h.; .. If ¨ s ¨ Ž t i .; popŽ ¨ , h.; w¨ 0 s W Ž h, h ly1 i

9. IMPLEMENTATION OF THE SQUARE HAAR WAVELET AND THE BOX SPLINE PYRAMID ALGORITHMS We now have the notation and results necessary to describe implementations of the interlaced quadtree algorithms to computations in the examples of Section 2. 9.1. The Edge Detecting Haar Wa¨ elets We wish to encode the level l triangle intensities in terms of the wavelet coefficients wlj, 1 and wlj, 3 for j F l. As we reach a2j i in the traversal of the leaf cycle we compute d ijy1 and a ijy1 from a2j i and a2j iq1 pushing m 3 Ž k . and h li s Ž d ijy1, a2j i , a2j iq1 . onto a stack, discarding a2j i and a2j iq1 , and saving a ijy1 in temporary storage for use in the level j y 1 calculations. When we reach the square again at index iU we pop d ijy1 from the stack U U using it along with d ijy1 to compute Ž wlj, 1 , wlj, 3 .. We also compute a ijy1 to put in temporary storage. This procedure is continued for each level j until level 0 is reached, and we must save a00 and a10 . The triangle intensities a ij may be reconstructed from the coefficients 0 a0 , a10 and wlj, 1 , wlj, 3 for j G 0 by the following procedure: Given a ijy1, wlj, 1 and wlj, 3 we push a ijy1 to a stack and pop it again at U to compute wlj, 2 and slj. Next we compute index iU using a ijy1 and a ijy1 the appropriate version of a s By1 w. Thus, the list aT s Ž a2j i , a2j iq1 , a2j iU , a2j iU q1 . is stored out of order in memory to be rearranged later. To compute the next level coefficients a ijq1 we traverse the level j q 1 leaf cycle in reverse order. When we reach a2j iU we push it to the stack for later jq1 computation of wljq1 but we also remove a2j iU from its memory , 2 and sl position and push it onto another stack to be popped again when we reach its place in the leaf cycle. 9.2. Box Spline Computations As mentioned in Section 2, the edge detection computations for the box spline wavelets of Ron and Shen w17x are the same as those of the square Haar wavelets above. We now look more closely at the simplicity of calculation of the four directional box spline coefficients in terms of the interlaced quadtree triangle vertex computations. In particular, we consider the representation of a function defined on a uniform rectangular mesh, in terms of coefficients of the well-known Zwart]Powell elements

126

D. J. HEBERT

scaled as in the multiresolution analysis studied, for example, in w14x and w29x. In the spirit of the previous examples, we consider again an image consisting of 2 m = 2 m pixel intensities which are identified with the clo points of depth m listed in the non-redundant ordering as  skl 4 , with l s 2 m y 1. We calculate for each level l, a set of coefficients skly1 which represent the coefficient of a discrete box spline centered at the kth lo point at level l y 1. At each level l we again identify the lo point of level l with the point m1Ž i . of the ith triangle of level l. Here is the algorithm: ALGORITHM 3. Ž1. Assign the value 0 to each of the vertices of triangles of level l. For each leaf cycle index i at level l, add the value Ž1r4. sil to the values assigned to the vertices ¨ 0 Ž i . and ¨ 1Ž i . of the triangle of index i. Ž2. If m1Ž k . s m 3 Ž i ., define 3

skly1 s

1 4

Ý bj , js0

where bj is the value assigned to the jth corner vertex of the quadtree block centered at m1Ž k . s m 3 Ž i . Žthe m1 point of the triangle of level l y 1 and index k .. The computations of this algorithm are accomplished by a quadtree interlace algorithm. Step Ž1. is performed by traversing the leaf cycle of level l. The vertices ¨ 0 Ž i . and ¨ 1Ž i . Žexcept on the edges of the square. are lo points representing triangles of levels greater than l y 1. Step Ž2. requires a second traversal of the leaf cycle at level l to complete the computation of skly1. For each i, if b and c are the values stored at the vertices ¨ 0 Ž i . and ¨ 1Ž i . then we add Ž1r4.Ž b q c . to the value at the skly1 at the point m 3 Ž i .. The algorithm computes 3

skly1 s

1 4

3

Ý bj s 161 Ý js0

s jq ,

j, qs0

where s jq s sil for some i. If skly1 is the intensity assigned to a qlo point, then the algorithm computes skly1 according to the mask of the discrete Zwart]Powell element given in Section 2 applied to the intensities sil corresponding to lo points which surround the point m1Ž k . at the center of the mask. When the algorithm is applied at successive levels starting at level 2 m y 1 we obtain coefficients sil for l - 2 m y 2 as discrete box spline masks which turn out to be the same masks which occur in the subdivision algorithm for the Zwart spline Žcf. w14x..

INTERLACED QUADTREE ALGORITHMS

127

REFERENCES 1. 2. 3. 4. 5.

6. 7. 8. 9. 10.

11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.

H. S. M. Coxeter, ‘‘Regular Polytopes,’’ Dover, New York, 1973. K. Georg, Generation of triangulations by reflections, Utilitas Math. 16 Ž1979., 123]129. H. Sagan, ‘‘Space Filling Curves,’’ Springer-Verlag, BerlinrNew York, 1994. M. C. Rivara, Mesh refinement processes based on the generalized bisection of simplices, SIAM J. Numer. Anal. 21 Ž1984., 604]613. J. M. Maubach, ‘‘Iterative Methods for Non-linear Partial Differential Equations,’’ N.W.O. Center for Mathematics and Computer Science, Amsterdam, 1991 ŽISBN 909004007-2.. J. M. Maubach, Local bisection refinement for n-simplicial grids generated by reflections, SIAM J. Sci. Comput. 16 Ž1995., 210]227. D. J. Hebert, Symbolic local refinement of tetrahedral grids, J. Symbolic Comput. 17 Ž1994., 457]472. H. Samet, ‘‘Applications of Spatial Data Structures: Computer Graphics, Image Processing, and GIS,’’ Addison-Wesley, Reading, MA, 1989. G. Dutton, Locational properties of quaternary triangular meshes, in ‘‘Proceedings of the Fourth International Symposium on Spatial Data Handling, Zurich, 1990,’’ pp. 901]910. M. F. Goodchild and Y. Shiren, A hierarchical spatial data structure for global geographic information systems, CVGIP: Graphical Models and Image Processing, 54 Ž1992., 33]41. G. Fekete, Rendering and Managing Spherical Data with Sphere Quadtrees, in ‘‘Proceedings of Visualization 90, San Francisco, 1990.’’ H. K.-C. Chang and C.-Y. Liou, Parallel implementation of linear quadtree codes using the n Cube 2 supercomputer system, Internat. J. Supercomputer Appl. 9 Ž1995., 220]231. P. Zwart, Multivariate splines with nondegenerate partitions, SIAM J. Numer. Anal. 10 Ž1973., 665]673. C. de Boor, K. Hollig, and S. Riemenschnieder, ‘‘Box Splines,’’ Springer-Verlag, ¨ BerlinrNew York, 1993. I. Daubechies, ‘‘Ten Lectures on Wavelets a61,’’ CBMS]NSF Conference Series in Applied Mathematics, SIAM, Philadelphia, 1992. A. Cohen and I. Daubechies, Non-separable bidimensional wavelet bases, Re¨ . Math. Iberoamericana 9 Ž1993., 51]137. A. Ron and Z. Shen, Compactly supported tight affine frames in L2 Ž R d ., Math. Comp., to appear. E. P. Simoncelli and E. H. Adelson, Non-separable extensions of quadrature mirror filters to multiple dimensions, Proc. IEEE 78 Ž1990., 652]663. J. Kovacevic and M. Vetterli, Nonseparable multidimensional perfect reconstruction filter banks and wavelet bases for R n , IEEE Trans. Inform. Theory. 28 Ž1992., 533]555. D. J. Hebert, A box spline subdivision pyramid algorithm, Appl. Math. Lett., to appear. D. J. Hebert and H. J. Kim, Image encoding with triangulation wavelets, Proc. SPIE 2569 Ž1995., 381]392. J. A. Provine and R. M. Ragayyan, Lossless compression of peanoscanned images, J. Electronic Imaging 3 Ž1994., 176]181. H. J. Kim, ‘‘Unified Image Compression Using Reversible and Fast Biorthogonal Wavelet Transforms,’’ thesis, University of Pittsburgh, 1996. G. Strang and T. Nguyen, ‘‘Wavelets and Filter Banks,’’ Wellesley-Cambridge Press, 1996. E. R. Dougherty and R. Giardina, ‘‘Image Processing}Continuous to Discrete,’’ Vol. 1, 1987.

128

D. J. HEBERT

26. D. J. Hebert and H. J. Kim, A fast-wavelet compass edge detector, Proc. SPIE 2825 Ž1996., 432]442. 27. A. P. Schroder ¨ and W. Sweldens, Spherical wavelets: Efficiently representing functions on the sphere, Computer Graphics Proceedings ŽSIGGRAPH 95. pp. 161]172, ACM, 1995. 28. M. Girardi and W. Sweldens, A new class of unbalanced Haar wavelets that form an unconditional basis for Lp on general measure spaces, J. Fourier Anal. Appl., to appear. 29. C. K. Chui, K. Jetter, and J. Stockler, Wavelets and frames on the four-directional mesh, ¨ in ‘‘Wavelets, Theory, Algorithms, and Applications,’’ ŽC. K. Chui, L. Montejusco, and L. Puccio, Eds.., pp. 213]230, 1994.