diff options
| -rw-r--r-- | .gitignore | 2 | ||||
| -rw-r--r-- | COPYING | 13 | ||||
| -rw-r--r-- | Makefile | 26 | ||||
| -rw-r--r-- | kd-forest.c | 302 | ||||
| -rw-r--r-- | kd-forest.h | 57 | ||||
| -rw-r--r-- | main.c | 281 | ||||
| -rw-r--r-- | util.c | 33 | ||||
| -rw-r--r-- | util.h | 20 |
8 files changed, 734 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f815932 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +/kd-forest +*.png @@ -0,0 +1,13 @@ + DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE + Version 2, December 2004 + + Copyright (C) 2004 Sam Hocevar <sam@hocevar.net> + + Everyone is permitted to copy and distribute verbatim or modified + copies of this license document, and changing it is allowed as long + as the name is changed. + + DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. You just DO WHAT THE FUCK YOU WANT TO. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..77bc5c9 --- /dev/null +++ b/Makefile @@ -0,0 +1,26 @@ +##################################################################### +# kd-forest # +# Copyright (C) 2014 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. # +##################################################################### + +CC ?= gcc +CFLAGS ?= -std=c99 -pipe -O2 -Werror -Wall -Wpedantic -Wextra -Wno-sign-compare -Wno-unused-parameter -Wunreachable-code -Wshadow -Wpointer-arith -Wwrite-strings -Wcast-align -Wstrict-prototypes +LDFLAGS ?= -lm -lpng +RM ?= rm -f + +kd-forest: kd-forest.c kd-forest.h util.c util.h main.c + $(CC) $(CFLAGS) $(LDFLAGS) -o kd-forest kd-forest.c util.c main.c + +image: kd-forest + ./kd-forest + +clean: + $(RM) kd-forest + +.PHONY: image clean diff --git a/kd-forest.c b/kd-forest.c new file mode 100644 index 0000000..4a5bf37 --- /dev/null +++ b/kd-forest.c @@ -0,0 +1,302 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 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 "kd-forest.h" +#include "util.h" +#include <stdbool.h> +#include <stdlib.h> +#include <string.h> + +void +kd_node_init(kd_node_t *node, unsigned int x, unsigned int y) +{ + node->left = node->right = NULL; + node->x = x; + node->y = y; + node->added = node->removed = false; +} + +void +kd_node_set_color(kd_node_t *node, uint32_t color) +{ + node->coords[0] = (color >> 16) & 0xFF; + node->coords[1] = (color >> 8) & 0xFF; + node->coords[2] = color & 0xFF; +} + +static size_t +kd_collect_nodes(kd_node_t *root, kd_node_t **buffer, bool include_removed) +{ + size_t count = 0; + if (include_removed || !root->removed) { + buffer[0] = root; + ++count; + } + if (root->left) { + count += kd_collect_nodes(root->left, buffer + count, include_removed); + } + if (root->right) { + count += kd_collect_nodes(root->right, buffer + count, include_removed); + } + return count; +} + +typedef int kd_comparator(const void *a, const void* b); + +static int kd_compare0(const void *a, const void *b) { + int aval = (*(const kd_node_t **)a)->coords[0]; + int bval = (*(const kd_node_t **)b)->coords[0]; + return aval - bval; +} + +static int kd_compare1(const void *a, const void *b) { + int aval = (*(const kd_node_t **)a)->coords[1]; + int bval = (*(const kd_node_t **)b)->coords[1]; + return aval - bval; +} + +static int kd_compare2(const void *a, const void *b) { + int aval = (*(const kd_node_t **)a)->coords[2]; + int bval = (*(const kd_node_t **)b)->coords[2]; + return aval - bval; +} + +static kd_comparator *kd_comparators[KD_DIMEN] = { + kd_compare0, + kd_compare1, + kd_compare2, +}; + +// When building k-d trees, we use KD_DIMEN sorted arrays of nodes plus one +// extra array for scratch space +#define KD_BUFSIZE (KD_DIMEN + 1) + +static kd_node_t * +kd_build_tree_recursive(kd_node_t **buffers[KD_BUFSIZE], size_t size, + unsigned int coord) +{ + if (size == 0) { + return NULL; + } + + size_t split = size/2; + size_t left_size = split, right_size = size - left_size - 1; + kd_node_t *root = buffers[coord][split]; + for (size_t i = 0; i < size; ++i) { + buffers[coord][i]->is_left = i < left_size; + } + + kd_node_t **right_buffers[KD_BUFSIZE]; + for (int i = 0; i < KD_DIMEN; ++i) { + right_buffers[i] = buffers[i] + left_size + 1; + } + + kd_node_t **scratch = buffers[KD_DIMEN]; + right_buffers[KD_DIMEN] = scratch; + + for (size_t i = 0; i < KD_DIMEN; ++i) { + if (i == coord) { + continue; + } + + kd_node_t **buffer = buffers[i]; + kd_node_t **right_buffer = right_buffers[i]; + for (size_t j = 0, k = 0, skip = 0; j < size; ++j) { + if (buffer[j]->is_left) { + buffer[j - skip] = buffer[j]; + } else { + if (buffer[j] != root) { + scratch[k] = buffer[j]; + ++k; + } + ++skip; + } + } + for (size_t j = 0; j < right_size; ++j) { + right_buffer[j] = scratch[j]; + } + } + + coord = (coord + 1)%KD_DIMEN; + root->left = kd_build_tree_recursive(buffers, left_size, coord); + root->right = kd_build_tree_recursive(right_buffers, right_size, coord); + + return root; +} + +static kd_node_t * +kd_build_tree(kd_node_t **buffers[KD_BUFSIZE], size_t size) +{ + for (int i = 1; i < KD_DIMEN; ++i) { + memcpy(buffers[i], buffers[0], size*sizeof(kd_node_t *)); + } + for (int i = 0; i < KD_DIMEN; ++i) { + qsort(buffers[i], size, sizeof(kd_node_t *), kd_comparators[i]); + } + return kd_build_tree_recursive(buffers, size, 0); +} + +static int +kd_distance_sq(kd_node_t *a, kd_node_t *b) +{ + int result = 0; + for (int i = 0; i < KD_DIMEN; ++i) { + int d = a->coords[i] - b->coords[i]; + result += d*d; + } + return result; +} + +static void +kd_find_nearest_recursive(kd_node_t *root, kd_node_t *target, kd_node_t **best, + int *limit, unsigned int coord) +{ + int dist = target->coords[coord] - root->coords[coord]; + int dist_sq = dist*dist; + + if (!root->removed) { + int root_dist_sq = kd_distance_sq(root, target); + if (!*best || root_dist_sq < *limit) { + *best = root; + *limit = root_dist_sq; + } + } + + coord = (coord + 1)%KD_DIMEN; + + if (root->left && (dist <= 0 || !*best || dist_sq <= *limit)) { + kd_find_nearest_recursive(root->left, target, best, limit, coord); + } + if (root->right && (dist >= 0 || !*best || dist_sq <= *limit)) { + kd_find_nearest_recursive(root->right, target, best, limit, coord); + } +} + +static void +kd_find_nearest(kd_node_t *root, kd_node_t *target, kd_node_t **best, + int *limit) +{ + kd_find_nearest_recursive(root, target, best, limit, 0); +} + +void +kdf_init(kd_forest_t *kdf) +{ + kdf->roots = NULL; + kdf->size = kdf->size_est = 0; + kdf->roots_size = kdf->roots_capacity = 0; +} + +void +kdf_destroy(kd_forest_t *kdf) +{ + free(kdf->roots); +} + +static size_t +kdf_collect_nodes(kd_forest_t *kdf, kd_node_t **buffer, unsigned int slot, + bool include_removed) +{ + size_t count = 0; + for (unsigned int i = 0; i < slot; ++i) { + if (kdf->roots[i]) { + count += kd_collect_nodes(kdf->roots[i], buffer + count, include_removed); + } + } + return count; +} + +static void +kdf_balance(kd_forest_t *kdf, kd_node_t *new_node, bool force) +{ + ++kdf->size; + + size_t slot, buffer_size; + if (force) { + buffer_size = kdf->size_est = kdf->size; + slot = kdf->roots_size; + } else { + ++kdf->size_est; + for (slot = 0; slot < kdf->roots_size; ++slot) { + if (!kdf->roots[slot]) { + break; + } + } + buffer_size = 1 << slot; + } + + kd_node_t **buffer = xmalloc(buffer_size*sizeof(kd_node_t *)); + buffer[0] = new_node; + kdf_collect_nodes(kdf, buffer + 1, slot, !force); + + kd_node_t **buffers[KD_BUFSIZE]; + for (int i = 1; i < KD_BUFSIZE; ++i) { + buffers[i] = xmalloc(buffer_size*sizeof(kd_node_t *)); + } + + if (slot >= kdf->roots_capacity) { + kdf->roots_capacity = slot + 1; + kdf->roots = xrealloc(kdf->roots, kdf->roots_capacity*sizeof(kd_node_t *)); + } + + size_t i, offset; + for (i = 0, offset = 0; offset < buffer_size; ++i) { + size_t chunk_size = 1 << i; + if (buffer_size & chunk_size) { + buffers[0] = buffer + offset; + kdf->roots[i] = kd_build_tree(buffers, chunk_size); + offset |= chunk_size; + } else { + kdf->roots[i] = NULL; + } + } + if (force || i > kdf->roots_size) { + kdf->roots_size = i; + } + + free(buffer); + for (i = 1; i < KD_BUFSIZE; ++i) { + free(buffers[i]); + } +} + +void +kdf_insert(kd_forest_t *kdf, kd_node_t *node) +{ + node->added = true; + + // If half or more of the nodes are removed, force a complete rebalance + bool force = (kdf->size_est + 1) >= 2*(kdf->size + 1); + kdf_balance(kdf, node, force); +} + +void +kdf_remove(kd_forest_t *kdf, kd_node_t *node) +{ + --kdf->size; + node->removed = true; +} + +kd_node_t * +kdf_find_nearest(kd_forest_t *kdf, kd_node_t *target) +{ + int limit; + kd_node_t *best = NULL; + + for (unsigned int i = 0; i < kdf->roots_size; ++i) { + kd_node_t *root = kdf->roots[i]; + if (root != NULL) { + kd_find_nearest(root, target, &best, &limit); + } + } + + return best; +} diff --git a/kd-forest.h b/kd-forest.h new file mode 100644 index 0000000..f10f982 --- /dev/null +++ b/kd-forest.h @@ -0,0 +1,57 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 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 KD_FOREST_H +#define KD_FOREST_H + +#include <stdbool.h> +#include <stddef.h> +#include <stdint.h> + +#define KD_DIMEN 3 + +// Single node in a k-d tree +typedef struct kd_node_t { + // Node coordinates + int coords[KD_DIMEN]; + // Sub-trees + struct kd_node_t *left, *right; + // Used to keep track of which sub-tree a node is in during construction + bool is_left; + // State flags + bool added, removed; + + // Corresponding image position for this node + unsigned int x, y; +} kd_node_t; + +void kd_node_init(kd_node_t *node, unsigned int x, unsigned int y); +void kd_node_set_color(kd_node_t *node, uint32_t color); + +// A forest of k-d trees +typedef struct { + // Array of k-d tree roots + kd_node_t **roots; + // Size and capacity of the roots array + unsigned int roots_size, roots_capacity; + // The actual size of this tree + size_t size; + // The size estimate for this tree, counting removed nodes + size_t size_est; +} kd_forest_t; + +void kdf_init(kd_forest_t *kdf); +void kdf_destroy(kd_forest_t *kdf); +void kdf_insert(kd_forest_t *kdf, kd_node_t *node); +void kdf_remove(kd_forest_t *kdf, kd_node_t *node); +kd_node_t *kdf_find_nearest(kd_forest_t *kdf, kd_node_t *target); + +#endif // KD_FOREST_H @@ -0,0 +1,281 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 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. * + *********************************************************************/ + +// Number of trailing zero bits on each color chanel, set to zero for all +// 24-bit images +#define ZERO_BITS 0 +// Whether to sort by hue +#define HUE_SORT 1 + +#if HUE_SORT +# define PASSES 16 +# define RANDOMIZE 0 +#else +# define PASSES 1 +# define RANDOMIZE 1 +#endif + +#include "kd-forest.h" +#include "util.h" +#include <errno.h> +#include <math.h> +#include <png.h> +#include <setjmp.h> +#include <stdio.h> +#include <stdlib.h> + +#if __unix__ +#include <unistd.h> +#endif + +unsigned int +rand_in(unsigned int range) +{ + // Compensate for bias if (RAND_MAX + 1) isn't a multiple of range + unsigned int limit = RAND_MAX + 1U - ((RAND_MAX + 1U)%range); + unsigned int res; + do { + res = rand(); + } while (res >= limit); + return res%range; +} + +kd_node_t * +try_neighbor(kd_node_t *node, unsigned int width, unsigned int height, + int which) +{ + int dx = which%3 - 1; + int dy = which/3 - 1; + + if (dx < 0 && node->x < -dx) { + return NULL; + } else if (dx > 0 && node->x + dx >= width) { + return NULL; + } else if (dy < 0 && node->y < -dy) { + return NULL; + } else if (dy > 0 && node->y + dy >= height) { + return NULL; + } + + return node + (int)width*dy + dx; +} + +kd_node_t * +next_neighbor(kd_node_t *node, unsigned int width, unsigned int height) +{ + unsigned int first = rand_in(9); + + for (unsigned int i = first; i < first + 9; ++i) { + int which = i%9; + if (which == 4) { + // Skip self + continue; + } + + kd_node_t *neighbor = try_neighbor(node, width, height, which); + if (neighbor && !neighbor->added) { + return neighbor; + } + } + + return NULL; +} + +void +remove_if_surrounded(kd_forest_t *kdf, kd_node_t *node, unsigned int width, + unsigned int height) +{ + if (node->added && !node->removed + && next_neighbor(node, width, height) == NULL) { + kdf_remove(kdf, node); + } +} + +void +remove_non_boundary(kd_forest_t *kdf, kd_node_t *node, unsigned int width, + unsigned int height) +{ + for (int i = 0; i < 9; ++i) { + kd_node_t *neighbor = try_neighbor(node, width, height, i); + if (neighbor) { + remove_if_surrounded(kdf, neighbor, width, height); + } + } +} + +#if HUE_SORT +#define PI 3.1415926535897932 + +static double +hue(uint32_t color) +{ + int R = (color >> 16) & 0xFF; + int G = (color >> 8) & 0xFF; + int B = color & 0xFF; + + double hue = atan2(sqrt(3.0)*(G - B), 2*R - G - B); + if (hue < 0.0) { + hue += 2.0*PI; + } + return hue; +} + +static int +hue_comparator(const void *a, const void *b) +{ + double ahue = hue(*(uint32_t *)a); + double bhue = hue(*(uint32_t *)b); + return (ahue > bhue) - (ahue < bhue); +} + +#endif + +int +main(void) +{ + const unsigned int jump = 1U << ZERO_BITS; + const unsigned int width = 1U << ((24 - 3*ZERO_BITS + 1)/2); // Round up + const unsigned int height = 1U << ((24 - 3*ZERO_BITS)/2); // Round down + const unsigned int size = width*height; + + printf("Generating a %ux%u image (%u pixels)\n", width, height, size); + + // Generate all the colors + uint32_t *colors = xmalloc(size*sizeof(uint32_t)); + for (unsigned int b = 0, i = 0; b < 0x100; b += jump) { + for (unsigned int g = 0; g < 0x100; g += jump) { + for (unsigned int r = 0; r < 0x100; r += jump, ++i) { + colors[i] = (r << 16) | (g << 8) | b; + } + } + } + srand(0); +#if RANDOMIZE + // Fisher-Yates shuffle + for (unsigned int i = size; i-- > 0;) { + unsigned int j = rand_in(i + 1); + uint32_t temp = colors[i]; + colors[i] = colors[j]; + colors[j] = temp; + } +#endif +#if HUE_SORT + qsort(colors, size, sizeof(uint32_t), hue_comparator); +#endif + + // Make a pool of potential k-d nodes + kd_node_t *nodes = xmalloc(size*sizeof(kd_node_t)); + for (unsigned int y = 0, i = 0; y < height; ++y) { + for (unsigned int x = 0; x < width; ++x, ++i) { + kd_node_init(nodes + y*width + x, x, y); + } + } + + // Make the forest + kd_forest_t kdf; + kdf_init(&kdf); + +#if __unix__ + bool tty = isatty(1); + const char *clear_line = tty ? "\033[2K\r" : ""; + const char *new_line = tty ? "" : "\n"; +#else + const char *clear_line = ""; + const char *new_line = "\n"; +#endif + + size_t max_size = 0; + + // Do multiple passes to get rid of artifacts in HUE_SORT mode + for (unsigned int i = 0, progress = 0; i < PASSES; ++i) { + for (unsigned int j = i; j < size; j += PASSES, ++progress) { + if (progress%width == 0) { + printf("%s%.2f%%\t| boundary size: %zu\t| max boundary size: %zu%s", + clear_line, 100.0*progress/size, kdf.size, max_size, new_line); + fflush(stdout); + } + + uint32_t color = colors[j]; + + kd_node_t target; + kd_node_set_color(&target, color); + + kd_node_t *new_node; + if (j == 0) { + // First node goes in the center + new_node = nodes + size/2 + width/2; + } else { + kd_node_t *nearest = kdf_find_nearest(&kdf, &target); + new_node = next_neighbor(nearest, width, height); + } + + kd_node_set_color(new_node, color); + kdf_insert(&kdf, new_node); + remove_non_boundary(&kdf, new_node, width, height); + + if (kdf.size > max_size) { + max_size = kdf.size; + } + } + } + printf("%s%.2f%%\t| boundary size: 0\t| max boundary size: %zu\n", + clear_line, 100.0, max_size); + + FILE *file = fopen("kd-forest.png", "wb"); + if (!file) { + abort(); + } + + png_structp png_ptr = + png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!png_ptr) { + abort(); + } + + png_infop info_ptr = png_create_info_struct(png_ptr); + if (!info_ptr) { + abort(); + } + + // libpng will longjmp here if it encounters an error from now on + if (setjmp(png_jmpbuf(png_ptr))) { + abort(); + } + + png_init_io(png_ptr, file); + png_set_IHDR(png_ptr, info_ptr, width, height, 8, + PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, + PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); + png_set_sRGB_gAMA_and_cHRM(png_ptr, info_ptr, PNG_sRGB_INTENT_ABSOLUTE); + + png_write_info(png_ptr, info_ptr); + + uint8_t *row = xmalloc(3*width*sizeof(uint8_t)); + for (unsigned int y = 0; y < height; ++y) { + for (unsigned int x = 0; x < width; ++x) { + kd_node_t *node = nodes + y*width + x; + row[3*x] = node->coords[0]; + row[3*x + 1] = node->coords[1]; + row[3*x + 2] = node->coords[2]; + } + png_write_row(png_ptr, (png_bytep)row); + } + + png_write_end(png_ptr, info_ptr); + + free(row); + png_destroy_write_struct(&png_ptr, &info_ptr); + fclose(file); + kdf_destroy(&kdf); + free(nodes); + free(colors); + return 0; +} @@ -0,0 +1,33 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 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 "util.h" +#include <stdlib.h> + +void * +xmalloc(size_t size) +{ + void *ret = malloc(size); + if (!ret) { + abort(); + } + return ret; +} + +void * +xrealloc(void *ptr, size_t size) +{ + void *ret = realloc(ptr, size); + if (!ret) { + abort(); + } + return ret; +} @@ -0,0 +1,20 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 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 UTIL_H +#define UTIL_H + +#include <stddef.h> + +void *xmalloc(size_t size); +void *xrealloc(void *ptr, size_t size); + +#endif // UTIL_H |