images/hilbert.h

   1 /*
   2  *      Image Library -- multidimensional Hilbert curves
   3  *
   4  *      (c) 2006 Pavel Charvat <pchar@ucw.cz>
   5  *
   6  *      This software may be freely distributed and used according to the terms
   7  *      of the GNU Lesser General Public License.
   8  *
   9  *
  10  *      References:
  11  *      - http://www.dcs.bbk.ac.uk/~jkl/mapping.c
  12  *        (c) 2002 J.K.Lawder
  13  *      - J.K. Lawder. Calculation of Mappings between One and n-dimensional Values
  14  *        Using the Hilbert Space-Filling Curve. Technical Report JL1/00, Birkbeck
  15  *        College, University of London, 2000.
  16  *
  17  *      FIXME:
  18  *      - the algorithm fails for some combinations of HILBERT_DIM and HILBERT_ORDER,
  19  *        but it should be safe for HILBERT_DIM = 2..8, HILBERT_ORDER = 8..32
  20  *      - clean and optimize the code
  21  */
  22
  23 #ifndef HILBERT_PREFIX
  24 #  error Undefined HILBERT_PREFIX
  25 #endif
  26
  27 #define P(x) HILBERT_PREFIX(x)
  28
  29 /*
  30  * HILBERT_DIM is the number of dimensions in space through which the
  31  * Hilbert Curve passes.
  32  * Don't use this implementation with values for HILBERT_DIM of > 31!
  33  * Also, make sure you use a 32 bit compiler!
  34  */
  35 #ifndef HILBERT_DIM
  36 #  define HILBERT_DIM 2
  37 #endif
  38
  39 #ifndef HILBERT_TYPE
  40 #  define HILBERT_TYPE u32
  41 #endif
  42
  43 #ifndef HILBERT_ORDER
  44 #  define HILBERT_ORDER (8 * sizeof(HILBERT_TYPE))
  45 #endif
  46
  47 typedef HILBERT_TYPE P(t);
  48
  49 /*
  50  * retained for historical reasons: the number of bits in an attribute value:
  51  * effectively the order of a curve
  52  */
  53 #define NUMBITS HILBERT_ORDER
  54
  55 /*
  56  * the number of bits in a word used to store an hcode (or in an element of
  57  * an array that's used)
  58  */
  59 #define WORDBITS HILBERT_ORDER
  60
  61 #ifdef HILBERT_WANT_ENCODE
  62 /*
  63  * given the coordinates of a point, it finds the sequence number of the point
  64  * on the Hilbert Curve
  65  */
  66 static void
  67 P(encode) (P(t) *dest, P(t) *src)
  68 {
  69   P(t) mask = (P(t))1 << WORDBITS - 1, element, temp1, temp2,
  70     A, W = 0, S, tS, T, tT, J, P = 0, xJ;
  71   uns i = NUMBITS * HILBERT_DIM - HILBERT_DIM, j;
  72
  73   for (j = 0; j < HILBERT_DIM; j++)
  74     dest[j] = 0;
  75   for (j = A = 0; j < HILBERT_DIM; j++)
  76     if (src[j] & mask)
  77       A |= (1 << HILBERT_DIM - 1 - j);
  78
  79   S = tS = A;
  80
  81   P |= S & (1 << HILBERT_DIM - 1);
  82   for (j = 1; j < HILBERT_DIM; j++)
  83     if( S & (1 << HILBERT_DIM - 1 - j) ^ (P >> 1) & (1 << HILBERT_DIM - 1 - j))
  84       P |= (1 << HILBERT_DIM - 1 - j);
  85
  86   /* add in HILBERT_DIM bits to hcode */
  87   element = i / WORDBITS;
  88   if (i % WORDBITS > WORDBITS - HILBERT_DIM)
  89     {
  90       dest[element] |= P << i % WORDBITS;
  91       dest[element + 1] |= P >> WORDBITS - i % WORDBITS;
  92     }
  93   else
  94     dest[element] |= P << i - element * WORDBITS;
  95
  96   J = HILBERT_DIM;
  97   for (j = 1; j < HILBERT_DIM; j++)
  98     if ((P >> j & 1) == (P & 1))
  99       continue;
 100     else
 101       break;
 102   if (j != HILBERT_DIM)
 103     J -= j;
 104   xJ = J - 1;
 105
 106   if (P < 3)
 107     T = 0;
 108   else
 109     if (P % 2)
 110       T = (P - 1) ^ (P - 1) / 2;
 111     else
 112       T = (P - 2) ^ (P - 2) / 2;
 113   tT = T;
 114
 115   for (i -= HILBERT_DIM, mask >>= 1; (int)i >= 0; i -= HILBERT_DIM, mask >>= 1)
 116     {
 117       for (j = A = 0; j < HILBERT_DIM; j++)
 118         if (src[j] & mask)
 119           A |= (1 << HILBERT_DIM - 1 - j);
 120
 121       W ^= tT;
 122       tS = A ^ W;
 123       if (xJ % HILBERT_DIM != 0)
 124         {
 125           temp1 = tS << xJ % HILBERT_DIM;
 126           temp2 = tS >> HILBERT_DIM - xJ % HILBERT_DIM;
 127           S = temp1 | temp2;
 128           S &= ((P(t))1 << HILBERT_DIM) - 1;
 129         }
 130       else
 131         S = tS;
 132
 133       P = S & (1 << HILBERT_DIM - 1);
 134       for (j = 1; j < HILBERT_DIM; j++)
 135         if( S & (1 << HILBERT_DIM - 1 - j) ^ (P >> 1) & (1 << HILBERT_DIM - 1 - j))
 136           P |= (1 << HILBERT_DIM - 1 - j);
 137
 138       /* add in HILBERT_DIM bits to hcode */
 139       element = i / WORDBITS;
 140       if (i % WORDBITS > WORDBITS - HILBERT_DIM)
 141         {
 142           dest[element] |= P << i % WORDBITS;
 143           dest[element + 1] |= P >> WORDBITS - i % WORDBITS;
 144         }
 145       else
 146         dest[element] |= P << i - element * WORDBITS;
 147
 148       if (i > 0)
 149         {
 150           if (P < 3)
 151             T = 0;
 152           else
 153             if (P % 2)
 154               T = (P - 1) ^ (P - 1) / 2;
 155             else
 156               T = (P - 2) ^ (P - 2) / 2;
 157
 158           if (xJ % HILBERT_DIM != 0)
 159             {
 160               temp1 = T >> xJ % HILBERT_DIM;
 161               temp2 = T << HILBERT_DIM - xJ % HILBERT_DIM;
 162               tT = temp1 | temp2;
 163               tT &= ((P(t))1 << HILBERT_DIM) - 1;
 164             }
 165           else
 166             tT = T;
 167
 168           J = HILBERT_DIM;
 169           for (j = 1; j < HILBERT_DIM; j++)
 170             if ((P >> j & 1) == (P & 1))
 171               continue;
 172             else
 173               break;
 174           if (j != HILBERT_DIM)
 175             J -= j;
 176
 177           xJ += J - 1;
 178           /* J %= HILBERT_DIM; */
 179         }
 180     }
 181   for (j = 0; j < HILBERT_DIM; j++)
 182     dest[j] &= ~(P(t))0 >> (8 * sizeof(P(t)) - WORDBITS);
 183 }
 184 #endif
 185
 186 #ifdef HILBERT_WANT_DECODE
 187 /*
 188  * given the sequence number of a point, it finds the coordinates of the point
 189  * on the Hilbert Curve
 190  */
 191 static void
 192 P(decode) (P(t) *dest, P(t) *src)
 193 {
 194   P(t) mask = (P(t))1 << WORDBITS - 1, element, temp1, temp2,
 195     A, W = 0, S, tS, T, tT, J, P = 0, xJ;
 196   uns i = NUMBITS * HILBERT_DIM - HILBERT_DIM, j;
 197
 198   for (j = 0; j < HILBERT_DIM; j++)
 199     dest[j] = 0;
 200
 201   /*--- P ---*/
 202   element = i / WORDBITS;
 203   P = src[element];
 204   if (i % WORDBITS > WORDBITS - HILBERT_DIM)
 205     {
 206       temp1 = src[element + 1];
 207       P >>= i % WORDBITS;
 208       temp1 <<= WORDBITS - i % WORDBITS;
 209       P |= temp1;
 210     }
 211   else
 212     P >>= i % WORDBITS; /* P is a HILBERT_DIM bit hcode */
 213
 214   /* the & masks out spurious highbit values */
 215   if (HILBERT_DIM < WORDBITS)
 216     P &= (1 << HILBERT_DIM) -1;
 217
 218   /*--- xJ ---*/
 219   J = HILBERT_DIM;
 220   for (j = 1; j < HILBERT_DIM; j++)
 221     if ((P >> j & 1) == (P & 1))
 222       continue;
 223     else
 224       break;
 225   if (j != HILBERT_DIM)
 226     J -= j;
 227   xJ = J - 1;
 228
 229   /*--- S, tS, A ---*/
 230   A = S = tS = P ^ P / 2;
 231
 232
 233   /*--- T ---*/
 234   if (P < 3)
 235     T = 0;
 236   else
 237     if (P % 2)
 238       T = (P - 1) ^ (P - 1) / 2;
 239     else
 240       T = (P - 2) ^ (P - 2) / 2;
 241
 242   /*--- tT ---*/
 243   tT = T;
 244
 245   /*--- distrib bits to coords ---*/
 246   for (j = HILBERT_DIM - 1; P > 0; P >>=1, j--)
 247     if (P & 1)
 248       dest[j] |= mask;
 249
 250
 251   for (i -= HILBERT_DIM, mask >>= 1; (int)i >= 0; i -= HILBERT_DIM, mask >>= 1)
 252     {
 253       /*--- P ---*/
 254       element = i / WORDBITS;
 255       P = src[element];
 256       if (i % WORDBITS > WORDBITS - HILBERT_DIM)
 257         {
 258           temp1 = src[element + 1];
 259           P >>= i % WORDBITS;
 260           temp1 <<= WORDBITS - i % WORDBITS;
 261           P |= temp1;
 262         }
 263       else
 264         P >>= i % WORDBITS;     /* P is a HILBERT_DIM bit hcode */
 265
 266       /* the & masks out spurious highbit values */
 267       if (HILBERT_DIM < WORDBITS)
 268         P &= (1 << HILBERT_DIM) -1;
 269
 270       /*--- S ---*/
 271       S = P ^ P / 2;
 272
 273       /*--- tS ---*/
 274       if (xJ % HILBERT_DIM != 0)
 275         {
 276           temp1 = S >> xJ % HILBERT_DIM;
 277           temp2 = S << HILBERT_DIM - xJ % HILBERT_DIM;
 278           tS = temp1 | temp2;
 279           tS &= ((P(t))1 << HILBERT_DIM) - 1;
 280         }
 281       else
 282         tS = S;
 283
 284       /*--- W ---*/
 285       W ^= tT;
 286
 287       /*--- A ---*/
 288       A = W ^ tS;
 289
 290       /*--- distrib bits to coords ---*/
 291       for (j = HILBERT_DIM - 1; A > 0; A >>=1, j--)
 292         if (A & 1)
 293           dest[j] |= mask;
 294
 295       if (i > 0)
 296         {
 297           /*--- T ---*/
 298           if (P < 3)
 299             T = 0;
 300           else
 301             if (P % 2)
 302               T = (P - 1) ^ (P - 1) / 2;
 303             else
 304               T = (P - 2) ^ (P - 2) / 2;
 305
 306           /*--- tT ---*/
 307           if (xJ % HILBERT_DIM != 0)
 308             {
 309               temp1 = T >> xJ % HILBERT_DIM;
 310               temp2 = T << HILBERT_DIM - xJ % HILBERT_DIM;
 311               tT = temp1 | temp2;
 312               tT &= ((P(t))1 << HILBERT_DIM) - 1;
 313             }
 314           else
 315             tT = T;
 316
 317           /*--- xJ ---*/
 318           J = HILBERT_DIM;
 319           for (j = 1; j < HILBERT_DIM; j++)
 320             if ((P >> j & 1) == (P & 1))
 321               continue;
 322             else
 323               break;
 324           if (j != HILBERT_DIM)
 325             J -= j;
 326           xJ += J - 1;
 327         }
 328     }
 329 }
 330 #endif
 331
 332 #undef P
 333 #undef HILBERT_PREFIX
 334 #undef HILBERT_DIM
 335 #undef HILBERT_TYPE
 336 #undef HILBERT_ORDER
 337 #undef HILBERT_WANT_DECODE
 338 #undef HILBERT_WANT_ENCODE
 339 #undef NUMBITS
 340 #undef WORDBITS