From 19fee44ed7e0d780267e0a049ec2eb82db00d46c Mon Sep 17 00:00:00 2001 From: Pavel Charvat Date: Wed, 3 May 2006 14:18:37 +0200 Subject: [PATCH] just playing with multidimensional Hilbert curves (trying to optimize duplicates search) --- images/Makefile | 4 +- images/hilbert-test.c | 59 ++++++++ images/hilbert-test.t | 3 + images/hilbert.h | 340 ++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 405 insertions(+), 1 deletion(-) create mode 100644 images/hilbert-test.c create mode 100644 images/hilbert-test.t create mode 100644 images/hilbert.h diff --git a/images/Makefile b/images/Makefile index b35ab2ec..ccc343a3 100644 --- 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 index 00000000..235c3ae4 --- /dev/null +++ b/images/hilbert-test.c @@ -0,0 +1,59 @@ +/* Tests for multidimensional Hilbert curves */ + +#define LOCAL_DEBUG + +#include "lib/lib.h" + +#include + +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 index 00000000..f8794e3c --- /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 index 00000000..4971826a --- /dev/null +++ b/images/hilbert.h @@ -0,0 +1,340 @@ +/* + * Image Library -- multidimensional Hilbert curves + * + * (c) 2006 Pavel Charvat + * + * 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 -- 2.39.2