just playing with multidimensional Hilbert curves

(trying to optimize duplicates search)

diff --git a/images/Makefile b/images/Makefile
index b35ab2ecaac92dba146f29b09cee87b7b80e4125..ccc343a323d481a64d240ae79e2cdc00211751ec 100644 (file)
--- a/images/Makefile
+++ b/images/Makefile
@@ -12,10 +12,12 @@ $(o)/images/image-idx: LIBS+=-lGraphicsMagick -ljpeg -lpng
 $(o)/images/image-test: $(o)/images/image-test.o $(o)/images/kd-tree.o $(LIBSH)
 
 $(o)/images/color-t: LIBS+=-lm
+$(o)/images/hilbert-test: LIBS+= $(LIBSH)
 
-TESTS+=$(addprefix $(o)/images/,color.test)
+TESTS+=$(addprefix $(o)/images/,color.test hilbert-test.test)
 
 $(o)/images/color.test: $(o)/images/color-t
+$(o)/images/hilbert-test.test: $(o)/images/hilbert-test
 
 # By :;DF
 $(o)/images/block_info.o $(o)/images/block_info.oo: CFLAGS+=-I/usr/include/GraphicsMagick

diff --git a/images/hilbert-test.c b/images/hilbert-test.c
new file mode 100644 (file)
index 0000000..235c3ae
--- /dev/null
+++ b/images/hilbert-test.c
@@ -0,0 +1,59 @@
+/* Tests for multidimensional Hilbert curves */
+
+#define LOCAL_DEBUG
+
+#include "lib/lib.h"
+
+#include <stdlib.h>
+
+static uns test1_dim;
+static uns test1_order;
+#define HILBERT_PREFIX(x) test1_##x
+#define HILBERT_DIM test1_dim
+#define HILBERT_ORDER test1_order
+#define HILBERT_WANT_DECODE
+#define HILBERT_WANT_ENCODE
+#include "images/hilbert.h"
+
+static void
+test1(void)
+{
+  uns a[32], b[32], c[32];
+  for (test1_dim = 2; test1_dim <= 8; test1_dim++)
+    for (test1_order = 8; test1_order <= 32; test1_order++)
+      for (uns i = 0; i < 1000; i++)
+        {
+          for (uns j = 0; j < test1_dim; j++)
+           a[j] = (uns)rand() >> (32 - test1_order);
+          test1_encode(b, a);
+          test1_decode(c, b);
+          for (uns j = 0; j < test1_dim; j++)
+           if (a[j] != c[j])
+             die("Error... dim=%d order=%d testnum=%d index=%d val1=0x%08x val2=0x%08x", test1_dim, test1_order, i, j, a[j], c[j]);
+        }
+}
+
+#if 0
+#include "images/hilbert-origin.h"
+static void
+test_origin(void)
+{
+  Hcode code;
+  Point pt, pt2;
+  pt.hcode[0] = 0x12345678;
+  pt.hcode[1] = 0x654321;
+  pt.hcode[2] = 0x11122233;
+  code = H_encode(pt);
+  pt2 = H_decode(code);
+  DBG("origin: [%08x, %08x, %08x] --> [%08x, %08x %08x] --> [%08x, %08x %08x]", 
+    pt.hcode[0], pt.hcode[1], pt.hcode[2], code.hcode[0], code.hcode[1], code.hcode[2], pt2.hcode[0], pt2.hcode[1], pt2.hcode[2]);
+}
+#endif
+
+int
+main(int argc UNUSED, char **argv UNUSED)
+{
+  test1();
+  //test_origin();
+  return 0;
+}

diff --git a/images/hilbert-test.t b/images/hilbert-test.t
new file mode 100644 (file)
index 0000000..f8794e3
--- /dev/null
+++ b/images/hilbert-test.t
@@ -0,0 +1,3 @@
+# Tests for multidimensional Hilbert curves
+
+Run:   obj/images/hilbert-test

