Utah Raster Toolkit  9999-git
URT Development version (post-3.1b)
hilbert.c
Go to the documentation of this file.
1 /*
2  * This software is copyrighted as noted below. It may be freely copied,
3  * modified, and redistributed, provided that the copyright notice is
4  * preserved on all copies.
5  *
6  * There is no warranty or other guarantee of fitness for this software,
7  * it is provided solely "as is". Bug reports or fixes may be sent
8  * to the author, who may or may not act on them as he desires.
9  *
10  * You may not include this software in a program or other software product
11  * without supplying the source, or without informing the end-user that the
12  * source is available for no extra charge.
13  *
14  * If you modify this software, you should include a notice giving the
15  * name of the person performing the modification, the date of modification,
16  * and the reason for such modification.
17  */
18 /*
19  * hilbert.c - Computes Hilbert curve coordinates from position and v.v.
20  *
21  * Author: Spencer W. Thomas
22  * EECS Dept.
23  * University of Michigan
24  * Date: Thu Feb 7 1991
25  * Copyright (c) 1991, University of Michigan
26  */
27 static char rcsid[] = "$Header: /l/spencer/src/urt/lib/RCS/hilbert.c,v 3.0.1.1 1992/04/30 14:07:11 spencer Exp $";
28 
29 /*
30  * Lots of tables to simplify calculations. Notation: p#i means bit i
31  * in byte p (high order bit first).
32  * p_to_s: Output s is a byte from input p such that
33  * s#i = p#i xor p#(i-1)
34  * s_to_p: The inverse of the above.
35  * p_to_J: Output J is "principle position" of input p. The
36  * principle position is the last bit s.t.
37  * p#J != p#(n-1) (or n-1 if all bits are equal).
38  * bit: bit[i] == (1 << (n - i))
39  * circshift: circshift[b][i] is a right circular shift of b by i
40  * bits in n bits.
41  * parity: Parity[i] is 1 or 0 depending on the parity of i (1 is odd).
42  * bitof: bitof[b][i] is b#i.
43  * nbits: The value of n for which the above tables have been
44  * calculated.
45  */
46 
47 typedef unsigned int byte;
48 
49 static int nbits = 0;
50 
51 static byte
52  p_to_s[512],
53  s_to_p[512],
54  p_to_J[512],
55  bit[9],
56  circshift[512][9],
57  parity[512],
58  bitof[512][9];
59 
60 /* Calculate the above tables when nbits changes. */
61 static void
63 int n;
64 {
65  register int i, b;
66  int two_n = 1 << n;
67 
68  if ( nbits == n )
69  return;
70 
71  nbits = n;
72  /* bit array is easy. */
73  for ( b = 0; b < n; b++ )
74  bit[b] = 1 << (n - b - 1);
75 
76  /* Next, do bitof. */
77  for ( i = 0; i < two_n; i++ )
78  for ( b = 0; b < n; b++ )
79  bitof[i][b] = (i & bit[b]) ? 1 : 0;
80 
81  /* circshift is independent of the others. */
82  for ( i = 0; i < two_n; i++ )
83  for ( b = 0; b < n; b++ )
84  circshift[i][b] = (i >> (b)) |
85  ((i << (n - b)) & (two_n - 1));
86 
87  /* So is parity. */
88  parity[0] = 0;
89  for ( i = 1, b = 1; i < two_n; i++ )
90  {
91  /* Parity of i is opposite of the number you get when you
92  * knock the high order bit off of i.
93  */
94  if ( i == b * 2 )
95  b *= 2;
96  parity[i] = !parity[i - b];
97  }
98 
99  /* Now do p_to_s, s_to_p, and p_to_J. */
100  for ( i = 0; i < two_n; i++ )
101  {
102  int s;
103 
104  s = i & bit[0];
105  for ( b = 1; b < n; b++ )
106  if ( bitof[i][b] ^ bitof[i][b-1] )
107  s |= bit[b];
108  p_to_s[i] = s;
109  s_to_p[s] = i;
110 
111  p_to_J[i] = n - 1;
112  for ( b = 0; b < n; b++ )
113  if ( bitof[i][b] != bitof[i][n-1] )
114  p_to_J[i] = b;
115  }
116 }
117 
118 /*****************************************************************
119  * TAG( hilbert_i2c )
120  *
121  * Convert an index into a Hilbert curve to a set of coordinates.
122  * Inputs:
123  * n: Number of coordinate axes.
124  * m: Number of bits per axis.
125  * r: The index, contains n*m bits (so n*m must be <= 32).
126  * Outputs:
127  * a: The list of n coordinates, each with m bits.
128  * Assumptions:
129  * n*m < (sizeof r) * (bits_per_byte), n <= 8, m <= 9.
130  * Algorithm:
131  * From A. R. Butz, "Alternative Algorithm for Hilbert's
132  * Space-Filling Curve", IEEE Trans. Comp., April, 1971,
133  * pp 424-426.
134  */
135 void
137 int n;
138 int m;
139 long int r;
140 int a[];
141 {
142  byte rho[9], rh, J, sigma, tau,
143  sigmaT, tauT, tauT1 = 0, omega, omega1 = 0, alpha[9];
144  register int i, b;
145  int Jsum;
146 
147  /* Initialize bit twiddle tables. */
148  calctables(n);
149 
150  /* Distribute bits from r into rho. */
151  for ( i = m - 1; i >= 0; i-- )
152  {
153  rho[i] = r & ((1 << n) - 1);
154  r >>= n;
155  }
156 
157  /* Loop over bytes. */
158  Jsum = 0;
159  for ( i = 0; i < m; i++ )
160  {
161  rh = rho[i];
162  /* J[i] is principle position of rho[i]. */
163  J = p_to_J[rh];
164 
165  /* sigma[i] is derived from rho[i] by exclusive-oring adjacent bits. */
166  sigma = p_to_s[rh];
167 
168  /* tau[i] complements low bit of sigma[i], and bit at J[i] if
169  * necessary to make even parity.
170  */
171  tau = sigma ^ 1;
172  if ( parity[tau] )
173  tau ^= bit[J];
174 
175  /* sigmaT[i] is circular shift of sigma[i] by sum of J[0..i-1] */
176  /* tauT[i] is same circular shift of tau[i]. */
177  if ( Jsum > 0 )
178  {
179  sigmaT = circshift[sigma][Jsum];
180  tauT = circshift[tau][Jsum];
181  }
182  else
183  {
184  sigmaT = sigma;
185  tauT = tau;
186  }
187 
188  Jsum += J;
189  if ( Jsum >= n )
190  Jsum -= n;
191 
192  /* omega[i] is xor of omega[i-1] and tauT[i-1]. */
193  if ( i == 0 )
194  omega = 0;
195  else
196  omega = omega1 ^ tauT1;
197  omega1 = omega;
198  tauT1 = tauT;
199 
200  /* alpha[i] is xor of omega[i] and sigmaT[i] */
201  alpha[i] = omega ^ sigmaT;
202  }
203 
204  /* Build coordinates by taking bits from alphas. */
205  for ( b = 0; b < n; b++ )
206  {
207  register int ab, bt;
208  ab = 0;
209  bt = bit[b];
210  /* Unroll the loop that stuffs bits into ab.
211  * The result is shifted left by 9-m bits.
212  */
213  switch( m )
214  {
215  case 9: if ( alpha[8] & bt) ab |= 0x01;
216  case 8: if ( alpha[7] & bt) ab |= 0x02;
217  case 7: if ( alpha[6] & bt) ab |= 0x04;
218  case 6: if ( alpha[5] & bt) ab |= 0x08;
219  case 5: if ( alpha[4] & bt) ab |= 0x10;
220  case 4: if ( alpha[3] & bt) ab |= 0x20;
221  case 3: if ( alpha[2] & bt) ab |= 0x40;
222  case 2: if ( alpha[1] & bt) ab |= 0x80;
223  case 1: if ( alpha[0] & bt) ab |= 0x100;
224  }
225  a[b] = ab >> (9 - m);
226  }
227 }
228 
229 /*****************************************************************
230  * TAG( hilbert_c2i )
231  *
232  * Convert coordinates of a point on a Hilbert curve to its index.
233  * Inputs:
234  * n: Number of coordinates.
235  * m: Number of bits/coordinate.
236  * a: Array of n m-bit coordinates.
237  * Outputs:
238  * r: Output index value. n*m bits.
239  * Assumptions:
240  * n*m <= 32, n <= 8, m <= 9.
241  * Algorithm:
242  * Invert the above.
243  */
244 void
246 int n;
247 int m;
248 int a[];
249 long int *r;
250 {
251  byte rho[9], J, sigma, tau,
252  sigmaT, tauT, tauT1 = 0, omega, omega1 = 0, alpha[9];
253  register int i, b;
254  int Jsum;
255  long int rl;
256 
257  calctables(n);
258 
259  /* Unpack the coordinates into alpha[i]. */
260  /* First, zero out the alphas. */
261  switch ( m ) {
262  case 9: alpha[8] = 0;
263  case 8: alpha[7] = 0;
264  case 7: alpha[6] = 0;
265  case 6: alpha[5] = 0;
266  case 5: alpha[4] = 0;
267  case 4: alpha[3] = 0;
268  case 3: alpha[2] = 0;
269  case 2: alpha[1] = 0;
270  case 1: alpha[0] = 0;
271  }
272 
273  /* The loop that unpacks bits of a[b] into alpha[i] has been unrolled.
274  * The high-order bit of a[b] has to go into alpha[0], so we pre-shift
275  * a[b] so that its high-order bit is always in the 0x100 position.
276  */
277  for ( b = 0; b < n; b++ )
278  {
279  register int bt = bit[b], t = a[b] << (9 - m);
280 
281  switch (m)
282  {
283  case 9: if ( t & 0x01 ) alpha[8] |= bt;
284  case 8: if ( t & 0x02 ) alpha[7] |= bt;
285  case 7: if ( t & 0x04 ) alpha[6] |= bt;
286  case 6: if ( t & 0x08 ) alpha[5] |= bt;
287  case 5: if ( t & 0x10 ) alpha[4] |= bt;
288  case 4: if ( t & 0x20 ) alpha[3] |= bt;
289  case 3: if ( t & 0x40 ) alpha[2] |= bt;
290  case 2: if ( t & 0x80 ) alpha[1] |= bt;
291  case 1: if ( t & 0x100 ) alpha[0] |= bt;
292  }
293  }
294 
295  Jsum = 0;
296  for ( i = 0; i < m; i++ )
297  {
298  /* Compute omega[i] = omega[i-1] xor tauT[i-1]. */
299  if ( i == 0 )
300  omega = 0;
301  else
302  omega = omega1 ^ tauT1;
303 
304  sigmaT = alpha[i] ^ omega;
305  /* sigma[i] is the left circular shift of sigmaT[i]. */
306  if ( Jsum != 0 )
307  sigma = circshift[sigmaT][n - Jsum];
308  else
309  sigma = sigmaT;
310 
311  rho[i] = s_to_p[sigma];
312 
313  /* Now we can get the principle position. */
314  J = p_to_J[rho[i]];
315 
316  /* And compute tau[i] and tauT[i]. */
317  /* tau[i] complements low bit of sigma[i], and bit at J[i] if
318  * necessary to make even parity.
319  */
320  tau = sigma ^ 1;
321  if ( parity[tau] )
322  tau ^= bit[J];
323 
324  /* tauT[i] is right circular shift of tau[i]. */
325  if ( Jsum != 0 )
326  tauT = circshift[tau][Jsum];
327  else
328  tauT = tau;
329  Jsum += J;
330  if ( Jsum >= n )
331  Jsum -= n;
332 
333  /* Carry forth the "i-1" values. */
334  tauT1 = tauT;
335  omega1 = omega;
336  }
337 
338  /* Pack rho values into r. */
339  rl = 0;
340  for ( i = 0; i < m; i++ )
341  rl = (rl << n) | rho[i];
342  *r = rl;
343 }
344 
345 
346 #ifdef test
347 #include <stdio.h>
348 
349 main()
350 {
351  int a[9], n, m, i;
352  long int r, r1;
353 
354  for (;;)
355  {
356  printf( "Enter n, m: " );
357  scanf( "%d", &n );
358  if ( n == 0 )
359  break;
360  scanf( "%d", &m );
361  while ( (i = getchar()) != '\n' && i != EOF )
362  ;
363  if ( i == EOF )
364  break;
365  for ( r = 0; r < 1 << (n*m); r++ )
366  {
367  hilbert_i2c( n, m, r, a );
368  if ( n*m <= 6 )
369  {
370  printf( "a = " );
371  for ( i = 0; i < n; i++ )
372  printf( "0x%0*x ", (m+3)/4, a[i] );
373  }
374  hilbert_c2i( n, m, a, &r1 );
375  if ( n*m <= 6 )
376  printf( "r = 0x%0*x\n", (n*m+3)/4, r1 );
377  else if ( r != r1 )
378  printf( "r = 0x%0*x; r1 = 0x%0*x\n", (n*m+3)/4, r,
379  (n*m+3)/4, r1 );
380  }
381  }
382 }
383 
384 p1t( tbl, n )
385 byte tbl[];
386 int n;
387 {
388  int i;
389 
390  for ( i = 0; i < n; i++ )
391  printf( "%02x ", tbl[i] );
392  putchar( '\n' );
393 }
394 
395 p2t( tbl, n, m )
396 byte tbl[][9];
397 int n;
398 {
399  int i;
400 
401  for ( i = 0; i < n; i++ )
402  {
403  printf( "%3d: ", i );
404  p1t( tbl[i], m );
405  }
406 }
407 
408 
409 #endif
static byte circshift[512][9]
Definition: hilbert.c:52
static byte p_to_J[512]
Definition: hilbert.c:52
void hilbert_i2c(int n, int m, long int r, a)
Definition: hilbert.c:136
void hilbert_c2i(int n, int m, a, long int *r)
Definition: hilbert.c:245
static byte s_to_p[512]
Definition: hilbert.c:52
static byte bitof[512][9]
Definition: hilbert.c:52
static byte bit[9]
Definition: hilbert.c:52
unsigned int byte
Definition: hilbert.c:47
static byte parity[512]
Definition: hilbert.c:52
static int nbits
Definition: hilbert.c:49
static void calctables(int n)
Definition: hilbert.c:62
static byte p_to_s[512]
Definition: hilbert.c:52
static char rcsid[]
Definition: hilbert.c:27