From 746fdfc554b5584476df93b6c8717d04814cc2d4 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Tue, 11 Mar 2014 19:51:47 -0400 Subject: Allow different color spaces to be used for similarity measurement. --- Makefile | 4 +- color.c | 119 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ color.h | 27 ++++++++++++++ kd-forest.c | 59 ++++++++++++------------------ kd-forest.h | 3 +- main.c | 68 +++++++++++++++++++++------------- 6 files changed, 216 insertions(+), 64 deletions(-) create mode 100644 color.c create mode 100644 color.h diff --git a/Makefile b/Makefile index 77bc5c9..582b209 100644 --- a/Makefile +++ b/Makefile @@ -14,8 +14,8 @@ CFLAGS ?= -std=c99 -pipe -O2 -Werror -Wall -Wpedantic -Wextra -Wno-sign-compare 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 +kd-forest: kd-forest.c kd-forest.h util.c util.h color.c color.h main.c + $(CC) $(CFLAGS) $(LDFLAGS) -o kd-forest kd-forest.c util.c color.c main.c image: kd-forest ./kd-forest diff --git a/color.c b/color.c new file mode 100644 index 0000000..34eab01 --- /dev/null +++ b/color.c @@ -0,0 +1,119 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 Tavian Barnes * + * * + * 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 "color.h" +#include + +void +color_unpack(uint8_t pixel[3], uint32_t color) +{ + pixel[0] = (color >> 16) & 0xFF; + pixel[1] = (color >> 8) & 0xFF; + pixel[2] = color & 0xFF; +} + +void +color_set_RGB(double coords[3], uint32_t color) +{ + uint8_t pixel[3]; + color_unpack(pixel, color); + for (int i = 0; i < 3; ++i) { + coords[i] = pixel[i]/255.0; + } +} + +// Inverse gamma for sRGB +double +sRGB_C_inv(double t) +{ + if (t <= 0.040449936) { + return t/12.92; + } else { + return pow((t + 0.055)/1.055, 2.4); + } +} + +static void +color_set_XYZ(double XYZ[3], uint32_t color) +{ + double RGB[3]; + color_set_RGB(RGB, color); + + RGB[0] = sRGB_C_inv(RGB[0]); + RGB[1] = sRGB_C_inv(RGB[1]); + RGB[2] = sRGB_C_inv(RGB[2]); + + XYZ[0] = 0.4123808838268995*RGB[0] + 0.3575728355732478*RGB[1] + 0.1804522977447919*RGB[2]; + XYZ[1] = 0.2126198631048975*RGB[0] + 0.7151387878413206*RGB[1] + 0.0721499433963131*RGB[2]; + XYZ[2] = 0.0193434956789248*RGB[0] + 0.1192121694056356*RGB[1] + 0.9505065664127130*RGB[2]; +} + +// CIE L*a*b* and L*u*v* gamma +static double +Lab_f(double t) +{ + if (t > 216.0/24389.0) { + return pow(t, 1.0/3.0); + } else { + return 841.0*t/108.0 + 4.0/29.0; + } +} + +// sRGB white point (CIE D50) in XYZ coordinates +static const double WHITE[] = { + [0] = 0.9504060171449392, + [1] = 0.9999085943425312, + [2] = 1.089062231497274, +}; + +void +color_set_Lab(double coords[3], uint32_t color) +{ + double XYZ[3]; + color_set_XYZ(XYZ, color); + + double fXYZ[] = { + [0] = Lab_f(XYZ[0]/WHITE[0]), + [1] = Lab_f(XYZ[1]/WHITE[1]), + [2] = Lab_f(XYZ[2]/WHITE[2]), + }; + + coords[0] = 116.0*fXYZ[1] - 16.0; + coords[1] = 500.0*(fXYZ[0] - fXYZ[1]); + coords[2] = 200.0*(fXYZ[1] - fXYZ[2]); +} + +void +color_set_Luv(double coords[3], uint32_t color) +{ + double XYZ[3]; + color_set_XYZ(XYZ, color); + + double uv_denom = XYZ[0] + 15.0*XYZ[1] + 3.0*XYZ[2]; + if (uv_denom == 0.0) { + coords[0] = 0.0; + coords[1] = 0.0; + coords[2] = 0.0; + return; + } + + double white_uv_denom = WHITE[0] + 16.0*WHITE[1] + 3.0*WHITE[2]; + + double fY = Lab_f(XYZ[1]/WHITE[1]); + double uprime = 4.0*XYZ[0]/uv_denom; + double unprime = 4.0*WHITE[0]/white_uv_denom; + double vprime = 9.0*XYZ[1]/uv_denom; + double vnprime = 9.0*WHITE[1]/white_uv_denom; + + coords[0] = 116.0*fY - 16.0; + coords[1] = 13.0*coords[0]*(uprime - unprime); + coords[2] = 13.0*coords[0]*(vprime - vnprime); +} diff --git a/color.h b/color.h new file mode 100644 index 0000000..6ce5168 --- /dev/null +++ b/color.h @@ -0,0 +1,27 @@ +/********************************************************************* + * kd-forest * + * Copyright (C) 2014 Tavian Barnes * + * * + * 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 COLOR_H +#define COLOR_H + +#include "kd-forest.h" +#include + +void color_unpack(uint8_t pixel[3], uint32_t color); + +// Use RGB coordinates +void color_set_RGB(double coords[3], uint32_t color); +// Use CIE L*a*b* coordinates +void color_set_Lab(double coords[3], uint32_t color); +// Use CIE L*u*v* coordinates +void color_set_Luv(double coords[3], uint32_t color); + +#endif // COLOR_H diff --git a/kd-forest.c b/kd-forest.c index 4a5bf37..5b230be 100644 --- a/kd-forest.c +++ b/kd-forest.c @@ -11,6 +11,7 @@ #include "kd-forest.h" #include "util.h" +#include #include #include #include @@ -24,14 +25,6 @@ kd_node_init(kd_node_t *node, unsigned int x, unsigned int 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) { @@ -52,21 +45,21 @@ kd_collect_nodes(kd_node_t *root, kd_node_t **buffer, bool include_removed) 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; + double aval = (*(const kd_node_t **)a)->coords[0]; + double bval = (*(const kd_node_t **)b)->coords[0]; + return (aval > bval) - (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; + double aval = (*(const kd_node_t **)a)->coords[1]; + double bval = (*(const kd_node_t **)b)->coords[1]; + return (aval > bval) - (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; + double aval = (*(const kd_node_t **)a)->coords[2]; + double bval = (*(const kd_node_t **)b)->coords[2]; + return (aval > bval) - (aval < bval); } static kd_comparator *kd_comparators[KD_DIMEN] = { @@ -80,8 +73,7 @@ static kd_comparator *kd_comparators[KD_DIMEN] = { #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) +kd_build_tree_recursive(kd_node_t **buffers[KD_BUFSIZE], size_t size, unsigned int coord) { if (size == 0) { return NULL; @@ -144,27 +136,26 @@ kd_build_tree(kd_node_t **buffers[KD_BUFSIZE], size_t size) return kd_build_tree_recursive(buffers, size, 0); } -static int +static double kd_distance_sq(kd_node_t *a, kd_node_t *b) { - int result = 0; + double result = 0.0; for (int i = 0; i < KD_DIMEN; ++i) { - int d = a->coords[i] - b->coords[i]; + double 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) +kd_find_nearest_recursive(kd_node_t *root, kd_node_t *target, kd_node_t **best, double *limit, unsigned int coord) { - int dist = target->coords[coord] - root->coords[coord]; - int dist_sq = dist*dist; + double dist = target->coords[coord] - root->coords[coord]; + double dist_sq = dist*dist; if (!root->removed) { - int root_dist_sq = kd_distance_sq(root, target); - if (!*best || root_dist_sq < *limit) { + double root_dist_sq = kd_distance_sq(root, target); + if (root_dist_sq < *limit) { *best = root; *limit = root_dist_sq; } @@ -172,17 +163,16 @@ kd_find_nearest_recursive(kd_node_t *root, kd_node_t *target, kd_node_t **best, coord = (coord + 1)%KD_DIMEN; - if (root->left && (dist <= 0 || !*best || dist_sq <= *limit)) { + if (root->left && (dist <= 0 || dist_sq <= *limit)) { kd_find_nearest_recursive(root->left, target, best, limit, coord); } - if (root->right && (dist >= 0 || !*best || dist_sq <= *limit)) { + if (root->right && (dist >= 0 || 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(kd_node_t *root, kd_node_t *target, kd_node_t **best, double *limit) { kd_find_nearest_recursive(root, target, best, limit, 0); } @@ -202,8 +192,7 @@ kdf_destroy(kd_forest_t *kdf) } static size_t -kdf_collect_nodes(kd_forest_t *kdf, kd_node_t **buffer, unsigned int slot, - bool include_removed) +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) { @@ -288,7 +277,7 @@ kdf_remove(kd_forest_t *kdf, kd_node_t *node) kd_node_t * kdf_find_nearest(kd_forest_t *kdf, kd_node_t *target) { - int limit; + double limit = INFINITY; kd_node_t *best = NULL; for (unsigned int i = 0; i < kdf->roots_size; ++i) { diff --git a/kd-forest.h b/kd-forest.h index f10f982..75d8666 100644 --- a/kd-forest.h +++ b/kd-forest.h @@ -21,7 +21,7 @@ // Single node in a k-d tree typedef struct kd_node_t { // Node coordinates - int coords[KD_DIMEN]; + double 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 @@ -34,7 +34,6 @@ typedef struct kd_node_t { } 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 { diff --git a/main.c b/main.c index 4a9bfcd..cd08967 100644 --- a/main.c +++ b/main.c @@ -15,16 +15,23 @@ // Whether to sort by hue #define HUE_SORT 1 +// Which color space to use +#define USE_RGB 0 +#define USE_LAB 1 +#define USE_LUV 0 + #define RANDOMIZE (!HUE_SORT) #include "kd-forest.h" #include "util.h" +#include "color.h" #include #include #include #include #include #include +#include #if __unix__ #include @@ -43,8 +50,7 @@ rand_in(unsigned int range) } kd_node_t * -try_neighbor(kd_node_t *node, unsigned int width, unsigned int height, - int which) +try_neighbor(kd_node_t *node, unsigned int width, unsigned int height, int which) { int dx = which%3 - 1; int dy = which/3 - 1; @@ -84,8 +90,7 @@ next_neighbor(kd_node_t *node, unsigned int width, unsigned int height) } void -remove_if_surrounded(kd_forest_t *kdf, kd_node_t *node, unsigned int width, - unsigned int height) +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) { @@ -94,8 +99,7 @@ remove_if_surrounded(kd_forest_t *kdf, kd_node_t *node, unsigned int width, } void -remove_non_boundary(kd_forest_t *kdf, kd_node_t *node, unsigned int width, - unsigned int height) +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); @@ -165,6 +169,12 @@ main(void) qsort(colors, size, sizeof(uint32_t), hue_comparator); #endif + // Make the actual bitmap image + png_bytepp rows = xmalloc(height*sizeof(png_bytep)); + for (unsigned int i = 0; i < height; ++i) { + rows[i] = xmalloc(3*width*sizeof(png_byte)); + } + // 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) { @@ -203,7 +213,15 @@ main(void) uint32_t color = colors[j]; kd_node_t target; - kd_node_set_color(&target, color); +#if USE_RGB + color_set_RGB(target.coords, color); +#elif USE_LAB + color_set_Lab(target.coords, color); +#elif USE_LUV + color_set_Luv(target.coords, color); +#else +# error "Pick one!" +#endif kd_node_t *new_node; if (j == 0) { @@ -211,16 +229,25 @@ main(void) new_node = nodes + size/2 + width/2; } else { kd_node_t *nearest = kdf_find_nearest(&kdf, &target); + if (!nearest) { + abort(); + } new_node = next_neighbor(nearest, width, height); + if (!new_node) { + abort(); + } } - kd_node_set_color(new_node, color); + memcpy(new_node->coords, target.coords, sizeof(target.coords)); kdf_insert(&kdf, new_node); remove_non_boundary(&kdf, new_node, width, height); if (kdf.size > max_size) { max_size = kdf.size; } + + png_bytep pixel = rows[new_node->y] + 3*new_node->x; + color_unpack(pixel, color); } } printf("%s%.2f%%\t| boundary size: 0\t| max boundary size: %zu\n", @@ -249,28 +276,19 @@ main(void) png_init_io(png_ptr, file); png_set_IHDR(png_ptr, info_ptr, width, height, 8, - PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, + PNG_COLOR_TYPE_RGB, PNG_INTERLACE_ADAM7, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); png_set_sRGB_gAMA_and_cHRM(png_ptr, info_ptr, PNG_sRGB_INTENT_ABSOLUTE); + png_set_rows(png_ptr, info_ptr, rows); + png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL); + png_destroy_write_struct(&png_ptr, &info_ptr); + fclose(file); - 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); + for (unsigned int i = 0; i < height; ++i) { + free(rows[i]); } + free(rows); - 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); -- cgit v1.2.3