diff options
| -rw-r--r-- | Makefile | 4 | ||||
| -rw-r--r-- | generate.c | 5 | ||||
| -rw-r--r-- | hilbert.c | 170 | ||||
| -rw-r--r-- | hilbert.h | 31 | ||||
| -rw-r--r-- | options.c | 4 | ||||
| -rw-r--r-- | options.h | 1 |
6 files changed, 212 insertions, 3 deletions
@@ -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) @@ -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 @@ -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; @@ -20,6 +20,7 @@ typedef enum { MODE_HUE_SORT, MODE_RANDOM, MODE_MORTON, + MODE_HILBERT, } mode_t; // Possible pixel selection modes |