diff --git a/images/hilbert.h b/images/hilbert.h
new file mode 100644 (file)
index 0000000..4971826
--- /dev/null
+++ b/images/hilbert.h
@@ -0,0 +1,340 @@
+/*
+ *     Image Library -- multidimensional Hilbert curves
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ *
+ *
+ *     References:
+ *     - http://www.dcs.bbk.ac.uk/~jkl/mapping.c
+ *       (c) 2002 J.K.Lawder
+ *     - J.K. Lawder. Calculation of Mappings between One and n-dimensional Values
+ *       Using the Hilbert Space-Filling Curve. Technical Report JL1/00, Birkbeck
+ *       College, University of London, 2000.
+ *
+ *     FIXME:
+ *     - the algorithm fails for some combinations of HILBERT_DIM and HILBERT_ORDER,
+ *       but it should be safe for HILBERT_DIM = 2..8, HILBERT_ORDER = 8..32
+ *     - clean and optimize the code
+ */
+
+#ifndef HILBERT_PREFIX
+#  error Undefined HILBERT_PREFIX
+#endif
+
+#define P(x) HILBERT_PREFIX(x)
+
+/*
+ * HILBERT_DIM is the number of dimensions in space through which the
+ * Hilbert Curve passes.
+ * Don't use this implementation with values for HILBERT_DIM of > 31!
+ * Also, make sure you use a 32 bit compiler!
+ */
+#ifndef HILBERT_DIM
+#  define HILBERT_DIM 2
+#endif
+
+#ifndef HILBERT_TYPE
+#  define HILBERT_TYPE u32
+#endif
+
+#ifndef HILBERT_ORDER
+#  define HILBERT_ORDER (8 * sizeof(HILBERT_TYPE))
+#endif
+
+typedef HILBERT_TYPE P(t);
+
+/*
+ * retained for historical reasons: the number of bits in an attribute value:
+ * effectively the order of a curve
+ */
+#define NUMBITS HILBERT_ORDER
+
+/*
+ * the number of bits in a word used to store an hcode (or in an element of
+ * an array that's used)
+ */
+#define WORDBITS HILBERT_ORDER
+
+#ifdef HILBERT_WANT_ENCODE
+/*
+ * given the coordinates of a point, it finds the sequence number of the point
+ * on the Hilbert Curve
+ */
+static void
+P(encode) (P(t) *dest, P(t) *src)
+{
+  P(t) mask = (P(t))1 << WORDBITS - 1, element, temp1, temp2,
+    A, W = 0, S, tS, T, tT, J, P = 0, xJ;
+  uns i = NUMBITS * HILBERT_DIM - HILBERT_DIM, j;
+
+  for (j = 0; j < HILBERT_DIM; j++)
+    dest[j] = 0;
+  for (j = A = 0; j < HILBERT_DIM; j++)
+    if (src[j] & mask)
+      A |= (1 << HILBERT_DIM - 1 - j);
+
+  S = tS = A;
+
+  P |= S & (1 << HILBERT_DIM - 1);
+  for (j = 1; j < HILBERT_DIM; j++)
+    if( S & (1 << HILBERT_DIM - 1 - j) ^ (P >> 1) & (1 << HILBERT_DIM - 1 - j))
+      P |= (1 << HILBERT_DIM - 1 - j);
+
+  /* add in HILBERT_DIM bits to hcode */
+  element = i / WORDBITS;
+  if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+    {
+      dest[element] |= P << i % WORDBITS;
+      dest[element + 1] |= P >> WORDBITS - i % WORDBITS;
+    }
+  else
+    dest[element] |= P << i - element * WORDBITS;
+
+  J = HILBERT_DIM;
+  for (j = 1; j < HILBERT_DIM; j++)
+    if ((P >> j & 1) == (P & 1))
+      continue;
+    else
+      break;
+  if (j != HILBERT_DIM)
+    J -= j;
+  xJ = J - 1;
+
+  if (P < 3)
+    T = 0;
+  else
+    if (P % 2)
+      T = (P - 1) ^ (P - 1) / 2;
+    else
+      T = (P - 2) ^ (P - 2) / 2;
+  tT = T;
+
+  for (i -= HILBERT_DIM, mask >>= 1; (int)i >= 0; i -= HILBERT_DIM, mask >>= 1)
+    {
+      for (j = A = 0; j < HILBERT_DIM; j++)
+        if (src[j] & mask)
+          A |= (1 << HILBERT_DIM - 1 - j);
+
+      W ^= tT;
+      tS = A ^ W;
+      if (xJ % HILBERT_DIM != 0)
+        {
+          temp1 = tS << xJ % HILBERT_DIM;
+          temp2 = tS >> HILBERT_DIM - xJ % HILBERT_DIM;
+         S = temp1 | temp2;
+         S &= ((P(t))1 << HILBERT_DIM) - 1;
+       }
+      else
+       S = tS;
+
+      P = S & (1 << HILBERT_DIM - 1);
+      for (j = 1; j < HILBERT_DIM; j++)
+       if( S & (1 << HILBERT_DIM - 1 - j) ^ (P >> 1) & (1 << HILBERT_DIM - 1 - j))
+         P |= (1 << HILBERT_DIM - 1 - j);
+
+      /* add in HILBERT_DIM bits to hcode */
+      element = i / WORDBITS;
+      if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+       {
+         dest[element] |= P << i % WORDBITS;
+         dest[element + 1] |= P >> WORDBITS - i % WORDBITS;
+       }
+      else
+       dest[element] |= P << i - element * WORDBITS;
+
+      if (i > 0)
+       {
+         if (P < 3)
+           T = 0;
+         else
+           if (P % 2)
+             T = (P - 1) ^ (P - 1) / 2;
+           else
+             T = (P - 2) ^ (P - 2) / 2;
+
+         if (xJ % HILBERT_DIM != 0)
+           {
+             temp1 = T >> xJ % HILBERT_DIM;
+             temp2 = T << HILBERT_DIM - xJ % HILBERT_DIM;
+             tT = temp1 | temp2;
+             tT &= ((P(t))1 << HILBERT_DIM) - 1;
+           }
+         else
+           tT = T;
+
+         J = HILBERT_DIM;
+         for (j = 1; j < HILBERT_DIM; j++)
+           if ((P >> j & 1) == (P & 1))
+             continue;
+           else
+             break;
+         if (j != HILBERT_DIM)
+           J -= j;
+
+         xJ += J - 1;
+          /* J %= HILBERT_DIM; */
+       }
+    }
+  for (j = 0; j < HILBERT_DIM; j++)
+    dest[j] &= ~(P(t))0 >> (8 * sizeof(P(t)) - WORDBITS);
+}
+#endif
+
+#ifdef HILBERT_WANT_DECODE
+/*
+ * given the sequence number of a point, it finds the coordinates of the point
+ * on the Hilbert Curve
+ */
+static void
+P(decode) (P(t) *dest, P(t) *src)
+{
+  P(t) mask = (P(t))1 << WORDBITS - 1, element, temp1, temp2,
+    A, W = 0, S, tS, T, tT, J, P = 0, xJ;
+  uns i = NUMBITS * HILBERT_DIM - HILBERT_DIM, j;
+
+  for (j = 0; j < HILBERT_DIM; j++)
+    dest[j] = 0;
+
+  /*--- P ---*/
+  element = i / WORDBITS;
+  P = src[element];
+  if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+    {
+      temp1 = src[element + 1];
+      P >>= i % WORDBITS;
+      temp1 <<= WORDBITS - i % WORDBITS;
+      P |= temp1;
+    }
+  else
+    P >>= i % WORDBITS;        /* P is a HILBERT_DIM bit hcode */
+
+  /* the & masks out spurious highbit values */
+  if (HILBERT_DIM < WORDBITS)
+    P &= (1 << HILBERT_DIM) -1;
+
+  /*--- xJ ---*/
+  J = HILBERT_DIM;
+  for (j = 1; j < HILBERT_DIM; j++)
+    if ((P >> j & 1) == (P & 1))
+      continue;
+    else
+      break;
+  if (j != HILBERT_DIM)
+    J -= j;
+  xJ = J - 1;
+
+  /*--- S, tS, A ---*/
+  A = S = tS = P ^ P / 2;
+
+
+  /*--- T ---*/
+  if (P < 3)
+    T = 0;
+  else
+    if (P % 2)
+      T = (P - 1) ^ (P - 1) / 2;
+    else
+      T = (P - 2) ^ (P - 2) / 2;
+
+  /*--- tT ---*/
+  tT = T;
+
+  /*--- distrib bits to coords ---*/
+  for (j = HILBERT_DIM - 1; P > 0; P >>=1, j--)
+    if (P & 1)
+      dest[j] |= mask;
+
+
+  for (i -= HILBERT_DIM, mask >>= 1; (int)i >= 0; i -= HILBERT_DIM, mask >>= 1)
+    {
+      /*--- P ---*/
+      element = i / WORDBITS;
+      P = src[element];
+      if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+       {
+         temp1 = src[element + 1];
+         P >>= i % WORDBITS;
+         temp1 <<= WORDBITS - i % WORDBITS;
+         P |= temp1;
+       }
+      else
+       P >>= i % WORDBITS;     /* P is a HILBERT_DIM bit hcode */
+
+      /* the & masks out spurious highbit values */
+      if (HILBERT_DIM < WORDBITS)
+        P &= (1 << HILBERT_DIM) -1;
+
+      /*--- S ---*/
+      S = P ^ P / 2;
+
+      /*--- tS ---*/
+      if (xJ % HILBERT_DIM != 0)
+       {
+         temp1 = S >> xJ % HILBERT_DIM;
+         temp2 = S << HILBERT_DIM - xJ % HILBERT_DIM;
+         tS = temp1 | temp2;
+         tS &= ((P(t))1 << HILBERT_DIM) - 1;
+       }
+      else
+       tS = S;
+
+      /*--- W ---*/
+      W ^= tT;
+
+      /*--- A ---*/
+      A = W ^ tS;
+
+      /*--- distrib bits to coords ---*/
+      for (j = HILBERT_DIM - 1; A > 0; A >>=1, j--)
+       if (A & 1)
+         dest[j] |= mask;
+
+      if (i > 0)
+       {
+         /*--- T ---*/
+         if (P < 3)
+           T = 0;
+         else
+           if (P % 2)
+             T = (P - 1) ^ (P - 1) / 2;
+           else
+             T = (P - 2) ^ (P - 2) / 2;
+
+         /*--- tT ---*/
+         if (xJ % HILBERT_DIM != 0)
+           {
+             temp1 = T >> xJ % HILBERT_DIM;
+             temp2 = T << HILBERT_DIM - xJ % HILBERT_DIM;
+             tT = temp1 | temp2;
+             tT &= ((P(t))1 << HILBERT_DIM) - 1;
+           }
+         else
+           tT = T;
+
+         /*--- xJ ---*/
+         J = HILBERT_DIM;
+         for (j = 1; j < HILBERT_DIM; j++)
+           if ((P >> j & 1) == (P & 1))
+             continue;
+           else
+             break;
+         if (j != HILBERT_DIM)
+           J -= j;
+         xJ += J - 1;
+       }
+    }
+}
+#endif
+
+#undef P
+#undef HILBERT_PREFIX
+#undef HILBERT_DIM
+#undef HILBERT_TYPE
+#undef HILBERT_ORDER
+#undef HILBERT_WANT_DECODE
+#undef HILBERT_WANT_ENCODE
+#undef NUMBITS
+#undef WORDBITS