summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorTavian Barnes <tavianator@tavianator.com>2015-10-04 13:36:25 -0400
committerTavian Barnes <tavianator@tavianator.com>2015-10-04 13:36:25 -0400
commitf0bc8b8ddbd1ef1041f1249d510bb379903777d5 (patch)
tree4b11e28a59986e5d0e5d83326c1e6dabbefc6a04
parentc6ae17f774721875a2833b1138bad9f516a72003 (diff)
downloadkd-forest-f0bc8b8ddbd1ef1041f1249d510bb379903777d5.tar.xz
Add support for Hilbert order.
-rw-r--r--Makefile4
-rw-r--r--generate.c5
-rw-r--r--hilbert.c170
-rw-r--r--hilbert.h31
-rw-r--r--options.c4
-rw-r--r--options.h1
6 files changed, 212 insertions, 3 deletions
diff --git a/Makefile b/Makefile
index 6041235..4e035e0 100644
--- a/Makefile
+++ b/Makefile
@@ -15,12 +15,12 @@ LDFLAGS = -Wl,-O1,--sort-common,--as-needed,-z,relro
LIBS = -lm -lpng
RM = rm -f
-DEPS = Makefile color.h generate.h kd-forest.h options.h util.h
+DEPS = Makefile color.h hilbert.h generate.h kd-forest.h options.h util.h
IMAGEFLAGS = -b 24 -s -l min -c Lab
ANIMFLAGS = -b 19 -s -l mean -c Lab
-kd-forest: color.o generate.o kd-forest.o main.o options.o util.o
+kd-forest: color.o generate.o hilbert.o kd-forest.o main.o options.o util.o
$(CC) $(CFLAGS) $(LDFLAGS) $^ $(LIBS) -o kd-forest
%.o: %.c $(DEPS)
diff --git a/generate.c b/generate.c
index 6b7e925..9a4eac8 100644
--- a/generate.c
+++ b/generate.c
@@ -11,6 +11,7 @@
#include "generate.h"
#include "color.h"
+#include "hilbert.h"
#include "util.h"
#include <stdlib.h>
@@ -38,6 +39,10 @@ generate_colors(const options_t *options)
}
break;
+ case MODE_HILBERT:
+ hilbert_point(3, grb_bits, n, grb);
+ break;
+
default:
for (unsigned int j = 0; j < 3; ++j) {
grb[j] = n & ((1 << grb_bits[j]) - 1);
diff --git a/hilbert.c b/hilbert.c
new file mode 100644
index 0000000..98cf247
--- /dev/null
+++ b/hilbert.c
@@ -0,0 +1,170 @@
+/*********************************************************************
+ * kd-forest *
+ * Copyright (C) 2015 Tavian Barnes <tavianator@tavianator.com> *
+ * *
+ * This program is free software. It comes without any warranty, to *
+ * the extent permitted by applicable law. You can redistribute it *
+ * and/or modify it under the terms of the Do What The Fuck You Want *
+ * To Public License, Version 2, as published by Sam Hocevar. See *
+ * the COPYING file or http://www.wtfpl.net/ for more details. *
+ *********************************************************************/
+
+#include "hilbert.h"
+#include <stdint.h>
+
+// These algorithms are described in "Compact Hilbert Indices" by Chris Hamilton
+
+// Right rotation of x by b bits out of n
+static uint32_t
+ror(uint32_t x, unsigned int b, unsigned int n)
+{
+ uint32_t l = x & ((1 << b) - 1);
+ uint32_t r = x >> b;
+ return (l << (n - b)) | r;
+}
+
+// Left rotation of x by b bits out of n
+static uint32_t
+rol(uint32_t x, unsigned int b, unsigned int n)
+{
+ return ror(x, n - b, n);
+}
+
+// Binary reflected Gray code
+uint32_t
+gray_code(uint32_t i)
+{
+ return i ^ (i >> 1);
+}
+
+// e(i), the entry point for the ith sub-hypercube
+uint32_t
+entry_point(uint32_t i)
+{
+ if (i == 0) {
+ return 0;
+ } else {
+ return gray_code((i - 1) & ~1U);
+ }
+}
+
+// g(i), the inter sub-hypercube direction
+unsigned int
+inter_direction(uint32_t i)
+{
+ // g(i) counts the trailing set bits in i
+ unsigned int g = 0;
+ while (i & 1) {
+ ++g;
+ i >>= 1;
+ }
+ return g;
+}
+
+// d(i), the intra sub-hypercube direction
+unsigned int
+intra_direction(uint32_t i)
+{
+ if (i & 1) {
+ return inter_direction(i);
+ } else if (i > 0) {
+ return inter_direction(i - 1);
+ } else {
+ return 0;
+ }
+}
+
+// T transformation inverse
+uint32_t
+t_inverse(unsigned int dimensions, uint32_t e, unsigned int d, uint32_t a)
+{
+ return rol(a, d, dimensions) ^ e;
+}
+
+// GrayCodeRankInverse
+void
+gray_code_rank_inverse(unsigned int dimensions, uint32_t mu, uint32_t pi, unsigned int r, unsigned int free_bits, uint32_t *i, uint32_t *g)
+{
+ // *i is the inverse rank of r
+ // *g is gray_code(i)
+
+ *i = 0;
+ *g = 0;
+
+ for (unsigned int j = free_bits - 1, k = dimensions; k-- > 0;) {
+ if (mu & (1 << k)) {
+ *i |= ((r >> j) & 1) << k;
+ *g |= (*i ^ (*i >> 1)) & (1 << k);
+ --j;
+ } else {
+ *g |= pi & (1 << k);
+ *i |= (*g ^ (*i >> 1)) & (1 << k);
+ }
+ }
+}
+
+// ExtractMask
+void
+extract_mask(unsigned int dimensions, const unsigned int extents[], unsigned int i, uint32_t *mu, unsigned int *free_bits)
+{
+ // *mu is the mask
+ // *free_bits is popcount(*mu)
+
+ *mu = 0;
+ *free_bits = 0;
+
+ for (unsigned int j = dimensions; j-- > 0;) {
+ *mu <<= 1;
+
+ if (extents[j] > i) {
+ *mu |= 1;
+ *free_bits += 1;
+ }
+ }
+}
+
+// CompactHilbertIndexInverse
+void
+hilbert_point(unsigned int dimensions, const unsigned int extents[], uint32_t index, uint32_t point[])
+{
+ unsigned int max_extent = 0, total_extent = 0;
+ for (unsigned int i = 0; i < dimensions; ++i) {
+ if (extents[i] > max_extent) {
+ max_extent = extents[i];
+ }
+ total_extent += extents[i];
+ point[i] = 0;
+ }
+
+ uint32_t e = 0;
+ unsigned int k = 0;
+
+ // Next direction; we use d instead of d + 1 everywhere
+ unsigned int d = 1;
+
+ for (unsigned int i = max_extent; i-- > 0;) {
+ uint32_t mu;
+ unsigned int free_bits;
+ extract_mask(dimensions, extents, i, &mu, &free_bits);
+ mu = ror(mu, d, dimensions);
+
+ uint32_t pi = ror(e, d, dimensions) & ~mu;
+
+ unsigned int r = (index >> (total_extent - k - free_bits)) & ((1 << free_bits) - 1);
+
+ k += free_bits;
+
+ uint32_t w, l;
+ gray_code_rank_inverse(dimensions, mu, pi, r, free_bits, &w, &l);
+
+ l = t_inverse(dimensions, e, d, l);
+
+ for (unsigned int j = 0; j < 3; ++j) {
+ point[j] |= (l & 1) << i;
+ l >>= 1;
+ }
+
+ e = e ^ ror(entry_point(w), d, dimensions);
+ d = (d + intra_direction(w) + 1)%dimensions;
+ }
+}
diff --git a/hilbert.h b/hilbert.h
new file mode 100644
index 0000000..4628ad6
--- /dev/null
+++ b/hilbert.h
@@ -0,0 +1,31 @@
+/*********************************************************************
+ * kd-forest *
+ * Copyright (C) 2015 Tavian Barnes <tavianator@tavianator.com> *
+ * *
+ * This program is free software. It comes without any warranty, to *
+ * the extent permitted by applicable law. You can redistribute it *
+ * and/or modify it under the terms of the Do What The Fuck You Want *
+ * To Public License, Version 2, as published by Sam Hocevar. See *
+ * the COPYING file or http://www.wtfpl.net/ for more details. *
+ *********************************************************************/
+
+#ifndef HILBERT_H
+#define HILBERT_H
+
+#include <stdint.h>
+
+/**
+ * Compute the point corresponding to the given (compact) Hilbert index.
+ *
+ * @param dimensions
+ * The number of spatial dimensions.
+ * @param extents
+ * The bit depth of each dimension.
+ * @param index
+ * The (compact) Hilbert index of the desired point.
+ * @param[out] point
+ * Will hold the point on the Hilbert curve at index.
+ */
+void hilbert_point(unsigned int dimensions, const unsigned int extents[], uint32_t index, uint32_t point[]);
+
+#endif // HILBERT_H
diff --git a/options.c b/options.c
index db48615..5d7c325 100644
--- a/options.c
+++ b/options.c
@@ -225,7 +225,7 @@ print_usage(FILE *file, const char *command, bool verbose)
#define usage(...) print_colorized(file, tty, __VA_ARGS__)
usage("Usage:\n");
usage(" !$! *%s* [-b|--bit-depth @DEPTH@]\n", command);
- usage(" %s [-s|--hue-sort] [-r|--random] [-M|--morton]\n", whitespace);
+ usage(" %s [-s|--hue-sort] [-r|--random] [-M|--morton] [-H|--hilbert]\n", whitespace);
usage(" %s [-l|--selection @min@|@mean@]\n", whitespace);
usage(" %s [-c|--color-space @RGB@|@Lab@|@Luv@]\n", whitespace);
usage(" %s [-w|--width @WIDTH@] [-h|--height @HEIGHT@]\n", whitespace);
@@ -305,6 +305,8 @@ parse_options(options_t *options, int argc, char *argv[])
options->mode = MODE_RANDOM;
} else if (parse_arg(argc, argv, "-M", "--morton", NULL, &i, &error)) {
options->mode = MODE_MORTON;
+ } else if (parse_arg(argc, argv, "-H", "--hilbert", NULL, &i, &error)) {
+ options->mode = MODE_HILBERT;
} else if (parse_arg(argc, argv, "-l", "--selection", &value, &i, &error)) {
if (strcmp(value, "min") == 0) {
options->selection = SELECTION_MIN;
diff --git a/options.h b/options.h
index 06ad69b..62038a4 100644
--- a/options.h
+++ b/options.h
@@ -20,6 +20,7 @@ typedef enum {
MODE_HUE_SORT,
MODE_RANDOM,
MODE_MORTON,
+ MODE_HILBERT,
} mode_t;
// Possible pixel selection modes